Atomic dynamics of alumina melt : A molecular dynamics simulation study

The atomic dynamics of Al2O3 melt are studied by molecular dynamics simulation. The particle interactions are described by an advanced ionic interaction model that includes polarization effects and ionic shape deformations. The model has been shown to reproduce accurately the static structure factors S(Q) from neutron and x-ray diffraction and the dynamic structure factor S(Q, ω) from inelastic x-ray scattering. Analysis of the partial dynamic structure factors shows inelastic features in the spectra up to momentum transfers, Q, close to the principal peaks of partial static structure factors. The broadening of the Brillouin line widths is discussed in terms of a frequency dependent viscosity η(ω).


Introduction
In recent years, substantial progress has been achieved in understanding the atomic dynamics of liquids.This is mainly due to two achievements, i.e., the development of inelastic x-ray scattering (IXS) as complementary technique to inelastic neutron scattering (INS) and the increasing power of computer simulations.It is well known that the interatomic potential has an increasing effect on the dynamics when time scales of atomic interactions and length scales of interatomic distances are approached [1,2].Experimentally, a detailed information about the interatomic forces can be obtained from measurements of the dynamic structure factor S(Q, ω), which represents the collective atomic dynamics, using inelastic neutron or x-ray scattering techniques.Q is the magnitude of the wave vector transfer and ω is the angular frequency, ω the respective energy transfer with being Planck's constant divided by 2π.Propagating acoustic-like modes have been observed in many liquids up to the terahertz frequency range, but they are much more pronounced and persistent in metallic than in molecular systems, which has been shown to be a consequence of different types of interaction [3].
For a long time the experimental observation of collective excitations in ionic liquids was very complicated or impossible due to the large sound velocity and kinematic restrictions of INS [4,5].These restrictions were lifted when IXS was developed [6].Since then the collective excitations in S(Q, ω) have been measured, e.g. in molten NaCl [7] and molten Al 2 O 3 [8], the latter in a containerless levitation experiment.On the other hand, the increasing computing power allows us to perform computer simulations of a new quality.The latter is essential since, especially in the field of liquids, the combination of experiment and simulation has always been a very fruitful symbiosis towards the understanding of atomistic processes.
If a realistic modelling of the atomic dynamics of liquids is attempted, an accurate description of the interatomic forces is crucial.The simplest representation of atomic interactions via pair potentials, which gives a good description of molecular systems in a wide range of thermodynamic states, cannot be easily applied to ionic materials.The reason is that the electronic structure of an ion in a condensed phase strongly depends on the interaction with its neighbours and hence manybody effects have to be considered.Rigid ion models of Born-Mayer-type, for example, subsume the many-body effects in some average sense and may be capable of reproducing some of the experimental results.However, they often do not provide a reliable description of both static and dynamic properties.Other effects such as polarization have to be included into an ionic model to obtain a more realistic description of disordered ionic systems, which has been shown for halide melts [9] or amorphous silica [10].Therefore, various polarizable ion models have been developed in order to model the properties of ionic solids and melts.However, most of them are optimized for a specific system and possess just a very limited transferability.Recently, a systematic approach to the construction of more transferable ionic potentials was described [11,12].The so-called Aspherical Ion Model (AIM) includes explicitly all components of the interaction that are considered essential.Its parameters are optimized using first-principles electronic structure calculations [11,12].
The work presented here is inspired by recent experimental studies of the atomic structure [13][14][15] and dynamics [8] of Al 2 O 3 melt.Most of the previous simulation studies of alumina melt (e.g.[16][17][18]) were concerned with structural properties at ambient pressure.Hoang and coworkers used a rigid ion model in order to study the self-diffusion [19] and the high pressure behaviour of the melt [20].An advanced ionic model of AIM-type [21] was recently used to perfectly reproduce both the experimental static and dynamic structure factors of alumina melt [22].Here, we focus on the relations between the spectra of collective modes, the vibrational spectra and the structural relaxation processes in the melt.

Interaction potential
Details of the AIM potential and its parametrization have been described earlier [21].The Al 2 O 3 potential is not only applicable to the melt [22] but it also satisfactorily describes the properties of solid alumina phases [21,23].The AIM interatomic potential V is constructed from of components The first two components, the charge-charge (V qq ) and dispersion (V disp ) interactions, are a purely pairwise additive: q i is the formal charge on ion i (+3 for Al and -2 for O), C ij 6 and C ij 8 are the dipole-dipole and dipole-quadrupole dispersion coefficients respectively, and f ij n are Tang-Toennies dispersion damping functions [24], which describe short-range corrections to the asymptotic dispersion term.
The overlap repulsion component (V rep ) is given by with and summation of repeated indexes is implied.The variable δσ i characterizes the deviation of the radius of oxide anion i from its default value, {ν i α } are a set of three variables describing the αβ = 3r ij,α r ij,β /r 2 ij − δ αβ are interaction tensors.The last summations include the self-energy terms, representing the energy required to deform the anion charge density, with β, ζ and η as effective force constants.The extent of each ion's distortion is determined at each molecular dynamics time-step by energy minimization.
The polarization part of the potential incorporates dipolar and quadrupolar contributions [25], αβ + where , α, B and C the dipole, dipole-dipole-quadrupole and quadrupole polarizabilities of the oxygen ion, respectively, and rij are the multipole interaction tensors [26].The instantaneous values of these moments are obtained by minimization of this expression.The charge-dipole and charge-quadrupole cation-anion asymptotic functions include terms which account for penetration effects at short-range by using Tang-Toennies damping functions [24] of the form, with D and Q standing for dipolar and quadrupolar parts.While the parameters b D and b Q determine the range at which the overlap of the charge densities affects the induced multipoles, the parameters c D and c Q determine the strength of the ion response to this effect.Only oxygen ions are considered polarizable.A compilation of the potential parameters is given in table 1.

Molecular dynamics simulation
Due to the "electronic" degrees of freedom of the AIM potential, the molecular dynamics simulations include additional annealing in the spirit of Born-Oppenheimer dynamics.Induced multipoles and ionic deformations are obtained at each time step by conjugate gradient minimization of the respective contributions to the total energy before forces acting on individual ions are calculated.The molecular dynamics simulations are performed with cubic simulation cells containing 2160 ions (432 formula units) and a constant temperature T of 2500 K.The cell size is significantly larger than in the previous study [22] to access smaller Q and to make better contact to the IXS experimental data.The smallest Q from the simulation is determined by the simulation box length, here it is 0.21 Å−1 .In the simulation, the system is first equilibrated for 10 ps in the NPT ensemble using an isotropic barostat coupled to a Nosé-Hoover thermostat [27].Thereafter the barostat is switched off and simulations continue in the NVT ensemble.The temperature is still controlled by the Nosé-Hoover thermostat [28] in order to correct for a small energy drift at long simulation times.Trajectories of the production runs are collected for 200 ps.

Vibrational spectra and self-diffusion coefficients
Single particle time correlation functions provide information regarding the distribution of vibrational frequencies in the melt.The velocity autocorrelation functions v X 1 (0)v X 1 (t) for X = Al and O are shown in figure 1.After a fast initial decay the typical oscillatory behaviour is observed.The correlation is essentially lost after a few tenth of picoseconds.Fourier transformation of the velocity autocorrelation functions defines a power spectrum of vibrational frequencies Both z Al (ω) and z O (ω) show a broad distribution with vibrational energies up to about 200 meV (see figure 1).Besides the principal peak at about 35 meV, there are two shoulders at around 75 and 100 meV.Around those frequencies, weakly dispersive peaks are observed at low Q in the transverse and longitudinal charge current correlation spectra, respectively.They are reminiscent of optic modes and seem to make a significant contribution to the total vibrational density of states.The zero frequency limit of z X (ω) yields the self-diffusion coefficient, D X , divided by π for species X.The obtained values for Al (D Al = 1.2 × 10 −5 cm 2 /s) and O (D O = 1.3 × 10 −5 cm 2 /s) compare well with our earlier simulation results using the slope of the mean square displacements to determine D X .Experimental data from tracer diffusion measurements suggest a slightly larger diffusion coefficient (D = 2.65 × 10 −5 cm 2 /s at 2475 K) [29], which is still in reasonable agreement, considering that tracer and self-diffusion may not be exactly equivalent.

Collective dynamics and viscosity
The collective atomic dynamics is best represented by the double Fourier transform of the van Hove correlation function [30], which is the dynamic structure factor S(Q, ω).For a multicomponent system, partial dynamic structure factors, S XY (Q, ω), may be defined to represent correlations between different species (X and Y ).S XY (Q, ω) is obtained via the time correlation function of density fluctuations in reciprocal space δn X (Q, t), which yields the partial intermediate scattering functions (N X being the number of particles of species X) Its time Fourier transform is the partial dynamic structure factor The spectral quantity accessible experimentally in inelastic neutron or x-ray scattering experiments is a sum of partials weighted by the respective concentrations and coherent scattering lengths.Since Al 3+ and O 2− contain the same number of electrons, the x-ray scattering (weighting) factors are almost equal at low Q and the x-ray dynamic structure factor probes the fluctuations in the ion number density.Comparison of the total S(Q, ω) calculated from the MD results using the respective x-ray scattering factors for Al and O with the experimental S(Q, ω) at Q = 0.21 Å−1 from inelastic x-ray scattering [8] shows excellent agreement and confirms the quality of the AIM potential (see figure 2).A similarly good agreement was achieved in our earlier study [22] using the smaller 640 ion cell.In that case the smallest accessible Q was 0.31 Å−1 .Generally, the total (x-ray or neutron weighted) S(Q, ω) is difficult to interpret and the partials cannot easily be extracted from the experimental data.However, at small Q (usually substantially smaller than 1 Å−1 ) models of simple liquids, such as generalized hydrodynamic theory [1,2,31], have been successfully applied to experimental data of binary liquids, including alumina melt [8]. Figure 2 shows that in the low-Q region all S XY (Q, ω) look very similar, which justifies the treatment as a one-component system.
The Q-range in which side peaks or shoulders are visible in S(Q, ω) of liquids strongly depends on the type of interaction [3].While in simple ionic liquids, such as alkali halide melts they are barely visible even at the lowest Q accessible to IXS [7], they seem to be much more pronounced in more complex ionic melts as in the present case of alumina melt.It is interesting to note that shoulders are still clearly visible in S OO (Q, ω) at Q = 2.1 Å−1 , which is close to the principal peak position of the respective S OO (Q) at 2.6 Å−1 .The total S(Q, ω) in figure 2 suggests that these 'collective' modes cannot be observed by IXS in alumina melt due to the cancellation with the partial Al-O spectrum.However, designing an experiment with increased weighting factor of S OO (Q, ω) compared to the Al-O and Al-Al partials could make these shoulders visible.Perhaps this could be a challenge for inelastic neutron scattering.Since clear peaks in S XY (Q, ω) are only visible at small Q, we can use the maxima ω m of the spectra of the partial longitudinal current correlation functions 2) to study the dispersion relations of the collective excitations.As shown in figure 3, the dispersions from J OO l (Q, ω) and J AlAl l (Q, ω) are very similar in the Q-range up to at least 2.1 Å−1 .Alternatively, the dynamic structure factor is often fitted by generalized hydrodynamic models, which in the limit of low Q can be written as a sum of three Lorentzians [31]: A c and A s are the amplitudes of the central and inelastic lines, respectively.Γ c and Γ s are the corresponding linewidths and ω s is the excitation frequency.ω s (Q) is consistent with the two ω m (Q) dispersions at the lowest Q (see figure 3).However, the Lorentzian model does not describe well the spectra at Q larger than about 0.6 Å−1 .Assuming a linear dispersion at low Q, an estimate of the longitudinal sound velocity is obtained from the initial slope.The value of 6800 m/s as indicated in figure 3 is in excellent agreement with the IXS data [8].The damping term Γ s (Q) increases approximately with Q 2 , which will be discussed in the following.
In the hydrodynamic limit (low Q), the inelastic (Brillouin) peak broadening is expected to be quadratic in Q with the proportionality constant related to the longitudinal viscosity η l and to the thermal diffusivity D T [31] Γ s = 1 2 with γ = C p /C v , the ratio of the specific heats at constant pressure and constant volume, estimated as ∼ 1.08 [8].This dependence on Q is actually observed experimentally [8] and also in the present simulations (see figure 3), despite the relatively large Q values at which the observations are made.However, this should not be taken as conclusive evidence that the hydrodynamic picture of viscously damped propagating modes is appropriate -even disregarding the probably small thermal contribution (second term in equation ( 13)) the proportionality constant is much smaller than the measured macroscopic viscosity.Knowing the mass density ρ of the melt, the kinematic viscosity η l = 4/3η + η b may be obtained.η and η b are the shear and bulk viscosities, respectively.
From the simulation we have ρ = 2.83 g/cm 3 and proportionality constants η l /2ρ (figure 3) of 9.6 × 10 −7 m 2 /s, which gives a kinematic viscosities of 5.4 mPa•s.This is reasonably consistent with the analysis of the experimental data, where the damping constant yields a kinematic viscosity of about 3.5 mPa•s at T = 2323 K [8].These values are considerably (here about one order of magnitude) lower than the macroscopic viscosities [32].If the macroscopic viscosity is used to estimate the Brillouin linewidth, through equation ( 13) it would predict that the Brillouin lines were heavily overdamped.A representation of the data with three generalized Lorentzians is also suggested by a generalized hydrodynamic approach of Sinn et al. [8].In the Mori-Zwanzig formalism [33] the width parameter Γ s would be frequency dependent and related to the spectrum of a memory function.Since the information on the width parameter comes from fitting the spectrum in the vicinity of the Brillouin lines (i.e.close to ω s ), Γ s should be thought of as reflecting the liquid dynamics at frequencies close to ω s , which is itself proportional to Q.The memory function of interest in the present context should closely resemble the correlation function of the elements of the stress tensor [34] σ αβ (t)σ αβ (0) , where σ αβ is an off-diagonal element of the stress tensor of the simulation cell.For an isotropic system, an average over different elements α, β = x, y, z can be performed.The resulting time correlation function is shown in figure 4. The time dependence of σ αβ (t)σ αβ (0) suggests two different relation channels, which are modelled by fitting two stretched exponential functions of the form A exp(−(t/τ ) β ).The relaxation times τ differ by more than one order of magnitude.While the fast relaxation channel has a relaxation time of only 13 fs, the slower process decays in about 0.55 ps.
The spectrum of this function defines a frequency-dependent viscosity which is also shown in figure 4, together with the spectra of the two stretched exponentials to indicate how the two relaxation channels contribute to the spectral shape.The spectrum associated with the fast component is seen to resemble the vibrational density of states.It extends over a similar frequency range to the density of states shown in figure 1 and exhibits similar features at around 100 meV.The slow component dominates the low frequency spectrum and is associated with the structural (configurational) relaxation of the fluid.The zero frequency viscosity, η(0) = 25 mPa•s (marked by a full circle in figure 1), is the macroscopic viscosity, which is again in excellent agreement with experimental data [21,32].We see that its value will be dominated by this structural relaxation component.The temperature dependence of the structural relaxation time determines the temperature dependence of the viscosity.In a glass, the structural relaxation becomes indefinitely slow, and the spectrum η(ω) reduces to a vibrational spectrum very similar to the vibrational contribution to the liquid spectrum.
The Brillouin frequencies probed in the 0.2−1.0Å−1 Q-range will sample this spectrum between about 10 and 50 meV, where the height of the spectrum is much lower than at ω=0.We can therefore understand the order of magnitude difference between the "viscosities" extracted from the Brillouin linewidths and the macroscopic value.In this frequency range, η(ω) has a weak frequency dependence, which the spectral decomposition suggests should be associated with the tail of the slow structural relaxation component with a flat background associated with the vibrational term.We note that Sinn et al. [8] fitted their IXS data with a three Lorentzian model which allowed for a frequency-dependent viscosity and extracted a structural relaxation time of τ ≈ 0.5 ps at T = 2500 K, which agrees well with our calculated value.The fact that we still see a contribution from the structural relaxation at 2500 K could be the reason why our calculated Brillouin widths are larger than those reported experimentally, as noted above.The experimental values pertain to T =2323 K and it is possible that at this lower temperature, the structural relaxation spectrum has narrowed and no longer makes as large a contribution to the frequency-dependent viscosity in the range of Brillouin frequencies as at 2500 K.
This analysis suggests that in Al 2 O 3 at 2500 K the damping of the high-Q acoustic modes observed in IXS is crossing over from the hydrodynamic, viscous damping associated with structural relaxation of the fluid, to damping associated with the vibrational dynamics of the atoms about some disordered configuration in which the atoms are trapped on the timescale of the acoustic oscillations.As the temperature is raised, the damping due to the structural relaxation should increase resulting in a broadening of the Brillouin lines.At lower temperatures, it would be more appropriate to think of the fluid as responding like a disordered solid, and expect the Brillouin linewidth to be only weakly temperature dependent.In this régime the Brillouin linewidth reflects the relaxation by dephasing [35] of density fluctuations of wavevector Q. Schirmacher et al. [36] have recently shown how this mechanism can also lead to a Q 2 -dependence of the Brillouin linewidth.

Conclusions
A realistic modelling of alumina melt is achieved by using an advanced ionic interaction potential and molecular dynamics.Atomic structure and dynamics as well as transport coefficients are quantitatively reproduced.The modelling approach provides a deeper insight into the relation between atomic-scale and macroscopic properties and facilitates the interpretation of experimental data.Inelastic features are visible in S OO (Q, ω) up to relatively large Q close to the principal peak of S(Q).The visibility and broadening of Brillouin peaks seem to be closely related to the frequency dependent viscosity, which is in agreement with recent theoretical developments.

Figure 1 .
Figure 1.Velocity autocorrelation functions and partial frequency spectra z O (ω) and z Al (ω) at T = 2500 K.The zero frequency limit yields the corresponding self-diffusion coefficients D X /π.

Figure 2 .
Figure 2.Top: Constant Q projections of the partial and total x-ray dynamic structure factor at T = 2500 K.For comparison with experimental data from Sinn et al.[8] (symbols) the calculated spectra are convoluted with a Gaussian resolution function with a width of 1.8 meV FWHM.Bottom: The same dynamic structure factors but not convoluted and multiplied by ω 2 , which is proportional to the partial and total spectra of the longitudinal current correlation functions, respectively.

Figure 3 .
Figure 3. Dispersion relations ω(Q) (full symbols) and line widths Γ(Q) (open symbols) obtained from the maxima of the partial longitudinal current correlation spectra (circles) and from fits of the total x-ray S(Q, ω) with the Lorentzian model (triangles).The initial slope of the dispersion relations assuming linear behaviour corresponds to a sound velocity of 6800 m/s.The Q 2dependence of the linewidth Γ at small Q is indicated.

Figure 4 .
Figure 4. Left: average stress tensor autocorrelation function (solid line) fitted with two stretched exponentials (dashed line).Right: corresponding spectrum (frequency dependent viscosity).The zero frequency value of η(ω) = 25 mPa.s, which corresponds to the macroscopic shear viscosity, is indicated by a full circle.At 10 meV, η(ω) has dropped by more than a factor of ten.Also shown are the spectra of the fitted curve and its fast and slow components.