{'fraction_non_alphanumeric': 0.03712499145766418, 'fraction_numerical': 0.012693911023030137, 'mean_word_length': 4.206635829923501, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heat transport and phonon localization in mass-disordered harmonic crystals

We investigate the steady state heat current in two and three dimensional disordered harmonic crystals in a slab geometry, connected at the boundaries to stochastic white noise heat baths at different temperatures.The disorder causes short wavelength phonon modes to be localized so the heat current in this system is carried by the extended phonon modes which can be either diffusive or ballistic. Using ideas both from localization theory and from kinetic theory we estimate the contribution of various modes to the heat current and from this we obtain the asymptotic system size dependence of the current. These estimates are compared with results obtained from a numerical evaluation of an exact formula for the current, given in terms of a frequency transmission function, as well as from direct nonequilibrium simulations. These yield a strong dependence of the heat flux on boundary conditions. Our analytical arguments show that for realistic boundary conditions the conductivity is finite in three dimensions but we are not able to verify this numerically, except in the case where the system is subjected to an external pinning potential. This case is closely related to the problem of localization of electrons in a random potential and here we numerically verify that the pinned three dimensional system satisfies Fourier's law while the two dimensional system is a heat insulator. We also investigate the inverse participation ratio of different normal modes.

We investigate the steady state heat current in two and three dimensional disordered harmonic crystals in a slab geometry, connected at the boundaries to stochastic white noise heat baths at different temperatures. The disorder causes short wavelength phonon modes to be localized so the heat current in this system is carried by the extended phonon modes which can be either diffusive or ballistic. Using ideas both from localization theory and from kinetic theory we estimate the contribution of various modes to the heat current and from this we obtain the asymptotic system size dependence of the current. These estimates are compared with results obtained from a numerical evaluation of an exact formula for the current, given in terms of a frequency transmission function, as well as from direct nonequilibrium simulations. These yield a strong dependence of the heat flux on boundary conditions. Our analytical arguments show that for realistic boundary conditions the conductivity is finite in three dimensions but we are not able to verify this numerically, except in the case where the system is subjected to an external pinning potential. This case is closely related to the problem of localization of electrons in a random potential and here we numerically verify that the pinned three dimensional system satisfies Fourier's law while the two dimensional system is a heat insulator. We also investigate the inverse participation ratio of different normal modes.

I. INTRODUCTION
Energy transport in dielectric crystals at low temperatures is limited by isotope mass disorder and by anharmonicities. For a crystal coupled to thermal heat baths, in a slab geometry, the total energy transported may in addition depend on how the bath is coupled to the boundaries of the bulk. In this paper we will ignore anharmonicities but study the two other effects in considerable depth.
Based on empirical evidence, one generically expects the validity of Fourier's law, i.e. , for a slab in contact with heat reservoirs at different temperatures the average energy current, J, should be proportional to 1/N , with N being the slab length. There have been many attempts to derive Fourier's law from microscopic dynamics. Using very reasonable physical assumptions of local equilibrium and the notion of the mean free path traveled by phonons between collisions yields a heuristic derivation of Fourier's law. There is however no fully convincing derivation. There have also been many computer simulations of the heat flux in the nonequilibrium stationary states of anharmonic crystals kept in contact with thermal reservoirs at different temperatures, and theoretical analysis based on the Green-Kubo formalism [1,2,3]. These studies suggest (but there is no proof despite some claims) that Fourier' law is not valid for one and two dimensional systems, even in the presence of anharmonic interactions, unless the system is also subjected to an external substrate pinning potential. Generically it is found that, for a system in contact with heat reservoirs at different fixed temperatures, the heat current density J scales anomalously with system length N as with µ = 1. The effective thermal conductivity behaves then as κ ∼ N α where α = 1 − µ. For two dimensional systems there are some analytic studies which suggest κ ∼ ln(N ).
Recent experiments on heat conduction in nanotubes and graphene flakes have reported observations which indicate such divergence of κ with system size [4,5]. For heat conduction in the ordered harmonic crystal there are exact results from which one has µ = 0 in all dimensions [6,7]. Heat conduction in a disordered harmonic crystal will be affected by Anderson localization [8] and by phonon scattering.
In this paper we report results of heat conduction studies in 2D and 3D disordered harmonic lattices with scalar displacements, connected to heat baths modeled by Langevin equations with white noise. We pay particular attention to the interplay between localization effects, boundary effects, and the role of long wavelength modes. The steady state heat current is given exactly as an integral over all frequencies of a phonon transmission coefficient. Using this formula and heuristic arguments, based on localization theory and kinetic theory results, we estimate the system size dependence of the current. The main idea behind our arguments is that the phonon states can be classified as ballistic modes, diffusive modes and localized modes. The classification refers both to the character of the eigenfunctions as well as to their transmission properties. Ballistic modes are spatially extended and approximately periodic; their transmission is independent of system size. Diffusive modes are extended but non-periodic and their transmission decays as 1/N . For localized modes transmission decays exponentially with N . In the context of kinetic theory calculations, the ballistic modes are the low frequency modes with phonon mean free path K (ω) > ∼ N , and their contribution to the current leads to divergence of the thermal conductivity. Here we will carefully examine the effect of boundary conditions on these modes.
Numerically we use two different approaches to study the nonequilibrium stationary state. The first is a numerical one which relies on the result that the current can be expressed in terms of a transmission coefficient. This transmission co-efficient can be written in terms of phonon Green's functions and we implement efficient numerical schemes to evaluate this. The second approach is through direct nonequilibrium simulations of the Langevin equations of motion and finding the steady state current and temperature profiles. We have also studied properties of the isolated system, i.e., of the disordered lattice without coupling to heat baths and looked at the normal mode frequency spectrum and the wavefunctions. One measure of the degree of localization of the normal modes of the isolated system is the so-called inverse participation ratio [IPR, defined in Eq. (10) below]. We have carried out studies of the IPR and linked these with the results from the transmission study.
Phonon localization: This is closely related to the electron localization problem. The effect of localization on linear waves in disordered media has been most extensively studied in the context of the Schrödinger equation for non-interacting electrons moving in a disordered potential. Looking at the eigenstates and eigenfunctions of the isolated system of a single electron in a disordered potential one finds that, in contrast to the spatially extended Bloch states in periodic potentials, there are now many eigenfunctions which are exponentially localized in space. It was argued by Mott and Twose [9] and by Borland [10], and proven rigorously by Goldsheid et al. [11], that in one dimension (1D) all states are exponentially localized. In two dimensions (2D) there is no proof but it is believed that again all states are localized. In three dimensions (3D) there is expected to be a transition from extended to localized states as the energy is moved towards the band edges [12]. The transition from extended to localized states, which occurs when the disorder is increased, changes the system from a conductor to an insulator. The connection between localization and heat transport in a crystal is complicated by the fact that phonons of all frequencies can contribute to energy transmission across the system. In particular account has to be taken of the fact that low frequency phonon modes are only weakly affected by disorder and always remain extended. The heat current carried by a mode which is localized on a length scale , decays with system length N as e −N/ . This depends on the phonon frequency and low frequency modes for which ∼ N will therefore be carriers of the heat current. The net current then depends on the nature of these low frequency modes and their scattering due to boundary conditions (BCs).
A renormalization group study of phonon localization in a continuum vector displacement model was carried out by John et al. [13]. They found that much of the predictions of the scaling theory of localization for electrons carry over to the phonon case. Specifically they showed that in one and two dimensions all non-zero frequency phonons are localized with the low frequency localization length diverging as ∼ ω −2 and ∼ e c/ω 2 , respectively (where c > 0 is some constant). This means that in 1D all modes with ω > ∼ ω L c = N −1/2 are localized while in 2D all modes with ω > ∼ ω L c = [log(N )] −1/2 are localized. In 3D the prediction is that there is an ω L c independent of N above which all modes are localized. However this study does not make any statements on the system size dependence of the conductivity.
Kinetic theory: If one considers the low frequency extended phonons, then the effect of disorder is weak and in dimensions d > 1 one expects that localization effects can be neglected and kinetic theory should be able to provide an accurate description. In this case one can think of Rayleigh scattering of phonons. This gives an effective mean free path where v, the sound velocity, can be taken to be a constant. For a finite system of linear dimension N we have D(ω) = vN for ω Kinetic theory then predicts where ρ(ω) ∼ ω d−1 is the density of states and we get κ ∼ N 1/(d+1) implying µ = d/(d + 1). The divergence of the phonon mean free path at low frequencies and the resulting divergence of the thermal conductivity of a disordered harmonic crystal has been discussed in the literature and it has been argued that anharmonicity is necessary to make κ finite [14,15].
Simulation results: There have been only few simulation studies of heat conduction in three dimensional disordered systems and none have been definitive concerning the validity of Fourier's law [16,17]. In two dimension a diverging thermal conductivity was reported in [18]. Some other studies have also looked at heat conduction in glassy systems at low temperatures where the harmonic approximation was used [19,20] but these did not address the questions of Ndependence of κ and the validity of Fourier's law.
In 1D it is well known from rigorous results and numerical studies that α = 0 and its value is strongly dependent on boundary conditions [21,22,23,24,25,26]. For fixed BCs one has α = −1/2 while for free BCs, α = 1/2. The precise definitions of the different BCs will be given later.
Here we explain the different physical situations they correspond to. If we model the heat reservoirs themselves by infinite ordered harmonic crystals then Langevin type equations for the system [27] are obtained on eliminating the bath degrees of freedom. The two different BCs then emerge naturally. Fixed BCs correspond to reservoirs with properties different from the system (e.g. different spring constants) and in this case one finds that effectively the particles at the boundaries (those coupled to reservoirs) experience an additional harmonic pinning potential. Free BCs correspond to the case where the reservoir is simply an extension of the system (without disorder) and in this case the end particles are unpinned. Free BCs have been studied in the literature in the context of heat conduction in one dimensional chains [22] and in studies on nanotubes [28,29]. In this paper we study lattices with both fixed and free BCs although we think that fixed BCs are more realistic.
In the presence of an external pinning potential low frequency modes are suppressed, hence one expects qualitative differences in transport properties. The pinned system has often been used as a model system to study the validity of Fourier's law. It has no translational invaraince and is thus more closely related to the problem of electrons moving in a random potential. Here we consider systems with and without external pinning potentials.
The rest of the paper is organized as follows. In Sec. (II) we define the specific model studied by us and present some general results for heat conduction in harmonic Hamiltonian systems connected to Langevin baths. We also give some details of the numerical and simulation methods used in the paper. The transfer matrix approach used in evaluating the phonon transmission function is explained in Appendix B. In Sec. (III) a brief review of results for the one dimensional case and the heuristic arguments for the higher dimensional cases are given. In Sec. (IV) we present results from both the numerical approach and from nonequilibrium simulations. The main results presented are for transmission functions, IPRs of normal modes and the system size dependence of the current in two and three dimensional disordered harmonic lattices. Along the way we also present results for the density of states ρ(ω). Finally we conclude with a discussion in Sec. (V).

II. MODELS AND METHODS
For simplicity we consider only the case where longitudinal and transverse vibration modes are decoupled and hence we can describe the displacement at each site by a scalar variable. Also we restrict our study to d-dimensional hypercubic lattices. Let us denote the lattice points by the vector n = {n 1 , n 2 , ..., n d } with n ν = 1, 2, ..., N . The displacement of a particle at the lattice site n is given by x n . In the harmonic approximation the system Hamiltonian is given by whereê refers to the 2d nearest neighbors of any site and we impose boundary conditions which will be specified later.
We have also included an external pinning harmonic potential with spring constant k o , which we will sometimes set equal to zero. We consider a binary mass disordered crystal. Specifically we set the masses of exactly half the particles at randomly chosen sites to bem − ∆ and the rest to bem + ∆. Thus ∆ gives a measure of the disorder. We couple all the particles at n 1 = 1 and n 1 = N to heat reservoirs, at temperatures T L and T R respectively, and use periodic boundary conditions in the other (d − 1) directions. The heat conduction takes place along the 1 direction. Each layer with constant n 1 consists of N = N d−1 particles. The heat baths are modeled by white noise Langevin equations of motion for the particles coupled to the baths. Using the notation n = (n 1 , n ), the equations of motion are given by: where the dissipative and noise terms are related by the usual fluctuation dissipation relations The particles at the surfaces n 1 = 1, N experience additional harmonic pinning potentials with spring constants k arising from coupling to the heat reservoirs. We consider two kinds of boundary conditions at the surfaces connected to reservoirs: (i) fixed BCs k > 0 and (ii) free BCs k = 0. As discussed in the introduction fixed BCs correspond to reservoirs with properties different from the system while free reservoirs correspond to the case where the reservoir is really an extension of the system but without disorder. For the pinned case we only consider fixed BC. A schematic of the models and the different boundary conditions that we have studied is given in Fig. (1). Henceforth we will use dimensionless variables: forceconstants are measured in units of k, masses in units of the average massm, time in units of the inverse frequency Ω −1 = (m/k) 1/2 , displacements are in units of the lattice spacing a, friction constant γ is in units ofmΩ, and finally temperature is measured in units ofma 2 Ω 2 /k B .
Driven by the reservoirs at two different temperatures T L and T R the system reaches a nonequilibrium steady state.
Our main interest will be in the steady state heat current in the system. Given the Langevin equations of motion Eq (4), one can find a formal general expression for the current. Let us denote by X a column vector with N d elements consisting of the displacements at all lattice sites. Similarly leṫ X represent the vector for velocities at all sites. Then we can write the Hamiltonian in Eq. (3) in the compact form H = 1 2Ẋ T MẊ + 1 2 X T VX , which defines the diagonal mass matrix M and the force constant matrix V. With this notation we have the following form for the steady state current per bond from the left to the right reservoir [23,27] : where and ∆T = T L − T R . The S + L , S + R represent terms arising from the coupling to the left and right baths respectively, and I L,R = Im[S + L,R ]. The specific form of S + L,R for our system described by Eqs. (4) is given in Appendix B. The matrix G + (ω) can be identified as the phonon Green's function of the system with self-energy corrections due to the baths [27]. The integrand in Eq. (6) T N (ω) can be thought of as the transmission coefficient of phonons at frequency ω from the left to the right reservoir. It will vanish, when N → ∞ , at values of ω for which the disorder averaged density of states is zero. Note that due to the harmonic nature of the forces the dependence of the heat flux on the reservoir temperatures enters only through the term ∆T in Eq. (6). The above expression for the current is of the Landauer form and has been derived using various other approaches such as scattering theory [30,31] and the nonequilibrium Green's function formalism [32,33].
Numerical approach: In Appendix B we describe how T N can be expressed in a form amenable to accurate numerical computation. The system sizes we study are sufficiently large so that T N (ω) has appreciable values only within the range of frequencies of normal modes of the isolated system, i.e., corresponding to γ = 0 in Eq. (4). Outside this range we find that the transmission rapidly goes to zero. By performing a discrete sum over the transmitting range of frequencies we do the integration in Eq. (6) to obtain the heat current density J . In evaluating the discrete sum over ω, step sizes of δω = 0.01−0.0001 are used and we verified convergence in most cases. With our choice of units we have k = 1,m = 1 and we fixed ∆T = 1. Different values of the mass variance ∆ and the on-site spring constant k o were studied for two and three dimensional lattices of different sizes. It is expected that the value of the exponent µ will not depend on γ and in our calculations we mostly set γ = 1, except when otherwise specified.
Simulation approach: The simulations of Eq. (4) are performed using a velocity-Verlet scheme as given in [34]. The current and temperature profiles in the system are obtained from the following time averages in the nonequilibrium steady state: n = 2, 3, ..., N , , n = 1, 2, ..., N .
We then obtained the average current J = ( In the steady state one has J n = J for all n and stationarity can be tested by checking how accurately this is satisfied. We chose a step size of ∆t = 0.005 and equilibrated the system for over The value of T N (ω), J and T n depend, of course, on the particular disorder realization. Mostly we will here be interested in disorder averages of these quantities which we will denote by [T ], J = [J ] and [T n ]. We also define the disorder averaged transmission per bond with the notation Numerical analysis of eigenmodes and eigenfunctions: We have studied the properties of the normal modes of the disordered harmonic lattices in the absence of coupling to reservoirs, again with both free and fixed boundary conditions. The d-dimensional lattice has p = 1, 2, ..., N d normal modes and we denote the displacement field corresponding to the p th mode by a n (p) and the corresponding eigenvalue by ω 2 p . The normal mode equation corresponding to the Hamiltonian in Eq. (3) is given by: where the a n satisfy appropriate boundary conditions. Introducing variables ψ n (p) = m 1/2 n a n (p), v n = (2d + k o )/m n and t n,l = 1/(m n m l ) 1/2 for nearest neighbour sites n, l the above equation transforms to the following form: This has the usual structure of an eigenvalue equation for a single electron moving in a d-dimensional lattice corresponding to a tight-binding Hamiltonian with nearest neighbour hopping t n,l and on-site energies v n . Note that t n,l and v n are correlated random variables, hence the disorder-energy diagram might differ considerably from a single band Anderson tight-binding model. We have numerically evaluated all eigenvalues and eigenstates of the above equation for finite cubic lattices of size upto N = 64 in 2D and N = 16 in 3D. One measure of the degree of localization of a given mode is the inverse participation ratio (IPR) defined as follows: For a completely localized state, i.e. a n = δ n,n0 , P −1 takes the value 1. On the other hand for a completely delocalized state, for which a n = N −d/2 e in.q where q is a wave vector, P −1 takes the value N −d . We will present numerical results for the IPR calculated for all eigenstates of given disorder realizations, in both 2D and 3D. Finally we will show some results for the density of states, ρ(ω), of the disordered system defined by: The density of states of disordered binary mass harmonic crystals was studied numerically by Payton and Visscher in 1967 [35] and reviewed by Dean in 1972 [36].

III. HEAT CONDUCTION IN DISORDERED HARMONIC CRYSTALS: GENERAL CONSIDERATIONS
Let us first briefly consider heat conduction in the one dimensional disordered harmonic chain. This has been extensively studied and is well understood [21,22,23,24,25,26]. The matrix formulation explained in the last section leads to a clear analytic understanding of the main results. The current is given by the general expression Eq. (6). From Eq. (B5) the transmission is given by There are three observations that enable one to determine the asymptotic system size dependence of the current. These are: is a complex number which can be expressed in terms of the product of N random 2 × 2 matrices. Using Furstenberg's theorem it can be shown that for almost all disorder realization, the large N behaviour of P (1,N ) for fixed ω > 0 is |P (1,N ) | ∼ e bN ω 2 , where b > 0 is a constant. This is to be understood in the sense that lim N →∞ (1/N ) log |P 1,N ) | ∼ bω 2 for ω → 0. Since T N (ω) ∼ |P (1,N ) | −2 ∼ e −2bN ω 2 , this implies that transmission is significant only for low frequencies ω < ∼ ω c (N ) ∼ 1/N 1/2 . The current is therefore dominated by the small ω behaviour of T N (ω).
(ii) The second observation made in [25] is that the transmission for ω < ω c (N ) is ballistic in the sense that T N (ω) is insensitive to the disorder.
(iii) The final important observation is that the form of the prefactors of e −bN ω 2 in T N (ω) for ω < ω c (N ) depends strongly on boundary conditions and bath properties [25,26]. For the white noise Langevin baths one finds T N (ω) ∼ ω 2 e −bN ω 2 for fixed BC and T N (ω) ∼ ω 0 e −bN ω 2 for free BC [26]. This difference arises because of the scattering of long wavelength modes by the boundary pinning potentials.
In Fig. (2) we plot numerical results showing T N (ω) for the 1D binary mass-disordered lattice with both fixed and free boundary conditions. One can clearly see the two features discussed above namely (i) dependence of frequency cut-off on system size and (ii) dependence of form of T N (ω) on boundary conditions. Using the three observations made above it is easy to arrive at the conclusion that J ∼ N −3/2 for fixed BC and J ∼ N −1/2 for free BC. In the presence of a pinning potential the low-frequency modes are suppressed and one obtains a heat insulator with J ∼ e −cN , with c a constant [24] (see also [37] and references there).
Higher dimensions. Let us try to extend the analysis of the 1D case to higher dimensions. For this we will use inputs from both kinetic theory and the theory of phonon localization. The main point of our arguments involves the assumption that normal modes can be classified as ballistic, diffusive or localized. Using localization theory we determine the frequency region where states are localized. The lowest frequency states with ω → 0 will be ballistic and we use kinetic theory to determine the fraction of extended states which are ballistic. We assume that at sufficiently low frequencies the effective disorder is always weak (even when the mass variance ∆ is large) and one can still use kinetic theory. Corresponding to the three observations made above for the 1D case we now make the following arguments: (i) From localization theory one expects all fixed non-zero frequency states in a 2D disordered system to be localized when the size of the system goes to infinity. As discussed in Sec. (I) localization theory gives us a frequency cut-off ω L c = (ln N ) −1/2 in 2D above which states are localized. In 3D one obtains a finite frequency cut-off ω L c independent of system size above which states are localized.
(ii) For the unpinned case with finite N there will exist low frequency states below ω L c , in both 2D and 3D, which are extended states. These states are either diffusive or ballistic. Ballistic modes are insensitive to the disorder and their transmission coefficient are almost the same as for the ordered case. To find the frequency cut-off below which states are ballistic we use kinetic theory results (see Appendix A). For the low-frequency extended states we expect kinetic theory to be reliable and this gives us a mean free path for phonons K ∼ ω −(d+1) . This means that for low frequencies ω < ∼ ω K c = N −1/(d+1) we have K (ω) > N and phonons transmit ballistically. We now proceed to calculate the contribution of these ballistic modes to the total current. This can be obtained by looking at the small ω form of T N (ω) for the ordered lattice.
(iii) For the ordered lattice T N (ω) is typically a highly oscillatory function with the oscillations increasing with system size. An effective transmission coefficient in the N → ∞ limit can be obtained by considering the integrated transmission. This asymptotic effective low-frequency form of T N (ω), for the ordered lattice can be calculated using meth-ods described in [26] and is given by: the result being valid for d = 1, 2, 3 [38].
Using the above arguments we then get the ballistic contribution to the total current density (for the unpinned case) as: We can now make predictions for the asymptotic system size dependence of total current density in two and three dimensions.
Two dimensions: From localization theory one expects that all finite frequency modes ω > ∼ ω L c = (ln N ) −1/2 are localized and their contribution to the total current falls exponentially with system size. Our kinetic theory arguments show that the low frequency extended states with The diffusive contribution to total current will then scale as J diff ∼ (ln N ) −1/2 N −1 . The ballistic contribution depends on BCs and is given by Eq. (13). This gives J ball ∼ N −4/3 for fixed BC and J ball ∼ N −2/3 for free BC. Hence, adding all the different contributions, we conclude that asymptotically: In the presence of an onsite pinning potential at all sites the low frequency modes get cut off and all the remaining states are localized, hence we expect: where b is some positive constant. Three dimensions: In this case localization theory tells us that modes with ω > ∼ ω L c are localized and ω L c is independent of N . From kinetic theory we find that the extended c are ballistic. The contribution to current from diffusive modes scales as J diff ∼ N −1 . The ballistic contribution (from states with ω < ∼ N −1/4 ) is obtained from Eq. (13) and gives J ball ∼ N −5/4 for fixed BC and J ball ∼ N −3/4 for free BC. Hence, adding all contributions, we conclude that asymptotically: In the presence of an onsite pinning potential at all sites the low frequency modes get cut off and, since in this case the remaining states form bands of diffusive and localized states, hence we expect: Thus in 3D both the unpinned lattice with fixed boundary conditions and the pinned lattice are expected to show Fourier type of behaviour as far as the system size dependence of the current is considered. Note that for free BC, the prediction for the current contribution from the ballistic part J ball ∼ N −d/(d+1) is identical to that from kinetic theory discussed in Sec. (I). This agreement can be traced to the small ω form of T (ω) ∼ ω d−1 for free BC [see Eq. (12)] which is identical to the form of the density of states ρ(ω) used in kinetic theory. The typical form of density of states for ordered and disordered lattices in different dimensions is shown in Fig. (3) and we can see that the low frequency form is similar in both cases and has the expected ω d−1 behaviour. However it seems reasonable to expect that, since the transport current phonons are injected at the boundaries, in kinetic theory one needs to use the local density of states evaluated at the boundaries. For fixed BC this will then give rise to an extra factor of ω 2 (from the squared wavefunction) and then the kinetic theory prediction matches with those given above.
We note that the density of states in Fig. (3) show apparent gaps in the middle ranges of ω for d = 2, 3. These might be expected to disappear when the size of the system goes to infinity when there should be large regions containing only masses of one type [23,24]. These regions will however be rare. In Fig. (4) we show plots of the density of states for the ordered and disordered harmonic lattices in the presence of pinning. In this case the gaps in the spectrum are more pronounced and, for large enough values of k o and ∆, may be present even in the thermodynamic limit.

IV. RESULTS FROM NUMERICS AND SIMULATIONS
We now present the numerical and simulation results for transmission coefficients, heat current density, temperature profiles and IPRs for the disordered harmonic lattice in various dimensions. The numerical scheme for calculating J is both faster and more accurate than nonequilibrium simulations. Especially, for strong disorder, equilibration times in nonequilibrium simulations become very large and in such cases only the numerical method can be used. However we also show some nonequilibrium simulation results. Their almost perfect agreement with the numerical results provides additional confidence in the accuracy of our results. In Sec. (IV A) we give the results for the 2D lattice for the unpinned case with both fixed and free boundary conditions and then for the pinned case. In Sec. (IV B) we present the results for the three dimensional case with and without substrate pinning potentials.

A. Results in two dimensions
In this section we consider N × N square lattices with periodic BCs in the ν = 2 direction and either fixed or free BCs in the conducting direction (ν = 1). One of the interesting questions here is as to how the three properties for the 1D case discussed in Sec. (III) get modified for the 2D case.

Disordered 2D lattice without pinning
Fixed BC : we have computed the transmission coefficients and the corresponding heat currents for different values of ∆ and for system sizes from N = 16 − 1024. The number of averages varied from over 100 samples for N = 16 to about two samples for N = 1024. In Figs. (5,6,7) we plot the disorder averaged transmission coefficient for three different disorder strengths, ∆ = 0.95, ∆ = 0.8 and ∆ = 0.2, for different system sizes. The corresponding plots of IPRs as a function of normal mode frequency ω p , for single disorder realizations, are also given. From the IPR plots we get an idea of the typical range of allowed normal mode frequencies and their degree of localization. Low IPR values which scale as N −2 imply extended states while large IPR values which do not change much with system size denote localized states. In Fig. (6) we also show typical plots of small IPR and large IPR wavefunctions. From Figs. (5,6,7) we make the following observations: (i) As expected we see significant transmission only over the range of frequencies with extended states. Thus in Fig. (5) for ∆ = 0.95 we see that, while there are normal modes in the range ω ≈ (0 − 12), transmission is appreciable only in the range ω ≈ (0 − 1.5) and this is also roughly the range where the IPR data shows a N −2 scaling behaviour. This can also be seen in Fig. (6) where the inset shows the decay of T (ω) in the localized region. Unlike the 1D case we see a very weak dependence on system size of the upper frequency cut-off ω L c beyond which states are localized and transmission is negligible. As discussed earlier, localization theory predicts ω L c ∼ (ln N ) −1/2 but this may be difficult to observe numerically. The overall transmission function T N (ω) decreases with increasing system size, with T (ω) ∼ 1/N at higher freqencies and T (ω) ∼ N 0 at the lowest frequencies.
(ii) In Fig. (7) we have also plotted T (ω) for the ordered binary mass case and we note that over a range of small frequencies, T (ω) for the disordered case is very close to the curve for the ordered case, which means that these modes are ballistic. As expected from the arguments in Sec. (III) we roughly find T (ω) ∼ ω 3 at small frequencies. The remaining transmitting states are either diffusive (with a 1/N scaling) or are in the cross-over regime between diffusive and ballistic and so do not have a simple scaling.
We next look at the the integrated transmission which gives the net heat current. The system size dependence of the disorder averaged current J for different values of ∆ is shown in Fig. (8). For the case ∆ = 0.2, we also show simulation results and one can see that there is excellent agreement with the numerical results. For ∆ = 0.2 we get an exponent µ ≈ 0.6 which is close to the value obtained earlier in [18] for a similar disorder strength. However with increasing disorder we see that this value changes and seems to settle to around µ ≈ 0.75. It seems reasonable to expect (though we have no rigorous arguments) that there is only one asymptotic exponent and for small disorder one just needs to go to very large system sizes to see the true value. In Fig. (9) we show temperature profiles obtained from simulations for lattices of different sizes with ∆ = 0.2. The jumps at the boundaries indicate that the asymptotic size limit has not yet been reached. This is consistent with our result that the exponent µ obtained at ∆ = 0.2 is different from what we believe is the correct asymptotic value (obtained at larger values of ∆). We do not have temperature plots at strong disorder where simulations are difficult.
Thus contrary to the arguments in Sec. (III) which predicted J ∼ (ln N ) −1/2 N −1 we find a much larger current scaling as J ∼ N −0.75 . It is possible that one needs to go to larger system sizes to see the correct scaling.
Free BC : In this case from the arguments in Sec. (III) we expect ballistic states to contribute most significantly to the current density giving J ∼ N −2/3 .
In Figs. (10,11) we plot the disorder averaged transmission coefficient for ∆ = 0.8 and ∆ = 0.2 for different system sizes. Qualitatively these results look very similar to those for fixed boundaries. However transmission is now significantly larger in the region of extended states. The behaviour at frequencies ω → 0 is also different and we now find T (ω) ∼ ω in contrast to T (ω) ∼ ω 3 for fixed boundaries. From the plots of IPRs in Fig. (10) we note that there is not much qualitative difference with the fixed boundary plots except in the low frequency region (see below).
The system size dependence of the disorder averaged current J for two different values of ∆ is shown in Fig. (12). For ∆ = 0.2 we get an exponent µ ≈ 0.5 while for the stronger disorder case ∆ = 0.8 we see a different exponent µ ≈ 0.6. Again we believe that the strong disorder value of µ = 0.6 is closer to the value of the true asymptotic exponent. This value is close to the expected µ = 2/3 for free BC and significantly different from the value obtained for fixed BC (µ ≈ 0.75). Thus the dependence of the value of σ on boundary conditions exists even in the 2D case.
For the case of free BCs, we find that the values of T (ω) in the diffusive regime matches with those for fixed BCs but are completely different in the ballistic regime. This is seen in Fig. (13) where we plot the effective mean free path l eff (ω) = N T (ω)/w d−1 in the low-frequency region [this is obtained by comparing Eq. (6) with the kinetic theory expression for conductivity Eq. (2)]. For free BC, l eff is roughly consistent with the kinetic theory prediction l −1 eff ∼ N −1 + K −1 (ω) but the behaviour for fixed BC is very different. The inset of Fig. (13) plots l eff for the equal mass ordered case and we find that in the ballistic regime it is very close to the disordered case, an input that we used in the heuristic derivation. The numerical data also confirms that for small ω, T (ω) ∼ ω for free BCs and as ω 3 for fixed BCs. The transmission for fixed BC shows rapid oscillations which increase with system size, and arise from scattering and interference of waves at the interfaces.

Disordered 2D lattice with pinning
We now study the effect of introducing a harmonic pinning potential at all sites of the lattice. It is expected that this will cut off low frequency modes and hence one should see strong localization effects. The localization length will decrease both with increasing ∆ and increasing k o (in 1D heuristic arguments give ∼ 1/(∆ 2 k o ) [37]). In Figs. (14,15) we plot the transmission coefficients for two cases with on-site potentials k o = 10.0 and k o = 2.0 respectively, and ∆ = 0.4. We also plot the IPR in Fig. (14). Unlike in the unpinned case we now find that the transmission coefficients are much smaller and fall more rapidly with system size.
From the plot of P −1 we find that for all the modes, the value of P −1 does not change much with system size which implies that all modes are localized. The allowed frequency bands correspond to the transmission bands. The two wavefunctions plotted in Fig. (14) correspond to one relatively small and one large P −1 value and clearly show that both states are localized.
The system size dependence of the integrated current is shown in Fig. (16) for the two parameter sets. The values of µ ≈ 1.6, 3.65 for the two sets indicate that at large enough length scales one will get a current falling exponentially with system size and hence we have an insulating phase. In Fig. (17) we plot the temperature profiles for the set with ∆ = 0.4, k o = 10.0 . In this case it is difficult to obtain steady state temperature profiles from simulations for larger system sizes. The reason is that the temperature (unlike current) gets contributions from all modes (both localized and extended) and equilibrating the localized modes takes a long time.

B. Results in three dimensions
In this section we mostly consider N × N × N lattices with periodic boundary conditions in the ν = 2, 3 directions. Some results for N × N 2 × N 3 lattices with N 2 = N 3 < N will also be described. Preliminary results for the case of free BCs are given and indicate that there is no dependence of the exponent µ on BCs. It is not clear to us whether this is related to the boundedness of the fluctuations in x n and the decay of the correlations between x n and x l (like |n − l| −1 ) in d = 3 and their growth (with N ) in d < 3.

Disordered 3D lattice without pinning
Fixed BC : we have used both the numerical approach and simulations for sizes up to 32×32×32 for which we have data for T (ω). For larger systems the matrices become too big and we have not been able to use the numerical approach. Hence, for larger system sizes we have only performed simulations, including some on N × N 2 × N 2 lattices. For these cases only the current J is obtained. The number of averages varies from over 100 samples for N = 16 to two samples for N = 64. In Figs. (18,19) we plot the disorder averaged transmission coefficient for two different disorder strengths, ∆ = 0.8 and ∆ = 0.2, for different system sizes. The corresponding plots of IPRs as a function of normal mode frequency ω p , for single disorder realizations, are also given. From the IPR plots we get an idea of the typical range of allowed normal mode frequencies and their degree of localization. Low IPR values which scale as N −3 imply extended states while large IPR values which do not change much with system size denote localized states. From Figs. (18,19) we make the following observations. (i) From the 3D data it is clear the effect of localization is weaker than in 1D and 2D. Both for ∆ = 0.2 and ∆ = 0.8 we find that there is transmission over almost the entire range of frequencies of the allowed normal modes. From the IPR plots we see that for ∆ = 0.2 most states are extended except for a small region in the high frequency band-edge. For ∆ = 0.8 the allowed modes form two bands and one finds significant transmission over almost the full range. At the band edges (except the one at ω = 0) there are again localized states. It also appears that there are some large IPR states interspersed within the high frequency band. As in the 2D case and unlike the 1D case, the frequency range over which transmission takes place does not change with system size, only the overall magnitude of transmission coefficient changes.
(ii) The plot of N T (ω) in Fig. (18) shows the nature of the extended states. The high frequency band and a portion of the lower frequency band have the scaling T (ω) ∼ N −1 and hence corresponds to diffusive states. In the lower-frequency band the fraction of diffusive states seems to be increasing with system size but it is difficult to verify the ω K c ∼ N −1/4 scaling. The ballistic nature of the low-frequency states is confirmed in Fig. (19) where we see that T (ω) for the binarymass ordered and disordered lattices match for small ω [with a T (ω) ∼ ω 4 dependence].
In Fig. (20) we show the system size dependence of the disorder averaged current density J for the two cases with weak disorder strength (∆ = 0.2) and strong disorder strength (∆ = 0.8). The results for cubic lattices of sizes up to N = 32 are from the numerical method while the results for larger sizes are from simulations. We find an exponent µ ≈ 0.6 at small disorder and µ ≈ 0.75 at large disorder strength. As in the 2D case here too we believe that at small disorder, the asymptotic system size limit will be reached at much larger system sizes and that the exponent obtained at large disorder strength is probably close to the true asymptotic value. The value (µ = 0.75) does not agree with the prediction (J ∼ N −1 ) made from the heuristic arguments in Sec. (III). A study of larger system sizes is necessary to confirm whether or not the asymptotic size limit has been reached.
The data point at N = 128 for the set with ∆ = 0.2 in Fig. (20) actually corresponds to a lattice of dimensions 128× 48 × 48 and we believe that the current value is very close to the expected fully 3D value. To see this point, we have plotted in Fig. (21) results from nonequilibrium simulations with N × N 2 × N 2 lattices with N 2 ≤ N .
Finally, in Fig. (22) we show temperature profiles (for single disorder realizations) obtained from simulations for lattices of different sizes and with ∆ = 0.2. The jumps at the boundaries again indicate that the asymptotic system size limit has not been reached even at the largest size.
Free BC : In this case from the arguments in Sec. (III) we expect ballistic states to contribute most significantly to the current density giving J ∼ N −3/4 . In Fig. (23) we plot the disorder averaged transmission coefficient for ∆ = 0.8 for different system sizes. The transmission function is very close to that for the fixed boundary case except in the frequency region corresponding to non-diffusive states. At ω → 0 we now expect, though it is hard to verify from the data, that T (ω) ∼ ω 2 in contrast to T (ω) ∼ ω 4 for fixed boundaries.
The system size dependence of the disorder averaged current J for two different values of ∆ is shown in Fig. (20). We find that the current values are quite close to the fixed BC case and the exponent obtained at the largest system size studied for this case is µ ≈ 0.71. This value is close to the expected µ = 3/4 for free BC.
We now compare the transmission coefficient for free and fixed BCs in the ballistic regime. This is plotted in Fig. (24) where we show the effective mean free path l eff (ω) = N T (ω)/w d−1 in the low-frequency region. As in the 2D case we again find that for free BCs, l eff is roughly consistent with the kinetic theory prediction l −1 eff ∼ N −1 + K −1 (ω) and the behaviour for fixed BCs is very different. The inset of Fig. (24) plots l eff for the equal mass ordered case and we find that in the ballistic regime it is very close to the disordered case. The numerical data confirms the input in our theory on the form of T (ω) for small ω, i.e. T (ω) ∼ ω 2 for free BCs and as ω 4 for fixed BCs. The transmission for fixed BC shows rapid oscillations which increase with system size, and arise from scattering and interference of waves at the interfaces. From the IPR plots we notice that the spectrum of the 3D disordered pinned chain has a similar interesting structure as in the 2D case with two bands and a gap which is seen at strong disorder. However unlike the 2D case where all states were localized, here the IPR data indicates that most states except those at the band edges are diffusive. We see localized states at the band edges and also there seem to be some localized states interspersed among the extended states within the bands. The insets in Figs. (25,26) show that there is a reasonable N −1 scaling of the transmission data in most of the transmitting region. This is clearer at the larger system sizes. Thus, unlike the unpinned case where low frequency extended states were ballistic or super-diffusive, here we find that there is no transmittance at small (ω → 0) frequencies and that all states are diffusive.

Disordered 3D lattice with pinning
From the above discussion we expect Fourier's law to be valid in the 3D pinned disordered lattice. The system size dependence of the disorder averaged current J for different disorder strengths is plotted in Fig. (27). For all the parameter sets the exponent obtained is close to µ = 1 corresponding to a finite conductivity and validity of Fourier's law. The temperature profiles plotted in Fig. (28) have small boundary temperature jumps and indicate that the asymptotic size limit has already been reached.
One might expect that at very strong disorder, all states should become localized and then one should get a heat insulator. The parameter set corresponding to Fig. (26) corresponds to strong disorder and for this we still find a signif- The numerical (and nonequilibrium simulation) results obtained in the paper are compared, in two and three dimensions, with the analytical predictions obtained from our heuristic arguments. The error bar for the numerically obtained exponent values is of the order ±0.02. This error is estimated from the errors in the last few points of the J-versus-N data. NB: The system sizes used may well be far from asymptotic.
icant fraction of extended states. Thus for the binary mass case it appears that there are always extended states. We have some results for the case with a continuous mass distribution ( masses are chosen from a uniform distribution between 1 − ∆ and 1 + ∆). In this case we find that the effect of disorder is stronger and the transmission at all frequencies is much reduced compared to the binary mass case. However we cannot see the exponential decrease in transmission with system size and so it is not clear if an insulating behaviour is obtained. Further numerical studies are necessary to understand the asymptotic behaviour.

V. DISCUSSION
We have studied heat conduction in isotopically disordered harmonic lattices with scalar displacements in two and three dimensions. The main question addressed is the system size dependence of the heat current, which is computed using Green's function based numerical methods as well as nonequilibrium simulations. We have tried to understand the size dependence by looking at the phonon transmission function T (ω) and examining the nature of the energy transport in different frequency regimes. We also described a heuristic analytical calculation based on localization theory and kinetic theory and compared their predictions with our numerical and simulation results. This comparison is summarized in Table (I). The most interesting findings of this work are: (i) For the unpinned system we find that in 2D there are a large number of localized modes for which phonon transmission is negligible. In 3D the number of localized modes is much smaller. The extended modes are either diffusive or ballistic. Our analytic arguments show that the contribution of ballistic modes to conduction is dependent on BCs and is strongly suppressed for the case of fixed BCs, the more realistic case. In 3D this leads to diffusive modes dominating for large system sizes and Fourier's law is satisfied. Thus a finite heat conductivity is obtained for the 3D disordered harmonic crystal without the need of invoking anharmonicity as is usually believed to be necessary [14,15]. This is similar to what one obtains when one adds stochasticity to the time evolution in the bulk as shown by [39]. Our numerical results verify the predictions for free BCs and we believe that much larger system sizes are necesary to verify the fixed BC results ( this is also the case in 1D [25,26]).
(ii) In two dimensions the pinned disordered lattice shows clear evidence of localization and we obtain a heat insulator with exponential decay of current with system size. (iii) Our result for the 3D pinned disordered lattice provides the first microscopic verification of Fourier's law in a three dimensional system. For the binary mass distribution we do not see a transition to insulating behaviour with increasing disorder. For a continuous mass distribution we find that the current is much smaller (than the binary mass case with the same value of ∆) but it is not clear whether all states get localized and if an insulating phase exists.
Acknowledgements: We thank G. Baskaran, Michael Aizenman, Tom Spencer and especially David Huse for useful discussions. We also thank Srikanth Sastry and Vishwas Vasisht for use of computational facilities. The research of J. L. Lebowitz was supported by NSF grant No. DMR0802120 and by AFOSR grant No. FA9550-07.

APPENDIX A: KINETIC THEORY
Kinetic theory becomes valid in the limit of small disorder. Its basic object is the Wigner function, f , which describes the phonon density in phase space and is governed by the transport equation Here r ∈ R d (boundary conditions could be imposed), k ∈ [−π, π] d is the wave number of the first Brioullin zone, ω is the dispersion relation of the constant mass harmonic crystal, and C is the collision operator. It acts only on wave numbers and is given by We refer to [40] for a derivation. In the range of validity of (A1), (A2) we can think of phonons as classical particles with energy ω and velocity ∇ω(k). They are scattered by the impurities from k to dk with the rate Collisions are elastic. We distinguish (i) no pinning potential. Then for small k one has ω(k) = |k| and |∇ω(k)| = 1. From (A2) the total scattering rate behaves as |k| d+1 . This is the basis for the discussion in connection with Eq. (2).
The prefactor in (A2) can be replaced by ω 2 0 . The velocity is k and the scattering is isotropic with rate |k| d−2 . Thus the diffusion coefficient results as D(k) ∼ = |k| −d+4 which vanishes as |k| → 0 for d = 2, 3. Hence there is no contribution to the thermal conductivity from the small k modes.

APPENDIX B: TRANSFER MATRIX APPROACH
We now outline steps by which T N can be expressed in forms which are amenable to accurate numerical evaluation.
We will give results whereby we express T N in terms of product of random matrices. These are related to the Green's function and transfer matrix methods used earlier in the calculation of localization lengths in disordered electronic systems [41]. Some related discussions for the phonon case can be found in [42]. For heat conduction in one dimensional disordered chains, the transfer matrix approach has been shown to be very useful in obtaining analytic as well as accurate numerical results and here we study the extension of this to higher dimensions.
The transmission coefficient is given by and we now specify the form of S + L,R corresponding to the equations of motion in Eqs. (4). Note that we have transformed to dimensionless variables ω → ω/Ω, M → M/m, V → V/k, γ → γ/(mΩ) where Ω = (k/m) 1/2 . We are considering heat conduction in the ν = 1 direction of a d-dimensional lattice with particles on the layers n 1 = 1 and n 1 = N being connected to heat baths at temperatures T L and T R respectively. The matrices S + L and S + R represent the coupling of the system to the left and right reservoirs respectively, and can be written as N × N block matrices where each block is a N × N matrix. The block structures are as follows: where I is a N × N unit matrix, and 0 is a N × N matrix with all elements equal to zero. Similarly the matrices M and V have the following block structure: where M n denotes the diagonal mass-matrix for the n 1 = n layer and Φ is a force-constant matrix whose off-diagonal terms correspond to coupling to sites within a layer. Hence the matrix G −1 = [−Mω 2 + V − S + L − S + R ] has the following structure: where a l = −M l ω 2 + Φ − δ l,1 Σ + L − δ l,N Σ + R . Now defining Γ L,R = Im[Σ + L,R ] and with the form of S + L,R given in Eqs. (B1),(B2), we find that the expression for the transmission coefficient reduces to the following form: We now show that G + N satisfies a simple recursion equation. We first introduce some notation. Let Y (l,l+n−1) with 1 ≤ n ≤ N − l + 1 denote a n × n tridiagonal block matrix whose diagonal entries are a l , a l+1 , ...a l+n−1 , where each a l is a N × N matrix. The off-diagonal entries are given by −I. For an arbitrary block matrix A (l,m) , A (l,m) (i,j) will denote the block sub-matrix of A (l,m) beginning with i th block row and column and ending with the j th block row and column, while A (l,m) i,j will denote the (i, j) th block element of A (l,m) . Also I n will denote a n × n block-diagonal matrix with diagonal elements I. The inverse of Y (1,N ) is denoted by [Y (1,N ) ] −1 = G (1,N ) and satisfies the equation: According to our notation we have G (1,N ) = G + and G (1,N ) N ) has the following structure: where W T N = (0, 0, ..., −I) is a 1 × N − 1 block vector. We then write Eq. (B6) in the form is a 1×N −1 block vector. From Eq (B7) we get the following four equations: Noting that [Y (1,N −1) ] −1 = G (1,N −1) we get, using the third equation above and the form of W N : The initial conditions for this recursion are: P (1,0) = I M and P (1,1) = a 1 . By proceeding similarly as before we can also obtain the following recursion relation: and P (1,N ) can be recursively obtained using the initial conditions P (N +1,N ) = I M and P (N,N ) = a N . Given the set {a i }, by iterating either of the above equations one can numerically find P (1,N ) and then invert it to find G N ) 1,N . However this scheme runs into accuracy problems since the numerical values of the matrix elements of the iterates grow rapidly. We describe now a different way of performing the recursion which turns out to be numerically more efficient. We first define From Eq. (B12) we immediately get: with the initial condition r 1 = a 1 . Then G is evaluated turns out to be numerically more accurate.
Finally we show that one can express G (1,N ) 1,N in the form of a product of matrices. The product form is such that the system and reservoir contributions are separated. First we note that the form of the matrices a l for our specific problem is: a l = c l − δ l,1 Σ 1 − δ l,N Σ N where c l = −M l ω 2 + Φ. We define system-dependent matrices Q (1,n) , Q (n,N ) by replacing a 1 , a N by c 1 , c N in the recursions for P s'. Thus Q (1,n) = P (1,n) (a 1 → c 1 , a N → c N ) and Q (n,N ) = P (n,N ) (a 1 → c 1 , a N → c N ). Clearly Qs' satisfy the same recursion as the P s' with a l replaced by c l . Then using Eqs. (B12),(B13), and similar equations for the Qs' we get: From the recursion relations for the Qs' it is easy to see that    (ii) Plot of N T (ω). The range of frequencies for which T (ω) ∼ 1/N is indicated by the dashed line. (iii) Plot of ρ(ω) for binary mass ordered and single disordered sample. (iv) Plot of N 2 P −1 for single samples (smoothed data). We see that even though the allowed normal modes occur over a large frequency band ≈ (0 − 12), transmission takes place in a small band ≈ (0 − 1.25) and is negligible elsewhere. The IPR plots confirm that the non-transmitting states correspond to localized modes. In (i) we see that ω L c is slowly decreasing with inrease of N .