Ballistic heat conduction and mass disorder in one dimension

It is commonly believed that in the absence of mass disorder, phonons propagate ballistically without dissipation in pristine one-dimensional (1D) harmonic systems, resulting in a linear relationship between the thermal conductivity κ and length L [1]. On the other hand, it is well-established from numerical and analytical results that the presence of mass disorder causes phonon localization [1–10], leading to a sublinear κ-L relationship $(\kappa\propto L^{\alpha}$ where α < 1) in the asymptotic (L → ∞) limit. This heat conduction scaling behaviour is commonly termed subballistic or superdiffusive in the literature [1, 5, 6]. It also implies that the presence of mass disorder and localization [3, 11] is incompatible with ballistic heat conduction although its mere presence is not sufficient for the normal heat diffusion behaviour [1, 5, 6].

However, recent measurements of the thermal conductivity in micrometer-long SiGe-alloy nanowires at room temperature by Hsiao and co-workers have yielded a linear κ-L relationship [12], suggesting the possibility of micrometer-scale ballistic phonon transport despite the presence of mass disorder. Thus, a fresh examination of the relationship between phonon localization and propagation in 1D systems is needed to understand this result i.e. the coexistence of ballistic heat conduction and mass disorder. This requires us to determine the conditions for ballistic conduction to occur in a disordered 1D structure. It also raises the question of how localization varies with the mass composition of the system and if there are any other experimentally detectable physical signatures, apart from the linear κ-L relationship, that can distinguish ballistic and subballistic heat conduction. This is important since the experimental determination of the κ-L relationship requires the fabrication and measurement of different samples of varying length [13], and the process of fabricating samples of different sizes may introduce variability in the measured thermal conductivity. Such measurement variability between samples can obscure the signs of ballistic heat conduction. Therefore, it would be useful to be able to verify this phenomenon in a single sample.

To understand the coexistence of ballistic heat conduction and mass disorder, we revisit in this paper the phenomenon of phonon propagation in the disordered 1D harmonic lattice, which we take to represent an idealization of real 1D systems with mass disorder [9, 14, 15]. Although this model may be simple, it is sufficient for our purpose of demonstrating and clarifying the coexistence of ballistic heat conduction and mass disorder. Thus, it can be applied to elucidate the essential physics of phonon transport in quasi-1D systems with disorder as well as for the analysis and interpretation of experiments [12, 16–18]. Its simplicity also enables us to relate our numerical results to known analytical results from conventional wave transmission theory and to pinpoint the experimentally relevant physical parameters affecting phonon transmission in the presence of mass disorder. The purpose of our paper is to analyze the coexistence of ballistic heat conduction and mass disorder within the established framework of phonon localization. We show how ballistic heat conduction can occur in the presence of disorder and how it can be distinguished from subballistic heat conduction. Although the asymptotic scaling theory of heat conduction [2, 4, 14] and the localization phenomenon are rather well understood, the implications for finite systems have not been clarified. Neither has there been any attempt to use these disordered 1D models to relate the heat conduction phenomenon to experimentally relevant variables such as impurity concentration and temperature.

The organization of our paper is as follows. We first analyze how mass disorder affects phonon transmittance. While it is known that mass disorder attenuates phonon propagation and that the transmittance attenuation length (λ) scales as the inverse square frequency (ω⁻²) for low-frequency modes in the L → ∞ and ω → 0 limit [2, 4], there is no straightforward numerical presentation of this result in the context of heat conduction. We use the standard nonequilibrium Green's function (NEGF) method to compute the frequency-dependent transmittance for chains of different lengths with different values of impurity mass. We compare the transmittance with known analytical estimates and find a reasonably good fit over a wide range of frequencies, allowing us to find the dependence on mass composition. A close analogy to the Beer-Lambert law is also found. This enables us to find a good closed form estimate of the thermal conductance, from which we are able to deduce its dependence on mass composition. We also derive from the closed form expression the dependence of the thermal conductance on length, mass composition and temperature. More importantly, we show and summarize how these dependences differ in the ballistic and subballistic regime. Given that large lattices are expected to exhibit a sublinear κ-L relationship in the L → ∞ limit, we identify the crossover length scale (L_C) below which phonon transport acquires a fully ballistic character independent of mass composition and above which phonon transport is affected by lattice disorder. We also show how the crossover length can be tuned by changing the impurity concentration and mass difference.

We choose as our model system the familiar 1D harmonic lattice first studied by Dyson [19]. Our system differs from that used in other papers in that the atomic mass is not treated as a continuous random variable. Rather, we have a random binary lattice (RBL), where the positions of the component atoms are randomly shuffled, which resembles more realistic systems with mass disorder and allows us to quantify the effect of experimentally relevant variables such as impurity concentration and mass difference on thermal conduction. For simplicity, only adjacent atoms are coupled. We do not expect our results to be qualitatively different with longer-range interactions since the interaction typically grows weaker with interatomic distance. The atoms are only permitted to move longitudinally and no attempt is made to include any anharmonic interaction in our model.

Our system consists of three parts. In the middle, there is a finite-size 'conductor' of N equally spaced atoms. On either side, there is a homogeneous lead i.e. a semi-infinite chain with no mass disorder. In effect, we have a finite disordered system embedded in an infinite homogeneous 1D lattice. Coupling between adjacent atoms is governed by the harmonic spring term $V(x_{i},x_{i+1})=\frac{1}{2}k(x_{i}-x_{i+1})^{2}$ where k is the spring constant and x_i is the displacement of the i-th atom from its equilibrium position. In the RBL, there are two species of atoms which we label 'A' and 'B'. Species 'A' is taken to be the substitutional impurity and exists only within the conductor. The rest of the atoms in the conductor and the leads are of species 'B'. We set the mass of the B atoms to m_B = 4.6 × 10⁻²⁶ kg (the mass of a Si atom), the spring constant to k = 32 Nm⁻¹ (the approximate strength of the Si–Si bond) and the interatomic spacing to a = 0.55 nm. Our choice of parameters follows that in [20]. We vary the ratio of the masses $(R=\frac{m_{\rm A}}{m_{\rm B}})$ and concentration (c_A = fraction of A atoms, c_B = 1 − c_A) in our simulations. The 1D lattice size is varied from N = 100 to 20 000, and the equivalent range of chain length values (L = Na) is L = 55 nm to 1.1 µm. A schematic of the system is shown in figure 1.

Zoom In Zoom Out Reset image size

To simulate phonon transport, we use the nonequilibrium Green's function (NEGF) method which is described in detail in [20, 21]. The method calculates the frequency-dependent energy transmittance function Ξ(ω) by modelling the conductor as an open quantum system and the semi-infinite leads as self-energies (open boundaries). The Green's function for the conductor is

$$\begin{equation*} {\bit G}(\omega)=\left[\omega^{2}+{\rm i}0^{+}-{\bit K_{\rm C}}-{\Sigma_{L}}(\omega)-{\Sigma_{R}}(\omega)\right]^{-1}, \end{equation*}$$

where Σ_L (Σ_R) is the so-called self-energy corresponding to the left (right) semi-infinite lead; K_C is the force constant matrix of the conductor, and its (i, j) matrix element which couples the i-th and j-th atom is given by

$$\begin{equation*} [{\bit K_{\rm C}}]_{ij}=(m_{i}m_{j})^{-\frac{1}{2}}\frac{\partial^{2}V}{\partial x_{i}\partial x_{j}}\ . \end{equation*}$$

To compute Ξ(ω), we use the Caroli transmission formula [22] Ξ(ω) = Tr(Γ_LGΓ_RG^†), where Γ_L (Γ_R) is the term coupling the conductor to the left (right) lead [20, 21]. As with any study of disordered systems, it is necessary to express the computed variables in terms of their ensemble average. In this paper, the ensemble-average transmittance function 〈Ξ(ω)〉 is obtained by averaging Ξ(ω) over 100 realizations of the RBL.

It is more useful to work with the thermal conductance σ which is related to the thermal conductivity via the relation κ ∝ σL. The thermal conductance is computed using the Landauer formula [20, 21]:

$$\begin{equation} \sigma(T)=\frac{1}{2\pi}\int_{0}^{\infty}\hbar\omega\frac{{\rm d}f(\omega)}{{\rm d}T}\langle\Xi(\omega)\rangle\ {\rm d}\omega \end{equation} \tag{ 1 }$$

where T and ℏ are respectively the temperature and reduced Planck constant, and $f(\omega)=[\exp(\frac{\hbar\omega}{k_{\rm B}T})-1]^{-1}$ is the Bose–Einstein occupation factor. The expression in equation (1) integrates the frequency-dependent transmittance spectrum and the differential heat capacity over the entire frequency range (0–∞). In the high temperature limit, we have $f(\omega)\approx\frac{k_{\rm B}T}{\hbar\omega}$ so that equation (1) becomes:

$$\begin{equation} \lim_{T\rightarrow\infty}\sigma(T)=\frac{k_{\rm B}}{2\pi}\int_{0}^{\infty}\langle\Xi(\omega)\rangle\ {\rm d}\omega, \end{equation} \tag{ 2 }$$

which is proportional to the total transmittance. In the case of the homogeneous chain (m_A = m_B), the transmittance is Ξ(ω) = Θ(ω_C − ω) where $\omega_{\rm C}=2\sqrt{k/m_{\rm B}}$ is the frequency cutoff, and equation (2) yields a length-independent thermal conductance i.e. $\lim_{T\rightarrow\infty}\sigma(T)=\frac{k_{\rm B}\omega{}_{\rm C}}{2\pi}$. This example implies that a length-dependent transmittance function is needed for the conductance to vary with L.

3.1. Length dependence of transmittance spectrum

Given that the transmittance spectrum is expected to vary with L, the question arises as to what its dependence on L should be. In general, the transmittance function can be interpreted as the ease with which a phonon traverses the conductor. To identify the functional relationship, we first compute 〈Ξ(ω)〉 for different values of chain length L at different values of R and c_A. Figure 2 shows 〈Ξ(ω)〉 decreasing with increasing ω and L for a particular combination of c_A and R(c_A = 0.5 and R = 2), with the rate of decrease faster for higher-frequency modes. From our calculations, we find that 〈Ξ(ω, L)〉 obeys the simple numerical relation:

$$\begin{equation} \langle\Xi(\omega,L)\rangle=\langle\Xi(\omega,L_{0})\rangle\exp\left[-\frac{L-L_{0}}{\lambda(\omega)}\right]\ , \end{equation} \tag{ 3 }$$

for ω < ω_C and L > L₀, analogous to the Beer-Lambert law in optics. The parameter λ(ω) is the frequency-dependent attenuation length and can be interpreted as the distance through which the phonon has its transmitted square amplitude reduced by a factor of 1/e. Numerical fitting yields the numerical relationship λ ∝ ω⁻².

Zoom In Zoom Out Reset image size

To determine the exact relationship between λ and ω, we use the recent result from Das and Dhar [23], who proved an intuitive and useful formula relating the energy transmittance Ξ(ω) and the plane-wave transmission function τ(ω) for the 1D harmonic chain:

$$\begin{equation} \Xi(\omega)=|\tau(\omega)|^{2}\ . \end{equation} \tag{ 4 }$$

In the disordered harmonic chain with N atoms, the ensemble average of its asymptotic amplitude has the form of an exponentially decaying function i.e. lim_N → ∞|τ(ω)| = exp[−γ(ω)N] where γ is the dimensionless Lyapunov exponent [2, 4]. This implies that the attenuation length in equation (3) is

$$\begin{equation} \lim_{L\rightarrow\infty}\lambda(\omega)=\frac{a}{2\gamma(\omega)}\ . \end{equation} \tag{ 5 }$$

In the low frequency (ω → 0⁺) limit, the closed form expression for γ is known [2] and given by

$$\begin{equation} \lim_{\omega\rightarrow0}\gamma(\omega)=\frac{\langle\delta m\rangle^{2}\omega^{2}}{8k\langle m\rangle}\ , \end{equation} \tag{ 6 }$$

where 〈δm〉 and 〈m〉 are respectively the standard deviation and the mean of the conductor atomic mass. Combining equations (6) and (5), we obtain the expression for the low-frequency attenuation length

$$\begin{equation} \lambda_{0}(\omega)=\frac{4ka\langle m\rangle}{\langle\delta m\rangle^{2}\omega^{2}} \end{equation} \tag{ 7 }$$

which we plot in figure 3 together with λ(ω) for R = 1/8, 1/2, 2 and 8 at c_A = 0.1 and 0.5. We find close agreement between the two at low frequencies (ω < 1 THz). At higher frequencies, λ₀(ω) is slightly smaller but still an excellent approximation to the numerically computed λ(ω). The approximation λ₀ however diverges increasingly from the extracted attenuation length λ as λ₀ approaches a, the interatomic spacing. This is more evident for the R = 8 curve in figure 3 when the impurity mass m_A is much greater than m_B.

Zoom In Zoom Out Reset image size

If we interpret λ₀(ω) as an impurity-limited mean free path (MFP) and assume a very dilute impurity concentration n_imp = c_A ≪ 1, then we have 〈δm〉² = c_A(1 − c_A)(m_A − m_B)² ≈ n_imp(m_A − m_B)² and obtain $\lambda_{0}\propto n_{{\rm imp}}^{-1}(m_{\rm A}-m_{\rm B})^{-2}\omega^{-2}$ i.e. the attenuation length varies as $n_{{\rm imp}}^{-1}$. This parallels the Beer-Lambert law for photons where the attenuation length is inversely proportional to the absorber concentration. In three dimensions, the impurity-limited MFP $(\lambda_{{\rm imp,3D}}\propto n_{{\rm imp}}^{-1}(m_{\rm A}-m_{\rm B})^{-2}\omega^{-4})$ has identical dependence on impurity concentration and mass difference but a very different frequency dependence [24].

3.2. Length, mass composition and temperature dependence for ballistic and subballistic thermal conductance

Having numerically verified the form of the transmittance function, we proceed to calculate the thermal conductance to determine its dependence on length, mass composition and temperature. From equations (1), (6) and (5), we arrive at an analytical expression for estimating the length-dependent thermal conductance:

$$\begin{equation} \sigma_{0}(T,L)=\frac{1}{2\pi}\int_{0}^{\omega_{\rm C}}\hbar\omega\frac{{\rm d}f}{{\rm d}T}\exp\left[-\frac{L\langle\delta m\rangle^{2}\omega^{2}}{4ka\langle m\rangle}\right]\ {\rm d}\omega. \end{equation} \tag{ 8 }$$

This approximation relies on the $\lim_{L\rightarrow\infty,\omega\rightarrow0^{+}}\lambda(\omega)$ expression in equation (5) i.e. it makes use of the low-frequency approximation for λ(ω) in the L → ∞ limit. At finite L, higher-frequency modes contribute more to the phonon transmittance and their attenuation behaviour may not follow that of equation (6).

To see how well equation (8) approximates the thermal conductance σ (calculated from NEGF) especially for finite L, we compute σ₀ for R = 2 and c_A = 0.1 at T = 300 K and plot it in figure 4 together with the corresponding σ. σ₀ scales as L^−0.5 and agrees well with σ but is consistently ∼20 percent smaller. We have also computed σ₀ for other values of R and c_A and found that it is also ∼20 percent smaller than σ computed from NEGF. This small discrepancy is expected given that equation (8) uses the low-frequency approximation λ₀ which overestimates the transmittance attenuation.

Zoom In Zoom Out Reset image size

Since figure 4 shows that the expression for σ₀ from equation (8) can fit the numerically computed σ, it means that we can use equation (8) to infer the asymptotic behaviour of the conductance in the L → ∞ and L → 0 limits. In the L → ∞ limit, only very low frequency modes contribute to the integral in equation (8), allowing us to use the high-temperature approximation i.e.

$$\begin{equation*} \frac{{\rm d}f}{{\rm d}T}=\frac{k_{\rm B}}{\hbar\omega}\ , \end{equation*}$$

so that equation (8) becomes

$$\begin{equation*} \lim_{L\rightarrow\infty}\sigma_{0}(T,L)\approx\frac{k_{\rm B}}{2\pi}\int_{0}^{\infty}\exp\left[-\frac{L\langle\delta m\rangle^{2}\omega^{2}}{4ka\langle m\rangle}\right]\ {\rm d}\omega\ . \end{equation*}$$

This yields

$$\begin{equation} \lim_{L\rightarrow\infty}\sigma_{0}=\left[\frac{ka\langle m\rangle k_{\rm B}^{2}}{4\pi\langle\delta m\rangle^{2}}\right]^{\frac{1}{2}}\frac{1}{\sqrt{L}} \end{equation} \tag{ 9 }$$

which gives the well-known $\lim_{L\rightarrow\infty}\sigma\propto L^{-0.5}$ relationship [2, 4, 9], indicating subballistic (superdiffusive) thermal transport. In the L → 0 limit, the conductor is transparent to all phonon modes with frequency below ω_C (the cutoff frequency determined by the atomic mass in the leads, m_B), i.e. Ξ(ω) ≈ Θ(ω_C − ω), and we recover the expression for the finite length-independent thermal conductance,

$$\begin{equation} \lim_{L\rightarrow0}\sigma_{0}=\frac{1}{2\pi}\int_{0}^{\omega_{\rm C}}\hbar\omega\frac{{\rm d}f}{{\rm d}T}\ {\rm d}\omega \end{equation} \tag{ 10 }$$

which is independent of the mass composition. We note that the ballistic conductance in equation (10) is proportional to the heat capacity of the system.

Although the model described by equation (8) is very simple, the variation in its thermal transport character in the different limits, as suggested by equations (9) and (10), has significant qualitative implications for more realistic experimental systems, apart from its $\sigma\propto L^{-0.5}$ property. We use them to deduce the difference in the dependence on the various physical variables. Firstly, in the ballistic limit, the thermal conductance is independent of its mass composition. On the other hand, in the subballistic limit, the conductance is inversely proportional to the standard deviation of the mass concentration, and in the dilute limit (c_A → 0), it should scale as $\sigma\propto n_{{\rm imp}}^{-0.5}|\Delta m|^{-1}$, where Δm = m_A − m_B is the mass difference, since 〈δm〉² = c_Ac_B(m_A − m_B)² ≈ n_impΔm² in equation (9). Secondly, the temperature dependence of the conductance is different for the subballistic and ballistic case. In the subballistic case (L → ∞), the range of conducting phonon modes decreases as L^−0.5 and only the low-frequency modes contribute to thermal transport. Hence, the temperature dependence of the thermal conductance is negligible. However, in the ballistic case (L → 0), the range of conducting modes (ω < ω_C) is determined by the phonon occupation in the leads. Thus, the conductance is more sensitive to temperature change and in the low-temperature limit (T ≪ ℏω_C/k_B), we obtain σ₀ ∝ T. We summarize the differences in the scaling dependence of the thermal conductance in table 1. These differences can be used to distinguish the ballistic and subballistic transport regimes.

Table 1. Scaling dependence of the thermal conductance σ on various physical variables (length, temperature, mass difference and impurity concentration) in the ballistic (L ≪ L_C) and subballistic (L ≫ L_C) regime in the 1D random binary lattice.

Physical variable	Scaling exponent α
	Ballistic (L ≪ LC)	Subballistic (L ≫ LC)
Length $(\sigma\propto L^{\alpha})$	0	−0.5 to −0.4
Temperature $(\sigma\propto T^{\alpha})$	1	0
Mass difference (σ ∝ |Δm|α)	0	−1
Impurity concentration $(\sigma\propto n_{{\rm imp}}^{\alpha})$	0	−0.5

3.3. Ballistic to subballistic crossover length

Given that thermal conduction is ballistic in the L → 0 limit and subballistic in the L → ∞ limit, there must exist an intermediate length scale that marks the crossover from ballistic to subballistic thermal transport although we do not observe one in figure 4. Since the thermal conductance depends on the mass composition in the L → ∞ limit (equation (9)), it implies that the crossover length varies with mass composition. This can be increased by weakening the amount of mass disorder so that ballistic heat conduction is detectable in the range of L values simulated.

For a dilute impurity concentration and small mass difference, the mass variance (〈δm〉²) is small and the mass disorder is weak. Thus, the crossover length between ballistic and subballistic thermal transport can be large. The crossover length scale L_C can be estimated by estimating the minimum possible attenuation length which is determined by the cutoff frequency ω_C. Using equation (7), we set L_C = λ₀(ω_C) which yields L_C = a〈m〉 m_B/[c_A(1 − c_A)Δm²]. This suggests that the ballistic regime can be extended by reducing the impurity concentration as well as mass ratio (R). Instead of using R = 2, we pick a small R so that the mass difference would be small (Δm ≪ m_B). Figure 5 shows the NEGF-computed σ and σ₀ from equation (8) as a function of L (56 nm to 88 µm) for (R, c_A) equal to (1.05, 0.1), (1.05, 0.5), (1.1, 0.1), (1.1, 0.5), (1.2, 0.1) and (1.2, 0.5) in the high temperature (T → ∞) limit. We find good agreement between σ and σ₀ over the range of L considered, with the agreement better for larger c_A and Δm. The corresponding crossover length ranges between L_C ≈ 0.05 to 2 µm. Below L < L_C, σ and σ₀ converge to the ballistic limit $\sigma=\frac{k_{\rm B}\omega_{\rm C}}{2\pi}$ while for L ≫ L_C, $\sigma\propto L^{-\alpha}$ with α < 0.5 and σ₀ ∝ L^−0.5, respectively. The discrepancy in the length dependence between σ and σ₀ may be a consequence of the small mass variance. In the 〈δm〉 → 0 limit, we expect fully ballistic thermal transport i.e. $\sigma\propto L^{0}$. Thus, the scaling exponent would deviate from −0.5 toward 0 for small 〈δm〉.

Zoom In Zoom Out Reset image size

In summary, we have studied heat conduction in the 1D random binary lattice model with mass disorder. We have shown that the energy transmittance attenuates exponentially with lattice size, similar to the Beer-Lambert law for photons. The attenuation length is shown numerically to have a frequency dependence of λ ∝ ω⁻² for low frequency modes (as predicted by localization theory) and also for high frequency modes (not predicted by localization theory). The existence of a maximum cutoff frequency implies a minimum attenuation length below which the transmittance of the phonons in the system converges to unity. This indicates that ballistic heat conduction can coexist with mass disorder for finite systems. The crossover length between ballistic and subballistic heat conduction can be estimated based on the cutoff frequency ω_C and the low-frequency expression for the attenuation length λ₀. This analytical estimate has been shown to agree closely with the numerical crossover length computed with NEGF. We have also derived the explicit length, temperature and mass composition dependence of the thermal conductance and shown how they differ in the ballistic and subballistic limit. The differences can be used to distinguish ballistic and subballistic heat conduction in experimental systems. By making use of the mass composition dependence of the attenuation length, we have shown how the crossover between ballistic and subballistic heat conduction can be manipulated by changing the impurity mass difference and concentration.

The authors gratefully acknowledge the financial support from the Agency for Science, Technology and Research (A*STAR), Singapore and the use of computing resources at the A*STAR Computational Resource Centre, Singapore.