Rep:Mod:js6313 mgo

Introduction

At any temperature above 0K, atoms within a crystal vibrate; the modes of these vibrations are called phonons. These can be modelled in a number of ways, including the quasi-harmonic approximation (QMA), where each bond is equivalent to a harmonic oscillator, and equilibrium bond lengths are temperature dependent. A more accurate model hostilities molecular dynamics (MD), which is based upon allowing the dynamic development of forces in agreement with Newton's second law. In order to ascertain a completely representative view of all physical properties of a lattice, the sum over all possible values of k would have to be taken - this would take infinitely long, so a compromise is made between accuracy and computation cost. For the QMA, this was found to be a grid size of 20.

When any crystal is heated, it expands. The extent to which this happens is encapsulated by the thermal expansion coefficient, $α_{V}$ , which relates how much the crystal expands when subjected to a change in temperature $δ V$ :

$α_{V} = \frac{1}{V} {(\frac{δ V}{δ T})}_{p}$

Magnesium oxide is a basic cubic lattice at temperatures below 2800 °C. Its high melting point can be attributed to the high lattice enthalpy that results from small ions with high charges. The unit cell is a parallelepiped with

$α = β = γ = 6 0^{\circ}$ and $a = b = c = 2.9783 Å$ which gives a conventional lattice constant of $2.10597 Å$ .

In general, both the QMA and MD models are in good agreement with experiment, although the thermal expansion of MgO at very low temperatures deviates from literature when using the former.

Overview of The First Brillouin Zone, Dispersion Curves and the Density of States (DOS)

In a crystal with lattice vector a, the wavevector $k = \frac{2 π}{λ}$ . Only wavevectors between $- \frac{π}{a} and \frac{π}{a}$ are unique; these are in what is known as the first Brillouin zone. Every k-point in this domain is unique, and has it's own energy level which under ideality is given by:

$E (k) = \frac{ℏ^{2} k^{2}}{2 m} + E_{p o t}$

Note that this only holds true within the stated boundary conditions for the first Brillouin zone.

Dispersion curves are a plot of frequency against k, within the first Brillouin zone. Each line is a phonon - at any point, a maximum of six phonons is observed. This is because there are 2 atoms in each unit cell, each with 3 vibrational degrees of freedom. At some points (where lines converge), phonons are degenerate.

**Fig 3.** Density of phonon states along k-point L (left) and along 50 k-points (right)

A density of states (DOS) graph shows the frequency of phonons along a given number of k-points; setting this to 1 gives a view along the L k-point shown below. Note that the convergence of the lower energy curves corresponds to a higher DOS. As the number of sampled k-points increases, a more representative view of the density of states is gained.

Computing the Free Energy - The Quasi-Harmonic Approximation vs Molecular Dynamics

For the quasi-harmonic approximation, Every MgO bond can be thought of as a mass on a spring, the motion of which is governed by Hooke's law:

$F = - k (x - \bar{x})$

This is a good approximation at very low temperatures, where bond lengths do not change significantly, however it breaks down once vibrations start to stretch bonds closer to their breaking points. The QMA is a variation on this model whereby k is dependent on the volume, and that the and that the harmonic approximation works for all volumes.

Molecular dynamics takes a set of starting conditions, and takes a series of snapshots at regular intervals (5 ps in this instance), until the properties reach equilibrium. A random selection of velocities is selected such that the overall temperature of the 32 unit cells is that specified beforehand; from here, the forces are computed, and then:

$a = \frac{F}{m}$

$v_{new} = v_{old} + a δ t$

$r_{new} = r_{old} + v_{new} δ t$

where r is the position, v is the velocity and a is the acceleration.

In order to ascertain the minimum shrinking values required to get an accurate (converged) value for the free energy, jobs were using values of 1, 2, 3, 4, 5, 10 and 20. All calculations here explicitly set the temperature of 300 K and the pressure to 0 GPa. Note how the computation time increases exponentially with shrinking value:

Grid size (n)	Number of Brillouin zone points sampled ( $\frac{n^{3}}{2}$ for n > 5)	Helmholtz free energy (eV)	Difference between FE at grid size = n and FE at grid size = 50 (meV)	Calculation time (ms)
1	1	-40.930301	3.818	2.0
2	4	-40.926609	0.126	2.8
3	18	-40.926432	0.051	6.3
4	32	-40.926450	0.033	9.1
5	75	-40.926463	0.020	20.9
10	500	-40.926480	0.003	136.7
20	4000	-40.926483	0	1479.8
50	62500	-40.926483	0	107751.6

From this, it was decided that a grid size of 20 was a good compromise between calculation time and accuracy; this was used for subsequent calculations using the quasi-harmonic approximation. A comparison of free energies can be seen in Fig 5. Accuracies to within 1 meV, 0.5 meV and 0.1 meV require grid sizes of 2, 2 and 3 respectively.

Thermal Expansion of MgO Using the Quasi-Harmonic Approximation and Molecular Dynamics

**Fig 5.** Graph showing temperature vs free energy for the QMA (green) and MD (red)

The calculations were run at temperatures ranging from 0 K to 6000 K; past 2000 K when using the QMA, DLV gave error messages, possibly due to the dissociation of the lattice. The results are summarised in Fig 5. The two models give similar results at low temperatures, but diverge at T > 450 K. This is because anharmonicity increaces with increacing temperature, so the QM model breaks down (it predicts that free energy will fall with tempreature, whereas the MD model uses the formula $E_{k} = \frac{k_{B} T}{2}$ . This means that kinetic energy of the individual ions and thus the free energy increases with temperature, as shown in the graph and proven by experiment.^[1]

The Thermal Expansion Coefficient

The thermal expansion coefficient, $α_{V}$ , was calculated using by comparison of the cell volumes over temperatures ranging from 100 K to 1000 K; this can be seen in Fig 6. From this, it can be seen that while MD shows linearity all the way down to 100 K, the QMA deviates at temperatures below 200 K. Because $α_{V} = \frac{1}{V} {(\frac{δ V}{δ T})}_{p}$ , the units of volume cancel out and the coefficient is the gradient of the graph - for the QMA this is found to be $4.47 \times 1 0^{- 4} K^{- 1}$ , and for MD it is $5.81 \times 1 0^{- 4} K^{- 1}$ . This is in good agreement with literature.^[2]

The MD model is, in this case, more accurate. The QMA model uses only a single unit cell, and assumes periodicity of vibrations which means that the amplitude of the phonons must not cause interraction with adjacent cells. The MD model makes no such assumptions, and does not in fact consider the unit cell at all; rather, it computes the motion and position of every ion individually with respect to their immediate surroundings, to give a more physically realistic picture.

Conclusions

Both of the models used were accurate for determining the thermal expansion coefficient, whereas the QMA deviated enormously from experiment (and theory) at temperatures above 500 K. This can be explained by the fact that at higher temperatures, anharmonicity increaces greatly which causes innacuracies in the QMA. In order to improve the models, the ions could take into account orbitals rather than being approximated to round charges.

References

↑ http://www.asminternational.org/documents/10192/3449965/ACFAAD6.pdf/2d574bfc-e104-48c5-8d8e-c33ffe91c3a0
↑ NIELSEN, T. H., and M. H. LEIPOLD. "Thermal Expansion In Air Of Ceramic Oxides To 2200Oc". Journal of the American Ceramic Society 46.8 (1963): 381-387.

[1] ttp://www.asminternational.org/documents/10192/3449965/ACFAAD6.pdf/2d574bfc-e104-48c5-8d8e-c33ffe91c3a0

[2] NIELSEN, T. H., and M. H. LEIPOLD. "Thermal Expansion In Air Of Ceramic Oxides To 2200Oc". Journal of the American Ceramic Society 46.8 (1963): 381-387.

[1]

[2]