Beating the amorphous limit in thermal conductivity by superlattices design

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 have discovered 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.

A crucial issue 4 is whether thermal conductivity can be lowered below the glass limit through nanoscale phononic design 3,11 . This possibility would allow to devise (meta-)materials which are excellent thermal insulators while preserving good electronic properties, as needed in many applications [1][2][3] . The most popular design to reach this goal is that of a lamellar superlattice [12][13][14][15] , often composed of two chemically different intercalated layers, e.g., Si-Ge 12,13 or GaAs-AlAs 14,15 (see also Fig. 1). In a superlattice, the thermal conductivity tensor is anisotropic, with the cross-plane component, κ CP , usually lower than the in-plane value, κ IP 16,17 .
In recent experiments [18][19][20] , ultra-low values of κ CP , suggested to be smaller than the amorphous limit, were measured. In particular, Costescu et al. 18 demonstrated that the presence of a high-density of interfaces decreases κ CP of W-Al 2 O 3 nanolaminates, below that of the amorphous Al 2 O 3 . An experiment by Chiritescu et al. 19 achieved ultra-low thermal conductivity in layered WSe 2 crystals, by disordering the crystalline WSe 2 sheets. Finally, Pernot et al. 20 also observed very low values of κ CP , below that of amorphous Si, in Ge nanodots multi-layers separated by Si crystals.
Although the above experimental works have demonstrated very low values of κ in superlattice systems, we note that these values have not been systematically compared to those assumed in the glasses with exactly the same chemical composition. Since different chemical species are expected to produce different effects on κ, it is therefore still not completely clear whether superlattice structuration alone can lower κ below the amorphous limit. Note, for instance, that very recent numerical work 21 , has shown that a superlattice composed by layers with randomized thicknesses can indeed show a κ below the value pertaining to the disordered-alloy with the same composition. This limit, however, is generally higher than that in the corresponding glass [5][6][7][8] , which should therefore be considered the true amorphous limit to be beaten.
In addition, a general framework to rationalize in a coherent single picture the previous observations of very low κ 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 8 , 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 point 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, A and B, both of thickness W/2. The competition between two length scales, the repetition period of the superlattice, W, and the mean free path of the superlattice phonons, , determines the coherent or incoherent character of phonon transport, as described in [22][23][24] and demonstrated by numerical simulations [25][26][27] and recent experiments 28 . The superlattice is composed of two FCC-lattice layers, A (red) and B (green). The two layers have identical thickness W/2, where W is the replication period. Here, we measure W as the number of monolayers of the lattice structure, e.g., W = 8 in the displayed case. The distance between adjacent monolayers is a/2 for the perfect FCC structure, where a is the lattice constant. For >  W , 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 29,30 , and the particle-like phonons are scattered within the layers (internal resistance) and at the interfaces (interfacial resistance or Kapitza resistance [31][32][33] ). The thermal conductivity in the cross-plane direction can be written as 22 Here, κ A and κ B are the thermal conductivities of materials A and B, and κ = ∞  R 2 K CP is the Kapitza length 34 . R 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) 31,32 , the interfacial resistance is relatively large (small) compared to the internal resistance. Both κ CP and κ IP (the in-plane thermal conductivity) increase with W, due to the decrease of the interfacial resistance density 29,30 . In the diffuse limit W → ∞, where the interfacial resistance can be neglected, κ CP and κ IP have the upper bounds κ ∞ CP and κ κ κ , respectively. When <  W , phonon transport is coherent [22][23][24][25][26][27][28] , and the wave nature of phonons cannot be neglected. In this regime, κ CP decreases with increasing W, in contrast with the incoherent case. The reduction of κ CP is explained with the emergence of a band gap at the Brillouin zone boundary, due to band-folding 35 : increasing W augments the frequency gap in the dispersion relation. This, in turn, decreases the average group velocity v of phonons, finally reducing κ CP . Mini-umklapp processes 36 , occurring at the mini-Brillouin zone, also contribute to the reduction of κ CP . At the crossover length ∼  W , between the incoherent and the coherent transport regimes, κ CP assumes a minimum value when plotted against W [22][23][24][25][26][27][28] . 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 [37][38][39][40][41][42][43][44][45][46] . It has been reported that interfacial roughness [37][38][39] or mixing 40,41 reduce both κ CP and κ IP , and can even suppress the coherent nature of phonons, with κ CP(IP) increasing monotonously at any W. The interface topology is also an important factor to determine the phonon transport 42,43 . 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 44,45 or tuning chemical bonding 46 , has significant effects on heat transport features, which will be demonstrated by the study of the S3 superlattice.  Table 1. Details of the three superlattice systems investigated in this work. They are based on the FCCcrystal lattice structure and are composed of: (S1) two intercalated crystalline layers (A and B) formed by point particles with different masses m A and m B ; (S2) ordered crystalline layers intercalated to massdisordered alloy layers; and (S3) identical crystalline layers with modified (weakened compared to those intra-layers) interactions across the interfaces. The control parameters are the mass ratio m B /m A in S1, the mass ratio m B2 /m B1 of the disordered alloy layer in S2, and the energy scale ε AB of the interactions across the interfaces in S3. Number density and temperature were fixed to the values ρ= . 1 015 (corresponding to a lattice constant a = 1.58) and T = 10 −2 , respectively. The quantities presented in the table are defined in the main text. In the last column we refer to the figure containing the data relative to the indicated system. Additional details about the investigated superlattices and parameters used are given in the Methods section.

Results
In Table 1, we present the details of the three superlattice systems studied in this work, with values of the important quantities: κ A and κ B are the thermal conductivities of layers A and B, respectively; κ ∞ CP and κ ∞ IP are the cross-and in-plane diffuse limits of κ CP and κ IP ; R is the interfacial resistance,  K the Kapitza length; κ glass and κ alloy are the thermal conductivities of the glass and disordered alloy with exactly the same composition as the indicated superlattice.
The number density of all systems was fixed at ρ = .
1 015, with a corresponding crystal lattice constant a = 1.58. In order to minimize anharmonic couplings and focus primarily on the contribution arising from the details of the nano-structuration, we considered a low temperature value T = 10 −2 .
Thermal conductivities have been estimated by the Green-Kubo formulation 47,48 . We have calculated both components of the superlattices thermal conductivities, κ CP and κ IP , by varying the pattern repetition period W. Note that W indicates the number of monolayers of the lattice structure, where the distance between adjacent monolayers is a/2 for the perfect FCC structure (see also Fig. 1). In the following, we systematically compare the value of κ CP to κ ∞ CP , κ glass , κ alloy , to evaluate the efficiency of the superlattice structures in minimizing heat transfer in the direction of the patterning. The in-plane behaviour has been similarly quantified by comparing κ IP to κ ∞ IP . In addition, we have also characterized the vibrational states by using a standard normal-modes analysis 8,49 .
All important information about the system models and methods used for the simulation production runs and analysis are given in the Methods section. Additional details about specific points are included in the Supplementary Materials. S1. Superlattice composed of two intercalated crystalline layers with different masses. In Fig. 2 we show the thermal conductivities, κ CP and κ IP (symbols), as functions of the replication period, W, for the layers mass ratios m B /m A = 2, 4, and 8. Note that we have chosen to fix the average mass, 〈 m〉 = (m A + m B )/2 = 1, rather than fixing a reference value m A = 1 and varying m B . This latter protocol has indeed an obvious drawback: the average mass would increase when considering different values of the ratio m B /m A , implying a trivial effect on the thermal conductivity which scales as 〈 m〉 −1/2 . This would therefore hinder the possibility to isolate the contribution to lowering the thermal conductivity which originates from the mass mismatch alone.
The values of the diffuse limits κ ∞ CP and κ ∞ IP as well as those of the glass and the disordered alloy constituted by the same species (see Table 1) are also shown as lines in Fig. 2. As expected, the relation holds for the pure materials. In the studied W-range, W = 2 to 40 (monolayers), the in-plane value κ IP shows a very weak dependence on W, as was observed for superlattices with perfect interfaces in Refs 26,39. The value of κ IP is close to, although lower than, κ ∞ IP , indicating that slight in-plane phonon scattering at the interface is still active.
More interestingly, as W increases, the cross-plane value κ CP decreases steeply, reaches a minimum value at  ⁎ W 20, and next increases mildly at larger W. This W-dependence is consistent with previous predictions [22][23][24][25][26][27][28] , and corresponds to the crossover at W * from coherent to incoherent phonon transport. In the incoherent regime, W > 20, from Eq. (1) and the data of κ CP (dashed line in Fig. 2) we can extract the values of the interfacial resistance, R, and the Kapitza length,  K , which are presented in Table 1. Note that for m B /m A = 8 ( Fig. 2(c)), we do not observe a clear thermal conductivity minimum. More precisely, even at the largest value W = 40, κ CP is still orders of magnitude lower than κ ∞ CP , indicating that the interfacial resistance R results in a strong reduction of κ CP in this range of W. Equivalently, the Kapitza length  K is significantly larger than the maximum period , the slope is so small that κ CP appears almost flat when plotted against W, for 20 ≤ W ≤ 40 ( Fig. 2(c)).
The data shown in Fig. 2 demonstrate that κ CP can be indeed lowered below the disordered alloy limit for m B /m A = 2, and even below the glass limit for higher mass heterogeneities, m B /m A = 4 and 8. These results are consistent with the experimental work of Ref. 18, and demonstrate that the interface formed between dissimilar materials effectively reduces κ CP . It is also worth noting that the thermal conductivity tensor is very strongly anisotropic in this case, with κ κ  CP IP . 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), g(ω), for m B /m A = 4 and W = 2 to 80. g A (ω) and g B (ω) of the bulk crystals of type A and B as well as the vDOS of the glass and of the disordered alloy are also shown for comparison. Note that ω ω ( B . At small W = 2, g(ω) 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 W increases, g(ω) generates features increasingly similar to those identifying g A (ω) and g B (ω), separately. In particular, in the low-ω region g(ω) follows g B (ω) (the heavy crystal B), whereas g A (ω) (the light crystal A) controls g(ω) 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 A (B). In this situation, phonon propagation is largely obstructed in the cross-plane direction, leading to the observed large reduction of κ CP . We remark that phonons propagate in the in-plane direction with small constraints, as shown by the large value of κ IP close to κ ∞ IP . 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 g(ω) becomes more clear when considering the vibrational amplitudes associated with the eigenstates k. In Fig. 4 we show the vibrational amplitudes, E A k and E B k (Eq. (6)), in the two layers A and B for each mode k, together with the binned average values (solid lines). Based on the relations , we can define a relative degree of excitation of particles in the two layers, by the threshold value 0.5: large excitations correspond to ≥ .
, particle vibrations in both layers are of the same degree and correlated.
At small W = 2,4 we find, particularly in the low-ω region, a large fraction of vibrational states with . . We therefore conclude that, for ω ω < B max , some phonons still propagate in the cross-plane direction, contributing to κ CP .  Table 1. The solid curve interpolating the κ CP data points in the entire Wrange is a guide for eye. The calculation of the displayed error bars is detailed in the Methods section.
We note that our observation of the vibrational separation in both the vDOS and vibrational amplitudes is consistent with results reported previously 27,41,50 . Indeed, the simulation work of Ref. 27 reported a separation in the vDOS of the Si isotopic-superlattice ( 28 Si-42 Si superlattice). A recent simulation work 41 focused on partial inverse participation ratios in a superlattice similar to the one considered here, reporting vibrational modes separation between layers. Ref. 50 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 W-dependence of κ CP . Indeed, in the range W = 2 to 20 identifying the coherent regime, the vibrational separation hinders the coherent phonon propagation in the cross-plane    = 4), this feature appears as soon as the vibrational separation is saturated, at the crossover point  ⁎ W 20. Thus, the saturation point of the vibrational separation identifies the minimum value of κ CP , which can be indeed below the glass limit.

S2. Superlattice composed of intercalated ordered crystalline layers and mass disordered alloy layers. This system consists of three components, with masses m A = 1 in the crystalline layer A,
and m B1 and m B2 in the disordered alloy layer B. Note that also in this case, on the basis of the same arguments discussed above for S1, the average mass in layer B and that of the entire system are fixed as , respectively, in order to eliminate the trivial contribution associated with different average values in the different cases.
In Fig. 5, we plot κ CP and κ IP for the mass ratios of the layer B, m B2 /m B1 = 2, 4, and 8. At small W ≤ 4, the values of both κ CP and κ IP are very close to those of the disordered bulk alloy formed by the same particles. As W increases, κ IP increases gradually toward κ ∞ IP . This increase is controlled by the development of in-plane phonon propagation in the ordered crystalline layer A. Indeed, the g(ω) of the superlattice, shown in Fig. 6, roughly follows that of the disordered bulk alloy at small W = 4, whereas at large W = 20,40 it is dominated by g A (ω). In particular, the longitudinal peak around ω .  14 5 becomes clear, corresponding to that of the crystalline layer A.
The cross-plane value κ CP also increases with W, but reaches the limit value κ ∞ CP already at  W 20. Since κ B of the disordered alloy layer B is low (see Table 1), κ ∞ CP remains low, typically less than twice the disordered alloy value. As a result, the variation of κ CP with W is small. This result indicates that scattering in the disordered alloy layer B dominates the thermal conduction in the cross-plane direction. Both experimental work 51 on Si(crystal)-SiGe(disordered alloy) nanowires and numerical simulations 40 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 [37][38][39] , or interfacial species mixing 40,41 destroy the coherent features of vibrational excitations present in the investigated superlattices. 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 κ CP just corresponds to the disordered alloy limit, i.e., κ CP cannot be reduced below the glass limit.
Finally it is worth to mention that the thermal conductivity tensor becomes increasingly anisotropic at larger W due to the increase of κ IP , showing a behaviour different than that observed in S1 where the anisotropy reaches the maximum at the crossover point  ⁎ W 20. Fig. 7 we show the W-dependences of κ CP and κ IP for the case where the energy scale associated with particles interactions across the interfaces (ε AB ) are lowered compared to those intra-layers, with ε AB = 0.5 and 0.1 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. Taking the one component system as reference is justified on the basis of earlier work 8 . There, it was shown that when a system of soft spheres is frozen in a disordered state, slightly modifying even a substantial fraction of the interactions does not have any appreciable influence on the vibrational properties, including thermal transport. As a consequence, in the present case we do not expect any relevant modification to κ originating from a limited fraction (scaling with the surface of the interfaces) of modified interactions in the glass sample. The in-plane value κ IP is almost independent of W, and is very close to the value pertaining to the crystal. In contrast, κ CP decreases monotonically by increasing W, and especially in the weaker case ε AB = 0.1, the observed reduction of κ CP is dramatic. At W = 10, κ CP equals the value obtained for the glassy sample, and it is almost two orders of magnitude lower than this value at W = 20. This extremely low κ CP is consistent with previous experimental work 19 .

S3. Superlattice composed by identical crystalline layers separated by weakly interacting interfaces. In
Some insight about the origin of this observation comes from the data shown in Fig. 8, where we display the average cross-plane distance δz between adjacent crystalline planes (monolayers), normalized to the value in the perfect lattice, a/2. For ε AB = 0.5 and W = 4, the system keeps the perfect lattice structure, with δz ≡ a/2 for all monolayers. In contrast, as ε AB decreases and for a large value W = 20, δz becomes substantially larger than a/2 at the interfaces, which therefore assumes a local density lower than the average. At the same time, slightly reduced δz 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 W = 20, the values of δz at the interfaces located at W/2 = 10 and W = 20 are different, with a large discrepancy for ε AB = 0.1. 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 δz at different interfaces.  The behaviour of κ CP can be further elucidated by inspection of the main features of the vibrational spectrum. In Fig. 9 we plot the g(ω) of superlattice S3, together with the vibrational amplitudes E A k and E B k . The transverse and longitudinal phonon branches are observed in g(ω) (top panels) for all cases, similar to the homogeneous bulk crystal. For ε AB = 0.5, W = 4, g(ω) shows an excess of lower-ω modes compared to those present in the one-component crystal, simply due to the weakened interactions at the interfaces. As ε AB decreases and W increases, g(ω) 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 the cases of ε AB = 0.5, W = 20 and ε AB = 0.1, W = 20, that the distance between monolayers far from the interfaces becomes smaller than a/2. The consequent larger mass density makes higher the frequency of phonon modes of given wavelength, leading to the shift of g(ω) towards higher frequencies.
We now focus on the vibrational amplitudes, E A k and E B k (Fig. 9, bottom panels). In the cases with ε AB = 0.5 and W = 4 and 20, the particles in the two layers A and B 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 κ CP . The situation changes drastically in the case ε AB = 0.1 and W = 20, where the ultra-low value of κ CP can be reached. Except for the low-ω modes, E A k and E B k are symmetrically randomly distributed around the average values ≡ .
( ) , indicating that particles in layers A and B 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 A and B, leading to extremely low value of κ CP , while keeping a high κ IP .
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.

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 κ CP below the corresponding glass limit, while keeping the in-plane contribution κ IP close to the pure crystalline value. More specifically, in the case of mass mismatch (S1), propagation of phonons with high frequencies (ω ω > B max ) is almost completely suppressed, whereas a fraction of low-frequency phonons (ω ω < B max ) are still able to propagate across the interfaces, contributing to κ CP (Fig. 4(d-f)). Also, the minimum in thermal conductivity as a function of the repetition period W (Fig. 2) corresponds to a maximum in the vibrational separation between the layers of type A and B. 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 a very wide range of frequencies (see Fig. 9(c)), which results into extremely low values of κ CP , 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 WSe 2 sheets in Ref. 19. Interfaces stiffness modification by controlling pressure 44,45 or chemical bonding 46 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 κ CP . 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 W (see Fig. 5). Also, as suggested in previous works, disorder in the interfacial roughness [37][38][39] or interfacial mixing 40,41 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 49,52 . 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. 53. A possibility to realize this task is embedding in the targeted material objects featuring larger typical sizes, including nano-particles 54,55 or nano(quantum)-dots 20,56 . Based on this strategy, very low thermal conductivity was achieved experimentally in a Si-Ge quantum-dot superlattice 20 , 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 57 . Here, values of thermal conductivity 10 4 times smaller than that of bulk Si were reached in Si phononic crystals with spherical pores.
An other important point must be underlined. In the present work we have investigated the different systems at low temperature, to both focus on plainly structural effects and keep contact with the well controlled harmonic limit. Anharmonicities, however, are expected to play a crucial role at temperatures higher than the Debye temperature, T D , which are those relevant for technological applications. This is a crucial aspect to be explored in extended future work.
In conclusion, we note that the three superlattice structures studied in the present work show totally different W-dependences of cross and in-plane thermal conductivities. Our results therefore not only contribute to a deeper understanding 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.

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 V = L 3 (L 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 where r is the distance between those two particles, and σ ij and ε ij are the interparticle diameter and interaction energy scale, respectively. The potential is cut-off and shifted at r c = 2.5σ ij . Particle i has mass m i , and we have used σ, ε/k B (k B is the Boltzmann constant), and m as units of length, temperature, and mass. As a reference, for Argon σ = 3.4 Å, ε/k B = 120 K, and m = 39.96 a.m.u. We considered the number density ρ = / = . N V 1 015, corresponding to a lattice constant ρ = ( / ) = . / a 4 158 1 3 . We prepared three superlattices, composed of intercalated FCC lattice layers, A and B, both of thickness W/2, as schematically illustrated in Fig. 1. The first superlattice (S1) consists of two crystalline layers formed by point particles with different masses, m A and m B . We have considered mass ratios m B /m A > 1, while keeping a constant average mass (m A + m B )/2 = 1. As an example, the case m B /m A = 4 corresponds to m A = 0.4 and m B = 1.6. We have dubbed A and B as the light and heavy layers, respectively. Note that a mass ratio of m B /m A = 2.5 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 σ ij = ε ij = 1. The energy scale of interactions between particles pertaining to different layers are, however, reduced to ε ij = ε AB < 1.

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 58 . The systems were first equilibrated at relatively low temperature T = 10 −2 by MD simulation in the NVT-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-T regime 9 . 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 59,60 , and we have therefore chosen to stay within a fully classical approach.
Following the equilibration stage, we performed the production runs in the NVE-ensemble. The Green-Kubo formulation 47,48 was next applied to calculate the thermal conductivities, in the cross-plane and in-plane directions, respectively: Here, J x,y and J z are the heat currents in the in-plane (x,y) and cross-plane (z) directions, and 〈 〉 denotes the ensemble average. Landry et al. 48 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. In the bulk glasses and disordered alloys, κ κ  CP IP , i.e., heat conduction is isotropic, whereas in the superlattices, they are expected to assume different values 16,17 .
More in details, the equations of motions were integrated numerically with a time step, δt = 5 × 10 −3 , for a total run-time t run = N run δt = 10 5 (N run = 2 × 10 7 steps) for S1 and t run = 10 4 (N run = 2 × 10 6 steps) for both S2 and S3. The systems snapshots extracted from the trajectory have been used to calculate the correlation functions, , , , , J t J 0 x y z xy z of Eq. (4), which have been subsequently integrated numerically over a finite time window Δ t. We have observed a clear convergence of the integrals for κ CP,IP , to the limits of Eq. (4) for Δ t = 10 4 for S1, and Δ t = 10 3 for both S2 and S3. For each considered superlattice, and all repetition periods W, we have performed 10 independent calculations starting from different initial system configurations. These has allowed us to obtain 10 independent sample values for κ CP and 2 × 10 = 20 (we have considered both the x and y components) for κ IP . These values were used to calculate averages and sample-to-sample fluctuations (standard deviations), which are shown as error bars in Figs. 2,5 and 7. More details on these calculations are reported in the Supplementary Materials. Analysis of finite system size effect on thermal conductivity. In the Green-Kubo calculations of thermal conductivities, one must be attentive to finite system size effects 47,48 . Indeed, long-wavelength phonons with λ > L are excluded from the simulation box due to the finite value L 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 Scientific RepoRts | 5:14116 | DOi: 10.1038/srep14116 anharmonic coupling (scattering) processes 47 . We note that the considered T = 10 −2 is low enough to substantially reduce anharmonic effects, but anharmonic couplings are still present. We can take care of finite size effects by increasing L to values where κ CP and κ IP become L-independent. For the glass and disordered alloy thermal conductivities, we have confirmed that a system size L = 10a (N = 4,000) is sufficiently large to obtain correct values of κ κ  CP IP , without any size effect 8 . In the superlattice cases, the appropriate L depends on the considered structure and the periodic repetition length W 48 . More in details, we paid particular attention to the number P of repetitions, defined from L = PW, 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 L = 10a (20 monolayers, N = 4,000) to 24a (48 monolayers, N = 55,296). We provide details of our analysis of the finite-size effects in the Supplementary Materials. In Figs. 2,5 and 7, we plot the values obtained by using the largest systems (the exact system size depends on both the superlattice type and W), which show the smallest finite-size effects.
For the S1 superlattice, we confirmed that the required number P of repetitions becomes larger for smaller W 48 : one period (L = W) only is adequate for W ≥ 20, whereas four periods or more (L ≥ 4W) are required for W ≤ 8 (see the Supplementary Materials). We have therefore employed four pattern repetitions (L = 4W) for 10 ≤ W ≤ 12 and two (L = 2W) for 14 ≤ W ≤ 18. This behaviour is simple to rationalize by inspecting the data in Fig. 2, where the crossover between incoherent and coherent phonon transport occurs around  ⁎ W 20. In the coherent regime W < 20, 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 P are needed (even P = 1) in the incoherent regime W > 20, 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 have confirmed that P = 1 or 2 (L = W or 2W) are sufficient for W ≥ 20, while two or more repetitions (L ≥ 2W) are appropriate for W < 20, for both S2 and S3 (see the Supplementary Materials). We can understand this behaviour by noticing that phonon transport 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 λ > L, play very little role in phonon transport and finite system size effects are consequently negligible.

Normal modes analysis.
We have characterized the superlattice vibrational states (superlattice phonons) by performing a standard normal-mode analysis 8,49 . We have diagonalized the dynamical (Hessian) matrix calculated at local minima of the potential energy landscape, and obtained eigenvalues  ), particles in layer A (B) contribute more to mode k than those in layer B (A). 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 T → 0 which, we believe, is an appropriate approximation for our case T = 10 −2 , where anharmonicities are weak.