Enhancing thermal transport in multilayer structures: A molecular dynamics study on Lennard−Jones solids

Abstract

We investigate the thermal transport properties of three kinds of multilayer structures: a perfect superlattice (SL) structure, a quasi-periodic multilayer structure consisted of two superlattice (2SL) structures with different periods, and a random multilayer (RML) structure. Our simulation results show that there exists a large number of aperiodic multilayer structures that have effective thermal conductivity higher than that of the SL counterpart, showing enhancement ratio in the effective thermal conductivity up to 193%. Surprisingly, some RML structures also exhibit enhanced thermal transport than the SL counterpart even in the presence of phonon localization. The detailed analysis on the underlying mechanism reveals that such peculiar enhancement is caused by the synergistic effect of coherent and incoherent phonon transport, which can be tuned by the structural configuration. Combined with molecular dynamics simulations and the machine learning technique, we further reveal that the enhancement effect of the effective thermal conductivity by 2SL structure is more significant when the period of SL structure is close to the critical transition period between the coherent and incoherent phonon transport regimes. Our study proposes a novel strategy to enhance the thermal transport in multilayer structures by regulating the wave-particle duality of phonons via the structure optimization, which might provide valuable insights to the thermal management in devices with densely packed interfaces.

Graphical abstract

Keywords

multilayer structures / thermal conductivity / machine learning / molecular dynamics simulation / wave-particle duality of phonon

1 Introduction

Heat conduction in low-dimensional materials [1-5] has attracted recent interests for the applications in thermal management [6,7] and thermoelectrics [8,9]. Among them, periodic structures have been extensively investigated recently due to the wave-particle duality of phonons [10-20]. For instance, Luckyanova et al. [10] found that the experimentally measured thermal conductivity κ of GaAs/AlAs superlattice (SL) increases almost linearly with length. The experimental study by Ravichandran et al. [11] further reveals the transition from incoherent to coherent phonon transport in SL structure when the interface density increases. Moreover, due to the coherent interference between phonons in SL structure, the total-transmission and total-reflection of individual phonons can be achieved in the presence of interfaces [14].

As certain level of mixing or roughness at the interface can often be introduced during the synthesis process of superlattice structure in experiment, the impact of interfacial mixing or roughness on the thermal transport in superlattice structure has also been investigated [21-23]. For instance, Huberman et al. [21] found that roughness or species mixing at the interface of superlattice structure leads to the reduction of thermal conductivity by suppressing the coherent phonon transport. In contrast to the non-monotonic dependence on the period length observed in perfect superlattice structure [11], thermal conductivity of the superlattice with roughness increases monotonically with the period length [22].

The miniaturization of nano-electronic devices has faced severe heat dissipation issue [24] due to the significantly enhanced interface density. Moreover, thermal management is also extremely important for the battery systems [25]. For the multilayer structures under the same system length and interface density, the periodic arrangement of interfaces in perfect SL structure seems to be the best option for heat conduction, since a suppressed thermal transport has been observed in various studies when introducing randomness to the perfect periodic structures due to the emergence of phonon localization [23,26-30]. Surprisingly, very limited counter-examples have been reported in recent studies [31,32]. For instance, Wei et al. [31] found that κ of a few aperiodic nanoporous graphene structures is abnormally higher than that of the periodic counterparts, although their study was based on Boltzmann transport equation method which ignores the wave characteristics of phonons. Another example is the simulation work by Chakraborty et al. [32] that they found κ of the one-dimensional (1D) multilayer structures does not decrease monotonically with increasing degree of randomness. These studies suggest that there might exist other competing mechanisms besides Anderson localization for the abnormally enhanced κ. However, only a small portion of the design space has been explored in previous studies. Besides, whether the observed enhanced thermal transport in aperiodic structure compared to the SL structure is a general feature in multilayer structure, and more importantly, the physical mechanism behind the enhanced thermal transport in aperiodic structure remains unclear.

To answer this question, it requires the search for the optimized configuration with maximized thermal transport ability among various 1D multilayer structures under the condition of the same interface density, which is a very challenging task for the traditional trial-and-error research procedures due to a huge amount of configurations in the design space. Fortunately, the machine learning (ML) technique has been demonstrated as a powerful tool to identify the optimized candidate for achieving targeted material properties [33-42].

In this work, we systematically examine the possibility for enhancing thermal transport in 1D multilayer structures via the structure optimization. Based on the molecular dynamics (MD) simulations, we observe that there exists a large number of aperiodic multilayer structures (even with randomly distributed interfaces) whose effective κ is higher than that of SL counterpart, under the condition of same interface density. The physical mechanism responsible for this abnormal enhancement is discussed in detail. Furthermore, the ML technique is applied to explore this abnormal feature in the whole design space and generalize the conclusion. Our work proposes a novel strategy to achieve enhanced thermal transport in multilayer structures by redistributing the interfaces and utilizing the coherent phonon transport.

2 Methods

2.1 Multilayer structures

Our target is to identify the optimized multilayer structure with maximized thermal transport ability better than the perfect SL structure, given the same interface density. Due to the transition from incoherent to coherent phonon transport, κ of SL structures [11,26] will first decrease and then increase with the decrease of SL’s period length P₀, exhibiting a minimum thermal conductivity at a critical transition period length P_t. An intuitive guess is that a combination of two SL structures with periods in two different transport regimes might have higher thermal conductivity than the perfect SL structure at P_t. Therefore, we consider three kinds of multilayer structures in this study: a perfect SL structure [Fig.1(a)], a combination of two SL structures (referred as 2SL in this work) with different periods [Fig.1(b)], and a random multilayer structure (RML) by introducing randomness to 2SL structure [Fig.1(c)]. For these three kinds of structures, we fix the total length L and the total number of interfaces so that the interface density is the same for all structures.

Fig.1 Schematic graphs for the three types of multilayer structures and NEMD setup. The dark and light regions represent material A and B, respectively. (a) SL structure with period length of P0 and number of layers of N0. (b) 2SL structure combined by two SL structures with different period lengths (P1 and P2) and number of layers (N1 and N2). R1 is the ratio of the left region. (c) RML structure. (d) Schematic setup for the NEMD simulations.

Full size|PPT slide

To build the SL structure, two conceptual materials A and B are considered and both of them have the conventional argon face-centered-cubic lattice structure with the same lattice constant of 5.23 Å. Here we define the unit cell (UC) length of argon as the length unit (UC = 5.23 Å). The value of P_t can be influenced by mass mismatch and the interatomic potential parameters between materials A and B. In this work, the masses of material A and B are equal to that of argon and the only difference between them is the interatomic potential parameters. The detailed selection basis for the parameters can be found in supporting information (SI) (see Figs. S5 and S6 in SI). All structures are constructed by stacking UC of material A and B along the [100] direction. The SL structure [Fig.1(a)] is constructed by alternately stacking material A and B with the same thickness, in which P₀is the period length and N₀ is the number of the layers (A or B). The 2SL structure [Fig.1(b)] is a combination of two SL structures with different period lengths (P₁ and P₂) and number of layers (N₁ and N₂). To guarantee the same interface density, we keep N₀=N₁+N₂. The RML structure is constructed based on the 2SL structure. Two layers are randomly selected from the left and right region of 2SL structure, respectively, and their positions are then swapped [22,30,32]. Note that these two layers selected in two regions must be of the same material to ensure that the total number of interfaces does not change after swapping. Such random swapping procedure is repeated S $\mathrm{\_}$num times to build the RML structure. In this work, S $\mathrm{\_}$num is set as 5000 to ensure the convergence of computed thermal transport properties of RML structure (see Section 2 in SI).

2.2 MD simulations

MD simulation is adopted in this work to study the thermal transport properties. Compared to the atomistic Green's function method and Boltzmann transport equation used in previous studies [31,34], MD simulations can consider both the full order of anharmonic phonon-phonon scatterings [42] and wave characteristics of phonons [43,44]. All non-equilibrium molecular dynamics (NEMD) simulations are performed using the LAMMPS [45]. The timestep is set as 1 fs. The interatomic interaction is described by the Lennard−Jones (LJ) potential $E(r_{ij})=4ϵ[{(\frac{\sigma }{r_{ij}})}^{12}-{(\frac{\sigma }{r_{ij}})}^6]$. The $ϵ$ for material A and B is set as 0.0416 eV and 0.1664 eV, respectively, which is chosen from a previous study [22] in order to ensure a large transition period between the coherent and incoherent transport regimes. The $\sigma$ is set as 0.34 nm for both material A and B, and the cutoff distance is 2.5 $\sigma$. Previous studies [26,30] reveal that the two-body LJ potential and many-body Tersoff potential actually capture essentially the same underlying physics regarding in the coherent and incoherent phonon transport in superlattice and disordered structures. Therefore, the merit of LJ potential is that it can capture the underlying physics with a simplified model. In this study, we focus on this LJ model system to study the control mechanism of κ for multilayer structures. The same strategy has been widely used in literature studies [22,30,32,46-48].

After constructing the structure, the whole system is first relaxed in the isothermal-isobaric ensemble with the periodic boundary condition applied in all directions. Then, NEMD simulations are performed according to the setup shown in Fig.1(d). At both ends of the simulation domain, a boundary layer of 2 nm in thickness is frozen as the fixed boundary. The structure under study (SL, 2SL, or RML) is sandwiched between a hot heat bath with temperature T_h and a cold bath with temperature T_c. Atoms in the heat baths are chosen to have the same atom type as the atoms in the neighboring layer next to the heat baths to minimize the temperature difference between the heat bath and central conduction region. Nosé−Hoover heat bath with a length of 4 nm is used at each end of the simulation domain. The temperatures of two heat bath are set atT_h = 35 K and T_c = 25 K, unless specified otherwise. Finally, we perform NEMD simulations long enough (~10 ns) to obtain the steady state heat flux and temperature profile.

As all three types of structures considered in our study is not homogeneous, and temperature discontinuity exists at the various interfaces inside each structure. Strictly speaking, a homogeneous thermal conductivity cannot be defined from the Fourier’s law for these inhomogeneous systems, especially for 2SL and RML structure. In order to compare the heat transport ability of these inhomogeneous structures, we select the same length for all structures and define the effective thermal conductivity for the inhomogeneous structure as

(1)${\kappa }_{eff}={\displaystyle \frac{JL}{S\mathrm{\Delta }T}},$

where J is the steady state heat current computed from the energy change rate in the heat bath, L is the length of the central conduction region, S is the cross-sectional area, and $\mathrm{\Delta }T$ is the temperature difference between the hot bath and cold bath. The same approach has been widely used in the study on thermal transport properties of multilayer structures [30,32,49].

2.3 Transmission calculation

With the mode analysis [50-54], the spectral contribution to thermal transport can be obtained. In this regard, we calculate the frequency-resolved phonon transmission across an imaginary interface in each structure as [51,52]

(2)$\Gamma \left(\omega \right)={\displaystyle \frac{2}{k_B\mathrm{\Delta }T}}\mathrm{R}\mathrm{e}{\displaystyle \sum _{k\in l}{\displaystyle \sum _{j\in r}{\int }_{\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{d}\tau {\mathrm{e}}^{\mathrm{i}\omega \tau }⟨{\mathit{F}}_{jk}\left(\tau \right)\cdot {\mathit{v}}_j\left(0\right)⟩,}}$

where k_B is the Boltzmann constant, F_jk is the interatomic force between atom j on left side of the imaginary interface and atom k on the right side of the imaginary interface, and v_j is the velocity of atom j. This method has the advantage to take into account the anharmonic effect on the transmission coefficient [51,52], which is ignored in the Green’s function approach. Therefore, it has been widely used to study the thermal transport mechanism in various systems [47,55-57].

3 Results and discussion

3.1 Impact of different structures

From the P₀−κ relationship in pure SL (see Fig. S1 in SI), we can see that phonon transport transits from incoherent diffusive phonon transport dominant regime to the coherent phonon wave transport dominant regime as P₀ decreases. Most of the phonons behave as coherent waves when P₀ is smaller than the critical transition period P_t. In the coherent phonon dominant heat transport regime, phonons can travel coherently without losing phase information [58]. In this regime, κ in SL structure decreases with increasing P₀ because of the reduced group velocity and the increase of band gaps [59,60]. Meanwhile, coherent phonons have the unique characteristics that they can travel through the interfaces as a wave rather than being scattered by the interfaces [30]. In the incoherent phonon dominant heat transport regime, the interface scattering leads to the phase breaking and decrease of phonon relaxation time [38,58], and κ in SL structure decreases with decreasing P₀ due to the enhanced interface density and consequently the enhanced incoherent phonon scattering.

In this section, we study three forms of multilayer structures with the same total length (L = 1024 UC) and number of interfaces (N₀ = N₁+ N₂) through the traditional MD simulations. Note that size effect has been tested in SI (see Figs. S3 and S4 in SI). In our model, we choose P₀, P₁ and P₂ as integer numbers of UC, and P₁ is always smaller than P₂. To maintain the same L and N₀, the following relationship must be satisfied: P₁ < P₀ < P₂. The P_t is a signature of the transition from the particlelike (incoherent) regime to the wavelike (coherent) regime. In this work, we select the model with P_t = 8 UC and the corresponding parameters are listed in SI (see Fig. S8 in SI). TheP₀of SL in this work varies from 4−1024 UC. In the following part, we discuss three different cases: (i)P₀ = P_t, (ii) P₀ < P_t and (iii) P₀ > P_t.

First, we study the most special case that P₀ = P_t = 8 UC. To ensure the same interface density as SL structure, the parameters in 2SL structure must satisfy the following condition: P₁ < P_t <P₂. In this case, the left region of 2SL is dominated by coherent phonons and the right region of 2SL is dominated by incoherent phonons. Here we use ${\kappa }_l$ and ${\kappa }_r$ to denote the thermal conductivity for the left and right region, respectively. Since P₀ = P_t, both ${\kappa }_l$ and ${\kappa }_r$ in 2SL are larger than that of SL because of the relationship betweenP₀and ${\kappa }_{eff}$ (see Fig. S8 in SI). The simulation results shown in Fig.2(a) indeed confirms that ${\kappa }_{eff}$ of most 2SL structures is larger than that of the corresponding SL (dashed line), especially for the case with small P₁ value. The enhancement ratio, $\mathrm{\Delta }=\frac{{\kappa }_{eff}^{2SL/RML}-{\kappa }_{eff}^{SL}}{{\kappa }_{eff}^{SL}}$, can even reach up to 193%. For P₁ = 6 UC, ${\kappa }_{eff}$ of 2SL is only slightly lower than that of the corresponding SL when the ratio of the left region R₁ is small, which suggests that the interfacial thermal resistance between the left and right region is negligible compared to the thermal resistance of each region in 2SL structure.

Fig.2 The relationship between the ratio of the left region R1 in 2SL and ${\kappa }_{eff}$ of multilayer structures with that the P0 of the corresponding SL is (a) 8 UC, (b) 4 UC, (c) 16 UC and (d) 16 UC. TheP1 of 2SL in (a), (b), (c) is smaller than the transition period lengthPt = 8 UC, while P1 in (d) is larger thanPt. The dashed lines in each plot denote the thermal conductivity of corresponding SL with the same interface density. Here Th= 35 K and Tc= 25 K are used in the simulations.

Full size|PPT slide

According to the constraint that N₀ and L are constant, the increase of R₁ for each P₁ leads to the increase of N₁ and decrease of N₂, and consequently the increase of P₂, which means that the difference between P₁and P₂ is enhanced. Fig.2(a) shows that ${\kappa }_{eff}$ of 2SL structure increases monotonically with R₁for each P₁. With the increase of R₁ in the left region of 2SL structure, more coherent phonons form and participate in phonon transport. At the same time, increasing R₁ results in the decrease of N₂ in the right region, which consequently reduces the scattering of incoherent phonons. Therefore, the synergistic effect of coherent phonons and incoherent phonons results in the increase of ${\kappa }_{eff}$ with R₁for each P₁ in 2SL structure. In other words, increasing the length ratio of coherent dominated left region can enhance the thermal transport in 2SL structure.

Furthermore, Fig.2(a) also shows that a larger P₁of 2SL leads to a smaller enhancement in ${\kappa }_{eff}$ compared to that of SL structure. The constraint of constant L and N₀ means that increasing P₁is achieved by decreasing P₂. Therefore, both ${\kappa }_l$ and ${\kappa }_r$ in 2SL decrease with increasing P₁ (see Fig. S8 in SI). Combined with the negligible interfacial thermal resistance between the left region and right region, 2SL structure with a large P₁ owns a relatively low ${\kappa }_{eff}$.

RML can be considered as a random structural disturbance, which can destroy the wave nature (coherence) of the phonons and cause the phonon localization phenomenon [22, 23, 26, 32, 38, 61]. The phonon localization effect will inhibit the phonon transport [28,62] and consequently cause the reduction of thermal conductivity. This point can be observed in Fig.2(a) that ${\kappa }_{eff}$ of RML structure (empty symbols) is obviously lower than that of the corresponding 2SL structure (solid symbols). Very interestingly, ${\kappa }_{eff}$ of certain RML structure (see for instance, P₁= 2 UC) is still larger than that of the corresponding SL structure, despite the localization induced reduction of thermal conductivity in RML. This peculiar behavior highlights the importance of coherent phonon transport for enhancing thermal transport compared to the SL structure even in the presence of randomness.

The second case is that P₀ < P_t, corresponding to the situation that the SL structure is in the coherent phonon dominant regime. This condition yields two possible situations: P₁<P₂<P_t or P₁ < P_t < P₂. For the situation P₁<P₂ < P_t, both the left and right regions in 2SL are dominated by coherent phonons. However, since there is only one such 2SL structure satisfying this condition (P₁ = 2 UC, N₁ = 256, P₂ = 6 UC, N₂ = 256) in our structure library, we do not consider this special case in our study, due to the rigor of statistical physics. Another situation is thatP₁<P_t < P₂. Here we consider SL with P₀ = 4 UC, and the corresponding 2SL and RML with P₁= 2 UC. As shown in Fig.2(b), ${\kappa }_{eff}$ of 2SL in this case also increases monotonically with the increase of R₁, since the same condition (P₁ < P_t < P₂) is satisfied as that in Fig.2(a). Moreover, ${\kappa }_{eff}$ of all 2SL structures in Fig.2(b) is larger than that of corresponding SL (dashed line), especially for structure with large R₁. The maximum enhancement ratio $\mathrm{\Delta }$ = 161%, which is smaller than that in Fig.2(a) ( $\mathrm{\Delta }$ = 193%). When random layers are introduced, ${\kappa }_{eff}$ of RML is reduced compared to that of 2SL, but there still exists a large number of RML structures (when R₁> 0.45) with ${\kappa }_{eff}$ larger than that of corresponding SL structure.

The last case is that P₀ > P_t, corresponding to the situation that the SL structure is in the incoherent phonon dominant regime. Two possible situations are P₁ < P_t< P₂or P_t < P₁ < P₂. Here we take samples with P₀= 16 UC as objects of study. When P₁ < P_t< P₂, we also observe in Fig.2(c) an increasing trend of ${\kappa }_{eff}$ in 2SL as R₁ increase for each P₁. The maximum enhancement ratio of ${\kappa }_{eff}$ relative to SL can reach by up to $\mathrm{\Delta }$ = 171%, which is smaller than that in Fig.2(a) ( $\mathrm{\Delta }$ = 193%). This is because there are fewer coherent phonons in 2SL for a largerP₀. Consequently, the reduction of ${\kappa }_{eff}$ from 2SL to RML structure in Fig.2(c) is smaller than that of Fig.2(a).

When the situation turns to that P_t < P₁<P₂, both the left and right regions in 2SL are dominated by incoherent phonons. In this configuration, there are still some 2SL structures with ${\kappa }_{eff}$ higher than SL counterpart (see Section 3 in SI). However, the maximum enhancement ratio in Fig.2(d) ( $\mathrm{\Delta }$ = 39%) is much smaller than that in Fig.2(c) ( $\mathrm{\Delta }$ = 171%). So far, we can conclude that the presence of coherent phonons is a key factor leading to the increase of ${\kappa }_{eff}$. Notice that the coherent phonon contribution in the incoherent regime might not be negligible when P₀ is close to P_t [26]. Therefore, when the randomness is introduced, ${\kappa }_{eff}$ of RML is still reduced compared to that of 2SL when P₀= 16 UC. When P₀is much larger than P_t, the difference in ${\kappa }_{eff}$ between RML and corresponding 2SL is very small due to the negligible coherent phonon transport (see Section 4 in SI).

To further understand the influence of structural configuration on thermal transport, we compute the phonon transmission across an imaginary interface for a typical set of SL, 2SL and RML structures. Here we consider P₀ = 8 UC for SL structure, and P₁ = 2 UC, N₁ = 255, P₂ = 1538 UC, N₂ = 2 for 2SL and RML structures. The imaginary interface in 2SL is between the left and right regions, while it is at approximately the interface near the center of the sample for SL and RML structures. The atomic velocity and force in the region with thickness of 1 nm on each side of the imaginary interface are recorded to compute the transmission coefficient. Although the location of the imaginary interface is different for different structures, this does not affect the result of the phonon transmission spectrum because the spectrum decomposition is insensitive to the exact position of the interface in the steady state [51].

Fig.3 shows that the phonon transmission in our system is dominated by phonons below 7 THz. There are two dominant peaks (around 3 THz and 4 THz) in the phonon transmission spectrum for all structures. The 2SL structure has the largest transmission coefficient for these two peaks, while the SL structure has the lowest transmission coefficient. This results are consistent with the independent NEMD calculations of ${\kappa }_{eff}$ for these three structures (SL: 3.86 W·m⁻¹·K⁻¹, 2SL: 11.30 W·m⁻¹·K⁻¹, and RML: 7.11 W·m⁻¹·K⁻¹). The phonon transmission calculation reveals that the introduction of coherent phonon transport in SL structure can notably promote the phonon transmission for the dominant peaks, leading to the enhanced ${\kappa }_{eff}$ in 2SL structure compared to that in SL structure.

Fig.3 The transmission spectra of a typical set of SL, 2SL and RML. The configurations of the three structures are set as follows: P0 = 8 UC for SL structure, and P1 = 2 UC, N1 = 255, P2 = 1538 UC, N2 = 2 for 2SL and RML structures. Here Th= 35 K and Tc= 25 K are used in the simulations.

Full size|PPT slide

Temperature can affect the thermal transport in crystalline materials through the anharmonic phonon-phonon interactions [63,64]. In order to further explore the validity of the above-mentioned enhancement of ${\kappa }_{eff}$ by aperiodic multilayer structures at different temperatures, we repeat the simulations for a typical set (P₀ = 8 UC and P₁ = 2 UC) of 2SL and RML structures at various temperatures. Fig.4 shows the calculation results of ${\kappa }_{eff}$ versus the ratio R₁ for this set of 2SL and RML structures at various temperatures. The ${\kappa }_{eff}$ of corresponding SL is 3.86 W·m⁻¹·K⁻¹ at 30 K, 3.54 W·m⁻¹·K⁻¹ at 50 K, and 3.06 W·m⁻¹·K⁻¹ at 80 K. A similar increasing trend of ${\kappa }_{eff}$ with the increase of R₁ is observed at different temperatures, which confirms that the enhancement of thermal conductivity by 2SL and RML structures is valid at different temperatures and the optimal structure with maximum enhancement ratio keeps the same (the large R₁ limit). Moreover, for a given R₁ ratio, the ${\kappa }_{eff}$ of 2SL structures decreases with increasing temperature, due to the enhanced anharmonic phonon scattering. Although the ${\kappa }_{eff}$ of corresponding SL is also reduced with increasing temperature, we find the maximum enhancement ratio $\mathrm{\Delta }=\frac{{\kappa }_{eff}^{2SL}-{\kappa }_{eff}^{SL}}{{\kappa }_{eff}^{SL}}$ for the optimal structure is actually smaller at higher temperatures (193% at 30 K, 157% at 50 K, 134% at 80 K), which suggests that such enhancement strategy by aperiodic structure works better at lower temperature due to the stronger coherent phonon transport.

Fig.4 The relationship between the ratio R1 and ${\kappa }_{eff}$ in a particular set of 2SL and RML multilayer structures at various temperatures. Here P0= 8 UC and P1= 2 UC are used in the simulations, and the temperature denotes the average temperature of two heat baths. The ${\kappa }_{eff}$ of corresponding SL is 3.86 W·m−1·K−1 at 30 K, 3.54 W·m−1·K−1 at 50 K, and 3.06 W·m−1·K−1 at 80 K. The maximum enhancement ratio $\mathrm{\Delta }$ for the optimal structure is 193% at 30 K, 157% at 50 K, and 134% at 80 K.

Full size|PPT slide

3.2 Predictions based on machine learning

So far, we have demonstrated for a limited number of structures that it is indeed possible to find alternative structures with higher ${\kappa }_{eff}$ than that of SL structure, given the same length and interface density. In the whole design space, there are 1105 and 16735 2SL structures corresponding to the situations that P₁ < 8 UC and P₁ > 8 UC, respectively (see Fig. S2 in SI). Accordingly, the number of the corresponding RML in the structure library is 17840. Therefore, there are 35680 structures in our structure library. Since the computational cost for each structure is about 160 core hours by NEMD simulations, such huge structure library makes it extremely challenging to explore the whole design space. In order to further verify whether the enhanced thermal transport is a general feature in the whole structure library, we resort to ML technique to predict ${\kappa }_{eff}$ of all structures in the whole design space.

The training set of ML needs to be large enough to ensure the accuracy of prediction and avoid the occurrence of overfitting, which can improve the generalization ability of ML predictions. To ensure enough 2SL structures, P_t should be sufficiently large, which can be achieved by selecting appropriate simulation parameters (see Section 1 in SI). Our targeted construction method enables us to consider the coherent phonon transport more effectively. Through optimal parameter selection (see Figs. S3−S7 in SI), we set the cross section as 3 UC × 3 UC andL = 1024 UC with a negligible 5% error tolerance for all structures in our calculations. The LJ parameters for Ar are ${ϵ}_{Ar}$ = 0.0104 eV and ${\sigma }_{Ar}$ = 0.34 nm, respectively [22]. The $ϵ$ of material A and B is 4 and 16 times of ${ϵ}_{Ar}$, and the $\sigma$ of material A and B is the same as that for argon. The cutoff radius is set as 2.5 ${\sigma }_{Ar}$. In all structures, the mass of materials A and B are both 40 g/mol, which is the mass of argon. The ambient temperature considered in our work is 30 K.

Prediction performance in ML depends not only on the completeness of the initial input data, but also on the selection of algorithms and descriptor dataset. Inspired by the recent successes of traditional ML algorithms [65-69], we examine various ML algorithms by comparing the coefficient of determination (R²), mean squared error (MSE), root mean squared error (RMSE). R², MSE and RMSE of Gradient Boosting Decision Tree (GBDT) [70] for the test dataset are 0.92, 0.39 and 0.63, which demonstate that GBDT is the best prediction algorithm (see Figs. S11 and S12 in SI). The descriptor dataset here consists of L, P₀, P₁, N₁, P₂, N₂, R₁, S $\mathrm{\_}$num, which can describe the various structures in detail (see Fig. S13 in SI). Therefore, the GBDT algorithm is used to predict ${\kappa }_{eff}$ of all structures, with detailed method described in SI (see Fig. S14 in SI).

To explore the general relationship between ${\kappa }_{eff}$ and structure arrangement, the ML predicted ${\kappa }_{eff}$ for all 2SL and RML structures with the configuration that P₁ < 8UC is plotted in Fig.5(a), while that with the configuration that P₁ > 8 UC is plotted in Fig.5(b). Fig.5 also shows that ${\kappa }_{eff}$ of RML is consistently lower than that of the corresponding 2SL structure, as a consequence of Anderson localization. Only when phonon coherence is weak, the introduction of randomness cannot significantly affect κ of nanostructures [13,30,32,38,49,71].

Fig.5 GBDT prediction of ${\kappa }_{eff}$ of 2SL and RML structures for (a) P1 < 8 UC, and (b) P1 > 8 UC. The structure id is arranged in order of increasing corresponding P0 for better visualization of the data.

Full size|PPT slide

To compare ${\kappa }_{eff}$ with SL structure, the configurations for the optimal 2SL structure with maximum enhancement ratio $\mathrm{\Delta }$ corresponding to each P₀ in Fig.5 are listed in Tab.1. It is noticeable that the optimal 2SL with a larger Δ in this situation all possess the characteristic that the difference between P₁ and P₂ is large. The importance analysis of GBDT also illustrates that these two factors (P₁ and P₂) are essential to determine ${\kappa }_{eff}$ (see Fig. S13 in SI). Meanwhile, we find that the deviation between P₀ and P_t (8 UC) diminishes the ability to improve ${\kappa }_{eff}$. With the increase of P₀, the maximum enhancement ratio $\mathrm{\Delta }$ in Tab.1 first increases and then decreases, which is due to the fact that the coherent phonon transport becomes weaker when the deviation between P₀ and P_t is larger. When P₀ = 1024 UC and P₁ < 8 UC, the maximum $\mathrm{\Delta }$ is negative because coherent phonon transport is negligible and the interface between the two regions in 2SL cannot be ignored. Except for this very large P₀ case, by comparing the impact on different region combinations in 2SL (P₁ < 8 UC and P₁ > 8 UC) on maximum $\mathrm{\Delta }$, we find that 2SL structures with the configuration P₁ < 8 UC consistently have a larger $\mathrm{\Delta }$ than those for the configuration P₁ > 8 UC to gain ideal improvement of ${\kappa }_{eff}$ (see Tab.1). This can also be explained by the fact that more coherent phonons are formed in 2SL structure with the configuration P₁ < 8 UC. Therefore, we can conclude that enhanced thermal transport relative to SL structure can be achieved in 2SL structure when the period of SL structure P₀ is near the critical transition period P_t. Since P_t typically has a very small value, our study provides an efficient way to promote the thermal transport in the multilayer structure by engineering the coherent phonon transport, despite the presence of very dense interfaces.

Tab.1 The machine learning predictions of ${\kappa }_{eff}$ in SL structure and the optimal 2SL structure with different configurations (P1 < 8 UC and P1 > 8 UC) versus the period P0of SL structure. Here, $\mathrm{\Delta }$ is the relative enhancement ratio defined as $\mathrm{\Delta }=\frac{{\kappa }_{eff}^{2SL}-{\kappa }_{eff}^{SL}}{{\kappa }_{eff}^{SL}}$.

P0(UC)	${\kappa }_{eff}$, SL (W·m−1·K−1)	P1 (UC)	P2(UC)	${\kappa }_{eff}$, 2SL (W·m−1·K−1)	Maximum $\mathrm{\Delta }$ (%)
4	4.25	2	1026	11.11	161
8	3.86	2	1538	11.30	193
16	4.32	2	1794	11.71	171
16	4.32	12	524	6.84	58
32	4.99	2	1922	11.65	134
32	4.99	10	1418	10.20	104
64	5.33	2	1986	11.56	117
64	5.33	16	1552	10.26	93
128	6.87	2	2018	11.55	68
128	6.87	18	1778	10.24	49
256	8.77	2	2034	11.55	32
256	8.77	16	1936	10.81	23
512	11.64	2	2042	11.84	2
512	11.64	368	560	12.30	6
1024	15.01	2	2046	12.39	−17
1024	15.01	1022	1026	15.01	0

4 Conclusion

To summarize, we have identified the aperiodic 2SL and its corresponding RML structures that have a higher ${\kappa }_{eff}$ than the SL counterpart by using both traditional MD simulations and a ML prediction model based on GBDT. For each SL period P₀, the enhancement ratio $\mathrm{\Delta }$ in ${\kappa }_{eff}$ increases monotonically with the increase of portion of left region with a smaller period P₁ in 2SL structure, due to the significant coherent phonon transport in the left region. Besides, the enhancement ratio $\mathrm{\Delta }$ can be further maximized with the decrease of P₁. When introducing randomness to the interface positions, ${\kappa }_{eff}$ of RML structure is reduced compared to that of 2SL structure, as a consequence of phonon localization. Surprisingly, there still exists RML structure with higher ${\kappa }_{eff}$ than that of the SL counterpart even in the presence of phonon localization, highlighting the importance of coherent phonon transport in enhancing the thermal transport. Furthermore, by exploring the whole design space via the machine learning technique, it is found that the enhancement effect of ${\kappa }_{eff}$ by 2SL structure is more significant when the period of SL structure P₀ is close to the critical transition period P_t between the coherent and incoherent phonon transport regimes. In addition, the 2SL structure with the configuration that P₁ < P_t is the best choice to gain maximum $\mathrm{\Delta }$, which can reach up to 193%. Our study provides a novel approach to enhance the thermal transport in multilayer structures by regulating the wave-particle duality of phonons via the structure optimization.

References

[1] Z.Zhang, Y.Ouyang, Y.Cheng, J.Chen, N.Li, G.Zhang. Size-dependent phononic thermal transport in low-dimensional nanomaterials. Phys. Rep. , 2020, 860 : 1

[2] W.Ren, Y.Ouyang, P.Jiang, C.Yu, J.He, J.Chen. The impact of interlayer rotation on thermal transport across graphene/hexagonal boron nitride van der Waals heterostructure. Nano Lett. , 2021, 21( 6): 2634

[3] C.Yu Y. Ouyang J.Chen, A perspective on the hydrodynamic phonon transport in two-dimensional materials, J. Appl. Phys. 130(1), 010902 ( 2021)

[4] G.Xie, D.Ding, G.Zhang. Phonon coherence and its effect on thermal conductivity of nanostructures. Adv. Phys. X , 2018, 3( 1): 1480417

[5] J.He, Y.Ouyang, C.Yu, P.Jiang, W.Ren, J.Chen. Lattice thermal conductivity of β₁₂ and χ₃ borophene. Chin. Phys. B , 2020, 29( 12): 126503

[6] A.L. Moore, L.Shi. Emerging challenges and materials for thermal management of electronics. Mater. Today , 2014, 17( 4): 163

[7] Y.Fu, J.Hansson, Y.Liu, S.Chen, A.Zehri, M.K. Samani, N.Wang, Y.Ni, Y.Zhang, Z.-B.Zhang. Graphene related materials for thermal management. 2D Mater. , 2019, 7 : 012001

[8] Y.Ouyang, Z.Zhang, D.Li, J.Chen, G.Zhang. Emerging theory, materials, and screening methods: New opportunities for promoting thermoelectric performance. Ann. Phys. , 2019, 531( 4): 1800437

[9] J.He, Y.Hu, D.Li, J.Chen. Ultra-low lattice thermal conductivity and promising thermoelectric figure of merit in borophene via chlorination. Nano Res. , 2022, 15( 4): 3804

[10] M.N. Luckyanova, J.Garg, K.Esfarjani, A.Jandl, M.T. Bulsara, A.J. Schmidt, A.J. Minnich, S.Chen, M.S. Dresselhaus, Z.Ren, E.A. Fitzgerald, G.Chen. Coherent phonon heat conduction in superlattices. Science , 2012, 338( 6109): 936

[11] J.Ravichandran, A.K. Yadav, R.Cheaito, P.B. Rossen, A.Soukiassian, S.Suresha, J.C. Duda, B.M. Foley, C.H. Lee, Y.Zhu, A.W. Lichtenberger, J.E. Moore, D.A. Muller, D.G. Schlom, P.E. Hopkins, A.Majumdar, R.Ramesh, M.A. Zurbuchen. Crossover from incoherent to coherent phonon scattering in epitaxial oxide superlattices. Nat. Mater. , 2014, 13( 2): 168

[12] T.Zhu, E.Ertekin. Phonon transport on two-dimensional graphene/boron nitride superlattices. Phys. Rev. B , 2014, 90( 19): 195209

[13] J.Maire, R.Anufriev, R.Yanagisawa, A.Ramiere, S.Volz, M.Nomura. Heat conduction tuning by wave nature of phonons. Sci. Adv. , 2017, 3( 8): e1700027

[14] P.Jiang, Y.Ouyang, W.Ren, C.Yu, J.He, J.Chen. Total-transmission and total-reflection of individual phonons in phononic crystal nanostructures. APL Mater. , 2021, 9( 4): 040703

[15] L.Yang, J.Chen, N.Yang, B.Li. Significant reduction of graphene thermal conductivity by phononic crystal structure. Int. J. Heat Mass Transf. , 2015, 91 : 428

[16] X.K. Chen, Z.X. Xie, W.X. Zhou, L.M. Tang, K.Q. Chen. Phonon wave interference in graphene and boron nitride superlattice. Appl. Phys. Lett. , 2016, 109( 2): 023101

[17] Z.Zhang, Y.Guo, M.Bescond, J.Chen, M.Nomura, S.Volz. Coherent thermal transport in nano-phononic crystals: An overview. APL Mater. , 2021, 9( 8): 081102

[18] Y.Zhou, X.Gong, B.Xu, M.Hu. First-principles and molecular dynamics study of thermoelectric transport properties of N-type silicon-based superlattice-nanocrystalline heterostructures. J. Appl. Phys. , 2017, 122( 8): 085105

[19] I.M. Felix, L.F. C. Pereira. Thermal conductivity of graphene−hBN superlattice ribbons. Sci. Rep. , 2018, 8( 1): 2737

[20] L.Razzaghi, F.Khoeini, A.Rajabpour, M.Khalkhali. Thermal transport in two-dimensional C₃N/C₂N superlattices: A molecular dynamics approach. Int. J. Heat Mass Transf. , 2021, 177 : 121561

[21] S.C. Huberman, J.M. Larkin, A.J. McGaughey, C.H. Amon. Disruption of superlattice phonons by interfacial mixing. Phys. Rev. B , 2013, 88( 15): 155311

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

[23] T.Juntunen, O.Vänskä, I.Tittonen. Anderson localization quenches thermal transport in aperiodic superlattices. Phys. Rev. Lett. , 2019, 122( 10): 105901

[24] E.Pop. Energy dissipation and transport in nanoscale devices. Nano Res. , 2010, 3( 3): 147

[25] C.Xiang, C.W. Wu, W.X. Zhou, G.Xie, G.Zhang. Thermal transport in lithium-ion battery: A micro perspective for thermal management. Front. Phys. , 2022, 17( 1): 13202

[26] S.Hu, Z.Zhang, P.Jiang, J.Chen, S.Volz, M.Nomura, B.Li. Randomness-induced phonon localization in graphene heat conduction. J. Phys. Chem. Lett. , 2018, 9( 14): 3959

[27] M.Luckyanova, J.Mendoza, H.Lu, B.Song, S.Huang, J.Zhou, M.Li, Y.Dong, H.Zhou, J.Garlow, L.Wu, B.J. Kirby, A.J. Grutter, A.A. Puretzky, Y.Zhu, M.S. Dresselhaus, A.Gossard, G.Chen. Phonon localization in heat conduction. Sci. Adv. , 2018, 4( 12): eaat9460

[28] S.Hu, Z.Zhang, P.Jiang, W.Ren, C.Yu, J.Shiomi, J.Chen. Disorder limits the coherent phonon transport in two-dimensional phononic crystal structures. Nanoscale , 2019, 11( 24): 11839

[29] T.Ma, C.-T.Lin, Y.Wang. The dimensionality effect on phonon localization in graphene/hexagonal boron nitride superlattices. 2D Mater. , 2020, 7 : 035029

[30] Y.Wang, H.Huang, X.Ruan. Decomposition of coherent and incoherent phonon conduction in superlattices and random multilayers. Phys. Rev. B , 2014, 90( 16): 165406

[31] H.Wei, H.Bao, X.Ruan. Genetic algorithm-driven discovery of unexpected thermal conductivity enhancement by disorder. Nano Energy , 2020, 71 : 104619

[32] P.Chakraborty, Y.Liu, T.Ma, X.Guo, L.Cao, R.Hu, Y.Wang. Quenching thermal transport in aperiodic superlattices: A molecular dynamics and machine learning study. ACS Appl. Mater. Interfaces , 2020, 12( 7): 8795

[33] A.Agrawal, A.Choudhary. Perspective: Materials informatics and big data: Realization of the “fourth paradigm” of science in materials science. APL Mater. , 2016, 4( 5): 053208

[34] S.Ju, T.Shiga, L.Feng, Z.Hou, K.Tsuda, J.Shiomi. Designing nanostructures for phonon transport via Bayesian optimization. Phys. Rev. X , 2017, 7( 2): 021024

[35] X.Wan, W.Feng, Y.Wang, H.Wang, X.Zhang, C.Deng, N.Yang. Materials discovery and properties prediction in thermal transport via materials informatics: A mini review. Nano Lett. , 2019, 19( 6): 3387

[36] Y.Ouyang, Z.Zhang, C.Yu, J.He, G.Yan, J.Chen. Accuracy of machine learning potential for predictions of multiple-target physical properties. Chin. Phys. Lett. , 2020, 37( 12): 126301

[37] S.Ju, S.Shimizu, J.Shiomi. Designing thermal functional materials by coupling thermal transport calculations and machine learning. J. Appl. Phys. , 2020, 128( 16): 161102

[38] P.R. Chowdhury, C.Reynolds, A.Garrett, T.Feng, S.P. Adiga, X.Ruan. Machine learning maximized Anderson localization of phonons in aperiodic superlattices. Nano Energy , 2020, 69 : 104428

[39] Y.Ouyang, C.Yu, G.Yan, J.Chen. Machine learning approach for the prediction and optimization of thermal transport properties. Front. Phys. , 2021, 16( 4): 43200

[40] L.Yang, X.Wan, D.Ma, Y.Jiang, N.Yang. Maximization and minimization of interfacial thermal conductance by modulating the mass distribution of the interlayer. Phys. Rev. B , 2021, 103( 15): 155305

[41] S.Arabha, Z.S. Aghbolagh, K.Ghorbani, S.M. Hatam-Lee, A.Rajabpour. Recent advances in lattice thermal conductivity calculation using machine-learning interatomic potentials. J. Appl. Phys. , 2021, 130( 21): 210903

[42] Y.Ouyang, C.Yu, J.He, P.Jiang, W.Ren, J.Chen. Accurate description of high-order phonon anharmonicity and lattice thermal conductivity from molecular dynamics simulations with machine learning potential. Phys. Rev. B , 2022, 105( 11): 115202

[43] Z.Zhang, Y.Guo, M.Bescond, J.Chen, M.Nomura, S.Volz. Generalized decay law for particlelike and wavelike thermal phonons. Phys. Rev. B , 2021, 103( 18): 184307

[44] Z.Zhang, Y.Guo, M.Bescond, J.Chen, M.Nomura, S.Volz. Heat conduction theory including phonon coherence. Phys. Rev. Lett. , 2022, 128( 1): 015901

[45] S.Plimpton. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. , 1995, 117( 1): 1

[46] P.Chakraborty, I.A. Chiu, T.Ma, Y.Wang. Complex temperature dependence of coherent and incoherent lattice thermal transport in superlattices. Nanotechnology , 2021, 32( 6): 065401

[47] A.Giri, J.L. Braun, P.E. Hopkins. Implications of interfacial bond strength on the spectral contributions to thermal boundary conductance across solid, liquid, and gas interfaces: A molecular dynamics study. J. Phys. Chem. C , 2016, 120( 43): 24847

[48] P.Chakraborty, L.Cao, Y.Wang. Ultralow lattice thermal conductivity of the random multilayer structure with lattice imperfections. Sci. Rep. , 2017, 7( 1): 8134

[49] B.Qiu, G.Chen, Z.Tian. Effects of aperiodicity and roughness on coherent heat conduction in superlattices. Nanoscale Microscale Thermophys. Eng. , 2015, 19( 4): 272

[50] Y.Zhou, X.Zhang, M.Hu. An excellent candidate for largely reducing interfacial thermal resistance: A nano-confined mass graded interface. Nanoscale , 2016, 8( 4): 1994

[51] K.Sääskilahti, J.Oksanen, S.Volz, J.Tulkki. Frequency-dependent phonon mean free path in carbon nanotubes from nonequilibrium molecular dynamics. Phys. Rev. B , 2015, 91( 11): 115426

[52] K.Sääskilahti, J.Oksanen, J.Tulkki, S.Volz. Role of anharmonic phonon scattering in the spectrally decomposed thermal conductance at planar interfaces. Phys. Rev. B , 2014, 90( 13): 134312

[53] Y.Zhou, X.Zhang, M.Hu. Quantitatively analyzing phonon spectral contribution of thermal conductivity based on nonequilibrium molecular dynamics simulations (I): From space Fourier transform. Phys. Rev. B , 2015, 92( 19): 195204

[54] Y.Zhou, M.Hu. Full quantification of frequency-dependent interfacial thermal conductance contributed by two- and three-phonon scattering processes from nonequilibrium molecular dynamics simulations. Phys. Rev. B , 2017, 95( 11): 115313

[55] Y.Ouyang, Z.Zhang, Q.Xi, P.Jiang, W.Ren, N.Li, J.Zhou, J.Chen. Effect of boundary chain folding on thermal conductivity of lamellar amorphous polyethylene. RSC Advances , 2019, 9( 57): 33549

[56] X.K. Chen, M.Pang, T.Chen, D.Du, K.Q. Chen. Thermal rectification in asymmetric graphene/hexagonal boron nitride van der Waals heterostructures. ACS Appl. Mater. Interfaces , 2020, 12( 13): 15517

[57] Y.Ma, Z.Zhang, J.Chen, K.Sääskilahti, S.Volz, J.Chen. Ordered water layers by interfacial charge decoration leading to an ultra-low Kapitza resistance between graphene and water. Carbon , 2018, 135 : 263

[58] G.Chen, Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons, Oxford University Press, 2005

[59] E.S. Landry, A.J. McGaughey. Effect of film thickness on the thermal resistance of confined semiconductor thin films. J. Appl. Phys. , 2010, 107( 1): 013521

[60] S.I. Tamura, Y.Tanaka, H.J. Maris. Phonon group velocity and thermal conduction in superlattices. Phys. Rev. B , 1999, 60( 4): 2627

[61] R.Hu, S.Iwamoto, L.Feng, S.Ju, S.Hu, M.Ohnishi, N.Nagai, K.Hirakawa, J.Shiomi. Machine-learning-optimized aperiodic superlattice minimizes coherent phonon heat conduction. Phys. Rev. X , 2020, 10( 2): 021050

[62] P.Jiang, S.Hu, Y.Ouyang, W.Ren, C.Yu, Z.Zhang, J.Chen. Remarkable thermal rectification in pristine and symmetric monolayer graphene enabled by asymmetric thermal contact. J. Appl. Phys. , 2020, 127( 23): 235101

[63] C.Yu, Y.Hu, J.He, S.Lu, D.Li, J.Chen. Strong four-phonon scattering in monolayer and hydrogenated bilayer BAs with horizontal mirror symmetry. Appl. Phys. Lett. , 2022, 120( 13): 132201

[64] S.Lu, W.Ren, J.He, C.Yu, P.Jiang, J.Chen. Enhancement of the lattice thermal conductivity of two-dimensional functionalized MXenes by inversion symmetry breaking. Phys. Rev. B , 2022, 105( 16): 165301

[65] T.Zhan, L.Fang, Y.Xu. Prediction of thermal boundary resistance by the machine learning method. Sci. Rep. , 2017, 7( 1): 7109

[66] Y.J. Wu, M.Sasaki, M.Goto, L.Fang, Y.Xu. Electrically conductive thermally insulating Bi–Si nanocomposites by interface design for thermal management. ACS Appl. Nano Mater. , 2018, 1( 7): 3355

[67] Y.J. Wu, L.Fang, Y.Xu. Predicting interfacial thermal resistance by machine learning. npj Comput. Mater. , 2019, 5 : 56

[68] Y.Liu, W.Hong, B.Cao. Machine learning for predicting thermodynamic properties of pure fluids and their mixtures. Energy , 2019, 188 : 116091

[69] Z.Hou, Y.Takagiwa, Y.Shinohara, Y.Xu, K.Tsuda. Machine-learning-assisted development and theoretical consideration for the Al₂Fe₃Si₃ thermoelectric material. ACS Appl. Mater. Interfaces , 2019, 11( 12): 11545

[70] J.H. Friedman. Greedy function approximation: A gradient boosting machine. Ann. Stat. , 2001, 29( 5): 1189

[71] M.R. Wagner, B.Graczykowski, J.S. Reparaz, A.El Sachat, M.Sledzinska, F.Alzina, C.M. Sotomayor Torres. Two-dimensional phononic crystals: Disorder matters. Nano Lett. , 2016, 16( 9): 5661

Declaration of interests

The authors declare no competing financial interest.

Electronic supplementary material

is available in the online version of this article at https://doi.org/10.1007/s11467-022-1170-5 and https://journal.hep.com.cn/fop/EN/10.1007/s11467-022-1170-5 and are accessible for authorized users.

Acknowledgements

This project was supported in part by the grants from the National Natural Science Foundation of China (Grant Nos. 12075168 and 11890703), the Science and Technology Commission of Shanghai Municipality (Grant Nos. 19ZR1478600 and 21JC1405600), and the Fundamental Research Funds for the Central Universities (Grant No. 22120220060).

RIGHTS & PERMISSIONS

2022 Higher Education Press