{"id":"d552a9a2-ccef-4737-8493-e93d61c470bf","arxiv_id":"2607.16783","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"Lattice initialisation—not domain length—changes the simulated Kapitza resistance of Cu-graphene interfaces by nearly 2x, via residual strain and a proposed damped interfacial layer.","lead":"This paper runs molecular-dynamics heat-flow simulations across copper–graphene–copper interfaces and finds that the computed interfacial thermal resistance changes by almost a factor of two depending only on how the atomic lattice is set up. The result matters because many published simulations use one such setup without checking the other.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Untested lateral-cell-size confound: the factor-two Case I vs Case II comparison changes Lx×Ly by ~24% and the relaxation path, not just lattice parameters.","rationale":"I am not challenging the empirical reproducibility of the two ITC values. The concern is about attribution: the abstract's 'despite differences of only a few per cent in the lattice parameters' implies the comparison isolates lattice parameters, but the two preparation protocols also differ in lateral cell dimensions and relaxation history. Since the proposed damping mechanism is wrinkling, which depends on lateral cell size, this is a load-bearing confound. The reader's local-equilibrium concern is legitimate but secondary: the Kapitza jump is extrapolated from bulk linear fits, and the factor-two gap is far larger than the reported fit uncertainties, so a local non-equilibrium near the interface would have to shift Case I and Case II very differently to erase the factor two. The lateral-size test is cheaper and more decisive. I would keep the conditional verdict: the empirical result may stand, but the mechanistic attribution needs this control before acceptance as stated.","tokens_in":25561,"tokens_out":10182,"duration_ms":106343,"concrete_test":"Run the Case II initialisation in a 50×50 Ų cross-section (same as Case I), using the nearest commensurate number of Cu and graphene unit cells and applying the same NPT relaxation; compute ITC at L_z ≈ 300 Å. If ITC stays near 300 MW/m²K, the lateral-size confound is immaterial. If it rises toward Case I's ≈550 MW/m²K, the factor-two claim must be re-attributed to lateral cell size / relaxation path rather than the few-per-cent lattice-parameter difference. A complementary arm—Case I with the 54.3×57.3 Ų cell—would bound the effect in the other direction.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central factor-of-two result rests on a comparison that is not a controlled lattice-parameter sweep. Case I (Sec. II C) uses a fixed 50×50 Ų cell with a_Cu=3.60 Å and a_G=2.46 Å; Case II uses an NPT-relaxed cell with a_Cu=3.632 Å, a_G=2.422 Å and final dimensions 54.3×57.3 Ų. The Sec. III C 'comparable' pair therefore changes (i) lattice constants by a few per cent, (ii) the lateral periodic cell area by ~24%, (iii) the number of graphene/copper unit cells, and (iv) the relaxation history. Lateral size directly sets the longest wavelength of graphene flexural modes and the possibility of wrinkling—the very mechanism invoked for Case II's damping layer—so the factor two cannot be attributed solely to 'initialisation-induced strain and density' without controlling L_x,L_y. The paper's own conclusion notes the interface area was kept constant, yet it differs between the two cases. No simulation in the paper varies the cross-section for a fixed initialisation, making this a live confound. If the factor two survives constant 50×50 Ų cross-section, the claim is secure; if not, the headline should be re-attributed to lateral finite-size/relaxation effects.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"This paper uses direct non-equilibrium molecular dynamics (NEMD) to study thermal transport across Cu–graphene–Cu interfaces. It compares two lattice-initialisation strategies: Case I uses experimental lattice constants (a_Cu = 3.60 Å, a_G = 2.46 Å) in a fixed 50×50 Ų lateral cell, while Case II uses separately relaxed and NPT-equilibrated lattices (a_Cu ≈ 3.632 Å, a_G ≈ 2.422 Å) with a larger, equilibrated lateral cell. The authors also vary domain length L_z and use two boundary-temperature protocols (ΔT = const. and δT = const.). The implementation is validated against Zhu et al. (2022), reproducing ITC = 630 ± 7 MW/m²K versus 640.2 MW/m²K. The central finding is that ITC is almost a factor of two larger in Case I (≈ 552 MW/m²K for L_z > 250 Å) than in Case II (≈ 289 MW/m²K), despite modest differences in lattice parameters. The paper argues that this difference is controlled by initialisation-induced strain and density, not by domain length or boundary enforcement, and interprets the lower ITC of Case II via a 'damping boundary layer' associated with graphene wrinkling and broadened interfacial phonon spectra. The copper lattice conductivity is found to increase with domain length, consistent with phonon mean-free-path effects.","tokens_in":25929,"tokens_out":6471,"duration_ms":61246,"significance":"If the factor-of-two initialisation sensitivity is confirmed, this is an important methodological result for the NEMD community: it would show that conventional lattice-construction choices can dominate over domain length and thermostat effects in interfacial thermal conductance prediction, and that bulk spectral overlap alone is insufficient. The paper has clear strengths: a successful literature validation, transparent error propagation, detailed supplementary descriptions of the relaxation protocol, and a data-availability statement. No target quantity is fitted, so circularity is not a concern. However, the central Case I versus Case II comparison is not a controlled lattice-parameter sweep: it simultaneously changes lateral periodic-cell dimensions, the number of unit cells, and the relaxation history. The explanatory 'damping boundary layer' mechanism is inferred from diagnostics on the same simulations and is entangled with this confound. The result is significant and worth pursuing, but the causal attribution to initialisation-induced strain/density is not yet established.","major_comments":[{"comment":"The two 'comparable' systems differ in more than lattice parameters: Case I uses (Lx,Ly)=(50.0,50.0) Å, while Case II uses (54.3,57.3) Å, i.e. an area increase of ~25%, plus a longer Lz and a completely different relaxation history (NVT fixed-cell vs NPT anisotropic relaxation). The concluding remark that 'the interface area was kept constant' (Sec. IV) is contradicted by these reported dimensions. Since lateral cell size directly sets the graphene flexural-mode cutoff and wrinkling propensity—the very mechanism invoked for Case II—the factor-of-two difference cannot be unambiguously attributed to initialisation-induced strain and density. A controlled test with identical Lx,Ly for both initialisations, or a lateral-size sweep within each case, is required to support the title claim.","section":"III C, Figs. 3 and 10; Table S1; Sec. IV"},{"comment":"The 'damping boundary layer' interpretation is built from post-hoc diagnostics of the same simulations: broadened regional VDOS, broader Cu–Cu nearest-neighbour distributions, and a nonlinear temperature profile in Case II. These diagnostics are consistent with the lateral-size/relaxation confound identified above and do not by themselves establish a causal mechanism. Moreover, the bulk spectral-overlap values (0.23 vs 0.25 in z) are presented without any uncertainty; a 0.02 difference is small relative to the factor-of-two ITC change, so the claim that 'spectral overlap alone cannot predict ITC' is not yet quantitatively supported. The paper should either provide an independent test (e.g., fixed lateral area, or artificially suppressing wrinkling) or explicitly present the damping layer as a hypothesis requiring further validation.","section":"III C, Figs. 7, 8, 11"},{"comment":"The central attribution of ITC variations to atomic density/strain is supported only by visual alignment of the ITC and density curves. No quantitative correlation, regression, or uncertainty interval is given. The density changes are a few per cent, and the ITC scatter is large for Lz < 100 Å (Case I) and throughout Case II. A simple Pearson/Spearman correlation with confidence intervals would substantially strengthen the claim that density is the primary controlling parameter. This point is secondary to the lateral-size confound but bears on the paper's mechanistic conclusion.","section":"III B, Fig. 6"}],"minor_comments":[{"comment":"The manuscript text as received contains many garbled symbols (e.g., Å rendered as '8A', missing overlines and subscripts in several equations). Please correct these throughout, as they significantly impede readability.","section":"General typesetting"},{"comment":"The text states that σ_Tbar_k is the standard error of the time-averaged temperature, but the least-squares formula does not show a weight 1/σ². Clarify whether the fit is weighted or unweighted.","section":"II B, Eq. (6)"},{"comment":"The values a_G = 2.46 Å (Case I) and a_G = 2.422 Å (Case II) are from different thermodynamic states (unrelaxed experimental vs 300 K relaxed). Please state this explicitly and consider a main-text table comparing the two initialisations: lattice parameters, cell dimensions, number of atoms, and residual strain.","section":"II C and Table S1"},{"comment":"The legend uses 'Atom/Volume' and the axis label 'N=LxLyLz' is unclear. Use 'number density (Å^{-3})' and define it explicitly in the caption.","section":"Fig. 6"},{"comment":"The data-availability statement says the data are in 4TU.ResearchData but gives no persistent identifier or link. Please provide a DOI or URL.","section":"Data availability"},{"comment":"The velocity-distribution check in SI SIII is performed for 'arbitrary slabs' in the copper domains. It does not directly address the interfacial regions I–V used for the VDOS analysis. A statement noting that equipartition is assumed rather than verified in the interfacial regions would be appropriate.","section":"III C, local equilibrium"}],"recommendation":"major_revision","confidential_remarks":"The paper addresses a question of real practical importance for NEMD simulations. The main obstacle is the lateral-size/relaxation confound in the central comparison. This is fixable with additional simulations: if the factor-of-two survives a constant (Lx,Ly) comparison, the claim would be solid. The spectral-overlap argument also needs error bars. The typesetting corruption in the submitted text should be corrected, although it may be an artifact of the submission pipeline."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Dear colleague,\n\nThe headline is that this paper delivers a clear warning: in NEMD of Cu-graphene interfaces, two conventional ways of setting up the lattice give ITC values that differ by nearly a factor of two (about 550 vs 290 MW/m2K at Lz around 300 Å). That difference is larger than many physical trends in the literature, so it matters. But before the causal claim is accepted, one confound needs closing: Case I is built on a 50 x 50 Å2 cross-section, while Case II is NPT-relaxed to about 54.3 x 57.3 Å2. That is a ~24% change in lateral area, plus a different relaxation history. Lateral size sets the longest graphene flexural wavelengths, so the factor two cannot be cleanly assigned to few-per-cent lattice-constant shifts without a fixed-area comparison. The paper's own \"interface area kept constant\" statement covers within-case sweeps, not the Case I/II contrast.\n\nWhat is genuinely good: the implementation is validated against Zhu et al. (630 ± 7 vs 640.2 MW/m2K), the domain-length sweep is systematic, and the VDOS and structural diagnostics are a reasonable first attempt to say what changes. The finding that greater bulk phonon overlap coincides with higher ITC, contrary to simple mismatch-model intuition, is worth taking seriously as a caution about spectral-overlap arguments. The damping-layer mechanism is clearly labelled as a suggestion, supported by local VDOS broadening and Cu nearest-neighbour distributions; indirect, but honest.\n\nSoft spots beyond the confound: \"no domain-length dependence\" is asserted from a figure rather than tested; that is a minor fix. The local-equilibrium temperature extraction is standard, and the supplementary velocity-distribution check helps, but the thin-slab jump remains definition-sensitive. The force-field combination is a known limitation, acknowledged. No target quantity is fitted, so there is no circularity problem.\n\nWho this is for: anyone doing NEMD of graphene-metal interfaces, and anyone interpreting small ITC differences from simulation as physical. It deserves a serious referee. The empirical factor-two difference is likely real, but the mechanism attribution needs one controlled cross-section sweep before it is established.","headline":"Useful warning that NEMD Kapitza resistance is very sensitive to lattice initialisation, but the factor-two claim is confounded by lateral cell size and relaxation path differences between the two cases.","tokens_in":26425,"tokens_out":2743,"would_cite":true,"duration_ms":28798,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Two lattice initialization choices produce a factor-of-two change in simulated graphene–copper Kapitza resistance.","keywords":["non-equilibrium molecular dynamics","Kapitza resistance","graphene-copper interface","lattice initialisation","residual strain","phonon spectral overlap","thermal interface conductance","finite-size effects"],"falsifier":"A concrete test is to take the relaxed (Case II) equilibrated structure, artificially flatten the graphene layer while keeping all lattice constants and densities fixed, and re-run the NEMD simulation: if the interfacial conductance jumps back to roughly 550 MW/m²K, the wrinkling/damping-layer mechanism is confirmed; if it stays near 290 MW/m²K, the cause of the factor-two difference lies elsewhere, such as in the density or the copper boundary layer alone. An independent check would be to replace the equipartition temperature estimator with a spectral or velocity-distribution-based local temp","tokens_in":25478,"feed_emoji":"🔥","tokens_out":4830,"duration_ms":47084,"temperature":0.7,"pith_summary":"This paper tries to establish that the computed thermal resistance across a graphene–copper interface in non-equilibrium molecular dynamics is controlled primarily by the strain and atomic density set during lattice initialization, not by the simulation domain length or boundary enforcement. It shows that two standard initialization strategies—one using experimental lattice constants, the other relaxing the lattices to minimize strain—produce lattices that differ by only a few per cent, yet give interfacial conductances of about 552 versus 289 MW/m²K, a factor of two. The lower-strain, relaxed case exhibits graphene wrinkling and a disordered, spectrally broadened copper boundary layer near the interface, which the authors interpret as a damping layer that raises thermal resistance. They further show that the relaxed case has slightly higher phonon spectral overlap between graphene and copper, yet lower conductance, indicating that spectral overlap alone cannot predict interfacial heat transport. A reader should care because it means absolute interfacial conductance values for this widely studied composite depend sensitively on seemingly minor simulation construction choices.","feed_headline":"Lattice choice doubles simulated graphene–copper heat barrier","feed_subtitle":"Relaxing the starting crystal creates a disordered copper layer that doubles Kapitza resistance despite nearly identical lattice constants.","key_machinery":"The key machinery is the direct NEMD setup: a Cu–graphene–Cu sandwich periodic in the transverse directions, with fixed atoms at the z-boundaries and Langevin thermostats enforcing hot and cold temperatures; from the resulting steady-state temperature profile, the Kapitza resistance is defined as the temperature jump δT at the graphene plane, obtained by extrapolating linear fits of the copper temperature, divided by the heat flux. The comparison hinges on two lattice initialization protocols: Case I sets lattice constants from experimental values (a_Cu = 3.6 Å, a_G = 2.46 Å) with no relaxation, while Case II relaxes each lattice separately at 0 K and 300 K, then combines them with minimal m","core_discovery":"The central discovery is that the Kapitza resistance of a Cu–graphene–Cu interface extracted from NEMD simulations is not a robust physical quantity under conventional domain construction: two initialization protocols, differing by only about 1–3% in lattice parameters, produce interfacial thermal conductances of 557.3 ± 3 MW/m²K (Case I, experimental lattice constants) versus 297.8 ± 1 MW/m²K (Case II, relaxed, minimal-strain lattices). The authors attribute this factor-of-two difference to residual strain and the resulting atomic density: higher strain in Case I keeps the graphene layer flat and the adjacent copper more crystalline, while the lower-strain Case II develops a wrinkled graphe","pith_inferences":["If this factor-two sensitivity generalizes, comparisons between published NEMD values for metal–graphene interfacial conductance should be treated with caution unless the initialization and equilibration protocol is reported in detail; even the same force field may yield answers that differ by up to a factor of two.","The damping-layer interpretation suggests a testable extension: computing the local thermal conductivity profile or mode-resolved transmission across the interface as a function of graphene corrugation amplitude would directly link wrinkling to the added resistance, and could be probed by intentionally introducing controlled ripples in an otherwise flat sheet.","The result implies that strain engineering of graphene–metal interfaces has an additional channel: residual strain affects not only phonon frequency shifts but also the structural order of the adjacent metal, which may dominate the interface conductance; this could guide experiments that anneal composites to modify interfacial thermal transport.","Since Case I and Case II differ in average copper–copper separation by only about 0.1 Å, the mechanism is subtle; a sensitivity study varying the copper–carbon Lennard-Jones parameters could test whether the disordered boundary layer is force-field specific or robust across potentials."],"forward_implications":["Simulated Kapitza resistance for graphene–copper is not converged by enlarging the simulation domain alone: once strain-related density effects are controlled, no residual length dependence appears up to ~500 Å, implying that reported length trends in earlier studies may be confounded by initialization artifacts.","The two initialization strategies imply that absolute interfacial conductance values for this interface carry an inherent configuration uncertainty of roughly a factor of two (about 290–560 MW/m²K), which should be quantified when comparing NEMD results to experiments or across different simulation studies.","Because bulk phonon spectral overlap increases while conductance decreases, predictions based on acoustic or diffusive mismatch models using bulk vibrational densities of states are insufficient; interfacial structure and local disorder must be included.","The copper lattice thermal conductivity in the same simulations shows clear domain-length and temperature dependence consistent with phonon mean-free-path limitation, meaning the interface resistance and the adjacent bulk conductivity have different sensitivities to simulation setup.","NEMD results can carry small statistical error bars yet be systematically controlled by domain configuration choices, so error propagation should include setup uncertainties rather than only sampling noise."],"fun_headline_variants":["Lattice prep doubles graphene–copper heat barrier","Strain choice doubles interface heat resistance","Two lattice setups, double Kapitza resistance","Configuration, not length, doubles graphene-Cu heat","Relaxed lattice doubles Cu-graphene heat barrier"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The central claim rests on the assumption that the local temperature extracted from kinetic energy via equipartition, together with the linear Fourier extrapolation near the interface, remains valid in the thin, strongly driven region next to graphene; if the velocity distribution near the interface is not quasi-equilibrium, the inferred temperature jump and hence the factor-of-two difference in Kapitza resistance would be artifacts of the post-processing.","fun_headline_variants_meta":{"raw":{"variants":["Lattice prep doubles graphene–copper heat barrier","Strain choice doubles interface heat resistance","Two lattice setups, double Kapitza resistance","Configuration, not length, doubles graphene-Cu heat","Relaxed lattice doubles Cu-graphene heat barrier"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000231,"raw_usage":{"total_tokens":1369,"prompt_tokens":836,"completion_tokens":533,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":580,"completion_tokens_details":{"reasoning_tokens":461}},"tokens_in":580,"tokens_out":533,"duration_ms":5490,"temperature":1.0,"reasoning_tokens":461,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-01T19:56:29.437348+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"A concrete test is to take the relaxed (Case II) equilibrated structure, artificially flatten the graphene layer while keeping all lattice constants and densities fixed, and re-run the NEMD simulation: if the interfacial conductance jumps back to roughly 550 MW/m²K, the wrinkling/damping-layer mechanism is confirmed; if it stays near 290 MW/m²K, the cause of the factor-two difference lies elsewhere, such as in the density or the copper boundary layer alone. An independent check would be to replace the equipartition temperature estimator with a spectral or velocity-distribution-based local temp","supporting_citations":[],"review_version":1}