Infrared optical properties of α quartz by molecular dynamics simulations

A subject of great current interest is that of a microscopic description of the ferroelectric transition. It is known that at the transition a divergence of the dielectric constant $\varepsilon(\omega)$ at $\omega=0$ occurs, which in most cases is understood as corresponding to the fact that the frequency of an infrared peak goes to zero. On the other hand, the infrared frequencies are usually studied via a linear analysis of phonon dispersion relations, while the nonlinear contribution to the dynamics should be relevant at the transition to the ferroelectric phase, as it should be near any phase transition. So the problem arises of how the infrared peaks should be described in a fully nonlinear setting. A description can actually be given through the study of the time autocorrelation of the polarization due to the ions, which can be computed in the classical approximation via molecular dynamics simulations. Computations of this type were indeed performed successfully in the case of LiF [1], with results that agree with the experimental data in a surprisingly good way.

At the moment we are unable to study ferroelectrics through molecular dynamics. So in this paper we limit ourselves to quartz, which is not a ferroelectric material, but however presents a divergence in the dielectric constant at the temperature of the $\alpha\text{-}\beta$ transition. We compute here its refractive-index curves in the infrared region. Our main concerns are the dynamical properties of the system, particularly at the $\alpha\text{-}\beta$ transition, and the quantitative agreement between calculated and experimental spectra in the infrared for α quartz.

As is well known [2], linear analysis shows that quartz is doubly refractive and that the peaks in its refractive-index curves correspond to the frequencies of the active normal modes. Such qualitative results are confirmed by the present MD simulations for α quartz at high temperatures, even at the transition to the β phase, notwithstanding the high nonlinearity of the system. This quantitative agreement of the nonlinear results with those of the linear analysis, in particular for the values of the frequencies, is quite surprising, in view of the large nonlinear contributions. On the other hand, both the linear analysis and our nonlinear study fail in perfectly reproducing the experimental data, as some systematic deviations are observed. We tried several procedures for choosing the parameters, both of the linear model and of the nonlinear one, in order to find a better agreement with the experimental curves. But there was no way of reducing the relative error below a threshold of the order of 6%. This fact, too, requires an explanation. These are the two main results of the present work.

It is known that the primitive cell of quartz contains nine atoms (three silicon and six oxygen atoms). Moreover, it has the form of a right prism of rhomboidal basis, corresponding to three basis vectors a, b, c with a and b forming an angle of 120° and c orthogonal to a and b. In table 1 we report the lengths of the three basis vectors at normal conditions of temperature and pressure (300 K and 1 bar), as given by [3], which we use in our simulations. We also report the fractional coordinates x, y, z (along the three basis vectors) of a silicon atom and of an oxygen atom, out of which all other coordinates can be generated by symmetry transformations¹ . The configuration corresponding to α quartz is thought of as being the more stable equilibrium configuration of the system. In order to simulate the crystal, we choose a domain $D\subset\mathbb{R}^3$ (fundamental box) constituted by $4\times 4\times 4$ primitive cells, with a number $N=9\times 4^3=576$ of point particles inside it. Due to the partially ionic character of the quartz crystal, the point particles have to be endowed with suitable effective charges, e_S for the silicon ion and e_O for the oxygen ion, with the neutrality constraint

$$\begin{equation} 2e_O+e_S=0. \end{equation} \tag{ 1 }$$

Thus, Coulomb long-range forces come into play and, in order to take them into account, working however with a small number N of particles, periodic boundary conditions are imposed. In addition to the electric forces, short-range two-body spherically symmetric potentials are introduced, one for each of the pairs Si-Si, O-O, Si-O. These potentials are taken of a form which is extensively used for quartz, namely (see [4,5]),

$$\begin{equation} U(r)= Ae^{-Br}-\frac{C}{r^{6}}, \end{equation} \tag{ 2 }$$

(r being the interatomic distance), with a triple of parameters A, B, C a priori different for each pair. In the numerical calculations, for the short-range interactions a cutoff of 9 Å was imposed, while the Coulomb interactions were dealt with through standard Ewald summations. This procedure is known to be necessary if the model has to reproduce the LO-TO splitting, namely, the splitting between longitudinal and transverse optical modes. A review of the short-range potentials employed for quartz can be found in [6,7].

The masses are taken from the literature, so that to fix the model there remains a total of 10 free parameters: one effective charge (for example that of oxygen) and the three parameters A, B, C of the short-range potential for each of the three pairs Si-O, O-O, Si-Si. The values of the potential parameters adopted are reported in table 2, while we used the values $e_S= 2.04191$ and $e_O=-1.02096$ for the effective charge (in unit of electron charge) of Si and O, respectively. The parameters were determined by optimization procedures aimed at obtaining the best possible agreement between the computed refractive-index curves and the empirical ones at 300 K, requiring in addition that the α structure be stable at that temperature and at larger (but not too much) temperatures. With the values thus determined for the parameters of the model it turns out, as will be shown below, that the $\alpha\text{-}\beta$ transition occurs at a temperature of about $T= 700\ \text{K}$, which is a rather low value.

The equations of motion were numerically solved using a Verlet algorithm, with integration step of 2 fs, at several values of temperature. Having chosen to work in a purely Hamiltonian frame, temperature was determined through the choice of the initial data. In principle this should be obtained by extracting the data according to a Gibbs distribution. This being impracticable, we followed the alternative standard procedure. Namely, one puts the particles at the α equilibrium point, while their velocities are extracted from a Maxwell-Boltzmann distribution at a suitable temperature. Then one lets the system thermalize, which usually takes a time of the order of 2 ps (1000 integration steps), and temperature is eventually identified through the mean kinetic energy of the ions. Then one starts computing means and correlations of the relevant quantities.

The refractive index is obtained by computing the electric permittivity tensor $\varepsilon_{ij}(\omega)$ as a function of frequency, and diagonalizing it at each given frequency. As expected, two eigenvalues are found to coincide, and the square root of such a value is precisely the refractive index of the ordinary ray. The refractive index of the extraordinary ray is instead the square root of the remaining eigenvalue².

The connection with dynamics is obtained through the susceptibility tensor $\chi_{ij}(\omega)$ due to the ions, which is related to permittivity by

$$\begin{equation} \varepsilon_{ij}(\omega)=\delta_{ij}+4\pi(\chi_{ij}\left(\omega)+\chi^{el}_{ij}\right)\!. \end{equation} \tag{ 3 }$$

Here, $\chi^{el}$ is the contribution of the electrons, which turns out to be constant in the infrared domain (see [8]). Instead, the ions' contribution $\chi_{ij}(\omega)$ is obtained numerically according to Green-Kubo linear response theory (see for example [9]) as follows. One considers the polarization P, which is defined in microscopic terms as

$$\begin{equation} \mathbf{P}= \frac 1{V} \sum_{l} e_l \mathbf{x}_{l}, \end{equation} \tag{ 4 }$$

where V is the volume of the simulation domain (or fundamental box), while $\mathbf{x}_{l}$ is the position vector of the l-th ion, of charge e_l.

Then at temperature T one has

$$\begin{equation} \chi_{ij}(\omega) = \frac V{k_BT}\int_0^{+\infty} e^{-i\omega t} \langle P_i(t) \dot P_j(0)\rangle \text{d} t, \end{equation} \tag{ 5 }$$

k_B being the Boltzmann constant. Here $\langle \cdot \rangle$ should in principle be the canonical average. Actually the averages were estimated as the mean of the time averages calculated along a certain number (usually 40) of different MD trajectories, calculated for 200 ps.

A first set of results is illustrated in figs. 1 and 2. In the first figure we report, vs. frequency, the real part of the refractive index for the ordinary ray at $T=300\ \text{K}$ (dark line) and at $T=1\ \text{K}$ (red line). The analogous spectra for the extraordinary ray are reported in fig. 2. At a temperature as low as $T=1\ \text{K}$ the spectrum is determined essentially by the linear approximation, so that the peaks correspond to the frequencies of the normal modes (actually, those of the so-called type A₂ and E; see [2]). The results show that the spectrum at $T=300\ \text{K}$ does not differ essentially from that corresponding to 1 K, apart from a consistent broadening of the peaks and some small shifts in their frequencies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A second set of results concerns the behavior of the spectra at the structural $\alpha\text{-}\beta$ phase transition which, in the present microscopic model, with the choice made for the parameters, turns out to occur in a region of temperatures roughly around $T=690\ \text{K}$. This will be shown in a moment. So we computed the refractive-index curves at $T=670\ \text{K}$, an $T=700\ \text{K}$, which are reported in fig. 3 for the ordinary ray and in fig. 4 for the extraordinary ray.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The figures show that, at the transition, the optical spectra are still dominated by the linear behavior. Indeed, in both figures the two curves relative to the two temperatures essentially superpose one another, and one can notice the same peaks of the previous figs. 1 and 2, just a little more broadened and noisy.

Some differences however show up. The most important one is the appearance of one more peak (see the inset) in the β phase at $\omega \simeq 455\ \text{cm}^{-1}$, that should correspond to a normal mode which is only Raman active in the α phase. In addition, a small peak appears at approximately 200 cm⁻¹ in the extraordinary ray, which should presumably be due to the nonlinearity.

So, on the one hand, the harmonic approximation essentially still dominates up to the transition temperature, at least for what concerns the refractive index. On the other hand we found that the transition has a relevant effect (in the model considered in this work) on a normal mode of the system, that with frequency $\omega=232\ \text{cm}^{-1}$ (corresponding in our model to the experimental frequency 207 cm⁻¹), which was predicted by Saksena [10] to be the soft mode triggering the transition. Such mode is Raman active, but not active in the infrared, and so does not show up in the refractive-index curves. In order to exhibit the impact of the nonlinearity on the dymamics of that mode in our model, we computed its power spectrum (i.e., the Fourier transform of the time autocorrelation of its amplitude), for two temperatures, 300 and 690 K. The results are reported in fig. 5. One sees that at 300 K one has, as expected, a rather sharp peak at about its characteristic linear frequency, with some small broadening, presumably due to an occuring of some chaoticity in the motion. Instead, at 690 K an overall large chaoticity shows up, extending up to zero frequency. Moreover, the original peak broadens and flattens, while a new one shows up at about 75 cm⁻¹. The latter feature seems to remind the form of the experimental Raman spectrum (see the curve at 800 K in fig. 1 of [11]).

Zoom In Zoom Out Reset image size

The occurrence of the transition is exhibited in terms of an "order parameter", which discriminates between the configurations of the two phases. Following essentially [12] and [13], we define it as follows.

Consider the representative position vectors of silicon and of oxygen, defined through their mean fractional coordinates (i.e., as the averages, over the elementary cells, of the fractional coordinates of such atoms). Consider also their time averages (actually a mean of such averages, taken over several independent simulations), which we denote by $\mathbf{x}_S$, $\mathbf{x}_O$. Then consider the position vectors $\mathbf{x}^\alpha_S$ and $\mathbf{x}^\alpha_O$, defined by the experimental fractional coordinates of silicon and of oxygen for α quartz, taken from table 1. Analogously define $\mathbf{x}^\beta_S$ and $\mathbf{x}^\beta_O$.

Thus the distance $d_{\alpha}$ of the mean configuration of the system from the equilibrium α configuration is naturally estimated as

$$\begin{equation*} d_{\alpha}=\sqrt{\|\mathbf{x}_S - \mathbf{x}^\alpha_S\|^2 + \|\mathbf{x}_O - \mathbf{x}^\alpha_O\|^2}, \end{equation*}$$

and analogously for the distance $d_{\beta}$ from the equilibrium β configuration. So one can introduce the variable u defined by

$$\begin{equation} u = \frac{d_{\alpha}-d_{\beta}}{d_{\alpha, \beta}} \end{equation} \tag{ 6 }$$

where ${d_{\alpha, \beta}}$ is a normalizing factor, the distance between the two equilibria, defined in the natural way.

A negative value of u clearly indicates that the atoms are, in the mean, near to the α configuration, while a positive value indicates that they are in the mean near to the β configuration. The graph of u vs. temperature is reported in fig. 6. One sees that a rather abrupt passage from a negative to a positive value occurs in a small region of temperatures about $T=690\ \text{K}$, at which a value of u very near to zero is obtained. Instead a value of about −0.6 is obtained at $T=670\ \text{K}$ and a value of about +0.25 is obtained at $T=700\ \text{K}$. This shows that at these temperatures the nonlinear effects become so important as to trigger a phase transition.

Zoom In Zoom Out Reset image size

In figs. 7 and 8 we report both the calculated refractive-index curves and the experimental ones, taken from [8,14], for the ordinary and the extraordinary rays respectively, at 300 K. For both types of rays the experimental and the calculated curves have the same general aspect: namely, the number of peaks is the same, and both the intensities and the broadening are of the same order. Actually the lowest peak in the theoretical curve corresponds to the normal mode at 140 cm⁻¹, which has a vanishingly small intensity in the data; on the contrary the lowest frequency peak in the experimental curve is at 265 cm⁻¹ and the corresponding mode in the theoretical curve has a very low intensity (see the normal modes frequencies in table 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3:.  Normal modes of α quartz grouped by symmetry. Experimental frequencies (from [15]) and calculated frequencies in cm⁻¹ are reported, together with the calculated components of the dipole moment μ for each mode.

Sym.	Exp. freq.	Calc. freq.	$\mu_x$	$\mu_y$	$\mu_z$
	1162	1125	−0.1947	0.1903	0
		1125	−0.1903	−0.1947	0
	1072	1083	0.3479	−0.4105	0
		1083	−0.4105	−0.3479	0
	795	727	0.3032	−0.1175	0
		727	−0.1175	−0.3032	0
	697	675	0.1058	0.1061	0
		675	−0.1061	0.1058	0
E	450	506	0.4470	0.2172	0
		506	0.2172	−0.4470	0
	394	359	−0.3075	0.0309	0
		359	−0.0309	−0.3075	0
	265	258	−0.0010	0.0010	0
		258	0.0010	0.0010	0
	128	141	0.0435	0.0072	0
		141	0.0072	−0.0435	0
	1080	1101	0	0	−0.5871
A2	778	732	0	0	0.3754
	495	544	0	0	0.4078
	364	358	0	0	−0.4294
	1085	1074	0	0	0
A1	464	474	0	0	0
	356	355	0	0	0
	207	232	0	0	0

However there is a clear quantitative disagreement, determined essentially by the positions of the peaks. A similar quantitative disagreement was met in the values of the elastic moduli (which we computed as a back check from the values of the parameters adapted to fit the frequencies) following Born and Huang [17]. The values are reported in table 4, together with the experimentally measured values taken from [16]. By inspection one can see that the computed values differ by a factor of order 2 from the experimental ones.

As in our model we have ten free parameters, one can investigate whether a better quantitative agreement can be obtained by optimizing them, or even by considering other types of models. We concentrated our attention on the frequencies of the spectral lines. However, notwithstanding a lot of effort, we were unable to significantly improve the agreement.

We now describe the strategy we followed for optimizing the parameters. The procedure is quite involved, because we have two objectives. On the one hand the system has to admit a global equilibrium configuration, periodic with respect to the primitive cell and furthermore reproducing the α quartz symmetries. On the other hand we require that the normal mode frequencies, calculated at the α equilibrium, reproduce the frequencies observed, both in the infrared spectra and in the Raman ones.

Our procedure was the following one. To start up, we consider the experimental crystallographic configuration of table 1, given by X-ray diffraction, and linearize the equations of motions at that point. So we can determine the normal modes frequencies at the corresponding equilibrium point, as functions of the parameters entering the potential. In this calculation, the symmetries of the crystal are automatically taken into account in the construction of the dynamical matrix, because only some components are directly calculated, all the others being derived by symmetry transformations. As a consequence, the modes are correctly grouped into the three irreducible representations of the symmetry group, namely 4 A₁ modes, 4 A₂ modes and 8 degenerate E modes, so that a total of 16 distinct frequencies are obtained. Such properties are reflected in the components of an electric dipole moment vector μ that we associate to each mode by multiplying the Cartesian displacement of each atom by its effective charge and summing all the vectors thus obtained (see table 3). This gives an easy criterion for the correct identification of frequencies in the minimization procedure.

Then, for any single set of parameters obtained we determine the corresponding equilibrium position and the corresponding set of normal mode frequencies. The set of parameters is accepted if the calculated equilibrium position is sufficiently near to the experimental one of α quartz, the frequencies are sufficiently near to the experimental ones, and furthermore the α structure is stable up to sufficiently high temperatures. No set of parameters found gave an agreement for the frequencies better than about 6–7%.

Other attempts were as follows. We started from the power n = 6 entering (2), letting n be a free parameter, different for each pair, adapting the cutoff parameter to each choice. No substantial improvement was obtained. Then we changed completely the form of the potentials, using Lennard–Jones ones. But this gave a drastic worsening of the results.

These facts show that the results depend in a very sensitive way on the form of the potentials. In order to bypass this problem we decided to restrict our studies to the linear model, assigning as parameters directly the elastic constants. In such a way one even eliminates the constraint that the elastic constants should be defined in terms of first and second derivatives of a given potential. However, no progress was obtained.

As a last resort, we eliminated the neutrality constraint (1) on the effective charge, by assuming both charges to be free parameters, but again without substantial improvement.

Actually, in the whole literature we were unable to find a paper in which the calculated frequencies agree with the experimental ones, in the mean, better than 3%, which is the result obtained in the old paper [18]. In such a paper a linear model is considered, which takes into account also the polarization of oxygen ions as a free parameter. However, the maximum error was larger than 7%.

For what concerns nonlinear models investigated by MD simulations (see the review [6]), the situation is worse. The error is of the order of 7%–8% for pair potentials (see, for instance, [19–21]), and does not decrease significantly when three-body potentials are considered (see [13,22]).

We interpret these facts as indicating that some structural deficiency is present in all models (including ours) that have been considered. Such a deficiency was sometimes acknowledged (see, for example, [23]) and is usually ascribed to some deficiency in the short-range potentials. Our opinion is instead that a relevant role is played by the long-range forces, since the main discrepancy concerns essentially the LO-TO splitting, which is due to the long-range forces. Our conjecture is that a significant improvement could be obtained if one takes into account retardation, which already proved to be an essential feature for a microscopic description of polaritons (see [24]).

We thank G. Grosso, G. Pastori Parravicini, N. Manini and G. Onida for useful discussions. The use of computing resources provided by CINECA is also gratefully acknowledged.