Rep:MgO:syl815

Introduction

Thermal properties of materials are described by their vibrational free energies, which can be described in terms of the relative motion of atoms or the motion of their centre-of-mass. These concepts give rise to different approaches in calculating vibrational free energy, and both methodologies will be explored in greater detail.

Quasi-harmonic Approximation (QHA)

Fundamentally, QHA invokes the description of a crystalline solid as a primitive unit cell. This is essential due to the impracticality of calculating all the vibrational degrees of freedom in a crystal—for a crystal of size $N$, there are $3N$ degrees of vibrational freedom, and in an infinitely large crystal lattice, 3N --> infinity. Nonetheless, the translational periodicity of the crystal lattice, where $f(r+T)=f(r)$, simplifies the dynamics of all atoms in the lattice into that of a unit cell. For such a simplification to be appropriate, the following assumptions are made.

The Adiabatic Approximation

The adiabatic approximation separates the motion of the ion cores from that of the electrons since former are much more massive than the latter. Hence, the ion cores can be assumed to be in their equilibrium positions and that their motion is dependent on the potential field generated from the average motion of electrons.

The Harmonic Approximation

The total potential energy of a crystal can be expressed as the sum of all interatomic potentials. A two-body system typically has an anharmonic potential energy surface $U\left(r\right)$, where $r$ is the interatomic separation. By considering a small displacement $x$, $$x=r-r_0\left(1\right)$$

Where $r_0$ is the equilibrium distance between the first and second atoms and is a minimum on the potential energy surface.

$U\left(r\right)$ can be expanded in a Taylor series about $r_0$. $$U\left(r\right)=U(r_0)+\frac{\partial U}{\partial x}x+\frac{{\partial }^{2}U}{\partial x^2}{x}^2+...$$ Since $U(r_0)$ is unimportant in dynamics, $\frac{\partial U}{\partial x}$ is a force term and must be 0 for an equilibrium configuration, and all higher order terms $x_n$ , where $n\ge 3$ are assumed to be close to 0. As such, only the quadratic term is considered in the harmonic approximation. The solutions are the normal modes of vibrations for a system of independent quantum oscillators.

A phonon is a quantum of vibrational energy, hw, associated with a wave vector k.

Hence, for a crystal, its potential energy $\mathrm{\Phi }$ is given in the following equation under the harmonic approximation $$\mathrm{\Phi }={\mathrm{\Phi }}_0+\sum _{lk}\sum _{\alpha }{\mathrm{\Phi }}_{\alpha }\left(lk\right)\cdot u+\alpha \left(lk\right)+\frac{1}{2}\sum _{l{l}^{\prime }k{k}^{\prime }}\sum _{\alpha \beta }{\mathrm{\Phi }}_{\alpha \beta }(lk,l^{\prime }{k}^{\prime })\cdot u_{\alpha }\left(lk\right)u_{\beta }\left(l^{\prime }{k}^{\prime }\right)+...$$ Where $l$ and $k$ are the labels of the unit cells and atoms in each unit cell respectively, $\alpha$ and $\beta$ the Cartesian coordinates. ${\mathrm{\Phi }}_0$, ${\mathrm{\Phi }}_{\alpha }$ and ${\mathrm{\Phi }}_{\alpha \beta }$ represent the zeroth, first and second order force constants respectively.

Limitations of Harmonic Approximation

The harmonic approximation predicts symmetric atomic vibrations about r0 at all temperatures, and is therefore incongruent with observed phenomena such as thermal expansion and heat conductivity. The QHA causes renormalisation of the phonon frequencies and atomic force constants as is appropriate for the thermal equation of state.

MD Simulation

MD considers the forces exerted on each atom and provides a classical description of an $N$-atom system. This is given by $$M\left(\frac{{\partial }^2{\mathbf{r}}_{\mathbf{i}}}{\partial t^2}\right)=\sum _{j=1,j\ne i}^N{\mathbf{F}}_{\mathbf{i}\mathbf{j}}$$

where ${\mathbf{r}}_{\mathbf{i}}$ is the distance between atoms where $i$ and $j$, ${\mathbf{F}}_{\mathbf{i}\mathbf{j}}$ is the force exerted on $i$ by $j$ and $M$ is the total mass of the system.

Methodology

Unless otherwise stated, all calculations were performed on a primitive unit cell of MgO with lattice parameters $a=2.9783\AA$, $a=b=c$; $\alpha =60$°, and $\alpha =\beta =\gamma$ with GULP version 1.4.43 and crystals visualised with DLV interface.

A phonon dispersion curve was computed by sampling 100 points within the first Brillouin zone. The phonon density of states (DOS) was calculated with various shrinking factors, and the graphs subsequently plotted with matplotlib. The free energy of MgO was calculated with different shrinking factors at 300 K, and a suitable shrinking factor selected for the subsequent investigation of the thermal expansion of MgO. For every run, the Gibbs free energy was optimised, and calculations were performed from 0 to 2960 K in temperature steps of 20 K.

All MD simulations were performed on an isothermal-isobaric ensemble of MgO supercell of 32 formula units, with the following cell parameters: $$a=8.4239\AA$$ $$a=b=c$$ $$\alpha =90^o$$ $$\alpha =\beta =\gamma$$

MD was performed over a temperature range of 20 K to 4000 K, with temperature steps of 20 K. All calculations were performed with a time step of 1 fs. From 20 K to 1680 K, the system was allowed to first equilibrate for 1 ps; this was increased to 5 ps from 1700 K to 4000 K. Following which, MD production was allowed to run for 5 ps for all temperatures.

All data was analysed with Python on Jupyter notebook, and all graphs plotted with matplotlib.

Results and Discussion

The lattice energy of MgO was calculated to be -41.0753 eV per primitive unit cell.

Phonon Modes of MgO

Figure 1 illustrates the phonon dispersion curve computed at 100 points for the primitive unit cell.

Figure 1 Phonon Dispersion Curve of MgO

A salient feature is the presence of 6 branches in the dispersion diagram. Assuming that the Born-von Karman boundary condition is satisfied, the edge effects of cells on dynamics can be ignored and $u_{N+1}=u_1$, where $u$ is the displacement and $N$ is the number of unit cells. This also implies the translational symmetry in $k$-space, such that all information of phonon dispersion can be derived by sampling in the first Brillouin zone (FBZ).

By considering a linear diatomic chain satisfying the periodic boundary condition, the solutions to the vibrational frequency can be expressed in the form $${\omega }^2=\mathrm{\Lambda }(\frac{1}{m}+\frac{1}{M})\pm {[{(\frac{1}{m}+\frac{1}{M})}^2-\frac{4}{mM}{\sin }^{2}ka]}^{\frac{1}{2}}$$

where $\omega$ is the angular frequency, $\mathrm{\Lambda }$ is the force constant of the bond, $(\frac{1}{m}+\frac{1}{M})$ is the reduced mass of the system, where $m$ is the mass of the lighter atom (O) and $M$ is the mass of the more massive atom (Mg), and $a$ is the length of the unit cell.

The equation highlights two possible solutions for each $k$-value in a linear chain. Moreover, when $m\ne M$, a gap is observed at$k=\frac{\pi }{2a}$, which is observed in Figure 1.

By extending the logic to a 3D crystal lattice, the number of branches observed is given by $3x$, where $x$ is the number of atoms per unit cell. This is in agreement with the observation in Figure 1.

By appraising the solutions for $k=0$ (long wavelength limit),

$${\omega }_1=2\mathrm{\Lambda }(\frac{1}{m}+\frac{1}{M})$$

$${\omega }_2=0$$

${\omega }_1$ corresponds to a high energy mode in which the atoms in the unit cell are moving out-of-phase, where frequency values are within the visible electromagnetic spectrum. The atoms are able to interact with an electric field of appropriate frequency due to the presence of both a positive and negative charge within the unit cell. It is hence naturally termed the optical mode.

On the other hand, ${\omega }_2$ corresponds to a low energy mode with the atoms moving in phase and the wave pattern is similar to sound waves—hence the term acoustic mode. For any crystal with $N$ atoms in the unit cell, there are only 3 acoustic—2 transverse and 1 longitudinal—and $3N-3$ optical branches. The transverse modes are perpendicular to $k$, while the longitudinal mode is parallel.

Computing Density of States (DOS)

The impracticality of sampling all $k$-points within the FBZ can be circumvented by the use of a commensurate grid of $k$-points. To determine this set of $k$-points, the Pack-Monkhorst (PM) shrinking factor was used to specify the number of equidistant $k$-points taken along each direction of $b_1$, $b_2$ and $b_3$ in one reciprocal lattice primitive unit cell.

A major advantage is its computational efficiency by restricting the number of $k$-points calculated to a finite value. Moreover, the accuracy obtained from calculations with a PUC can be comparable to that of a supercell as long as the shrinking factor is appropriate.

Table 1 illustrates the effect of modifying the PM shrinking factor on the number of k-points calculated.

Table 1. Grid size against number of $k$-points

Grid Size (n x n x n)	Number of k-points
1	1
2	4
3	18
4	32
5	75
6	108
8	256
10	500
16	2048
20	4000
64	>99 999

As the mesh of $k$-points increases, the number of $k$-points calculated increases as well. This is contrary to the prediction from the above equation, where we would expect $k_x\times k_y\times k_z$ number of points. This can be attributed to the mapping of equivalent $k$-points onto each other and thus the number of $k$-points calculated is reduced.

Determining Optimal Grid Size for MgO

The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of $k$-points is increased.

Figure 2 DOS plots against wavenumber when PM shrinking factor was varied.

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication which will be further discussed in the later section.

An initial plot of the density of states was obtained from a $1\times 1\times 1$ grid yielding six resultant modes. Sharp and distinct peaks are observed in the plot, since only one $k$-point was sampled.

Notably, only four unique peaks are observed even though we should observe 6 modes of vibrations. The final two modes are degenerate at $286c{m}^{-1}$ and $351c{m}^{-1}$. Compared to the non-degenerate acoustic and optical peaks ($676c{m}^{-1}$ and $806c{m}^{-1}$ respectively), the degenerate acoustic modes are higher in energy whereas the degenerate optical modes are lower in energy correspondingly. It can therefore be deduced that the degenerate acoustic and optical modes are transverse in nature.

The $k$-point used in the DOS calculation could be identified by comparing with the dispersion curve. Since point M contains all of the frequency values in Figure 1, it can be determined that the point represented in the DOS curve is M, where $k_x$ = 0.5, $k_y$ = 0.5 and $k_z$ = 0.5.

Relationship between the Dispersion Curve and DOS

The DOS curve illustrates the number of energy states per unit energy, demonstrating a mode at $414c{m}^{-1}$. This correlates well with Figure 1. By constructing a horizontal line at frequency = 414 cm-1, it can be observed that the branches intersect this line frequently. This implies that a significant proportion of k-points have vibrational modes of frequency $414c{m}^{-1}$. The DOS curve can thus be interpreted as the orthogonal of the dispersion curve.

The dispersion diagram is useful in locating the band gaps of the acoustic and optical modes - for electronic dispersion diagram, this is useful in identifying whether a material has a direct bandgap or an indirect one, which affects the properties of the material and its use.

However, the dispersion diagram only illustrates the energy values calculated at the special points chosen, interpolating the energies of the vibrational modes for the k-points which are not calculated. The DOS plot is in this respect more meaningful, the energy states for all of $k$ values are accounted in this representation.

Computing the Free Energy Using the Harmonic Approximation

Figure 3 demonstrates the relationship between the PM shrinking factor used and the computed Helmholtz free energy of the system.

Figure 3 Free energy vs. PM shrinking factor

The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased.

From the above figure, the free energy of MgO is observed to increase and converge to a value of -40.926 483 eV, and it is observed that this occurs for a grid size of 8x8x8.

A 2x2x2 grid is sufficient for calculating the free energy of MgO to 1 meV. A 4x4x4 grid is necessary for a precision to 0.5 meV and 0.1 meV.

Thermal Expansion

The Helmholtz free energy of a crystal is given by the sum of the energies of independent vibrational waves. The energy level $n$ of a quantum harmonic oscillator are given by:

$$E_n=(n+\frac{1}{2})h\nu$$

where $h$ is Planck's constant and $\nu$ is the frequency of energy level $n$. For $N$ atoms and $3N$ independent harmonic oscillators, the vibrational energy is given by:

$$E_{vib}=\sum _{n=1}^{3N}(n+\frac{1}{2})h\nu$$

For a canonical $NVT$ ensemble, the partition function is given by

$$Z=\sum _n{e}^{(-\beta {\epsilon }_n)}=\frac{e^{-\frac{\beta {\epsilon }_n}{2}}}{1-e^{-\beta {\epsilon }_n}}$$

where $\beta =\frac{1}{kT}$ and $E_n$ enumerates all vibrational energy states.

For $N$ atoms and $3N$ independent harmonic oscillators,

$$Z_N=\prod _n^{3N}{e}^{(-\beta {\epsilon }_n)}=\frac{e^{-\frac{\beta {\epsilon }_n}{2}}}{1-e^{-\beta {\epsilon }_n}}$$

The phonon entropy can then be expressed in terms of the partition function: $$S=-k_B\ln Z_N$$ where $k_B$ is the Boltzmann constant.

Given the relation $$F=U+TS$$ where $U$ is the internal energy of the system— for a crystal this is its electric potential energy $$U_E=\frac{1}{2}\sum _{i=1}^{3N}{q}_i\sum _{j=1,j\ne i}^{3N}\frac{q_j}{4\pi {\epsilon }_0{r}_{ij}}$$

where $i$ and $j$ are the indices of the ions, $r_{ij}$ is the distance between $i$ and $j$ and ${\epsilon }_0=8.8542\times 10^{-12}F\cdot m^{-1}$

The Helmholtz free energy of a crystal is thus given by

$$F=\frac{1}{2}\sum _{i=1}^{3N}{q}_i\sum _{j=1,j\ne i}^{3N}\frac{q_j}{4\pi {\epsilon }_0{r}_{ij}}+\sum _{n=1}^{3N}(n+\frac{1}{2})h\nu +k_{B}T\sum _n^{3N}\ln (1-e^{-\beta {\epsilon }_n})$$

This equation could be used to qualitatively rationalise the free energy dependence on temperature. The data obtained is plotted in Figure 4.

Figure 4 Free Energy of System against Temperature

Particularly, there are two salient regimes of interest. At low temperatures, T < 100 K, the graph is flat. However, at high temperatures, the behaviour is approximately linear. These observations are in agreement with the above equation, which highlights the temperature dependence of entropy $S$. At low temperatures, the term $k_{B}T\sum _n^{3N}\ln (1-e^{-\beta {\epsilon }_n})$ is extremely small, and hence the free energy term is dominated by the internal energy of the crystal. At high temperatures, the term $-k_{B}T\sum _n^{3N}\ln (1-e^{-\beta {\epsilon }_n})$ dominates and therefore the free energy of the system appears to have a dependence in temperature.

Variation of Lattice Parameter with Temperature

Figure 5 Cell parameter of MgO against Temperature

As the temperature increases, the lattice parameter increases. It can thus be observed that the cell volume $V$ has a dependence on temperature $T$, and the thermal expansion coefficient $\alpha$ is given by

$$\alpha =\frac{1}{3V}{\left(\frac{\partial V}{\partial T}\right)}_P=\frac{1}{3B}{\left(\frac{\partial P}{\partial T}\right)}_V$$ Where $B$ is the bulk modulus and $P$ is the pressure. At 300 K, $\alpha =2.06\times 10^{-5}{K}^{-1}$, compared to a literature value of $\alpha =3.06\times 10^{-5}{K}^{-1}$.

MD Simulation

The cell volume per formula unit of MgO was plotted against temperatures between 20 K to 4000 K.

Figure 5 Experimental Data MD

Under MD, the cell volume generally increases linearly with temperature throughout. By considering the mean kinetic energy of the crystal $$⟨E_k⟩=\frac{1}{2}M\sum _{i=1}^N{v}_i^2=\frac{3}{2}N{k}_B{T}_{MD}$$ Where $⟨E_k⟩$ is the average kinetic energy of the atoms, $M$ is the mass of the crystal lattice, and $v_i$ represents the velocity of the atom $i$. It can be observed that the cell energy is linearly dependent on temperature. In a constant pressure system, this would predict volume expansion as temperature increases under classical MD.

It can be observed that at high temperatures when $T\ge 2000K$, more noise is present in the data due to the large cell volume and the large kinetic energy of the atoms. A longer equilibration time might be necessary to minimise this effect.

Figure 6 Comparison of cell volume per formula unit under QHA and MD calculations.

At extremely low temperatures of $T\le \approx 200K$, QHA predicts a larger cell volume than MD. This can be attributed to the significant quantum effects when the temperature of the system is below Debye temperature. However, when performing MD, the average temperature of the system is only dependent on the mean kinetic energy of the atoms and neglects the zero point vibrations at that temperature. Consequently, classical MD predicts a smaller cell volume with the atoms closer together.

Nonetheless, this can be circumvented through quantum corrections. Wang et al. described one such approach involving a scaling correction of the system temperature.

The data obtained for MD and QHA demonstrate strong agreement for temperatures between 200 to 1000 K. At these temperatures, the thermal energy of the system is sufficiently large such that the motion of the particles can be described classically.