Molecular Dynamics Study on Structure, Vibrational Properties, and Transport Coefficients of Liquid Alumina

The structure, vibrational density of states, and transport coefficients of liquid alumina were studied using molecular dynamics simulations. At the temperature of 2500 K, 3000 K, 3500 K, and 4000 K, systems with three different densities were constructed, respectively, including the configurations with densities of 2.81 g/cm3 and 3.17 g/cm3, and the relaxed ones with nearly zero pressure at each temperature. With the changes in temperature or density, the transformations on the structural, vibrational and transport properties were discussed. The Born–Mayer–Huggins type of atomic interactions was used, with newly optimized parameters. The analysis of the interatomic correlations indicated that the short-range order of liquid alumina was mainly constructed by AlO4 tetrahedra, also a certain number of AlO3 and AlO5 was present. Meanwhile, the structural transitions on the elemental units occurred as either the temperature or density increased. Two primary frequency bands were observed in each vibrational density of states spectrum, with the higher frequency bands produced by the O atom vibrations, and the lower frequency ones generated by the Al atom vibrations. Self-diffusion coefficients were estimated using the linear behavior of the mean-squared displacement at long time, while by using the Green–Kubo relation during equilibrium molecular dynamics simulations, thermal conductivities and viscosities were calculated. Significantly, the viscosity at 2500 K with a density of 2.81 g/cm3 was equal to 25.23 mPa s, which was very close to the experimental finding.

Understandings of the origins of the outstanding macroscopic features can be promoted by the knowledge of the microscopic structures. The high melting point of alumina (up to 2327 K) [20], however, frequently caused contamination from the container throughout the experiment, reducing the accuracy of the obtained structural information. With the use of X-ray diffraction [21][22][23], neutron scattering measurements [24,25], nuclear magnetic resonance measurements (NMR) [26][27][28], and aerodynamic-levitation and laser-heating techniques [29], some processes in the structural investigations of molten alumina have been achieved. Most of the experiments obtained the coordination numbers of Al-O pairs between 4 and 5. However, there were also some disagreements; for example, the measurement conducted by Waseda et al. [22] using X-ray diffraction at 2363 K concluded that the coordination number Z Al-O was equal to 5.6, indicating that the elemental units for the construction of liquid alumina were octahedral networks. Besides the structural information, the physical properties of molten alumina, for example, electrical conductivity [30], heat capacity [31], surface tension [32,33], emissivity [34], and viscosity [32], were investigated by non-contract measurements.
However, the experimental data were usually inadequate and difficult to reproduce because of the challenges experienced during the experiment processes. An alternative method was to estimate the characteristics of liquid alumina using first-principle calculations or molecular dynamics (MD) simulations. While it was discovered that liquid alumina was mostly made up of AlO 4 tetrahedrons [29,[35][36][37][38][39], others suggested that AlO 5 may have also been present in significant amounts [40][41][42]. Even Hemmati [43] claimed that liquid alumina was mostly composed of AlO 5 units. It was also evident from the MD simulations that the pore aggregation underwent a continuous transition as the pressure increased [41,[44][45][46]. Additionally, by using the slopes of the mean-squared displacements, the self-diffusion coefficients were determined [36,37,40,41]. The specific heat, vibrational density of states, and ionic conductivity were examined at the temperature of 2600 K basing on an atomic interaction consisting of two-and three-body terms by Vashishta et al. [35]. Meanwhile, the charge-transport properties of liquid alumina were discussed in detail by Gheribi et al. [47].
It is worth noting that most computational research on liquid alumina was performed at high or ambient pressures, and the relaxation configurations at particularly zero pressure have been less well understood. Furthermore, the earlier studies on liquid alumina have primarily aimed at the structural characteristics, with only a small number of them being concerned with the physical properties, particularly given the fact that a detailed exploration of transport properties has rarely been reported.
Therefore, in this study, the configurations with nearly zero or a certain pressure were all investigated. Furthermore, besides the structure and vibrational density of states, the transport coefficients of liquid alumina were carefully examined in this research. Using the newly optimized parameters [48], the examinations were conducted by classical MD simulations based on Born-Mayer-Huggins (BMH) atomic interactions. Systems with three different densities were established at temperatures of 2500 K, 3000 K, 3500 K, and 4000 K, respectively. At each temperature, the examinations were performed on not only the configurations with densities of 2.81 g/cm 3 and 3.17 g/cm 3 (with some pressures), but also the relaxed ones with nearly zero pressure. Moreover, transformations in structural, vibrational, and transport characteristics were discussed in relation to the changes in both density and temperature.

Simulation Methods
2000 atoms (1200 O and 800 Al) were packed randomly into a basic cube, creating a system with periodic boundary conditions in x, y and z directions. Using the Verlet algorithm in the velocity form with a timestep of 1.0 × 10 −3 ps, Newton's equations of motion were integrated. The short-range interaction was cut off at 8.0 Å, and the Coulomb interaction was calculated by the Ewald summation method. After being equilibrated in the isobaric-isothermal ensemble (NPT) at zero pressure and the temperature of 4500 K, the initial configuration was cooled down to 2500 K with a cooling rate of 1 K/ps. The corresponding configuration was taken at intervals of 500 K, and it was equilibrated for 1000 ps at a constant temperature and zero pressure. Since the object of this study was liquid alumina, the simulation conducted temperature was required to be higher than the melting point of 2327 K [20]. Therefore, the simulations were carried out at 2500 K, 3000 K, 3500 K, and 4000 K, respectively. Furthermore, in order to compare the structure, vibrational density of states and transport coefficients of liquid alumina as a function of density, two configurations with densities of 2.81 g/cm 3 and 3.17 g/cm 3 were produced by adjusting the volume of the simulation box at each temperature. In particular, the density of 2.81 g/cm 3 was experimentally recorded from neutron scattering measurements [25], while many previous simulations of liquid alumina were also performed on the density of 3.17 g/cm 3 , which could be utilized to check the reliability of the simulation results in this paper by the comparisons with the reported results, including those obtained from both experiments and simulations. Moreover, at each temperature, configurations with densities of 2.81 g/cm 3 and 3.17 g/cm 3 as well as those with anticipated densities at nearly zero pressure were investigated. The snapshots of the considered configurations at different densities and temperatures are shown in Figure 1. Classical MD simulations were carried out with the LAMMPS code [49] (version: stable_3Mar2020) by applying the Born-Mayer-Huggins (BMH) potential, and atomic interactions between O and Al particles were provided by with α or β = O or Al, and r αβ being the distance between the centers of the particle α and the particle β. q α and q β are the effective charges in the first term of the pair potential, which is associated with the long-range Coulomb. The Born repulsive is shown in the second term, while the dipolar expansion corresponds to the third and the fourth ones, though just the van der Waals term was calculated. Recently, Bouhadja et al. [48] adjusted the parameters of A αβ , σ αβ , ρ αβ, and C αβ in Equation (1), as displayed in Table 1.
The atomic correlations were used to determine the short-range order of the structural characteristics of liquid alumina, including the pair distribution function (g(r)) and the coordination number as well as the bond-angle distribution. These methods were commonly employed for the analysis of microscopic structures and have been proven to be efficient [50][51][52][53][54]. The partial pair distribution function (g αβ (r), PPDF) is shown as Equation (2): where N β is the number of the particle β, V represents the volume of the system, n αβ (r, r + ∆r) represents the average number of the β particles spanning the space from r to r+∆r. The total pair distribution function (PDF) was calculated as Equation (3): By using the minimum following the first peak of g αβ (r) as a cutoff radius (R), which was taken as R Al-O = 2.38 Å, R Al-Al = 3.76 Å and R O-O = 3.65 Å, the coordination numbers Z αβ (R), characterizing the mean number of β particles surrounding an α particle, were calculated as Equation (4): Regarding the dynamics, the vibrational density of states (VDOS) and transport coefficients (self-diffusion coefficient, thermal conductivity and viscosity) were estimated. For each configuration, five independent simulations with various beginning momenta were performed to decrease the statistical errors, since the transport coefficients are sensitive to the initial conditions. Additionally, the results were calculated by averaging the standard deviations over the individual simulations.
By conducting a Fourier transform on the velocity autocorrelation function (Z(t)), VDOS was determined as Equation (5): Figure 1. Snapshots of the considered configurations in this study at the temperatures of 2500 K, 3000 K, 3500 K, and 4000 K, with three different densities at each temperature (the blue particles: Al atoms, the red particles: O atoms). Table 1. Potential parameters for liquid alumina [48]. The self-diffusion coefficient (D), in this study, was computed from the linear behavior of the mean-squared displacement (MSD) over a long period, as displayed in Equation (6). This method has been widely used in the calculations of self-diffusion coefficients for many different systems, such as drugs [50], gases [55,56] nanoparticles [57,58], and liquids [59,60]. Moreover, the obtained results often showed to be stable due to the long correlation time.
where r α (t) denotes the vector of the position related to particle α when at time t, and N denotes the number of the particle α.
Since liquid alumina was a homogeneous system, thermal conductivity (κ) and viscosity (η) could be estimated by the equilibrium molecular dynamics (EMD) simulations with the Green-Kubo relation. The integral of time of the heat-current autocorrelation function (HCACF) was used to calculate the thermal conductivity, while that of the stress autocorrelation function (SACF) was used to evaluate the viscosity. The definition of the heat-current (J) is: where v α denotes the velocity of the particle α, e α represents the potential and kinetic energy of the particle α, and S α denotes the atomic stress of the particle α. Subsequently, using the Green-Kubo relation, EMD simulations were carried out to compute κ: where k B represents the Boltzmann constant, and T means the equilibrium temperature. The component of molecular stress tensor P ij for SACF (ij = xy, xz, and yz directions) is defined as: where m α is the mass of the particle α, v αi and v αj are the velocities of the particle α, r αβ denotes the relative position of the particle α about particle β, and F αβ represents the force operating on the particle α as a result of the interaction with particle β. Then η was computed by EMD simulation as Equation (10):

Results and Discussion
In this study, the calculated zero-pressure densities were equal to 2.86 g/cm 3 , 2.77 g/cm 3 , 2.67 g/cm 3 , and 2.56 g/cm 3 , respectively, at 2500 K, 3000 K, 3500 K, and 4000 K, decreasing with the increasing temperature. Notably, the predicted density values at 2500 K and 3000 K were consistent with the experimental values of 2.81 g/cm 3 and 2.78 g/cm 3 , which were measured, respectively, by the neutron scattering measurements [25] and aerodynamiclevitation technique [29], indicating the efficiency of the optimized parameters of the BMH potential.

Structural Properties
At 2500 K, 3000 K, 3500 K, and 4000 K, the position (r αβ ) and the full width at half maximum (FWHM) of the initial peak in each PPDF for the configurations with three different densities are shown in Table 2 Table 2. Simultaneously, for the configurations at zero pressure or with a constant density (ρ = 2.81 g/cm 3 or ρ = 3.17 g/cm 3 ), changes in r αβ also happened as the temperature increased, with r Al-Al increasing and r Al-O decreasing. Regarding r O-O , it increased along with the temperature for the zero-pressure densities, whereas it was found to be declined when at a constant density. Additionally, for all the Al-O, Al-Al, and O-O correlations, the FWHM of the initial peaks grew along with the temperature. Compared to those at 3000 K, the intensities of the initial peaks of all the g αβ (r) were shown to be lower at 4000 K.   Table 3 demonstrated that for the atom pairs of Al-O, the average coordination numbers (Z αβ ) were approximately located at~4, consistent with the experimental, ab initio and MD simulation findings. For all the atom pairs of α-β, Z αβ increased along with the density at each temperature. However, when at zero pressure or a fixed density, Z αβ roughly decreased as the temperatures rose from 2500 K to 4000 K.   (Figure 3a) were distributed between 2 and 6. Besides AlO 4 , there was also a certain amount of AlO 3 and AlO 5 present. With the temperature increasing, the fractions of AlO 2 and AlO 3 increased, whereas those of AlO 4 , AlO 5 and AlO 6 decreased, demonstrating a structural transition on AlO n . Moreover, when at a fixed temperature, some changes in the coordination number on AlO n also can be observed as the density rose, with a general decrease in the fractions of AlO 3 and AlO 4 , accompanying an increase in those of AlO 5 , indicating the trend of structural transition on the elemental units from a tetrahedral network to an octahedral network. These variations in the percentages of AlO n enabled the Z Al-O to decrease with increasing temperature, and to increase along with density, in accordance with the findings in Table 3. The coordination numbers of O-Al (Figure 3b) also had a narrow distribution, with a range between 1 and 5. In contrast, the coordination numbers for O-O (Figure 3c) and Al-Al (Figure 3d) exhibited broader ranges, with 4 to 15, and 3 to 13, respectively. The bond-angle distributions, as seen in Figure 4, revealed details about the polyhedral geometry and packing within the structure of liquid alumina. The bond-angel distributions of O-Al-O and Al-O-Al, respectively, can be applied to examine the intraand inter-polyhedral connectivity. According to Skinner's experiment [29], the presence of AlO 4 was associated with peaks around 101 • and 106 • in the bond-angle distribution of O-Al-O, while AlO 5 corresponded to the main peak at 86 • , following a minor one over the degrees of 140-170 • . As shown in Figure 4b (Figure 3a). Meanwhile, as the temperature increased, the initial peaks also changed to lower degrees, which might be due to the relatively high AlO 3 percentages. The Al-O-Al bond-angle distributions at 3000 K and 4000 K are shown in Figure 4a. At both 3000 K and 4000 K, there were two main peaks in each Al-O-Al bond-angle distribution: the one approximately distributed at 93 • corresponded to the linkages of AlO 5 -AlO 5 or AlO 5 -AlO 4 by edge-sharing [29], while the other one located at 119 • was attributed to the oxygen atoms that were three-fold coordinated, and linked by their corners to AlO 4 or AlO 5 [29]. The intensity ratio between the peak around 93 • and the one at 119 • (I 93 • /I 119 • ) increased along with the density or temperature, indicating the increasing fractions of the edge-sharing networks. Figure 5 displays the partial and total VDOSs for the configurations of liquid alumina with three densities, respectively, at 3000 K and 4000 K. It is clear that the VDOS spectra at 3000 K ranged from 0 to 50 THz, while those at 4000 K had wider ranges, up to about 55 THz. There were two principal bands in each VDOS spectrum, divided at frequencies between 15 and 20 THz. In accordance with the findings of Vashishta et al. [35], the Al atom vibrations largely produced the band with lower frequencies, whereas the O atom vibrations provided the band with higher frequencies. Additionally, similar to those for other amorphous networks that were predominantly made up of basic tetrahedrons [63], both the modes of bond-bending and vibrations of inter-tetrahedra were related to the lower frequency bands, while the bonding-stretching modes and vibrations of intra-tetrahedra were associated with the higher frequency bands [64].

Self-Diffusion Coefficient
The time dependence of MSDs and the density dependence of self-diffusion coefficients at temperatures of 2500 K, 3000 K, 3500 K, and 4000 K are illustrated in Figure 6. Al atoms moved more quickly than O atoms when they were present at the same density and temperature, as evidenced by the self-diffusion coefficients of Al atoms being slightly larger than that of O atoms, as shown in Table 4. The distinct diffusion mechanisms that occurred between Al and O particles could be able to explain this phenomenon. In addition to being in the free states, the breakdown or reformation of Al-O bonds influenced the diffusions of Al particles, while the transitions between the bridging and the nonbridging, often needing higher energies, affected the diffusions of O particles [65]. Table 4 displays the numerical values for the self-diffusion coefficients. Significantly, the results for the configuration at 2500 K and 2.81 g/cm 3 agreed well with those calculated by Jahn et al. [66] using MD simulations based on an advanced ionic model type potential. For the models at zero pressure or with a fixed density, the self-diffusion coefficients increased along with the temperature. Meanwhile, it also can be observed that the self-diffusion coefficients generally decreased with density increasing. Table 4. Self-diffusion coefficients for liquid alumina at the temperature ranging from 2500 K to 4000 K.

Thermal Conductivity
The normalized heat-current autocorrelation functions (HCACFs) and thermal conductivities (κ) of liquid alumina are given as functions of correlation time in Figures 7 and 8, respectively, at temperatures of 2500 K, 3000 K, 3500 K and 4000 K. Arima et al. reported that the thermal conductivities would be better estimated when the corresponding normalized HCACFs settled to zero or a constant value [67]. The normalized HCACFs in this investigation gradually decreased to negative until ultimately decaying to zero (Figure 7), confirming the reliability of the simulated thermal conductivity results.  It took a few picoseconds, as shown in Figure 8, for the curves of κ to stabilize at a specified value, which was accepted to be the value of thermal conductivity. As the influence of phonons grew with the temperature, the curves of κ required less time to reach the constant values, and the correlations dissipated more quickly as well. These findings matched those reported by Alexander et al. [68], who investigated the thermal conductivities for the actinide dioxides of UO 2 , and ThO 2 using the EMD calculations. Thus, by averaging the related curves from 6 to 12 ps at 2500 K, from 3 to 11 ps at 3000 K, from 3 to 9 ps at 3500 K and from 3 to 8 ps at 4000 K, the thermal conductivity values were determined. When at a fixed temperature, κ enhanced with density increasing. Meanwhile, it raised along with temperature while at the density of 2.81g/cm 3 or 3.17 g/cm 3 , but declined as the temperature increased for the zero-pressure densities.

Viscosity
The normalized stress autocorrelation functions (SACFs) and viscosities (η) for the configurations of liquid alumina at temperatures from 2500 K to 4000 K are shown in Figures 9 and 10, respectively. The normalized SACF curves fell monotonically and quickly declined to zero, taking just 6 ps at 2500 K, 2.5 ps at 3000 K, 1.5 ps at 3500 K, and even less than 1 ps at 4000 K. Similar to the cases of the thermal conductivities, the time taken by the viscosities curves to reach a stable value decreased with the temperature increasing, along with the decline of the duration at this stable value. Therefore, the viscosity values of η were achieved by averaging the values between 11 and 18 ps at 2500 K, between 3 and 9 ps at 3000 K, between 2 and 7 ps at 3500 K, and between 1 and 5 ps at 4000 K. As can be seen, the viscosity values increased along with density at each temperature. Moreover, for the configurations at zero pressure or with a constant density, they declined quickly when the temperature increased. Moreover, it was worth noting that the predicted viscosity value of 25.23 mPa·s at 2500 K with the density of 2.81g/cm 3 , was in good accordance with the experimental value of 25.6 mPa·s given by Paradis et al. [32], who tested the viscosities of molten alumina by the electrostatic levitation and multi-beam radiative heating techniques. Furthermore, the calculated viscosities of the models at 2500 K with densities of 2.81 g/cm 3 and 2.86 g/cm 3 , as well as those at 3000 K with densities of 2.77 g/cm 3 and 2.81 g/cm 3 (Figure 10), were all very close to the results obtained by the density-function theory based on electronic structure calculations [69].

Conclusions
Using molecular dynamics simulations, the structural, vibrational, and transport characteristics of liquid alumina were investigated, with the temperature increasing from 2500 K to 4000 K. Configurations with densities of 2.81 g/cm 3 and 3.17 g/cm 3 , as well as the relaxed ones with nearly zero pressure, were explored at intervals of 500 K. Transformations in structural, vibrational, and transport characteristics were discussed in relation to the changes in temperatures or densities.
The newly optimized parameters associated with the Born-Mayer-Huggins potential with were taken to calculate the atomic interactions. Although they were originally developed for aluminosilicate (AS), the findings of this paper proved their efficiency in characterizing the structure and some dynamical features of the pure liquid alumina.
The structure of liquid alumina was primarily made up of the slightly disordered AlO 4 tetrahedral units, with tiny fractions of AlO 2 and AlO 6 , as well as a certain amount of AlO 3 and AlO 5 . Structural transitions on the elemental units occurred as either the temperature or density increased. The explorations of the intermediate-range order revealed the two primary approaches that the elemental units being connected to each other: sharing edge or corner.
By applying the Fourier transform to the velocity autocorrelation function, the vibrational density of states (VDOS) was calculated, which presented two principal frequency bands, one with the higher frequency bands from the vibrations of O atoms, relating to the bonding-stretching modes and vibrations of intra-tetrahedra, while the other with the lower frequency ones from the vibrations of Al atoms, associating with the modes of bond-bending and vibrations of inter-tetrahedra.
The self-diffusion coefficient (D) was estimated by the linear behavior of the meansquared displacement over a long time, and by using the Green-Kubo relation, the thermal conductivity (κ) and viscosity (η) were calculated through the EMD simulations. When at a fixed temperature, the thermal conductivities and viscosities increased with the density, although the self-diffusion coefficients for both Al and O particles decreased. Meanwhile, changes in temperatures also affected the transport coefficients. For the models with the density of 2.81g/cm 3 or 3.17 g/cm 3 , with the increasing temperature, D and κ were discovered to be enhanced, whereas η declined, while for the zero-pressure models, D also increased along with the temperature, but κ and η decreased. Notably, the viscosity at 2500 K with the density of 2.81 g/cm 3

Data Availability Statement:
The parameters for the Born-Mayer-Huggins potential used in this manuscript were reported in Bouhadja's paper, defined in the references section as [48]. New data are available from the corresponding authors upon reasonable request.