Comparison of linear and quadratic dispersion models for phonon transport in one-dimensional mass-disordered systems

Acoustic phonons in crystals have a linear dispersion relation, ω_q = sq, and travel at a finite speed, s, even in the low frequency limit. This leads to a variety of anomalies in phonon transport phenomena. In the low frequency limit, the mean free path of the phonons becomes very large^1,2 because of the vanishing scattering rate of a typical scattering mechanism, such as phonon-defect and phonon–phonon scattering.^3–5 This forces us to introduce the boundary scattering whose characteristic length is, however, hard to be predicted from first-principles calculations.⁵ For one-dimensional mass-disorder systems, the phonon thermal conductance K shows anomalous diffusion and does not obey the Fourier law.^6–16 For a free boundary condition in a system consisting of a channel (with random impurities) of N atoms with attached reservoirs (without impurities) of the same properties as the channel, the thermal conductance exhibits the super-diffusive nature⁶ of K ∝ N^−1/2. On the other hand, for a fixed boundary condition where the reservoirs with properties different from those of the channel, the thermal conductance is predicted to be the sub-diffusive nature¹² of K ∝ N^−3/2.

For a quadratic dispersion ω_q = cq²/2, the phonon speed, cq, converges to zero in the long wavelength limit, which does not contribute to the thermal transport. It is, therefore, interesting to investigate the consequences of the quadratic dispersion on the anomalous diffusion. The phonon mode with the quadratic dispersion is present in two-dimensional materials, such as graphene^17,18 and transition metal dichalcogenides.¹⁹ It is also important in this sense to investigate the effects of the quadratic dispersion on the thermal transport. We calculated²⁰ the thermal conductance of the in-plane mode (linear dispersion) and that of the out-of-plane mode (quadratic dispersion) in an armchair-edge graphene nanoribbon of length L and found that the thermal conductance due to the in-plane mode exhibits the super-diffusive nature of K ∝ L^−1/2. On the other hand, the thermal conductance due to the out-of-plane mode did not show a clear power-law behavior in the calculated L range of up to L = 430 µm. To obtain a definite answer to the consequences of the quadratic dispersion on the anomalous diffusion, we have to calculate a longer system so that we can identify a clear power-law behavior.

In the present study, we use a one-dimensional atomic chain model that shows the quadratic dispersion and investigate the thermal conductance for a longer system (N ∼ 10⁹) to obtain a clear power-law behavior. We assume a reservoir–channel–reservoir-type of the system and adopt the non-equilibrium Green function (NEGF) method^21–29 to calculate the thermal conductance. The efficiency of the numerical algorithm is of utmost importance for simulating larger systems, and we introduce the R-matrix propagation algorithm^30–32 to reduce the computational burden. We focus on the impact of the dimensionality of the dispersion relation on the thermal conductance and consider only a homogeneous system (except for the presence of impurity atoms in the channel) with the free boundary condition. We calculate a single phonon transmission function through the channel, and the many-body interaction between phonons is neglected.

This paper is organized as follows. In Sec. II, we first introduce the models used for the present study: a one-dimensional linear dispersion model and a one-dimensional quadratic dispersion model. After that, the calculation method for the thermal conductance and the R-matrix propagation algorithm are briefly summarized for the sake of the completeness. Section III shows the results and discussion. We first present the results on the linear dispersion model and confirm that the present method properly reproduces the results of the previous research. The results on the quadratic dispersion model are then presented and discussed. A conclusion is given in Sec. IV.

We consider a one-dimensional linear dispersion model whose schematic diagram is given in Fig. 1(a). It consists of three regions: the finite channel region (C) of N atoms and two semi-infinite reservoirs (L and R) on the left and right sides. The atom numeration is chosen as n ≤ 0 for the left reservoir, 1 ≤ n ≤ N for the channel, and n ≥ N + 1 for the right reservoir. The mass of the nth atom is designated as M_n. The spring constant between the nearest neighboring atoms is constant throughout the system and is written as k. The lattice constant is a. We consider the binary-mass-disorder with the mass of the host atom [open circles in Fig. 1(a)] of m and the mass of the impurity atom [hatched circles in Fig. 1(a)] of M,

The impurity atoms exist only in the channel region with the impurity concentration of C = N_imp/N (N_imp is the number of impurity atoms). The reservoirs consist of the host atom; M_n = m for n ≤ 0 or n ≥ n + 1.

The equation of motion in the frequency domain can be written as

where $UL=(…,u−2,u−1,u0)T$⁠, $UC=(u1,u2,…,uN)T$⁠, and $UR=(uN+1,uN+2,uN+3,…)T$ are the column vectors of the mass-normalized displacement u_n of the nth atom for the left reservoir, channel, and right reservoir, respectively. The mass-normalized dynamical matrix Φ⁽¹⁾ is given by

with d_n = 2 k/M_n, $sn=−k/(MnMn+1)1/2$ for n = 1, 2, …, N and d₀ = 2 k/m, s₀ = −k/m. Here, we assume that there is no impurity on the channel edges at n = 1 and N. This condition can always be satisfied by an appropriate choice of the interface positions between the channel and the reservoirs. Note that in the one-dimension model, s_n is just a real number, but we keep the symbol of transpose T so that our equations are valid in a more general case.

The quadratic dispersion model can be constructed by introducing a second nearest neighbor interaction with the negative spring constant of −k/4 [see Fig. 1(b)]. In the following, we assume that N is an even number. The mass-normalized dynamical matrix Φ⁽²⁾ is then given by

where

for j = 1, 2, …, N/2, and

Here, we assume that there is no impurity on the channel edges at n = 1, 2, N − 1, N.

The dispersion relations of the linear and quadratic dispersion models can be obtained by considering a uniform system without impurities and are given by

for the linear dispersion model and

for the quadratic dispersion model. Here, we define ω₀ = (k/m)^1/2. In the low frequency limit (or long wavelength limit), the dispersions become

The dispersion relations are plotted in Fig. 2. Both models have the same spectral boundary: $ωq(1)=ωq(2)=2ω0$ at q = π/a.

The thermal conductance is calculated within the non-equilibrium Green function (NEGF) method. The channel retarded Green function of the N × N matrix is given by

where the superscript d specifies the dispersion model (d = 1 for the linear dispersion model and d = 2 for the quadratic dispersion model), $ΦC(d)$ is the mass-normalized dynamical matrix of the channel, and $Πα(d)(ω)$ (α = L, R) are the contact self-energies (whose dimension is ω²). The self-energies $Πα(d)(ω)$ have non-zero elements only in the upper-left d × d matrix for α = L or lower-right d × d matrix for α = R. The non-zero part of the contact self-energies, $πα(d)(ω)$⁠, can be written as

for the linear dispersion model, where q is a solution of $ω=ωq(1)$⁠, and

for the quadratic dispersion model, where q and γ are solutions of $ω=ωq(2)$ and ω = 2ω₀ sinh²(γa/2), respectively.

The phonon transmission function, T^(d)(ω), through the channel is given in terms of D^(d)(ω) and $Πα(d)(ω)$ as follows:

where $Γα(d)(ω)$ is defined as

When the left and right reservoirs are maintained at constant temperatures T_L and T_R, respectively, the thermal current, Q^(d), from the left reservoir to the right reservoir is given by

Here, $N(T;ω)=[exp(ℏω/kBT)−1]−1$ is the Bose–Einstein distribution function and k_B is the Boltzmann constant. By setting T_L = T + ΔT and T_R = T and taking the limit of ΔT → 0, we obtain the thermal conductance, K^(d), at temperature T as follows:

where $w(x)=(3/π2)(x2ex)/(ex−1)2$⁠. In the present study, we consider the high temperature limit [w(x) → 3/π²], and we have

When random impurities are considered, the thermal conductance is calculated by taking the sample average of the transmission functions with different impurity configurations.

We use the R-matrix method to evaluate the transmission function T(ω) from the retarded Green function D(ω). Note that in this subsection, we omit the superscript d to keep the notation simple. If we use a direct matrix inversion of the whole channel matrix in Eq. (12) whose size is N × N, then the computational time scales as ∼N³, which prevents us from simulating a larger system. The computational burden can be substantially reduced in the R-matrix method, where we first define an N × N matrix,

which represents Green’s function for the isolated channel. Here, ω is a real number without an infinitesimal shift +i0. From R(ω), we can compute D(ω) as follows:

with Π(ω) = Π_L(ω) + Π_R(ω). Since the coupling with the reservoirs is localized at the channel ends, the surface part of R(ω) is sufficient to perform phonon transport simulation. For the dispersion model d (=1, 2), the surface part consists of four d × d matrices $Rαβ(ω)$ (α, β = L, R), defined as

We then define the R-matrix $R(ω)≡R(ω)surface$ by

Since the dynamical matrices Φ^(d) (d = 1, 2) have the block-diagonal form with the block size d × d [see Eqs. (3) and (4)], the following propagation algorithm can be used to compute the channel R-matrix $R(ω)$⁠.

1. Compute a d × d matrix of $r1(ω)=[ω2−D1]−1$⁠, and define the R-matrix for the isolated first block with the elements $RLL=RLR=RRL=RRR=r1(ω)$⁠.

2. Repeat the iterations

   $RRR=ω2−Dk+1−SkTRLLSk−1,$

   (28)

   $RRL=RRRSkTRRL,$

   (29)

   $RLL=RLL+RLRSnRRL,$

   (30)

   $RLR=RRLT$

   (31)

   for k = 1, 2, …, N/d − 1 to obtain the R-matrix for the channel of N/d blocks. Note that D_k and S_k should be regarded as d_k and s_k for the linear dispersion model [see Eq. (3)].

We first present the results of the linear dispersion model. The channel length N dependence of the thermal conductance K⁽¹⁾ is shown in Fig. 3 for the mass ratio M/m = 1.1 and the impurity concentration C = 0.1. The thermal conductance is normalized by K₀ = k_Bω₀/π, which is the thermal conductance of the system without impurity. We calculate two cases in the spatial distribution of impurities: periodically arranged impurities (open circles) and randomly arranged impurities (closed circles). In the case of the periodic systems, the impurity atoms were placed in every C⁻¹ atoms in the channel region. In the case of random systems, we simulated ∼100 samples with different impurity configurations. For each sample, CN impurities were randomly placed in the channel region using computer-generated pseudorandom numbers.³³ The thermal conductance K depends on the exact locations of impurities, and the sample-averaged values are plotted in Fig. 3. The normalized standard deviation δK/K is found to be less than 1% for the linear dispersion model and 2% for the quadratic dispersion model. For the periodically arranged impurities, the thermal conductance K⁽¹⁾ is nearly independent of N. This implies that the wave nature of phonons is properly incorporated into the present simulation method. For the randomly arranged impurities, the thermal conductance decreases with the system size as K⁽¹⁾ ∝ N^−1/2 for N ≫ 1, which can be explained in terms of the existence of low frequency non-localized modes.^6,8

Figure 4 shows the color plot of the sample-averaged transmission function ⟨T⁽¹⁾(ω)⟩ on the (ω, N)-plane for the mass ratio M/m = 0.1 and the impurity concentration C = 0.1. We see that higher frequency modes are cut off as the channel length becomes longer. The demarcation frequency, ω_d, between high frequency localized modes and low frequency non-localized modes is given by⁸ 

where ⟨m⟩ = (1 − C)m + CM is the average mass and ⟨δm²⟩ = C(1 − C)(M − m)² is the variance. The demarcation frequencies are plotted as the white dashed line in Fig. 4. Since $ωd(1)$ is inversely proportional to N^1/2, the thermal conductance in the high temperature limit exhibits the super-diffusive nature of K⁽¹⁾ ∝ N^−1/2.

Matsuda and Ishii⁶ derived an analytical expression of K⁽¹⁾ for N ≫ 1 based on the Kubo formula,

Figure 5 shows a comparison between the simulated K⁽¹⁾ (marks) and the analytical expression $KMI(1)$ (lines) for the mass ratio 0.1 ≤ M/m ≤ 10 and the impurity concentration C = 0.01, 0.1, and 0.5. The analytical formula $KMI(1)$ is systematically about 20% smaller than the simulated results, as was observed in Ref. 15 (see also Fig. 7).

We next present the results on the quadratic dispersion model. Figure 6 shows the color plot of the sample-averaged transmission function ⟨T⁽²⁾(ω)⟩ on the (ω, N)-plane for the mass ratio M/m = 0.1 and the impurity concentration C = 0.1. As in the linear dispersion model, the quadratic model shows that higher frequency modes are cut off as the channel length increases. It can be seen that the demarcation frequency steeply depends on the channel length, N, as ∝N⁻¹. By calculating the sample-averaged transmission functions ⟨T⁽²⁾(ω)⟩ for a parameter range 0.1 ≤ M/m ≤ 10 and 0.01 ≤ C ≤ 0.5, we find that the demarcation frequency is well approximated by

which is shown by the white dashed line in Fig. 6.

Figure 7 shows the channel length N dependence of the thermal conductance K^(d) for the mass ratio M/m = 1.1 and the impurity concentration C = 0.1. Only the random impurity case is plotted. Open squares show the thermal conductance K⁽²⁾ in the quadratic dispersion model, and closed circles show K⁽¹⁾ in the linear dispersion model. We find that K⁽²⁾ is inversely proportional to N for N ≫ 1 showing the normal-diffusion. Since the analytical formula $KMI(1)$ of Eq. (33) can be written as $KMI(1)=(kB/8π)ωd(1)$⁠, we expect that an analytical formula for K⁽²⁾ could be written as

As shown by the dotted line in Fig. 7, this analytical formula is found to reasonably well reproduce the simulated results. (It is about 10% larger than the simulated results.)

We calculated the thermal conductance of the one-dimensional mass-disordered atomic chains with the linear and quadratic dispersion relation under the free boundary condition. We used the R-matrix method to reduce the computational burden, which enables us to compute the thermal conductance for a longer system of the size of up to N ∼ 10⁹. For the linear dispersion model, the thermal conductance exhibits the super-diffusive nature of K⁽¹⁾ ∝ N^−1/2. The NEGF results are consistent with the analytical results for the demarcation frequency and the thermal conductance. For the quadratic dispersion, we find that the thermal conductance shows the normal-diffusion of K⁽²⁾ ∝ N⁻¹. An analytical formula of the demarcation frequency $ωd(2)$ is proposed, which agrees with the simulated transmission functions ⟨T(ω)⟩ for parameter ranges of 0.1 ≤ M/m ≤ 10 and 0.01 ≤ C ≤ 0.5. In the present study, we calculated the single phonon transmission function, and the many-body interaction between phonons is neglected. Since the considered system length is very long, anharmonic phonon–phonon scattering^{14,25,34–40} would become relevant. It has been suggested that one-dimensional systems with linear dispersion exhibit anomalous diffusion even in the presence of anharmonic interactions.¹⁴ However, no solid proof has yet been obtained for linear and quadratic dispersion models. Such a study remains to be performed in the future.

This work was supported by the JSPS KAKENHI under Grant No. 20H00250.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In this appendix, we show that the quadratic dispersion model with the dynamical matrix Φ⁽²⁾ of Eq. (4) can be considered as a discretized version of the Euler–Bernoulli equation of an unloaded thin beam.

The dynamic beam equation for the transverse displacement, y(x, t), is given by

where κ(x) is the flexural rigidity and ρ(x) is the mass density. By applying the finite difference method on the grid points x_j = ja (j = …, −1, 0, 1, 2, …), Eq. (A1) can be written in the frequency domain as

By introducing a mass-normalized displacement defined as $uj=(aρj)1/2yj$⁠, Eq. (A2) becomes

where Φ_ij are defined as follows:

If we consider aρ_j as an atomic mass, M_j, of a one-dimensional atomic chain model and 4κ_j/a³ as a spring constant k_j, the dynamical matrix is written as

For a uniform spring constant case (k_j = k for all j), this equation reduces to Eq. (4).