Finite-size effects on phonon-mediated thermal transport across Si–Ge interfaces: Spectral analysis and parameter optimization for molecular dynamics simulations

The rapid development of miniaturized electronic devices, high-efficiency energy systems, and advanced composite materials has highlighted the critical role of interfacial thermal management in modern technology.^[1–4] Interfaces, atomic-scale boundaries between materials or phases with distinct properties (such as lattice structures and chemical compositions), exhibit thermal resistance that can dominate overall heat transfer. As the density of interfaces increases in nanoscale systems, the interfacial thermal resistance (ITR) often exceeds the intrinsic thermal resistance of the bulk materials, significantly impacting applications such as heat dissipation in high-power devices, thermoelectric energy conversion, and the thermal stability of heterostructures.^[5–9] Early research on ITR mainly focused on microscopic characterization and empirical regulation.^[10–12] However, advances in nanotechnology, microfabrication, and characterization techniques have shifted the focus to atomic-scale mechanisms, including the roles of phonon/electron coupling, interfacial atomic arrangement, and non-equilibrium thermal transport.^[13–20]

Historically, ITR was often neglected in heat transfer analyses. Keesom^[21] first emphasized its importance in 1936, but systematic study began only after Kapitza’s seminal work on copper/liquid nitrogen interfaces,^[22] leading to the term “Kapitza resistance.” Today, ITR is investigated through complementary approaches, including experimental methods (e.g., time-domain thermoreflectance (TDTR),^[23] frequency-domain thermoreflectance (FDTR),^[24–26] and scanning thermal microscopy (SThM)^[27,28]), theoretical models (e.g., the acoustic mismatch model (AMM),^[29] and the diffusive mismatch model (DMM)^[30]), and computational methods (e.g., non-equilibrium molecular dynamics (NEMD),^[31] non-equilibrium Green’s function (NEGF),^[32–34] and the Boltzmann transport equation (BTE)^[35]).

The AMM, developed in the mid-20th century, was among the first theoretical frameworks for ITR.^[29] This model treats phonons as elastic waves, predicting thermal resistance based on the mismatch in acoustic impedance (sound velocity and density) between materials. It assumes purely specular phonon reflection and neglects inelastic scattering, making it accurate only at low temperatures. As a result, AMM predictions typically underestimate experimental values of ITR. To address these limitations, the DMM was proposed.^[30] This method assumes complete diffuse scattering of phonons at the interface, with transmission probabilities governed by the mismatch in vibrational density of states between materials. While this approach better describes high-temperature or rough interfaces, it tends to overestimate thermal resistance due to its assumption of full phonon randomization. In practice, the true ITR lies between the AMM and DMM predictions, as real interfaces exhibit a combination of elastic and inelastic phonon scattering mechanisms.

With the rapid development of computer technology, numerical methods have become essential for studying interfacial thermal transport. In particular, they can accurately consider the influence of interface microstructure on interfacial heat transfer. Among them, NEGF^[32–34] and lattice dynamics (LD)^[36] calculations can precisely quantify phonon transmission coefficients while accounting for atomic-scale interface details, thereby improving prediction accuracy and playing an important role in understanding and regulating ITR.^[37] However, the NEGF and LD methods can only handle small systems and find it difficult to capture anharmonic phonon scattering in practical situations.^[38] Molecular dynamics (MD) simulation,^[39] which naturally incorporates any order of anharmonicity, can simulate systems with millions of atoms. Consequently, it serves as one of the most important methods for studying ITR.^[40]

Among various MD simulation methods, NEMD is the most widely adopted approach for investigating ITR. In this method, thermal baths are applied at both ends of the simulation domain, and the ITR is calculated from the steady-state heat flow and temperature drop at the interface according to Fourier’s law. However, NEMD results are known to be sensitive to simulation parameters, including the thermostat coupling constant τ, thermal bath length L_bath, and system length L_sys. Previous studies have demonstrated significant variations in predicted ITR depending on these parameters. For example, different research groups have reported disparate values for Si/Ge ITR when using different system lengths.^[41,42] Comprehensive reviews have confirmed that thermostat type and system dimensions substantially influence NEMD results, with particular sensitivity to system length due to finite-size effects.^[43,44] Additional studies have shown that thermal bath length affects phonon absorption and reflection at interfaces,^[45] though the underlying mechanisms remain incompletely understood.

In this study, we present a systematic investigation of parameter effects in NEMD simulations of Si–Ge heterojunctions. To evaluate the influence of thermostat properties on ITR, we designed two modified thermostats beyond the standard defect-free (pristine) configuration. These consisted of a structure with 20% random atomic vacancies and a second structure featuring a periodic array of triangular prism-shaped voids. The atomic vacancies intensify phonon scattering, thereby promoting rapid phonon thermalization inside the thermostat. In contrast, the geometrically ordered voids create a phononic crystal structure that fundamentally alters the phonon dispersion. Consequently, these tailored thermostats can simulate the behavior of materials distinct from the main simulation domain. Our study reveals that L_bath and L_sys significantly impact the calculated ITR, which exhibits a linear relationship with the reciprocals of L_bath and L_sys. While the thermostat coupling constant τ has minimal effect on resistance values, excessively large values fail to accurately control the temperature of the thermal baths. Spectral thermal conductance analysis shows that L_bath and L_sys mainly affect low-frequency phonons, which maintain ballistic transport characteristics over long distances due to their long mean free paths. When L_bath is small, these phonons reflect back to the interface rather than being absorbed by the thermal bath, thereby reducing the overall heat flux.

We modeled a Si–Ge heterojunction system with diamond cubic crystal structures for both materials (Fig. 1(a)). Atomic interactions were described using the Tersoff potential optimized by Lindsay and Broido,^[46,47] which has proven accurate for calculating thermal properties in Si and Ge systems.^[18,48,49] Because there is a small lattice mismatch between Si and Ge, to simplify the interface structure while maintaining physical relevance, we employed the Si potential parameters for both materials, preserving only the mass difference (28 amu for Si vs. 72.6 amu for Ge) across the interface. This approximation is justified by the small lattice mismatch and similar potential functions between Si and Ge, an approach successfully implemented in previous studies.^[50]

All simulations were performed using LAMMPS with the Verlet integration algorithm and a timestep of 1 fs. The system dimensions were 4.3 nm in both cross-sectional directions (y and z) with periodic boundary conditions. Along the heat transport direction (x-axis), the total system length was L_x = 2L_fix + 2L_bath + L_sys, where L_fix = 1.1 nm represented fixed regions at both ends to prevent rigid-body motion (Fig. 1(a)), while L_bath and L_sys were variable parameters corresponding to the thermal bath and central system lengths, respectively. The Langevin thermostat^[51] was adopted to control the temperatures at the two ends of the simulation domain, as previous studies have shown that the Langevin thermostat performs better than other types in heat transport simulations.^[52,53] The hot thermal bath was set at the Si end with a temperature T + ΔT = 315 K, while the cold bath was at the Ge region with a temperature T – ΔT = 285 K. All systems were first equilibrated for 1 ns under NPT ensemble conditions (300 K, 0 Pa). Then the systems were coupled with the thermal baths for NEMD simulations. After an initial 1 ns equilibration period to reach a steady state, we conducted 5 ns production runs while recording the heat flux and temperature profiles.

From the steady-state temperature profile (Fig. 1(b)), we extracted the interfacial temperature drop ΔT_int. The ITR, R, was then calculated using Fourier’s law

where J represents the heat flux obtained from the NEMD simulations.

To analyze frequency-dependent phonon contributions to thermal conductance, the spectral decomposition of heat flow was carried out in this work. In NEMD simulations, the heat flux q_i→j(ω) between atoms i and j in the left (L) and right (R) regions of the interface (Fig. 1(a)) can be calculated as follows:^[54]

where v^iα and v^jβ are the Fourier transforms of the velocities of atoms i and j in the α and β directions, respectively. Kijαβ represents the force constant between atom i in the α direction and atom j in the β direction, t_simu is the simulation time, and ω is the phonon frequency. The total interfacial spectral heat flux was obtained by summing over all atom pairs in the designated interfacial regions L and R (Fig. 1(a))

from which the spectral thermal conductance was calculated as

In the spectral decomposition calculation of heat flow, the sizes of the regions L and R should be larger than the potential cutoff radius to ensure that no interactions occur between atoms outside the designated regions.

To investigate how system dimensions influence the ITR, we first examined the relationship between the scattering region length L_sys and ITR. Considering that thermal bath materials may differ from the scattering region in practical systems (e.g., superlattices), three bath configurations were constructed (Fig. 1(c)): (1) defect-free Si/Ge, (2) baths containing 20% atomic vacancies, and (3) baths incorporating triangular prism-shaped voids (1 nm equilateral triangle base, occupying approximately 20% volume). The introduction of these two types of defects effectively modifies the phonon spectra,^[55] thereby simulating cases in which the bath and scattering regions have distinct material properties.

With L_bath fixed at 20 nm and a thermostat coupling constant of 1 ps, we explored the effect of L_sys on ITR. The results are shown in Fig. 2(a). The ITR of the Si/Ge heterojunction decreases as L_sys increases from 10 nm to 200 nm, eventually converging to similar values (2.29 m²·K·GW⁻¹) regardless of bath type. This convergence indicates that the bath material or structure becomes irrelevant once L_sys is sufficiently large. However, at small L_sys, defective baths yield higher resistance than pristine ones — a phenomenon associated with phonon ballistic transport. Introducing defects in the thermal baths creates additional interfaces between the baths and the scattering region, thereby altering the phonon spectra in these regions.

For large L_sys, phonons emitted from the thermal baths gradually convert to phonon modes corresponding to pristine Si or Ge. In contrast, when L_sys is small, phonons with long mean free paths (MFPs) do not have sufficient distance to fully convert to pristine Si or Ge modes through anharmonic scattering. These long-MFP phonons reach the Si/Ge interface without scattering, behaving as if the interface were directly formed between the two thermal baths. Consequently, the temperature jump experienced by such phonons equals the temperature difference between the two baths, effectively amplifying the calculated ITR through this interface-only scattering mechanism.

To further elucidate the relationship between ITR and L_sys, we plotted ITR as a function of 1/L_sys, as shown in Fig. 2(b). The results reveal a linear relationship between ITR and 1/L_sys, consistent with findings in AlN/GaN systems.^[44] This linear dependence allows extrapolation of ITR to infinite L_sys (1/L_sys = 0), yielding an intrinsic Si/Ge ITR of 2.285 m²·K·GW⁻¹. The reduction in ITR with increasing L_sys also stems from ballistic phonon effects. In short systems, ballistic phonons reflected by the thermostats can re-enter the simulation domain, decreasing the net heat flux and artificially inflating the measured ITR.^[44] As L_sys increases, a larger portion of phonons undergo thermalization before reaching the interface, reducing this spurious backward heat flow. The ITR stabilizes once L_sys exceeds the maximum phonon MFP in both materials.

Since the dependence of ITR on system length originates from the ratio between L_sys and the phonon MFP, we propose that the linear relationship between ITR and 1/L_sys is a universal phenomenon, independent of the specific material system. This universality is further supported by simulations employing the Stillinger–Weber (SW) potential,^[56] which offers a more realistic description of Si–Si, Si–Ge, and Ge–Ge interactions. Although the absolute ITR values differ between the Tersoff and SW potentials, both models exhibit the same decreasing trend with increasing L_sys and the same linear dependence on 1/L_sys, enabling consistent extrapolation to infinite system size (Figs. 2(a) and 2(b)).

Interestingly, the slope of the ITR vs. 1/L_sys plot is steeper for the SW potential than for the Tersoff potential. This difference likely arises because the SW potential yields a smaller intrinsic ITR. A lower ITR indicates weaker phonon scattering at the interface, resulting in a longer effective phonon mean free path (MFP) throughout the system. Consequently, finite-size effects persist over longer length scales, making the size dependence more pronounced.

Previous studies have demonstrated that the size of thermal baths in NEMD simulations is crucial for accurately determining thermal conductivity,^[57] yet the effect of thermostat size on ITR remains unclear. To examine the influence of thermostat length, we fixed L_sys at 100 nm and the thermostat coupling constant at 1 ps, while varying the thermostat length from 3 nm to 30 nm. As shown in Fig. 2(c), the resistance decreases with increasing L_bath and converges when L_bath ≥ 10 nm. This establishes a minimum bath length requirement (approximately 10 nm for Si/Ge) to obtain intrinsic resistance values. Remarkably, the resistance also scales linearly with 1/L_bath (Fig. 2(d)), allowing extrapolation to infinite bath length. The extrapolated ITR values shown in Fig. 2(d) are 2.288 m²·K·GW⁻¹, 2.307 m²·K·GW⁻¹, and 2.307 m²·K·GW⁻¹ for pristine thermostats, thermostats with atomic vacancies, and thermostats with triangular prism-shaped voids, respectively. These extrapolated values exhibit excellent consistency, with variations below 1% across different bath types. The dependence of ITR on L_bath originates from phonon absorption dynamics: short baths fail to fully absorb long-MFP phonons, leading to partial reflections that reduce the net heat flux and artificially increase the apparent resistance. The corresponding temperature and heat flux data are provided in Figs. S1–S3 of the supporting information.

The coupling parameter τ of the thermal bath plays a dual role in molecular dynamics simulations, influencing both temperature control accuracy and phonon dynamics within the bath region. Strong coupling (small τ values) ensures stricter temperature control but introduces larger perturbations to atomic vibrations. To systematically assess these effects, we examined τ values ranging from 0.1 ps to 10 ps while maintaining constant system dimensions (L_bath = 20 nm, L_sys = 100 nm) for the three bath configurations. The results in Fig. 3(a) show that ITR is remarkably insensitive to variations in τ. Across all three bath types, the calculated ITR varies by less than 0.02 m²·K·GW⁻¹, corresponding to relative differences below 3% throughout the tested τ range.

Temperature control analysis revealed important practical considerations for simulation design (Fig. 3(b)). The inset in Fig. 3(b) clearly shows that smaller τ values (≤ 1 ps) maintain precise temperature control in the bath regions, whereas larger τ values (particularly 10 ps) cause noticeable deviations from the target temperatures. Quantitative analysis of the temperature difference between hot and cold baths indicates a linear decrease with increasing τ, with pristine baths exhibiting a more rapid reduction than defective configurations. This behavior stems from two key factors: the intrinsically higher thermal conductivity of pristine baths facilitates more efficient temperature equilibration, while the effective ITR between defective baths and the scattering region reduces phonon transmission efficiency. When the temperature bias between the two thermostats remains reasonably small, the obtained ITR is independent of the temperature bias. This inherent property also explains why the thermostat coupling constant exerts a negligible influence on ITR.

Based on these findings, we recommend maintaining τ values at or below 5 ps to ensure adequate temperature control in thermal transport simulations. For an optimal balance between temperature regulation and minimal disturbance of phonon dynamics, moderate τ values in the range of 0.5–2 ps appear most suitable. These recommendations should be considered alongside bath structure effects, particularly when comparing systems with different defect configurations or material compositions. The demonstrated insensitivity of ITR to variations in τ provides valuable flexibility in selecting simulation parameters while preserving physical accuracy.

The preceding results established the significant influence of system length L_sys and bath length L_bath on ITR, underscoring the necessity of sufficiently large dimensions to achieve convergence. To elucidate the phonon-level mechanisms underlying these size effects, we performed spectral decomposition of heat flow and computed the interfacial spectral thermal conductance. To reduce statistical noise in the spectral decomposition, the spectral conductance was averaged over five independent simulations (Fig. S4). The analysis of L_sys dependence employed fixed bath parameters (L_bath = 20 nm and τ = 1 ps) across all three bath configurations. Here, we focused on spectral thermal conductance rather than resistance to avoid mathematical singularities at frequencies with zero heat flow. Figure 4 presents the spectral distribution of thermal conductance for all three thermal baths under various L_sys values. The results reveal a consistent trend: interfacial thermal conductance increases with L_sys and eventually converges, mirroring the behavior observed for total thermal conductance. To quantify these frequency-dependent variations, we calculated the relative change η(ω) at frequency ω between systems with L_sys = 10 nm and 100 nm

where G(ω, 10 nm) and G(ω, 100 nm) denote the thermal conductance at frequency ω for L_sys = 10 nm and 100 nm, respectively.

The spectral analysis demonstrates that the observed size dependence primarily originates from mid-to-low frequency phonons (ω < 6 THz). Across the three thermal bath systems, the thermal conductance of these phonons increases by approximately 20%–30% as L_sys increases from 10 nm to 100 nm. The physical origin of this behavior lies in the characteristic length scales of phonon transport: mid-to-low frequency phonons possess substantially longer mean free paths than their high-frequency counterparts, making their contribution to thermal transport more sensitive to system size variations.

Figure 5 shows the interfacial thermal conductance spectra for various L_bath values across the three thermal bath types. Similar to the system length dependence, the spectral thermal conductance increases with L_bath and eventually converges. However, bath length exerts a comparatively weaker influence on thermal conductance than system length, consistent with the trends observed in total thermal resistance. A detailed analysis of relative conductance variations reveals that L_bath predominantly affects low-frequency phonon modes. In the pristine bath system, phonons in the 2–4.5 THz range exhibit approximately 15% conductance enhancement as L_bath increases from 3 nm to 15 nm, while modes below 1.5 THz show a more substantial increase of about 70%. This frequency-dependent behavior is also observed in the triangular prism-shaped-void system. The atomic vacancy system displays reduced overall sensitivity to L_bath variations, yet it follows the same general trend, with maximum changes occurring in the low-frequency regime.

These results demonstrate that thermal bath length primarily influences interfacial thermal conductance through its impact on low-frequency phonons. The physical origin of this behavior lies in the characteristic length scales of phonon transport: low-frequency phonons possess significantly longer mean free paths than their high-frequency counterparts. When L_bath is insufficient, these long-wavelength phonons cannot be fully absorbed by the thermal bath and instead reflect back into the system. This reflection process effectively reduces their contribution to the net heat flux,^[44,58] thereby lowering the measured thermal conductance. The particularly strong effects observed for sub-1.5 THz phonons confirm their dominant role in establishing size-convergence criteria for thermal bath dimensions in molecular dynamics simulations.

In this study, we systematically investigated the effects of key simulation parameters — system length L_sys, thermal bath length L_bath, and thermostat coupling strength τ — on the ITR of Si/Ge heterojunctions using NEMD simulations. Our results demonstrate that both L_sys and L_bath significantly influence the calculated ITR, with resistance values decreasing as these parameters increase before eventually converging. In contrast, the thermostat coupling parameter τ shows negligible impact on ITR, although excessively large values (τ > 5 ps) compromise temperature control accuracy. We observed linear relationships between ITR and both L_sys and L_bath, providing a robust extrapolation method for determining bulk-like ITR from finite-size simulations. Spectral analysis reveals that these size effects primarily originate from mid-to-low-frequency phonons (< 6 THz), whose long mean free paths make their transport particularly sensitive to system dimensions. This work establishes practical guidelines for NEMD parameter selection in interfacial thermal transport studies while advancing fundamental understanding of phonon-mediated heat transfer across material boundaries. The identified convergence criteria and extrapolation methods enable more reliable prediction of intrinsic ITR values from molecular dynamics simulations. Furthermore, the mechanistic insights into frequency-dependent phonon transport provide a foundation for targeted interface engineering strategies in thermal management applications.