Beating the amorphous limit in thermal conductivity by superlattices design

# Beating the amorphous limit in thermal conductivity by superlattices design

## Abstract

The value measured in the amorphous structure with the same chemical composition is often considered as a lower bound for the thermal conductivity of any material: the heat carriers are strongly scattered by disorder, and their lifetimes reach the minimum time scale of thermal vibrations. An appropriate design at the nano-scale, however, may allow one to reduce the thermal conductivity even below the amorphous limit. In the present contribution, using molecular-dynamics simulation and the Green-Kubo formulation, we study systematically the thermal conductivity of layered phononic materials (superlattices), by tuning different parameters that can characterize such structures. We discover that the key to reach a lower-than-amorphous thermal conductivity is to block almost completely the propagation of the heat carriers, the superlattice phonons. We demonstrate that a large mass difference in the two intercalated layers, or weakened interactions across the interface between layers result in materials with very low thermal conductivity, below the values of the corresponding amorphous counterparts.

## I Introduction

Materials with low thermal conductivity, , are employed in many modern technologies, such as thermal management in electronic devices or thermoelectric energy conversion et al. (2001, 2009a); Maldovan (2013). In general, low values of are observed in disordered solids Goodson (2007), including topologically disordered systems (glasses) and crystalline solids with size or mass disorder (disordered alloys) Cahill and Pohl (1988); et al. (1992); Allen and Feldman (1993); et al. (2013a, ). This behaviour can be rationalized by considering the phenomenological kinetic theory expression Kittel (1996) , which relates average velocity, , and lifetime, , (and therefore the mean free path ) of phonons to ( is the specific heat per unit volume). In good crystals, phonons lifetime is primarily controlled by anharmonic interactions. In contrast, in disordered solids, the disorder (or the elastic heterogeneity et al. (2013b)) reduces (or ) and, as a result, .

In early experimental investigations Cahill and Pohl (1988); et al. (1992), Cahill et al. have studied the disordered alloys, e.g., , and shown that can be reduced to the glass value by controlling the relative composition . In our works et al. (2013a, ) we in turn demonstrated that, in size-disordered crystal, progressively decreases with increasing size mismatch, eventually converging to the corresponding glass value. When this limit is reached, is comparable to the time scale of thermal vibrations ( to the particle size), i.e., to the minimum time (length) scale et al. (2013a, ). Heat propagation can therefore be described as a random walk of vibrational energies Cahill and Pohl (1988); et al. (1992), or in terms of non-propagating delocalized modes, the diffusons Allen and Feldman (1993). For this reason, the value in the glass is generally considered as a lower bound for of materials with homogeneous chemical composition Cahill and Pohl (1988); et al. (1992).

A crucial issue Goodson (2007) is whether thermal conductivity can be lowered below the glass limit through nanoscale phononic design et al. (2011a); Maldovan (2013). This possibility would allow to devise (meta-)materials which are excellent thermal insulators while preserving good electronic properties, as needed in many applications et al. (2001, 2009a); Maldovan (2013). The most popular design to reach this goal is that of a lamellar superlattice et al. (1997, 2000a, 1999); Daly and Maris (2002); Venkatasubramanian (2000), often composed of two chemically different intercalated layers, e.g., Si-Ge et al. (1997, 2000a) or GaAs-AlAs et al. (1999); Daly and Maris (2002) (see also Fig. 1). In a superlattice, the thermal conductivity tensor is anisotropic, with the cross-plane component, , usually lower than the in-plane value,  et al. (2002a); et al. (2007a). In recent experiments  et al. (2004, 2007b, 2010), ultra-low values of , suggested to be smaller than the amorphous limit, were measured. In particular, Costescu et al. et al. (2004) demonstrated that the presence of a high-density of interfaces decreases of W- nanolaminates, below that of the amorphous . An experiment by Chiritescu et al. et al. (2007b) achieved ultra-low thermal conductivity in layered crystals, by disordering the crystalline sheets. Finally, Pernot et al. et al. (2010) also observed very low values of , below that of amorphous Si, in Ge nanodots multi-layers separated by Si crystals.

Although the above works have demonstrated very low values of in superlattice systems, we note that these have not been systematically compared to the values assumed in the glasses with exactly the same chemical composition. Also, a general framework to rationalize in a coherent single picture all these observations is, to the best of our knowledge, still lacking.

In this work, we address these two issues. Building on the comparison of the superlattice with the corresponding amorphous structure, we clarify the mechanisms allowing for ultra-low thermal conductivity in the former. We have studied by computer simulation a numerical model that allows one to exactly compare ordered and disordered systems with identical chemical composition and access detailed information on the entire normal modes spectrum, providing, as a consequence, a complete understanding of the heat transfer process. As the lifetime of heat carriers is already minimum in glasses et al. (2013a, ), we demonstrate that the key to even lower thermal conductivities is to suppress their propagation across the interfaces between the constituent layers.

More in details, we have focused on three distinct design principles for superlattices, mimicking similar configurations actually employed in experiments. These are based on the face-centered-cubic (FCC) lattice structure, and are composed of: (S1) two intercalated crystalline layers formed by sphere particles with different masses; (S2) ordered crystalline layers intercalated to mass-disordered alloy layers; and (S3) identical crystalline layers with modified (weakened) interactions across the interfaces (see the Methods section and Table 1). We show that a large mass difference between layers (S1) and weakened interactions between layers (S3) efficiently obstruct the propagation of phonons, resulting in a very large reduction of the superlattice thermal conductivity, even below the values pertaining to the glass phases with identical composition. Based on our results, we conclude with a discussion of the optimal strategy to follow towards very low thermal conductivity materials.

In Fig. 1 we show a schematic illustration of a superlattice composed by two intercalated layers, and , both of thickness . The competition between two length scales, the repetition period of the superlattice, , and the mean free path of the superlattice phonons, , determines the coherent or incoherent character of phonon transport, as described in Simkin and Mahan (2000); Yang and Chen (2003); Garg and Chen (2013) and demonstrated by numerical simulations et al. (2005, 2007c, 2008a) and recent experiments et al. (2014a).

For , the incoherent phonon transport is independent in the different layers, and phonons can be effectively treated as particles. In this case, the Boltzmann transport equation applies Chen and Neagu (1997); Chen (1998), and the particle-like phonons are scattered within the layers (internal resistance) and at the interfaces (interfacial resistance) et al. (2000b, 2012a). The thermal conductivity in the cross-plane direction can be written as

 κCP=2κ−1A+κ−1B+4Rw−1=κ∞CP(11+ℓKw−1), (1)

where

 κ∞CP=limw→∞κCP=2κAκBκA+κB. (2)

Here, and are the thermal conductivities of materials and , and is the Kapitza length Nan and Birringer (1998); Barrat and Chiaruttini (2003). is the interfacial resistance, which exists even at a perfect interface and depends on the nature of the contacting materials (e.g., crystal-crystal, crystal-glass) et al. (2000b, 2012a). For (), the interfacial resistance is relatively large (small) compared to the internal resistance. Both and (the in-plane thermal conductivity) increase with , due to the decrease of the interfacial resistance density Chen and Neagu (1997); Chen (1998). In the diffuse limit , where the interfacial resistance can be neglected, and have the upper bounds and , respectively.

When , phonon transport is coherent Simkin and Mahan (2000); Yang and Chen (2003); Garg and Chen (2013); et al. (2005, 2007c, 2008a); et al. (2014a), and the wave nature of phonons cannot be neglected. In this regime, decreases with increasing , in contrast with the incoherent case. The reduction of is explained with the emergence of a band gap at the Brillouin zone boundary, due to band-folding et al. (1988); Mizuno and i. Tamura (1992): increasing augments the frequency gap in the dispersion relation. This, in turn, decreases the average group velocity of phonons, finally reducing . Mini-umklapp processes Ren and Dow (1982), occurring at the mini-Brillouin zone, also contribute to the reduction of . At the crossover length , between the incoherent and the coherent transport regimes, assumes a minimum value when plotted against  Simkin and Mahan (2000); Yang and Chen (2003); Garg and Chen (2013); et al. (2005, 2007c, 2008a); et al. (2014a). We have encountered this situation in the case of superlattice S1, as we will see below.

Details of the structure of the interface between layers are also known to significantly affect phonon transport et al. (2002b, 2003a, 2003b); Landry and McGaughey (2009); et al. (2013c); et al. (2009b, 2011b); et al. (2011c, d, 2012b, 2013d). It has been reported that interfacial roughness et al. (2002b, 2003a, 2003b) or mixing Landry and McGaughey (2009); et al. (2013c) reduce both and , and can even suppress the coherent nature of phonons, with increasing monotonously at any . The interface topology is also an important factor to determine the phonon transport et al. (2009b, 2011b). While we will not address precisely this situation in detail here, the superlattice S2 of our study bears some similarities with it.

Finally, the stiffness of interfacial bondings, which can be controlled by applying pressure et al. (2011c, d) or tuning chemical bonding et al. (2012b), has significant effects on heat transport features, which will be demonstrated by the study of the S3 superlattice.

## Ii Results

In Table 1, we present the details of the three superlattice systems studied in this work, with values of the important quantities: and are the thermal conductivities of layers and , respectively; and are the cross- and in-plane diffuse limits of and ; is the interfacial resistance, the Kapitza length; and are the thermal conductivities of the glass and disordered alloy with exactly the same composition as the indicated superlattice. Thermal conductivities have been estimated by molecular-dynamics (MD) simulation and the Green-Kubo formulation McGaughey and Kaviany (2006); et al. (2008b). The number density and the temperature are fixed at (the corresponding crystal lattice constant is ) and , respectively. Vibrational states were also characterized by using a standard normal-modes analysis Ashcroft and Mermin (1976). Details about the systems and the methods used for the simulation production runs and analysis are given in the Methods section.

S1. Superlattice composed of two intercalated crystalline layers with different masses. In Fig. 2 we show the thermal conductivities, and (symbols), as functions of the replication period, , for the layers mass ratios , and . The values of the diffuse limits and as well as those of the glass and the disordered alloy constituted by the same species (see Table 1) are also shown as lines. As expected, the relation holds for the pure materials. In the studied -range, to (monolayers), the in-plane value shows a very weak dependence on , as was observed for superlattices with perfect interfaces in Refs. et al. (2007c, 2003b). The value of is close to, although lower than, , indicating that slight in-plane phonon scattering at the interface is still active.

More interestingly, as increases, the cross-plane value decreases steeply, reaches a minimum value at , and next increases mildly at larger . This -dependence is consistent with previous predictions Simkin and Mahan (2000); Yang and Chen (2003); Garg and Chen (2013); et al. (2005, 2007c, 2008a); et al. (2014a), and corresponds to the crossover at from coherent to incoherent phonon transport. In the incoherent regime, , from Eq. (1) and the data of (dashed line in Fig. 2) we can extract the values of the interfacial resistance, , and the Kapitza length, , which are presented in Table 1. Note that for (Fig. 2(c)), we do not observe a clear thermal conductivity minimum. More precisely, even at the largest value , is still orders of magnitude lower than , indicating that the interfacial resistance results in a strong reduction of in this range of . Equivalently, the Kapitza length is significantly larger than the maximum period . The data shown in Fig. 2 demonstrate that can be indeed lowered below the disordered alloy limit for , and even below the glass limit for higher mass heterogeneities, and . These results are consistent with the experimental work of Ref. et al. (2004), and demonstrate that the interface formed between dissimilar materials effectively reduces . It is also worth noting that the thermal conductivity tensor is very strongly anisotropic in this case, with .

The vibrational modes of the structure, i.e., the superlattice phonons, are key to understand the above behaviour of thermal conductivity. In Fig. 3 we show the vibrational density of states (vDOS), , for and to . and of the bulk crystals of type and . The vDOS of the glass and of the disordered alloy are also shown for comparison. Note that . At small , of the superlattice roughly follows that of the disordered alloy, implying that the vibrational states in the two layers are strongly mixed. In this situation, phonons are able to propagate in both the cross- and in-plane directions. On the other hand, as increases, generates features increasingly similar to those identifying and , separately. In particular, in the low- region follows (the heavy crystal ), whereas (the light crystal ) controls in the high- region. This result indicates that different parts of the vibrational spectrum are active in the two layers, with high(low)- modes preferentially excited in the light (heavy) layer (). In this situation, phonon propagation is largely obstructed in the cross-plane direction, leading to the observed large reduction of . We remark that phonons propagate in the in-plane direction with few constraints, as shown by the large value of close to . This implies that phonons, whose propagations are blocked in the cross-plane direction, are actually specularly reflected at the interface and confined in the in-plane direction.

The separation of the vibrational states found in the becomes more clear when considering the vibrational amplitudes associated with the eigenstates . In Fig. 4 we show the vibrational amplitudes, and (Eq. (6)), in the two layers and for each mode , together with the binned average values (solid lines). Based on the relations and , we can define a relative degree of excitation of particles in the two layers, by the threshold value : large excitations correspond to , small excitations to . If , particle vibrations in both layers are of the same degree and correlated.

At small we find, particularly in the low- region, a large fraction of vibrational states with . As increases, in the high frequency region , where is the high-frequency boundary in , only particles in the light layer vibrate (), whereas those in the heavy layer are almost immobile, as indicated by . In this -region, phonon propagation in the cross-plane direction is therefore almost completely suppressed. On the other hand, for , particles pertaining to the heavy layer show large vibrational amplitudes (), while vibrations in layer tend to be small (). More in details, for , we see that the averaged amplitudes are much larger in the layer () than in the layer () in the range. Contrary to the case of , however, a significant number of modes are excited in both layers and , even with . We therefore conclude that, for , some phonons still propagate in the cross-plane direction, contributing to .

We note that our observation of the vibrational separation in both the vDOS and vibrational amplitudes is consistent with results reported previously et al. (2008a, 2013c, 2012c). Indeed, the simulation work of Ref. et al. (2008a) reported a separation in the vDOS of the Si isotopic-superlattice (- superlattice). A recent simulation work et al. (2013c) focused on partial inverse participation ratios in a superlattice similar to the one considered here, reporting vibrational modes separation between layers. Ref. et al. (2012c) attributed the reduction of thermal conductivity to a mechanism described as phonon localization, which we consider to be essentially the same phenomenon as the vibrational separation described here.

We believe that this concept of vibrational separation is a simple and accurate framework to rationalize the behaviour of thermal conductivity in superlattices. In particular, it provides a complete characterization of the minimum in the -dependence of . Indeed, in the range to identifying the coherent regime, the vibrational separation hinders the coherent phonon propagation in the cross-plane direction, leading to the large reduction of . In contrast, in-plane phonon propagation is very mildly affected by the vibrational separation and, therefore, keeps high values. Also, by considering and (solid lines), we conclude that the separation saturates to its maximum level at . Upon further increase , although averaged values show no significant changes, we recognize an increasing fraction of modes with and for (panels (e) and (f) in Fig. 4). This observation indicates that the separation tendency for modes with and becomes weaker, i.e., the correlation of vibrational features in the two layers decreases, which corresponds exactly to the incoherent transport picture, and leads to the increase of . Although transport becomes completely incoherent only for values of of the order of the Kapitza length (note that for ), this feature appears as soon as the vibrational separation is saturated, at the crossover point . Thus, the saturation point of the vibrational separation identifies the minimum value of , which can be indeed below the glass limit.

S2. Superlattice composed of intercalated ordered crystalline layers and mass disordered layers. This system consists of three components, with masses in the crystalline layer , and and in the disordered alloy layer . In Fig. 5, we plot and for the mass ratios of the layer , , , and . At small , the values of both and are very close to those of the disordered bulk alloy formed by the same particles. As increases, increases gradually toward . This increase is controlled by the development of in-plane phonon propagation in the ordered crystalline layer . Indeed, the of the superlattice, shown in Fig. 6, roughly follows that of the disordered bulk alloy at small , whereas at large it is dominated by . In particular, the longitudinal peak around becomes clear, corresponding to that of the crystalline layer .

The cross-plane value also increases with , but reaches the limit value already at . Since of the disordered alloy layer is low (see Table 1), remains low, typically less than twice the disordered alloy value. As a result, the variation of with is small. This result indicates that scattering in the disordered alloy layer dominates the thermal conduction in the cross-plane direction. Both experimental work et al. (2003c) on Si(crystal)-SiGe(disordered alloy) nanowires and numerical simulations Landry and McGaughey (2009) have reported similar observations. We also note that the coherent nature of the superlattice phonons in the cross-plane direction, which we observed in the S1 system, breaks down in S2. This is essentially equivalent to previous findings that disorder in interfacial roughness et al. (2002b, 2003a, 2003b), or interfacial species mixing Landry and McGaughey (2009); et al. (2013c) destroy the coherent features of vibrational excitations present in the investigated superlattices. In addition, the thermal conductivity tensor becomes increasingly anisotropic at larger due to the increase of , showing a behaviour different than that observed in S1. As a consequence of these features, in superlattices of type S2 the variability of the cross-plane heat transport is strongly bounded, and the minimum limit of just corresponds to the disordered alloy limit, i.e., cannot be reduced below the glass limit.

S3. Superlattice composed by identical crystalline layers separated by weakly interacting interfaces. In Fig. 7 we show the -dependences of and for the case where the energy scale associated to particles interactions across the interfaces () are lowered compared to those intra-layers, with and in the two panels. In the figure, we also plot as lines the data for the corresponding one-component crystal and glass with unmodified interactions. The in-plane value is almost independent of , and is very close to the value pertaining to the crystal. In contrast, decreases monotonically by increasing , and especially in the weaker case , the observed reduction of is dramatic. At , equals the value obtained for the glassy sample, and it is almost two orders of magnitude lower than this value at . This extremely low is consistent with previous experimental work et al. (2007b).

Some insight about the origin of this observation comes from the data shown in Fig. 8, where we display the average cross-plane distance between adjacent crystalline planes (monolayers), normalized to the value in the perfect lattice, . For and , the system keeps the perfect lattice structure, with for all monolayers. In contrast, as decreases and for a large value , becomes substantially larger than at the interfaces, which therefore assumes a local density lower than the average. At the same time, slightly reduced are also observed for the other intra-monolayers, leading to an increase of the local density compared to the average. This heterogeneity hinders energy propagation across the interface and, as a result, phonons are specularly reflected and confined in the in-plane direction. We remark that in the cases with , the values of at the interfaces located at and are different, with a large discrepancy for . We rationalize this behaviour by observing that, during the preparation stage of the sample, the applied selective weakening of the interactions destabilizes the global equilibrium of the superlattice, with a concentration of mechanical stress close to the interfaces. Lattice planes far from the boundaries easily recover mechanical equilibrium by coherently reducing their mutual distance. In contrast, particles in monolayers adjacent to the interfaces move both out-of-plane and in-plane, to optimize the local effective spring constants. The optimal solution found depends in general on the details of the local environment, explaining the observed discrepancy in at different interfaces.

The behaviour of can be further elucidated by inspection of the main features of the vibrational spectrum. In Fig. 9 we plot the of superlattice S3, together with the vibrational amplitudes and . The shows transverse and longitudinal phonon branches for all cases, similar to the homogeneous bulk crystal. As increases deforms, following the appearance of an increasing fraction of modes at increasing higher frequencies. This behaviour is certainly correlated to the observation made above (see Fig. 8) for , that the distance between monolayers far from the interfaces becomes smaller than . The consequent larger mass density makes higher the frequency of phonon modes of given wavelength, leading to the shift of towards higher frequencies. This global shift has as a consequence a mild increase of with , as it is clear from Fig. 7(b) ( case). Note that for and (Fig. 9(a)), shows an excess of lower- modes compared to those present in the one-component crystal, simply due to the weakened interactions at the interfaces.

We now focus on the vibrational amplitudes, and (Fig. 9, bottom panels). In the cases with and and , the particles in the two layers and show completely equivalent and correlated vibrations for the vast majority of the modes, as indicated by . This result implies that phonons indeed propagate across the weakened interfaces in the cross-plane direction, but they are also partially reflected at the interface, causing the observed reduction of . The situation changes drastically in the case and , where the ultra-low value of can be reached. Except for the low- modes, and are symmetrically randomly distributed around the average values , indicating that particles in layers and vibrate independently, in an uncorrelated manner. As a consequence, a very large fraction of vibrational modes do not cross at all the interfaces, but rather undergo a perfect specular reflection. In this situation, heat is not transferred between two adjacent layers and , leading to extremely low value of , while keeping a high . We conclude by noticing that although specular reflection was also observed in the system S1, the physical mechanism behind this phenomenon is different in the two cases: vibrational separation causes reflection in the former, whereas weakened interactions across the interface, with the resulting augmented spacing between the layers, completely block cross-plane phonon propagation in the latter.

## Iii Discussion

We have provided numerically, for the first time to our knowledge, a clear demonstration of very low thermal conductivities in superlattices, below the glassy limit of the corresponding amorphous structures. Blocking phonon propagation in ordered structures via interfaces design is the key principle. We have identified two possible strategies to achieve this goal: imposing a large mass heterogeneity in the intercalated layers (as in system S1) or degrading inter-layers interactions compared to those intra-layers (as in S3). We have found that in both cases phonons are specularly reflected at the interface and confined in the in-plane direction. This reduces the cross-plane thermal conductivity below the corresponding glass limit, while keeping the in-plane contribution close to the pure crystalline value.

More specifically, in the case of mass mismatch (S1), propagation of phonons with high frequencies () is almost completely suppressed, whereas a fraction of low-frequency phonons () are still able to propagate across the interfaces, contributing to (Fig. 4(d)). Also, the minimum in thermal conductivity as a function of the repetition period (Fig. 2) corresponds to a maximum in the vibrational separation between the layers of type and . These therefore act as true filters in complementary regions of the vibrational spectrum, suppressing significantly phonons transport in the direction of the replication pattern. On the other hand, attenuated interactions across the interfaces (S3) are able to block phonons at almost all frequencies (see Fig. 9(c)), which results into extremely low values of , even orders of magnitude lower than the corresponding glass limit (Fig. 7(b)). In this sense, directly modifying the interfaces seems to be the most effective strategy to obtain very low heat transfer. Note that this is a practically feasible route, since attenuated interfaces can be designed by exploiting materials with weak van der Waals forces among adjacent crystalline planes, as demonstrated in the case of sheets in Ref. et al. (2007b). Interfaces stiffness modification by controlling pressure et al. (2011c, d) or chemical bonding et al. (2012b) are additional possible routes to directly tune the strength of interfaces.

Our data also suggest that intercalating disordered alloy layers in ordered crystalline layers (S2) is not effective in lowering . Indeed, we have demonstrated that in this case disorder is not sufficient to block the propagation of vibrational excitations, even though it makes phonons lifetimes short. The intercalated disordered alloy layer dominates phonon transport in the entire superlattice, notwithstanding the presence of the crystalline layers. As a result, thermal conductivity is very similar to the one of the disordered alloy and is only marginally modified by modulation of the period (see Fig. 5). Also, as suggested in previous works, disorder in the interfacial roughness et al. (2002b, 2003a, 2003b) or interfacial mixing Landry and McGaughey (2009); et al. (2013c) seems to already dominate over phonon transport, and destroy the coherent nature of phonons.

In addition, as we understand from our analysis of vibrational amplitudes (Figs. 4 and 9), it is much more problematic to block low- (long wavelength, ) phonons propagation, than those with high- (short ). This situation is similar to what has been observed in bulk glasses, where the long- acoustic waves are not scattered by the disorder and can propagate over long distances by carrying heat energy Monaco and Mossa (2009); et al. (2014b). Therefore, blocking or efficiently scattering the long- phonons is also a key factor to achieve very low thermal conductivities, as was pointed out in Ref. et al. (2012d). A possibility to realize this task is embedding in the targeted material objects featuring larger typical sizes, including nano-particles et al. (2006); Zhang and Minnich (2014) or nano(quantum)-dots et al. (2010, 2011e). Based on this strategy, very low thermal conductivity was achieved experimentally in a Si-Ge quantum-dot superlattice et al. (2010), even below the amorphous Si value. The additional possibility of introducing large size defects by the porous structuring of materials has also been explored in a recent numerical work et al. (2014c). Here, values of thermal conductivity times smaller than that of bulk Si were reached in Si phononic crystals with spherical pores.

In conclusion, we note that the three superlattice structures studied in the present work show totally different -dependences of cross and in-plane thermal conductivities. Our results therefore not only contribute to a deeper comprehension of the physical mechanisms behind very-low thermal conductivity, they also provide insight for developing new design concepts for materials with controlled heat conduction behaviour.

## Iv Methods

System description. In this Section we provide details on the numerical models we have used for the superlattices. The corresponding amorphous structures (glasses) and disordered alloys with exactly the same composition were also prepared, for the sake of comparison with superlattice phases. We have considered in all cases a 3-dimensional cubic box, of volume ( being the linear box size), with periodic boundary conditions in all directions. In the superlattice and disordered alloy cases, particles were distributed on the FCC lattice sites. In the glass phases, they were frozen in topologically random positions following a rapid quench from the normal liquid phase below the glass transition temperature , avoiding crystallization (see, for instance, Ref. Monaco and Mossa (2009) for details on the preparation of glasses). Particles, and , interact via soft-sphere (SS) or Lennard-Jones (LJ) potentials:

 vijSS(r) =ϵij(σijr)12, (3) vijLJ(r) =4ϵij[(σijr)12−(σijr)6],

where is the distance between those two particles, and and are the interparticle diameter and interaction energy scale, respectively. The potential is cut-off and shifted at . Particle has mass , and we have used , ( is the Boltzmann constant), and as units of length, temperature, and mass. As a reference, for Argon ,  K, and  a.m.u. We considered the number density , corresponding to a lattice constant .

We prepared three superlattices, composed of intercalated FCC lattice layers, and , both of thickness , as schematically illustrated in Fig. 1. The first superlattice (S1) consists of two crystalline layers formed by sphere particles with different masses, and . We have considered mass ratios , while keeping a constant average mass . As an example, the case corresponds to and . We have dubbed and as the light and heavy layers, respectively. Note that a mass ratio of corresponds to the case of the realistic Si-Ge superlattice. Except for the above mass difference in the different layers, all particles are characterized by the same properties. In particular, they interact via the SS potential , with .

The second superlattice (S2) is composed of an ordered crystalline layer intercalated to a disordered alloy layer . in , whereas in half of the particles have mass , the others, and are randomly distributed on the lattice sites. Again, and are determined by the mass ratio , keeping a constant average value . All particles in both layers interact via the SS potential with .

The third superlattice (S3) is composed of identical crystalline layers and , but the interactions among particles in different layers (i.e., across the interfaces) are modified (weakened) compared to those intra-layers. All particles have mass , and interact via the LJ potential , with . The energy scale of interactions between particles pertaining to different layers are, however, reduced to .

MD simulation and the Green-Kubo method for the calculation of thermal conductivity. In the present study, all simulations have been realized by using the large-scale, massively parallel molecular dynamics simulation tool LAMMPS Plimpton (1995); lam (). The systems were first equilibrated at relatively low temperature by MD simulation in the -ensemble. This choice was dictated by the need to reduce anharmonic effects, in order to primarily focus on the contribution of the structural features of the superlattices on thermal conductivity. We must note that our approach is classical, and does not take into account the quantum mechanisms active in the low- regime Kittel (1996). These effects have important implications, increasing the contribution to the thermal conductivity coming from low- vibrational excitations. At present, however, it is not obvious and still under debate how to effectively include quantum effects into a classical system et al. (2009c, 2014d), and we have therefore chosen to stay within a fully classical approach.

Following the equilibration stage, we performed the production runs in the -ensemble. The Green-Kubo formulation McGaughey and Kaviany (2006); et al. (2008b) was next applied to calculate the thermal conductivities, in the cross-plane and in-plane directions, respectively:

 κCP =1VT2∫∞0⟨Jz(t)Jz(0)⟩dt, (4) κIP =12VT2∫∞0⟨Jx(t)Jx(0)+Jy(t)Jy(0)⟩dt.

Here, , and are the heat currents in the in-plane () and cross-plane () directions, and denotes the ensemble average. In the bulk glasses and disordered alloys, , i.e., heat conduction is isotropic, whereas in the superlattices, they are expected to assume different values et al. (2002a); et al. (2007a). Landry et al. et al. (2008b) have carefully confirmed the validity of the Green-Kubo method for the calculation of superlattices thermal conductivity, by comparison with the direct method based on non-equilibrium simulation. Also, in the Green-Kubo calculations, one must be attentive to finite system size effects McGaughey and Kaviany (2006); et al. (2008b). Indeed, long-wavelength phonons with are excluded from the simulation box due to the finite value of the box size, which imposes important size effects on the numerical determination of . The box size therefore needs to be large enough to include a vibrational spectrum sufficient to establish an accurate description of anharmonic coupling (scattering) processes McGaughey and Kaviany (2006). We note that the considered is low enough to substantially reduce anharmonic effects, but anharmonic couplings are still present.

We can take care of finite size effects by increasing to values where and become -independent. For the glass and disordered alloy thermal conductivities, we have confirmed that a system size () is sufficiently large to obtain correct values of , without any size effect et al. (2013a, ). In the superlattice cases, the appropriate depends on the considered structure and the periodic repetition length  et al. (2008b). More in details, we paid particular attention to the number of repetitions, defined from , necessary to produce sufficient anharmonic couplings of phonons in the cross-plane direction. We have therefore investigated the presence of finite size effects by analyzing different systems with sizes ranging from ( monolayers, ) to ( monolayers, ). In Figs. 2, 5, and 7, we show multiple data points at some -values, obtained for different system sizes. For the S1 superlattice, we confirmed that the required number of repetitions becomes larger for smaller  et al. (2008b): one period () only is adequate for , whereas four periods or more () are required for . We have therefore employed four pattern repetitions () for and two () for . This behaviour is simple to rationalize by inspecting the data in Fig. 2, where the crossover between incoherent and coherent phonon transport occurs around . In the coherent regime , the wave character of the phonons becomes important, and therefore a larger number of repetitions is necessary to produce the coherent wave interference processes correctly. In contrast, smaller values of are needed (even ) in the incoherent regime , where the incoherent particle nature of the phonons appears.

For the S2 and S3 superlattices the system size effects issue is much less pronounced than in the S1 case. We can understand this behaviour by noticing that phonon tranport is mainly determined by the scattering processes in the disordered alloy layer in S2, and the blocking at the weak interface for S3. In both cases the missing long wavelength phonons, with , play very little role in phonon transport and finite system size effects are consequently negligible. We therefore used () for and one or more repetitions () for , for both S2 and S3.

Normal modes analysis. We have characterized the superlattice vibrational states (superlattice phonons) by performing a standard normal-mode analysis Ashcroft and Mermin (1976) with ARPACK arp (). We have diagonalized the dynamical (Hessian) matrix calculated at local minima of the potential energy landscape, and obtained eigenvalues and eigenvectors (polarization vectors) . Here, is the particle index, and is the eigenmode number, where we have disregarded the three vanishing Goldstone modes. The eigenvectors are normalized such that , where is the Kronecker delta function. The eigenfrequencies are next calculated as , and the associated probability distribution (normalized histogram) directly provides the vDOS:

 g(ω)=13N−33N−3∑k=1δ(ω−ωk). (5)

In addition, from the eigenvector we have defined the vibrational amplitudes of mode for layers and :

 EkA(B)=∑j∈layerA(B)(ekj⋅ekj). (6)

Note that for each and, therefore, . Based on the values of and , one can determine in which layer particles are more displaced (excited) according to the associated eigenvector . In particular, if (), particles in layer () contribute more to mode than those in layer (). In the case , particles in both layers contribute equivalently, and in a correlated manner. Note that the normal mode analysis provides us with the system vibrational states in the harmonic limit which, we believe, is an appropriate approximation for our case , where anharmonicities are weak.

###### Acknowledgements.
We thank P. Keblinski for helpful correspondence. This work was supported by the Nanosciences Foundation of Grenoble. J.-L. B is supported by the Institut Universitaire de France. Most of the computations presented in this work were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

### References

1. R. Venkatasubramanian et al., “Thin-film thermoelectric devices with high room-temperature figures of merit,” Nature 413, 597–602 (2001).
2. A. J. Minnich et al., “Bulk nanostructured thermoelectric materials: current research and future prospects,” Energy Environ. Sci. 2, 466–479 (2009a).
3. M. Maldovan, “Sound and heat revolutions in phononics,” Nature 503, 209–217 (2013).
4. K. E. Goodson, “Ordering up the minimum thermal conductivity of solids,” Science 315, 342–343 (2007).
5. D. G. Cahill and R. O. Pohl, “Lattice vibrations and heat transport in crystals and glasses,” Annual Review of Physical Chemistry 39, 93–121 (1988).
6. D. G. Cahill et al., “Lower limit to the thermal conductivity of disordered crystals,” Phys. Rev. B 46, 6131–6140 (1992).
7. P. B. Allen and J. L. Feldman, “Thermal conductivity of disordered harmonic solids,” Phys. Rev. B 48, 12581–12588 (1993).
8. H. Mizuno et al., “Elastic heterogeneity, vibrational states, and thermal conductivity across an amorphisation transition,” EPL 104, 56001 (2013a).
9. H. Mizuno et al., in preparation .
10. C. Kittel, Introduction to Solid State Physics, 7th ed. (John Wiley and Sons, New York, 1996).
11. H. Mizuno et al., “Measuring spatial distribution of the local elastic modulus in glasses,” Phys. Rev. E 87, 042306 (2013b).
12. P. E. Hopkins et al., “Reduction in the thermal conductivity of single crystalline silicon by phononic crystal patterning,” Nano Letters 11, 107–112 (2011a).
13. S.-M. Lee et al., “Thermal conductivity of Si-Ge superlattices,” Applied Physics Letters 70, 2957–2959 (1997).
14. S. Volz et al., “Computation of thermal conductivity of Si/Ge superlattices by molecular dynamics techniques,” Microelectronics Journal 31, 815–819 (2000a).
15. W. S. Capinski et al., “Thermal-conductivity measurements of GaAs/AlAs superlattices using a picosecond optical pump-and-probe technique,” Phys. Rev. B 59, 8105–8113 (1999).
16. B. C. Daly and H. J. Maris, “Calculation of the thermal conductivity of superlattices by molecular dynamics simulation,” Physica B: Condensed Matter 316–317, 247–249 (2002).
17. R. Venkatasubramanian, “Lattice thermal conductivity reduction and phonon localizationlike behavior in superlattice structures,” Phys. Rev. B 61, 3091–3097 (2000).
18. B. Yang et al., “Measurements of anisotropic thermoelectric properties in superlattices,” Applied Physics Letters 81, 3588–3590 (2002a).
19. A. Mavrokefalos et al., “In-plane thermal conductivity of disordered layered and superlattice films,” Applied Physics Letters 91, 171912 (2007a).
20. R. M. Costescu et al., “Ultra-low thermal conductivity in W/ nanolaminates,” Science 303, 989–990 (2004).
21. C. Chiritescu et al., “Ultralow thermal conductivity in disordered, layered crystals,” Science 315, 351–353 (2007b).
22. G. Pernot et al., “Precise control of thermal conductivity at the nanoscale through individual phonon-scattering barriers,” Nature Mater. 9, 491–495 (2010).
23. M. V. Simkin and G. D. Mahan, “Minimum thermal conductivity of superlattices,” Phys. Rev. Lett. 84, 927–930 (2000).
24. B. Yang and G. Chen, “Partially coherent phonon heat conduction in superlattices,” Phys. Rev. B 67, 195311 (2003).
25. J. Garg and G. Chen, “Minimum thermal conductivity in superlattices: A first-principles formalism,” Phys. Rev. B 87, 140302 (2013).
26. Y. Chen et al., “Minimum superlattice thermal conductivity from molecular dynamics,” Phys. Rev. B 72, 174302 (2005).
27. T. Kawamura et al., “An investigation of thermal conductivity of nitride-semiconductor nanostructures by molecular dynamics simulation,” Journal of Crystal Growth 298, 251–253 (2007c).
28. N. Yang et al., “Ultralow thermal conductivity of isotope-doped silicon nanowires,” Nano Letters 8, 276–280 (2008a).
29. J. Ravichandran et al., “Crossover from incoherent to coherent phonon scattering in epitaxial oxide superlattices,” Nature Mater. 13, 168–172 (2014a).
30. G. Chen and M. Neagu, “Thermal conductivity and heat transfer in superlattices,” Applied Physics Letters 71, 2761–2763 (1997).
31. G. Chen, “Thermal conductivity and ballistic-phonon transport in the cross-plane direction of superlattices,” Phys. Rev. B 57, 14958–14973 (1998).
32. E.-K. Kim et al., “Thermal boundary resistance at interface,” Applied Physics Letters 76, 3864–3866 (2000b).
33. E. Lampin et al., “Thermal boundary resistance at silicon-silica interfaces by molecular dynamics simulations,” Applied Physics Letters 100, 131906 (2012a).
34. C.-W. Nan and R. Birringer, “Determining the kapitza resistance and the thermal conductivity of polycrystals: A simple model,” Phys. Rev. B 57, 8264–8268 (1998).
35. J.-L. Barrat and F. Chiaruttini, “Kapitza resistance at the liquid-solid interface,” Molecular Physics 101, 1605–1610 (2003).
36. S. Tamura et al., “Acoustic-phonon propagation in superlattices,” Phys. Rev. B 38, 1427–1449 (1988).
37. S. Mizuno and S. i. Tamura, “Theory of acoustic-phonon transmission in finite-size superlattice systems,” Phys. Rev. B 45, 734–741 (1992).
38. S. Y. Ren and J. D. Dow, “Thermal conductivity of superlattices,” Phys. Rev. B 25, 3750–3755 (1982).
39. B. C. Daly et al., “Molecular dynamics calculation of the thermal conductivity of superlattices,” Phys. Rev. B 66, 024301 (2002b).
40. K. Imamura et al., “Lattice thermal conductivity in superlattices: molecular dynamics calculations with a heat reservoir method,” Journal of Physics: Condensed Matter 15, 8679–8690 (2003a).
41. B. C. Daly et al., “Molecular dynamics calculation of the in-plane thermal conductivity of GaAs/AlAs superlattices,” Phys. Rev. B 67, 033308 (2003b).
42. E. S. Landry and A. J. H. McGaughey, “Effect of interfacial species mixing on phonon transport in semiconductor superlattices,” Phys. Rev. B 79, 075316 (2009).
43. S. C. Huberman et al., “Disruption of superlattice phonons by interfacial mixing,” Phys. Rev. B 88, 155311 (2013c).
44. K. Termentzidis et al., “Nonequilibrium molecular dynamics simulation of the in-plane thermal conductivity of superlattices with rough interfaces,” Phys. Rev. B 79, 214307 (2009b).
45. K. Termentzidis et al., “Cross-plane thermal conductivity of superlattices with rough interfaces using equilibrium and non-equilibrium molecular dynamics,” International Journal of Heat and Mass Transfer 54, 2014–2020 (2011b).
46. W.-P. Hsieh et al., “Pressure tuning of the thermal conductance of weak interfaces,” Phys. Rev. B 84, 184107 (2011c).
47. M. Shen et al., “Bonding and pressure-tunable interfacial thermal conductance,” Phys. Rev. B 84, 195432 (2011d).
48. M. D. Losego et al., “Effects of chemical bonding on heat transport across interfaces,” Nature Mater. 11, 502–506 (2012b).
49. Z. Wei et al., “Negative correlation between in-plane bonding strength and cross-plane thermal conductivity in a model layered material,” Applied Physics Letters 102, 011901 (2013d).
50. A. J. H. McGaughey and M. Kaviany, Advances in Heat Transfer, edited by G. Greene, Y. Cho, J. Hartnett, and A. Bar-Cohen, Vol. 39 (Elsevier, New York, 2006) pp. 169–255.
51. E. S. Landry et al., “Complex superlattice unit cell designs for reduced thermal conductivity,” Phys. Rev. B 77, 184302 (2008b).
52. N. W. Ashcroft and N. D. Mermin, Solid State Physics (Harcourt College Publishers, New York, 1976).
53. L. Yang et al., “Reduction of thermal conductivity by nanoscale 3d phononic crystal,” Nature Scientific Reports 3, 1143 (2012c).
54. D. Li et al., “Thermal conductivity of Si/SiGe superlattice nanowires,” Applied Physics Letters 83, 3186–3188 (2003c).
55. G. Monaco and S. Mossa, “Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale,” Proc. Natl. Acad. Sci. USA 106, 16907–16912 (2009).
56. H. Mizuno et al., “Acoustic excitations and elastic heterogeneities in disordered solids,” Proc. Natl. Acad. Sci. USA 111, 11949–11954 (2014b).
57. M. N. Luckyanova et al., “Coherent phonon heat conduction in superlattices,” Science 338, 936–939 (2012d).
58. W. Kim et al.Phys. Rev. Lett. 96, 045901 (2006).
59. H. Zhang and A. J. Minnich, “The best nanoparticle size distribution for minimum thermal conductivity,” arXiv:1404.1438  (2014).
60. D. L. Nika et al., “Reduction of lattice thermal conductivity in one-dimensional quantum-dot superlattices due to phonon filtering,” Phys. Rev. B 84, 165415 (2011e).
61. L. Yang et al., “Extreme low thermal conductivity in nanoscale 3D Si phononic crystal with spherical pores,” Nano Letters 14, 1734–1738 (2014c).
62. S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” Journal of Computational Physics 117, 1–19 (1995).
63. Http://lammps.sandia.gov.
64. J. E. Turney et al., “Assessing the applicability of quantum corrections to classical thermal conductivity predictions,” Phys. Rev. B 79, 224305 (2009c).
65. O. N. Bedoya-Martinez et al., “Computation of the thermal conductivity using methods based on classical and quantum molecular dynamics,” Phys. Rev. B 89, 014303 (2014d).
66. Http://www.caam.rice.edu/software/ARPACK.
104262