Many-body quantum Monte Carlo study of 2D materials: cohesion and band gap in single-layer phosphorene

Quantum Monte Carlo (QMC) is applied to obtain the fundamental (quasiparticle) electronic band gap, $\Delta_f$, of a semiconducting two-dimensional (2D) phosphorene whose optical and electronic properties fill the void between graphene and 2D transition metal dichalcogenides. Similarly to other 2D materials, the electronic structure of phosphorene is strongly influenced by reduced screening, making it challenging to obtain reliable predictions by single-particle density functional methods. Advanced GW techniques, which include many-body effects as perturbative corrections, are hardly consistent with each other, predicting the band gap of phosphorene with a spread of almost 1 eV, from 1.6 to 2.4 eV. Our QMC results, from infinite periodic superlattices as well as from finite clusters, predict $\Delta_f$ to be about 2.4 eV, indicating that available GW results are systematically underestimating the gap. Using the recently uncovered universal scaling between the exciton binding energy and $\Delta_f$, we predict the optical gap of 1.75 eV that can be directly related to measurements even on encapsulated samples due to its robustness against dielectric environment. The QMC gaps are indeed consistent with recent experiments based on optical absorption and photoluminescence excitation spectroscopy. We also predict the cohesion of phosphorene to be only slightly smaller than that of the bulk crystal. Our investigations not only benchmark GW methods and experiments, but also open the field of 2D electronic structure to computationally intensive but highly predictive QMC methods which include many-body effects such as electronic correlations and van der Waals interactions explicitly.


I. INTRODUCTION
Two-dimensional (2D) materials have already revolutionized science, and have the potential to revolutionize technology due to their unique electronic, optical, thermal, spin, and magnetic properties. [1][2][3][4][5][6][7] Remarkably, 2D materials cover a wide range of electronic structures. The electronic properties range from metallic single atom layers of palladium and rhodium, 8 semimetallic graphene, 3 semiconducting transition metal dichalcogenides, 5 to insulating wide-gap h-BN. 6 Crucial for device applications are materials with a proper band gap. Layered black phosphorus (BP) features fundamental band gaps in the range of 0.3-2 eV, bridging semimetallic graphene and transition metal dichalcogenides. 9 This range is specifically important for optoelectronic, photovoltaic, photocatalytic, fiber optic telecommunications, and thermal imaging applications. 10 A single layer of black phosphorus-phosphorenecomprises sp 3 bonded phosphorus atoms forming an anisotropic puckered honeycomb lattice, see Figure 1. The three-fold bonding coordination implies that each phosphorus atom has a lone pair orbital which makes phosphorene reactive to air. 9 This oxidation degradation is eliminated by capping or encapsulating phosphorene with(in) an insulator.
Black phosphorus features a direct band gap at the Γ point, all the way from single-layer phosphorene up to bulk. However, and this makes phosphorene particularly attractive, there are several ways to manipulate the gap. (i) First, the gap varies with the number of TABLE I. Comparison of selected experimental and computed fundamental (∆ f ) and optical (∆o) gaps and exciton binding energies, in [eV], in unstrained single-layer phosphorene. STS, PL, PLES, and OA stand for scanning tunneling spectroscopy, photoluminiscence, photoluminescence excitation spectroscopy, and optical absorption, respectively. The experimental samples correspond to exfoliation onto different substrates (Si, SiO2, sapphire) and were studied either freshly prepared or capped by h-BN layer. G0 and W0 imply that the Green's function and screened Coulomb repulsion in the GW approach are calculated non-self-consistently. In GW calculations, the exciton binding energies, where quoted, were calculated in the GW-BSE (Bethe-Salpeter equation) approach. 11 GGA means generalized gradient approximation, H hybrid DFT functional. The present results appear in bold: DFT with GGA (PBE 12 ) and H functional (B3LYP 13 ) and DMC (upper two entries correspond to the periodic system with B3LYP and PBE nodal hypersurfaces, respectively, while the bottom entry to the cluster system with B3LYP treatment). The curly brace indicates compound value due to the cluster upper limit property, see the text. The star at the exciton binding energy means estimation as ∆ b = 0.27∆ f . 14 ∆ f ∆o  29 This strain engineering is predicted to affect the band gap by ≈ ±50%. 30 Finally, (iii) the gap of phosphorene is predicted to be strongly susceptible to the dielectric environment. 2,10,26 Assuming that the fundamental (quasiparticle) band gap ∆ f depends on the dielectric that protects it against degradation, can we infer from experiments some key characteristics about pristine phosphorene providing thus the benchmark for both experiment and theory? Two established facts make the answer positive. (a) The dielectric environment affects both the fundamental gap, as well as the exciton binding energy, ∆ b . 31 Remarkably, the difference, ∆ o = ∆ f − ∆ b , which is the optical gap, is essentially unaffected by the dielectric. 23,32-36 (b) Recently, a universal linear scaling between the exciton binding energy and the fundamental gap, ∆ b ≈ 0.27∆ f was predicted based on many examples from the 2D realm 14 (see also the predecessor 18 ), including phosphorene. Combining these two observations, (a) and (b), allows to estimate ∆ o from a calculation of ∆ f on a pristine 2D material, namely ∆ o ≈ 0.73∆ f , and compare with experimental ∆ o obtained from an encapsulated or capped sample.
Let us now turn to the existing experimental and theoretical state-of-the-art in determining the electronic gap(s) of single-layer phosphorene. Table I shows a rather comprehensive selection of measured and calculated fun-damental and optical gaps, as well as exciton binding energies for phosphorene. Photoluminescence is often contaminated by defect emission, as is also evidenced by the scatter of the measured values for ∆ o . Most reliable is optical absorption. A recent experiment 26 has reported ∆ o ≈ 1.73 eV, for phosphorene encapsulated in hBN. This value would lead to ∆ f ≈ 2.4 eV for pristine phosphorene, according to the above mentioned linear scaling. Similarly, it would suggest the exciton binding energy, ∆ b , of about 0.65 eV. Certainly both (a) and (b) observations are not exact, so the above estimates of ∆ f and ∆ b would carry a scatter of perhaps 10% or so. If we next look at the photoluminescence excitation spectroscopy data for phosophorene on silicon oxide, 37 the obtained fundamental gap is ∆ f ≈ 2.2 ± 0.1 eV. Considering the influence of the oxide, it is reasonable to deduce from this experiment that the fundamental gap of free-standing phosphorene is 0.1-0.2 eV greater, that is, 2.3-2.4 eV. 37 On the theory side, we see from Tab. I that singleparticle density functional methods predict, unsurprisingly, too low and strongly method-dependent values for ∆ f . Inclusion of GW corrections 11 is essential to bring the gaps closer to 2 eV. However, various GW approximants (G 0 W 0 or GW 0 ) give different values coming from different implementations, ranging from 1.6 eV to 2.4 eV. The largest value, 2.4 eV, which would be consistent with the aforementioned optical absorption experiments, is obtained by using a hybrid functional. 37 However, the same implementation predicts 0.6 eV band gap for bulk BP 37 (experimental value is 0.3 eV), making it clear that there is a limited predicting power from this calculation. As for the exciton binding energies, predictions (see Tab. I) range from 0.5 to 1 eV, again with little consensus in both theory and experiment. The experimental value of 0.9 eV 37 is likely affected by the optical edge of 1.3 eV of the emission peak (compare to 1.73 eV of the absorption experiment 26 ).
To obtain accurate bounds and reliable estimates of the band gap of phosphorene, we propose to employ the quantum Monte Carlo (QMC) method. QMC is an efficient, albeit computationally demanding 38,39 way to benchmark electronic structure calculations in condensed matter. [40][41][42] Indeed, this method has been applied to compute band gaps in three-dimensional systems, 40 clusters, 43 and nanoparticles. 44 In the 2D realm it was already used to obtain reference binding energies of 2D bilayers, 28,45 but thus far has not been systematically employed to obtain electronic structure parameters.
Here we report QMC calculations for the fundamental band gap of single-layer phosphorene. We use both the periodic lattice as well as supercell approaches, to demonstrate convergence and consistency. The accuracy of the ground-state properties is evidenced by calculating cohesion, which differs little from the bulk value. We stress that our QMC calculation is the full method, not relying on phenomenological interactions. In fact, we solve the Schrödinger equation for hundreds of electrons interacting mutually as well as with the lattice ions, to obtain the ground and excited states needed for the band gap calculation. Finally, we note that the knowledge of the band gap of phosphorene is important not only on its own right as a fundamental electronic quantity of a potentially technologically relevant material, but it is crucial also for building effective theories such as tight-binding 21 and k·p models. 46 We believe that our adaptation of QMC methods will open the way for this powerful technique to investigations of electronic structures of 2D systems, which are inherently prone to strong interactions and require careful considerations.

II. SIMULATION TECHNIQUES
The band gap ∆ f was determined from both extended and cluster approximants with lattice parameters fixed to the experimental values in the black phosphorus crystal. In fact, it was shown that the experimental lattice parameters agree with the lattice parameters determined by QMC methods within the error bars 28 . The gap ∆ f was extracted as the singlet-singlet vertical excitation energy. Here ∆ f ≈ E ss v = E s 1 − E s 0 , with E 0 and E 1 being, respectively, the ground-and the first excitedstates obtained by fixed-node QMC 38 not allowing any relaxation of the DFT nodal hypersurfaces due to the HOMO→LUMO electron excitation; no vibronic effects are included. The fixed-node approximation is the only fundamental approximation in the electronic structure QMC. 38 In the periodic setup the E 0 and E 1 were computed from DMC (diffusion Monte Carlo) energies in the fixednode approximation using the VMC (variational Monte Carlo) trial wave functions with the nodal hypersurfaces determined by two different sets of DFT orbitals: the generalized gradient approximation PBE 12 (PBE/DMC) and hybrid B3LYP 13 (B3LYP/DMC), at the Γ-point of the Brillouin zone, optimizing the short-range correlations of the Jastrow factor. 38 The consistency check using both PBE and B3LYP DFT nodal hypersurfaces was deemed important as at the DFT level the HOMO-LUMO gaps of the two DFT functionals differ by ≈1 eV, see Table I. The Yeh-Berkowitz 47 modification in the 3D Ewald summation technique for systems with a slab geometry that are periodic in two dimensions and have a finite length in the third dimension, was adopted. We cross checked in detail that this agrees with an alternative derivation of Ewald sums for 2D slab geometries. 48 In the cluster setup we used the B3LYP 13 (B3LYP/DMC) nodal hypersurfaces.
Finite-size scaling towards the thermodynamic limit was performed for a series of 1 × 1 to 6 × 6 series of L × L periodic approximants, see Fig. 1 a), and for 4 × 4 and 5 × 5 H-terminated cluster approximants, assuming a linear scaling with 1/N , where N = 4×L×L is the number of P atoms. Our approach corresponds to quasi-exact many-body treatment, to within the fixed-node approximation, of the 2D electron polarizability entering the equations for ∆ f . 14 The ground-state energy E 0 was also used to determine the cohesion energy. More details can be found in Supplementary Material (SM).

III. RESULTS AND DISCUSSION
Periodic supercells. Our finite-size diffusion Monte-Carlo (DMC) scaling study of ∆ f is shown in Fig. 2 a). Although in 3D periodic calculations typical QMC extrapolations from supercells to bulk behave as 1/N where N is the number of atoms in the supercell, in the 2D slab systems the issue is more complicated. 48 The strong periodic Coulomb interactions unscreened in the direction perpendicular to the layer make the renormalization of electron interactions huge and the electrostatic energies slowly convergent. Therefore, the supercells have to be sufficiently large to marginalize the periodicity effects and to allow for reliable finite-size extrapolations. To this end, we have carried out extensive calculations of ∆ f for the sequence of 1 × 1 up to 6 × 6 periodic supercell systems with two DFT nodal hypersurfaces. For sizes 3×3 or larger the linear scaling describes the results very well. For smaller supercells, the different components of the total energies, such as the potential energy, are too biased to make the trends transparent.
The infinite-size extrapolation from the extendedsystem calculations yields fixing the nodal hypersurfaces by the hybrid B3LYP 13 DFT functional. Recall that the corresponding DFT value is 1.7 eV, see Tab. I. If we now fix the nodal hypersurfaces by PBE 12 DFT single-particle orbitals, see Figure 2 a), we get the extrapolation to which overlaps with ∆ ext,1 f within the error bars. The PBE DFT gap is 0.8 eV. Remarkably, despite the significant difference of almost 1 eV in the DFT gaps, the QMC method that starts with the corresponding DFT wave functions (and fixing their nodal hypersurfaces), gives a consistent output for both! Surprisingly, the convergence of ∆ f in the periodic setting is determined mainly by the Hartree-Fock energy components, with the rest, including the correlation energies, converging much faster. For example, the ground-state correlation energy is essentially converged at the 3 × 3 supercell size, whereas the excited-state at a slightly larger size of 5 × 5. This provides a way of estimating the phosphorene properties from the more slowly converging Hartree-Fock results, simply by adding the corresponding correlation energies; see SM for more details. We also plot the HOMO/LUMO orbitals for the extended calculations at the Γ point in Fig. 2 b) and c), respectively. The orbital's distribution and bonding type can give insight into Coulomb finite size effects; see SM for a discussion.
Clusters. Due to intricacies of extrapolations of the periodic supercells, and for consistency, we have also calculated ∆ f in a finite cluster setting. The clusters corresponded to supercell sizes in the periodic method above, with terminations of unsaturated bonds with H atoms, see Fig. 1 c) and d). Finite-size scaling constructed from these cluster approximants, given in Fig. 2 a), extrapolates to reasonably close to ∆ ext f considering the fact that we compare periodic systems versus isolated clusters in vacuum and that these two models actually exhibit opposite trends as the functions of the system size. Due to finitesize confinement, the computed excitations in finite clusters as a rule overestimate the gaps. The DMC values for the 5 × 5 clusters clearly illustrate this, see Fig. 2 a). However, as explained in SM, the cluster gaps are bound to converge to the fundamental gap in the infinite size limit and also provide upper bounds for the estimation of fundamental gaps as indicated above. The upper bound property of cluster gaps enables us to probe usefulness of alternative extrapolations to bulk schemes, such as, for instance, 1/ √ N scaling, 49 see SM. Since both our estimates of ∆ f from the periodic supercells are greater than the cluster estimate, there is likely a small bias in the evaluation of excitations in the periodic supercells. In particular, the fixed-node errors in the excited and the ground states are most likely not identical and could be also intertwined with the remnants of the finite size errors even at larger sizes.
The convergence of the excited states is much more challenging in extended 2D systems than in bulk. The key reason is that the non-periodic direction enables changes in the single-particle densities that have dominant contribution to the periodic Coulomb interactions, a problem that is naturally avoided in cluster geometries; for details see SM. Therefore, we estimate that a systematic uncertainty of 0.1 to 0.2 eV could be present in our periodic-supercells gap calculations. The multiple calculations of the excited states that we present clearly illustrate this aspect.
The error is certainly greater for the excited-state energy, than for the ground state. Since the variational principle behind QMC gives the upper bound on the energy, it follows that the cluster value is a superior estimate for the fundamental gap. We believe that the true value of ∆ f for intrinsic phosphorene is about 2.4 eV, as also indicated in Tab. I. The GW values, defining the current state-of-the-art, are, compared to QMC, widely scattered and systematically underestimating the gap.
Exciton binding energy, optical gap, and comparison with experiment. Using the universal linear scaling, 14 ∆ b ≈ 0.27∆ f , we can also estimate the excitonic binding energy of phosphorene as ∆ b ≈ 0.65 eV. The optical gap is then ∆ o = ∆ f − ∆ b ≈ 1.75 eV. This is consistent with the optical absorption experiment 26 which reports 1.73 eV. As we discussed in the introduction, the optical gap of 2D semiconductors should be insensitive to the dielectric environment, justifying the consistency claim.
In terms of the fundamental gap obtained in experiment, (see Tab. I), photoluminescence emission spectroscopy (silicon substrate, no capping), 23 optical absorption (sapphire substrate, h-BN capping), 26 and scanning tunneling spectroscopy (no capping), 16 yield values of ∆ f of 2.2, 1.8, and 2.0 eV, respectively. The optical absorption value was attributed to the h-BN capping which also seems to reduce the exciton binding energy to just 0.1 eV. 26 The photoluminescence value of 2.2 eV is closest to our QMC prediction. Considering that the sample was on a dielectric substrate which lowers ∆ f by perhaps 0.1-0.2 eV, 26 our result is also consistent with this experiment.
Cohesion energy. One important quantity which, to the best of our knowledge, has not been determined experimentally yet, is the free-standing phosphorene cohesion energy. The cohesion energy of bulk black phosphorus is 3.43 eV/atom. 50 Compared to ordinary semiconducting 3D materials the phosphorus crystal is quite different being a van der Waals system of covalently bonded slabs. Even though the system shows moderately large cohesion, the stability of black phosphorus is rather low mainly due to an easy detachment of P 4 molecular units that per atom are bonded almost as strongly as the bulk material. Having calculated the phosphorene DMC ground state energies for a variety of sizes we can estimate the cohesion energy in the thermodynamic limit. Since the P 4 molecule is the dominant sublimation product from phosphorus solids, to estimate the cohesion we follow the thermodynamic path P 4 molecule → bulk black phosphorus → phosphorene using also the QMC calcu-lation of the van der Waals binding energy estimated recently. 28 The finite-size scaling, Fig. 3, yields a fixednode corrected cohesion energy of 3.268(4) in DMC/PBE and 3.284(6) eV/atom in the DMC/B3LYP periodic calculations, while from the clusters the estimated value is 3.26(1) eV/atom, see the SM for more details. Note that our estimated DMC values of cohesion energy of phosphorene from the two models are almost identical and very close to the estimation of 3.35 eV/atom for the 2D system obtained by using the van der Waals interlayer bonding from another DMC calculation. 28 The corresponding DFT values are 3.09 eV/atom 51 and 3.45 eV/atom. 52

IV. CONCLUSION AND PERSPECTIVES
We have performed systematic fixed-node QMC calculations of the quasiparticle band gap of free-standing single-layer phosphorene in both periodic and cluster settings. Using the universal linear scaling between the gap and the optical binding energy, we have also extracted the optical gap, which can be compared with experiments done on phosphorene on dielectric substrates or on encapsulated samples. Our results are consistent with available optical absorption and photoluminescence emission spectroscopy experiments. We argue that previous calculations based on GW underestimate the quasiparticle gap and do not give consistent predictions for phosphorene. Our ground state is essentially exact, evidenced from the calculated cohesion and its agreement with available (indirect) experimental data.
Our calculated quantities, band gap and cohesion, are key inputs into more qualitative and approximate theories, such as tight-binding and k·p, as well as into experimental interpretations. In particular, there is a clear path how our results can be modified to accommodate dielectric environments which enter the experiments, leaving the core of our QMC results unchanged. This can be done by Bethe-Salpeter modeling or by use of model dielectric screening 34 to find the appropriate value of the exciton binding energy. Hence, our explicitly many-body QMC results provide a reference ground for further studies on phosphorene based on strain and layer engineering as well as chemical doping and structural defects and indeed in any other 2D material. We have also demonstrated that these cutting-edge calculations are now feasible for a range of 2D systems and we expect the QMC methods to find a place at the top of the list of the toolkit for studying 2D systems. This is underscored by the fact that the scatter of predicted values for electronic parameters is significant not only in phosphorene, but in other 2D materials as well. 34 Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de). R.D., K.T., and I.S. acknowledge that the results of this research have been obtained using the DECI resources Beskow at the PDC Center for High Performance Computing at the KTH Royal Institute of Technology, Stockholm, Sweden from the PRACE project NaM2D and Sisu based in Finland at CSC IT Center for Science with support from the PRACE DoCSiNaP project. We also gratefully acknowledge use of the Hitachi SR16000/M1 supercomputer system at CCMS/IMR, Tohoku University, Japan. We also acknowledge useful discussions with Martin Gmitra and Alexey Chernikov.

I. SIMULATION DETAILS
The electronic fundamental band gaps ∆ f , were determined for both, extended and cluster systems with B3LYP 1 DFT orbitals and also extended systems with PBE 2,3 DFT orbitals. Lattice parameters were fixed to the experimental values of black phosphorus crystal 4,5 .
The fixed-node QMC 6,7 electronic gap calculations were performed in three steps: 1) the trial wave functions for ground-and excited-states were constructed from DFT B3LYP 1 and PBE 2,3 singleparticle DFT wave functions. Note that, unlike the GGA-type PBE, the hybrid B3LYP functional yields a fairly realistic band gap of 1.7 eV already at the DFT level (see Table I of the main text). Note also that no charge transfer is associated with the HOMO→LUMO transition, indicating that range-separated hybrids, such as CAM-B3LYP functional 8 , may not be needed. 2) the ground state trial wave function constructed from the DFT wave functions was optimized using VMC (variational Monte Carlo) techniques, and 3) finally fundamental gaps were computed from DMC (diffusion Monte Carlo) energies of ground-and excitedstates as first singlet-singlet vertical excitation energies, where always the Jastrow factor of the ground state was used: with E 0 , E 1 being, respectively, the ground-and the first excited-states obtained by fixed-node DMC 6 (PBE/DMC and B3LYP/DMC) not allowing any relaxation of the DFT nodal hypersurfaces due to the HOMO→LUMO electron excitation. Our QMC gap estimates use the vertical transition approximation, i.e. they neglect all the adiabatic, vibronic and zero-point vibrational energy effects. Such approximations tend to increase the gap value compared to the experiments, but on much smaller scale than that of importance here. * jaroslav.fabian@ur.de † ivan.stich@savba.sk For DFT modeling we use the CRYSTAL code 9 , while for all VMC and DMC calculations the QWalk code 10 was used. In all calculations the atomic cores were replaced by effective core potentials (ECP) 11 used in combination with a VTZ basis set.
Dynamical correlations were explicitly built in into the trial wave functions via isotropic Schmidt-Moskowitz Jastrow factors 7,12 , including electron-electron and electronnucleus correlations. Parameters of the trial wave functions 7 were optimized by minimizing linear combination of energy and variance 13 . The T-moves scheme 14 was used to keep the calculations with ECPs variational. Our systems are 2D+h in nature, i.e. while being 2D, they also have a finite thickness h. In order to take that structural characteristics into account, the Yeh-Berkowitz 15 modification of 3D Ewald summation technique for systems with a slab geometry that are periodic in two dimensions and have a finite length in the third dimension was adopted.
Two different types of finite-size scaling were performed: 1) a series of L × L supercell approximants where L goes from 1 to 6, see Figure 1 of the main text, and 2) 4 × 4 and 5 × 5 H-terminated cluster approximants, see Figure 1 of the main text. More details of the finite-size scaling are given in Sec. III.
The choice of these two types of models (periodic and finite) was driven by difficulties in estimating the band gap in the periodic setting. As explained and demonstrated below on solid data, the periodic calculations for excited states converge to the thermodynamic limit very slowly. Remarkably, we find that the cluster models converge to the infinite limit much faster and with much smaller variations of extrapolated estimators.

II. COHESION ENERGY
We start the estimate of the cohesion energy from the P 4 molecule that by bonding patterns and energetics is very close to white and red phosphorus bulk crystals. P 4 is also the dominant sublimation product from phosphorus solids. Therefore we will use the following thermodynamic path molecule P 4 ↔ black phosphorus bulk ↔ phosphorene to find the cohesion. For P 4 binding we initially use very accurate calculations and experimental data from the W4 testing set 16 . Coupled Cluster extrapolation to infinite basis and also experiment give an atomization energy of 285.03 kcal/mol (experiment being within 1 kcal/mol) and a zero point energy (ZPE) of 0.98 kcal/mol per atom. Our best Coupled Cluster extrapolated estimation is within 1.5 kcal/mol of the W4 value, see Table I. We also estimate the core-valence correlation and relativity corrections of −0.4 kcal/mol per atom 17 from calculations of P 2 . Therefore the P 4 atomization energy (bottom of the well, infinitely heavy nuclei) at T = 0 K is 286.95 kcal/mol or 3.110(9) eV/atom.
Cohesion of the bulk black phosphorus can be estimated from the sublimation energy of black phosphorus (P 4 ) → gas(P 4 ). We take a value of the heat of sublimation of 25.5(2.5) kJ/mol per atom or 0.26 eV/atom 18 (or almost the same value from an alternative experiment that estimates the cohesion energy difference between white and black bulks to be 21.2(2.5) kJ/mol and the heat of sublimation from white phosphorus to gas of 3.4 kcal/mol per P 4 19 so that we sum it to 25.5(2.5) kJ/mol = 0.26 eV/atom). Therefore the bulk black phosphorus cohesive energy with infinitely heavy nuclei is about 3.37(9) eV/atom (we assume that ZPE/atom is the same as in P 4 and this assumption also increases the error bar).
Bulk black phosphorus consists of stacked layers of phosphorene that are bounded by ≈ 0.08 eV/atom 5 . Therefore a reasonable estimate for phosphorene cohesion that we can infer from this data is E phosphorene coh,T ≈0,∞nuclei = 3.29(9) eV/atom. We neglected relaxation of phosphorene layer in vacuum as compared to the layers in black phosphorus bulk but that is probably within the approximately estimated error bar.
Let us now analyze the fixed node diffusion Monte Carlo (FNDMC) calculations and results.
Single FNDMC atom with PBE DFT orbitals, Burkatzki-Filippi-Dolg pseudopotentials 20 , and single reference is  with B3LYP nodal hypersurfaces. Both energies have fixed-node errors. We approximately estimate them from results of exact CCSD(T) in complete basis set limit (CBS) for the P atom and P 4 cluster since its cohesion/binding and bonding patterns are similar to phosphorene. The results of Burkatzki-Filippi-Dolg pseudopotential 20 calculations for a P atom and a P 4 cluster are in Supplementary Table I where the uncertainty is coming fully from the fixed-node correction since the statistical error bars are well below this. These fixed-node corrected cohesion energies compare also favorably with the estimate of 3.29(9) eV/atom given in the introduction.

III. FINITE-SIZE SCALING
A. HOMO-LUMO and ∆ f gap The fundamental gap can be determined from total energies as 6 : with N e being the number of electrons. In the GGA approximation, ∆ f converges to the DFT HOMO-LUMO gap in the N e → ∞ limit and demonstrates that in the infinite limit at the DFT GGA level there is no interaction between the electron and the hole in the exciton and the HOMO-LUMO gap corresponds to the fundamental gap, ∆ f . Such a conclusion is also valid in the generalized Kohn-Sham theory 21 , i.e. also for hybrid functionals used in the main text. This is a key property for understanding our QMC cluster results in the main text. However, since the LUMO level is constructed identically for finite and infinite systems, the DMC cluster results will, in the infinite system size limit, converge to the true DMC value of ∆ f . In addition, due to the confinement of the excitation, which can only increase the band gap, the cluster results will provide a strict upper limit for our extended system results.

B. Finite size scaling of QMC gaps
An analytical formula for the asymptotic finite-size scaling of the band gap is not known. In the main text we use the scaling versus 1/N . However, our calculated gap values could equally well be fitted with other fitting formulas, such as for instance with any 1/ n √ N formula. our data versus 1/ 2 √ N , Supplementary Figure 2, yields a fit that extrapolates to ∆ f of 3.08 eV. This, though, is in disagreement with the upper bound provided by the cluster calculations which exclude any higher root extrapolation than 1, i.e., the linear finite-size extrapolation versus 1/N is used the main text. Further analysis of the behavior for larger sizes supports these conclusions and provides therefore a rather solid estimation. This is further supported by the analysis of the cohesion energy that shows that the ground state is converged at sizes 4×4 or larger (see main text Fig. 3). The convergence of the excited state in periodic setting is slower due to the presence of vacuum in the direction orthogonal to the phosphorene slab. This enables uninhibited restructuring of the charge in this direction and thus generates longrange (artificial) interactions. In particular, the excited state, as it is clear from Figure 2 of the main text, shows basically an antibonding pattern along one of the P-P bonds that is almost orthogonal to the phosphorene plane (xy-direction) and therefore possesses a significant tail in the z−direction. The excited state is then a repetition of the same bonding pattern at all symmetry-related bonds since the excitation corresponds to the Γ-point symmetry. This results in arrays of dipoles on both sides of the slab and contributes significantly to the potential energy. Indeed, this is amply visible in the corresponding HF energies for the excited states (see Sec. C). The convergence to the thermodynamic limit is therefore very slow and extrapolation is quite difficult and very demanding on computational resources.

C. Hartree-Fock band gaps
The main contribution to the band gap, the Hartree-Fock part, is plotted in Supplementary Figure 3. This was calculated from the VMC wave function with B3LYP orbitals, setting the Jastrow correlation factor to 1. The Hartree-Fock band gap is extrapolated to 4.720(80) eV in the infinite size limit.

IV. EXTRAPOLATIONS BASED ON CALCULATED CORRELATION ENERGIES
Having shown in Sec. III D that the correlation energies tend to saturate at smaller system sizes than the Hartree-Fock energies, we now show that this fact may be used to estimate energies and band gaps for larger approximants. We applied this idea to the 6 × 6 approximant in Supplementary Figure 4. Finite-size scaling of the correlation energies of the ground-state S0 (green line) and the first excited-state S1 (purple line) versus 1/N . the cluster B3LYP/DMC treatment. From the smaller 4 × 4 and 5 × 5 DMC cluster calculations the (converged) correlation energies per P and H atom can be estimated, provided the Hartree-Fock energies are also known, Supplementary Table II. The relevant numbers are compiled in Supplementary Table II. From the energies compiled in the Supplementary Table II we estimate a fundamental gap of ∆ 6×6, clst, est f = 2.580(150) eV, which is almost within the error bar of the extrapolated value. It confirms that the band gaps of the clusters beyond 4 × 4 are much more weakly dependent on the size. This enables us to cross check the values from periodic calculations that show much stronger variation with the system size.
Consistency of cluster and periodic calculations can be ascertained also for the total energies per atom and the subsequent cohesion estimation. For clusters L × L we can use the following total energy model E L = N P E P + N H E H that adds energies of corresponding numbers of P and H atoms (see Table II). Based on this we can find linear relationships between different E L such as, for example, E 6 = (12/5)E 5 − (3/2)E 4 , etc. Influence of finite sizes is as we mentioned much smaller for clusters than for the periodic calculations. The data enables us therefore to estimate the corresponding total energy per P atom in the thermodynamic limit also in the cluster setting giving the following value It compares very favorably with the result from our periodic calculations (compare Sec. II, difference ≈ 0.02 eV/atom) and shows thus remarkably close agreement between these estimations despite significant difference in boundary conditions of the two models.

V. TABULATED DMC GAPS
All computed gap values are compiled in Supplementary Table III.