Diffusion Quantum Monte Carlo Benchmark of Interlayer Binding and Charge Redistribution in Chemically Distinct Two-Dimensional Van Der Waals Bilayers

Introduction

The emergence of 2D materials has revolutionized condensed matter physics and materials science, offering unprecedented opportunities for both fundamental research and technological applications [1, 2, 3]. These materials, which include graphene, hexagonal boron nitride (h-BN), and transition metal dichalcogenides (TMDs), exhibit a range of novel electronic, optical, and mechanical properties that make them suitable for applications in nanoelectronics, optoelectronics, and quantum information science [4, 5, 6]. Beyond monolayers, layered structures introduce an additional degree of freedom in material design through interlayer interactions, stacking arrangements, and electronic band hybridization, leading to tunable properties such as superconductivity, ferroelectricity, and topological states [7, 8]. Interlayer sliding and stacking can break inversion symmetry and produce switchable out-of-plane polarization, as in sliding ferroelectric bilayers, while magnetic stacking and crystal symmetry likewise control exchange pathways in 2D magnets and in related phenomena such as altermagnetism. Accurate interlayer separations, binding curves, and charge redistribution are therefore prerequisites for modeling these stacking-tunable polar and magnetic responses.

The accurate theoretical characterization of layered materials remains a challenging task because several distinct physical ingredients coexist. Non-local correlation is the electron correlation that depends on spatially separated regions of the density, hence missing in the local and semilocal exchange–correlation functionals. van der Waals (vdW) interactions are the long-range, asymptotically decaying part of that non-local correlation and provides the majority of the attraction between weakly interacting layers. These interactions are introduced to standard density functional theory (DFT) functionals via the addition of pairwise interactions[9, 10] up to self-consistent non-local vdW functionals with minimal empirical input[11, 12, 13, 14]. More sophisticated, beyond-DFT, approaches based on random-phase approximation (RPA) provide an explicit treatment of long-range correlation effects [15]; however, RPA atomization energies for solids are known to be systematically overestimated by up to 15% [15], and its quantitative reliability for layered and weakly bonded 2D systems remains to be validated using other many-body methods. An additional challenge is interlayer charge transfer at short distances in certain bilayers. As a consequence, the accuracy of widely used methods for interlayer interactions in 2D bilayers is not well-established, and the absence of systematically controlled reference data hinders both methodological assessment and further theoretical development.

This has motivated growing interest in quantum Monte Carlo (QMC) applications, in particular diffusion Monte Carlo (DMC), which enables highly accurate total energy calculations in molecular and periodic systems through an explicit treatment of electron–electron correlations [16] with few approximations. By describing intra- and interlayer interactions on an equal footing without empirical correction, DMC provides a natural framework for establishing a many-body benchmark for interlayer binding in layered 2D materials [17]. Improvements in the methods, their implementation, and computational resources, now enables the relatively small binding energies to be well resolved.

To move beyond isolated case studies and expose systematic trends, we selected 10 bilayer systems that span the main families of layered 2D materials while remaining accessible to a consistent high-throughput workflow based on bulk crystal structures from Inorganic Crystal Structure Database (ICSD) [18]. Bilayer graphene and hexagonal boron nitride provide canonical, weakly bound $sp^2$-bonded benchmarks with extensive prior many-body and experimental data, including prior DMC studies of bilayer graphene [19, 20]. Six transition-metal dichalcogenides, $\mathrm{MoX}_2$ and $\mathrm{WX}_2$ with $\mathrm{X} = \mathrm{S}, \mathrm{Se},$ and $\mathrm{Te}$, cover technologically relevant semiconductors and introduce systematic variation in atomic mass, polarizability, and interlayer binding strength. Chromium trihalides, $\mathrm{CrI}_3$ and $\mathrm{CrBr}_3$, in both antiferromagnetic and ferromagnetic interlayer alignments, extend the benchmark to magnetically ordered, correlated insulators for which mean-field DFT plus empirical dispersion corrections are expected to be the least reliable.

Together, these systems probe a wide range of interlayer separations and binding energies providing a stringent and chemically diverse test bed for assessing whether DFT-based approaches can reproduce many-body reference interlayer energetics. Fixed-node DMC has become a standard reference for weak and strong interlayer coupling in 2D systems when DFT dispersion treatments disagree [17]. While bilayer graphene [19, 20], bilayer hexagonal boron-nitride [21] and their heterostructures [22] and bilayer $\mathrm{MoS}_2$ [23] were previously studied using DMC, the study of layered materials with DMC has been often limited to monolayers, layered bulk materials and their defects [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42].

In this work, we use fixed-node DMC to systematically benchmark interlayer physics in chemically distinct 2D bilayers. We benchmarked the binding energy and interlayer separation at its minima ($E_b$, $d_0$), out-of-plane breathing mode frequencies $\omega$ and the layer-resolved dipole moments extracted from plane-averaged charge densities $\rho(z)$. We aim to establish a reliable benchmark dataset that can aid in the validation of approximate computational methods. This database will not only serve as a reference for researchers exploring new bilayered materials but also provide insights into the fundamental interactions governing their stability and functionality. Tabulated energetics and the raw charge densities are available in the supplementary information [43] and in associated data repositories.

Methods

Utilized geometries

We use a top-down approach to generate layered structures from bulk materials in the ICSD [18]. Therefore, the stacking order in the bilayer is transferred from the stacking order in the bulk. Detailed information on the structures used and their stackings are provided in the Supporting Information (SI) [43]. The monolayer and bilayer geometries used for the interlayer binding curves are consistently obtained from the bulk experimental coordinates. For direct comparison, the same geometries were used for the DFT binding curves: in-plane lattice constants and intralayer atomic coordinates are taken from the bulk-derived structures and held fixed while the interlayer separation $d$ is scanned rigidly. For DFT calculations, interlayer separations were studied up to 8 Å, while in DMC they were studied up to 5 Å. The use of periodic boundary conditions in z-direction means that image interactions are possible. We used a vacuum layer of at least 12 Å to minimize the image interactions.

Calculated quantities

We defined the interlayer separation ($d$) of the bilayer as the minimum vertical distance between the atoms in each layer of the bilayer. Using this definition, we fitted

to the binding energies calculated in DFT and DMC. We have extensively tested different functional forms including exponentials and other polynomials (SI Fig. S2 and Table S4). We found that using the 6–12 Lennard-Jones potential produced reasonable fits ($R^2 > 0.9$) for more repulsive DFT functionals (PBE, $\mathrm{r}^2$SCAN and PBEsol with no dispersion correction), but failed ($R^2 < 0.7$) for attractive DFT functionals (LDA, DFT+$vdW$) and DMC. Additional terms were added to improve the quality of the fit. We find that using these four parameters yields observable quantities that are robust with respect to the choice of fitting function. The form above was used to extract all interlayer separations $d_0$, binding energies $E_b$, and out-of-plane breathing-mode frequencies, $\omega$, reported in this work. This functional form is related to earlier DMC fits for bilayer graphene [19], where $E_b'(d) = A_2d^{-4} + A_4d^{-8} + A_6d^{-12} + A_8d^{-16}$ was found to give reliable interlayer separations and phonon breathing-mode frequencies. In our tests, the (1) also produces high quality fits ($R^2 > 0.99$ for most materials and functionals; mean $R^2 \approx0.99$ for DMC). Previous DMC and RPA calculations on graphene predicted that for the long-range portion of the binding-energy curve the dominant terms have $d^{-\alpha}$ scaling with $\alpha< 4$, which agrees with our choice of fitting function [20, 44].

We calculated the binding energy as $E_b = E_{\mathrm{bilayer}} - 2 \times E_{\mathrm{monolayer}}$, where $E_{\mathrm{bilayer}}$ and $E_{\mathrm{monolayer}}$ are the total energies of the bilayer and the corresponding isolated monolayer, respectively. Tabulated values in the SI are reported per atom (meV/atom) [43]. For cross-material comparison in Figure 1(b) we further normalize $E_b$ by the in-plane cell area (meV/Å$^2$). To obtain phonon breathing modes, the second derivative of the fitted binding-energy curve, $k$, is evaluated at the equilibrium separation $d_0$. The reduced mass of the bilayer reduces to $\mu= \langle m \rangle/ 2$, where $\langle m \rangle$ is the average mass per atom, which can be written as $\langle m \rangle= (nw_A + mw_B)/(n + m)$ for $\mathrm{A}_n\mathrm{B}_m$. The breathing-mode frequency (spectroscopic wavenumber) is calculated as

where $c$ is the speed of light, so that $\omega$ is reported in cm$^{-1}$ for comparison with Raman experiments. In DMC, the mean and the uncertainty of all quantities are obtained through Monte Carlo bootstrapping. The equilibrium $d_0$, $E_b$, and $\omega$ summarized in Figure 1 are extracted from these fitted curves; representative examples are shown in Figure 2, with the full benchmark set in SI Fig. S1 [43].

Diffusion Monte Carlo calculations

Fixed-node many-body DMC calculations [16] were performed as implemented in QMCPACK v4.0.0 [45, 46] via the Nexus workflow automation package [47]. For the QMC calculations, single-determinant Slater-Jastrow type trial wavefunctions [48, 49] were generated using correlation-consistent effective core potential (ccECP) potentials [50, 51, 52, 53] via the Quantum ESPRESSO package v7.4.0 [54] using the Perdew–Burke–Ernzerhof (PBE) [55] exchange-correlation functional for all materials. We used T-moves approximation for variational estimation of non-local pseudopotential terms [56]. For magnetic Cr-halides, the DMC ground state was obtained from PBE+U trial wavefunctions [57] with the Hubbard correction applied to Cr $3d$ states. We used a Hubbard $U$ of 1.0 eV for CrI$_3$ and CrBr$_3$ which minimizes the fixed-node DMC energy in their monolayers [43]. The single particle orbitals were represented in a hybrid basis set [58] to reduce memory consumption. We used two and three body Jastrow terms [59] to optimize the trial wavefunction. These Jastrow parameters were optimized via mixed energy/reweighted variance minimization [60]. A DMC time step of 0.005 Ha$^{-1}$ was used to minimize imaginary-time-step errors, and DMC total energies were extrapolated to the infinite-size limit using linear extrapolations over $N^{-5/4}$ in the number of atoms $N$, as previously demonstrated to be appropriate for 2D materials [61]. We used supercells constructed from $2 \times2$ scalar matrices ($N \times N$) where $N$ is an integer in the 2–4 range acting upon the primitive cells of each material to extrapolate the finite size errors. For bilayer graphene and hBN we used supercells containing up to 256 atoms, 1024 electrons, for TMDs up to 96 atoms, 832 electrons and for Cr-trihalides we used up to 144 atoms and 1260 electrons. We used twist averaged boundary conditions to reduce residual one-body finite-size errors in the Brillouin zone. Regular twist mesh grids were converged in the primitive cells using DFT calculations, then commensurate twist grids were used for all supercells. In this dataset the resulting statistical errors are typically $\sim0.01$ Å on $d_0$, $\sim0.1$–$1$ meV/atom on $E_b$, and $1$–$5$ cm$^{-1}$ on $\omega$, comparable to or smaller than the DFT–DMC discrepancies discussed below.

Density Functional Theory Calculations

DFT calculations for the structural and energetic benchmarks were performed in VASP v5.4.4 [62, 63] using projector augmented-wave (PAW) potentials [64]. For charge density benchmarks, in order to be consistent with the DMC calculations, we used the same ccECP potentials [50, 51, 52, 53] and settings as for the DMC trial wavefunction generation, changing only the exchange–correlation functional and any dispersion correction. Geometries follow the rigid $d$ scan described above; only the exchange–correlation functional and dispersion correction are varied. We used LDA [65], PBE [55], PBEsol [66] and r$^2$SCAN [67, 68] functionals and D3 [9], TS [10], MBD [13], and rVV10 [14] vdW correction methods. D3, TS and MBD corrections were only paired with PBE, which otherwise lacks long range dispersion, while the rVV10 correction was only paired with r$^2$SCAN. A kinetic energy cutoff of 520 eV was used, with a reciprocal-grid density of 350 per Å$^{-3}$ of the Brillouin zone (BZ) (a $12 \times12 \times1$ grid for graphene).

Results and Discussion

We summarize the results for equilibrium interlayer separation $d_0$, binding energy $E_b$, and out-of-plane breathing-mode frequency $\omega$ in Figure 1 and Figure 2. All individual results are given in SI Tables S1–S3 [43]. For every quantity below, we also provide a comparison between DMC and high-level ab initio calculations and experimental measurements from the literature. We find that our DMC results are in good agreement with the prior high-level ab initio calculations. Care must be taken when comparing any of the theoretical results to experimental data due to, e.g., the use of fixed lattice constants and coordinates.

For each observable $X \in\{d_0, E_b, \omega\}$, we define per-bilayer deviations $\Delta X_i = X_i^{\mathrm{DFT}} - X_i^{\mathrm{DMC}}$ over ten bilayers ($N = 10$). The signed mean deviation (MD) is

; the mean absolute deviation (MAD) is

Values quoted in Results are from SI Tables S1–S3 unless noted [43].

Equilibrium interlayer separation

Semilocal functionals exhibit a systematic ordering of equilibrium interlayer separation $d_0$: LDA systematically underestimates by MD/MAD of $-0.26/0.26$ Å, PBE systematically overestimates by MD/MAD of $0.86/0.86$ Å, and PBEsol lies between LDA and PBE, tracking DMC more closely with MD/MAD of $0.10/0.16$ Å. We also find that PBEsol performs better for multi-atom-thick layers (e.g., TMDs) than for single-sheet graphene and h-BN, consistent with its design for slowly varying densities in solids and slabs [66]. The standard deviation of the signed errors $\Delta d_0$ remains similar for all three semilocal functionals ($\sim0.14$–$0.17$ Å; fig:1a). This can indicate that for PBEsol, while the systematic bias is cancelled out, the spread of the errors is still significant and does not indicate a significant improvement over LDA and PBE. While r$^2$SCAN also yields improved $d_0$ values, compared to LDA and PBE, it performs worse than PBEsol in terms of the MD error. However, the standard deviation of the signed errors in r$^2$SCAN is smaller ($\sim0.06$ Å) than that of PBEsol ($\sim0.16$ Å), as indicated with the sharper peak in fig:1a.

The average results in Figure 1 show that dispersion-corrected functionals markedly improve $d_0$ for all the DFT functionals they are paired with. Most importantly, r$^2$SCAN–rVV10 achieves an impressive agreement with DMC with MAD of 0.05 Å, while PBE-based dispersion functionals perform qualitatively similarly (MAD $\sim0.15$ Å). Dispersion corrections based on PBE significantly improve the mean and mean-absolute deviation error metrics; however, the standard deviation of the signed errors in all dispersion-corrected PBE functionals is only minimally improved compared to bare PBE, and remains larger than that of bare r$^2$SCAN. In comparison, the addition of the rVV10 correction to r$^2$SCAN significantly improves every error metric we investigated compared to bare r$^2$SCAN.

Our DMC results for $d_0$ are in good agreement with high-level ab initio calculations and experimental measurements. For bilayer graphene, we find $d_0 = 3.38(1)$ Å, which compares well with previous DMC of bilayer graphene [3.43 Å] [19] and bulk graphene [3.36 Å] [20]. RPA predicts similar $d_0$ values for bilayer graphene (3.39 Å [69]) and bulk graphene (3.34 Å [70]), close to the experimental bulk spacing of 3.36 Å [71]. Similar agreement is found for bilayer h-BN with $d_0 = 3.36(1)$ Å, consistent with previous DMC (3.25–3.50 Å [21]) and Møller–Plesset perturbation theory (3.34 Å [72]) and with the experimental bulk spacing of 3.35 Å [73]. For TMDs and Cr-halides, our DMC vdW gaps lie within $\sim0.1$ Å of experimental nearest-layer spacings [74, 75, 76, 77, ?, 79]. Direct comparison with RPA bilayer TMDs requires a definition caveat: He et al. report the metal–metal (or c/2) distance, 6.27 Å for AA′ MoS$_2$ [80], whereas our $d_0$ is the chalcogen–chalcogen vdW gap. Subtracting the S–S thickness of one MoS$_2$ layer in our ICSD-derived geometry [75] ($t = 3.130$ Å) converts that RPA value to a gap of 3.14 Å, in agreement with our 3.13(1) Å. Qualitative trends are also recovered: larger terminating atoms give larger $d_0$ in the TMD series, and ferromagnetic Cr-halide stacking has a larger DMC $d_0$ than antiferromagnetic

ends_mid=1 stacking.

Interlayer binding energy

Figure 2 shows representative $E_b(d)$ curves and fits used to extract the benchmark observables in Figure 1. Panel (a) for bilayer graphene shows a shallow semilocal minimum and deep dispersion-corrected wells that overbind relative to DMC at the fitted minimum. Panel (b) for $\mathrm{MoS}_2$ illustrates the contrasting curvature of the binding curves: PBE yields a very shallow well, whereas LDA and DMC produce deeper minima at shorter separations. Panels (c) and (d) for antiferromagnetic and ferromagnetic $\mathrm{CrI}_3$ show that magnetic stacking order shifts the DMC minimum to larger $d$ in the FM case, whereas semilocal and dispersion-corrected DFT curves fail to reproduce this ordering (curves for all ten bilayers plus two magnetic configurations are given in SI fig:S1 [43]).

The semilocal ordering of $d_0$ and $E_b$ can be understood from how each functional treats exchange–correlation in regions of weak interlayer density overlap. LDA, based on the uniform electron gas, underestimates exchange repulsion in inhomogeneous systems and favors interlayer density overlap, yielding the shortest $d_0$. PBE enhances exchange repulsion through gradient corrections designed for molecular atomization, suppressing interlayer attraction and producing the largest $d_0$ and the weakest binding. PBEsol partially restores the gradient expansion appropriate for slowly varying solid-state densities, reducing PBE’s excessive repulsion while remaining between LDA and PBE on both observables. Even where LDA produces short separations and deeper wells than PBE, the fitted $E_b$ of LDA (and of other semilocal functionals) still typically falls below DMC for TMDs and Cr trihalides because nonlocal correlation is missing (SI tbl:S2 [43]). Graphene and h-BN are partial exceptions where LDA $E_b$ can approach DMC by coincidence of overestimated overlap attraction and underestimated dispersion.

Comparison with experimental quantities must be made with care: detailed visual representations of the differences between cleavage, exfoliation and cohesive energy can be found in the Supplementary Information of Ref. [81]. While these quantities are usually of similar magnitude, these terms will differ when the interactions between second and further neighbor planes are significant, as in the case of graphene/graphene [82]. Bulk interlayer binding of graphene—the energy to separate the crystal fully into monolayers—was reported as 56(5) meV/atom in DMC [20] and 48 meV/atom in RPA [70]. Cleavage energy (the energy to separate two semi-infinite stacks) is typically $\sim15$–$20\%$ larger than the binding energy in model estimates [82, 83, ?]ref-84, while Wang et al. later measured $\sim 64$ meV/atom for AB-stacked graphene {cite:p}ref-85`. Exfoliation and related cohesion/adhesion esti-

mates for graphene—$43 \pm5$ [83], $52 \pm5$ [84], and $\sim54$ meV/atom from G/G adhesion [86]—fall in the same $\sim45$–$55$ meV/atom range as the bulk DMC/RPA binding energies above, supporting that multilayer cohesion is much stronger than freestanding-bilayer binding. Our bilayer graphene $E_b$ is $15.4(4)$ meV/atom, consistent with bilayer DMC [$17.7(9)$ meV/atom [19]]. For MoS$_2$, RPA bilayer binding in the AA′ stacking is $\sim27$ meV/atom [80], close to our $24.3(1)$ meV/atom. Although we expect the changes to be small, and not yet practical to perform for all the materials, a full relaxation of all the geometries in DMC could further improve the agreement with experiment.

Interlayer breathing modes

Out-of-plane breathing-mode frequencies $\omega$ probe the curvature of $E_b(d)$ at $d_0$ and therefore provide a complementary test to the location and depth of the minimum (Methods; Figure 1; SI tbl:S3 [43]). As for $d_0$, bare semilocal functionals systematically order $\omega$: LDA yields the stiffest wells, PBE the softest, and PBEsol/r$^2$SCAN intermediate values. PBE underestimates $\omega$ by tens of cm$^{-1}$ across the set (MAD $24$ cm$^{-1}$), consistent with its shallow binding curves (Figure 2). LDA, dispersion-corrected PBE, and r$^2$SCAN–rVV10 all yield similar $\omega$ frequencies with MD less than $\sim2.2$ and MAD less than $\sim6$ cm$^{-1}$. In comparison, PBE, PBEsol, and r$^2$SCAN have systematic errors with MD of $\sim-24$, $\sim-13$, and $\sim-12$ cm$^{-1}$, respectively.

Dispersion corrections restore much of the missing curvature: PBE-D3, PBE-MBD, and r$^2$SCAN–rVV10 reduce the MAD on $\omega$ to $\sim3$–$4$ cm$^{-1}$, comparable to LDA’s fortuitous agreement. Improvements are not uniform, however—PBE-TS overestimates graphene $\omega$ ($89$ vs. DMC $73$ cm$^{-1}$), reflecting the same overbinding seen in $E_b$. In Cr trihalides, large DFT–DMC gaps in well depth and shape propagate into $\omega$: semilocal functionals underestimate AFM CrBr$_3$ frequencies, and magnetic stacking (AFM vs. FM) changes DMC $\omega$ more than most dispersion-corrected functionals predict. Overall, $\omega$ is more challenging to predict and functionals that partially correct $d_0$ can still fail for curvature-dependent observables as chemistry or magnetism changes.

Where Raman data exist, DMC agrees reasonably with other high-level theories and experiments. For graphene, we predict $\omega= 73(1)$ cm$^{-1}$ compared to a previous DMC result of $83(7)$ cm$^{-1}$ [19] and experiment at $80(2)$ cm$^{-1}$ [87]. For MoS$_2$ our DMC value $\omega= 36(4)$ cm$^{-1}$ agrees well with experiment at $40$–$41$ cm$^{-1}$ [88, 89, 90]. WS$_2$ agrees within the DMC uncertainty ($31(5)$ vs. $33$ cm$^{-1}$ [88]), whereas MoSe$_2$, MoTe$_2$, and WSe$_2$ are several cm$^{-1}$ below the reported Raman values. We did not find reports of rigid-layer breathing frequencies for h-BN, WTe$_2$, or the Cr trihalides to compare with our DMC results. Ultralow-frequency Raman of few-layer h-BN maps the interlayer shear mode ($52.5$ cm$^{-1}$ in bulk; $46$ cm$^{-1}$ at

3L) but not the breathing branch, which is Raman silent in the bulk [91]; the graphene $80$ cm$^{-1}$ value was obtained from resonant overtones of the silent ZO′ mode [87], a route closed for insulating h-BN. Bulk $T_d$-WTe$_2$ shows an $A_1$ interlayer peak at $8$–$9$ cm$^{-1}$ that overlaps a shear mode [92, 93, 94] but not a rigid-layer analog of our $\omega= 21(1)$ cm$^{-1}$. Raman of CrI$_3$ and CrBr$_3$ reports intramolecular $A_g/E_g$ phonons from $\sim50$–$280$ cm$^{-1}$ [95, 96] and the magnetism-coupled $\sim129$ cm$^{-1}$ $A_g$ Davydov splitting in few-layer CrI$_3$ [97], well above the rigid-layer well.

Interlayer dipole moments

We quantify interlayer-induced charge redistribution through the layer-resolved dipole moment

where $\rho(z)$ is the plane-averaged total bilayer density, $z_{\mathrm{vac}}$ is the bilayer vacuum midpoint, and $z_{\mathrm{layer}}$ is an individual layer center. The vacuum-bounded integral is evaluated on continuous sub-intervals $[z_{\min},z_{\mathrm{vac}})$ and $[z_{\mathrm{vac}},z_{\max}]$ using linear interpolation of $\rho(z)$. Bilayers are centered in the simulation cell along the $z$-axis, and periodic boundary contributions are neglected. Figure 3 summarizes $|\mu_z^{\mathrm{layer}}|$ as a function of nominal interlayer separation.

At representative short separations ($d = 2.5$–$3.0$ Å), interlayer polarization is strongest in Cr-trihalides and weakest in graphene/h-BN stacks, with TMDs occupying an intermediate regime. Graphene and h-BN are grouped as single-sheet layers—each repeat unit is one atomic plane (graphene or monolayer h-BN)—whereas TMD and Cr-halide layers are multi-atomic-thick.

For example, at $d = 2.5$ Å using PBE, $|\mu_z^{\mathrm{layer}}| \approx0.73\ e$ Å (CrI$_3$), $0.40\ e$ Å (CrBr$_3$), $0.24\ e$ Å (MoTe$_2$), $0.10\ e$ Å (MoS$_2$), and $\approx0.07\ e$ Å (h-BN and graphene). CrI$_3$ exceeds CrBr$_3$ by nearly a factor of two at fixed $d$, consistent with the larger polarizability and more diffuse halogen $p$ states of iodine. Among TMDs, heavier chalcogens give larger $|\mu_z^{\mathrm{layer}}|$ (e.g. MoS$_2 <$ MoSe$_2 <$ MoTe$_2$ and WS$_2 <$ WSe$_2$) at $d = 2.5$ Å, paralleling the halogen trend in CrX$_3$. For CrX$_3$, all semilocal functionals and DMC agree within a few percent at short range; the material hierarchy is unchanged across these methods. Magnetic stacking (AFM vs FM) changes $|\mu_z^{\mathrm{layer}}|$ by less than 1% and does not affect the trends below [43].

The $|\mu_z^{\mathrm{layer}}|$ decreases monotonically with increasing $d$ from 2.5 to 5 Å for all materials, reflecting the expected reduction of interlayer charge overlap as the layers decouple. The systematic increase of interlayer polarization along Br $\to$ I (CrX$_3$) and S $\to$ Se $\to$ Te (MX$_2$) suggests that heavier, more polarizable anions enhance interlayer charge coupling at fixed $d$. If such redistribution accompanies stronger interlayer exchange pathways, it may contribute, together with other mechanisms, to the stronger interlayer magnetic coupling observed in CrI$_3$ relative to CrBr$_3$.

For CrX$_3$ at short separations ($d \lesssim3.5$ Å), PBE and DMC agree within a few percent (e.g. CrI$_3$ at 2.5 Å: 0.73 vs. 0.77 $e$ Å). At larger $d$, where PBE predicts $|\mu_z^{\mathrm{layer}}| \to0$ for all classes, DMC retains a long-ranged tail in both Cr-halides independent of the magnetic configuration: at 5.0 Å, $|\mu_z^{\mathrm{layer}}| \approx0.007$ $e$ Å (PBE) vs. 0.082 $e$ Å (DMC) for CrBr$_3$, and 0.014 vs. 0.134 $e$ Å for CrI$_3$, whereas representative TMD and single-sheet bilayers have $|\mu_z^{\mathrm{layer}}| \lesssim0.015$ $e$ Å in both methods. Thus, while both approaches predict rapid decay of interlayer coupling with separation, many-body correlations in DMC support weak but non-negligible interlayer charge redistribution throughout Cr-halide family.

Conclusion

We have constructed a DMC benchmark of interlayer binding for ten 2D bilayers spanning sp$^2$, TMD, and magnetic chromium-trihalide chemistries with two magnetic configurations each. From fitted $E_b(d)$ curves we report equilibrium separations, binding energies, and breathing-mode frequencies, and we use them to assess widely used DFT approximations. No single functional simultaneously matches DMC for $d_0$, $E_b$, and $\omega$ across the full set (Figure 1); exact numerical agreement is not expected, but the pattern of failures is systematic. As expected, functionals without any van der Waals correction perform poorly. Of the tests functionals that include these corrections, r$^2$SCAN–rVV10 gives the best overall agreement. The tabulated bilayer curves and charge densities provide a transferable reference for developing and validating approximate interlayer methods.

Acknowledgments

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility (ALCF), which is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357. This research also used resources of the Oak Ridge Leadership Computing Facility (OLCF), which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725.

Supporting Information

Detailed numerical tables for equilibrium interlayer separation (Table S1), binding energy (Table S2), and out-of-plane breathing-mode frequency (Table S3), the full set of binding-energy curves (Fig. S1), and fit-quality diagnostics (fig:S2, tbl:S4) are provided in the separate SI document (SI.pdf).

Data Availability

The tabulated data and binding curves in the SI, together with ab initio electron densities available at DOI 10.13139/OLCF/2565002, provide a transferable reference for developing and validating approximate interlayer methods.