Rep:Ga714MgO

Introduction

Fig 1. Thermal expansion of the tracks on a railway

Thermal expansion in an extremely important property to consider in a variety of fields, but few more important than in engineering. Depending on a material's thermal expansion coefficient ($\alpha$), it may expand more or less than another material (see Eqn 1). Selecting an inappropriate material for a particular application or environment can have catastrophic consequences. For example, railway tracks can bend to the extent that trains would be derailed (Fig 1). Furthermore, choosing materials of differing thermal expansion coefficients for the construction of a bridge can result in extremely high stress at the interface between the two materials as one has expanded more than the other. The result of this is at best distortion of the structure but can also lead to fracturing and failure of the design. It is clear that understanding the thermal expansion of materials is critical to the design of some of the world's greatest projects.

This experiment involved simulating the thermal expansion of the face centred cubic MgO lattice. For the purpose of this experiment, it is assumed that the crystal is perfect with no defects and infinitely large. Initially, the phonon dispersion curves and the phonon density of states are simulated for a variety of grid sizes in order to select an optimal grid size on which the following calculations are based. This is a compromise between accuracy of the result and time as with increased grid size more accurate results are obtained but the calculation is more time expensive. Subsequently, by setting the temperature, the Helmholtz free energy and volume of the lattice under investigation were calculated using the Quasi-Harmonic Approximation and Molecular Dynamics. The two different calculations were then compared and the thermal expansion coefficient was calculated

Thermal Expansion Coefficient

The thermal expansion coefficient is an intrinsic property of a material. It is the change in volume of a material for an incremental change in temperature and is normalised by dividing by the initial volume.

${\alpha }_v=\frac{1}{V_0}{\left(\frac{\delta V}{\delta T}\right)}_P$ Eqn. 1

Experimental

Initially, the structure of MgO was optimised using the at 0K and the structure and properties of the MgO lattice of the output file was analysed. Subsequently, the phonon dispersion curves were calculated. Then an optimal grid size was found by calculating the phonon density of states using various shrinking factors along the A, B and C directions. Using the Quasi-Harmonic Approximation, the free energy of the lattice was calculated at 300K for a range of grid sizes, following the convergence of the Free Energy value. The thermal expansion of the MgO lattice was simulated using the same method with a 19x19x19 shrinking factor for the temperatures in the range of 0-1000K at 100K increments. Finally, the thermal expansion of MgO was investigated using Molecular Dynamics simulations. These were performed for the temperatures in the range of 100-1000K for every hundred Kelvin. Finally, both simulations were run at 1500, 2000 and 2852 K in order to understand what happens as the melting point of the material is approached.

Computational Methods

All the simulations were run in DL Visualize, which is a graphical user interface that allowed for the visualisation of the crystal lattice in question. The calculations were performed using the General Utility Lattice Program (GULP), which is a general purpose code which takes atoms into account, rather than electrons and nuclei. The atoms or ions are prevented from collapsing by a repulsive wall which is placed between the ions. It is possible to run both Quasi-Harmonic Approximations and Molecular Dynamics simulations in this program.

Fig 2. The harmonic oscillator superimposed on the Morse potential

Quasi-Harmonic Approximation (QHA)

The Quasi-Harmonic Approximation is an improvement on a purely harmonic approximation as it allows for the volume-dependence of the harmonic phonon. A purely harmonic phonon model would imply that no expansion occurs as the interatomic distance would be temperature independent. Consequently, the Quasi-Harmonic Approximation can be viewed as a series of independent harmonic phonons that stands true for every lattice parameter. Despite the approximation being based on the classical mechanics of the harmonic oscillator, the Quasi-Harmonic Approximation is a quantum mechanical problem made up of discrete vibrational energy levels. The large number of phonon energy levels contributes to the band structure of the energy and so the eigenstates can be described by a Bloch function where the wave function and thus energy varies continuously with reciprocal lattice vector k. As a result, quantum mechanical phenomena are observed, for example the zero-point energy of the system, as the lowest energy vibration (a translation) is not zero. The Quasi-Harmonic Approximation does not account for an anharmonicity as can be seen in Fig. 2. In order to calculate the Helmholtz free energy, a temperature is input and the harmonic vibrational frequencies are calculated. This allows for the entropy of vibration to be calculated for which a function for the Helmholtz free energy is obtained. The minimum value of this function is obtained and that is the Helmholtz free energy of the structure at the given temperature.

Molecular Dynamics (MD)

Molceular Dynamics is a classical approach to the same problem that is being solved by the Quasi-Harmonic Approximation, it is based on the Newtonian mechanics of forces and his equations of motion. A Molecular Dynamics simulation run at a specific temperature randomly assigns atomic velocities periodically from a Maxwell-Boltzmann distribution. From here the atoms are allowed to interact for this periodic time and the dynamics and evolution of the system are observed. This process is repeated many times in order to obtain an average. For systems that obey the Ergodic hypothesis, the time average properties allow for the macroscopic thermodynamic properties to be calculated and in this way the free energy of the system and thus the volume and thermal expansion coefficient can be calculated.

Results and Discussion

MgO Structure

MgO has a face-centred cubic lattice structure. This is most easily visualised in a conventional cell with ions on each vertex and centre of face of the cube, giving 4 Mg and 4 O atoms per cell. This, however, is not the simplest way to describe the lattice. Therefore, the output of the optimisation is the primitive cell with lattice parameters 2.9783 Å and 60.0000°. This is a rhombohedral structure with atoms on each vertex and one at the centre, giving 1 Mg and 1 O atom per cell. Four of these primitive cells exist in the conventional cell and one primitive cell can be seen stretching from opposite vertices inside the conventional cell. So the simulated conventional cell length can be calculated to give 4.212 Å. This compares very well to the literature value of 4.211 Å and it can be concluded that this initial simulation provides a suitable approximation for the MgO lattice.

Phonon Dispersion Curves

Fig 3. Simulated Phonon Dispersion Curve for MgO Lattice

Phonon dispersion curves arise due to the variations in vibrational frequency with the reciprocal lattice vector. The graph is made up of branches and the number of branches in the curve depends on the number of atoms in the primitive cell and the number of dimensions. These branches can be divided into two types: acoustic modes, lower energy branches arising from the in-phase vibration of neighbouring atoms and optical, higher energy branches due to the out-of-phase vibration of neighbouring atoms. At the centre of the Brillouin zone ($\mathrm{\Gamma }$), the energy of acoustic vibrations tend to zero and optical branches are non-zero. Instead, at $\mathrm{\Gamma }$ the out-of-phase vibrations result in a time-dependent electric dipole moment. In 3D,there are 3 acoustic branches and 3N-3 optical branches. Therefore, for the MgO crystal, there are 3 optical branches. Branches may merge at specific points in the Brillouin zone due to symmetry. The total number of branches is always the number of dimensions multiplied by the number of atoms in the primitive cell. Therefore, the equation below renders itself clear using the aforementioned information, where D is the number of dimensions and N the number of atoms in the primitive cell.

$TotalBranches=D(1_{AcousticContribution}+(N-1{)}_{OpticalContribution})$ Eqn. 2

Fig 3 shows the simulated phonon dispersion curve for the MgO lattice. As expected, 6 branches can be seen at certain points in the Brillouin zone. At L we can see that 2 of the acoustic branches merge and so too another acoustic and one optical. Then as we move along the Brillouin zone towards the centre the 3 acoustic branches tend to zero, the two merged acoustic branches continue merged and the other merged branch splits as the acoustic branch tends to zero and the optical remains non-zero.

Phonon Density of States (PDOS)

The PDOS is a plot of the phonon frequency vs intensity. Thus, it shows the number of vibrational states for a given frequency. The integral of the entire plot is 1 showing that entire system is described by the vibrational states enclosed in the integral. The density of states plays a key role in determining the physical properties of a solid. However, it is time expensive for the simulation to cover every point in the Brillouin zone and give the complete PDOS. Therefore, a compromise must be made between the accuracy and time taken to run the calculation. The calculation is based on the grid size or shrinking factor. With increasing grid size more values of the reciprocal lattice vector (k) are covered, and so a more accurate representation of the PDOS is obtained.

Grid Size - Accuracy Vs Time Compromise

Fig 4. The DOS calculation for the grid size 1x1x1

The first grid size selected was the smallest and therefore the quickest to run. With a 1x1x1 grid size the DOS for a single value of k was returned. It can be seen that there are two large peaks of intensity 0.35 at approximately 300 and 350 cm^-1 and two peaks of intensity 0.15 at around 700 and 800 cm^-1. Comparing the peaks in Fig 4 to the phonon dispersion curve in Fig 3, it is clear that the k-point at which these phonons are observed are at symmetry point L (0.5, 0.5, 0.5) and this is the primary symmetry point used by the code and not the origin $\mathrm{\Gamma }$. This observation also explains the difference in intensity between the lower energy and higher energy peaks. On the phonon dispersion plot, it can be seen that two branches feed into each of the lower energy frequencies at k-point L. The resulting phonons are doubly degenerate and therefore it is expected that the DOS is greater for these frequencies than for those with only one branch.

Fig 4. The DOS calculations for (starting top left and going clockwise) 5x5x5, 20x20x20, 30x30x30, 40x40x40

Obviously, one value of k does not give a fair representation for the DOS across the crystal lattice. Subsequently, the changes in the shape of the DOS graph changes with increasing grid size (see Fig 5). With increasing grid size the time for the calculation to run increased significantly. However, it can be seen that a true understanding of the density of states was obtained for larger grid sizes as the curves get smoother and the sharp peaks disappear. From 20x20x20 upwards, the shape of the curve remains roughly the same, as does the range. There are small changes in intensity for a few frequencies, notably around 300 cm^-1 but not significant enough to warrant the extra processing time. Therefore, a suitable grid size to give a reasonable approximation for the PDOS in this experiment was 20x20x20.

The changes in the Helmholtz Free Energy (A) with increasing grid size were calculated using the Quasi-Harmonic Approximation for the MgO lattice at 300k. Below is a table of the data obtained. It can be seen that the Helmholtz Free Energy rapidly converges to -40.926483 eV. This data concurs with the previous simulations in the sense that a 20x20x20 grid size is suitable for the calculations on the MgO lattice. The 20x20x20 grid size was selected over larger grid sizes as it is less time expensive to run simulations for this grid size. It appears that the free energy converges significantly faster than the PDOS graph.

Grid Size	A (eV)	ΔA (eV)
1x1x1	-40.930301	0.003818
2x2x2	-40.926609	0.000126
4x4x4	-40.926450	-0.000033
8x8x8	-40.926478	-0.000005
16x16x16	-40.926482	-0.000001
17x17x17	-40.926482	-0.000001
19x19x19	-40.926483	0.000000
20x20x20	-40.926483	N/A

This chosen 20x20x20 grid size is only suitable for the MgO lattice being considered. Changing the lattice parameters changes the smallest grid size required to give reasonable results. Consequently, the lattices of calcium oxide, faujasite and lithium will be compared to that for MgO.

$a^{*}=\frac{2\pi }{a}$ Eqn 3

Given that Ca and Mg are in the same group it can be assumed that CaO also adopts a face-centred cubic lattice structure. The lattice parameter for calcium oxide is slight larger than that of MgO (4.800 Å). This is due to the Ca²⁺ ion being larger than the magnesium cation. Eqn 3 shows us the relationship between reciprocal space lattice parameter a^* and a is multiplicatively inverse. This larger lattice parameter in real space results in a smaller Brillouin zone in reciprocal space. Thus, it would theoretically be ok to use a smaller grid size in order to obtain accurate results. However, the difference in size of the lattice parameters of MgO and CaO is small and so little time would be saved in making the grid size marginally smaller. Therefore, using the same 20x20x20 grid size would be recommendable.

Faujasite also has a cubic structure and, despite showing more covalent character in its bonding, it will be assumed that it is directly comparable to MgO. The lattice parameter for faujasite is much larger than that of MgO and CaO. Therefore, a significantly smaller grid size would suffice to render a smooth PDOS graph, and thus be appropriate for running calculations.

Lithium, on the other hand, demonstrates metallic bonding. This means that the delocalised sea of electrons would shift to mitigate repulsive interactions between lithium cations when they are displaced. This means the phonon energies are lower and there is a smaller phonon dispersion. Consequently, the PDOS is over a smaller range and so a smaller grid size is required to cover the less disperse range.

QHA Thermal Expansion Calculations MgO

Free Energy Vs Temperature

Fig 5. Plot of Helmholtz Free Energy Vs Temperature for the QHA calculation

Fig 5 shows that using the Quasi-Harmonic Approximation, the Helmholtz Free Energy of the lattice gets more negative with increasing temperature. The gradient of the graph also gets more negative until an approximately linear region is observed from 1000K upwards. These phenomena can be explained by Eqn 3 below.

$A=U-TS$ Eqn. 3

It can be seen that for a material with an exothermic lattice energy and given that absolute entropy is always positive the Helmholtz Free Energy will always be negative. At 0K, the zero-point energy of the system is observed, a quantum mechanical phenomenon that indicates there is residual free energy even at 0K. At low temperatures U dominates the value of the free energy. The non-linear region is the part where both the entropic and internal energy components have significant contributions to the free energy. The linear regions occurs where the entropic component dominates the value of the free energy as the temperature is high. The changes in U and S with increasing temperature are insignificant in comparison so A becomes approximately proportional to T.

Lattice Parameter and Cell Volume Calculations

With increasing temperature, more energy is available to the system. This results in higher energy phonons being available for population. These phonons tend to have larger displacements as neighbouring atoms will be moving out-of-phase with respect to eachother. This is essentially increasing the kinetic energy of the molecules so that they vibrate more and thus increase the average internuclear separation. The result is thermal expansion which is most common consequence of increasing temperature. However, although rare, it is possible that a material can contract upon heating and this phenomenon is known as negative thermal expansion. Fig 6 demonstrates the increase in the lattice parameter with increasing temperature for the QHA. Therefore, MgO exhibits normal thermal expansion with increasing temperature. If the calculation were based on purely harmonic vibrations then there would be no change in the lattice parameter. This is because despite increasing vibrational energy level the equilibrium internuclear separation will never change. See Fig 2 to see that the mid point of each vibrational energy level is the equilibrium position. The QHA corrects the Harmonic Approximation at each vibrational energy level to give a more Morse-like potential curve. This means with increasing vibrational energy level the equilibrium position for the internuclear separation increases (see Fig 2). This explains the trend we see in Fig 6 for the plot of Lattice Parameter vs Temperature. However, it is important to note that using the QHA the Harmonic component determines that no dissociation will occur.

Fig 6. Plot of Lattice Parameter Vs Temperature for the QHA calculation

Fig 7. Plot of Lattice Parameter Vs Temperature for both QHA and MD calculations

As the lattice parameter has increased, the volume of the primitive cell has increased. In Fig 7, the increase in volume with temperature is plotted. This plot allows us to calculate the thermal expansion coefficient of MgO for different temperature ranges and the results are discussed below.

Molecular Dynamics Calculations

Molecular Dynamics simulations are initiated by giving an or some atom(s) velocities and allowing the system to equilibrate. A primitive cell is too small to perform this sort of calculation as it would not allow for many vibrational modes and every primitive cell in the large lattice would be moving in phase. Therefore, a larger cell is needed to emulate the classical mechanics that would occur in an infinite lattice. However, increasing the cell size makes the calculation more time expensive. So a compromise must be reached between accuracy of results and the run-time. It was decided that a 64 atom supercell would yield sufficiently accurate results. The 64 atom super cell has 32 times the volume of the primitive cell and, hence, contains 32 MgO units.

Molecular Dynamics Cell Volume Vs Temperature

Fig 7 shows the effect of temperature on the cell volume simulated using Molecular Dynamics. This simulation agrees with the QHA simulation in that it predicts that MgO exhibits normal thermal expansion. Furthermore, in the temperature range of 500-1000 K, the values are fairly similar. However, outside of this region the two curves diverge. The reasons for this observation are discussed with the thermal expansion coefficient below. However, it is clear that, as the Molecular Dynamics simulation is based on classical mechanics, the cell volume decreases linearly as the temperature approaches zero. This is because the quantum mechanical zero-point energy is not account for. This is also the reason that no data was collected at 0K because, according to classical mechanics, the system would have no kinetic energy, and so the Molecular Dynamics simulation fails. Above 1000K the curves diverge as the QHA predicts increasing infinite expansion and Molecular Dynamics, as it only relies on the mechanics of motion, does not increase as significantly as the melting point temperature is approached.

Molecular Dynamics Internal Energy Vs Temperature

Fig 8. Plot of Internal Energy Vs Temperature for MD simulation

The most striking feature of Fig 8 is that the gradient is the opposite sign to that of Helmholtz Free Energy Vs Temperature as calculated by the QHA. However, both are true, whilst the QHA computed the Helmholtz Free Energy, which gets more negative with increasing temperature, Molecular Dynamics computes the internal energy (U), as it does not calculate the entropy of the system.

$\mathrm{\Delta }U=\mathrm{\Delta }q+\mathrm{\Delta }w$ Eqn. 4

It has already been seen that with increasing temperature the cell volume increases. This expansion is work done by the lattice and therefore this makes the internal energy more positive, see Eqn 4 (+q is heat transfer out the system and +w is work done by the system). This explains the positive gradient shown in Fig 8.

Molecular Dynamics at the Melting Point

The Molecular Dynamics calculation begins to deteriorate as the experiment approaches the melting point. The standard deviation of the instantaneous Cell Volumes at 2852 K is 12.97 K, whereas at 1000 K the standard deviation is 4.11 K. This indicates that as the simulation is run closer to the melting point the simple Molecular Dynamics simulation begins to break down. But it can be seen below, that despite the increase in the standard deviation of the instantaneous cell volumes, the MD simulation still yields accurate results at high temperatures approaching the melting point.

Thermal Expansion Coefficient Results

Temperature Range (K)	Experimental $\alpha$ * 10-6(K-1) (QHA)	Experimental $\alpha$ * 10-6(K-1) (MD)	Literature $\alpha$ * 10-6(K-1)
0-300	9.47	27.9	10.4
300-1000	27.9	28.8	31.1
1500-2000	41.6	36.2	37.8

Upon comparing the calculated thermal expansion coefficients with the literature values, in the temperature range of 0-300 K, it is clear that the QHA is the far more suitable approximation. This is due to the fact that the QHA accounts for the zero-point energy of a material and so at low temperature, where thermal motion has little effect on the expansion of the lattice. The MD result for this temperature range is far from accurate as it ignores this dominating factor at low temperatures. However, with increasing temperature it is clear that the MD calculations become more accurate and also more appropriate than the QHA calculations. At higher temperatures, classical mechanics dominates the thermal expansion property of a material and so the zero-point energy has little influence. The QHA becomes unsuitable as it does not correctly simulate that anharmonicity at higher temperatures and predicts infinite expansion without dissociation. Therefore, at lower temperatures it would be recommendable to use a QHA to determine the thermal expansion coefficient of a material and for higher temperatures the MD simulation is more appropriate.

Conclusion

It is clear that each simulation method has its advantages and limitations. The QHA is far more accurate at lower temperatures and thus is more suitable for studies into designs of projects that will be exposed to every-day weather. However, if a material like a ceramic were to be used at high temperatures, a Molecular Dynamics simulation would be more suitable for modelling its thermal expansion properties. Furthermore, both simulations were extremely quick and facile to run. So, in comparison to measuring the properties using neutron-scattering, they are fantastic tools that give good approximations for properties that would usually take significantly longer to obtain values. Therefore, more research should be placed into refining the simulations in order to allow for the QHA to better factor in the anharmonicity and, potentially a correction to be made to the MD simulation to account for the zero-point energy.

It would also be interesting to not how these simulations fare upon encountering a material that exhibits negative thermal expansion. The ideas upon which this calculations are based may not account for this phenomenon, and so a different simulation may be required.