# Beating the amorphous limit in thermal conductivity by superlattices design

Hideyuki Mizuno, Stefano Mossa, Jean-Louis Barrat

► To cite this version:

Hideyuki Mizuno, Stefano Mossa, Jean-Louis Barrat. Beating the amorphous limit in thermal conductivity by superlattices design. Scientific Reports, 2015, 5, pp.14116. 10.1038/srep14116. cea-01734605

HAL Id: cea-01734605

https://cea.hal.science/cea-01734605

Submitted on 4 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_1.png" /></figure>

OPEN

# Beating the amorphous limit in thermal conductivity by superlattices design

Received: 05 May 2015

Accepted: 03 August 2015

Published: 16 September 2015

Hideyuki Mizuno<sup>1,2,†</sup>, Stefano Mossa<sup>3,4,†</sup> & Jean-Louis Barrat<sup>1,2,6</sup>

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.

Materials with low thermal conductivity, $\kappa$, are employed in many modern technologies, such as thermal management in electronic devices or thermoelectric energy conversion<sup>1-3</sup>. In general, low values of $\kappa$ are observed in disordered solids<sup>4</sup>, including topologically disordered systems and crystalline solids with size or mass disorder<sup>5-8</sup>. We refer to the former as glasses, to the latter as disordered alloys. While in glasses disorder originates both from molecules size or mass heterogeneity and the topologically amorphous structure, in disordered alloys molecules characterized by heterogeneous attributes are located at ordered lattice sites. The low $\kappa$ in disordered solids can be rationalized by considering the phenomenological kinetic theory expression<sup>9</sup> $\kappa = (1/3)Cv^2\tau$, which relates average velocity, $v$, and lifetime, $\tau$, (and therefore the mean free path $\ell = v\tau$) of phonons to $\kappa$ (C 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<sup>10</sup>) reduces $\tau$ (or $\ell$) and, as a result, $\kappa$.

In early experimental investigations<sup>5,6</sup>, Cahill et al. have studied the disordered alloys, e.g., $(KBr)_{1-x}(KCN)_x$, $(NaCl)_{1-x}(NaCN)_x$, and shown that $\kappa$ can be reduced to the glass value by controlling the relative composition $x$. In our work<sup>8</sup> we in turn demonstrated that, in size-disordered crystal (alloy), $\kappa$ progressively decreases with increasing size mismatch, eventually converging to the corresponding glass value. When this limit is reached, $\tau$ is comparable to the time scale of thermal vibrations ($\ell$ to the particle size), i.e., to the minimum time (length) scale<sup>8</sup>. Heat propagation can therefore be described as a random walk of vibrational energies<sup>5,6</sup>, or in terms of non-propagating delocalized modes, the diffusions<sup>7</sup>. For this reason, the value in the glass is generally considered as a lower bound for $\kappa$ of materials with homogeneous chemical composition<sup>5,6</sup>.

<footnote n="1"><sup>1</sup>Univ. Grenoble Alpes, LIPHY, F-38000 Grenoble, France. <sup>2</sup>CNRS, LIPHY, F-38000 Grenoble, France. <sup>3</sup>Univ. Grenoble Alpes, INAC-SPRAM, F-38000 Grenoble, France. <sup>4</sup>CNRS, INAC-SPRAM, F-38000 Grenoble, France. <sup>5</sup>CEA, INAC-SPRAM, F-38000 Grenoble, France. <sup>6</sup>Institut Laue-Langevin - 6 rue Jules Horowitz, BP 156, 38042 Grenoble, France. <sup>†</sup>Present address: Institut für Materialphysik im Weltraum, Deutsches Zentrum für Luft- und Raumfahrt (DLR), 51170 Köln, Germany. Correspondence and requests for materials should be addressed to S.M. (email: stefano.mossa@cea.fr)</footnote>
SCIENTIFIC REPORTS

<figure><img src="image_2.png" /></figure>

OPEN

# Beating the amorphous limit in thermal conductivity by superlattices design

Received: 05 May 2015

Accepted: 03 August 2015

Published: 16 September 2015

Hideyuki Mizuno<sup>1,2,†</sup>, Stefano Mossa<sup>3,4,†</sup> & Jean-Louis Barrat<sup>1,2,6</sup>

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.

Materials with low thermal conductivity, $\kappa$, are employed in many modern technologies, such as thermal management in electronic devices or thermoelectric energy conversion<sup>1-3</sup>. In general, low values of $\kappa$ are observed in disordered solids<sup>4</sup>, including topologically disordered systems and crystalline solids with size or mass disorder<sup>5-8</sup>. We refer to the former as glasses, to the latter as disordered alloys. While in glasses disorder originates both from molecules size or mass heterogeneity and the topologically amorphous structure, in disordered alloys molecules characterized by heterogeneous attributes are located at ordered lattice sites. The low $\kappa$ in disordered solids can be rationalized by considering the phenomenological kinetic theory expression<sup>9</sup> $\kappa = (1/3)Cv^2\tau$, which relates average velocity, $v$, and lifetime, $\tau$, (and therefore the mean free path $\ell = v\tau$) of phonons to $\kappa$ (C 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<sup>10</sup>) reduces $\tau$ (or $\ell$) and, as a result, $\kappa$.

In early experimental investigations<sup>5,6</sup>, Cahill et al. have studied the disordered alloys, e.g., $(KBr)_{1-x}(KCN)_x$, $(NaCl)_{1-x}(NaCN)_x$, and shown that $\kappa$ can be reduced to the glass value by controlling the relative composition $x$. In our work<sup>8</sup> we in turn demonstrated that, in size-disordered crystal (alloy), $\kappa$ progressively decreases with increasing size mismatch, eventually converging to the corresponding glass value. When this limit is reached, $\tau$ is comparable to the time scale of thermal vibrations ($\ell$ to the particle size), i.e., to the minimum time (length) scale<sup>8</sup>. Heat propagation can therefore be described as a random walk of vibrational energies<sup>5,6</sup>, or in terms of non-propagating delocalized modes, the diffusions<sup>7</sup>. For this reason, the value in the glass is generally considered as a lower bound for $\kappa$ of materials with homogeneous chemical composition<sup>5,6</sup>.

<footnote n="1"><sup>1</sup>Univ. Grenoble Alpes, LIPHY, F-38000 Grenoble, France. <sup>2</sup>CNRS, LIPHY, F-38000 Grenoble, France. <sup>3</sup>Univ. Grenoble Alpes, INAC-SPRAM, F-38000 Grenoble, France. <sup>4</sup>CNRS, INAC-SPRAM, F-38000 Grenoble, France. <sup>5</sup>CEA, INAC-SPRAM, F-38000 Grenoble, France. <sup>6</sup>Institut Laue-Langevin - 6 rue Jules Horowitz, BP 156, 38042 Grenoble, France. <sup>†</sup>Present address: Institut für Materialphysik im Weltraum, Deutsches Zentrum für Luft- und Raumfahrt (DLR), 51170 Köln, Germany. Correspondence and requests for materials should be addressed to S.M. (email: stefano.mossa@cea.fr)</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic illustration of the considered superlattice structures. 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.</figcaption></figure>

A crucial issue⁴ is whether thermal conductivity can be lowered below the glass limit through nano-scale phononic design³,¹¹. This possibility would allow to devise (meta-)materials which are excellent thermal insulators while preserving good electronic properties, as needed in many applications¹⁻³. The most popular design to reach this goal is that of a lamellar superlattice¹²⁻¹⁵, often composed of two chemically different intercalated layers, e.g., Si-Ge¹²,¹³ or GaAs-AlAs¹⁴,¹⁵ (see also Fig. 1). In a superlattice, the thermal conductivity tensor is anisotropic, with the cross-plane component, $\kappa_{CP}$, usually lower than the in-plane value, $\kappa_{IP}$¹⁶,¹⁷.

In recent experiments¹⁸⁻²⁰, ultra-low values of $\kappa_{CP}$ suggested to be smaller than the amorphous limit, were measured. In particular, Costescu *et al.*¹⁸ demonstrated that the presence of a high-density of interfaces decreases $\kappa_{CP}$ of W-Al₂O₃ nanolaminates, below that of the amorphous Al₂O₃. An experiment by Chiritescu *et al.*¹⁹ achieved ultra-low thermal conductivity in layered WSe₂ crystals, by disordering the crystalline WSe₂ sheets. Finally, Pernot *et al.*²⁰ also observed very low values of $\kappa_{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 $\kappa$ 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 $\kappa$, it is therefore still not completely clear whether superlattice structuration alone can lower $\kappa$ below the amorphous limit. Note, for instance, that very recent numerical work²¹, has shown that a superlattice composed by layers with randomized thicknesses can indeed show a $\kappa$ below the value pertaining to the disordered-alloy with the same composition. This limit, however, is generally higher than that in the corresponding glass⁵⁻⁸, 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 $\kappa$ 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⁸, 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, $l$, determines the coherent or incoherent character of phonon transport, as described in²²⁻²⁴ and demonstrated by numerical simulations²⁵⁻²⁷ and recent experiments²⁸.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic illustration of the considered superlattice structures. 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.</figcaption></figure>

A crucial issue⁴ is whether thermal conductivity can be lowered below the glass limit through nano-scale phononic design³,¹¹. This possibility would allow to devise (meta-)materials which are excellent thermal insulators while preserving good electronic properties, as needed in many applications¹⁻³. The most popular design to reach this goal is that of a lamellar superlattice¹²⁻¹⁵, often composed of two chemically different intercalated layers, e.g., Si-Ge¹²,¹³ or GaAs-AlAs¹⁴,¹⁵ (see also Fig. 1). In a superlattice, the thermal conductivity tensor is anisotropic, with the cross-plane component, $\kappa_{CP}$, usually lower than the in-plane value, $\kappa_{IP}$¹⁶,¹⁷.

In recent experiments¹⁸⁻²⁰, ultra-low values of $\kappa_{CP}$ suggested to be smaller than the amorphous limit, were measured. In particular, Costescu et al.¹⁸ demonstrated that the presence of a high-density of interfaces decreases $\kappa_{CP}$ of W-Al₂O₃ nanolaminates, below that of the amorphous Al₂O₃. An experiment by Chiritescu et al.¹⁹ achieved ultra-low thermal conductivity in layered WSe₂ crystals, by disordering the crystalline WSe₂ sheets. Finally, Pernot et al.²⁰ also observed very low values of $\kappa_{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 $\kappa$ 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 $\kappa$, it is therefore still not completely clear whether superlattice structuration alone can lower $\kappa$ below the amorphous limit. Note, for instance, that very recent numerical work²¹, has shown that a superlattice composed by layers with randomized thicknesses can indeed show a $\kappa$ below the value pertaining to the disordered-alloy with the same composition. This limit, however, is generally higher than that in the corresponding glass⁵⁻⁸, 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 $\kappa$ 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⁸, 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, $l$, determines the coherent or incoherent character of phonon transport, as described in²²⁻²⁴ and demonstrated by numerical simulations²⁵⁻²⁷ and recent experiments²⁸.
<table><thead><tr><th>System</th><th>Control Parameter</th><th>κ<sub>A</sub></th><th>κ<sub>B</sub></th><th>κ<sub>CP</sub><sup>∞</sup></th><th>κ<sub>IP</sub><sup>∞</sup></th><th>R</th><th>l<sub>K</sub></th><th>κ<sub>glass</sub></th><th>κ<sub>alloy</sub></th><th>Fig.</th></tr></thead><tbody><tr><td rowspan="3">(S1) Mass mismatch</td><td>m<sub>B</sub>/m<sub>A</sub> = 2</td><td>488.6</td><td>335.4</td><td>397.8</td><td>412.0</td><td>0.5</td><td>398</td><td>5.7</td><td>20.4</td><td>Fig. 2(a)</td></tr><tr><td>m<sub>B</sub>/m<sub>A</sub> = 4</td><td>625.9</td><td>306.8</td><td>411.8</td><td>466.3</td><td>1.9</td><td>1564</td><td>4.2</td><td>9.9</td><td>Fig. 2(b)</td></tr><tr><td>m<sub>B</sub>/m<sub>A</sub> = 8</td><td>843.8</td><td>291.1</td><td>432.8</td><td>567.5</td><td>—</td><td>—</td><td>3.3</td><td>7.7</td><td>Fig. 2(c)</td></tr><tr><td rowspan="3">(S2) Crystal-alloy</td><td>m<sub>B</sub>/m<sub>B<sub>1</sub></sub> = 2</td><td>381.6</td><td>20.4</td><td>38.7</td><td>201.0</td><td>—</td><td>—</td><td>5.7</td><td>33.2</td><td>Fig. 5(a)</td></tr><tr><td>m<sub>B</sub>/m<sub>B<sub>1</sub></sub> = 4</td><td>381.6</td><td>9.9</td><td>19.3</td><td>195.8</td><td>—</td><td>—</td><td>4.5</td><td>14.3</td><td>Fig. 5(b)</td></tr><tr><td>m<sub>B</sub>/m<sub>B<sub>1</sub></sub> = 8</td><td>381.6</td><td>7.7</td><td>15.1</td><td>194.6</td><td>—</td><td>—</td><td>4.0</td><td>8.2</td><td>Fig. 5(c)</td></tr><tr><td rowspan="2">(S3) Weakened interface</td><td>ε<sub>AB</sub> = 0.5</td><td>587.3</td><td>587.3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>10.6</td><td>—</td><td>Fig. 7(a)</td></tr><tr><td>ε<sub>AB</sub> = 0.1</td><td>587.3</td><td>587.3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>10.6</td><td>—</td><td>Fig. 7(b)</td></tr></tbody></table>

**Table 1.** Details of the three superlattice systems investigated in this work. They are based on the FCC-crystal 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 mass-disordered 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 $\epsilon_{AB}$ of the interactions across the interfaces in S3. Number density and temperature were fixed to the values $\hat{\rho}=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.

For $W > \ell$, 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-33}$). The thermal conductivity in the cross-plane direction can be written as$^{22}$

$$ \kappa_{CP} = \frac{2}{\kappa_A^{-1} + \kappa_B^{-1} + 4RW^{-1}} = \kappa_{CP}^{\infty} \left( \frac{1}{1 + \ell_K W^{-1}} \right), \quad (1) $$

where

$$ \kappa_{CP}^{\infty} = \lim_{W \to \infty} \kappa_{CP} = \frac{2\kappa_A \kappa_B}{\kappa_A + \kappa_B}. \quad (2) $$

Here, $\kappa_A$ and $\kappa_B$ are the thermal conductivities of materials A and B, and $\ell_K=2R\kappa_{CP}^{\infty}$ 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}$. For $W < \ell_K$ ($W > \ell_K$), the interfacial resistance is relatively large (small) compared to the internal resistance. Both $\kappa_{CP}$ and $\kappa_{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 \to \infty$, where the interfacial resistance can be neglected, $\kappa_{CP}$ and $\kappa_{IP}$ have the upper bounds $\kappa_{CP}^{\infty}$ and $\kappa_{IP}^{\infty} = (\kappa_A + \kappa_B)/2$, respectively.

When $W < \ell$, phonon transport is coherent$^{22-28}$, and the wave nature of phonons cannot be neglected. In this regime, $\kappa_{CP}$ decreases with increasing W, in contrast with the incoherent case. The reduction of $\kappa_{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 $\kappa_{CP}$. Mini-umklapp processes$^{36}$, occurring at the mini-Brillouin zone, also contribute to the reduction of $\kappa_{CP}$. At the crossover length $W \sim \ell$, between the incoherent and the coherent transport regimes, $\kappa_{CP}$ assumes a minimum value when plotted against $W^{22-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-46}$. It has been reported that interfacial roughness$^{37-39}$ or mixing$^{40,41}$ reduce both $\kappa_{CP}$ and $\kappa_{IP}$, and can even suppress the coherent nature of phonons, with $\kappa_{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><caption>**Table 1.** Details of the three superlattice systems investigated in this work. They are based on the FCC-crystal 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 mass-disordered 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 $\epsilon_{AB}$ of the interactions across the interfaces in S3. Number density and temperature were fixed to the values $\hat{\rho}=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.</caption><thead><tr><th>System</th><th>Control Parameter</th><th>κ<sub>A</sub></th><th>κ<sub>B</sub></th><th>κ<sub>CP</sub><sup>∞</sup></th><th>κ<sub>IP</sub><sup>∞</sup></th><th>R</th><th>ℓ<sub>K</sub></th><th>κ<sub>glass</sub></th><th>κ<sub>alloy</sub></th><th>Fig.</th></tr></thead><tbody><tr><td rowspan="3">(S1) Mass mismatch</td><td>m<sub>B</sub>/m<sub>A</sub> = 2</td><td>488.6</td><td>335.4</td><td>397.8</td><td>412.0</td><td>0.5</td><td>398</td><td>5.7</td><td>20.4</td><td>Fig. 2(a)</td></tr><tr><td>m<sub>B</sub>/m<sub>A</sub> = 4</td><td>625.9</td><td>306.8</td><td>411.8</td><td>466.3</td><td>1.9</td><td>1564</td><td>4.2</td><td>9.9</td><td>Fig. 2(b)</td></tr><tr><td>m<sub>B</sub>/m<sub>A</sub> = 8</td><td>843.8</td><td>291.1</td><td>432.8</td><td>567.5</td><td>—</td><td>—</td><td>3.3</td><td>7.7</td><td>Fig. 2(c)</td></tr><tr><td rowspan="3">(S2) Crystal-alloy</td><td>m<sub>B</sub>/m<sub>B<sub>1</sub></sub> = 2</td><td>381.6</td><td>20.4</td><td>38.7</td><td>201.0</td><td>—</td><td>—</td><td>5.7</td><td>33.2</td><td>Fig. 5(a)</td></tr><tr><td>m<sub>B</sub>/m<sub>B<sub>1</sub></sub> = 4</td><td>381.6</td><td>9.9</td><td>19.3</td><td>195.8</td><td>—</td><td>—</td><td>4.5</td><td>14.3</td><td>Fig. 5(b)</td></tr><tr><td>m<sub>B</sub>/m<sub>B<sub>1</sub></sub> = 8</td><td>381.6</td><td>7.7</td><td>15.1</td><td>194.6</td><td>—</td><td>—</td><td>4.0</td><td>8.2</td><td>Fig. 5(c)</td></tr><tr><td rowspan="2">(S3) Weakened interface</td><td>ε<sub>AB</sub> = 0.5</td><td>587.3</td><td>587.3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>10.6</td><td>—</td><td>Fig. 7(a)</td></tr><tr><td>ε<sub>AB</sub> = 0.1</td><td>587.3</td><td>587.3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>10.6</td><td>—</td><td>Fig. 7(b)</td></tr></tbody></table>

For $W > \ell$, 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–33}$). The thermal conductivity in the cross-plane direction can be written as$^{22}$

$$ \kappa_{CP} = \frac{2}{\kappa_A^{-1} + \kappa_B^{-1} + 4RW^{-1}} = \kappa_{CP}^{\infty} \left( \frac{1}{1 + \ell_K W^{-1}} \right), \quad (1) $$

where

$$ \kappa_{CP}^{\infty} = \lim_{W \to \infty} \kappa_{CP} = \frac{2\kappa_A \kappa_B}{\kappa_A + \kappa_B}. \quad (2) $$

Here, $\kappa_A$ and $\kappa_B$ are the thermal conductivities of materials A and B, and $\ell_K=2R\kappa_{CP}^{\infty}$ 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}$. For $W < \ell_K$ ($W > \ell_K$), the interfacial resistance is relatively large (small) compared to the internal resistance. Both $\kappa_{CP}$ and $\kappa_{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 \to \infty$, where the interfacial resistance can be neglected, $\kappa_{CP}$ and $\kappa_{IP}$ have the upper bounds $\kappa_{CP}^{\infty}$ and $\kappa_{IP}^{\infty} = (\kappa_A + \kappa_B)/2$, respectively.

When $W < \ell$, phonon transport is coherent$^{22–28}$, and the wave nature of phonons cannot be neglected. In this regime, $\kappa_{CP}$ decreases with increasing W, in contrast with the incoherent case. The reduction of $\kappa_{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 $\kappa_{CP}$. Mini-umklapp processes$^{36}$, occurring at the mini-Brillouin zone, also contribute to the reduction of $\kappa_{CP}$. At the crossover length $W \sim \ell$, between the incoherent and the coherent transport regimes, $\kappa_{CP}$ assumes a minimum value when plotted against $W^{22–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–46}$. It has been reported that interfacial roughness$^{37–39}$ or mixing$^{40,41}$ reduce both $\kappa_{CP}$ and $\kappa_{IP}$, and can even suppress the coherent nature of phonons, with $\kappa_{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.
## Results

In Table 1, we present the details of the three superlattice systems studied in this work, with values of the important quantities: $\kappa_A$ and $\kappa_B$ are the thermal conductivities of layers A and B, respectively; $\kappa_{CP}^\infty$ and $\kappa_{IP}^\infty$ are the cross- and in-plane diffuse limits of $\kappa_{CP}$ and $\kappa_{IP}$; R is the interfacial resistance, $\ell_K$ the Kapitza length; $\kappa_{glass}$ and $\kappa_{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 $\hat{\rho} = 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, $\kappa_{CP}$ and $\kappa_{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 $\kappa_{CP}$ to $\kappa_{CP}^\infty$, $\kappa_{glass}$, $\kappa_{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 $\kappa_{IP}$ to $\kappa_{IP}^\infty$. 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, $\kappa_{CP}$ and $\kappa_{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, $\langle m \rangle = (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 $\langle m \rangle^{-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 $\kappa_{CP}^\infty$ and $\kappa_{IP}^\infty$ 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 $\sqrt{\overline{m_A}}\kappa_A = \sqrt{\overline{m_B}}\kappa_B$ holds for the pure materials. In the studied W-range, $W=2$ to 40 (monolayers), the in-plane value $\kappa_{IP}$ shows a very weak dependence on $W$, as was observed for superlattices with perfect interfaces in Refs 26,39. The value of $\kappa_{IP}$ is close to, although lower than, $\kappa_{IP}^\infty$, indicating that slight in-plane phonon scattering at the interface is still active.

More interestingly, as $W$ increases, the cross-plane value $\kappa_{CP}$ decreases steeply, reaches a minimum value at $W^* \simeq 20$, and next increases mildly at larger $W$. This $W$-dependence is consistent with previous predictions$^{22-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 $\kappa_{CP}$ (dashed line in Fig. 2) we can extract the values of the interfacial resistance, $R$, and the Kapitza length, $\ell_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$, $\kappa_{CP}$ is still orders of magnitude lower than $\kappa_{CP}^\infty$, indicating that the interfacial resistance $R$ results in a strong reduction of $\kappa_{CP}$ in this range of $W$. Equivalently, the Kapitza length $\ell_K$ is significantly larger than the maximum period $W=40$. Also, if $\ell_K \gg W$, Eq. (1) can be approximated as $\kappa_{CP} \simeq (\kappa_{CP}^\infty/\ell_K)W$, i. e., $\kappa_{CP}$ increases linearly with $W$, with a slope, $\kappa_{CP}^\infty/\ell_K$, which decreases by increasing the ratio $m_B/m_A$. For instance, it is $\simeq 1.0$ for $m_B/m_A=2$ and $\simeq 0.26$ for $m_B/m_A=4$. For $m_B/m_A=8$, the slope is so small that $\kappa_{CP}$ appears almost flat when plotted against $W$, for $20 \le W \le 40$ (Fig. 2(c)).

The data shown in Fig. 2 demonstrate that $\kappa_{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 $\kappa_{CP}$. It is also worth noting that the thermal conductivity tensor is very strongly anisotropic in this case, with $\kappa_{CP} \ll \kappa_{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(\omega)$, for $m_B/m_A=4$ and $W=2$ to 80. $g_A(\omega)$ and $g_B(\omega)$ 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 $g_A(\sqrt{\overline{m_A}}\omega)/\sqrt{\overline{m_A}} = g_B(\sqrt{\overline{m_B}}\omega)/\sqrt{\overline{m_B}}$. At small $W=2$, $g(\omega)$ 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(\omega)$ generates features increasingly similar to those identifying $g_A(\omega)$ and $g_B(\omega)$, separately. In particular, in the low-$\omega$ region $g(\omega)$ follows $g_B(\omega)$ (the heavy crystal B), whereas $g_A(\omega)$ (the light crystal A) controls $g(\omega)$ in the high-$\omega$ region. This result indicates that different parts of the vibrational spectrum are active in the two layers, with high (low)-$\omega$ modes preferentially excited in the light (heavy)
## Results

In Table 1, we present the details of the three superlattice systems studied in this work, with values of the important quantities: $\kappa_A$ and $\kappa_B$ are the thermal conductivities of layers A and B, respectively; $\kappa_{CP}^\infty$ and $\kappa_{IP}^\infty$ are the cross- and in-plane diffuse limits of $\kappa_{CP}$ and $\kappa_{IP}$; R is the interfacial resistance, $\ell_K$ the Kapitza length; $\kappa_{glass}$ and $\kappa_{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 $\hat{\rho} = 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, $\kappa_{CP}$ and $\kappa_{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 $\kappa_{CP}$ to $\kappa_{CP}^\infty$, $\kappa_{glass}$, $\kappa_{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 $\kappa_{IP}$ to $\kappa_{IP}^\infty$. 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, $\kappa_{CP}$ and $\kappa_{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, $\langle m \rangle = (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 $\langle m \rangle^{-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 $\kappa_{CP}^\infty$ and $\kappa_{IP}^\infty$ 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 $\sqrt{\overline{m_A}}\kappa_A = \sqrt{\overline{m_B}}\kappa_B$ holds for the pure materials. In the studied W-range, $W=2$ to 40 (monolayers), the in-plane value $\kappa_{IP}$ shows a very weak dependence on $W$, as was observed for superlattices with perfect interfaces in Refs 26,39. The value of $\kappa_{IP}$ is close to, although lower than, $\kappa_{IP}^\infty$, indicating that slight in-plane phonon scattering at the interface is still active.

More interestingly, as $W$ increases, the cross-plane value $\kappa_{CP}$ decreases steeply, reaches a minimum value at $W^* \simeq 20$, and next increases mildly at larger $W$. This $W$-dependence is consistent with previous predictions$^{22-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 $\kappa_{CP}$ (dashed line in Fig. 2) we can extract the values of the interfacial resistance, $R$, and the Kapitza length, $\ell_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$, $\kappa_{CP}$ is still orders of magnitude lower than $\kappa_{CP}^\infty$, indicating that the interfacial resistance $R$ results in a strong reduction of $\kappa_{CP}$ in this range of $W$. Equivalently, the Kapitza length $\ell_K$ is significantly larger than the maximum period $W=40$. Also, if $\ell_K \gg W$, Eq. (1) can be approximated as $\kappa_{CP} \simeq (\kappa_{CP}^\infty/\ell_K) W$, i. e., $\kappa_{CP}$ increases linearly with $W$, with a slope, $\kappa_{CP}^\infty/\ell_K$, which decreases by increasing the ratio $m_B/m_A$. For instance, it is $\simeq 1.0$ for $m_B/m_A=2$ and $\simeq 0.26$ for $m_B/m_A=4$. For $m_B/m_A=8$, the slope is so small that $\kappa_{CP}$ appears almost flat when plotted against $W$, for $20 \le W \le 40$ (Fig. 2(c)).

The data shown in Fig. 2 demonstrate that $\kappa_{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 $\kappa_{CP}$. It is also worth noting that the thermal conductivity tensor is very strongly anisotropic in this case, with $\kappa_{CP} \ll \kappa_{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(\omega)$, for $m_B/m_A=4$ and $W=2$ to 80. $g_A(\omega)$ and $g_B(\omega)$ 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 $g_A(\sqrt{\overline{m_A}}\omega)/\sqrt{\overline{m_A}} = g_B(\sqrt{\overline{m_B}}\omega)/\sqrt{\overline{m_B}}$. At small $W=2$, $g(\omega)$ 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(\omega)$ generates features increasingly similar to those identifying $g_A(\omega)$ and $g_B(\omega)$, separately. In particular, in the low-$\omega$ region $g(\omega)$ follows $g_B(\omega)$ (the heavy crystal B), whereas $g_A(\omega)$ (the light crystal A) controls $g(\omega)$ in the high-$\omega$ region. This result indicates that different parts of the vibrational spectrum are active in the two layers, with high (low)-$\omega$ modes preferentially excited in the light (heavy)
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Thermal conductivity in superlattice SI composed of two intercalated crystalline layers with different masses. The cross-plane, $\kappa_{CP}$, and in-plane, $\kappa_{IP}$, components of thermal conductivity are plotted as functions of the repetition period $W$. The ratio $m_B/m_A$ of the masses in layers A and B is 2 in panel (a), 4 in (b), and 8 in (c). The values $\kappa_{CP}^\infty$ and $\kappa_{IP}^\infty$ of the diffuse limits ($W \to \infty$), as well as those in the glass and the disordered alloy with the same constituent species are indicated by the horizontal lines. In panels (a) and (b) we also show (dashed black lines), the prediction of Eq. (1) for $\kappa_{CP}$ in the incoherent regime, $W > 20$, with the values of $R$ and $l_K$ included in Table 1. The solid curve interpolating the $\kappa_{CP}$ data points in the entire $W$-range is a guide for eye. The calculation of the displayed error bars is detailed in the **Methods** section.</figcaption></figure>

layer A (B). In this situation, phonon propagation is largely obstructed in the cross-plane direction, leading to the observed large reduction of $\kappa_{CP}$. We remark that phonons propagate in the in-plane direction with small constraints, as shown by the large value of $\kappa_{IP}$ close to $\kappa_{IP}^\infty$. 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(\omega)$ 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 $E_k^A + E_k^B = 1$ and $0 \le E_k^k \le 1$, we can define a relative degree of excitation of particles in the two layers, by the threshold value 0.5: large excitations correspond to $E_{A,B}^k \ge 0.5$, small excitations to $E_{A,B}^k < 0.5$. If $E_A^k = E_B^k = 0.5$, particle vibrations in both layers are of the same degree and correlated.

At small $W=2,4$ we find, particularly in the low-$\omega$ region, a large fraction of vibrational states with $E_A^k \approx E_B^k \approx 0.5$. As $W$ increases, in the high frequency region $\omega > \omega_B^{\max}$, where $\omega_B^{\max} \approx 13$ is the high-frequency boundary in $g_B(\omega)$, only particles in the light layer A vibrate ($E_A^k \approx 1$), whereas those in the heavy layer B are almost immobile, as indicated by $E_B^k \approx 0$. In this $\omega$-region, phonon propagation in the cross-plane direction is therefore almost completely suppressed. On the other hand, for $\omega < \omega_B^{\max}$, particles pertaining to the heavy layer B show large vibrational amplitudes ($E_B^k > 0.5$), while vibrations in layer A tend to be small ($E_A^k < 0.5$). More in details, for $W \ge 20$, we see that the averaged amplitudes are much larger in the B layer ($\langle E_B^k \rangle > 0.8$) than in the A layer ($\langle E_A^k \rangle < 0.2$) in the $2 < \omega < 7.5$ range. Contrary to the case of $\omega > \omega_B^{\max}$, however, a significant number of modes are excited in both layers A and B, even with $E_A^k \approx E_B^k \approx 0.5$. We therefore conclude that, for $\omega < \omega_B^{\max}$, some phonons still propagate in the cross-plane direction, contributing to $\kappa_{CP}$.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Thermal conductivity in superlattice SI composed of two intercalated crystalline layers with different masses. The cross-plane, $\kappa_{CP}$, and in-plane, $\kappa_{IP}$, components of thermal conductivity are plotted as functions of the repetition period $W$. The ratio $m_B/m_A$ of the masses in layers A and B is 2 in panel (a), 4 in (b), and 8 in (c). The values $\kappa_{CP}^\infty$ and $\kappa_{IP}^\infty$ of the diffuse limits ($W \to \infty$), as well as those in the glass and the disordered alloy with the same constituent species are indicated by the horizontal lines. In panels (a) and (b) we also show (dashed black lines), the prediction of Eq. (1) for $\kappa_{CP}$ in the incoherent regime, $W > 20$, with the values of $R$ and $\ell_K$ included in Table 1. The solid curve interpolating the $\kappa_{CP}$ data points in the entire $W$-range is a guide for eye. The calculation of the displayed error bars is detailed in the **Methods** section.</figcaption></figure>

layer A (B). In this situation, phonon propagation is largely obstructed in the cross-plane direction, leading to the observed large reduction of $\kappa_{CP}$. We remark that phonons propagate in the in-plane direction with small constraints, as shown by the large value of $\kappa_{IP}$ close to $\kappa_{IP}^\infty$. 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(\omega)$ 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 $E_k^A + E_k^B = 1$ and $0 \le E_k^k \le 1$, we can define a relative degree of excitation of particles in the two layers, by the threshold value 0.5: large excitations correspond to $E_{A,B}^k \ge 0.5$, small excitations to $E_{A,B}^k < 0.5$. If $E_A^k = E_B^k = 0.5$, particle vibrations in both layers are of the same degree and correlated.

At small $W=2,4$ we find, particularly in the low-$\omega$ region, a large fraction of vibrational states with $E_A^k \approx E_B^k \approx 0.5$. As $W$ increases, in the high frequency region $\omega > \omega_B^{\max}$, where $\omega_B^{\max} \approx 13$ is the high-frequency boundary in $g_B(\omega)$, only particles in the light layer A vibrate ($E_A^k \approx 1$), whereas those in the heavy layer B are almost immobile, as indicated by $E_B^k \approx 0$. In this $\omega$-region, phonon propagation in the cross-plane direction is therefore almost completely suppressed. On the other hand, for $\omega < \omega_B^{\max}$, particles pertaining to the heavy layer B show large vibrational amplitudes ($E_B^k > 0.5$), while vibrations in layer A tend to be small ($E_A^k < 0.5$). More in details, for $W \ge 20$, we see that the averaged amplitudes are much larger in the B layer ($\langle E_B^k \rangle > 0.8$) than in the A layer ($\langle E_A^k \rangle < 0.2$) in the $2 < \omega < 7.5$ range. Contrary to the case of $\omega > \omega_B^{\max}$, however, a significant number of modes are excited in both layers A and B, even with $E_A^k \approx E_B^k \approx 0.5$. We therefore conclude that, for $\omega < \omega_B^{\max}$, some phonons still propagate in the cross-plane direction, contributing to $\kappa_{CP}$.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Vibrational density of states in superlattice S1, for a mass ratio $m_B/m_A = 4$ with $m_A = 0.4$ and $m_B = 1.6$. In panels (a–f) we show the data corresponding to the repetitions period values $W = 2,4,10,20,40,80$. For comparison, we also plot $g_{A(B)}(\omega)$ for the homogeneous bulk crystal composed by light (heavy) $m_{A(B)}$ masses only, together with the data for the glass and the disordered alloy formed by the same constituent species.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Vibrational amplitudes of normal modes in superlattice S1, for a mass ratio $m_B/m_A = 4$ with $m_A = 0.4$ and $m_B = 1.6$. The vibrational amplitudes of the eigenvectors, $E_A^k$ and $E_B^k$, in layers A (light) and B (heavy) for all normal modes $k$, are plotted as functions of the corresponding eigenfrequency $\omega^k$. $E_A^k$ and $E_B^k$ are defined in Eq. (6). The repetitions period values are $W = 2,4,10,20,40,80$ in panels (a–f). The solid line represents the average values $(E_A^k)$ and $(E_B^k)$ calculated in bins of the form $\omega^k \pm \delta\omega^k/2$, with $\delta\omega^k = 0.5$. The horizontal dotted lines represent the threshold value $E_A^k = E_B^k = 0.5$, the vertical lines indicate $\omega = \omega_B^{\max} \approx 13$, corresponding to the high frequency edge of $g_B(\omega)$.</figcaption></figure>

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 $\kappa_{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
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Vibrational density of states in superlattice S1, for a mass ratio $m_B/m_A = 4$ with $m_A = 0.4$ and $m_B = 1.6$. In panels (a–f) we show the data corresponding to the repetitions period values $W = 2,4,10,20,40,80$. For comparison, we also plot $g_{A(B)}(\omega)$ for the homogeneous bulk crystal composed by light (heavy) $m_{A(B)}$ masses only, together with the data for the glass and the disordered alloy formed by the same constituent species.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Vibrational amplitudes of normal modes in superlattice S1, for a mass ratio $m_B/m_A = 4$ with $m_A = 0.4$ and $m_B = 1.6$. The vibrational amplitudes of the eigenvectors, $E_A^k$ and $E_B^k$, in layers A (light) and B (heavy) for all normal modes $k$, are plotted as functions of the corresponding eigenfrequency $\omega^k$. $E_A^k$ and $E_B^k$ are defined in Eq. (6). The repetitions period values are $W = 2,4,10,20,40,80$ in panels (a–f). The solid line represents the average values $(E_A^k)$ and $(E_B^k)$ calculated in bins of the form $\omega^k \pm \delta\omega^k/2$, with $\delta\omega^k = 0.5$. The horizontal dotted lines represent the threshold value $E_A^k = E_B^k = 0.5$, the vertical lines indicate $\omega = \omega_B^{\max} \approx 13$, corresponding to the high frequency edge of $g_B(\omega)$.</figcaption></figure>

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 $\kappa_{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
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Thermal conductivity in superlattice S2 composed of ordered crystalline layers intercalated with mass disordered alloy layers. The two components $\kappa_{CP}$ and $\kappa_{IP}$ are plotted as functions of W. The mass ratio of the disordered alloy layer $m_{B2}/m_{B1}$ is (a) 2, (b) 4, and (c) 8. The values $\kappa_{CP}^{\infty}$ and $\kappa_{IP}^{\infty}$ of the diffuse limits, together with those in the glass and the disordered alloy are indicated by lines.</figcaption></figure>

direction, leading to the large reduction of $\kappa_{CP}$. In contrast, in-plane phonon propagation is very mildly affected by the vibrational separation and, therefore, $\kappa_{IP}$ keeps high values. Also, by considering $\langle E_A^k \rangle$ and $\langle E_B^k \rangle$ (solid lines), we conclude that the separation saturates to its maximum level at $W \simeq 20$. Upon further increase $W > 20$, although averaged values show no significant changes, we recognize an increasing fraction of modes with $E_A^k > 0.5$ and $E_B^k < 0.5$ for $\omega < \omega_B^{\max}$ (panels (e) $W=40$ and (f) $W=80$ in Fig. 4). This observation indicates that the separation tendency for modes with $E_A^k < 0.5$ and $E_B^k > 0.5$ 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 $\kappa_{CP}$. Although transport becomes completely incoherent only for values of W of the order of the Kapitza length (note that $\ell_z \simeq 1600$ for $m_B/m_A = 4$), this feature appears as soon as the vibrational separation is saturated, at the crossover point $W^* \simeq 20$. Thus, the saturation point of the vibrational separation identifies the minimum value of $\kappa_{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 $\langle m_B \rangle = (m_{B1} + m_{B2})/2 = 1$ and $\langle m \rangle = (m_A + \langle m_B \rangle)/2 = 1$, respectively, in order to eliminate the trivial contribution associated with different average values in the different cases.

In Fig. 5, we plot $\kappa_{CP}$ and $\kappa_{IP}$ for the mass ratios of the layer B, $m_{B2}/m_{B1}=2$, 4, and 8. At small $W \le 4$, the values of both $\kappa_{CP}$ and $\kappa_{IP}$ are very close to those of the disordered bulk alloy formed by the same particles. As $W$ increases, $\kappa_{IP}$ increases gradually toward $\kappa_{IP}^{\infty}$. This increase is controlled by the development of in-plane phonon propagation in the ordered crystalline layer A. Indeed, the $g(\omega)$ 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(\omega)$. In particular, the longitudinal peak around $\omega \simeq 14.5$ becomes clear, corresponding to that of the crystalline layer A.

The cross-plane value $\kappa_{CP}$ also increases with $W$, but reaches the limit value $\kappa_{CP}^{\infty}$ already at $W \simeq 20$. Since $\kappa_B$ of the disordered alloy layer B is low (see Table 1), $\kappa_{CP}^{\infty}$ remains low, typically less than twice the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Thermal conductivity in superlattice S2 composed of ordered crystalline layers intercalated with mass disordered alloy layers. The two components $\kappa_{CP}$ and $\kappa_{IP}$ are plotted as functions of W. The mass ratio of the disordered alloy layer $m_{B2}/m_{B1}$ is (a) 2, (b) 4, and (c) 8. The values $\kappa_{CP}^{\infty}$ and $\kappa_{IP}^{\infty}$ of the diffuse limits, together with those in the glass and the disordered alloy are indicated by lines.</figcaption></figure>

direction, leading to the large reduction of $\kappa_{CP}$. In contrast, in-plane phonon propagation is very mildly affected by the vibrational separation and, therefore, $\kappa_{IP}$ keeps high values. Also, by considering $\langle E_A^k \rangle$ and $\langle E_B^k \rangle$ (solid lines), we conclude that the separation saturates to its maximum level at $W \simeq 20$. Upon further increase $W > 20$, although averaged values show no significant changes, we recognize an increasing fraction of modes with $E_A^k > 0.5$ and $E_B^k < 0.5$ for $\omega < \omega_B^{\max}$ (panels (e) $W=40$ and (f) $W=80$ in Fig. 4). This observation indicates that the separation tendency for modes with $E_A^k < 0.5$ and $E_B^k > 0.5$ 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 $\kappa_{CP}$. Although transport becomes completely incoherent only for values of W of the order of the Kapitza length (note that $\ell_z \simeq 1600$ for $m_B/m_A = 4$), this feature appears as soon as the vibrational separation is saturated, at the crossover point $W^* \simeq 20$. Thus, the saturation point of the vibrational separation identifies the minimum value of $\kappa_{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 $\langle m_B \rangle = (m_{B1} + m_{B2})/2 = 1$ and $\langle m \rangle = (m_A + \langle m_B \rangle)/2 = 1$, respectively, in order to eliminate the trivial contribution associated with different average values in the different cases.

In Fig. 5, we plot $\kappa_{CP}$ and $\kappa_{IP}$ for the mass ratios of the layer B, $m_{B2}/m_{B1}=2$, 4, and 8. At small $W \le 4$, the values of both $\kappa_{CP}$ and $\kappa_{IP}$ are very close to those of the disordered bulk alloy formed by the same particles. As $W$ increases, $\kappa_{IP}$ increases gradually toward $\kappa_{IP}^{\infty}$. This increase is controlled by the development of in-plane phonon propagation in the ordered crystalline layer A. Indeed, the $g(\omega)$ 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(\omega)$. In particular, the longitudinal peak around $\omega \simeq 14.5$ becomes clear, corresponding to that of the crystalline layer A.

The cross-plane value $\kappa_{CP}$ also increases with $W$, but reaches the limit value $\kappa_{CP}^{\infty}$ already at $W \simeq 20$. Since $\kappa_B$ of the disordered alloy layer B is low (see Table 1), $\kappa_{CP}^{\infty}$ remains low, typically less than twice the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Vibrational density of states in superlattice S2. We show our results for the mass ratio of the disordered alloy layer $m_{B2}/m_{B1} = 4$, with $m_{B1} = 0.4$ and $m_{B2} = 1.6$. The period $W$ is 4, 20, and 40 for (a–c), respectively. For comparison, we plot $g_A(\omega)$ of the bulk crystal formed by particles of mass $m_A = 1$ (crystal A), $g_B(\omega)$ of the disordered alloy with masses $m_{B1} = 0.4$ and $m_{B2} = 1.6$ (alloy B), and the vDOS of the glass and the disordered alloy formed by the same constituent species.</figcaption></figure>

disordered alloy value. As a result, the variation of $\kappa_{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⁵¹ on Si(crystal)-SiGe(disordered alloy) nanowires and numerical simulations⁴⁰ 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³⁷⁻³⁹, or interfacial species mixing⁴⁰,⁴¹ 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 $\kappa_{CP}$ just corresponds to the disordered alloy limit, i.e., $\kappa_{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 $\kappa_{IP}$ showing a behaviour different than that observed in S1 where the anisotropy reaches the maximum at the crossover point $W^* \simeq 20$.

**S3. Superlattice composed by identical crystalline layers separated by weakly interacting interfaces.** In Fig. 7 we show the $W$-dependences of $\kappa_{CP}$ and $\kappa_{IP}$ for the case where the energy scale associated with particles interactions across the interfaces ($\varepsilon_{AB}$) are lowered compared to those intra-layers, with $\varepsilon_{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⁸. 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 $\kappa$ originating from a limited fraction (scaling with the surface of the interfaces) of modified interactions in the glass sample.

The in-plane value $\kappa_{IP}$ is almost independent of $W$, and is very close to the value pertaining to the crystal. In contrast, $\kappa_{CP}$ decreases monotonically by increasing $W$, and especially in the weaker case $\varepsilon_{AB} = 0.1$, the observed reduction of $\kappa_{CP}$ is dramatic. At $W = 10$, $\kappa_{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 $\kappa_{CP}$ is consistent with previous experimental work¹⁹.

Some insight about the origin of this observation comes from the data shown in Fig. 8, where we display the average cross-plane distance $\delta z$ between adjacent crystalline planes (monolayers), normalized to
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Vibrational density of states in superlattice S2. We show our results for the mass ratio of the disordered alloy layer $m_{B2}/m_{B1} = 4$, with $m_{B1} = 0.4$ and $m_{B2} = 1.6$. The period $W$ is 4, 20, and 40 for (a–c), respectively. For comparison, we plot $g_A(\omega)$ of the bulk crystal formed by particles of mass $m_A = 1$ (crystal A), $g_B(\omega)$ of the disordered alloy with masses $m_{B1} = 0.4$ and $m_{B2} = 1.6$ (alloy B), and the vDOS of the glass and the disordered alloy formed by the same constituent species.</figcaption></figure>

disordered alloy value. As a result, the variation of $\kappa_{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⁵¹ on Si.crystal-SiGe(disordered alloy) nanowires and numerical simulations⁴⁰ 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³⁷⁻³⁹, or interfacial species mixing⁴⁰,⁴¹ 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 $\kappa_{CP}$ just corresponds to the disordered alloy limit, i.e., $\kappa_{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 $\kappa_{IP}$ showing a behaviour different than that observed in S1 where the anisotropy reaches the maximum at the crossover point $W^* \simeq 20$.

**S3. Superlattice composed by identical crystalline layers separated by weakly interacting interfaces.** In Fig. 7 we show the $W$-dependences of $\kappa_{CP}$ and $\kappa_{IP}$ for the case where the energy scale associated with particles interactions across the interfaces ($\epsilon_{AB}$) are lowered compared to those intra-layers, with $\epsilon_{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⁸. 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 $\kappa$ originating from a limited fraction (scaling with the surface of the interfaces) of modified interactions in the glass sample.

The in-plane value $\kappa_{IP}$ is almost independent of $W$, and is very close to the value pertaining to the crystal. In contrast, $\kappa_{CP}$ decreases monotonically by increasing $W$, and especially in the weaker case $\epsilon_{AB} = 0.1$, the observed reduction of $\kappa_{CP}$ is dramatic. At $W = 10$, $\kappa_{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 $\kappa_{CP}$ is consistent with previous experimental work¹⁹.

Some insight about the origin of this observation comes from the data shown in Fig. 8, where we display the average cross-plane distance $\delta z$ between adjacent crystalline planes (monolayers), normalized to
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Thermal conductivity in superlattice S3 composed of identical crystalline layers with weakened interface. The two components $κ_{CP}$ and $κ_{IP}$ are plotted as functions of W. The interface interaction $ε_{AB}$ is 0.5 in (a) and 0.1 in (b). We also show, by the horizontal lines, the thermal conductivities of the corresponding one-component homogeneous bulk crystal and glass with unmodified interactions.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Distance between adjacent monolayers in superlattice S3. The average cross-plane distance $δz$ between adjacent crystalline planes plotted for each monolayer, identified by the corresponding order index. We present the value of $δz$ normalized to $a/2$, the horizontal line $δz/(a/2) = 1$ therefore indicates the value in the perfect crystalline lattice. The displacements observed in the cases $W = 20$ are discussed in the main text.</figcaption></figure>

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.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Thermal conductivity in superlattice S3 composed of identical crystalline layers with weakened interface. The two components $\kappa_{CP}$ and $\kappa_{IP}$ are plotted as functions of W. The interface interaction $\epsilon_{AB}$ is 0.5 in (a) and 0.1 in (b). We also show, by the horizontal lines, the thermal conductivities of the corresponding one-component homogeneous bulk crystal and glass with unmodified interactions.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Distance between adjacent monolayers in superlattice S3. The average cross-plane distance $\delta z$ between adjacent crystalline planes plotted for each monolayer, identified by the corresponding order index. We present the value of $\delta z$ normalized to $a/2$, the horizontal line $\delta z/(a/2) = 1$ therefore indicates the value in the perfect crystalline lattice. The displacements observed in the cases $W = 20$ are discussed in the main text.</figcaption></figure>

the value in the perfect lattice, $a/2$. For $\epsilon_{AB} = 0.5$ and $W = 4$, the system keeps the perfect lattice structure, with $\delta z \equiv a/2$ for all monolayers. In contrast, as $\epsilon_{AB}$ decreases and for a large value $W = 20$, $\delta 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 $\delta 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 $\delta z$ at the interfaces located at $W/2 = 10$ and $W = 20$ are different, with a large discrepancy for $\epsilon_{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 $\delta z$ at different interfaces.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Vibrational density of states and vibrational amplitudes in superlattice S3. We report data corresponding to the indicated values of the interfacial interaction energy $\epsilon_{AB}$ and the repetition period $W$: (a) $\epsilon_{AB}=0.5$, $W=4$, (b) $\epsilon_{AB}=0.5$, $W=20$, and (c) $\epsilon_{AB}=0.1$, $W=20$. In panels at the top, we show the vDOS $g(\omega)$ for the superlattices of type S3, together with those of the corresponding one-component homogeneous crystal and glass with unmodified interactions. In panels at the bottom we show the vibrational amplitudes, $E_A^k$ and $E_B^k$, in layers A and B, plotted as functions of the eigenfrequency $\omega^k$. The solid line represents the average values $\langle E_A^k \rangle$ and $\langle E_B^k \rangle$, calculated in bins of the form $\omega^k \pm \delta\omega^k/2$, with $\delta\omega^k = 0.5$. The horizontal dotted lines indicate $E_A^k = E_B^k = 0.5$.</figcaption></figure>

The behaviour of $\kappa_{CP}$ can be further elucidated by inspection of the main features of the vibrational spectrum. In Fig. 9 we plot the $g(\omega)$ 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(\omega)$ (top panels) for all cases, similar to the homogeneous bulk crystal. For $\epsilon_{AB}=0.5$, $W=4$, $g(\omega)$ shows an excess of lower-$\omega$ modes compared to those present in the one-component crystal, simply due to the weakened interactions at the interfaces. As $\epsilon_{AB}$ decreases and $W$ increases, $g(\omega)$ 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 $\epsilon_{AB}=0.5$, $W=20$ and $\epsilon_{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(\omega)$ 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 $\epsilon_{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 $E_A^k \equiv E_B^k \equiv 0.5$. 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 $\kappa_{CP}$. The situation changes drastically in the case $\epsilon_{AB}=0.1$ and $W=20$, where the ultra-low value of $\kappa_{CP}$ can be reached. Except for the low-$\omega$ modes, $E_A^k$ and $E_B^k$ are symmetrically randomly distributed around the average values $\langle E_{A(B)}^k \rangle \equiv 0.5$, 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 $\kappa_{CP}$, while keeping a high $\kappa_{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
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Vibrational density of states and vibrational amplitudes in superlattice S3. We report data corresponding to the indicated values of the interfacial interaction energy $\epsilon_{AB}$ and the repetition period $W$: (a) $\epsilon_{AB}=0.5$, $W=4$, (b) $\epsilon_{AB}=0.5$, $W=20$, and (c) $\epsilon_{AB}=0.1$, $W=20$. In panels at the top, we show the vDOS $g(\omega)$ for the superlattices of type S3, together with those of the corresponding one-component homogeneous crystal and glass with unmodified interactions. In panels at the bottom we show the vibrational amplitudes, $E_A^k$ and $E_B^k$, in layers A and B, plotted as functions of the eigenfrequency $\omega^k$. The solid line represents the average values $\langle E_A^k \rangle$ and $\langle E_B^k \rangle$, calculated in bins of the form $\omega^k \pm \delta\omega^k/2$, with $\delta\omega^k = 0.5$. The horizontal dotted lines indicate $E_A^k = E_B^k = 0.5$.</figcaption></figure>

The behaviour of $\kappa_{CP}$ can be further elucidated by inspection of the main features of the vibrational spectrum. In Fig. 9 we plot the $g(\omega)$ 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(\omega)$ (top panels) for all cases, similar to the homogeneous bulk crystal. For $\epsilon_{AB}=0.5$, $W=4$, $g(\omega)$ shows an excess of lower-$\omega$ modes compared to those present in the one-component crystal, simply due to the weakened interactions at the interfaces. As $\epsilon_{AB}$ decreases and $W$ increases, $g(\omega)$ 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 $\epsilon_{AB}=0.5$, $W=20$ and $\epsilon_{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(\omega)$ 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 $\epsilon_{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 $E_A^k \equiv E_B^k \equiv 0.5$. 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 $\kappa_{CP}$. The situation changes drastically in the case $\epsilon_{AB}=0.1$ and $W=20$, where the ultra-low value of $\kappa_{CP}$ can be reached. Except for the low-$\omega$ modes, $E_A^k$ and $E_B^k$ are symmetrically randomly distributed around the average values $\langle E_{A(B)}^k \rangle \equiv 0.5$, 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 $\kappa_{CP}$, while keeping a high $\kappa_{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 $\kappa_{CP}$ below the corresponding glass limit, while keeping the in-plane contribution $\kappa_{IP}$ close to the pure crystalline value.

More specifically, in the case of mass mismatch (S1), propagation of phonons with high frequencies ($\omega > \omega_B^{max}$) is almost completely suppressed, whereas a fraction of low-frequency phonons ($\omega < \omega_B^{max}$) are still able to propagate across the interfaces, contributing to $\kappa_{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 $\kappa_{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 $\kappa_{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-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-$\omega$ (long wavelength, $\lambda$) phonons propagation, than those with high-$\omega$ (short $\lambda$). This situation is similar to what has been observed in bulk glasses, where the long-$\lambda$ 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-$\lambda$ 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(quan- tum)-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 $T_g$, avoiding crystallization (see, for instance, Ref. 49 for details on the preparation of glasses). Particles, i and j, interact via soft-sphere (SS) or Lennard-Jones (LJ) potentials:

$$v_{SS}^{ij}(r) = \epsilon^{ij} \left( \frac{\sigma^{ij}}{r} \right)^{12}, \\ v_{LJ}^{ij}(r) = 4\epsilon^{ij} \left[ \left( \frac{\sigma^{ij}}{r} \right)^{12} - \left( \frac{\sigma^{ij}}{r} \right)^{6} \right], \quad (3)$$
in the in-plane direction. This reduces the cross-plane thermal conductivity $\kappa_{CP}$ below the corresponding glass limit, while keeping the in-plane contribution $\kappa_{IP}$ close to the pure crystalline value.

More specifically, in the case of mass mismatch ($S1$), propagation of phonons with high frequencies ($\omega > \omega_B^{max}$) is almost completely suppressed, whereas a fraction of low-frequency phonons ($\omega < \omega_B^{max}$) are still able to propagate across the interfaces, contributing to $\kappa_{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 $\kappa_{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 $\kappa_{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-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-$\omega$ (long wavelength, $\lambda$) phonons propagation, than those with high-$\omega$ (short $\lambda$). This situation is similar to what has been observed in bulk glasses, where the long-$\lambda$ 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-$\lambda$ 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(quan-tum)-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 $T_g$, avoiding crystallization (see, for instance, Ref. 49 for details on the preparation of glasses). Particles, i and j, interact via soft-sphere (SS) or Lennard-Jones (LJ) potentials:

$$v_{SS}^{ij}(r) = \epsilon^{ij} \left( \frac{\sigma^{ij}}{r} \right)^{12}, \\ v_{LJ}^{ij}(r) = 4\epsilon^{ij} \left[ \left( \frac{\sigma^{ij}}{r} \right)^{12} - \left( \frac{\sigma^{ij}}{r} \right)^{6} \right], \qquad (3)$$
where *r* is the distance between those two particles, and $\sigma^{ij}$ and $\varepsilon^{ij}$ are the interparticle diameter and interaction energy scale, respectively. The potential is cut-off and shifted at $r_c = 2.5\sigma^{ij}$. Particle *i* has mass $m^i$, and we have used $\sigma$, $\varepsilon/k_B$ ($k_B$ is the Boltzmann constant), and *m* as units of length, temperature, and mass. As a reference, for Argon $\sigma = 3.4$ Å, $\varepsilon/k_B = 120$ K, and $m = 39.96$ a.m.u. We considered the number density $\hat{\rho} = N/V=1.015$, corresponding to a lattice constant $a = (4/\hat{\rho})^{1/3} = 1.58$.

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 $v_{SS}^{ij}(r)$, with $\sigma^{ij} = \varepsilon^{ij} = 1$.

The second superlattice (S2) is composed of an ordered crystalline layer A intercalated to a disordered alloy layer B. $m_A = 1$ in A, whereas in B half of the particles have mass $m_{B1}$, $m_{B2}$ the others, and are randomly distributed on the lattice sites. Again, $m_{B1}$ and $m_{B2}$ are determined by the mass ratio $m_{B2}/m_{B1} > 1$, keeping a constant average value ($m_{B1} + m_{B2})/2 = 1$. All particles in both layers interact via the SS potential $v_{SS}^{ij}(r)$ with $\sigma^{ij} = \varepsilon^{ij} = 1$.

The third superlattice (S3) is composed of identical crystalline layers A and B, 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 $m_A = m_B = 1$, and interact via the LJ potential $v_{LJ}^{ij}(r)$, with $\sigma^{ij} = \varepsilon^{ij} = 1$. The energy scale of interactions between particles pertaining to different layers are, however, reduced to $\varepsilon^{ij} = \varepsilon_{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-$\omega$ 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:

$$ \begin{aligned} \kappa_{CP} &= \frac{1}{VT^2} \int_0^\infty \langle J_z(t) J_z(0) \rangle dt, \\ \kappa_{IP} &= \frac{1}{VT^2} \int_0^\infty \langle J_{x,y}(t) J_{x,y}(0) \rangle dt. \end{aligned} \quad (4) $$

Here, $J_{x,y}$ and $J_z$ are the heat currents in the in-plane $(x,y)$ and cross-plane $(z)$ directions, and $\langle \rangle$ 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, $\kappa_{CP} \simeq \kappa_{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, $\delta t = 5 \times 10^{-3}$, for a total run-time $t_{run} = N_{run}\delta t = 10^5$ ($N_{run} = 2 \times 10^7$ steps) for S1 and $t_{run} = 10^4$ ($N_{run} = 2 \times 10^6$ steps) for both S2 and S3. The systems snapshots extracted from the trajectory have been used to calculate the correlation functions, $\langle J_{x,y,z}(t) J_{x,y,z}(0) \rangle$ of Eq. (4), which have been subsequently integrated numerically over a finite time window $\Delta t$. We have observed a clear convergence of the integrals for $\kappa_{CP,IP}$ to the limits of Eq. (4) for $\Delta t = 10^4$ for S1, and $\Delta 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 $\kappa_{CP}$ and $2 \times 10 = 20$ (we have considered both the $x$ and $y$ components) for $\kappa_{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 $\lambda > 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 $\kappa$. The box size therefore needs to be large enough to include a vibrational spectrum sufficient to establish an accurate description of
where r is the distance between those two particles, and $\sigma^{ij}$ and $\varepsilon^{ij}$ are the interparticle diameter and interaction energy scale, respectively. The potential is cut-off and shifted at $r_c = 2.5\sigma^{ij}$. Particle i has mass $m^i$, and we have used $\sigma$, $\varepsilon/k_B$ ($k_B$ is the Boltzmann constant), and m as units of length, temperature, and mass. As a reference, for Argon $\sigma = 3.4\AA$, $\varepsilon/k_B = 120$ K, and $m = 39.96$ a.m.u. We considered the number density $\hat{\rho} = N/V=1.015$, corresponding to a lattice constant $a = (4/\hat{\rho})^{1/3} = 1.58$.

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 $v_{SS}^{ij}(r)$, with $\sigma^{ij} = \varepsilon^{ij} = 1$.

The second superlattice (S2) is composed of an ordered crystalline layer A intercalated to a disordered alloy layer B. $m_A = 1$ in A, whereas in B half of the particles have mass $m_{B1}$, $m_{B2}$ the others, and are randomly distributed on the lattice sites. Again, $m_{B1}$ and $m_{B2}$ are determined by the mass ratio $m_{B2}/m_{B1} > 1$, keeping a constant average value ($m_{B1} + m_{B2})/2 = 1$. All particles in both layers interact via the SS potential $v_{SS}^{ij}(r)$ with $\sigma^{ij} = \varepsilon^{ij} = 1$.

The third superlattice (S3) is composed of identical crystalline layers A and B, 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 $m_A = m_B = 1$, and interact via the LJ potential $v_{LJ}^{ij}(r)$, with $\sigma^{ij} = \varepsilon^{ij} = 1$. The energy scale of interactions between particles pertaining to different layers are, however, reduced to $\varepsilon^{ij} = \varepsilon_{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-$\omega$ 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:

$$ \begin{aligned} \kappa_{CP} &= \frac{1}{VT^2} \int_0^\infty \langle J_z(t) J_z(0) \rangle dt, \\ \kappa_{IP} &= \frac{1}{VT^2} \int_0^\infty \langle J_{x,y}(t) J_{x,y}(0) \rangle dt. \end{aligned} \quad (4) $$

Here, $J_{x,y}$ and $J_z$ are the heat currents in the in-plane $(x,y)$ and cross-plane $(z)$ directions, and $\langle \rangle$ 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, $\kappa_{CP} \simeq \kappa_{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, $\delta t = 5 \times 10^{-3}$, for a total run-time $t_{run} = N_{run}\delta t = 10^5$ ($N_{run} = 2 \times 10^7$ steps) for S1 and $t_{run} = 10^4$ ($N_{run} = 2 \times 10^6$ steps) for both S2 and S3. The systems snapshots extracted from the trajectory have been used to calculate the correlation functions, $\langle J_{x,y,z}(t) J_{x,y,z}(0) \rangle$ of Eq. (4), which have been subsequently integrated numerically over a finite time window $\Delta t$. We have observed a clear convergence of the integrals for $\kappa_{CP,IP}$ to the limits of Eq. (4) for $\Delta t = 10^4$ for S1, and $\Delta 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 $\kappa_{CP}$ and $2 \times 10 = 20$ (we have considered both the x and y components) for $\kappa_{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 $\lambda > 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 $\kappa$. 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⁴⁷. 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 $\kappa_{CP}$ and $\kappa_{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 $\kappa_{CP} \simeq \kappa_{IP}$, without any size effect⁸.

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 \ge 20$, whereas four periods or more ($L \ge 4W$) are required for $W \le 8$ (see the Supplementary Materials). We have therefore employed four pattern repetitions ($L=4W$) for $10 \le W \le 12$ and two ($L=2W$) for $14 \le W \le 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^* \simeq 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 \ge 20$, while two or more repetitions ($L \ge 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 $\lambda > 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⁸,⁴⁹. We have diagonalized the dynamical (Hessian) matrix calculated at local minima of the potential energy landscape, and obtained eigenvalues $\lambda^k$ and eigenvectors (polarization vectors) $\mathbf{e}^k = \{\mathbf{e}_1^k, \dots, \mathbf{e}_j^k, \dots, \mathbf{e}_N^k\}$. Here, $j$ is the particle index, and $k=1, 2, \dots, 3N-3$ is the eigenmode number, where we have disregarded the three vanishing Goldstone modes. The eigenvectors are normalized such that $\mathbf{e}^k \cdot \mathbf{e}^l = \sum_{j=1}^N (\mathbf{e}_j^k \cdot \mathbf{e}_j^l) = \delta_{kl}$, where $\delta_{kl}$ is the Kronecker delta function. The eigenfrequencies are next calculated as $\omega^k = \sqrt{\lambda^k}$, and the associated probability distribution (normalized histogram) directly provides the vDOS:

$$g(\omega) = \frac{1}{3N-3} \sum_{k=1}^{3N-3} \delta(\omega - \omega^k). \qquad (5)$$

In addition, from the eigenvector $\mathbf{e}^k$ we have defined the vibrational amplitudes of mode $k$ for layers A and B:

$$E_A^k(B) = \sum_{j \in \text{layer } A(B)} (\mathbf{e}_j^k \cdot \mathbf{e}_j^k). \qquad (6)$$

Note that $E_A^k + E_B^k = \mathbf{e}^k \cdot \mathbf{e}^k = 1$ for each $k$ and, therefore, $0 \le E_A^k, E_B^k \le 1$. Based on the values of $E_A^k$ and $E_B^k$, one can determine in which layer particles are more displaced (excited) according to the associated eigenvector $\mathbf{e}^k$. In particular, if $E_A^k \ge 0.5, E_B^k < 0.5$ ($E_A^k < 0.5, E_B^k \ge 0.5$), particles in layer A (B) contribute more to mode $k$ than those in layer B (A). In the case $E_A^k = E_B^k = 0.5$, 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 \to 0$ which, we believe, is an appropriate approximation for our case $T = 10^{-2}$, where anharmonicities are weak.

## References

1. Venkatasubramanian, R., Siivola, E., Colpitts, T. & O'Quinn, B. Thin-film thermoelectric devices with high room-temperature figures of merit. *Nature* **413**, 597–602 (2001).

2. Minnich, A. J., Dresselhaus, M. S., Ren, Z. F. & Chen, G. Bulk nanostructured thermoelectric materials: current research and future prospects. *Energy Environ. Sci.* **2**, 466–479 (2009).

3. Maldovan, M. Sound and heat revolutions in phononics. *Nature* **503**, 209–217 (2013).

4. Goodson, K. E. Ordering up the minimum thermal conductivity of solids. *Science* **315**, 342–343 (2007).

5. Cahill, D. G. & Pohl, R. O. Lattice vibrations and heat transport in crystals and glasses. *Ann. Rev. Phys. Chem.* **39**, 93–121 (1988).

6. Cahill, D. G., Watson, S. K. & Pohl, R. O. Lower limit to the thermal conductivity of disordered crystals. *Phys. Rev. B* **46**, 6131–6140 (1992).
anharmonic coupling (scattering) processes⁴⁷. 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 $\kappa_{CP}$ and $\kappa_{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 $\kappa_{CP} \simeq \kappa_{IP}$, without any size effect⁸.

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 \ge 20$, whereas four periods or more ($L \ge 4W$) are required for $W \le 8$ (see the Supplementary Materials). We have therefore employed four pattern repetitions ($L=4W$) for $10 \le W \le 12$ and two ($L=2W$) for $14 \le W \le 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^* \simeq 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 \ge 20$, while two or more repetitions ($L \ge 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 $\lambda > 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⁸,⁴⁹. We have diagonalized the dynamical (Hessian) matrix calculated at local minima of the potential energy landscape, and obtained eigenvalues $\lambda^k$ and eigenvectors (polarization vectors) $\mathbf{e}^k = \{\mathbf{e}_1^k, \dots, \mathbf{e}_j^k, \dots, \mathbf{e}_N^k\}$. Here, $j$ is the particle index, and $k=1, 2, \dots, 3N-3$ is the eigenmode number, where we have disregarded the three vanishing Goldstone modes. The eigenvectors are normalized such that $\mathbf{e}^k \cdot \mathbf{e}^l = \sum_{j=1}^N (\mathbf{e}_j^k \cdot \mathbf{e}_j^l) = \delta_{kl}$, where $\delta_{kl}$ is the Kronecker delta function. The eigenfrequencies are next calculated as $\omega^k = \sqrt{\lambda^k}$, and the associated probability distribution (normalized histogram) directly provides the vDOS:

$$g(\omega) = \frac{1}{3N-3} \sum_{k=1}^{3N-3} \delta(\omega - \omega^k). \qquad (5)$$

In addition, from the eigenvector $\mathbf{e}^k$ we have defined the vibrational amplitudes of mode $k$ for layers A and B:

$$E_A^k(B) = \sum_{j \in \text{layer} A(B)} (\mathbf{e}_j^k \cdot \mathbf{e}_j^k). \qquad (6)$$

Note that $E_A^k + E_B^k = \mathbf{e}^k \cdot \mathbf{e}^k = 1$ for each $k$ and, therefore, $0 \le E_A^k, E_B^k \le 1$. Based on the values of $E_A^k$ and $E_B^k$, one can determine in which layer particles are more displaced (excited) according to the associated eigenvector $\mathbf{e}^k$. In particular, if $E_A^k \ge 0.5, E_B^k < 0.5$ ($E_A^k < 0.5, E_B^k \ge 0.5$), particles in layer A (B) contribute more to mode $k$ than those in layer B (A). In the case $E_A^k = E_B^k = 0.5$, 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 \to 0$ which, we believe, is an appropriate approximation for our case $T = 10^{-2}$, where anharmonicities are weak.

## References

1. Venkatasubramanian, R., Siivola, E., Colpitts, T. & O'Quinn, B. Thin-film thermoelectric devices with high room-temperature figures of merit. *Nature* **413**, 597–602 (2001).

2. Minnich, A. J., Dresselhaus, M. S., Ren, Z. F. & Chen, G. Bulk nanostructured thermoelectric materials: current research and future prospects. *Energy Environ. Sci.* **2**, 466–479 (2009).

3. Maldovan, M. Sound and heat revolutions in phononics. *Nature* **503**, 209–217 (2013).

4. Goodson, K. E. Ordering up the minimum thermal conductivity of solids. *Science* **315**, 342–343 (2007).

5. Cahill, D. G. & Pohl, R. O. Lattice vibrations and heat transport in crystals and glasses. *Ann. Rev. Phys. Chem.* **39**, 93–121 (1988).

6. Cahill, D. G., Watson, S. K. & Pohl, R. O. Lower limit to the thermal conductivity of disordered crystals. *Phys. Rev. B* **46**, 6131–6140 (1992).
7. Allen, P. B. & Feldman, J. L. Thermal conductivity of disordered harmonic solids. *Phys. Rev. B* **48**, 12581–12588 (1993).

8. Mizuno, H., Mossa, S. & Barrat, J.-L. Elastic heterogeneity, vibrational states, and thermal conductivity across an amorphisation transition. *EPL* **104**, 56001 (2013).

9. Kittel, C. *Introduction to Solid State Physics* (John Wiley and Sons, New York, 1996), 7th edn.

10. Mizuno, H., Mossa, S. & Barrat, J.-L. Measuring spatial distribution of the local elastic modulus in glasses. *Phys. Rev. E* **87**, 042306 (2013).

11. Hopkins, P. E. et al. Reduction in the thermal conductivity of single crystalline silicon by phononic crystal patterning. *Nano Lett.* **11**, 107–112 (2011).

12. Lee, S.-M., Cahill, D. G. & Venkatasubramanian, R. Thermal conductivity of Si-Ge superlattices. *Appl. Phys. Lett.* **70**, 2957–2959 (1997).

13. Volz, S., Saulnier, J. B., Chen, G. & Beauchamp, P. Computation of thermal conductivity of Si/Ge superlattices by molecular dynamics techniques. *Microelectronics Journal* **31**, 815–819 (2000).

14. Capinski, W. S. et al. Thermal-conductivity measurements of GaAs/AlAs superlattices using a picosecond optical pump-and-probe technique. *Phys. Rev. B* **59**, 8105–8113 (1999).

15. Daly, B. C. & Maris, H. J. Calculation of the thermal conductivity of superlattices by molecular dynamics simulation. *Phys. B: Condens. Matter* **316–317**, 247–249 (2002).

16. Yang, B., Liu, W. L., Liu, J. L., Wang, K. L. & Chen, G. Measurements of anisotropic thermoelectric properties in superlattices. *Appl. Phys. Lett.* **81**, 3588–3590 (2002).

17. Mavrokefalos, A., Nguyen, N. T., Pettes, M. T., Johnson, D. C. & Shi, L. In-plane thermal conductivity of disordered layered WSe₂ and (W₃)(WSe₂)₃ superlattice films. *Appl. Phys. Lett.* **91**, 171912 (2007).

18. Costescu, R. M., Cahill, D. G., Fabreguette, F. H., Sechrist, Z. A. & George, S. M. Ultra-low thermal conductivity in W/Al₂O₃ nanolaminates. *Science* **303**, 989–990 (2004).

19. Chiritescu, C. et al. Ultralow thermal conductivity in disordered, layered WSe₂ crystals. *Science* **315**, 351–353 (2007).

20. Pernot, G. et al. Precise control of thermal conductivity at the nanoscale through individual phonon-scattering barriers. *Nature Mater.* **9**, 491–495 (2010).

21. Wang, Y., Gu, C., Ruan, X. Optimization of the random multilayer structure to break the random-alloy limit of thermal conductivity. *Appl. Phys. Lett.* **106**, 073104 (2015).

22. Simkin, M. V. & Mahan, G. D. Minimum thermal conductivity of superlattices. *Phys. Rev. Lett.* **84**, 927–930 (2000).

23. Yang, B. & Chen, G. Partially coherent phonon heat conduction in superlattices. *Phys. Rev. B* **67**, 195311 (2003).

24. Garg, J. & Chen, G. Minimum thermal conductivity in superlattices: A first-principles formalism. *Phys. Rev. B* **87**, 140302 (2013).

25. Chen, Y., Li, D., Lukes, J. R., Ni, Z. & Chen, M. Minimum superlattice thermal conductivity from molecular dynamics. *Phys. Rev. B* **72**, 174302 (2005).

26. Kawamura, T., Kangawa, Y. & Kakimoto, K. An investigation of thermal conductivity of nitride-semiconductor nanostructures by molecular dynamics simulation. *J. Crystal Growth* **298**, 251–253 (2007).

27. Yang, N., Zhang, G. & Li, B. Ultralow thermal conductivity of isotope-doped silicon nanowires. *Nano Lett.* **8**, 276–280 (2008).

28. Ravichandran, J. et al. Crossover from incoherent to coherent phonon scattering in epitaxial oxide superlattices. *Nature Mater.* **13**, 168–172 (2014).

29. Chen, G. & Neagu, M. Thermal conductivity and heat transfer in superlattices. *Appl. Phys. Lett.* **71**, 2761–2763 (1997).

30. Chen, G. Thermal conductivity and ballistic-phonon transport in the cross-plane direction of superlattices. *Phys. Rev. B* **57**, 14958–14973 (1998).

31. Kim, E.-K., Kwun, S.-I., Lee, S.-M., Seo, H. & Yoon, J.-G. Thermal boundary resistance at Ge₂Sb₂Te₅/ZnS:SiO₂ interface. *Appl. Phys. Lett.* **76**, 3864–3866 (2000).

32. Lampin, E., Nguyen, Q.-H., Francioso, P. A. & Cleri, F. Thermal boundary resistance at silicon-silica interfaces by molecular dynamics simulations. *Appl. Phys. Lett.* **100**, 131906 (2012).

33. Chen, J., Walther, J. H. & Koumoutsakos, P. Strain Engineering of Kapitza Resistance in Few-Layer Graphene. *Nano Lett.* **14**, 819–825 (2014).

34. Barrat, J.-L., & Chiaruttini, F. Kapitza resistance at the liquid-solid interface. *Mol. Phys.* **101**, 1605–1610 (2003).

35. Tamura, S., Hurley, D. C. & Wolfe, J. P. Acoustic-phonon propagation in superlattices. *Phys. Rev. B* **38**, 1427–1449 (1988).

36. Ren, S. Y. & Dow, J. D. Thermal conductivity of superlattices. *Phys. Rev. B* **25**, 3750–3755 (1982).

37. Daly, B. C., Maris, H. J., Imamura, K. & Tamura, S. Molecular dynamics calculation of the thermal conductivity of superlattices. *Phys. Rev. B* **66**, 024301 (2002).

38. Imamura, K., Tanaka, Y., Nishiguchi, N., Tamura, S. & Maris, H. J. Lattice thermal conductivity in superlattices: molecular dynamics calculations with a heat reservoir method. *J. Phys.: Condens. Matter* **15**, 8679–8690 (2003).

39. Daly, B. C., Maris, H. J., Tanaka, Y. & Tamura, S. Molecular dynamics calculation of the in-plane thermal conductivity of GaAs/AlAs superlattices. *Phys. Rev. B* **67**, 033308 (2003).

40. Landry, E. S. & McGaughey, A. J. H. Effect of interfacial species mixing on phonon transport in semiconductor superlattices. *Phys. Rev. B* **79**, 075316 (2009).

41. Huberman, S. C., Larkin, J. M., McGaughey, A. J. H. & Amon, C. H. Disruption of superlattice phonons by interfacial mixing. *Phys. Rev. B* **88**, 155311 (2013).

42. Termentzidis, K., Chantrenne, P. & Keblinski, P. Nonequilibrium molecular dynamics simulation of the in-plane thermal conductivity of superlattices with rough interfaces. *Phys. Rev. B* **79**, 214307 (2009).

43. Termentzidis, K., Merabia, S., Chantrenne, P. & Keblinski, P. 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 (2011).

44. Hsieh, W.-P., Lyons, A. S., Pop, E., Keblinski, P. & Cahill, D. G. Pressure tuning of the thermal conductance of weak interfaces. *Phys. Rev. B* **84**, 184107 (2011).

45. Shen, M., Evans, W. J., Cahill, D. & Keblinski, P. Bonding and pressure-tunable interfacial thermal conductance. *Phys. Rev. B* **84**, 195432 (2011).

46. Losego, M. D., Grady, M. E., Sottos, N. R., Cahill, D. G. & Braun, P. V. Effects of chemical bonding on heat transport across interfaces. *Nature Mater.* **11**, 502–506 (2012).

47. McGaughey, A. J. H. & Kaviani, M. Phonon transport in molecular dynamics simulations: formulation and thermal conductivity prediction. *Advances in Heat Transfer* **39**, 169–255 (2006).

48. Landry, E. S., Hussein, M. I. & McGaughey, A. J. H. Complex superlattice unit cell designs for reduced thermal conductivity. *Phys. Rev. B* **77**, 184302 (2008).

49. Monaco, G. & Mossa, S. Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale. *Proc. Natl. Acad. Sci. USA* **106**, 16907–16912 (2009).

50. Yang, L., Yang, N. & Li, B. Reduction of thermal conductivity by nanoscale 3d phononic crystal. *Sci. Rep.* **3**, 1143 (2012).

51. Li, D., Wu, Y., Fan, R., Yang, P. & Majumdar, A. Thermal conductivity of Si/SiGe superlattice nanowires. *Appl. Phys. Lett.* **83**, 3186–3188 (2003).
7. Allen, P. B. & Feldman, J. L. Thermal conductivity of disordered harmonic solids. *Phys. Rev. B* **48**, 12581–12588 (1993).

8. Mizuno, H., Mossa, S. & Barrat, J.-L. Elastic heterogeneity, vibrational states, and thermal conductivity across an amorphisation transition. *EPL* **104**, 56001 (2013).

9. Kittel, C. *Introduction to Solid State Physics* (John Wiley and Sons, New York, 1996), 7th edn.

10. Mizuno, H., Mossa, S. & Barrat, J.-L. Measuring spatial distribution of the local elastic modulus in glasses. *Phys. Rev. E* **87**, 042306 (2013).

11. Hopkins, P. E. et al. Reduction in the thermal conductivity of single crystalline silicon by phononic crystal patterning. *Nano Lett.* **11**, 107–112 (2011).

12. Lee, S.-M., Cahill, D. G. & Venkatasubramanian, R. Thermal conductivity of Si-Ge superlattices. *Appl. Phys. Lett.* **70**, 2957–2959 (1997).

13. Volz, S., Saulnier, J. B., Chen, G. & Beauchamp, P. Computation of thermal conductivity of Si/Ge superlattices by molecular dynamics techniques. *Microelectronics Journal* **31**, 815–819 (2000).

14. Capinski, W. S. et al. Thermal-conductivity measurements of GaAs/AlAs superlattices using a picosecond optical pump-and-probe technique. *Phys. Rev. B* **59**, 8105–8113 (1999).

15. Daly, B. C. & Maris, H. J. Calculation of the thermal conductivity of superlattices by molecular dynamics simulation. *Phys. B: Condens. Matter* **316–317**, 247–249 (2002).

16. Yang, B., Liu, W. L., Liu, J. L., Wang, K. L. & Chen, G. Measurements of anisotropic thermoelectric properties in superlattices. *Appl. Phys. Lett.* **81**, 3588–3590 (2002).

17. Mavrokefalos, A., Nguyen, N. T., Pettes, M. T., Johnson, D. C. & Shi, L. In-plane thermal conductivity of disordered layered WSe₂ and (Wₓ)(WSe₂)₃ superlattice films. *Appl. Phys. Lett.* **91**, 171912 (2007).

18. Costescu, R. M., Cahill, D. G., Fabreguette, F. H., Sechrist, Z. A. & George, S. M. Ultra-low thermal conductivity in W/Al₂O₃ nanolaminates. *Science* **303**, 989–990 (2004).

19. Chiritescu, C. et al. Ultralow thermal conductivity in disordered, layered WSe₂ crystals. *Science* **315**, 351–353 (2007).

20. Pernot, G. et al. Precise control of thermal conductivity at the nanoscale through individual phonon-scattering barriers. *Nature Mater.* **9**, 491–495 (2010).

21. Wang, Y., Gu, C., Ruan, X. Optimization of the random multilayer structure to break the random-alloy limit of thermal conductivity. *Appl. Phys. Lett.* **106**, 073104 (2015).

22. Simkin, M. V. & Mahan, G. D. Minimum thermal conductivity of superlattices. *Phys. Rev. Lett.* **84**, 927–930 (2000).

23. Yang, B. & Chen, G. Partially coherent phonon heat conduction in superlattices. *Phys. Rev. B* **67**, 195311 (2003).

24. Garg, J. & Chen, G. Minimum thermal conductivity in superlattices: A first-principles formalism. *Phys. Rev. B* **87**, 140302 (2013).

25. Chen, Y., Li, D., Lukes, J. R., Ni, Z. & Chen, M. Minimum superlattice thermal conductivity from molecular dynamics. *Phys. Rev. B* **72**, 174302 (2005).

26. Kawamura, T., Kangawa, Y. & Kakimoto, K. An investigation of thermal conductivity of nitride-semiconductor nanostructures by molecular dynamics simulation. *J. Crystal Growth* **298**, 251–253 (2007).

27. Yang, N., Zhang, G. & Li, B. Ultralow thermal conductivity of isotope-doped silicon nanowires. *Nano Lett.* **8**, 276–280 (2008).

28. Ravichandran, J. et al. Crossover from incoherent to coherent phonon scattering in epitaxial oxide superlattices. *Nature Mater.* **13**, 168–172 (2014).

29. Chen, G. & Neagu, M. Thermal conductivity and heat transfer in superlattices. *Appl. Phys. Lett.* **71**, 2761–2763 (1997).

30. Chen, G. Thermal conductivity and ballistic-phonon transport in the cross-plane direction of superlattices. *Phys. Rev. B* **57**, 14958–14973 (1998).

31. Kim, E.-K., Kwun, S.-I., Lee, S.-M., Seo, H. & Yoon, J.-G. Thermal boundary resistance at Ge₂Sb₂Te₅/ZnS:SiO₂ interface. *Appl. Phys. Lett.* **76**, 3864–3866 (2000).

32. Lampin, E., Nguyen, Q.-H., Francioso, P. A. & Cleri, F. Thermal boundary resistance at silicon-silica interfaces by molecular dynamics simulations. *Appl. Phys. Lett.* **100**, 131906 (2012).

33. Chen, J., Walther, J. H. & Koumoutsakos, P. Strain Engineering of Kapitza Resistance in Few-Layer Graphene. *Nano Lett.* **14**, 819–825 (2014).

34. Barrat, J.-L., & Chiaruttini, F. Kapitza resistance at the liquid-solid interface. *Mol. Phys.* **101**, 1605–1610 (2003).

35. Tamura, S., Hurley, D. C. & Wolfe, J. P. Acoustic-phonon propagation in superlattices. *Phys. Rev. B* **38**, 1427–1449 (1988).

36. Ren, S. Y. & Dow, J. D. Thermal conductivity of superlattices. *Phys. Rev. B* **25**, 3750–3755 (1982).

37. Daly, B. C., Maris, H. J., Imamura, K. & Tamura, S. Molecular dynamics calculation of the thermal conductivity of superlattices. *Phys. Rev. B* **66**, 024301 (2002).

38. Imamura, K., Tanaka, Y., Nishiguchi, N., Tamura, S. & Maris, H. J. Lattice thermal conductivity in superlattices: molecular dynamics calculations with a heat reservoir method. *J. Phys.: Condens. Matter* **15**, 8679–8690 (2003).

39. Daly, B. C., Maris, H. J., Tanaka, Y. & Tamura, S. Molecular dynamics calculation of the in-plane thermal conductivity of GaAs/AlAs superlattices. *Phys. Rev. B* **67**, 033308 (2003).

40. Landry, E. S. & McGaughey, A. J. H. Effect of interfacial species mixing on phonon transport in semiconductor superlattices. *Phys. Rev. B* **79**, 075316 (2009).

41. Huberman, S. C., Larkin, J. M., McGaughey, A. J. H. & Amon, C. H. Disruption of superlattice phonons by interfacial mixing. *Phys. Rev. B* **88**, 155311 (2013).

42. Termentzidis, K., Chantrenne, P. & Keblinski, P. Nonequilibrium molecular dynamics simulation of the in-plane thermal conductivity of superlattices with rough interfaces. *Phys. Rev. B* **79**, 214307 (2009).

43. Termentzidis, K., Merabia, S., Chantrenne, P. & Keblinski, P. 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 (2011).

44. Hsieh, W.-P., Lyons, A. S., Pop, E., Keblinski, P. & Cahill, D. G. Pressure tuning of the thermal conductance of weak interfaces. *Phys. Rev. B* **84**, 184107 (2011).

45. Shen, M., Evans, W. J., Cahill, D. & Keblinski, P. Bonding and pressure-tunable interfacial thermal conductance. *Phys. Rev. B* **84**, 195432 (2011).

46. Losego, M. D., Grady, M. E., Sottos, N. R., Cahill, D. G. & Braun, P. V. Effects of chemical bonding on heat transport across interfaces. *Nature Mater.* **11**, 502–506 (2012).

47. McGaughey, A. J. H. & Kaviani, M. Phonon transport in molecular dynamics simulations: formulation and thermal conductivity prediction. *Advances in Heat Transfer* **39**, 169–255 (2006).

48. Landry, E. S., Hussein, M. I. & McGaughey, A. J. H. Complex superlattice unit cell designs for reduced thermal conductivity. *Phys. Rev. B* **77**, 184302 (2008).

49. Monaco, G. & Mossa, S. Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale. *Proc. Natl. Acad. Sci. USA* **106**, 16907–16912 (2009).

50. Yang, L., Yang, N. & Li, B. Reduction of thermal conductivity by nanoscale 3d phononic crystal. *Sci. Rep.* **3**, 1143 (2012).

51. Li, D., Wu, Y., Fan, R., Yang, P. & Majumdar, A. Thermal conductivity of Si/SiGe superlattice nanowires. *Appl. Phys. Lett.* **83**, 3186–3188 (2003).
52. Mizuno, H., Mossa, S. & Barrat, J.-L. Acoustic excitations and elastic heterogeneities in disordered solids. *Proc. Natl. Acad. Sci. USA* **111**, 11949–11954 (2014).

53. Luckyanova, M. N. et al. Coherent phonon heat conduction in superlattices. *Science* **338**, 936–939 (2012).

54. Kim, W. et al. Thermal conductivity reduction and thermoelectric figure of merit increase by embedding nanoparticles in crystalline semiconductors. *Phys. Rev. Lett.* **96**, 045901 (2006).

55. Zhang, H. & Minnich, A. J. The best nanoparticle size distribution for minimum thermal conductivity. *Sci. Rep.* **5**, 8995 (2015).

56. Nika, D. L. et al. Reduction of lattice thermal conductivity in one-dimensional quantum-dot superlattices due to phonon filtering. *Phys. Rev. B* **84**, 165415 (2011).

57. Yang, L., Yang, N. & Li, B. Extreme low thermal conductivity in nanoscale 3D Si phononic crystal with spherical pores. *Nano Lett.* **14**, 1734–1738 (2014).

58. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. *J. Comput. Phys.* **117**, 1–19 (1995). Http://lammps.sandia.gov.

59. Turney, J. E., McGaughey, A. J. H. & Amon, C. H. Assessing the applicability of quantum corrections to classical thermal conductivity predictions. *Phys. Rev. B* **79**, 224305 (2009).

60. Bedoya-Martinez, O. N., Barrat, J.-L. & Rodney, D. Computation of the thermal conductivity using methods based on classical and quantum molecular dynamics. *Phys. Rev. B* **89**, 014303 (2014).

## 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.

## Author Contributions

H.M., S.M. and J.-L.B. designed research, performed research, analyzed data, and wrote the paper.

## Additional Information

Supplementary information accompanies this paper at http://www.nature.com/srep

Competing financial interests: The authors declare no competing financial interests.

How to cite this article: Mizuno, H. et al. Beating the amorphous limit in thermal conductivity by superlattices design. *Sci. Rep.* **5**, 14116; doi: 10.1038/srep14116 (2015).

<figure><img src="image_11.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/
52. Mizuno, H., Mossa, S. & Barrat, J.-L. Acoustic excitations and elastic heterogeneities in disordered solids. *Proc. Natl. Acad. Sci. USA* **111**, 11949–11954 (2014).

53. Luckyanova, M. N. et al. Coherent phonon heat conduction in superlattices. *Science* **338**, 936–939 (2012).

54. Kim, W. et al. Thermal conductivity reduction and thermoelectric figure of merit increase by embedding nanoparticles in crystalline semiconductors. *Phys. Rev. Lett.* **96**, 045901 (2006).

55. Zhang, H. & Minnich, A. J. The best nanoparticle size distribution for minimum thermal conductivity. *Sci. Rep.* **5**, 8995 (2015).

56. Nika, D. L. et al. Reduction of lattice thermal conductivity in one-dimensional quantum-dot superlattices due to phonon filtering. *Phys. Rev. B* **84**, 165415 (2011).

57. Yang, L., Yang, N. & Li, B. Extreme low thermal conductivity in nanoscale 3D Si phononic crystal with spherical pores. *Nano Lett.* **14**, 1734–1738 (2014).

58. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. *J. Comput. Phys.* **117**, 1–19 (1995). Http://lammps.sandia.gov.

59. Turney, J. E., McGaughey, A. J. H. & Amon, C. H. Assessing the applicability of quantum corrections to classical thermal conductivity predictions. *Phys. Rev. B* **79**, 224305 (2009).

60. Bedoya-Martinez, O. N., Barrat, J.-L. & Rodney, D. Computation of the thermal conductivity using methods based on classical and quantum molecular dynamics. *Phys. Rev. B* **89**, 014303 (2014).

## 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.

## Author Contributions

H.M., S.M. and J.-L.B. designed research, performed research, analyzed data, and wrote the paper.

## Additional Information

Supplementary information accompanies this paper at http://www.nature.com/srep

Competing financial interests: The authors declare no competing financial interests.

How to cite this article: Mizuno, H. et al. Beating the amorphous limit in thermal conductivity by superlattices design. *Sci. Rep.* **5**, 14116; doi: 10.1038/srep14116 (2015).

<figure><img src="image_12.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/