Full counting statistics of phonon transport in disordered systems

Abstract

The coherent potential approximation (CPA) within full counting statistics (FCS) formalism is shown to be a suitable method to investigate average electric conductance, shot noise as well as higher order cumulants in disordered systems. We develop a similar FCS-CPA formalism for phonon transport through disordered systems. As a byproduct, we derive relations among coefficients of different phonon current cumulants. We apply the FCS-CPA method to investigate phonon transport properties of graphene systems in the presence of disorders. For binary disorders as well as Anderson disorders, we calculate up to the 8-th phonon transmission moments and demonstrate that the numerical results of the FCS-CPA method agree very well with that of the brute force method. The benchmark shows that the FCS-CPA method achieves 20 times more speedup ratio. Collective features of phonon current cumulants are also revealed.

Graphical abstract

Keywords

phonon transport / disordered systems / coherent potential approximation / full counting statistics

1 Introduction

Due to the advancement of nanotechnology, the feature size of electronic devices has reached nanoscale. Power dissipation of nano-devices increases as they become smaller which may limit device size and density. Hence it is important to study how power dissipation especially Joule heat is dissipated in nano-devices. Since phonon transport is one of the primary dissipation mechanisms, investigating phonon transport through nano-devices [1, 2] can obtain important information in order to design low power devices. From scientific point of view, phonon quantum transport has shown a number of interesting phenomena such as the quantized thermal conductance [3], breakdown of Fourier’s law in nanoscale systems [4], phonon filtering [5], and topological phononic crystals [6–8], to name just a few.

Impurities are always present in nano-devices giving rise to disorders for phonon transport, which has been studied for various disordered systems. For topological phononic systems, it was found that the topological edge state is intact in the presence of uncorrelated disorder and is gradually destroyed when the disorder is spatially correlated [9]. For one-dimensional harmonic chain with mass disorders, the effect of long-range interaction on phonon transport was examined and found to enhance (reduce) the transmission of high (low) frequency phonon [10]. To avoid huge computational burden of the brute force calculation of disorder average, the non-equilibrium vertex correction (NVC) theory [11, 12] within coherent potential approximation (CPA) [13, 14] was developed to perform analytic average over phonon transmission coefficient [15]. When both mass disorder and force constant disorder are present, it was found that the interplay of two disorders gives rise to an anomalous transparency in low dimensional phononic systems [16]. To release the computational burden, a CPA like theory was developed to deal with phonon transmission coefficient with both mass and force constant disorder present [17]. A dynamic cluster approximation formalism was also developed to investigate density of states (DOS) of periodic binary systems with both diagonal and off-diagonal disorders [18]. So far, the CPA formalism has been used to study DOS and phonon transmission coefficient which corresponds to taking disorder average over one and two Green’s functions.

In this work, we develop a CPA formalism based on the full counting statistics (FCS) framework that is capable of calculating disorder averaging over 2n Green’s functions needed for studying n-th cumulant of phonon current in the presence of disorders. The central idea of this formalism is to make a nonlinear transformation so that the problem of disorder averaging 2n Green’s function is mapped into that of averaging one generalized Green’s function with one additional parameter, i.e., the counting field [19]. After this transformation, it is purely a CPA problem and no NVC and other higher order NVC are needed. To obtain the n-th cumulant of phonon current, one performs n-th order numerical derivative respect to the counting field. We use this FCS-CPA formalism on graphene systems in the presence of binary disorders and Anderson disorders. The phonon transmission moments are calculated up to 8-th order using FCS-CPA approach. Compared with the brute force method, very good agreement is found. In addition, relations among various coefficients of phonon current cumulants are derived from the symmetry relation of the cumulant generating function.

The paper is organized as follows. In Section 2, a brief introduction of CPA method is given and the theoretical formalism of FCS-CPA for phonon transport is derived. Two methods of calculating the n-th moment of phonon transmission are given: (i) numerical derivatives with respect to the counting field; (ii) solving a series of NVC equations. Finally, relations among various phonon current cumulants are derived. In Section 3, the FCS-CPA formalism is applied to study phonon transport properties of a honeycomb-lattice system and numerical results are presented and analyzed. A short summary is given in Section 4.

2 Theoretical formalism

2.1 Coherent potential approximation

Before introducing the FCS-CPA method, it is necessary to briefly review the coherent potential approximation (CPA). The CPA method can be summarized as: given a complex matrix A\in {ℂ}^{nm\times nm}, and a block diagonal random matrix B of the same size with each block defined by some independent random distributions B_{ii}\sim X^{m\times m}, the average of the matrix inversion can be approximated by

\begin{array}{}\text{(1)}& \langle \frac{1}{A-B}\rangle =\frac{1}{A-C}\mathrm{,}\end{array}

\begin{array}{}\text{(2)}& C_{ij}=\langle B_{ii}{\left(I-{\left(\frac{1}{A-C}\right)}_{ii}\left(B_{ii}-C_{ii}\right)\right)}^{-1}\rangle {\delta }_{ij}\mathrm{,}\end{array}

where the matrix C is the solution of the above self-consistent equations. It should be noted that we use i\mathrm{,}j to label the matrix block, and two more subscripts are required if we want to retrieve one specific matrix element (e.g., B_{i\alpha \mathrm{,}j\beta }). Such an indexing convention is similar to the tight-binding hamiltonian where the system contains n sites and each site has several orbits (e.g., s/p/d/f-orbit for electron and x/y/z degrees of freedom for phonon).

To derive the above equations, we define G=(A-C{)}^{-1} which is called the renormalized Green’s function and introduce the T matrix such that {(A-B)}^{-1}=G+GTG. It is easy to find

\begin{array}{}\text{(3)}& \begin{array}{ll}T\hfill & =(B-C){\displaystyle \sum _{i=0}^{\infty }}{[G(B-C)]}^i\hfill \\ \hfill & =\left(B-C\right)\left(I+GT\right)\hfill \\ \hfill & =\left(I+TG\right)\left(B-C\right)\mathrm{.}\hfill \end{array}\end{array}

Obviously, \langle T\rangle =0 gives rise to Eq. (1) and corresponds to CPA. To facilitate the derivation, it is convenient to define the following three matrices B_i\mathrm{,}{C}_i\mathrm{,}{T}_i as

\begin{array}{}\text{(3)}& \begin{array}{l}{\left(B_i\right)}_{jk}=\{\begin{array}{l}{B}_{ii}\mathrm{,}j=k=i\mathrm{,}\hfill \\ \mathrm{0,}\text{else}\mathrm{,}\hfill \end{array}\\ {\left(C_i\right)}_{jk}=\{\begin{array}{l}{C}_{ik}\mathrm{,}j=i\mathrm{,}\hfill \\ \mathrm{0,}\text{else}\mathrm{,}\hfill \end{array}\\ {\left(T_i\right)}_{jk}=\{\begin{array}{l}{T}_{ik}\mathrm{,}j=i\mathrm{,}\hfill \\ \mathrm{0,}\text{else}\mathrm{,}\hfill \end{array}\end{array}\end{array}

such that {\displaystyle {\sum }_i}{B}_i=B, {\displaystyle {\sum }_i}{C}_i=C, and {\displaystyle {\sum }_i}{T}_i=T. It should be noted that all these matrices are of the same size as A and the subscript i here denotes that all rows except the i-th block row are set to zero. In general, the matrix C is a full matrix. But we will show below, under the single site approximation (SSA), the matrix C is structurally similar to B so that all off-diagonal blocks of C are zero. A visual representation of these matrices look like the following.

\begin{array}{}\text{(3)}& B_i\sim \left(\begin{array}{lll}\square \hfill & \square \hfill & \square \hfill \\ \square \hfill & ■\hfill & \square \hfill \\ \square \hfill & \square \hfill & \square \hfill \end{array}\right);C_i\text{ or }{T}_i\sim \left(\begin{array}{lll}\square \hfill & \square \hfill & \square \hfill \\ ■\hfill & ■\hfill & ■\hfill \\ \square \hfill & \square \hfill & \square \hfill \end{array}\right)\mathrm{.}\end{array}

From Eq. (3) we have

\begin{array}{}\text{(3)}& T_{ij}={\displaystyle \sum _k}{(B-C)}_{ik}{(I+GT)}_{kj}\mathrm{,}\end{array}

which is equivalent to

\begin{array}{}\text{(3)}& T_i=(B_i-C_i)(I+G{T}_i+G{\displaystyle \sum _{j\ne i}}{T}_j)\mathrm{.}\end{array}

Solving for T_i, we find

\begin{array}{}\text{(3)}& T_i=t_i+t_{i}G{\displaystyle \sum _{j\ne i}}{T}_j\mathrm{,}\end{array}

where

\begin{array}{}\text{(4)}& \begin{array}{ll}{t}_i\hfill & ={\left[I-(B_i-C_i)G\right]}^{-1}(B_i-C_i)\hfill \\ \hfill & =(B_i-C_i){\left[I-G(B_i-C_i)\right]}^{-1}\mathrm{.}\hfill \end{array}\end{array}

So far, no approximation has been made. Now we assume \langle t_{i}G{T}_j\rangle \approx \langle t_i\rangle G\langle T_j\rangle for i\ne j. This is the so-called single site approximation (SSA), which is good at weak disorders where multi-scattering can be neglected. Under this approximation, we have

\begin{array}{}\text{(4)}& \langle T\rangle ={\displaystyle \sum _i}\langle t_i\rangle +{\displaystyle \sum _{i\mathrm{,}j\ne i}}\langle t_i\rangle G\langle T_j\rangle \mathrm{.}\end{array}

From the above equation, we see that \langle t_i\rangle =0 is a sufficient condition for Eq. (1). From \langle t_i\rangle =0, we obtain a useful identity,

\begin{array}{}\text{(5)}& \begin{array}{ll}\langle t_i\rangle \hfill & =\langle {(I-(B_i-C_i)G)}^{-1}(B_i-C_i)\rangle \hfill \\ \hfill & =\langle {(I-(B_i-C_i)G)}^{-1}-I\rangle G^{-1}\mathrm{,}\hfill \end{array}\end{array}

or

\begin{array}{}\text{(6)}& \langle {\left(I-\left(B_i-C_i\right)G\right)}^{-1}\rangle =I\mathrm{.}\end{array}

Inserting this identity back to the first line of Eq. (5), we arrive at

\begin{array}{}\text{(6)}& C_i=\langle {\left(I-\left(B_i-C_i\right)G\right)}^{-1}{B}_i\rangle \mathrm{.}\end{array}

Since all columns except the i-th block of the rightmost matrix B_i are zero, we can deduce that the only nonzero block of C_i is on the diagonal and equivalently the matrix C is block diagonal. So the self-consistent equation for the i-th diagonal block of C can be written as

\begin{array}{}\text{(6)}& C_{ii}=\langle {\left(I-\left(B_{ii}-C_{ii}\right){\left(\frac{1}{A-C}\right)}_{ii}\right)}^{-1}{B}_{ii}\rangle \mathrm{.}\end{array}

Similarly, if we use the first line of Eq. (4), we will arrive at Eq. (2).

Eqs. (1) and (2) frequently appear in calculating transport properties of disordered systems. To evaluate the average of a matrix such as the left hand side of Eq. (1), the most intuitive way is the brute force (BF) method. The procedure of brute force calculations is as follows: (i) generating thousands of random samples of matrix B; (ii) calculating matrix inversions; (iii) averaging over the whole ensemble. Apparently, the brute force method is extremely time-consuming. On the contrary, in the CPA method, only a reasonable number of inverse operations are required to solve the self-consistent equations. Although SSA has been made in the derivation, we will show in the numerical results that the accuracy of the CPA simulation is comparable to that of the brute force method.

2.2 Tight-binding Hamiltonian and Green’s function formalism

The system we considered is shown in Fig. 1, which is a zigzag graphene nanoribbon and consists of three parts: the central scattering region, the left and right leads. Driven by the temperature difference between two leads, the phonon flows from the left lead to the right one. For convenience of discussion, we assume that the temperature in the left lead is always higher, T_L>T_R. Since both leads are of semi-infinite size, which are also called phonon reservoirs, a steady phonon current is established in the long time limit. For phonon transport, the Hamiltonian is expressed as

\begin{array}{}\text{(6)}& H=H_C+H_L+H_R+H_T\mathrm{,}\end{array}

\begin{array}{}\text{(6)}& H_C={\displaystyle \sum _{i\in C}}\frac{p_i^2}{2m_i}+{\displaystyle \sum _{i\mathrm{,}j\in C}}\frac{1}{2}{K}_{ij}{x}_i{x}_j\mathrm{,}\end{array}

\begin{array}{}\text{(6)}& H_L={\displaystyle \sum _{k\in L}}\frac{p_k^2}{2m_k}+{\displaystyle \sum _{k\mathrm{,}{k}^{\prime }\in L}}\frac{1}{2}{K}_{k{k}^{\prime }}{x}_k{x}_{k^{\prime }}\mathrm{,}\end{array}

\begin{array}{}\text{(6)}& H_R={\displaystyle \sum _{k\in R}}\frac{p_k^2}{2m_k}+{\displaystyle \sum _{k\mathrm{,}{k}^{\prime }\in R}}\frac{1}{2}{K}_{k{k}^{\prime }}{x}_k{x}_{k^{\prime }}\mathrm{,}\end{array}

\begin{array}{}\text{(6)}& \begin{array}{ll}{H}_T=\hfill & \frac{1}{2}{\displaystyle \sum _{i\in C\mathrm{,}k\in L}}(K_{ik}{x}_i{x}_k+K_{ki}{x}_k{x}_i)\hfill \\ \hfill & +\frac{1}{2}{\displaystyle \sum _{i\in C\mathrm{,}k\in R}}(K_{ik}{x}_i{x}_k+K_{ki}{x}_k{x}_i)\mathrm{.}\hfill \end{array}\end{array}

The above four terms are Hamiltonians of the central region H_C, the left lead H_L, the right lead H_R, and the coupling H_T between the central region and leads, respectively. The operator x_i describes the displacement from the equilibrium position of the i-th atom and p_i is the corresponding momentum operator. x_i and p_j satisfy the commutation relation, [x_i\mathrm{,}{p}_j]=i\hslash {\delta }_{ij}. The mass term m_i is usually absorbed into the displacement operator x_i for simplification. But this strategy is not available here, since we will treat the mass term as a random variable in disordered systems. We adopt the force constant matrix K from Ref. [20] and only consider the first nearest neighbor coupling. The force constant between the neighboring atoms is a 3×3 matrix since each atom has three degrees of freedom (x/y/z), which can be written as

\begin{array}{}\text{(6)}& K_{ij}=\left(\begin{array}{lll}-t_s{{\displaystyle \cos }}^2\theta -t_i{{\displaystyle \sin }}^2\theta \hfill & \left(t_i-t_s\right)\cos \theta \sin \theta \hfill & \hfill \\ \left(t_i-t_s\right)\cos \theta \sin \theta \hfill & -t_s{{\displaystyle \sin }}^2\theta -t_i{{\displaystyle \cos }}^2\theta \hfill & \hfill \\ \hfill & \hfill & t_o\hfill \end{array}\right)\mathrm{,}\end{array}

where the stretching factor t_s, the in-plane bending factor t_i, and the out-of-plane bending factor t_o are obtained by fitting the phonon dispersion relations from experimental data [21]. To apply the acoustic sum rule, we subtract the summations of those force constants from the diagonal elements as

\begin{array}{}\text{(7)}& K_{ii}=-{\displaystyle \sum _{j\ne i\mathrm{,}\alpha =x\mathrm{,}y\mathrm{,}z}}\text{diag}\{K_{ixj\alpha }\mathrm{,}{K}_{iyj\alpha }\mathrm{,}{K}_{izj\alpha }\}\mathrm{.}\end{array}

In terms of the non-equilibrium Green’s function (NEGF), the phonon current is given by the Landauer formula [22–25]

\begin{array}{}\text{(8)}& J_E={\displaystyle {\int }_0^{\infty }}\frac{\text{d}\omega }{2\text{π}}\text{Tr}[\widehat{T}](f_L-f_R)\hslash \omega \mathrm{,}\end{array}

where \widehat{T} is the transmission matrix given by

\begin{array}{}\text{(8)}& \widehat{T}={\mathrm{\Gamma }}_L{G}^r{\mathrm{\Gamma }}_R{G}^a\mathrm{,}\end{array}

and other quantities are defined as

\begin{array}{}\text{(9)}& \begin{array}{l}{G}^r=\underset{\eta \to 0^{+}}{{\displaystyle \lim }}\frac{I}{M{\left(w+\mathrm{i}\eta \right)}^2-K_C-{\mathrm{\Sigma }}_L^r-{\mathrm{\Sigma }}_R^r}\mathrm{,}\\ G^a={\left(G^r\right)}^{†}\mathrm{,}\\ M=\text{diag}\{m_1\mathrm{,}{m}_1\mathrm{,}{m}_1\mathrm{,}\cdots \mathrm{,}{m}_N\mathrm{,}{m}_N\mathrm{,}{m}_N\}\mathrm{.}\end{array}\end{array}

In Eq. (8), f_{L/R}=[{\text{e}}^{\hslash \omega /(k_B{T}_{L/R})}-{1]}^{-1} is the Bose–Einstein distribution of the left/right lead, where k_B is the Boltzmann constant, \hslash is the Plank constant, and ω is the angular frequency of phonon modes. K_C is the force constant matrix of the central region. M is a diagonal matrix and its element m_i represents the mass of the i-th atom. The same three m_is in M corresponds to three space degrees of freedoms of each atom. For an open system, the retarded (advanced) Green’s function G^r (G^a) for the central scattering region can be obtained by absorbing leads’ contribution via the self-energy ({\mathrm{\Sigma }}_L^r and {\mathrm{\Sigma }}_R^r) and these self-energies can be numerically calculated using the recursive Green’s function algorithm [26]. Also, the linewidth function {\mathrm{\Gamma }}_{L/R}=\text{i}({\mathrm{\Sigma }}_{L/R}^r-{\mathrm{\Sigma }}_{L/R}^a) gives the information about the propagating wavefunction in the lead. The phonon current can be easily understood as follows. Firstly, phonons are propagating from one lead to the other side and only the phonon population difference between two leads f_L-f_R makes the net contribution to the phonon current. Secondly, the transmission coefficient \text{Tr}(\widehat{T}) describes the propagation probability for the phonon with a specific energy, and the larger \text{Tr}(\widehat{T}) is, the larger the phonon current will be. Lastly, the phonon in all frequency domains will contribute to the current, therefore the integration with respect to the frequency gives rise to the total phonon current.

For phonon transport in disordered systems, we assume that the disorder only exists in the central scattering region. A commonly used disorder type in phonon systems is the isotopic disorder, where the mass of each atom becomes a random variable. If two different isotopes C₁ and C₂ of the same element exist in the system, the mass of the i-th atom m_i is random and has the following distribution,

\begin{array}{}\text{(10)}& p(m_i)=\{\begin{array}{l}p\mathrm{,}{m}_i=C_1\mathrm{,}\hfill \\ 1-p\mathrm{,}{m}_i=C_2\mathrm{.}\hfill \end{array}\end{array}

Therefore, the mass matrix M becomes a diagonal random matrix.

In the following, we will mainly focus on the phonon transmission coefficient \text{Tr}(\widehat{T}) since it contains all disorder effects in the system. Theoretically, the exact average transmission coefficient \langle \text{Tr}(\widehat{T})\rangle can be evaluated by enumerating all possible configurations,

\begin{array}{}\text{(11)}& \langle \text{Tr}(\widehat{T})\rangle =\text{}\text{}{\displaystyle \sum _{m_1\mathrm{,}{m}_2\mathrm{,}\cdots \mathrm{,}{m}_N}}\text{}\text{}\text{Tr}(\widehat{T})(m_1\mathrm{,}{m}_2\mathrm{,}\cdots \mathrm{,}{m}_N){\displaystyle \prod _i}p(m_i)\mathrm{.}\end{array}

Obviously, such a method is infeasible for large systems, since it requires 2^N samples for a N-atom system with binary isotopic diorders. A more practical way is the Monte Carlo method as shown in Eq. (12), which we refer as the brute force (BF) method throughout the paper,

\begin{array}{}\text{(12)}& \langle \text{Tr}(\widehat{T})\rangle \approx \frac{1}{K}{\displaystyle \sum _{i=1}^K}\text{Tr}(\widehat{T})({\text{sample}}_i\{m_1\mathrm{,}{m}_2\mathrm{,}\cdots \mathrm{,}{m}_N\})\mathrm{.}\end{array}

Here K is the number of random configurations which is far more less than 2^N. Different from the exact one, the brute force method only takes average over reasonable amount of random samples and still gives a quite accurate result.

2.3 Phonon transmission moments

In the previous subsection, we introduced average phonon transport properties of disordered systems. Since the phonon current alone is not enough to fully characterize phonon transport properties, we need higher order cumulants or moments of the phonon current operator and the average value is simply the first order cumulant. Full counting statistics (FCS), which has been widely used in transport study such as energy current [27, 28] and transient dynamics [29, 30], is suitable for evaluating high-order cumulants and moments. Different from the previous work [19] which focused on the transmission coefficient cumulants, we will show that the transmission moments will give more precise results.

We start with the cumulant generating function [31] (CGF) \ln \mathcal{Z} as shown in Eq. (13),

\begin{array}{}\text{(13)}& \begin{array}{l}\ln \mathcal{Z}={\displaystyle {\int }_0^{\infty }}\frac{\text{d}\omega }{2\text{π}}\mathcal{X}\mathrm{,}\\ \mathcal{X}=-\text{Trln}(I-\widehat{T}Y)\mathrm{,}\\ Y=({\text{e}}^{\text{i}\lambda \hslash \omega }-1)f_L(1+f_R)+({\text{e}}^{-\text{i}\lambda \hslash \omega }-1)f_R(1+f_L)\mathrm{.}\end{array}\end{array}

The n-th cumulant of the phonon current {\mathcal{C}}_n is obtained as the n-th derivative of the CGF with respect to \text{i}\lambda at λ=0,

\begin{array}{}\text{(14)}& {\mathcal{C}}_n={\frac{{\partial }^n\ln \mathcal{Z}}{{(\partial \text{i}\lambda )}^n}|}_{\lambda =0}={\displaystyle {\int }_0^{\infty }}\frac{\text{d}\omega }{2\text{π}}{\frac{{\partial }^n\mathcal{X}}{{(\partial \text{i}\lambda )}^n}|}_{\lambda =0}\mathrm{.}\end{array}

Expanding \mathcal{X} in terms of λ and taking the derivative, we find

\begin{array}{}\text{(15)}& {\frac{{\partial }^n\mathcal{X}}{{(\partial \text{i}\lambda )}^n}|}_{\lambda =0}={\displaystyle \sum _{m=1}}\frac{1}{m}\text{Tr}({\widehat{T}}^m)Y_{mn}\mathrm{,}\end{array}

where

\begin{array}{}\text{(15)}& Y_{mn}={\frac{{\partial }^n{Y}^m}{{(\partial \text{i}\lambda )}^n}|}_{\lambda =0}\mathrm{.}\end{array}

Eq. (15) can be used to perform brute force calculations of \langle {\mathcal{C}}_n\rangle where Y_{mn} does not depend on disorders. We now formulate the FCS-CPA method and use it to calculate \langle {\mathcal{C}}_n\rangle or equivalently \langle \text{Tr}({\widehat{T}}^m)\rangle. From Eq. (15) we see that the calculation of \langle {\mathcal{C}}_n\rangle amounts to calculate disorder average of multi-Green’s functions. The key advantage of FCS-CPA is to transform the representation of multi-Green’s functions to the representation of single generalized Green’s function so that CPA can be directly applied [32]. Specifically, one can transform the CGF so that it contains only one generalized Green’s function. We first use the following relation with the fact that \ln \text{Det}G=\text{Tr}\ln G

\begin{array}{}\text{(16)}& \mathcal{X}=\ln \text{Det}(I+\tilde{G}{\mathrm{\Gamma }}_0\alpha )=\text{Tr}\ln (I+\tilde{G}{\mathrm{\Gamma }}_0\alpha )\mathrm{,}\end{array}

where \alpha =\sqrt{Y} and {\mathrm{\Gamma }}_0=\left(\begin{array}{ll}\hfill & {\mathrm{\Gamma }}_R\hfill \\ {\mathrm{\Gamma }}_L\hfill & \hfill \end{array}\right)\mathrm{,}\tilde{G}=\left(\begin{array}{ll}{G}^r\hfill & \hfill \\ \hfill & G^a\hfill \end{array}\right)\mathrm{.}

Writing Eq. (16) in an integral form, we have

\begin{array}{}\text{(16)}& \begin{array}{ll}\mathcal{X}\hfill & =-\text{Tr}{\displaystyle {\int }_0^{\alpha }}\text{d}x\frac{1}{{\mathrm{\Gamma }}_{0}x+{\tilde{G}}^{-1}}{\mathrm{\Gamma }}_0\hfill \\ \hfill & =-\text{Tr}{\displaystyle {\int }_0^{\alpha }}\text{d}x\frac{1}{{\mathrm{\Gamma }}_{0}x+G^{-1}-B}{\mathrm{\Gamma }}_0\mathrm{,}\hfill \end{array}\end{array}

where

\begin{array}{}\text{(17)}& \begin{array}{l}{G}^{-1}=-\left(\begin{array}{ll}{K}_C+{\mathrm{\Sigma }}_L^r+{\mathrm{\Sigma }}_R^r\hfill & \hfill \\ \hfill & K_C+{\mathrm{\Sigma }}_L^a+{\mathrm{\Sigma }}_R^a\hfill \end{array}\right)\mathrm{,}\\ B=-\left(\begin{array}{ll}M{(\omega +\text{i}\eta )}^2\hfill & \hfill \\ \hfill & M{(\omega -\text{i}\eta )}^2\hfill \end{array}\right)\mathrm{.}\end{array}\end{array}

The denominator has been divided into two parts, one is a constant matrix ({\mathrm{\Gamma }}_{0}x+G^{-1}) and the other is a random diagonal matrix B. Then we can apply the CPA to obtain the average quantity

\begin{array}{}\text{(18)}& \langle \mathcal{X}\rangle =-\text{Tr}{\displaystyle {\int }_0^{\alpha }}\text{d}x\frac{1}{{\mathrm{\Gamma }}_{0}x+G^{-1}-\mathrm{\Delta }(x)}{\mathrm{\Gamma }}_0\mathrm{,}\end{array}

\begin{array}{}\text{(19)}& \begin{array}{ll}{\mathrm{\Delta }}_{ii}(x)\hfill & =\langle B_{ii}\{I\hfill \\ \hfill & -{\left[\frac{1}{{\mathrm{\Gamma }}_{0}x+G^{-1}-\mathrm{\Delta }\left(x\right)}\right]}_{ii}[B_{ii}-{\mathrm{\Delta }}_{ii}(x)]{\}}^{-1}\rangle \mathrm{.}\hfill \end{array}\end{array}

We note that {\mathrm{\Delta }}_{ii}(x) is a 6×6 matrix not 3×3, since the random matrix M is duplicated in Eq. (17). In order to obtain various cumulants, we need to take derivatives of Eq. (18). The first derivative of \langle \mathcal{X}\rangle with respect to \text{i}\lambda gives

\begin{array}{}\text{(20)}& \frac{\partial \langle \mathcal{X}\rangle }{\partial \text{i}\lambda }=-\frac{\partial \alpha }{\partial \text{i}\lambda }\text{Tr}\left[\frac{1}{{\mathrm{\Gamma }}_0\alpha +G^{-1}-\mathrm{\Delta }(\alpha )}{\mathrm{\Gamma }}_0\right]\mathrm{.}\end{array}

Apparently, the integration is avoided. Therefore, a direct way of calculating the n-th cumulant of phonon current is by taking (n-1)-th derivative on Eq. (20) numerically with small enough λ. Alternatively, we can apply the Taylor expansion on Eq. (19) and derive the following equations

\begin{array}{}\text{(21)}& \mathrm{\Delta }(x)={\mathrm{\Delta }}_0+{\displaystyle \sum _{k=1}^{\infty }}{\mathrm{\Delta }}_k{x}^k\mathrm{,}\end{array}

\begin{array}{}\text{(22)}& \frac{1}{{\mathrm{\Gamma }}_{0}x+G^{-1}-\mathrm{\Delta }\left(x\right)}=N_0+{\displaystyle \sum _{k=1}^{\infty }}{N}_k{x}^k\mathrm{,}\end{array}

\begin{array}{}\text{(23)}& \begin{array}{l}{\left\{I-{\left[\frac{1}{{\mathrm{\Gamma }}_{0}x+G^{-1}-\mathrm{\Delta }\left(x\right)}\right]}_{ii}[B_{ii}-{\mathrm{\Delta }}_{ii}\left(x\right)]\right\}}^{-1}\\ ={\left(K_0\right)}_{ii}+{\displaystyle \sum _{k=1}^{\infty }}{\left(K_k\right)}_{ii}{x}^k\mathrm{.}\end{array}\end{array}

Here matrices N_k and K_k are Taylor expansion coefficients and detailed expressions are shown in Appendix A. Inserting Eq. (21) and Eq. (23) into Eq. (19), we can get

\begin{array}{}\text{(24)}& {({\mathrm{\Delta }}_k)}_{ii}=\langle B_{ii}{(K_k)}_{ii}\rangle \mathrm{.}\end{array}

When k=0, Eq. (24) can be decoupled into two self-consistent CPA equations of block size 3 and they are conjugate to each other, which means only one self-consistent CPA equation is needed to solve. Such a property can be easily derived from Eq. (19) by setting x=0. When k=1, Eq. (24) actually corresponds to the NVC equation [19] which is necessary to calculate the average transmission coefficient. For k=2, Eq. (24) is solved to find the average second cumulant of the phonon current when {\mathrm{\Delta }}_0 and {\mathrm{\Delta }}_1 are provided. Hence for k>1, Eq. (24) will be referred as the higher order NVC equation. Once {\mathrm{\Delta }}_0\mathrm{,}{\mathrm{\Delta }}_1\mathrm{,}\cdots \mathrm{,}{\mathrm{\Delta }}_{k-1} have been established, Eq. (24) becomes linear with respect to {\mathrm{\Delta }}_k (see Appendix A) and only one more matrix inversion is required to solve for {\mathrm{\Delta }}_k. With these Taylor expansion coefficients, we can carry out the derivatives on Eq. (18) with respect to \text{i}\lambda

\begin{array}{}\text{(25)}& \begin{array}{ll}{\frac{{\partial }^n\langle \mathcal{X}\rangle }{{(\partial \text{i}\lambda )}^n}|}_{\lambda =0}\hfill & =-{\frac{{\partial }^n}{{(\partial \text{i}\lambda )}^n}|}_{\lambda =0}\left[{\displaystyle \sum _{k=0}^{\infty }}\frac{1}{k\text{}+\text{}1}\text{Tr}(N_k{\mathrm{\Gamma }}_0){\alpha }^{k\text{}+\text{}1}\right]\hfill \\ \hfill & =-{\displaystyle \sum _{m=1}}\frac{1}{2m}\text{Tr}(N_{2m-1}{\mathrm{\Gamma }}_0)Y_{mn}\mathrm{.}\hfill \end{array}\end{array}

In the second line, we drop the even terms of N_{2m} since the trace of its product with {\mathrm{\Gamma }}_0 is equal to zero. Compared with Eq. (15), we arrive at

\begin{array}{}\text{(26)}& \langle \text{Tr}({\widehat{T}}^m)\rangle =-\text{Tr}(N_{2m-1}{\mathrm{\Gamma }}_0)/2.\end{array}

Eq. (26) provides us two ways to calculate phonon transmission moments. The left hand side of Eq. (26) is for the brute force method via averaging over thousands of random samples, and the right hand side corresponds to the FCS-CPA method which is efficient and time-saving.

2.4 Relations among coefficients of phonon current cumulants

In the previous subsection, we have introduced phonon current cumulants defined as the derivative of the cumulant generating function (CGF) with respect to \text{i}\lambda. In terms of v=1/(k_B{T}_L)-1/(k_B{T}_R), we define the coefficient of phonon current cumulants as

\begin{array}{}\text{(27)}& C_n={\frac{{\partial }^n\mathcal{X}}{{(\partial \text{i}\lambda )}^n}|}_{\lambda =0}={\displaystyle \sum _{m=0}^{\infty }}\frac{v^m}{m!}{C}_n^{(m)}\mathrm{.}\end{array}

In this subsection, we will examine the relations among various coefficients C_n^{(m)}. We start from the symmetry relation of the CGF

\begin{array}{}\text{(28)}& \mathcal{X}[-\text{i}\lambda \mathrm{,}v]=\mathcal{X}[\text{i}\lambda +v\mathrm{,}v]\mathrm{.}\end{array}

Taking the derivative of Eq. (28) with respect to v, we have

\begin{array}{}\text{(29)}& \frac{\partial \mathcal{X}[-\text{i}\lambda \mathrm{,}v]}{\partial v}=\frac{\partial \mathcal{X}[\text{i}\lambda +v\mathrm{,}v]}{\partial v}+\frac{\partial \mathcal{X}[\text{i}\lambda +v\mathrm{,}v]}{\partial \text{i}\lambda }\mathrm{.}\end{array}

From the definition of cumulants and Eq. (27), we have the Taylor expansion of \mathcal{X} in terms of both \text{i}\lambda and v,

\begin{array}{}\text{(30)}& \mathcal{X}={\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \sum _{m=0}^{\infty }}\frac{{\left(\text{i}\lambda \right)}^n}{n!}\frac{v^m}{m!}{C}_n^{\left(m\right)}\mathrm{.}\end{array}

Inserting Eq. (30) into Eq. (29), the left hand side of Eq. (29) becomes

\begin{array}{}\text{(31)}& \text{L}\text{.H}\text{.S}\mathrm{.}={\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \sum _{m=0}^{\infty }}\frac{{\left(\text{i}\lambda \right)}^n}{n!}\frac{v^m}{m!}{\left(-1\right)}^n{C}_n^{\left(m+1\right)}\mathrm{.}\end{array}

After rearranging the summation order (see Appendix B), the right hand side of Eq. (29) becomes

\begin{array}{}\text{(32)}& \begin{array}{ll}\text{R}\text{.H}\text{.S}\text{.}\hfill & ={\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \sum _{m=0}^{\infty }}{\displaystyle \sum _{r=0}^n}\frac{{\left(\text{i}\lambda \right)}^r}{n!}\left(\begin{array}{c}n\\ r\end{array}\right)\frac{v^{m\text{}+\text{}n\text{}-\text{}r}}{m!}(C_n^{\left(m\text{}+\text{}1\right)}\text{}+\text{}{C}_{n+1}^{\left(m\right)})\hfill \\ \hfill & ={\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \sum _{m=0}^{\infty }}{\displaystyle \sum _{s=0}^{m+1}}\frac{{\left(\text{i}\lambda \right)}^n}{n!}\frac{v^m}{m!}\left(\begin{array}{c}m+1\\ s\end{array}\right)C_{n+s}^{\left(m+1-s\right)}.\hfill \end{array}\end{array}

Comparing Eq. (31) with Eq. (32), their corresponding coefficients should be equal which gives

\begin{array}{}\text{(33)}& {\left(-1\right)}^n{C}_n^{\left(m+1\right)}={\displaystyle \sum _{s=0}^{m+1}}\left(\begin{array}{c}m+1\\ s\end{array}\right)C_{n+s}^{\left(m+1-s\right)}\mathrm{,}\end{array}

where m=\mathrm{0,1,2,}\cdots.

For m=0, we have

\begin{array}{}\text{(34)}& \{\begin{array}{l}{C}_{2n-1}^{\left(0\right)}=\mathrm{0,}\hfill \\ 2C_{2n-1}^{\left(1\right)}+C_{2n}^{\left(0\right)}=0.\hfill \end{array}\end{array}

From Eq. (34) we see that C_{2n-1} contains no equilibrium contribution and vanishes when the temperature gradient is zero. In contrast, C_{2n} has both contributions from equilibrium and non-equilibrium processes. The results for m=\mathrm{1,2,3,4,5} are listed here,

\begin{array}{}\text{(35)}& \begin{array}{ll}m=1:\hfill & C_{2n-1}^{\left(2\right)}+C_{2n}^{\left(1\right)}=\mathrm{0,}\hfill \\ m=2:\hfill & 2C_{2n-1}^{\left(3\right)}+3C_{2n}^{\left(2\right)}+C_{2n+1}^{\left(1\right)}=\mathrm{0,}\hfill \\ m=3:\hfill & C_{2n-1}^{\left(4\right)}+2C_{2n}^{\left(3\right)}+C_{2n+1}^{\left(2\right)}=\mathrm{0,}\hfill \\ m=4:\hfill & 2C_{2n-1}^{\left(5\right)}+5C_{2n}^{\left(4\right)}+4C_{2n+1}^{\left(3\right)}+C_{2n+2}^{\left(2\right)}=\mathrm{0,}\hfill \\ m=5:\hfill & C_{2n-1}^{\left(6\right)}+3C_{2n}^{\left(5\right)}+3C_{2n+1}^{\left(4\right)}+C_{2n+2}^{\left(3\right)}=0.\hfill \end{array}\end{array}

These five equations can be summarized as

\begin{array}{}\text{(35)}& \begin{array}{ll}m=2r-1:\hfill & {\displaystyle \sum _{i=0}^r}\left(\begin{array}{c}r\\ i\end{array}\right)C_{2n-1+i}^{\left(2r-i\right)}=\mathrm{0,}\hfill \\ m=2r:\hfill & {\displaystyle \sum _{i=0}^{r+1}}\left\{\begin{array}{c}r+1\\ i\end{array}\right\}{C}_{2n-1+i}^{\left(2r+1-i\right)}=\mathrm{0,}\hfill \end{array}\end{array}

where \left\{\begin{array}{c}r+1\\ i\end{array}\right\} are coefficients of xⁱ in the Taylor expansion of (2x+1)(x+{1)}^r given by

Using these relations, we can express the coefficients of even cumulants C_{2n} in terms of coefficients of odd cumulants C_{2n-1}.

3 Numberical resutls

3.1 Phonon transmission moments

With the developed FCS-CPA method for phonon transport, we numerically study phonon transmission moments and cumulants of the graphene system shown in Fig. 1 to demonstrate the accuracy and efficiency of this method. In the calculation, a 4×4 honeycomb lattice has 16 atoms in the central scattering region, while a larger 16×16 lattice hosts 256 sites. The force constant parameters of graphene are adopted from Ref. [21] that the stretching factor t_s=365.0 J/m², the in-plane bending factor t_i=245.0 J/m², and the out-of-plane bending factor t_o=98.2 J/m². The brute force (BF) method, used as a comparison, always takes average over 1000 configurations.

First, we calculate the phonon transmission moments using both BF and FCS-CPA methods for a small 4×4 honeycomb lattice system. The atomic mass in left and right leads is 12 while the mass in the central region is chosen to be 12 or 14 with equal probability to simulate two carbon isotopes, which is called the binary disorder. The first four orders of transmission moments by two methods are presented in Fig. 2 with FCS-CPA results plotted in solid line while BF results in dots. We use various colors to represent different order of moments. For BF results, the estimated standard deviation is plotted as the error-bar in the figure. The transmission coefficients for the system without disorders are also plotted for reference, which shows clear step behavior. For such a small lattice, its transmission spectrum is quite flat and simple, so we can investigate the disorder influence closely. The small range of error-bar in the figure indicates that the BF method achieve a good convergence. Also, the solid lines go through the dotted data quite well, which suggests that the FCS-CPA results are in good agreement with the BF results in this frequency range. Because the CPA method is only valid in the weak disorder regime while the effective disorder strength in the phonon system \mathrm{\Delta }m{\omega }^2 grows rapidly, the FCS-CPA results at higher angular frequency, which is not plotted in the figure, are found to have large error or even diverge. Similar issue is also observed in electron systems. Compared with the quantized transmission in the pure system, all average transmission moments with disorders become much smoother. Transmission moment \langle \text{Tr}({\widehat{T}}^m)\rangle decreases when m increases due to the fact that all transmission eigenvalue are less than one, except the near-zero frequency region where transmission moments stay at the value of 3 which is enforced by the acoustic sum rule.

Then, we conduct the same calculation on a binary-disordered 16×16 lattice system and the results are shown in Fig. 3. Similar to the 4×4 system, the FCS-CPA results are in agreement with the BF results in the wide range of phonon frequency. One of the major differences between the disordered electronic and phononic systems is that the disorder strength for the electron system does not depend on the electron energy, while for the phonon system it is proportional to the square of phonon energy (equivalently the frequency). To demonstrate such difference, we calculate the transmission moments of a 16×16 lattice system with Anderson type disorder. The atomic mass in two leads is 12 as before while in the central scattering region the mass takes value from a uniform distribution. The first transmission moment \langle \text{Tr}(\widehat{T})\rangle for different disorder distributions are shown in Fig. 4. Uniform distributions U(\mathrm{11.8,12.2}), U(\mathrm{11,13}) and U(\mathrm{9,15}) with disorder strength \mathrm{\Delta }m=\mathrm{0.2,1,3}, respectively, are shown in the orange, gray and black lines. Also, the transmission coefficient for the system without disorder is plotted in blue line. The regions where the moments with different disorder strength deviate from the pure transmission spectrum with similar amount are marked with ellipses which are near \omega =\mathrm{40,70,150}, respectively. For these three regions (or three different disorder strength), the quantity \mathrm{\Delta }m{\omega }^2=\mathrm{4500,4900,4800} are close suggesting that the effective disorder strength is indeed proportional to \mathrm{\Delta }m{\omega }^2. Numerical results on the second transmission moment \langle \text{Tr}({\widehat{T}}^2)\rangle shows similar behaviors.

To benchmark the performance of the FCS-CPA method, we calculate transmission moments of a much larger system. For the phonon transport, we use a 48×48 lattice model and other parameters are kept the same as above. For such a system, the matrix involved in the FCS-CPA calculation is of order 20736 when solving Eq. (24). Since each site has only one orbit, a larger lattice 100×100 is used to benchmark for the electronic system and other parameters are adopted from the previous work [19]. The benchmark is conducted on a workstation of Intel Xeon Gold 5118 CPU with 4 cores and processor frequency 2.30 GHz. The results are listed in Table 1 where the time usage is for one frequency point (energy point for electron). We see that the FCS-CPA method achieves a great speedup ratio, around 20 times for both phonon and electron systems. The speedup ratio could be even larger if we consider the following aspects:

1) Sampling 1000 configurations is not accurate enough for the BF method especially for the second and higher order moments and the previous work [19] suggests that at least 10 thousands samples should be used.

2) When calculating transmission moments for a series of neighboring frequency points, the solved {\mathrm{\Delta }}_0 from the previous phonon frequency could be a good starting point for the next frequency calculation instead of using a random initial value as is done in this benchmark.

3.2 Phonon current cumulants

Once transmission moments are obtained efficiently using the FCS-CPA method, we can evaluate phonon current cumulants by integrating with respect to the frequency through Eq. (14), and numerical results are shown in this subsection.

First, we investigate the influence of binary disorders on phonon current cumulants. The calculation is performed on a 16×16 graphene lattice with binary isotopic disorder. The atomic mass in two leads is set to m=12, while in the central region two isotopes with masses m=12 and m=14 exist with different probabilities, which follows the distribution in Eq. (10). We set the temperatures in two leads as T_L=101 K and T_R=99 K. Results of the first eight cumulants are shown in Fig. 5. The probability of the m=12 isotope is labeled as the concentration for m=12 in the horizontal axis. In the figure, all even cumulants plotted are multiplied by a constant a=\frac{\mathrm{\Delta }T}{2k_B{\overline{T}}^2}=1.16 eV^–1 which can be deduced from Eq. (34). Here \mathrm{\Delta }T=T_L-T_R is the temperature difference and \overline{T}=(T_L+T_R)/2 is the average temperature. As expected in Section 2.4, the odd order cumulants are closely matched with the even order cumulants. One can see that, all cumulants decrease as the concentration for m=12 increases, reach their minima around concentration 0.4, and grow rapidly with further increasing of the m=12 concentration. Such phenomena also appear in the 12×12 lattice system, suggesting a common behavior.

Since all cumulants exhibit similar shapes, we reprocess the odd cumulants in Fig. 5 and plot the scaled results in Fig. 6. To investigate the underlying connection between these cumulants, we apply the first order polynomial fitting {\mathcal{C}}_i^{\prime }=a_i{\mathcal{C}}_i+b_i to find the best approximation. Higher order cumulants are first normalized to the range [0,1] and then linearly fitted to the first cumulant {\mathcal{C}}_1. The fitting results are shown in Fig. 6 with different colors and coefficients a_i and b_i are displayed in the inset table. Clearly, scaled odd cumulants {\mathcal{C}}_3\prime, {\mathcal{C}}_5\prime, and {\mathcal{C}}_7\prime are all in good agreement with each other and their corresponding curves almost overlap. They are also close to the first cumulant {\mathcal{C}}_1 and slight deviation appears only in the small concentration region around [0,0.1]. Due to the close relation between cumulant coefficients C_{2n} and C_{2n-1}, one can expect that even cumulants {\mathcal{C}}_4, {\mathcal{C}}_6, and {\mathcal{C}}_8 have the same linear connection. If we accept the tolerant error in this linear fitting, one can calculate only two points of higher order cumulants with respect to the concentration, and apply curve fitting to obtain other data points. This approximation provides us an efficient way to calculate high order phonon current cumulants.

Then, we study the effect of Anderson disorders on phonon current cumulants. The system under investigation is a 16×16 zigzag graphene nanoribbon with random atomic mass. The randomness follows the Anderson disorder where the mass m is uniformly distributed in U(12-\mathrm{\Delta }m\mathrm{,12}+\mathrm{\Delta }m) and \mathrm{\Delta }m is referred as the mass disorder strength. Other system parameters are kept unchanged. Cumulants up to the 8-th order are shown in Fig. 7. Different from the binary-disorder case where cumulants show parabolic behavior, all cumulants monotonically decrease with the increasing of \mathrm{\Delta }m in the presence of Anderson disorder. Similarly, we see that the odd cumulant {\mathcal{C}}_{2n-1} matches the even cumulant {\mathcal{C}}_{2n}. Since all cumulants have close shapes, again we use linear fitting {\mathcal{C}}_i^{\prime }=a_i{\mathcal{C}}_i+b_i to reveal their collective features. Odd cumulants from Fig. 7 are processed and numerical results are demonstrated in Fig. 8. One can find that the scaled odd cumulants {\mathcal{C}}_3\prime, {\mathcal{C}}_5\prime, and {\mathcal{C}}_7\prime are close to {\mathcal{C}}_1 in general, but the similarity is not good as the binary-disorder situation. Small deviations between {\mathcal{C}}_i\prime \text{s} exist in the whole disorder strength range, especially at the beginning and ending regions of \mathrm{\Delta }m. From Fig. 7, the even cumulants are expected to behave likely.

4 Conclusion

In this work, we developed a theoretical formalism to investigate phononic transport in disordered systems. Our formalism is based on full counting statistics theory combined with the coherent potential approximation. This formalism is capable of calculating disorder averaging of the 2n-th Green’s function, as well as disorder averaging of Green’s functions with different energies. We applied our theory to study the full counting statistics of phonon current in the presence of binary or Anderson disorders. The benchmark of our method shows 20 times speedup ratio compared with the brute force method in calculating phonon transmission moments, and the numerical accuracy is kept. Specifically, we have calculated phonon current cumulants up to the 8-th order. We have also derived a general relation among different coefficients of phonon current cumulants from the symmetry of the generating function, which fits well with numerical results.

References

[1] Z. Liu, X. Yu, and K. Chen, Thermal transport associated with ballistic phonons in asymmetric quantum structures, Front. Phys. China 4(3), 420 (2009)

[2] R. Su, Z. Yuan, J. Wang, and Z. G. Zheng, Interface-facilitated energy transport in coupled Frenkel–Kontorova chains, Front. Phys. 11(2), 114401 (2016)

[3] T. Yamamoto, S. Watanabe, and K. Watanabe, Universal features of quantized thermal conductance of carbon nanotubes, Phys. Rev. Lett. 92(7), 075502 (2004)

[4] C. Chang, D. Okawa, H. Garcia, A. Majumdar, and A. Zettl, Breakdown of Fourier’s law in nanotube thermal conductors, Phys. Rev. Lett. 101(7), 075903 (2008)

[5] X. Chen, Y. Xu, J. Wang, and H. Guo, Valley filtering effect of phonons in graphene with a grain boundary, Phys. Rev. B 99(6), 064302 (2019)

[6] P. Wang, L. Lu, and K. Bertoldi, Topological phononic crystals with one-way elastic edge waves, Phys. Rev. Lett. 115(10), 104302 (2015)

[7] S. H. Mousavi, A. B. Khanikaev, and Z. Wang, Topologically protected elastic waves in phononic metamaterials, Nat. Commun. 6(1), 8682 (2015)

[8] R. Süsstrunk and S. D. Huber, Observation of phononic helical edge states in a mechanical topological insulator, Science 349(6243), 47 (2015)

[9] Z. Y. Ong and C. H. Lee, Transport and localization in a topological phononic lattice with correlated disorder, Phys. Rev. B 94(13), 134203 (2016)

[10] H. B. Zhou, G. Zhang, J. S. Wang, and Y. W. Zhang, Phonon transport in a one-dimensional harmonic chain with long-range interaction and mass disorder, Phys. Rev. E 94(5), 052123 (2016)

[11] Y. Q. Ke, K. Xia, and H. Guo, Disorder scattering in magnetic tunnel junctions: Theory of nonequilibrium vertex correction, Phys. Rev. Lett. 100(16), 166805 (2008)

[12] Y. Q. Ke, K. Xia, and H. Guo, Oxygen-vacancy-induced diffusive scattering in Fe/MgO/Fe magnetic tunnel junctions, Phys. Rev. Lett. 105(23), 236801 (2010)

[13] P. Soven, Coherent-potential model of substitutional disordered alloys, Phys. Rev. 156(3), 809 (1967)

[14] D. W. Taylor, Vibrational properties of imperfect crystals with large defect concentrations, Phys. Rev. 156(3), 1017 (1967)

[15] X. X. Ni, M. L. Leek, J. S. Wang, Y. P. Feng, and B. W. Li, Anomalous thermal transport in disordered harmonic chains and carbon nanotubes, Phys. Rev. B 83(4), 045408 (2011)

[16] J. X. Zhai, Q. Y. Zhang, Z. H. Cheng, J. Ren, Y. Q. Ke, and B. W. Li, Anomalous transparency induced by cooperative disorders in phonon transport, Phys. Rev. B 99(19), 195429 (2019)

[17] Z. H. Cheng, J. X. Zhai, Q. Y. Zhang, and Y. Q. Ke, Auxiliary coherent medium theory for lattice vibrations in random binary alloys with mass and force-constant disorders, Phys. Rev. B 99(13), 134202 (2019)

[18] W. R. Mondal, T. Berlijn, M. Jarrell, and N. S. Vidhyadhiraja, Phonon localization in binary alloys with diagonal and off-diagonal disorder: A cluster Green’s function approach, Phys. Rev. B 99(13), 134203 (2019)

[19] B. Fu, L. Zhang, Y. Wei, and J. Wang, Full counting statistics of conductance for disordered systems, Phys. Rev. B 96(11), 115410 (2017)

[20] J. A. Young and J. U. Koppel, Phonon spectrum of graphite, J. Chem. Phys. 42(1), 357 (1965)

[21] R. A. Jishi, L. Venkataraman, M. S. Dresselhaus, and G. Dresselhaus, Phonon modes in carbon nanotubules, Chem. Phys. Lett. 209(1–2), 77 (1993)

[22] L. G. C. Rego and G. Kirczenow, Quantized thermal conductance of dielectric quantum wires, Phys. Rev. Lett. 81(1), 232 (1998)

[23] M. P. Blencowe, Quantum energy flow in mesoscopic dielectric structures, Phys. Rev. B 59(7), 4992 (1999)

[24] J. S. Wang, J. Wang, and J. T. Lü, Quantum thermal transport in nanostructures, Eur. Phys. J. B 62(4), 381 (2008)

[25] X. Chen, Y. Liu, and W. Duan, Thermal engineering in low-dimensional quantum devices: A tutorial review of nonequilibrium Green’s function methods, Small Methods 2(6), 1700343 (2018)

[26] M. P. Lopez Sancho, J. M. Lopez Sancho, J. M. L. Sancho, and J Rubio, Highly convergent schemes for the calculation of bulk and surface Green functions, J. Phys. F: Met. Phys. 15, 851 (1985)

[27] Z. Yu, G. M. Tang, and J. Wang, Full-counting statistics of transient energy current in mesoscopic systems, Phys. Rev. B 93(19), 195419 (2016)

[28] G. M. Tang, Z. Yu, and J. Wang, Full-counting statistics of energy transport of molecular junctions in the polaronic regime, New J. Phys. 19(8), 083007 (2017)

[29] G. M. Tang and J. Wang, Full-counting statistics of charge and spin transport in the transient regime: A nonequi- librium Green’s function approach, Phys. Rev. B 90(19), 195422 (2014)

[30] G. M. Tang, Y. Xing, and J. Wang, Short-time dynamics of molecular junctions after projective measurement, Phys. Rev. B 96(7), 075417 (2017)

[31] J. S. Wang, B. K. Agarwalla, H. Li, and J. Thingna, Nonequilibrium Green’s function method for quantum thermal transport, Front. Phys. 9(6), 673 (2014)

[32] L. Zhang, B. Fu, B. Wang, Y. Wei, and J. Wang, Frequency-dependent transport properties in disordered systems: A generalized coherent potential approximation approach, Phys. Rev. B 99(15), 155406 (2019)

Rights & permissions

Higher Education Press 2021