🔬 materials science

Lattice initialisation and finite-size effects of non-equilibrium molecular dynamics simulations for heat transfer across graphene-copper interfaces

This study demonstrates that Kapitza resistance at graphene-copper interfaces in non-equilibrium molecular dynamics simulations is highly sensitive to lattice initialization and residual strain, which can cause a two-fold variation in thermal resistance through strain-dependent phonon spectral shifts and the formation of damping boundary layers, rather than being primarily determined by domain size or simple phonon spectral overlap.

Original authors: L. A. van Goor, W. N. Edeling, D. Jafari, H. Lee, E. Luesink, W. W. Wits, A. V. Lyulin, B. J. Geurts

Published 2026-07-21

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: L. A. van Goor, W. N. Edeling, D. Jafari, H. Lee, E. Luesink, W. W. Wits, A. V. Lyulin, B. J. Geurts

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to build the ultimate heat-sink for a super-fast computer chip. You want to pull heat away from a tiny, scorching processor and dump it into a giant, cool block of metal. To do this, you might sandwich a super-thin, super-strong sheet of carbon (called graphene) between layers of copper. It sounds like a perfect plan: graphene is a heat-transfer superstar, and copper is a great conductor. But here's the catch: when you stick two different materials together, they don't always play nice. The heat often gets stuck at the boundary, like a traffic jam at a toll booth. This "traffic jam" is called thermal resistance, and if you can't figure out how to clear it, your computer chip might overheat and crash. Scientists use a powerful computer trick called Molecular Dynamics to simulate how atoms vibrate and move to see how heat flows across these boundaries. They want to know: does the size of our simulation matter? Does the way we set up the atoms at the start change the results? And most importantly, how do we get the heat to flow smoothly from the copper into the graphene?

This paper dives into that exact puzzle, simulating a sandwich of copper, a single layer of graphene, and more copper. The researchers were trying to solve a mystery: why do different computer simulations of the same setup sometimes give wildly different answers about how well heat crosses the interface? They focused on two main suspects. First, they looked at the "size" of the simulation box—does making the copper blocks longer or shorter change the heat flow? Second, and more importantly, they looked at how the atoms were "dressed" before the simulation started. Did the researchers stretch the atoms tight like a drum, or did they let them relax into a slightly wrinkled, natural state?

The team ran thousands of simulations using a method called Non-Equilibrium Molecular Dynamics (NEMD). Think of this as setting one side of the copper block on fire (hot) and the other side in an ice bath (cold), then watching how the heat tries to crawl through the graphene layer. They tested two very different ways of setting up the atomic lattices. In the first setup (Case I), they built the atoms using standard, textbook measurements for copper and graphene, which left the graphene slightly stretched and tight. In the second setup (Case II), they let the atoms relax and settle into their most comfortable, low-stress positions, which caused the graphene layer to wrinkle and ripple like a crumpled piece of paper.

Here is the big surprise they found: the way they set up the atoms at the very beginning changed the results by almost a factor of two! When the graphene was tight and stretched (Case I), the heat flowed across the interface much faster. When the graphene was relaxed and wrinkled (Case II), the heat got stuck, and the resistance doubled. This happened even though the difference in the atomic spacing was only a few percent. It turns out that the "wrinkled" graphene created a chaotic, disordered zone right next to the copper. This zone acted like a dampening layer, scattering the heat-carrying vibrations (phonons) and slowing them down. It's as if the tight graphene was a smooth highway for heat, while the wrinkled graphene was a bumpy, pothole-ridden road that made the heat cars crash and slow down.

The researchers also checked if the size of the simulation box mattered. They found that for the copper itself, the size did matter: shorter boxes meant the heat-carrying vibrations hit the walls more often, reducing the copper's ability to conduct heat. But for the interface between the copper and the graphene, the size of the box didn't really change the resistance, as long as the box was big enough. This suggests that the "traffic jam" at the interface is a local problem caused by the atomic arrangement, not a global problem caused by the size of the simulation.

One of the most interesting things the paper argues against is a common belief in the field: that heat transfer is mostly about how well the "music" of the copper atoms matches the "music" of the graphene atoms. Scientists often look at the overlap of their vibration frequencies to predict how well heat will flow. In this study, the "relaxed" (wrinkled) case actually had better frequency overlap than the "stretched" case. Yet, the heat flow was worse! This proves that just looking at the frequency match isn't enough. The local disorder and the "damping" effect of the wrinkled layer are the real bosses here.

In short, the paper suggests that if you want to predict how well heat moves across a copper-graphene interface, you can't just pick any starting point for your simulation. The tiny amount of leftover stress (strain) in the atoms, and whether the graphene is flat or wrinkled, makes a massive difference. While the size of the simulation box didn't change the interface results, the initial "outfit" of the atoms did. The authors conclude that we need to be very careful about how we set up these simulations to get reliable answers, because a small change in how we build the atomic model can lead to a huge change in the predicted heat flow. They didn't prove this is the only way nature works, but their simulations strongly suggest that the "wrinkles" and the resulting disorder are the hidden culprits behind the resistance.

Technical Summary: Lattice Initialisation and Finite-Size Effects in NEMD Simulations of Graphene–Copper Heat Transfer

Problem Statement
Non-Equilibrium Molecular Dynamics (NEMD) simulations are widely used to predict interfacial thermal conductance (ITC) and Kapitza resistance in copper–graphene composites, which are critical for electronic cooling applications. However, the reliability of these predictions is often compromised by uncertainties in computational setup choices. Specifically, the literature lacks a systematic understanding of how finite computational domain lengths and, more critically, the initialisation of atomic lattices (including residual strain and lattice constants) influence predicted thermal transport properties. While NEMD is known to be sensitive to domain size due to phonon mean-free-path limitations, the extent to which lattice initialisation strategies—often chosen arbitrarily to match experimental lattice parameters or to minimize strain—affect interfacial dynamics remains underexplored.

Methodology
The authors employed Direct NEMD simulations using the LAMMPS package to model heat transfer across a single-layer graphene interface embedded between copper domains. The study utilized a combination of force fields: the Embedded Atom Method (EAM) for copper, the Adaptive Intermolecular Reactive Bond Order (AIREBO) potential for carbon, and a Lennard-Jones potential for copper–carbon interactions.

To isolate specific sources of uncertainty, the study restricted variations to domain length ( $L_z$ ), lattice initialisation strategy, and temperature enforcement methods, while keeping the force field, thermostat, and cross-sectional area constant. Two distinct lattice initialisation strategies were compared:

Case I: Initialisation based on experimentally derived lattice parameters ( $a_{Cu} = 3.6$ Å, $a_G = 2.46$ Å), resulting in a system with residual strain and a flat graphene layer.
Case II: Initialisation based on a domain relaxation procedure (NPT equilibration) to minimize mismatch strain, resulting in slightly different lattice parameters ( $a_{Cu} \approx 3.61$ Å, $a_G \approx 2.41$ Å) and a visibly wrinkled graphene layer.

Two temperature enforcement schemes were tested to assess boundary artefacts:

Constant $\Delta T$ : Fixed temperature difference across the entire domain.
Constant $\delta T$ : Fixed temperature difference specifically across the interface, with extrapolated boundary temperatures.

Post-processing involved calculating local temperatures via the equipartition theorem, fitting linear temperature profiles to derive heat flux and Kapitza resistance ( $R_K$ ), and analyzing vibrational density of states (VDOS) to understand phonon spectra and overlap.

Key Contributions and Results

Sensitivity to Lattice Initialisation: The most significant finding is the extreme sensitivity of Kapitza resistance to lattice initialisation choices. Despite the two cases differing by only a few percent in lattice parameters and atomic density, Case I (higher residual strain, flat graphene) yielded an ITC of approximately 557 MW/m²K, while Case II (lower residual strain, wrinkled graphene) yielded an ITC of approximately 298 MW/m²K—a difference of nearly a factor of two.
Mechanism of Resistance: The study challenges the conventional view that spectral overlap alone dictates interfacial heat transfer. Although Case II exhibited slightly higher phonon-mode overlap between copper and graphene, it displayed significantly higher thermal resistance. The authors attribute this to the formation of a "damping boundary layer" in the lower-strain Case II. This layer is characterized by:
- Structural disorder and wrinkling of the graphene interface.
- Increased local structural disorder in adjacent copper atoms (broader nearest-neighbour distance distribution).
- Spectral broadening and increased scattering in the interfacial copper regions, which suppresses heat transfer.
Finite-Size Effects on ITC vs. Bulk Conductivity:
- Interfacial Conductance (ITC): Once strain-related effects (density variations) were accounted for, the ITC showed no significant dependence on domain length or the specific temperature enforcement method ( $\Delta T$ vs. $\delta T$ ). This suggests that ITC is robust to boundary artefacts and thermostat choices within the investigated range.
- Bulk Copper Conductivity: In contrast, the thermal conductivity of the copper domains adjacent to the interface exhibited clear dependence on both domain length and temperature. Smaller domains showed reduced conductivity due to boundary scattering limiting the phonon mean free path. Additionally, conductivity varied between hot and cold sides due to temperature-dependent scattering mechanisms (e.g., Umklapp scattering), consistent with established lattice thermal conductivity theory.
Phonon Spectra Analysis: The study demonstrated that bulk phonon spectra shifts (due to strain) are insufficient to explain the ITC differences. Instead, local interfacial dynamics, specifically the broadening of phonon modes and increased disorder in the immediate vicinity of the interface, are the dominant factors governing thermal resistance.

Significance and Claims
The paper claims that while NEMD simulations can yield precise numerical values for thermal properties, the reliability of the interfacial results is heavily contingent on the initialisation protocol. The authors argue that the "strain-dependent shift" in lattice constants is a primary driver of ITC variability, often overshadowing finite-size effects.

The study concludes that spectral overlap is an insufficient metric for predicting interfacial heat transfer in this system; the local structural and dynamic disorder at the interface plays a more critical role. Consequently, the authors assert that to reliably determine thermal transport properties for graphene–copper composites, researchers must carefully control and report lattice initialisation strategies and residual strain, rather than relying solely on domain size convergence or spectral overlap arguments. The findings highlight a need for further exploration of NEMD sensitivity to continuously varied input parameters to distinguish genuine interfacial physics from simulation artefacts.

More like this