# Correlations in optical phonon spectra of complex solids

\[

## Abstract

Spectral correlations in the optical phonon spectrum of a solid with a complex unit cell are analysed using the Wigner-Dyson statistical approach. Despite the fact that all force constants are real, we find that the statistics are predominantly of the GUE type depending on the location within the Brillouin zone of a crystal and the unit cell symmetry. Analytic and numerical results for the crossover from GOE to GUE statistics are presented.

\]

Wigner-Dyson statistical analysis  has become a widespread approach to characterise quantum spectra of complex dynamical systems, such as nuclei, atoms and molecules, disordered quantum electron systems and electromagnetic resonances in chaotic cavities . Both experimental and theoretical studies of such systems point to the fact that the quantum energy levels universally obey the same mutual correlations and spectral rigidity as the eigenvalues of random Gaussian-distributed matrices. In the present paper, we apply such a statistical approach to characterise vibrational spectra of multi-component solids with a complex unit cell structure, in order to assess the correlation properties of their optical phonon spectra taken at various positions in the Brillouin zone (BZ).

In random matrix theory (RMT), different ensembles of matrices, which are defined by their fundamental symmetry, result in distinct level correlation properties. The two most common in condensed matter physics are the Gaussian ensemble of real symmetric matrices (GOE) which adequately describes the spectral correlations in quantum chaotic electron systems with time reversal symmetry, and the Gaussian ensemble of Hermitian matrices (GUE) which describes chaotic electron billiards with time-reversal symmetry broken by a magnetic field. One may expect a set of $`3N`$ coupled classical oscillators to be an example of a system with time-reversal symmetry, thus obeying the GOE spectral statistics. This is indeed the case with the acoustic spectroscopy of irregularly shaped solid resonators , or with the spectrum of a regularly shaped system consisting of coupled oscillators whose masses are random . Below we demonstrate the counter-intuitive result that, on the whole, the spectral correlations of the optical phonon modes associated with the same (albeit arbitrary) point of the BZ of a complex solid obey the GUE statistics. This complies with an earlier observation  on the electronic structure of highly excited bands in solids. In particular, we report a study of the distribution function $`P(s)`$ of the nearest-level-spacing in vibrational spectra of a complex solid based on numerical simulations of crystalline structures, both with and without mirror reflection symmetry in the unit cell. We also analyse the dependence of statistical properties on the phonon wave number $`𝐐`$ within the BZ, and obtain a detailed description of the crossover between the limiting regimes of GOE-type correlations specific to the phonon frequencies exactly at the center of the BZ where $`𝐐=0`$, and of GUE-type for sufficiently large values of $`Q`$.

To simulate a complex solid, we have adopted the following model. The unit cell of a crystal was taken in the form of a parallelepiped consisting of $`N=8\times 10\times 12`$ atoms with equal pair interactions but randomly chosen masses arranged on an fcc lattice. The unit cell size $`L`$ determines the periodicity of the Bravais lattice of the entire solid, $`𝐋=n_1𝐋_1+n_2𝐋_2+n_3𝐋_3`$. The spectrum of optical phonons in it can be found by determining all $`3N`$ normal modes corresponding to the linearised set of equations for coupled harmonic oscillators, $`m_𝐣\ddot{u}_𝐣^\alpha =_{\mathrm{𝐢𝐋}}K_{\mathrm{𝐢𝐣}}^{\alpha \beta }(𝐋)u_{𝐢+𝐋}^\beta `$, where $`𝐢,𝐣`$ are atomic positions within the unit cell and $`\alpha ,\beta `$ denote cartesian components. Fourier transforming we obtain for each given point in the BZ (wave number $`𝐐`$),

$$K_{\mathrm{𝐢𝐣}}^{\alpha \beta }(𝐐)=\underset{𝐋}{}k(l_{𝐋+𝐢,𝐣}^\alpha l_{𝐋+𝐢,𝐣}^\beta e^{ı\mathrm{𝐐𝐋}}4\delta _{𝐋+𝐢,𝐣}\delta ^{\alpha \beta }),$$
(1)   
where $`k`$ is the interatomic force constant, $`l_{𝐋+𝐢,𝐣}^\alpha =(𝐣𝐢𝐋)^\alpha /\left|𝐣𝐢𝐋\right|`$, and $`(𝐣𝐢𝐋)`$ belongs to the first coordination sphere of an fcc lattice. Note that the latter implies non-zero matrix elements only for nearest neighbors and applies restrictions to the sum over $`𝐋`$ in Eq. (1). The sets of phonon frequencies $`\{\omega _k(𝐐)\}`$ (which result in a spaghetti of $`3N`$ dispersion curves, each corresponding to a particular phonon branch) can be obtained by solving numerically the eigenvalue equation

$$\mathrm{det}(D\omega ^2I)=0,D(𝐐)=M^{1/2}K(𝐐)M^{1/2}.$$
(2)   
This equation also defines the dynamical matrix, $`D(𝐐)`$, in which the complexity in the composition of a solid is introduced by the random diagonal matrix $`M_{\mathrm{𝐢𝐣}}^{\alpha \beta }=m_𝐢\delta ^{\alpha \beta }\delta _{\mathrm{𝐢𝐣}}`$. For the masses $`m_𝐣`$ we have used a box distribution over the interval of masses $`[m\delta m,m+\delta m]`$ with $`\delta m/m=0.3`$. This corresponds to an r.m.s. of mass disorder $`\delta m^2^{1/2}/m0.17`$. The mean value $`m`$ defines the cut-off frequency $`\omega _c=\sqrt{8k/m}`$ of the vibrational spectrum in the disorder-free limit and sets a characteristic scale to measure the eigenfrequencies, whereas lengths are measured in units of the lattice constant. As explained below, such a model also allows us to exploit an analogy between the numerical results obtained by us and the properties of spectra of systems exhibiting quantum diffusion. To avoid complications brought into the problem by the localisation effects related to the vibrations of light atoms in a heavy matrix, we restrict our statistical analysis to the frequency range $`\omega /\omega _c<0.95`$.

At the middle of the BZ ($`𝐐=0`$) and at other symmetry points, such as $`𝐐=(0,\pi /L_2,\pi /L_3)`$, the dynamical matrix $`D`$ is real and symmetric. For an arbitrary $`𝐐\mathrm{𝟎}`$, the matrix elements of $`D(𝐐)`$ related to the sites on the edges of the unit cell (we use the nearest-neighbor interaction) acquire complex phase factors which make the whole dynamical matrix complex and Hermitian, $`D(𝐐)=D_S(𝐐)+ıD_A(𝐐)`$, where $`D_S(𝐐)`$, $`D_A(𝐐)`$ are real symmetric and antisymmetric matrices, respectively. Without imposing any spatial symmetries onto the unit cell, these two scenarios span the two limiting cases for the D-matrix symmetry; the difference between these two regimes manifests itself in the form of the normalized distribution function $`P(s)`$ of the nearest-level-spacing, $`s=(\omega _{k+1}\omega _k)/\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ is the mean level spacing.

We have constructed $`P(s)`$ employing a two-step averaging procedure. The first step is to use an ensemble of 50 random realizations of the distribution of masses in the sample. A further averaging of $`P(s)`$ is applied over a broad frequency range, namely, $`0.45<\omega /\omega _c<0.95`$, for each of the calculated spectra, by using an observation that they become statistically indistinguishable after having been rescaled by $`\mathrm{\Delta }`$ . Note that we do not analyze the lowest part of the spectrum because of poor statistics. The nearest-level-spacing distribution function for the optical phonon spectrum at $`𝐐=0`$ is shown in Fig. 1(a). It coincides with the random matrix theory prediction for real symmetric matrices given by Wigner-Dyson distribution function for the GOE . In contrast, optical phonons with $`𝐐0`$ exhibit a nearest-level-spacing statistics which is best fitted by the Wigner-Dyson distribution function for the GUE. A typical $`P(s)`$-histogram is plotted in Fig. 1(b), and is compared with the GUE analytical result . A similar observation has been made earlier by Mucciolo et al  in relation to the electronic band structures in crystals. For completeness, the inset in Fig. 1(c) shows the $`P(s)`$-histogram for a corner of the BZ, which is of GOE-type. Fig. 1(c) also illustrates the form of $`P(s)`$ in the intermediate regime between two distinct statistical classes. It is compared with the result of a fit based on the RMT prediction  for an interpolating ensemble between the GOE and GUE which contains a single fitting parameter.

The inter-ensemble crossover takes place at relatively small values of $`𝐐`$ (the parameter responsible for the imaginary part of the dynamical matrix), as it happens with a similar crossover in the energy spectra of chaotic electronic systems in the presence of a weak magnetic field . Note also that for small $`Q`$’s, the function $`P(s)`$ has the following asymptotic behavior: it resembles the GUE distribution function at $`s0`$, whereas it follows the GOE analytical result for large s. This is because the splitting of levels with $`s<1`$ is more sensitive to a small antisymmetric addition to a symmetric dynamical matrix than the splitting of rare pairs with $`s1`$. The above fact also implies that the crossover $`P(s)`$ has a form which cannot be simply reduced to a trivial mixing between two typical GOE and GUE distributions . The crossover studies were based upon the analysis of spectra of 600 random realizations of atomic masses in the unit cell for various values of $`𝐐`$ along $`𝐋_1`$. The result of numerical analysis of the GOE-GUE crossover in different intervals of vibrational spectra is shown at the bottom of Fig. 2 in the form of a gray scale shading of the parametric plane of $`q=𝐐𝐋_1=QL_1`$ and frequency $`\omega `$, where a darker colour stresses higher similarity of $`P(s)`$ to that of the orthogonal symmetry class, and white indicates the dominance of the GUE spectral statistics.

Below we present semiclassical arguments which provide the crossover value of the rescaled $`Q`$, $`q=QL`$, as function of the frequency of the mode considered. Our approach consists of viewing the dynamics as that of a wave-packet of lattice vibrations (in our case optical modes), spreading over a unit cell treated here as a region of disordered medium. This treatment is analogous to the semiclassical treatment of parametric spectral correlations developed in the studies of quantum disordered electron systems subjected to a weak magnetic flux . In the present analysis, it is assumed that the participation of each atom in a given optical phonon mode can be described semiclassically for time intervals shorter than $`t_H1/\mathrm{\Delta }`$, after which the discreteness of the spectrum for each value of $`𝐐`$ starts to dominate. The spread of vibrations over the unit cell and the role of the individual atoms in this dynamics is considered as a diffusion of waves through an fcc lattice with mass disorder and is determined by the interference pattern of a variety of equally probable diffusion paths . These paths are independent of the exact value of $`𝐐`$, whereas the phases of diffusive waves involved in such an analysis are large and random, so that one obtains correlated spectra .

However, the cyclic (albeit non-periodic) boundary conditions impose a torus geometry in the unit cell. If $`W`$ is the total number of windings around the torus then the type of correlations depends on the phase factor $`e^{i\delta \phi }`$ determined by the phase difference between a path with a positive W and its time-reversed counterpart with -W. Therefore, $`\delta \phi `$ gives us a measure of time-reversal symmetry breaking and is controlled by the ’external’ parameter $`q1`$. A relevant crossover parameter can then be found by estimating the r.m.s. value of the symmetry breaking phase $`\delta \phi `$, acquired by waves whose propagation is followed along a diffusive path with the maximal length allowed by the limit set by the time, $`t_H1/\mathrm{\Delta }`$ , where $`\mathrm{\Delta }1/L^3\nu `$ and $`\nu (\omega )`$ is the acoustic phonon density of states.

The parameter $`q`$ affects the phase of the partial amplitudes only for paths that cross the unit cell, i.e., paths whose ’winding’ number is non-zero. Hence, the total symmetry breaking phase is $`\delta \phi qW`$, where $`W=w_i`$ is a winding number made of approximately $`t_H/t_D`$ random contributions $`w_i=\pm 1`$. The time $`t_DL^2/D`$ with $`D`$ the diffusion coefficient, is typical of a diffusive spread over the unit cell of a vibration with frequency $`\omega `$. For our model, $`D=\frac{1}{3}s^2\tau `$ where $`s^2(\omega )`$ is the angular and polarization average of the squared phonon group velocity in the clean limit and $`\tau ^1`$ is the scattering rate, caused by the variation of atomic masses. The latter can be estimated  as $`\tau ^1(\omega )\omega ^2\nu \left(\delta m^2/m^2\right)`$. It follows from our considerations that an estimation of the r.m.s. value of $`\delta \phi `$ yields $`\delta \phi ^2^{1/2}q\sqrt{t_H/t_D}`$. The crossover can be assigned to $`\delta \phi ^2^{1/2}1`$. Therefore, the form of the crossover line is determined by

$$q_c\sqrt{t_D/t_H}\sqrt{\mathrm{\Delta }L^2/s^2\tau }\omega /\sqrt{s^2(\omega )},$$
(3)   
which is shown in the top of Fig. 2 and is in agreement with the numerically obtained gray scale plot.

The correlation properties of spectra of a solid are sensitive to the geometrical symmetries  of its unit cell structure, such as the existence of a mirror plane in it. Below, we extend the numerical analysis to a solid with $`n3`$ mirror symmetry planes in the unit cell, each characterised by unit vector $`\widehat{\eta }_i`$. Numerical simulations similar to the ones described above (with statistics collected from 100 random realizations of masses and averaging extended over the spectral window $`[0.45,0.95]`$) show that the effect of geometrical symmetries on the spectral statistics depends on the point in the BZ. This is because the phonon momentum allows for breaking both the orthogonal and point-group symmetries. Typical forms of the function $`P(s)`$ for each symmetry case are shown in Fig. 3, for various wave vectors $`𝐐`$. The solid lines in these figures illustrate the RMT result for overlapping sequences of GOE or GUE spectra with equal fractional densities . When $`𝐐=0`$, the spectra split into sequences of levels corresponding to different parities, so that $`P(s)`$ coincides with what one would expect for $`2^n`$ overlapping GOE’s (e.g., $`P(s0)=12^n`$).

For finite $`Q`$, spectral statistics is determined by the orientation of $`𝐐`$ with respect to the mirror symmetry planes, falling into one of the GOE or GUE classes. A summary of all distinct statistical limits is given in Table 1. Note that for $`𝐐`$ with all non-zero components perpendicular to a symmetry plane, correlations are of GOE-type. This effect is due to invariance of $`D(𝐐)`$ under the combination of reflection and complex conjugation operations, which results in a real representation of the dynamical matrix . Thus, for a unit cell with 3 mirror symmetry planes only the orthogonal ensemble statistics is realised for an arbitrary $`𝐐`$. The crossover parameter $`q_c`$ which determines which class from Table 1 should be expected, can be estimated using Eq. (3).

The authors thank I. Lerner and J. Pendry for discussions. This work was supported by EPSRC, and by European Union RTN and TMR programmes. Y.G. acknowledges support by the Centre of Excellence of the Israeli Academy of Sciences and Humanities.