Mod:MgOjh6215

Thermal Expansion of MgO

Written by James Hale.

Introduction

The Structure of MgO

Figure 2 - 2x1x2 Primitive cell of MgO

Figure 1 - Conventional Cell of MgO

The MgO lattice is a face centred cubic (fcc) ionic lattice, which is represented in three different ways during this experiment. The conventional cell (figure 1), the primitive cell (figure 2) and the super cell. The conventional cell contains 4 MgO units and has the lattice structure: a=b=c, α=β=γ=90°. It has the same density but 4 times the volume of the primitive cell, which only contains 1 MgO unit and has the lattice structure: a=b=c, α=β=γ=60°. Even though the primitive cell contains fewer MgO units, the conventional cell is more widely used due to the simplicity of its symmetry. The supercell is used for the molecular dynamics calculations and is the conventional cell doubled in all 3 dimensions to make a 2x2x2 cell. This has 8 times the number of atoms and the lattice parameters defining the cell are twice as large (2a_Conventional = a_Supercell).

It is important to note that the computer will output the results for a primitive cell on completing calculations using the quasi-harmonic approximation. Therefore, the results will have to be converted to those of a conventional for direct comparison with other data.

Quasi-harmonic Approximation

The quasi-harmonic approximation is an improvement on the harmonic approximation for modelling lattice dynamics. For the harmonic approximation the atoms in the lattice only interact with their nearest neighbours. The energy of a 1 dimensional monatomic chain can be modelled by the following Taylor expansion of the interactions between the nearest neighbours:

$E = N\phi + \sum_{s\geq 1} \frac{1}{s!} \frac{\partial^s \phi}{\partial u^s} \sum_n (u_n - u_{n+1})^s ~~ Equation~1$

The equation states that the energy of a 1 dimensional monatomic chain, $E$, is equal to: the rest energy of every atom; where $N$ is the number of atoms and $\phi$ is the energy between 2 neighbours as a function of the distance between them, and the energy of the vibrations within the monatomic chain. The vibrational energy is equal to the Taylor series for the displacement of each atom from its rest position, where: $s$ is the order of the Taylor series, $u_n$ is the displacement of the nth atom and $u_{n+1}$ is the displacement of the nth +1 atom. This represents the correction of the energy between 2 neighbours from the rest energy where the distance between two neighbours is equal to the lattice constant, $a$ . For this equation to be formed the displacement of each atom, $u$ , must be small in comparison to the lattice constant, $a$ otherwise the series may not converge. Usually $u$ is much smaller than $a$, and so a quadratic approximation can be applied ($s=2$). This approximation can be made not just because the higher order terms are small and can be ignored but because harmonic equations have exact solutions, whereas anharmonic equations give numerical solutions. The additional level of accuracy is not required in this investigation and so using the anharmonic method would lead to an unnecessarily high computational cost. Therefore, the harmonic energy of a 1 dimensional monatomic chain is given by equation 2:

$E^{harmonic} = \frac{1}{2} \frac{\partial^2 \phi}{\partial u^2} \sum_n (u_n - u_{n+1})^2 ~~ Equation~2$

To solve the problem of having open ends of the monatomic chain, the Born-von Karman periodic boundary condition is imposed. Harmonic motion along the continuous chain is sinusoidal and the motion of the entire system can be modelled by travelling waves. The time dependent motion of the nth atom is dictated by the linear superposition of the traveling waves and is shown by equation 3:

$u_n(t) = \sum_k u_k exp(i[kx-\omega_kt]) ~~ Equation~3$

Where $u_n$ is the displacement of the nth atom, $t$ is time, $k$ is the wave vector and is equivalent to $\frac{2\pi}{Wavelength}$, $\omega_k$ is the angular frequency and $u_k$ is the amplitude. $x$ is the displacement of the wave can only be equal to $na$ where $n$ is an integer. The wave equation is for each wave can be solved if the wave vector corresponding to that wave is considered individually. This can be simplified to equation 4:

$\omega_k = \left(\frac{4J}{m}\right)^{1/2} |sin(ka/2)| ~~ Equation~4~~~~$

Where $J$ is the force constant $J = \frac{\partial^2 \phi}{\partial u^2}$ and m is the atomic mass. This is the equation which is used to the plot dispersion curve, which is the frequencies of phonon vibration at each k-point (figure 3). However, reciprocal space, which is equivalent the Fourier transform of the original lattice, must be used to plot the dispersion curve. The relation between the lattice constant in real and reciprocal space is:

$a* = \frac{2\pi}{a} ~~ Equation ~5$

Using reciprocal space allows for the wave vectors in higher dimensional lattices to be calculated. The density of states (table 1) at each k-point on the dispersion curve can be found and summing over all of the k-points, all of the vibrational contributions to the free energy can be found. Combined with the quasi-harmonic theory, reciprocal space can therefore facilitate the calculations of force constants and energy between atoms in a three dimensional lattice, ultimately leading to the calculation of the Helmholtz free energy of the lattice in equation 6 below:

$A = E_0 + \frac{1}{2}\sum_{k,j}\hbar\omega_{k,j}+ k_bT\sum_{k,j}ln[1-exp(-h\omega_{k,j}/ k_bT)]~~ Equation~6$

Where $A$, is the Helmholtz free energy, $E_0$is the internal energy of the lattice, $\frac{1}{2}\sum_{k,j}\hbar\omega_{k,j}$ is the vibrational ground state energy of the lattice and the final term is the contribution from higher vibrational states. The quasi-harmonic approximation is able to calculate the free energy of the lattice accurately because it is able to consider the contributions from the internal energy, which is trying to make the crystal smaller and the entropy, which is trying to increase the distance between the atoms to reduce the frequency of phonon vibrations. The quasi-harmonic approximation makes the assumption that all of the anharmonicity is contained within the resulting thermal expansion and therefore can assign volume dependent harmonic approximations to each lattice volume. This works well at lower temperatures because the higher order terms from equation 1 remain small. However, at high temperatures the anharmonic effects become larger and the higher order terms begin to make the quasi-harmonic approximation diverge from reality.

For diatomic lattices, the dispersion plots contain two types of branches. The first is the acoustic branch, where the different neighbouring atoms vibrate in phase and the second is the optic branch, where the different neighbouring atoms vibrate out of phase. The optic vibrations are of higher energy and are close to the optical part of the electromagnetic spectrum.

Thermal Expansion

Thermal expansion is an entropy driven process and is the tendency of a substance to increase in volume with increasing temperature. It is modelled by equation 7 below:

$\alpha_V = \frac{1}{V_0} \left(\frac{\partial V}{\partial T}\right)_P ~~ Equation~ 7$

Where $\alpha_V$ is the coefficient of thermal expansion, $V_0$ is the initial volume and $\left(\frac{\partial V}{\partial T}\right)_P$ is the partial derivative of volume with respect to temperature at constant pressure. This equation can be evaluated graphically with the gradient of a temperature against volume graph equalling the expansion coefficient after it is correct by the initial volume term. Thermal expansion occurs because a larger volume reduces the frequency and thus the energy of vibrations. The vibrations increase in frequency at higher temperatures and the volume of the substance increases to account for their contribution.

Thermal expansion is very important when considering the construction of buildings and machines. The thermal expansion may cause a structure to collapse or be needed to allow a mechanism to work more efficiently, like in combustion engines.

Experimental

All parts of the investigation were run on a linux operating system, using a DLV interface to observe the structure and properties of the MgO lattice. The lattice structures to be investigated were downloaded from 2 premade files: one conventional cell and one supercell. All calculations were run using GULP. To prevent excessive computational cost the shells were modelled as hard spheres of charge rather than shells of charge, which can be polarised and so Include Shells option was not selected.

Calculating the Phonon Modes of MgO

The conventional cell, which can be simplified to the primitive cell, was used for this part of the investigation because calculations run with the conventional cell require much less computational power than with the supercell. The phonon dispersion curve was calculated first. The phonon energies were found at 50 points along the k-space between the k-points W-L-Γ-X-W-K. To achieve this the phonon dispersion was computed, with phonon eigenvectors calculated. The temperature was set to 300 K and the Npoints were set to 50. Some phonons were then visualized to check a successful calculation had been completed.

The phonon density of states was then found for k-point grid sizes ranging from 1x1x1 to 64x64x64 at 300 K. The grid sizes were adjusted using shrinking factors A,B and C. For example, when A=1, B=1 and C=1 a k-point grid is formed with dimensions 1x1x1.

Computing the Free Energy - The Quasi-harmonic Approximation

The log files produced during the density of states calculations were examined to find the free energy given by each grid size. The larger the grid size the more accurate the free energy calculation. The ideal grid size is a balance between computational cost and the appropriate level of accuracy. The converged value for the free energy and the ideal lattice size were found.

The Thermal Expansion of MgO

The thermal expansion of MgO was investigated using the quasi-harmonic approximation with the conventional cell and then the molecular dynamics method with the supercell.

For the quasi-harmonic approximation a 24x24x24 grid size was found to be the ideal grid size to provide the best balance between accuracy and computational cost. The density of dates calculation was still active in order to allow for the larger grid size. The conventional cell was used to maintain a low computational cost. The Gibbs energy of the MgO lattice was then optimised over the temperature range 0 K to 2900 K in 100 K steps. Even though the melting point of MgO is 3125 K, the quasi-harmonic approximation would give an error for temperature at 3000 K or above.

For the molecular dynamics method the supercell was used. Even though it requires greater computational cost, the molecular dynamics approximation would not work with a conventional cell because the every cell in the lattice would be moving in phase. This is not a valid representation of the MgO structure and therefore a larger starting cell must be used. The molecular dynamics calculations were run with pressure, the number of particles and temperature all fixed. The time step between each recorded position of the individual atoms was set to 1 fs and the number of equilibration steps, the number of steps allowed for the system to reach equilibrium, was set to 500. The production steps, the number of steps run after the system has reached equilibrium, was also set to 500. The number of sampling steps, the number of steps over which the means are taken, and trajectory write steps, the number of steps between which the geometry is saved in order to produce animations, were both set to 5. The temperature range for the calculations was 100 K - 3300 K with 100 K steps. A 50 K reading was also taken because the simulation did not work at 0 K. The vibrations were animated after each calculation to check that it had run correctly.

Results and Discussion

Calculating the Phonon Modes of MgO

Figure 3 - The phonon dispersion curve of MgO

Table 1

Grid Size	Phonon Density of States Plot
1x1x1
3x3x3
24x24x24
32x32x32
64x64x64

Computing the Free Energy - The Quasi-harmonic Approximation

Figure 4

Table 2

Grid Size	Free Energy / eV	Difference Between Measured and Converged Values
1x1x1	-40.930301	0.003818
2x2x2	-40.926609	0.000126
3x3x3	-40.926432	0.000051
4x4x4	-40.92645	0.000033
5x5x5	-40.926463	0.000020
6x6x6	-40.926471	0.000012
8x8x8	-40.926478	0.000005
10x10x10	-40.92648	0.000003
12x12x12	-40.926481	0.000002
16x16x16	-40.926482	0.000001
20x20x20	-40.926483	0.000000
24x24x24	-40.926483	Converged
32x32x32	-40.926483	Converged
48x48x48	-40.926483	Converged
64x64x64	-40.926483	Converged

Figure 5

The Thermal Expansion of MgO

Table 3

Temperature / K	300	800	1300	1800
Thermal Expansion Coefficient for the Quasi-harmonic Approximation x10-6 / K-1	26.7	32.7	38.5	47.9
Thermal Expansion Coefficient for the Molecular Dynamics Method x10-6 / K-1	30.1	32.6	23.8	39.7
Literature Values for the Thermal Expansion Coefficient x10-6 / K-1	31.2	42.6	47.1	51.3

Figure 6

Figure 7

Figure 8

Figure 9