Thermal conductance, *K*, of one-dimensional mass-disordered systems is calculated. Two types of dispersion models are considered: linear and quadratic dispersion model. For the linear dispersion model, the thermal conductance shows the super-diffusive nature, which is consistent with the previous research. For the quadratic dispersion model, the thermal conductance exhibits the normal-diffusion in which *K* is inversely proportional to the system size.

## I. INTRODUCTION

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.^{5} 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^{6} 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^{12} of *K* ∝ *N*^{−3/2}.

For a quadratic dispersion *ω*_{q} = *cq*^{2}/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.^{19} It is also important in this sense to investigate the effects of the quadratic dispersion on the thermal transport. We calculated^{20} 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^{9}) 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.

## II. THEORY

### A. Model

#### 1. Linear dispersion model

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 *n*th 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=(\u2026,u\u22122,u\u22121,u0)T$, $UC=(u1,u2,\u2026,uN)T$, and $UR=(uN+1,uN+2,uN+3,\u2026)T$ are the column vectors of the mass-normalized displacement *u*_{n} of the *n*th atom for the left reservoir, channel, and right reservoir, respectively. The mass-normalized dynamical matrix Φ^{(1)} is given by

with *d*_{n} = 2 *k*/*M*_{n}, $sn=\u2212k/(MnMn+1)1/2$ for *n* = 1, 2, …, *N* and *d*_{0} = 2 *k*/*m*, *s*_{0} = −*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.

#### 2. Quadratic dispersion model

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 Φ^{(2)} 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*.

### B. Dispersion relation

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 *ω*_{0} = (*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: $\omega q(1)=\omega q(2)=2\omega 0$ at *q* = *π*/*a*.

### C. Thermal conductance

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), $\Phi C(d)$ is the mass-normalized dynamical matrix of the channel, and $\Pi \alpha (d)(\omega )$ (*α* = L, R) are the contact *self-energies* (whose dimension is *ω*^{2}). The self-energies $\Pi \alpha (d)(\omega )$ 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, $\pi \alpha (d)(\omega )$, can be written as

for the linear dispersion model, where *q* is a solution of $\omega =\omega q(1)$, and

for the quadratic dispersion model, where *q* and *γ* are solutions of $\omega =\omega q(2)$ and *ω* = 2*ω*_{0} sinh^{2}(*γa*/2), respectively.

The phonon transmission function, *T*^{(d)}(*ω*), through the channel is given in terms of *D*^{(d)}(*ω*) and $\Pi \alpha (d)(\omega )$ as follows:

where $\Gamma \alpha (d)(\omega )$ 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;\omega )=[exp(\u210f\omega /kBT)\u22121]\u22121$ 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/\pi 2)(x2ex)/(ex\u22121)2$. In the present study, we consider the high temperature limit [*w*(*x*) → 3/*π*^{2}], 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.

### D. *R*-matrix method

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*^{3}, 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\alpha \beta (\omega )$ (*α*, *β* = L, R), defined as

We then define the *R*-matrix $R(\omega )\u2261R(\omega )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(\omega )$.

Compute a

*d*×*d*matrix of $r1(\omega )=[\omega 2\u2212D1]\u22121$, and define the*R*-matrix for the isolated first block with the elements $RLL=RLR=RRL=RRR=r1(\omega )$.- Repeat the iterations(28)$RRR=\omega 2\u2212Dk+1\u2212SkTRLLSk\u22121,$(29)$RRL=RRRSkTRRL,$(30)$RLL=RLL+RLRSnRRL,$for(31)$RLR=RRLT$
*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)].

## III. RESULTS AND DISCUSSION

### A. Linear dispersion model

We first present the results of the linear dispersion model. The channel length *N* dependence of the thermal conductance *K*^{(1)} 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*_{0} = *k*_{B}*ω*_{0}/*π*, 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*^{−1} 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.^{33} 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*^{(1)} 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*^{(1)} ∝ *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*^{(1)}(*ω*)⟩ 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^{8}

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

Matsuda and Ishii^{6} derived an analytical expression of *K*^{(1)} for *N* ≫ 1 based on the Kubo formula,

Figure 5 shows a comparison between the simulated *K*^{(1)} (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).

### B. Quadratic dispersion model

We next present the results on the quadratic dispersion model. Figure 6 shows the color plot of the sample-averaged transmission function ⟨*T*^{(2)}(*ω*)⟩ 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*^{−1}. By calculating the sample-averaged transmission functions ⟨*T*^{(2)}(*ω*)⟩ 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*^{(2)} in the quadratic dispersion model, and closed circles show *K*^{(1)} in the linear dispersion model. We find that *K*^{(2)} 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\pi )\omega d(1)$, we expect that an analytical formula for *K*^{(2)} 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.)

## IV. CONCLUSION

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^{9}. For the linear dispersion model, the thermal conductance exhibits the super-diffusive nature of *K*^{(1)} ∝ *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*^{(2)} ∝ *N*^{−1}. An analytical formula of the demarcation frequency $\omega 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.^{14} However, no solid proof has yet been obtained for linear and quadratic dispersion models. Such a study remains to be performed in the future.

## ACKNOWLEDGMENTS

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

## DATA AVAILABILITY

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

### APPENDIX: QUADRATIC DISPERSION MODEL

In this appendix, we show that the quadratic dispersion model with the dynamical matrix Φ^{(2)} 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\rho 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*^{3} 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).