Thermal conductance of one dimensional disordered harmonic chains

TL;DR

This paper investigates how disorder affects heat conduction in one-dimensional harmonic chains, revealing universal scaling laws and conditions under which Fourier's law is satisfied.

Contribution

It introduces a universal scaling parameter linking disorder, temperature, and localization length, providing new insights into phonon-mediated heat conduction in disordered systems.

Findings

Thermal conductance scales non-trivially with system size.

Strong disorder with heavy-tailed distributions can satisfy Fourier's law.

A universal behavior of conductance emerges when considering the scaling parameter.

Abstract

We study heat conduction mediated by longitudinal phonons in one dimensional disordered harmonic chains. Using scaling properties of the phonon density of states and localization in disordered systems, we find non-trivial scaling of the thermal conductance with the system size. Our findings are corroborated by extensive numerical analysis. We show that a system with strong disorder, characterized by a `heavy-tailed' probability distribution, and with large impedance mismatch between the bath and the system satisfies Fourier's law. We identify a dimensionless scaling parameter, related to the temperature scale and the localization length of the phonons, through which the thermal conductance for different models of disorder and different temperatures follows a universal behavior.

Tables

Table 1. Table 1: Summary of the scaling behavior for thermal conductance, G 𝐺 G (from this work), density of states (DOS), ρ 𝜌 \rho and localization length, ξ 𝜉 \xi (from ref  Ziman ( 1982 ) ) under different impedance mismatch (coupling k ) k) and disorder strengths.

Disorder	Localization	DOS	Impedance Mismatch	Conductance
Uniform	$\xi \sim {\omega }^{-2}$	$\rho \sim 1$	Any	$G\sim L^{-\frac{1}{2}}$
Power law $ϵ\ge 2$ (weak)	$\xi \sim {\omega }^{-2}$	$\rho \sim 1$	Any	$G\sim L^{-\frac{1}{2}}$
Power law, $1\le ϵ\le 2$	$\xi \sim {\omega }^{-ϵ}$	$\rho \sim 1$	Any	$G\sim L^{-\frac{1}{ϵ}}$
Power law $ϵ\le 1$ (strong)	$\xi \sim {\omega }^{-\frac{2ϵ}{ϵ+1}}$	$\rho \sim {\omega }^{\frac{ϵ-1}{ϵ+1}}$	Low	$G\sim L^{-\frac{1+ϵ}{2ϵ}}$
			High	$G\sim L^{-1}$

Equations

$$G(L,T)\approx {\int }_0^{\infty }\frac{\mathrm{d}\omega }{(2\pi {)}^2}h\omega \frac{\partial f_T(\omega )}{\partial T}\tau (\omega ).$$

$$G(L,T)=\frac{3g_{qm}}{{\pi }^2}{\int }_0^{\infty }\mathrm{d}x\frac{x^2{\mathrm{e}}^x}{{\left({\mathrm{e}}^x-1\right)}^2}\tau (x{\omega }_T).$$

$$G_{qm}(L,T)\approx \frac{3}{{\pi }^2}\left(\frac{{\omega }_L}{{\omega }_T}\right)+\mathcal{O}{\left(\frac{{\omega }_L}{{\omega }_T}\right)}^2.$$

$$G(L)\approx \Sigma {\int }_0^{{\omega }_L}\mathrm{d}\omega \rho (\omega ),$$

$$G$$

$$M\frac{{\mathrm{d}}^2{u}_n}{\mathrm{d}{t}^2}=K_n(u_{n+1}-u_n)-K_{n-1}(u_n-u_{n-1})$$

$$-M{\omega }^2{x}_n=-K_{n-1}(x_n-x_{n-1})+K_n(x_{n+1}-x_n)$$

$$M{\omega }^2\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 0& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_N\end{array}\right]$$

$$=\left[\begin{array}{c}-k{x}_0\\ 0\\ ⋮\\ 0\\ -k{x}_{N+1}\end{array}\right]\equiv \left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_{N-1}\\ b_N\end{array}\right]$$

$$\text{or, }\left(M{\omega }^{2}I+A\right)|X⟩$$

$$\text{or, }|X⟩$$

$$|X⟩=\sum _{i=1}^N\frac{1}{M{\omega }^2+{\lambda }_i}|V_i⟩⟨V_i|B⟩$$

$$|X⟩\approx \frac{1}{\delta {\lambda }_i}|V_i⟩⟨V_i|B⟩$$

$$|X⟩\approx \frac{1}{\delta {\lambda }_i}|V_i⟩⟨V_i|B⟩+\frac{1}{\delta {\lambda }_{N-1}}|V_{N-1}⟩⟨V_{N-1}|B⟩+\frac{1}{\delta {\lambda }_N}|V_N⟩⟨V_N|B⟩.$$

$$x_1$$

$$\chi ={\cos }^{-1}\left(1-\frac{{\omega }^2}{2{\omega }_0^2}\right).$$

$$t$$

$$\tau =|t{|}^2$$

$$=\frac{4k^2(v_1{v}_N{)}^2{\chi }^2}{\left[{\chi }^2+(\frac{k}{\delta {\lambda }_N}+1-{\chi }^2{)}^2\right]\left[\delta {\lambda }^2{\chi }^2+{\left\{k\left(\frac{\delta \lambda }{\delta {\lambda }_N}+(v_1^2+v_N^2)\right)-\delta \lambda ({\chi }^2-1)\right\}}^2\right]}.$$

$$\frac{k}{\delta {\lambda }_N}=\frac{k}{{\lambda }_N+M{\omega }^2}\approx \frac{k}{{\lambda }_N}\approx -1.$$

$${\chi }^2+{\left(\frac{k}{\delta {\lambda }_N}+1-{\chi }^2\right)}^2$$

$$\delta {\lambda }^2{\chi }^2+{\left\{k\left(\frac{\delta \lambda }{\delta {\lambda }_N}+(v_1^2+v_N^2)\right)-\delta \lambda ({\chi }^2-1)\right\}}^2$$

$$\tau$$

$$\approx \frac{4k^2(v_1{v}_N{)}^2/{\chi }^2}{\left[(\delta \lambda -s{)}^2+\frac{s^2}{{\chi }^2}\right]}.$$

$$\tau (\omega )$$

$$\gamma =\frac{2\pi (v_1{v}_N{)}^{2}k}{M\tilde{\omega }(v_1^2+v_N^2)\chi }=\frac{2\pi (v_1{v}_N{)}^{2}k}{M{\tilde{\omega }}^2(v_1^2+v_N^2)}\sqrt{\frac{k}{m}},$$

$$\sigma =\frac{k(v_1^2+v_N^2)}{M\tilde{\omega }\chi }=\frac{k(v_1^2+v_N^2)}{M{\tilde{\omega }}^2}\sqrt{\frac{k}{m}}.$$

$${\gamma }_s=\frac{2\pi (v_1{v}_N{)}^{2}k\chi }{M\tilde{\omega }(v_1^2+v_N^2)}=\frac{2\pi (v_1{v}_N{)}^2}{M(v_1^2+v_N^2)}\sqrt{km},$$

$${\sigma }_s=\frac{k\chi (v_1^2+v_N^2)}{M\tilde{\omega }}=\frac{(v_1^2+v_N^2)}{M}\sqrt{km}.$$

$$A$$

$$A^b$$

$$A^d$$

$$\tilde{A}$$

$$⟨1⟩V_n^{(1)}\approx \sum _i\frac{⟨V_i^{(0)}⟩A^d{V}_n^{(0)}}{{\lambda }_n^{(0)}-{\lambda }_i^{(0)}}⟨1⟩V_i^{(0)}$$

$$⟨1⟩V_n\approx ⟨1⟩V_n^{(1)}$$

Full text

Introduction- The study of heat transport via phonons in low dimensional (spatial dimension $d<3$) classical and quantum mechanical systems has attracted considerable theoretical and experimental attention in recent years Tamaki et al. (2017); Cividini et al. (2017); Xu (2016); Dhar (2008); Bonetto et al. (2000); Li and Wang (2003); Monasterio et al. (2005); Narayan and Ramaswamy (2002); Mai and Narayan (2006); Dhar (2001); Moghaddasi Fereidani and Segal (2019). One of the main objectives of these studies is to understand the scaling of heat flux $J$ which, according to Fourier’s law Dhar (2008), should scale with the system size $L$ as $J\propto L^{-1}$ ($L$ is measured along the direction of heat propagation). But extensive numerical and analytical studies in the past few decades have revealed the possible violation of Fourier’s law for low dimensional systems Casher and Lebowitz (1971); Rubin and Greer (1971); Casati et al. (1984); Dhar (2001); Lepri et al. (2003); Hurtado and Garrido (2016). These studies show that $J\propto L^{{(\gamma-1)}}$ with $\gamma\neq 0$ which in turn implies $L$-dependent thermal conductivity, $\kappa=\lim\limits_{L\to\infty}\lim\limits_{\Delta T\to 0}\frac{JL}{\Delta T}\propto L^{\gamma}$ ($\Delta T$ being the temperature difference across the system) Dhar (2008). The violation of Fourier’s law in low-dimensional systems is also observed experimentally in the case of carbon nanotubes Chang et al. (2008), nanowires Yang et al. (2010) and graphene Xu et al. (2014).

For systems of finite size, instead of thermal conductivity $\kappa$, it is useful to study thermal conductance, $G=\kappa L^{d-2}$. Thus, for one dimensional systems $(d=1)$, according to Fourier’s law, we expect $G(L)\propto L^{-\beta}$, with $\beta=1$ for normal heat transport, while $\beta\neq 1$ implies anomalous heat transport. Note that the scaling exponents $\gamma$ and $\beta$ characterizing the thermal conductivity and thermal conductance, respectively, are related by $\beta=1-\gamma$. One interesting question is under what conditions $\beta=1$ (or, $\gamma=0$)?

Various aspects, such as disorder Casher and Lebowitz (1971); Herrera-González et al. (2015); Amir et al. (2018), phonon-phonon interaction Dhar and Lebowitz (2008); Hu et al. (1998), presence of pinning potential Roy and Dhar (2008); Bodyfelt et al. (2013), nature of the heat baths Dhar (2001) and the coupling between the system and the heat bath Amir et al. (2018), have been shown to affect heat transport. Particularly, theoretical studies for one dimensional isotopically (mass) disordered harmonic chains show that $J\propto L^{-\frac{1}{2}}$ with free boundary conditions Rubin and Greer (1971) while $J\propto L^{-\frac{3}{2}}$ with fixed boundary condition Casher and Lebowitz (1971), implying that $\beta$ can be $\frac{1}{2}$ or $\frac{3}{2}$. For this particular model, it was also shown that normal scaling (i.e. $\beta=1$) can be observed only under specific choices of the thermal bath Dhar (2001). It was also argued, under free boundary conditions, that one-dimensional harmonic chains with spatially correlated disorder may exhibit normal heat conduction asymptotically Herrera-González et al. (2015).

Can one have normal heat transport in one dimensional disordered (uncorrelated) harmonic chains even within free boundary condition? A recent theoretical study Amir et al. (2018) predicts that for a weakly coupled disordered harmonic chain one may observe normal heat transport in the presence of strong disorder, when disorder is characterized by a heavy-tailed distribution. While it is important to verify this theoretical prediction, it is equally interesting to ask: How does thermal conductance scale with $L$ if the coupling between the system and the heat bath is not weak? For a given coupling, how does $\beta$ depend on the nature of the disorder? In this Letter, we address these questions by studying, analytically as well as numerically, the scaling of thermal conductance in one dimensional disordered harmonic chains for different types of disorder and coupling between the system and the heat bath.

Heat conduction by phonons is similar to electrical conduction, but with a crucial difference: the presence of a localization threshold at zero frequency. This leads to a diverging localization length, $\xi(\omega)$, for $\omega\to 0$ Azbel (1983); Ishii (1973) and has strong consequences on the scaling of thermal conductance. Specifically, for a given $L$ and disorder strength, one can define a cut-off frequency $\omega_{L}$, for which $\xi(\omega_{L})=L$. All phonons with $\omega\leq\omega_{L}$ are effectively delocalized, i.e. $\xi(\omega)\geq L$, and contribute to the heat transport.

*Model & background.- *We consider a one dimensional disordered system consisting of $N$ particles, each of mass $M$, connected by harmonic springs with spring constants ${K_{i}}(i=1,2,\cdots N-1)$, chosen randomly from a given distribution, cf. Fig. 1. $K_{i}$ is the spring constant of the spring connecting particles $i$ and $i+1$ in the disordered chain. The two ends $(i=1$ and $i=N)$ of the chain of length $L$ are coupled to two heat baths at temperatures $T_{1}$ (left bath) and $T_{2}(<T_{1}$; right bath), respectively. Here, $L=(N-1)r_{0}$ with $r_{0}$ being the average interparticle distance. Heat baths are modeled as ordered harmonic chains consisting of an infinite number of equal masses $(m)$, and connected by identical springs $(k)$. The system is coupled to two heat baths via two springs each having spring constant $k$. If $k$ is much smaller (larger) compared to the typical spring constant in the disordered chain, we refer to the system as weakly (strongly) coupled to the reservoir. Note that our setup corresponds to the case of ‘free boundary condition’ considered in the literature Rubin and Greer (1971); Dhar (2001). Below we work in units where the mass $M$ of the system’s particles, the natural frequency of the bath $\omega_{0}=\sqrt{k/m}$, and Boltzmann’s constant are all set to unity. We express the stiffness of the springs in units of $M\omega_{0}^{2}=1$.

In the current study, we consider two models of disorder: (1) Uniform distribution: $K_{i}=(1+R_{i})$ where $R_{i}$ follows a uniform distribution of width $W$, i.e., $R_{i}\in\left[-W/2,W/2\right]$. Large values of $W$ correspond to stronger disorder and $W=2$ is the strongest possible disorder strength. (2) Power-law distribution: $K_{i}$ follow a power-law probability distribution, $P(K)\propto K^{\epsilon-1}$, where $0<K\leq 1$ and disorder strength is quantified by the dimensionless parameter $\epsilon(>0)$ FN2 (a). This situation arises naturally if $K$ decays exponentially with interparticle separation which follows a Poisson process Amir et al. (2013, 2010). Small $\epsilon(<1)$ corresponds to strong disorder.

As noted above, transport is mediated by effectively delocalized low-frequency phonons. Thus, it will be crucial to understand the scaling behavior of the localization length, $\xi(\omega)$ and density of states (DOS), $\rho(\omega)$, in the limit $\omega\to 0$. Here, we briefly summarize earlier theoretical predictions related to these scaling behaviors. With uniform disorder, $\rho(\omega)$ approaches a constant (equivalent to the Debye scaling for 1D ordered systems Ashcroft and Mermin (1976)) while $\xi(\omega)\propto\omega^{-2}$ Alexander et al. (1981); Ziman (1982); Amir et al. (2010); Azbel (1983); Ishii (1973). These results were predicted for weak disorder, but we will show that these hold for any value of $W$. With power law disorder, we define three qualitatively different regimes: in the weak disorder regime $(\epsilon>2)$, $\rho$ and $\xi$ have the same scaling behavior as that of the uniform disorder; In the intermediate disorder regime $(1<\epsilon\leq 2)$, where the variance of the compressibility diverges (but its mean does not), the localization length has a non-trivial scaling $\xi(\omega)\propto\omega^{-\epsilon}$, but the DOS still exhibits Debye scaling, and in the strong disorder regime $(\epsilon\leq 1)$, the mean of the system’s compressibility diverges, and both $\rho$ and $\xi$ feature anomalous scaling Ziman (1982): $\rho(\omega)\propto\omega^{\frac{\epsilon-1}{\epsilon+1}}$, $\xi(\omega)\propto\omega^{-\frac{2\epsilon}{\epsilon+1}}$. These results, which are crucial for the discussion below, are summarized in Table 1 and numerically demonstrated in Fig. 2(a-b) and Fig. 3.

Analytical results - To study the heat transport we follow the Landauer scattering approach in which propagation of a phonon of given frequency $\omega$ through the disordered chain is characterized by a transmission coefficient, $\tau(\omega)$. For a one dimensional system the thermal conductance $G(L,T)$ is Rego and Kirczenow (1998):

[TABLE]

where $f_{T}$ is the Bose-Einstein distribution function and we also assume $\Delta T=T_{1}-T_{2}\ll T=\frac{T_{1}+T_{2}}{2}$.

Note that the system features two competing frequency scales: the disorder-related $\omega_{L}$ and the thermal frequency $\omega_{T}=T/\hbar$. Phonons with $\omega>\omega_{L}$ do not contribute to conductance because they are localized. Phonons with $\omega>\omega_{T}$ do not contribute because they are not populated. Therefore, the integral in Eq. (1) is better represented in terms of the non-dimensional frequency $x=\omega/\omega_{T}$,

[TABLE]

Here $g_{qm}=\pi^{2}T/(3h)$ is the quantum of thermal conductance Rego and Kirczenow (1998); Pendry (1983) which is the universal value of $G(L,T)$ in the limit $T\to 0$. To see this, note that $\tau(\omega)\to 1$ for $\omega\to 0$ due to the existence of a Goldstone mode, related to the translational invariance of the system. Using this fact, it is straightforward to show that for a given $L$ Eq. (2) yields $G(L,T)\to g_{qm}$ for very small $T~{}(\omega_{T}\ll\omega_{L})$, regardless of any other property of $\tau(\omega)$. $g_{qm}$ is thus the natural unit of conductance for our system and below we express all results in these units by defining $G_{qm}(L,T)=G(L,T)/g_{qm}$.

Of course, the limit of $L\to\infty$ and finite $T$ is of more interest, but an exact evaluation of the integral in the general case is not feasible. Nonetheless, much insight can still be gained in some interesting cases.

We first consider the situation where the stiffness of the coupling spring $k$ is comparable to that of the chain, i.e. $k\approx 1$ FN2 (b). In this case there is relatively little impedance mismatch between the chain and the bath, implying less reflectance of the incident phonons from the bath-system boundary. In this situation, phonons get transmitted even when their frequency is not close to an eigenfrequency of the chain. Therefore, from a scaling perspective we can approximate that $\tau(\omega)=1$ for all phonons with $\omega\leq\omega_{L}$ and zero otherwise. A numerical calculation of $\tau(\omega)$, shown in Fig. 2c, demonstrates that this approximation is crude but reasonable. As shown below, it quantitatively captures the scaling behavior.

With this approximation, Eq. (1) depends only on the dimensionless combination $\omega_{L}/\omega_{T}$ (which is the upper integration limit) implying that thermal conductance for $k\approx 1$ should follow a universal curve, independent of temperature and disorder, when expressed in terms of $\omega_{L}/\omega_{T}$. In fact the integral can be carried out in closed form, and for large $L$ (or large $T$ .i.e., $\omega_{T}\gg\omega_{L})$ it reads

[TABLE]

In order to get the explicit dependence on system size in this limit, we use the known scaling $\xi(\omega)\propto\omega^{-\alpha}$, see Amir et al. (2018) and Fig. 3. Straightforward manipulation shows that this implies $G_{qm}\propto L^{-\frac{1}{\alpha}}T^{-1}$, that is, $\beta=\alpha^{-1}$ in this limit of small impedance mismatch and large $L$. Also, for a given $L$, $G(T)\propto T$ for small $T$ and $G(T)\sim$ constant (saturates) for high $T$ sup .

For uniform disorder, theory predicts $\alpha=2$ implying $G(L)\propto L^{-1/2}$, in accord with previously reported results for mass-disordered chain under free boundary condition Rubin and Greer (1971); Dhar (2001). For power-law disorder, as $\alpha$ depends on disorder strength $\epsilon$, $\beta$ also depends on $\epsilon$ with $G\propto L^{-1/2},L^{-1/\epsilon}$ and $L^{-\frac{1+\epsilon}{2\epsilon}}$ in the weak, intermediate and strong disorder regimes, respectively. These results are summarized in Table. 1.

This concludes the case of $k\approx 1$, where transmission is approximately constant for all $\omega$ below a certain cutoff. How does the picture change in the case of strong impedance mismatch, $k\gg 1$ or $k\ll 1$? In this case transmission is negligible for almost all frequencies, except those which are close to an eigenfrequency of the disordered chain. In previous work Amir et al. (2018), it was shown that in the weak coupling limit, $k\ll 1$, $\tau$ has a structure of non-overlapping Lorentzians for phonons with $\omega<\omega_{L}$, cf. Fig. 2c. Each Lorentzian is centered around an eigenfrequency of the disordered chain and the area of each Lorentzian, i.e. it’s integrated contribution to the thermal conductance, was shown to be $\omega$-independent for the delocalized modes Amir et al. (2018). Calculating the integral in general for any $\omega_{L}/\omega_{T}$ is difficult, but if we are only interested in the scaling behavior for large $L$, the integral essentially counts the number of eigenmodes of the disordered chain within the frequency range $0<\omega\leq\omega_{L}$ FN (3):

[TABLE]

where $\Sigma$ is the area of each Lorentzian. Considering $\rho(\omega)=D\omega^{s}$ (see Amir et al. (2018) and Fig. 2a-b), where $D$ depends on disorder, for large $L$ we get

[TABLE]

Thus, in the weak coupling regime $G_{qm}\propto L^{-\frac{s+1}{\alpha}}T^{-1}$ for a fixed $T$. For uniform disorder, as well as power-law disorder with $\epsilon>1$ (i.e. the weak and intermediate regimes), we have $s=0$ and thus the scaling of thermal conductance with $L$ remains the same as in case of an impedance-matched bath $k\approx 1$. Interestingly, for strong disorder, $\epsilon<1$, we have $s=\frac{\epsilon-1}{\epsilon+1}$ and $\alpha=\frac{2\epsilon}{1+\epsilon}$, which together cancel out exactly to yield normal Fourier-like heat conduction $\beta=1$. Also, note that, when expressed in terms of $\omega^{s+1}_{L}/\omega_{T}$, thermal conductance should follow a universal curve in the weak coupling regime and the large $L$ limit.

Lastly, we deal with the case of very large $k$, i.e. the strong coupling regime. A careful analysis, presented fully in the supplementary material sup , shows that this limit is equivalent to a system with the first and last particles excluded, i.e., effectively a system of $(N-2)$ particles. Therefore, like in the case of weak coupling, $\tau(\omega)$ is composed of non-overlapping peaks with an $\omega$-independent area, cf. Fig. 2. Since the density of states and the localization length are independent of the coupling $k$, all our predictions for $k\ll 1$ hold also for $k\gg 1$. That is, the same scaling exponents emerge in the case of strong impedance mismatch, regardless of whether $k$ is very small or very large. All our theoretical predictions for different disorder types and coupling strengths are summarized in Table. 1.

Numerical results.- We test our theoretical predictions by numerically computing thermal conductance and other properties, such as density of states and localization lengths, for all cases considered above. The density of states, $\rho(\omega)$, for different disorder strengths is shown in Fig. 2(a-b) for the uniform and power-law disorder, respectively, with $N=2000$ sup . For uniform disorder, $\rho(\omega)$ approaches a constant as $\omega\rightarrow 0$, that is, $s=0$. For power law disorder, $\rho(\omega)$ diverges with an exponent consistent with the theoretical prediction, $s=\frac{\epsilon-1}{\epsilon+1}$. For weak and intermediate disorder theory predicts $s=0$ but a weak divergence is observed for $\epsilon=1$, the origin of which in not clear to us at present.

To calculate $G$, we compute $\tau(\omega)$ directly for different $k$ and disorder types using a transfer matrix method sup . For a single realization of the disorder, the dependence of $\tau(\omega)$ on $k$ is shown in Fig. 2(c) (with $W=0.50$ and $N=11$).

For a given disorder and a fixed $\omega$, we find that $\tau(\omega)$ decays exponentially with $L$ sup . This defines a length scale which we interpret as the localization length $\xi$, i.e., $\tau(\omega,L)\propto\exp[-L/\xi(\omega)]$ sup . We find that $\xi$ diverges like $\xi(\omega)\propto\omega^{-\alpha}$, consistent with theoretical predictions (see Fig. 3).

Finally, we compute $G$ for different $L$, $T$, disorder types and coupling strengths, using Eq. (1). This is presented in Fig. 4(a) (for $k=1.0$) and 4(b) (for $k=0.01$), which shows that the numerical results for all cases agree with the theoretical predictions (for $k=100$, see supplemental material sup ). Heat transport is anomalous $(\beta\neq 1)$ for all cases except for the power-law disorder in the weak/strong coupling regime and also when $\epsilon=1$ (Fig.4(b)). In addition, panels (c) and (d) demonstrate that when expressed in terms of $\omega_{L}/\omega_{T}$ (for $k=1$) and $\omega_{L}^{s+1}/\omega_{T}$ (for $k=0.01)$, respectively, all data collapse on a single curve, following Eq. (3) and Eq. (5).

*Conclusions and discussion.- *In this Letter we studied analytically the thermal conductance of one-dimensional disordered harmonic chains and corroborated the theoretical predictions with extensive numerical simulations. We found a non-trivial scaling behavior of the conductance, which depends both on the nature of disorder and the coupling between the system and the heat baths. This dependence is mediated by the scaling of the localization length and density of states. In addition, we identified the dimensionless scaling parameter with which one has a unified description for all temperatures and systems sizes.

Specifically, we found that the conductance features anomalous scaling with $L$ for uniform disorder and for weak and intermediate power-law disorder (i.e. with a well-defined mean). Interestingly, for strong power-law disorder and strong impedance mismatch, $k\ll 1$ or $k\gg 1$, normal scaling $G\propto L^{-1}$ is observed. For strong disorder and low impedance mismatch, i.e. $k\approx 1$, the scaling exponent $\beta$ can be greater than unity. While most of the previous works on anomalous heat transport focused on the classical regime, our study encompasses both, quantum and classical regimes. It would be interesting to study in the future the fluctuations in the thermal conductance, and in particular ask if the universal fluctuations observed in electronic systems have a counterpart here.

*Acknowledgements.- *We thank David Mross for valuable discussions. BA acknowledges the computational facilities at the Weizmann Institute of Science and Harvard University. Financial support from the Weizmann Institute of Science is gratefully acknowledged. YBS was supported by the JSMF post-doctoral fellowship for the study of complex systems.

Supplemental material for ‘Thermal conductance for one dimensional disordered harmonic chains’

I Analytical prediction of the transmission coefficient in the regime of strong impedance mismatch

The model of a one dimensional disordered chain with $N$ particles coupled to two heat baths at the two ends is shown in Fig. 1 of the main article. $u_{n}$ denotes the displacement of the $n$-th particle in the disordered chain from its equilibrium position. The equation of motion of the $n$-th particle is

[TABLE]

Considering $u_{n}=x_{n}\exp[-i\omega t]$ (i.e. looking at normal mode solution) in Eq. (S.1), we get

[TABLE]

Note that $K_{0}=k,$ and $K_{N}=k$. We can write the above set of equations in a compact form using the following matrix notation:

[TABLE]

where, $I$ is a unit matrix of order $N\times N$, and $A$ is a tri-diagonal matrix of order $N\times N$ with $A_{i,i}=-K_{i-1}-K_{i},A_{i,i+1}=K_{i}$ and $A_{i,i-1}=K_{i-1}$. Also, $|X\rangle$ and $|B\rangle$ are column vectors and all components of $|B\rangle$ are zero except for the first and last elements; $b_{1}=-kx_{0}$ and $b_{N}=-kx_{N+1}$. We are using the quantum-mechanical bra-ket notation for convenience.

As the normalized eigenvectors $\{|V_{i}\rangle\}$, corresponding to the eigenvalues $\{\lambda_{i}\}(i=1,2,\cdots,N)$ of $A$ form a complete basis, we can expand the right hand side of Eq. (S.3) in terms of $\{|V_{i}\rangle\}$ and write the Eq. (S.3) as

[TABLE]

Whenever there is strong impedance mismatch between the bath and the system (which occurs in case of weak and strong coupling), transmission is appreciable only for frequencies close to the eigenfrequencies of $A$. In this case the sum in Eq. (S.4) is dominated by a small number of summands.

An eigenmode can have a non negligible contribution to the sum of Eq. (S.4) in one of two cases: either

1. (a)

it is close to resonance, i.e. the denominator $M\omega^{2}+\lambda_{i}$ is small, or 2. (b)

the projection of the eigenmode on the chain ends, $\left\langle V_{i}|B\right\rangle$ is large.

As was shown in ref. Amir et al. (2018), in the case of very weak coupling ($k\ll 1$) only case (a) occurs and near an eigenfrequency of the chain Eq. (S.4) is well approximated by

[TABLE]

where, $\delta\lambda_{i}=M(\omega^{2}-\omega_{i}^{2})$ as $\lambda_{i}=-M\omega_{i}^{2}$. Using Eq. (S.5), one can derive an expression for transmission coefficient as described in ref. Amir et al. (2018), which is valid for weak coupling.

The case of very large coupling $(k\gg 1)$ requires more care as both options (a) and (b) occur. In addition to resonant modes, the strong springs in the ends of the chains give rise to two special modes, each localized in one of the ends of the chain. These modes’ contribution to the sum of Eq. (S.4) is not negligible for any frequency (with respect to the contributions from other modes).

To see this, let us order the eigenvalues of $A$ in ascending order of their magnitude: $\lambda_{1}<\lambda_{2}<\cdots<\lambda_{i}<\cdots<\lambda_{N-1}<\lambda_{N}$. In Sec. II we show, within a rigorous degenerate perturbation theory, that for $k\gg 1$ these two special modes are $|V_{N-1}\rangle\approx[1~{}0~{}0~{}\cdots~{}0]^{T}$ and $|V_{N}\rangle\approx[0~{}0~{}0~{}\cdots~{}1]^{T}$, and their associated eigenvalues are $\lambda_{N-1}\approx\lambda_{N}\approx-k$. All other eigenvalues are much smaller.

Thus, in the strong coupling limit, instead of Eq. (S.5), we have

[TABLE]

Writing $\omega_{i}=\bar{\omega}$ and $\delta\lambda_{i}=\delta\lambda$, thus, the displacements $x_{1}$ and $x_{N}$ are given by

[TABLE]

where, $v_{j}\equiv\braket{j}{v}$ is the amplitude of the eigenmode $\bar{\omega}$ at site $j$. Assuming that a phonon of frequency $\omega$ (with unit amplitude) and wave number $q$ coming from the left side of the disordered chain, we have $x_{n}=\mathrm{e}^{in\chi}+r\mathrm{e}^{-in\chi}(n<1)$ and $x_{m}=t\mathrm{e}^{im\chi}(m>N)$, where $r$ and $t$ represent the reflection and transmission amplitudes and $\chi=qa$, $a$ being the lattice constant in the ordered region. Using phonon dispersion relation Ashcroft and Mermin (1976), we have

[TABLE]

Note that $x_{N+1}=t$ and $x_{0}=1+r$. Substituting these expressions for $x_{0}$ and $x_{N+1}$ in Eq. (S.7) and using Eq. (S.2), we can find the transmission amplitude. In order to simplify the algebra, we assume $\lambda_{N}=\lambda_{N-1}\approx-k$, so that $\delta\lambda_{N}=M\omega^{2}+\lambda_{N}=\delta\lambda_{N-1}$. With this, we get the following expression for transmission amplitude:

[TABLE]

For large $N$, as only low frequency modes contribute to the transmission coefficient $\chi\ll 1$ leading to $\chi=qa\approx\omega\sqrt{\frac{m}{k}}$. For very small $\chi$, the transmission coefficient takes the following form

[TABLE]

As $\omega$ is very small and $\lambda_{N}\approx-k$, we can write

[TABLE]

Using this, we have

[TABLE]

where, $s=k(v_{1}^{2}+v_{N}^{2})$. Using above simplifications, finally, we can write

[TABLE]

We define a frequency $\widetilde{\omega}$ such that $\delta\lambda-s=M(\omega^{2}-\widetilde{\omega}^{2})$. As transmission is appreciable for frequencies close to the eigenfrequencies of the chain, we have $\omega\approx\widetilde{\omega}\approx\bar{\omega}$ and consequently $\delta\lambda-s\approx 2M\widetilde{\omega}(\omega-\widetilde{\omega})$. Substituting this into Eq. (S.13) shows that the transmission coefficient is approximately a Lorentzian in terms of $\omega:$

[TABLE]

where the area $(\gamma)$ and width $(\sigma)$ are given by

[TABLE]

For large $k$, we compute $v_{1}$ and $v_{N}$ using perturbation theory and interestingly, find (see sec. II) that for the low-frequency delocalized modes $|v_{1}|\approx|v_{N}|\sim\frac{\omega}{k\sqrt{N}}$ and subsequently $\gamma\sim\frac{2\pi}{\sqrt{km}MN}$, i.e., for the delocalized modes area associated with each Lorentzian is independent of the frequency. Thus, their contribution to the thermal conductance is constant, i.e., independent of the frequency.

We also note that the equations for $\gamma$ and $\sigma$ for large $k$ are slightly different from those obtained for small $k(\ll 1)$ where the area $(\gamma_{s})$ and the width $(\sigma_{s})$ of each Lorentzian are given by (see eq. 15 in ref. Amir et al. (2018)):

[TABLE]

For small $k,\chi$ appears in the numerator of the expressions for the area and the width. In the weak coupling scenario $|v_{1}|=|v_{N}|=\frac{1}{\sqrt{N}}$ for the delocalized modes and thus one finds $\gamma_{s}=\frac{2\pi\sqrt{km}}{MN}$ Amir et al. (2018), i.e., the area of each Lorentzian associated with the delocalized modes are independent of the frequency.

To verify that $\tau(\omega)$ is indeed Lorentzian in $\omega$ with area and width given by Eq. (S.15) and Eq. (S.16) for large $k$, we fit the numerical data for $\tau(\omega)$ (see Fig. 2c in the main article) with a Lorentzian given in Eq. (S.14). As shown in Fig. S.1, we find that the area and the width given by Eq. (S.15) and Eq. (S.16), respectively, agree well with those obtained by fitting the numerical data.

II Characteristics of the eigenmodes and the eigenfrequencies of $A$ in the strong coupling limit

In Sec. I we used two properties of the eigenmodes in the $k\gg 1$ limit:

• •

There are two boundary modes (that is, eigenmodes localized at one of the ends of the system) with eigenvalue $\approx k$ that does not scale with $N$.

• •

For eigenmodes with $\omega\to 0$ the vibration amplitude of the first and last particles is proportional to $\omega$.

The former was used to approximate Eq. (S.6) from Eq. (S.4). The latter was used to show that the contribution of each delocalized mode to the condactuance is $\omega$-independent in the $\omega\to 0$ limit, Eq. (S.16). Here we derive both results within first order degenerate perturbation theory.

To set the ground, we write $A$ as:

[TABLE]

where,

[TABLE]

In the strong coupling limit $(k\gg 1)$ as the typical value of $K_{i}$ is much smaller than $k$, we consider $A^{d}$ as a perturbation and compute the first order correction to the eigenvalues and eigenvectors of the unperturbed matrix $A^{b}$. We see that $A^{b}$ has two eigenspaces (“degenerate bands”): The first consists of two modes, $\ket{1}=[1~{}0~{}0~{}\cdots~{}0]^{T}$ and $\ket{N}=[0~{}0~{}0~{}\cdots~{}1]^{T}$, localized at the first and last particle, respectively. Both have eigenvalue $-k$. The second eigenspace consists of $(N-2)$ modes with $\omega=0$. This eigenspace consists of all vibrations that keep the first and the last particles fixed.

First, we can already conclude that the two boundary modes with eigenvalue $\approx-k$ indeed exist. To zeroth order their eigenvalue is exactly $-k$ and they are entirely localized at the first and last particles. Any corrections to this picture become smaller as $k$ increases and thus the perturbation becomes smaller.

We turn to analyze the modes in the $(N-2)$-fold degenerate band. Here the zeroth order picture is trivial and we need the first order correction in order to see the structure. Degenerate perturbation theory Sakurai (1994) tells us that the relevant eigenmodes of the perturbed system are those which diagonalize the perturbation when projected to the degenerate eigenspace. In our case, we need to find the eigenmodes of the matrix that is obtained by deleting the first and last row and column from $A^{d}$, namely

[TABLE]

This matrix describes a chain of $N-2$ particles, connected to two fixed particles at $i=1$ and $i=N$ with the springs $K_{1}$ and $K_{N-1}$ respectively. The eigenmodes of this matrix, $|V_{n}^{(0)}\rangle$, when padded with zeros at both ends to make them $N$-dimensional, are also eigenmodes of $A$, to zeroth order. Note that $\widetilde{A}$ is identical to $A$ of Eq. (S.19), with the difference that it involves $N-2$ particles instead of $N$, and that the first and last springs are $K_{1}$ and $K_{N-1}$ instead of $k$.

Our goal is to understand how $\Braket{1}{V_{n}}$ and $\braket{N}{V_{n}}$ scale with $N$ and $\omega$. In the present context, degenerate perturbation theory Sakurai (1994) tells us that

[TABLE]

where the summation is performed only on $\Ket{V_{i}}$ outside the degenerate eigenspace of $\Ket{V_{n}}$. In our case there are exactly two modes outside the degenerate subspace – the boundary modes $\ket{1}$ and $\ket{N}$ – and consequently the amplitude of the first particle of the $n$-th eigenmode $\Ket{V_{n}}$ is, to first order,

[TABLE]

That is, to first order, the vibration amplitude of the first particle is proportional to the zeroth order vibration amplitude of the second particle.

In an ordered chain, when all $K_{i}$ are equal, $\Ket{V_{n}^{(0)}}$ is a sinusoidal mode which vanishes at $i=1$ and $i=N$,

[TABLE]

Therefore, for large $N$ we have $\braket{2}{V_{n}}\propto n$ and since for small $\omega$ dispersion is linear, $n$ is proportional to $\omega$. Combining this with Eq. (S.22), we get

[TABLE]

as utilized after Eq. (S.16). An identical argument shows that the vibration amplitude $\braket{N}{V_{n}}$ of the last particle follows the same scaling.

Lastly, we argue that Eq. (S.24) holds for disordered chains as the low-frequency modes are effectively delocalized and follow approximately the same dispersion relation as the ordered chain. This is demonstrated in Fig. S.2(a) where we plot the eigenvectors associated with the lowest four eigenfrequencies of a disordered chain with $N=2000,W=1.95$ and $k=100$. Panel (a) shows that the eigenmodes indeed correspond approximately to those of an ordered chain with fixed boundaries. We also verify that the scaling of $|v_{1}|$ and $|v_{N}|$ with $\omega$ and $N$ follow Eq. (S.24) as shown in Fig. S.3(a-b). Note that the linear $\omega$ dependence of $|v_{1}|$ and $|v_{N}|$ holds only for the low-frequency modes and fails for $\omega\gg\omega_{L}$ as the high frequency modes get strongly affected by disorder and become localized (i.e., does not follow the dispersion relation of an ordered harmonic chain).

It is interesting to note that the scaling relation (S.24), which holds for $k\gg 1$, is very different from what one would have in case of the weak coupling $(k\ll 1)$. In that case the roles of $A^{b}$ and $A^{d}$ in Eq. (S.19) are interchanged since $A^{b}$ is much smaller and is treated as a perturbation to $A^{d}$. Note that $A^{d}$ describes to a chain whose end particles are free (i.e., not connected to “external” springs like $K_{1}$ and $K_{N-1}$ in $\tilde{A}$), as can be seen from the fact that all rows and columns sum to zero. That is, a global uniform translation is an eigenmode (actually, a Goldstone mode) of the system. In this scenario the modes do not vanish at the boundary but rather their gradient does. This is demonstrated in Fig. S.2(b), showing eigenmodes of the same disordered chain as in Fig. S.2(a) but with $k=0.01$.

Lastly, we clarify that these notions of “free” and “fixed” boundary conditions are distinct from the way these notions are used in the literature Casher and Lebowitz (1971); Rubin and Greer (1971); Dhar (2001, 2008). In their jargon our system always has free boundary conditions, irrespective of the coupling strength, since the two end particles of the disordered chain are not pinned, i.e., these particles are not connected to some external springs other than the ones that couple the system with the heat baths. We only use the notions of “free” and “fixed” boundary conditions in the context of the eigenmode behavior near the boundaries under weak and strong coupling, respectively.

III Density of states

The density of states (DOS), $\rho(\omega)$, gives the number of states per unit frequency interval at the frequency $\omega$. It is defined as

[TABLE]

where $N$ is the number of particles in the disordered chain and $\left\langle\right\rangle$ indicates averaging with respect to the independent realizations of the disorder. Thus, we have $\int{\rm d}\omega\rho(\omega)=1$. We compute $\rho(\omega)$ for the isolated disordered chain by constructing the histogram of the frequencies $\omega(=\sqrt{\lambda})$ obtained by diagonalizing the following $N\times N$ matrix:

[TABLE]

where $\lambda$ are the eigenvalues of $H$. For an ordered chain, we know Ashcroft and Mermin (1976)

[TABLE]

where $\omega_{m}=2\omega_{0}$ is the maximum allowed frequency for the phonons, and $\omega_{0}$ is the natural frequency of an ordered chain. Thus, for an ordered chain $\rho(\omega)$ diverges as $\omega\rightarrow\omega_{m}$, as discussed in the main text and consequently $\rho(\omega)=0$ for $\omega>2\omega_{0}$. The divergence of density of states at $\omega=\omega_{m}$ for the ordered chain $(W=0)$ corresponds to the van-Hove singularity Ashcroft and Mermin (1976).

IV Details of the transfer matrix method to compute transmission coefficient, $\tau(\omega)$

To compute the transmission coefficient $\tau(\omega)$ for a phonon of frequency of $\omega$ as it passes through the disordered chain of length $L$, we construct the transfer matrix for the disordered chain as described below. We can recast Eq. (S.2) in the following matrix form:

[TABLE]

where, $m_{n}=M$ for $(n=1,2,\cdots,N)$ and $m_{n}=m$ otherwise. We can relate the two ends of the disordered chain using the matrix $M_{n}$ iteratively, and thus we have

[TABLE]

Eq. (S.29) relates the amplitudes of the waves on the two sides (i.e. in the ordered region) of the disordered chain. In the ordered region, we can express the amplitude of the displacement at any lattice site as the superposition of plane waves. For example, in the region $n>N$, for phonons with wave number $q$, we can write:

[TABLE]

where, $A$ and $B$ are two constants. The left hand side of Eq. (S.29) can be written as

[TABLE]

Now, multiplying both sides of Eq. (S.29) by $Q^{-1}$ from the left, we get

[TABLE]

where,

[TABLE]

and we have introduced the transfer matrix $T^{(N)}$, connecting the waves on two sides of the disordered chain, as

[TABLE]

Thus, from Eq. (S.32), we have

[TABLE]

We can simplify Eq. (S.35) further by considering the following situation. Let us assume a phonon with wave number by $q$ coming from the right side $(n>N)$ of the disordered chain. There is only a transmitted wave on the left side $(n<1)$ of the disordered chain while on the right side of the system we can have a reflected wave along with the incident wave. Therefore the amplitudes $x_{0}$ and $x_{-1}$ of the wave with frequency $\omega$ can be taken as

[TABLE]

Here, wave number $q$ and frequency $\omega$ are related by the phonon dispersion relation Ashcroft and Mermin (1976)

[TABLE]

From Eq. (S.36), we obtain

[TABLE]

Thus, Eq. (S.35) reduces to

[TABLE]

where $T^{(N)}_{ij}$ represents the $(i,j)$ element of the transfer matrix $T^{(N)}$ defined in Eq. (S.34). Finally, the transmission coefficient, $\tau(\omega)$, for an incoming wave with frequency, $\omega$, is given by Markos and Soukoulis (2008)

[TABLE]

Since, in our case the disorder does not break the time reversal symmetry, the transmission coefficient for phonons coming from the right side of the disordered system is same as those coming from the left side of the system Markos and Soukoulis (2008).

Using Eq. (S.29), one can obtain $x_{N+2}$ and $x_{N+1}$ and then calculate the transmission coefficient, $\tau(\omega)$, using Eq. (S.41). For a given $L$, once we obtain the transmission coefficient for different $\omega$, we can compute the thermal conductance for the disordered chain using the Landauer formula given in main article. Finally, we average the results over $500$ independent realizations of the disorder.

V Dependence of transmission coefficient on $L$ for a fixed frequency

For a given strength of disorder, the transmission coefficient $\tau(\omega)$ depends on the length $L$ of the disordered chain. The $L$-dependence of $\ln\tau(\omega)$ for two arbitrarily chosen frequencies: (a) $\omega=0.5$ and (b) $\omega=0.8$ are shown in Fig. S.4(a-b). We set $k=1$ and $m=1$. In general, we find that $\langle\ln\tau(\omega)\rangle$ decays linearly with $L:\langle\ln\tau(\omega)\rangle\propto{-\frac{L}{\xi(\omega)}}$ for different different values of $W$. Thus, for a given $W$, the inverse of the slope of $\langle\ln\tau(\omega)\rangle-L$ curve gives us the localization length associated with the particular mode of frequency $\omega$.

VI Scaling of thermal conductance with $L$ for strong coupling, $k(\gg 1)$

In the main article, we discuss results for $k=0.01$ (representing weak coupling between the system and the heat bath) and $k=1.0$ (intermediate coupling). Here, we show the effect of large $k(\gg 1)$ on the transmission coefficient and the scaling of thermal conductance. The transmission coefficient for different frequencies with $k=100$ is shown in Fig. S.5(a), where we consider a disordered chain (power law disorder) with $N=100$ and disorder strength $\epsilon=2.50$. For $k=100$, $\tau(\omega)$ develops well separated Lorentzians having area and width given by Eq. (S.15)-S.16. We have also verified that the scaling exponents for thermal conductance with $k=100$ and $k=0.01$ remain same in different disorder regimes as demonstrated in Fig. S.5(b) (for uniformly disordered chain) and Fig. S.5(c) (for power-law disorder). When expressed in terms of the variable $\omega_{L}^{s+1}/\omega_{T}$, thermal conductance for one dimensional disordered harmonic chains with different models of disorder collapse approximately for small $\omega_{L}^{s+1}/\omega_{T}$ (see Fig. S.5(d)). For a given $T$, we expect this collapse to work for large $L$ (.i.e., small $\omega_{L}$ and fixed $\omega_{T}$).

VII Scaling of thermal conductance with $T$ for different coupling $k$

The dependence of thermal conductance $G_{qm}(T)$ on temperature $T$ for a given $L$ for different disorder and coupling strength $(k)$ is shown in Fig. S.6(a-f). We first analyze the case $k=1.0$. We find that $G_{qm}(T)$ is close to unity for small $T(\ll\hbar\omega_{L})$ and it decays as $T^{-1}$ for $T>\hbar\omega_{L}$ in all disorder regimes. Here, $\xi(\omega_{L})=L$. For a given $L$, we can understand the dependence of $G_{qm}(T)$ on $T$ by considering the following observations. Phonons with frequencies such that the associated localization length is less than the system size, $L$, are localized and consequently do not contribute to the thermal conductance. If $T$ is such that $\omega_{T}<\omega_{L}$ then the modes which are delocalized are thermally excited. For a given $L$, as $T$ is increased from very small value to $T^{*}_{W}=\hbar\omega_{L}$, more and more delocalized (ballistic) phonons contribute to the thermal conductance $G(L,T)$ and consequently $G(T)$ should initially increase with $T$ for a given $L$. In Fig. S.6(a), we find that for a given $W$, $G_{qm}(T)$ remains unity for $T\ll T^{*}_{W}$ implying that $G(T)$ increases linearly with $T$ for $T\ll T^{*}_{W}$ (note that $T$ is included in the definition of $G_{qm}(L,T)$). For $T>T^{*}_{W}$, the modes which are localized $(\xi(\omega)<L)$ are thermally excited and these modes do not contribute to the heat transport, causing $G(T)$ to saturate for $T>T^{*}_{W}$. Such a saturation is reflected as the linear decay of $G_{qm}(T)(\propto T^{-1})$ in Fig. S.6(a), which is independent of the disorder strength, $W$. Note that for a given $L,T^{*}_{W}$ decreases with increasing $W$ and subsequently the temperature regime for $G_{qm}\approx 1$ shrinks to smaller values with increasing $W$.

We also find that the scaling $G_{qm}(T)\propto T^{-1}$ for $T>\hbar\omega_{L}$ holds irrespective of nature of disorder and coupling strength (see Fig. S.6(a-f)), following the theoretical prediction made in the main article.