Relation of vibrational excitations and thermal conductivity to elastic heterogeneities in disordered solids

In crystals, molecules thermally vibrate around the periodic lattice sites. Vibrational motions are well understood in terms of phonons, which carry heat and control heat transport. The situation is notably different in disordered solids, where vibrational excitations are not phonons and can be even localized. Recent numerical work has established the concept of elastic heterogeneity: disordered solids show inhomogeneous local mechanical response. Clearly, the heterogeneous nature of elastic properties strongly influences vibrational and thermal properties, and it is expected to be the origin of anomalous features, including boson peak, vibrational localization, and temperature dependence of thermal conductivity. These are all crucial long-standing problems in materials physics, which we address in the present work. We have considered a toy model able to stabilize different states of matter, by introducing an increasing amount of size disorder. The phase diagram generated by molecular dynamics simulations encompasses the perfect crystalline state with a spatially homogeneous elastic moduli distribution, multiple defective phases with increasing moduli heterogeneities, and eventually a series of amorphous states. We have established clear correlations among the heterogeneous local mechanical response, vibrational states, and thermal conductivity. We provide evidence that elastic heterogeneity controls both vibrational and thermal properties, and is a key concept to understand the anomalous puzzling features of disordered solids.

Figure 1

(a) Mean-squared displacement $〈\mathrm{\Delta }{r}^2〉$ and (b) $Q_6$ vs the disorder parameter $\lambda$, at $T={10}^{-2}$. The system is initialized in the perfect FCC-crystal state $\lambda =1$, where $〈\mathrm{\Delta }{r}^2〉\simeq 0$ and $Q_6\simeq 0.57$. Next, $\lambda$ is decreased from 1 to 0.7 in the fully developed amorphous state, as discussed in the text. Data for two independent instances of the size disorder are shown by open circles and triangles. Both samples undergo the amorphization transition at $\lambda ={\lambda }^{*}\simeq 0.81$, indicated by the vertical line. We also show in (b) by closed circles data for $Q_6$ obtained by following the reverse path, increasing $\lambda$ from 0.7 to 1. In this case, $Q_6$ shows no significant changes for $\lambda ={\lambda }^{*}$, keeping the value pertaining to the amorphous state. This hysteresis effect is discussed in the text.

Figure 2

Schematic illustration of (a) bulk, (b) pure shear, and (c) simple shear deformations. Bulk ($K^m$), pure shear ($G_p^m$), and simple shear ($G_s^m$) local moduli corresponding to these deformations can be determined as discussed in the text.

Figure 3

$\lambda$ dependence of the macroscopic (average) values of (a) bulk $K$, (b) pure shear $G_p$, and (c) simple shear $G_s$ moduli, together with the corresponding standard deviations, (d) $\delta K$, (e) $\delta G_p$, and (f) $\delta G_s$. Circles, squares, and triangles indicate the values of the total modulus $C^m$, the affine term $C_A^m$, and the nonaffine term $C_N^m$, respectively. The vertical lines indicate the transition point ${\lambda }^{*}$, where $G_p$ vanishes and $\delta G_p\simeq 5$. We show with filled circles the data obtained by increasing $\lambda$ from 0.7 to 1 (only total values are shown). In this case, the local moduli distributions are insensitive to the size disorder $\lambda$ and show no significant variations in both average values and standard deviations.

Figure 4

$\lambda$ dependence of the macroscopic (average) values of (a) pressure $p$, pure shear stress ${\sigma }_p$, simple shear stress ${\sigma }_s$; (b) mass density $\rho$, volume fraction $\varphi$; and (c) bond-order $Q_6$ and centrosymmetry $CS$ order parameters. We also plot the corresponding standard deviations in (d) $\delta p, \delta {\sigma }_p, \delta {\sigma }_s$, (e) $\delta \rho , \delta \varphi$, and (f) $\delta Q_6, \delta CS$. The vertical lines indicate the transition point ${\lambda }^{*}$. Note that the density is kept constant $\rho =1.015$ at all values of $\lambda$, while the volume fraction mildly varies in the range of $\varphi =53\%$ to $56\%$.

Figure 5

$\lambda$ dependence of the correlation parameters ${\mathrm{\Psi }}_{C^m{X}^m}$ between local moduli ($C^m$ total modulus, $C_A^m$ affine term only), (a) $K^m$, (b) $G_p^m$, (c) $G_s^m$, (d) $K_A^m$, (e) $G_{pA}^m$, (f) $G_{sA}^m$, and local quantities, $X^m=p^m, {\rho }^m, {\varphi }^m, {\sigma }_p^m, {\sigma }_s^m, Q_6^m, C{S}^m$. The vertical lines indicate the transition point ${\lambda }^{*}$. A detailed discussion of these data is included in the main text.

Figure 6

Vibrational densities of states, $g\left(\omega \right)$, at the indicated values of $\lambda$. These data are discussed in depth in the main text. The vertical line indicates $\omega =5$ for reference. The inset shows the details of the low-$\omega$ ($\omega <5$) region on double-logarithmic scale. We show with lines $g\left(\omega \right)\propto {\omega }^2$ for the case of the perfect crystal at $\lambda =1$, and $g\left(\omega \right)\propto {\omega }^{3/2}$ at the amorphization transition point $\lambda \simeq {\lambda }^{*}$.

Figure 7

Participation ratios, ${\mathcal{P}}_k$, versus the eigenfrequencies ${\omega }_k$ at the indicated values of $\lambda$. The solid lines represent the averaged values $〈{\mathcal{P}}_k〉$ calculated by smoothing the data in bins of width $\mathrm{\Delta }\omega =0.015$. A detailed discussion of these data is given in the main text.

Figure 8

(a) Reduced density of state $\tilde{g}\left(\omega \right)=g\left(\omega \right)/{\omega }^2$ at the indicated values of $\lambda$. (b) $\widehat{g}\left(\omega \right)=g\left(\omega \right)/g_D\left(\omega \right)$ plotted vs $\omega /{\omega }_D$, where $g_D\left(\omega \right)$ is the Debye-model prediction, and ${\omega }_D$ the Debye frequency. For $\lambda \le 0.78$, we observe the boson peak, at the frequency ${\mathrm{\Omega }}^{\text{BP}}\sim 1$, typical of many studies on Lennard-Jones type of glasses [18, 20].

Figure 9

The Debye frequency ${\omega }_D$ (left axis) and the Debye level $3/{\omega }_D^3$ (right axis), plotted as functions of $\lambda$. The vertical line indicates the transition point ${\lambda }^{*}$, where ${\omega }_D$ and $3/{\omega }_D^3$ assume minimum and maximum values, respectively.

Figure 10

Autocorrelation function of the mode energy fluctuations, $C_{E_k}\left(t\right)$, for (a) $\lambda =1$, with ${\omega }_k=4.93$ and ${\tau }_k=34.6$ and (b) $\lambda =0.7$, with ${\omega }_k=3.27$ and ${\tau }_k=32.4$. Correlation functions for total $\delta E_k\left(t\right)$, potential $\delta E_k^{\text{P}}\left(t\right)$, and kinetic $\delta E_k^{\text{K}}\left(t\right)$ energies are shown. Total energy $\left[C_{E_k}\left(t\right)\right]$ shows an exponential decay with a relaxation time ${\tau }_k/2$, where ${\tau }_k$ is the lifetime of mode $k$. Potential $\left[C_{E_k^{\text{P}}}\left(t\right)\right]$ and kinetic $\left[C_{E_k^{\text{K}}}\left(t\right)\right]$ energies exhibit a dumped oscillating decay of frequency $2{\omega }_k$, where ${\omega }_k$ is the mode frequency.

Figure 11

Schematic illustration of the considered acoustic plane waves in the three directions $\left(100\right), \left(110\right)$, and $\left(111\right)$. The wave vector $\mathbf{q}$, and the longitudinal and transverse ($X=L,T$) polarization vectors ${\mathbf{p}}_X$ are also shown [1, 2].

Figure 12

(Left) Lifetimes of the normal modes ${\tau }_k$, and of the acousticlike waves ${\tau }_{X\mathbf{q}}$, plotted as functions of ${\omega }_k$ and ${\omega }_{X\mathbf{q}}$, respectively, at the indicated values of $\lambda$. Here we consider the longitudinal ($L$) and transverse ($T$) acoustic waves propagating in the $\left(100\right), \left(110\right)$, and $\left(111\right)$ directions (see Fig. 11). Specifications of the different data sets are shown in the key. The vertical lines indicate $\omega =5$ for reference. In the perfect crystalline state $\lambda =1, {\tau }_k\simeq {\tau }_{X\mathbf{q}}$, while in defective and disordered states, lifetimes of acoustic-like modes strongly deviate from those pertaining to normal modes. For $\lambda \le 0.86$, the ${\tau }_{X\mathbf{q}}$ at high frequencies are of the order of the minimum time scale set by the Einstein period ${\tau }_E$. This is the typical time scale of the thermal motion of particles, and was estimated as ${\tau }_{\text{E}}\sim \mathcal{O}\left({10}^{-1}\right)$ for the present soft-sphere system [101, 102]. (Right) The averaged lifetimes of the normal vibrational modes $〈{\tau }_k〉$ (j), and of acoustic plane waves $〈{\tau }_{X\mathbf{q}}〉$ (k), versus the frequency ${\omega }_k$ and ${\omega }_{X\mathbf{q}}$, respectively, at the indicated values of $\lambda$. These data are the same as those shown in the left panels, averaged in bins of width $\mathrm{\Delta }\omega =1$, and over all considered propagation directions and polarizations. A comprehensive discussion of these data is included in the main text.

Figure 13

Parametric plot of the lifetimes ${\tau }_k$ and the participation ratios ${\mathcal{P}}_k$ for normal modes at the indicated values of $\lambda$. The solid lines represent the averaged values $〈{\tau }_k〉$ calculated by smoothing the data in bins of width $\mathrm{\Delta }{\mathcal{P}}_k=0.02$. Although modes with larger ${\mathcal{P}}_k$ indeed tend to show larger ${\tau }_k$, the overall correlations are quite weak.

Figure 14

(a) The lifetimes ${\tau }_{\text{ac}}, {\tau }_{\text{ac}}^h, {\tau }_{\text{ac}}^l$ of the acoustic waves, averaged over the entire, high- ($\omega >5$) and low- ($\omega <5$) frequency ranges, respectively, as detailed in the text. In the main panel, the data of ${\tau }_{\text{ac}}$ and ${\tau }_{\text{ac}}^h$ are plotted as functions of the extent of the elastic heterogeneities, $\delta K+\delta G_p+\delta G_s$, for $\lambda \ge {\lambda }^{*}$ (open symbols) and $\lambda <{\lambda }^{*}$ (closed symbols). The line is an exponential fit ${\tau }_{\text{ac}}\simeq {\tau }_{\text{ac}}^h\sim exp[-(\delta K+\delta G_p+\delta G_s)/g_{\tau }]$, with $g_{\tau }\simeq 0.4$. In the inset, we plot ${\tau }_{\text{ac}}^l$ as a function of $\delta G=\delta G_p$ for $\lambda \ge {\lambda }^{*}$ (open symbols) and $\delta G=\delta G_p+\delta G_s$ for $\lambda <{\lambda }^{*}$ (closed symbols). The solid line is a fit of the form ${\tau }_{\text{ac}}^l\sim exp[-\delta G/g_{{\tau }^l}]$ with $g_{{\tau }^l}\simeq 0.5$. (b) The thermal conductivity $\kappa$ shown as a function of the extent of the elastic heterogeneity, together with an exponential fit $\kappa \sim exp[-(\delta K+\delta G_p+\delta G_s)/g_{\kappa }]$, with $g_{\kappa }\simeq g_{\tau }\simeq 0.4$. Note that in the highly disordered states with large values of $\delta K+\delta G_p+\delta G_s$, both ${\tau }_{\text{ac}}$ and $\kappa$ reach the minimum allowed values, where the lifetime ${\tau }_{\text{ac}}\sim \mathcal{O}\left({10}^{-1}\right)$ is the Einstein period [101, 102], and the mean-free-path of the acoustic waves is of the order of the particles diameter.

Figure 15

(a) $\lambda$ dependence of the thermal conductivity at $T={10}^{-2}$ (open circles). We also show (filled circles) the data for the case where $\lambda$ is increased from 0.7 to 1 which show no significant variations with $\lambda$. (The system keeps the amorphous state in this case, see Fig. 1.) (b) $T$ dependencies at the indicated values of $\lambda$. The vertical lines indicate the transition point ${\lambda }^{*}$ in (a), and the temperature $T={10}^{-2}$ in (b). For the perfect crystal ($\lambda =1$), $\kappa$ decreases with $T$ as $\kappa \sim T^{-1}$, due to the anharmonic effects. As $\lambda$ decreases, the value of $\kappa$ decreases and saturates to a $\lambda$-independent value for $\lambda \le 0.86$. Note that $\kappa$ also becomes very mildly dependent on $T$ in the same $\lambda$ range. We also show the values of the simple model of Eq. (21) by squares in (a) and lines in (b), which capture the overall variation of the simulation results, as detailed in the text. We recall that the melting temperature is $T_m\simeq 0.6$, the glass transition temperature is $T_g\simeq 0.2$, and observe that in the liquid state, Eq. (21) is certainly not valid. Indeed, for $T>T_m$ our simulation data for $\kappa$ converge to a value independent of $\lambda$ [24], which cannot be accounted for by the model.

Figure 16

(a) $T$ dependence of the anharmonic term, $1/{\tau }_{\text{anh}}$ [Eq. (23)] at $\lambda =1$, and (b) $\lambda$ dependence of the disorder term, $1/{\tau }_{\text{dis}}$ [Eq. (24)] at the indicated values of $T$. The vertical line indicates the melting temperature $T_m\simeq 0.6$ at $\lambda =1$ in (a), and the transition point ${\lambda }^{*}$ in (b). In (b) we also plot (closed black circles) the values extracted from the lifetimes ${\tau }_{\text{ac}}$ shown in Fig. 14 by using Eq. (25). We can see that $1/{\tau }_{\text{dis}}$ does not depend on temperature for $T\le {10}^{-1}$, and there is a good agreement with the values extracted from ${\tau }_{\text{ac}}$. This observation supports the validity of the additive decomposition of the attenuation of Eq. (22).

Figure 17

Temperature dependence of (a) $C_{\text{QM}}/T^3$ and (b) ${\kappa }_{\text{QM}}$ at the indicated values of $\lambda$. The specific heat $C_{\text{QM}}$ [Eq. (27)] is calculated from the $g\left(\omega \right)$ data of Fig. 6. The thermal conductivity ${\kappa }_{\text{QM}}$ is obtained by applying Eq. (28) to the MD data of Fig. 15. For those calculations, we have used physical Argon units, $\sigma =3.405\AA , \epsilon /k_B=125.2$ K, and $\tau =2.11$ ps. In the figures, the temperature, specific heat, and thermal conductivity are measured in units of K, $k_B$, and $\text{W}{\text{m}}^{-1}{\text{K}}^{-1}$, respectively. For $\lambda \le 0.86, {\kappa }_{\text{QM}}$ increases as ${\kappa }_{\text{real}}\propto T^{\gamma }$ with $\gamma =1.8\simeq 2$, and a plateau value is observed around $T\sim 20$ K. The excess peak in $C_{\text{QM}}/T^3$, however, is found in the lower temperature region $T\le 1$ K. A discussion of this point is included in the text.