Application of the random matrix theory to the boson peak in glasses

TL;DR

This paper applies random matrix theory to analyze vibrational states in glasses, deriving equations that explain the boson peak and match numerical results, providing insights into the vibrational properties of amorphous materials.

Contribution

The study introduces a new approach using correlated random matrix theory to derive the vibrational density of states in glasses, accurately modeling the boson peak.

Findings

Vibrational density follows Debye law at low frequencies.

The boson peak appears as an additional feature at higher frequencies.

Derived equations match numerical simulations well.

Abstract

The density of vibrational states $g(\omega)$ of an amorphous system is studied by using the random-matrix theory. Taking into account the most important correlations between elements of the random matrix of the system, equations for the density of vibrational states $g(\omega)$ are obtained. The analysis of these equations shows that in the low-frequency region the vibrational density of states has the Debye behavior $g(\omega) \sim \omega^2$. In the higher frequency region, there is the boson peak as an additional contribution to the density of states. The obtained equations are in a good agreement with the numerical results and allow us to find an exact shape of the boson peak.

Figures

Equations

$$⟨B_{ij}{B}_{kl}⟩=\frac{1}{N}{C}_{ik}{\delta }_{jl},$$

$$\begin{gathered} \rho (\epsilon )=\frac{1}{8{\pi }^3}{\int }_{-\pi }^{\pi }{\int }_{-\pi }^{\pi }{\int }_{-\pi }^{\pi }\delta \left(\epsilon -8{\sin }^2\frac{k_x}{2}-8{\sin }^2\frac{k_y}{2}-8{\sin }^2\frac{k_z}{2}\right)d{k}_{x}d{k}_{y}d{k}_z \\ ={\int }_0^{\infty }{J}_0^3(\xi )\cos (3\xi \epsilon )d\xi , \end{gathered}$$

$$\int \frac{\rho (\epsilon ){\epsilon }^2}{(X-\epsilon {)}^2+Y^2}d\epsilon =1+\kappa ,$$

$${\omega }^2=(X^2+Y^2)\int \frac{\rho (\epsilon )\epsilon }{(X-\epsilon {)}^2+Y^2}d\epsilon ,$$

$$g(\omega )=\frac{1}{\pi }\frac{2\omega Y}{X^2+Y^2}.$$

$$2X{I}_1+\frac{(X^2-Y^2)}{16\pi Y}\sqrt{X+\sqrt{X^2+Y^2}}=\kappa ,$$

$${\omega }^2=(X^2+Y^2)\left(I_1+\frac{X}{16\pi Y}\sqrt{X+\sqrt{X^2+Y^2}}\right),$$

$$I_1=\int d\epsilon \frac{\rho (\epsilon )}{\epsilon }=\frac{1}{4}{\int }_0^{\infty }{J}_0(\zeta )\sin (3\zeta )d\zeta \approx \mathrm{0.126365.}$$

$$2X{I}_1+\frac{X^2}{16\pi Y}\sqrt{2X}=\kappa ,$$

$${\omega }^2=X^2\left(I_1+\frac{X}{16\pi Y}\sqrt{2X}\right),$$

$$g_l(\omega )=\frac{1}{\pi }\frac{2\omega Y}{X^2}.$$

$$g_l(\omega )=\frac{\omega }{8{\pi }^2\sqrt{I_1\kappa }}\sqrt{\frac{1-\sqrt{1-{\omega }^2/{\omega }_c^2}}{1-{\omega }^2/{\omega }_c^2}},$$

$$g_D(\omega )=\frac{{\omega }^2}{2{\pi }^2(2\kappa {)}^{3/2}}+\mathcal{O}({\omega }^4).$$

$${\omega }^2=(X_c^2+Y^2)I_1,$$

$$g_r(\omega )=\frac{1}{\pi }\frac{2\omega Y}{X_c^2+Y^2}.$$

$$g_r(\omega )=\frac{2\sqrt{I_1}}{\pi }\sqrt{1-\frac{{\omega }_c^2}{{\omega }^2}}.$$

$$g_D=\frac{1}{2{\pi }^2}{E}^{-3/2}{\omega }^2.$$

$$E=2\kappa =2\left(\frac{K}{N}-1\right).$$

$${\omega }_b=\sqrt{\frac{3}{2}}{\omega }_c=\sqrt{\frac{3}{8I_1}}\kappa .$$

$${\omega }_b=\sqrt{\frac{3}{32I_1}}E.$$

Full text

1 Introduction

The finding of common vibrational properties in amorphous dielectrics (glasses) is one of the key problems in the physics of disordered systems. One of these properties is a boson peak in the reduced density of vibrational states [1]. It characterizes a presence of the excess vibrational density of state compared to the Debye law.

The boson peak is observed in many experiments, such as the Raman scattering, X-ray scattering, an inelastic neutron scattering, and the measurements of the heat capacity $C(T)$ as a maximum in $C(T)/T^{3}$. The boson peak frequency $\omega_{b}$ is correlated with Ioffe-Regel crossover frequency $\omega_{\textsc{ir}}$ when the phonon mean free path becomes of the order of the phonon wavelength [2, 3, 4, 5, 6]. Despite of a number of articles about the boson peak, still there is no common physical interpretation. In the present paper, to study the nature of the boson peak we apply an approach based on the use of random matrices.

2 Random matrix model

In amorphous solids, due to the local disorder, the symmetric force constant matrix $\hat{\Phi}$ has a random nature. Therefore, the dynamical matrix $M_{ij}=(m_{i}m_{j})^{-1/2}\Phi_{ij}$ [7] has some disorder as well. However, the dynamical matrix $\hat{M}$ has several physical constraints. First of all, there is a translational invariance (or a sum rule), which means that the sum of all elements in each row or column of the matrix $\hat{M}$ must be equal to zero (we consider unit masses $m_{i}=1$). Secondly, eigenvalues of the dynamical matrix $\hat{M}$ determine the eigenfrequencies squared and must be nonnegative [8] due to the requirement of the mechanical stability.

Taking into account the above properties, the dynamical matrix $\hat{M}$ can be written as $\hat{M}=\hat{A}\hat{A}^{T}$ [9]. The matrix $\hat{A}$ is some non-symmetric random matrix with a constraint $\sum_{i}A_{ij}=0$. In general, the matrix $\hat{A}$ can be a rectangular $N\times K$ random matrix.

Let us consider a disordered system on a simple cubic lattice with the lattice constant $a_{0}=1$. In this model each atom is placed in a site of the lattice, but the interaction between atoms is random. Such a system can be built by the following structure of the matrix $\hat{A}$. The non-diagonal elements $A_{ij}$ are Gaussian random numbers if atoms with indices $i$ and $j$ are nearest neighbors in the lattice. Other non-diagonal elements are zeros. The diagonal elements are defined as $A_{ii}=-\sum_{j\neq i}A_{ji}$ to satisfy the condition of the translational invariance.

As it was shown in [9, 14], the system, described by the square matrix $\hat{A}$, has a zero Young modulus. The point is that each column of $\hat{A}$ corresponds to one bond of the system, each row of $\hat{A}$ corresponds to one degree of freedom of the system. Therefore, for the square matrix $\hat{A}$, the number of bonds and the number of degrees of freedom are equal. This is the so-called Maxwell rule. In this case, the system is very soft and do not have a macroscopic rigidity. In such a system, phonons do not propagate, and the boson peak is absent too [9, 14].

In order to introduce the macroscopic rigidity and the boson peak into the matrix model $\hat{M}=\hat{A}\hat{A}^{T}$, we build a new random rectangular matrix $\hat{B}$ with size $N\times K$ with $K>N$. We take two independent realizations of the square random matrix $\hat{A}$ with size $N\times N$: $\hat{A}$ and $\hat{A^{\prime}}$. In order to generate the matrix $\hat{B}$, we randomly insert some number of columns of the matrix $\hat{A^{\prime}}$ into the matrix $\hat{A}$. This random insertion of the new columns corresponds to a random addition of new bonds to the vibrational system. The relative number of new bonds is $\kappa=(K-N)/N$, where $K$ is the total number of columns in the resulting matrix $\hat{B}$. Therefore, the new rectangular matrix $\hat{B}$ satisfies conditions of translational invariance and mechanical stability and describes the vibrational system as $\hat{M}=\hat{B}\hat{B}^{T}$.

The analysis of the boson peak with a help of the random matrix theory can shed light on its universal nature. However, the nontrivial structure of the matrix $\hat{B}$ significantly complicates the problem. Nevertheless, the model under consideration has simpler correlations between matrix elements than three similar random matrix models, which were considered in [9, 14]. In the next section, we present equations for finding the vibrational density of states (VDOS) $g(\omega)$, obtained with an account of the most important correlations between elements of the matrix $\hat{B}$.

3 Density of vibrational states

The pairwise correlations between matrix elements in the matrix $\hat{B}$ has a form

[TABLE]

where the averaging is performed by different realizations of the matrix $\hat{B}$ with the fixed parameter $\kappa$. The matrix $\hat{C}$ has the following structure. The non-diagonal elements $C_{ij}=-2$ if atoms with indices $i$ and $j$ are nearest neighbors in the lattice. Otherwise, $C_{ij}=0$. Diagonal elements are $C_{ii}=12$. One can see that the matrix $\hat{C}$ can be considered as a dynamical matrix of a perfect crystal.

The distribution of eigenvalues of the correlation matrix $\hat{C}$ can be written as

[TABLE]

where $J_{0}(\xi)$ is the Bessel function of zero order.

If we take into account the pairwise correlations only (Eq. (1)), the VDOS can be obtained from the following system of equations [10, 11]:

[TABLE]

Equation (3) defines a so-called critical horizon, that is a contour in the complex plane $Z=X+iY$. Equations (4) and (5) yield an implicit dependency between $g(\omega)$ and $\omega$ from the critical horizon. Figure 1 shows the VDOS of the ensemble $\hat{M}=\hat{B}\hat{B}^{T}$. Solid lines show a numerical solution of Eqs. (2)–(5). One can see a crossover between the low-frequency Debye law $g(\omega)\sim\omega^{2}$ and approximately constant VDOS at high frequencies. The crossover frequency $\omega_{c}$ depends on the value of the parameter $\kappa$.

In Fig. 1 dotted lines show the result of the Kernel Polynomial Method (KPM) [15, 16] applied to a randomly generated matrix $\hat{B}$ for a system with $N=400^{3}$ atoms. It has a qualitative agreement with the solution of Eqs. (2)–(5) for the same values of $\kappa$. The quantitative difference can be caused by neglecting high-order correlations between elements of the matrix $\hat{B}$.

4 Low-frequency asymptotics

Let us consider the low-frequency asymptotics of the VDOS $g(\omega)$ for $\kappa\ll 1$. One can show that this region is defined by the region $X\ll 1$ and $Y\ll 1$ of the critical horizon (3). In this case we can assume that $\rho(\varepsilon)=\sqrt{\varepsilon/128\pi^{4}}+{\cal O}(\varepsilon^{3/2})$ for $\varepsilon\ll 1$. In this case, we can write Eqs. (3) and (4) in the following form:

[TABLE]

where

[TABLE]

Below we consider two cases: $Y\ll X\ll 1$ and $X\sim Y\ll 1$. As we will see, the first case is the left part of the boson peak, including low-frequency Debye behavior. The second case is the right part of the boson peak, which includes the maximum of $g(\omega)/\omega^{2}$.

4.1 Left part of the boson peak

In the case $Y\ll X\ll 1$, we neglect $Y^{2}$ compared to $X^{2}$ in Eqs. (6), (7) and (5) and obtain the following system of equations:

[TABLE]

From the above equations, we find the VDOS which corresponds to the left part of the boson peak:

[TABLE]

where $\omega_{c}=\kappa/(2\sqrt{I_{1}})$ is the crossover frequency. One can show that $|g_{l}(\omega)-g(\omega)|\ll g(\omega)$ in the region $\omega<\omega_{c}-\delta$, where $\delta\sim\kappa^{3/2}$.

In the low-frequencies region $\omega\ll\omega_{c}$, the VDOS (12) has the Debye behavior:

[TABLE]

4.2 Right part of the boson peak

In the case $X\sim Y\ll 1$, the critical horizon (6) is approximately vertical line $X\approx X_{c}=\kappa/(2\sqrt{I_{1}})$. Therefore, Eqs. (7) and (5) can be written as

[TABLE]

From the above equations, we find the VDOS, which corresponds to the right part of the boson peak:

[TABLE]

One can show that $|g_{r}(\omega)-g(\omega)|\ll g(\omega)$ in the region $\omega>\omega_{c}+\delta$, where $\delta\sim\kappa^{3/2}$. As a result, the asymptotics $g_{l}(\omega)$ and $g_{r}(\omega)$ almost coincide with $g(\omega)$ in the full low-frequency region $\omega\ll 1$ except a narrow crossover region $\omega_{c}-\delta<\omega<\omega_{c}+\delta$ (Fig. 2a).

As can be seen in Fig.2b, there exists a boson peak in the reduced VDOS $g(\omega)/g_{D}(\omega)$. For small parameters $\kappa$, two asymptotic expressions (12) and (16) have a good agreement with the precise solution of Eqs. (2)–(5). In the narrow region $\omega_{c}-\delta<\omega<\omega_{c}+\delta$ there is a crossover between $g_{l}(\omega)$ and $g_{r}(\omega)$.

5 Young modulus and boson peak frequency

It was found that the boson-peak frequency $\omega_{b}$ (frequency, where $g(\omega)/g_{D}(\omega)$ has a maximum) depends on the Young modulus $E$ of the system as $\omega_{b}\sim E$ [12, 13, 14]. We can find the relation between $\omega_{b}$ and $E$ from the obtained asymptotics (12) and (16).

In the simple cubic lattice with unit lattice constant $a_{0}=1$ and unit masses $m_{i}=1$, the Debye’s law has a form

[TABLE]

Comparing Eqs. (13) and (17), we find that

[TABLE]

As it follows from this form, if the matrix $\hat{B}$ is square, Young modulus $E$ is equal to zero. Indeed, according to the Maxwell rule, a system has no macroscopical rigidity when the number of bonds is equal to the number of degrees of freedom. The result for the Young modulus (18) is valid only for small $\kappa$, i.e. sides of matrix $\hat{B}$ are not much different. In other words, the system should be soft enough, but have a non-zero stiffness, which can be calculated from (18) in the random matrix model.

Figure 2 shows that the boson-peak frequency $\omega_{b}$ lies on the right asymptotic $g_{r}(\omega)$. In order to find expressions for $\omega_{b}$, we find the derivative of $g_{r}(\omega)/g_{D}(\omega)$ with respect to $\omega$ and equate it to zero. As a result, we obtain the boson-peak frequency

[TABLE]

From Eqs. (18) and (19) we find the linear relation between the boson peak on the Young modulus:

[TABLE]

6 Conclusion

Summarizing, in the present paper we derived closed analytical expressions to find the VDOS $g(\omega)$ of the random matrix model of an amorphous system. The solution of these equations for $g(\omega)$ agrees with the result of the direct calculation of the VDOS by the KPM. Also, from these equations, we found the crossover, which determines the shape of the boson peak. The boson peak frequency $\omega_{b}$ is close to the crossover frequency $\omega_{c}=\sqrt{2/3}\,\omega_{b}$. The VDOS to the left and to the right of $\omega_{c}$ are described by the asymptotic expressions, $g_{l}(\omega)$ and $g_{r}(\omega)$, respectively. The low-frequency behavior of $g_{l}(\omega)$ exhibits the Debye law $g(\omega)\sim\omega^{2}$ and gives the Young modulus $E\sim\kappa\sim\omega_{b}$. The linear relation between the Young modulus and the boson peak has been found in different systems [12, 13, 14] but was not explained so far.

\ack

The authors are thankful to V. I. Kozub for his constructive suggestions and comments. This work is supported by the Russian Federation President Grant no. MK-3052.2019.2.