Fundamental gap of fluorographene by many-body GW and fixed-node diffusion Monte Carlo methods

Fluorographene^1,2 (FG) is a two-dimensional (2D) stoichiometric graphene derivative (C₁F₁) material exhibiting a large direct bandgap. This 2D material has attracted a significant amount of research activity over the past years due its promising physical and chemical properties.^3,4 Some of the key properties of interest include its large bandgap, dielectric characteristics, and potential surface physics and chemistry applications. Despite this attention, several crucial quantities that determine the electronic properties of FG has not been fully and unequivocally established. In the first place, this involves the fundamental and optical bandgaps as well as related derived properties.

Motivated by the existing discrepancies between theory and experiment^5–7 described below, we focus on fundamental electronic and optical FG properties. Among our key interests is the fundamental (energy) gap,

also interpreted as a quasiparticle gap (an energy to create a non-interacting electron–hole pair), where I is a first ionization potential and A is an electron affinity of a studied periodic solid, while E_N(±1) are total energies of N(±1)-electron ground states. Another quantity that is important in practical applications is an optical gap,

where $EN*$ is total energy of the first optically active excited state. The optically accessible Δ_opt state is, by definition, lower in energy than Δ_f that corresponds to the continuum limit of the exciton series, and the difference between the two is referred to as the exciton binding energy,

which may range from a few meV in bulk solid-state semiconductors such as Si to a few eV in materials exhibiting strong excitonic effects, such as gapped 2D nanomaterials.^8,9

At present, the onset of the FG optical absorption spectrum, Δ_opt, has been estimated by experiments to lie between 3 eV and 5 eV.^1,10–12 This range is simply too large to provide a reliable information for applications as well as for cross-validating existing theoretical estimations. On the theoretical side, Δ_f of FG predicted by the many-body perturbational GW approximation for the electron self-energy¹³ points to the interval of 7 eV–8 eV,^5,7,14 which indicates the presence of strong excitonic effects. Indeed, the Bethe–Salpeter equation (BSE) building upon the GW level of theory for the quasi-particle spectrum confirmed that E_b of FG is quite generally close to 2 eV,^5,7,14 irrespective of input sets of orbitals and technical parameters. However, the best available theoretical estimates of Δ_opt for FG [Eq. (3)] lead to a range ∼5 eV to 6 eV, i.e., they still embody a significant uncertainty that complicates a reliable cross-check with respect to the available experiments.^5,6

The main goal of the present study is to establish a reference quality estimation of Δ_f for a freestanding FG free of defects at a desired benchmark accuracy level of ∼0.1 eV. We employ two independent state-of-the-art and complementary many-body approaches,¹⁵ namely, GW perturbation theory pushed to its limits and real-space continuum fixed-node diffusion Monte Carlo (FNDMC) method that gained momentum in (correlated) materials research in past years.^16–21 In order to settle the exciton binding energy E_b as well, we build upon our high accuracy GW inputs and we solve the Bethe–Salpeter equation (BSE) so as to settle the value of Δ_opt at the mentioned accuracy level.

An important methodological part of this study is a discussion of an option to extract Δ_f estimates in solids (or periodic systems in general) from energy differences of neutral state total energies. The total energies are provided by FNDMC using single-reference/single-determinant trial wave functions, Ψ^T, with the symmetry of the periodic (super)cells Bloch orbitals. We present arguments that using periodic Bloch orbitals with built-in symmetries of ideal crystal structure under consideration very significantly restrict the possibilities for the occurrence of significant, localized excitonic effects. The basic argument is straightforward. The wave function built from Bloch orbitals is translation invariant. It is well-known that the corresponding nodal surface share the same property.^22,23 However, the localized exciton with significant binding would nominally break this translation symmetry. The nodal surface is simply too restrictive²⁴ for a significant charge restructuring that would enable appreciable electron–hole interactions, i.e., E_b = 0 in the thermodynamic limit. In such a situation, the energy difference of Eq. (2) that is formally equivalent to $ΔoptFN$ nominally corresponds to $ΔfFN$ [by Eq. (3)]. That is, in the thermodynamic limit, we arrive at an expression

that can be realized in calculations by constructing proper symmetry Ψ^T (see Sec. II C for details). Note that this is, in particular, pertinent to materials with large values of dielectric constants that also naturally suppress the formation of excitons. Here, we show that a similar conclusion can be reached also in the system where significant excitonic effects are present. Note that we emphasize the fact that the system is neutral. This simplifies the treatment of the finite size effects as well as the influence of charges that could decrease symmetries of the considered state and further complicate reliable estimations of the total energies that are consistent in the thermodynamical limit.²⁵ We also specify the limits when the equality above can become problematic, and we offer some insight into how this can be potentially treated. Experience^26–30 and rigorous arguments³¹ suggest that the FNDMC method provides unique information about bandgaps of periodic systems in a fully many-body framework, which is typically out of range of other correlated wave function methods. In most cases, the results are of benchmark quality so that we can write

and existing imperfections in this regard can be ascribed to current technical limitations, for example, sizes of supercells that truly reach the scaling regime with respect to the thermodynamic limit. The key FNDMC ingredients seem to be its following of proper constraints on E_N as a continuous function of N expected for an exact electronic structure theory,^32–34 linearity of E_N, and discontinuity of $∂EN∂N$ at integer N values.³¹ 

Equation (4) resembles a similar (but not quite the same) result derived for Kohn–Sham (KS) density functional theory (DFT) of solids,³⁵ 

where g^DFT is a one-electron KS valence band maximum (VBM)—conduction band minimum (CBM) gap. This equation holds for generic exchange-correlation (xc) functionals missing E_N derivative discontinuity when continuously crossing integer N.³⁵ The physics and accuracy of these two approximations to Δ_f [Eq. (4) vs Eq. (6)], however, fundamentally differ [see Sec. II C for a rigorous derivation of Eq. (4) valid for FNDMC; see Refs. 35 and 36 for a discussion of KS gaps in solids].

This paper is organized as follows: Sec. II introduces details on the used FG structure model (Sec. II A), GW and BSE setup (Sec. II B), FNDMC and promotion approach (Sec. II C), and quantum Monte Carlo (QMC) simulation setup used in production calculations (Sec. II D). The individual results for GW/BSE and FNDMC are presented in Secs. III and IV, respectively. Discussion of the agreement of methods vs each other in the context of experiment (Sec. IV) is followed by conclusions (Sec. V).

The unit cell of FG containing four atoms, C₂F₂, was assumed to adopt a chair-like conformation shown to be the most stable FG structure.³ Input geometrical parameters of FG were obtained by Perdew, Burke, and Ernzerhof (PBE) optimization from our previous studies^5,14 and were as follows: a lattice constant of a = 2.6074 Å, a C–C bond length of 1.5825 Å, a C–F bond length of 1.3833 Å, and a C–C–C bond angle of 110.9°. Reported DFT gaps have been only marginally dependent on minor FG structure variations³⁷ (as seen in structures optimized with various DFT functionals).^5,14 FG is a direct material with the electronic bandgap at the Γ point (for more details on the band structure, see Sec. III).

The Vienna ab initio simulation package (VASP)³⁸ implementing the projector augmented-wave (PAW) method^38,39 was used in all GW, BSE, and underlying DFT calculations. Since these calculations serve here as the cross-validating reference, it is worth noting that the VASP plane wave calculations based on PAWs can be considered as highly accurate, as shown in Ref. 40 using the all-electron full-potential linearized–augmented plane waves (FLAPW) method, which is regarded as an established reference method in condensed matter physics.

We have used the GW approximation (where G is the Green’s function and W is the screened Coulomb potential) for the electron self-energy Σ with input orbitals ϕ_nk(r) from DFT and Perdew, Burke, and Ernzerhof (PBE) xc functional⁴¹ electronic structure calculations. Therefore, the quasi-particle energies $ϵnkGW$ are calculated as first-order corrections to the DFT single-particle energies ϵ_nk,

and the quasiparticle gap as $ΔfGW=ϵCBMGW−ϵVBMGW$ [analogous to Eq. (6)], where Z is the normalization factor,^42,T is the kinetic energy operator, V_n−e is the potential of the nuclei, and V_H is the Hartree potential. Recently, it was shown on layered BN,⁴³ monolayer MoS₂,⁴⁴ and bulk ZnO⁴⁵ that well converged GW with starting orbitals from PBE agrees very well with experiment and that there is no necessity to use higher level of GW theory by iterating G or W. In order to reach the desired accuracy of our calculations, we used GW set of PAWs, cut-off energy E_cut = 600 eV, and strict electronic step convergence criterion (1 × 10⁻⁷ eV). Convergence of the critical parameters (plane wave cutoff, number of bands considered in GW, size of the vacuum region in the simulation cell, and k-point grid) is elaborated in more detail in Sec. III.

For insulating FG, the materials with completely unoccupied conduction bands and occupied valence bands (c and $v$ indices, respectively), the Bethe–Salpeter equation can be rewritten as the eigenvalue problem^46–48 

where ν is the Coulomb kernel 1/|r − r′|, the eigenvectors $Acvkλ$ correspond to the amplitudes of free electron–hole pair configurations composed of electron states |ϕ_ck⟩ and hole states |ϕ_$v$k⟩, and the eigenenergies $Eexcλ$ correspond to the excitation energies (with $ΔoptBSE≡Eexcλ=1$⁠).

FNDMC is a stochastic many-body real-space projector ab initio method^49–52 that can be effectively applied to systems with tens to hundreds of valence electrons. Therefore, it has been gaining traction in electronic structure calculations of real solids due to its accuracy, scalability, and versatility.^{26,27,29,53–58} For a Hamiltonian $Ĥ$⁠, it projects out the most optimal ground state within the fixed-node constraint, which has non-zero overlap with the antisymmetric trial state Ψ^T. Using imaginary time τ, we can write

where E_T is an offset energy that keeps the norm of Ψ^FN asymptotically constant. The superscript FN indicates that Ψ^FN and Ψ^T share the same fixed node F = {R; Ψ(R) = 0}, where R is an electron position vector. FNDMC simulation provides a stochastic estimation of the total energy expectation value, $EFN±σK$⁠, where σ is a local energy [$EL(R)=ĤΨT(R)ΨT(R)$] dispersion and K is the number of independent samples. E^FN is an upper bound to the exact ground state⁵⁹ energy of $Ĥ$ within the same symmetry class as Ψ^T, thus enabling the study of excited states by picking Ψ^T orthogonal to the lower states.^60–63 The related FN bias grows quadratically with the increasing inaccuracy of F.⁵² Viable sophistication in constructions of trial states is used to marginalize the FN bias as much as possible.^64,65 Practical FNDMC computations that are feasible for large systems nevertheless require computationally efficient Ansätze, such as the widely used Slater–Jastrow trial wave functions,⁶⁶ Ψ^T = Ψ^SJ. Here, Ψ_S is a single-reference Slater determinant part and J is an explicitly correlated Jastrow term⁶⁷ that is appropriately parametrized. The trial wave function is variationally refined, and then, it serves as a starting point for the FNDMC calculation.

In periodic systems, it is a usual practice to pick Ψ^S of Bloch orbitals^20,25,26,68 for a ground state and then use the same orbital set for a construction of an electron-promoted state,

and charged states,

where a and a^† are usual creation and annihilation operators in second quantization of electrons in VBM and CBM, respectively. Note that in the thermodynamic limit, orbital relaxation is irrelevant.⁶⁹ Clearly, a one-electron perturbation (electron addition or removal) of an infinite crystal is negligible so that the ith orbital in N and N ± 1-electron system is indistinguishable, ϕ_i(N ± 1) ≈ ϕ_i(N). We continue with a textbook-level derivation of symbolic expressions leading to Δ_f and Δ_opt so that the explicit difference is apparent [Eq. (3)] in order to stress that in special cases, theoretical estimates of these quantities may be nominally equivalent.

Consider a system of N electrons where two electrons in two unrelaxed orbitals, ϕ₁(r) and ϕ₂(r) (single doubly occupied and single unoccupied orbital, e.g., VBM and CBM), will be considered explicitly. The remaining unspecified part (other N − 2 electrons, external potentials, etc.) is denoted as “rest” (R). The total ground state energy of a neutral state in one-determinant approximation with two electrons in ϕ₁ amounts to⁷⁰ 

where ϵ_i is an energy of an electron in orbital i, E_R is a total energy of a subsystem R, and J_ij (K_ij) is a Coulomb (exchange) interaction of electrons in orbitals ϕ_i(r) and ϕ_j(r), respectively, or an interaction of an electron in orbital ϕ_i(r) with R (tilde indicates an approximation). The total energy of a system where one of the electrons has been “promoted” from ϕ₁ to ϕ₂ amounts to

An optical gap approximation follows immediately [Eq. (2)]

The total energy of an anion state with N + 1 electrons (one electron added to ϕ₂) may be expressed as

whereas the total energy of a cation state with N − 1 electrons (single electron removed from ϕ₁) amounts to

The corresponding ionization potential and electron affinity $Ĩ$ and $Ã$ read

respectively, enabling us to express

Finally, we arrive at the desired explicit approximation to the exciton binding energy [using Eqs. (3), (14), and (19)],

that is independent of R. It depends on the initially selected orbitals occupied by the explicitly considered electrons. J₁₂ is a Coulomb repulsion of the electrons in orbitals ϕ₁ and ϕ₂, respectively. Since ϕ₁ is also a hole orbital left behind by the promoted electron, energy −J₁₂ can be interpreted as an electron–hole attraction (in electron–hole picture, $Ẽb=K12−J12$⁠). K₁₂ is an exchange energy that can be alternatively computed from the total energy difference involving triplet (spin–flip) promotion

leading to

Note that the expression of Eq. (20) is valid in a system of arbitrary size, including the thermodynamic limit since in a unit cell containing N electrons, N − 2 of them are swept into R. In the same way, for supercells containing κ unit cells (used to approach the thermodynamic limit by taking κ → ∞), R contains κN − 2 electrons.

In addition, Eq. (20), based on one-determinant total energies, is a very accurate approximation of the exact E_b in systems with a large number of electrons N approaching the size-extensive limit, if the corresponding exact wave functions are all dominated by the one-determinant and the selected orbitals are eigenstates of the one-particle reduced density matrix of the exact solution. In such a case, one can observe that analogous derivation using

where C_N, etc., denote the correlation energies missing at the corresponding mean-field level, will result in

where $C=CN−1+CN+1−CN−CN*$⁠. In the size-extensive limit, it follows that

and, together with an assumption $CN≈CN*$ that is reasonable for large N, that implies

Therefore, the correlation energy contributions largely cancel out. We can therefore state that to a very good approximation in systems with a large N,

Equation (20) is therefore very relevant for a vast class of orbital-based methods, including many-body approaches such as FNDMC. This result reaches beyond just formal importance. One can check that the FNDMC correlation energy converges rapidly with respect to the supercell size, reaching a constant already for sizes of practical relevance such as 2D ∼4 × 4, as has been observed in phosphorene.²⁰ 

Before focusing on FNDMC, we are interested in general $Ẽb$ of an extended system. Let us identify three possible scenarios with respect to the nature of the chosen orbitals ϕ₁ and ϕ₂:

• Both orbitals are of Bloch type (e.g., VBM and CBM).

• One of the orbitals is of Bloch type, whereas the other one is localized (translation symmetry is broken).

• Both orbitals are localized.

For the scenarios (i) and (ii), it is possible to demonstrate that

and

so that $Ẽb=0$⁠, whereas for the scenario (iii), one gets finite non-zero results in the κ → ∞ limit, allowing also nonzero $Ẽb$⁠. An important and perhaps obvious conclusion is that the result is orbital dependent. If any of the orbitals is a Bloch state obeying the crystal translation symmetry, in the thermodynamic limit, an electron and hole effectively do not interact, $Ẽb=0$⁠, and $Δ̃f=Δ̃opt$⁠.

As an example, let us sketch the critical steps of the proof for Eq. (28) in 2D and scenario (i). Let $ϕi(r→1)$ and $ϕj(r→2)$ be one-particle Bloch orbitals, expanded in a plane wave basis set,

normalized on region Σ = (κL) × (κL), where κ is an integer, $k→=πκL[n1,n2]$⁠, and $w→=πκL[m1,m2]$⁠, where n₁, n₂, m₁, and m₂ are integers, respectively, and L is a constant unit cell size. The Coulomb integral of our interest may be expressed as

Substituting Eqs. (30) and (31) into (32) leads to an expression that contains the integrals of type

where $λ→=k→+w→,μ→=k→′+w→′$ and $Wmlλ→μ→$ and $Fmlλ→μ→$ are expansion coefficients. The solution leads to

where $Ak→k→′w→w→′$ is an expansion coefficient independent of κ (see the supplementary material for more details). The limit of our interest [Eq. (28)] clearly vanishes, since the logarithmic singularity from the Coulomb interaction operator in Q is suppressed by the normalization factor. Similar conclusions hold for the case of K₁₂ [Eq. (29)]. In addition, both conclusions were also verified analytically for 1D and scenario (i) (see the supplementary material for details) and numerically (using dense real-space grids) for scenarios (i) and (ii) with a few specific cases of sine orbitals in 1D and 2D.

We can finally discuss consequences of the presented conclusions for FNDMC simulations using the Ψ^T determinant constructed of Bloch one-particle orbitals in combination with Eq. (2). What do we obtain from a single-electron promotion approach within the commonly used single-reference (single-determinant) FNDMC?

As mentioned above, FNDMC projects out an exact ground state consistent with the nodal boundary condition fixed by F(Ψ^T), which enforces antisymmetry and periodicity of the simulated state Ψ^FN [Eq. (8)]. In this way, F imposes a hard bound for the lowering of FNDMC energy in such calculation. The orbital promotion leads to another F with a different shape if the state is of different symmetry, but still translation invariant, and therefore, it locks in corresponding bound for the excitation. In case the excited state symmetry is the same as that of the ground state and there is no degeneracy, the same holds to a significant degree since any overlap is usually negligible, and the systematic errors are dominated by the fixed-node bias. It is also known that if there is a degeneracy in excited orbitals, this might not hold, but this is a different issue, which has little to do with excitonic effects considered here. Clearly, in a system where an exact exciton many-body state $ΨNexciton$ is an eigenstate of $Ĥ$⁠, for $F(ΨNT*)$ based on Eq. (9) where VBM and CBM are Bloch states obeying the lattice translation symmetry, we have $ΨNT*|ΨNexciton≈0$⁠. Therefore, FNDMC projects out the total energy solution that corresponds to the promotion in the considered band structure. Of course, this would change if the promotion would be built from orbitals that would enable translation symmetry breaking so as to approximate the electron–hole pair (for a broken-symmetry state, $ΨNT*,bs|ΨNexciton≠0$⁠). There is an additional issue of finite sizes that forces the electron and the hole to be “closer” than in the true thermodynamic limit. However, this does not change the above argument based on the fixed boundary stemming from the nodal surface. It would affect only the Mott–Wannier limit where the exciton is so delocalized that within the given supercell, it appears as translation invariant; however, the corresponding effects are tiny and usually fall below the error bars.⁷¹ 

Let us partially summarize possible outcomes of the single-reference single-electron promotion FNDMC calculations. For molecules, energy differences based on the described promotion approach are expected to approximate optical transition energies since all involved orbitals are localized in space and the translation symmetry is absent.¹⁵ Any inaccuracies originate from mundane issues such as virtual orbital imperfections, energetic ordering, and amount of multi-reference mixing. At least partial charge restructuring is common since the system has no boundary, and the charge can shift from inside to outside and vice versa as the energetics dictates, to the extent that the nodal surface allows it.⁷² This is straightforward to check in anions⁷³ where it is common that the tails of the trial wave functions might not be correct but can be repaired by the FNDMC solution. As we argued already before, in the low density limit, the fixed-node errors become less significant and provide exceedingly accurate energy differences.⁵⁷ 

In solids, nodal surfaces of trial wave functions with Bloch orbitals impose new restriction²⁴ so that the promotion into another Bloch orbital will nominally result in $ΔfFN$⁠. At the same time, we have to clearly state that this is not universal. Let us consider cases where at least part of the orbitals (occupied and virtuals) remain localized on constituents even in the thermodynamic limit (molecular crystal) or where periodicity is broken to a significant extent (defects, impurities, etc.). In such systems, the promotion FNDMC method result will be somewhere between the corresponding optical and fundamental gaps. For example, in van der Waals molecular crystals, the periodicity affects the ground states only marginally, and even the lowest excitations might not be changed too significantly. Therefore, this would apply provided that (i) there is “enough space” for local charge restructuring similar to an isolated molecule, and (ii) that there are low-lying excited orbitals that mimic the isolated molecular states very closely even within the Bloch periodicity. In such cases, one can expect that the resulting energy difference would be a better approximation to the optical transition(s), i.e., $ΔoptFN$ will result. It is therefore important to probe for the character of corresponding states in order to find out the relevant limit.

In the case of “in between” situations, one has to work on the trial wave function to build better approximation to one or both limits, and the interpretation of the results has to be verified accordingly. Note that one can obtain some indications, although for the time being a bit crude, when the band limit applies. For example, in another project (to be published elsewhere), the band energies/gaps do not change regardless of how the promotion is built in the spin symmetry channels. In particular, in the Si solid, one can check that the promotion energies remain essentially constant for the following three types of trial wave functions:⁷⁴ 

• promotion of an electron in a given spin channel (i.e., really a mix of true singlet and triplet states),

• promotion of an electron with simultaneous spin flip (true triplet state), and

• two-determinant single-electron excitation in the proper singlet state (i.e., open-shell singlet).

Some of such tests have been carried out also in the present calculations, as explained in Sec. IV. These results indicate that the exchange and correlation effects are indeed delocalized, and the Bloch periodicity provides essentially a full lock on the corresponding energy difference. Note that the electron correlation still matters since the obtained bandgap differs from the mean-field. Although this is most probably far from being universal and it requires further study on more systems, it provides an illustration that the presented analysis above is essentially correct. It is also intuitively clear that in a contrasting case of an isolated molecule, reasonably localized singlet and triplet excitations would differ rather significantly.

It is highly desirable to reach beyond the current status and to build accurate trial wave functions for both limits. At present, however, we still encounter significant difficulties in achieving such a goal. The complications are mostly on the technical side due to restricted availability of efficient tools for multi-reference trial functions in periodic setting. This can change in the near future so that this is clearly a promising subject for further research.

The QMC total energy computations consisted of the commonly used stages:^20,75

• DFT one-particle Bloch orbitals were obtained using Gaussian G09⁷⁶ or Quantum Espresso (QE),⁷⁷ DFT codes with PBE or PBE0⁷⁸ xc functionals, and aug-VTZ basis⁷⁹ or plane wave basis set (a kinetic energy plane wave cutoff of 100 Ry), respectively. Brillouin zone sampling meshes ranged from 2 × 2 × 1 to 5 × 5 × 1, consistent with the target real-space supercell sizes used within real-space Γ-point QMC computations (direct Γ − Γ transition). All computations used norm-conserving Burkatzki–Filippi–Dolg (BFD) effective core potentials (ECPs).⁷⁹ 

• Variational Monte Carlo (VMC) optimizations of the parametric Jastrow explicit correlation terms entering Ψ^T used two-particle (electron–electron and electron–nucleus) terms and 2 × 2 supercell (QWalk code⁸⁰) or up to three-particle terms (electron–electron, electron–nucleus, and electron–electron–nucleus) within the 3 × 3 supercell (QMCPACK⁸¹).

• FNDMC total energies, corresponding to Ψ^T of ground-states and single-electron VBM-CBM-promoted states, were obtained from supercells containing κ = 4, 9, 16 (QWalk) or κ = 9, 16, 25 (QMCPACK) unit cells with VMC-optimized ground-state Jastrows, a time step of 0.01 a.u., and T-moves approximation for the treatment of ECPs.⁸² 

• Finite size effects were suppressed by linear extrapolations of the series of production results to the thermodynamic limit by taking 1/κ → 0, where κ is the total number of unit cells present in the simulation supercell. We have verified that misfit of the ground state and promoted state FNDMC total energies normalized to κ amounts to 0.002 a.u. in thermodynamic limit so that the total energy consistency condition was not enforced.²⁵ 

The computations of the two sets of the FNDMC results for a series of FG supercells using different codes and setups, QWalk (Gaussian basis sets, Coulomb interaction, two-particle Jastrows, and a vacuum size of 40 Å) and QMCPACK [plane waves, model periodic Coulomb (MPC) interaction,⁸³ three-particle Jastrows, and a vacuum size of 20 Å], show consistency over the overlapping system sizes (3 × 3 and 4 × 4). Therefore, we rule out some of the possible method biases: one-particle basis set bias, slab-model vacuum size, and modified electron interaction.

Additional tests confirming robustness of the chosen production QMC setup for total energy differences were performed as well. Time step error tested by considering 0.01 a.u. vs 0.005 a.u. resulted in a difference of the order of error bars ∼0.1 ± 0.1 eV in energy differences. Shape effects of one-particle orbitals were ruled out by considering the hybrid xc functional (PBE vs PBE0 in a 4 × 4 supercell within QMCPACK), and the results were comparable within the error bar (6.11 ± 0.05 eV vs 6.03 ± 0.06 eV). Similarly, CBM orbital replacement by the equivalent orbital from DFT computation using triplet spin multiplicity resulted in no detectable bias within FNDMC energy differences. BFD ECPs transferability was checked vs a new generation of correlation-consistent ccECPs:⁸⁴ in a 4 × 4 supercell, the promotion energy differences were indistinguishable within the statistical error (6.11 ± 0.05 eV vs 6.12 ± 0.1 eV). The Ewald summation consistency test for our 2D FG slab model was tested by considering slab-modified Ewald interaction (Yeh–Berkowitz⁸⁵) within the appropriately adapted QWalk code; no deviation from to the usual 3D Ewald formulation was observed for promotion energy differences within the statistical error of ∼0.06 eV.

In this section, we present refined and extrapolated GW quasiparticle gap values, significantly improving upon the previous calculations.^5,14 The currently available best GW values rely on an old generation PAW potentials and do not consider the converged cutoff and/or vacuum size and/or sufficient Brillouin zone integrations and spread between 6.98 eV and 8.28 eV,⁵ depending on input sets of orbitals and level of GW calculations (e.g., 7.01 eV,⁷ 6.99 eV, and 7.82 eV¹⁴).

Here, using the new generation of GW-optimized PAW potentials, we obtained ultimate convergence of the simulation parameters, namely, an inter-sheet distance L_z → ∞ and k-point grid (N_k → ∞). The number of bands N_B and energy cutoff in GW calculation $EcutGW$ were investigated in detail (Fig. 1) and determined the production setup: $EcutGW$ = 400 eV and N_B = 1152. Subsequently, the effect of inter-sheet distance L_z was investigated for various k-point grids. Because of the nonlocal nature of the GW approximation, the quasiparticle gaps depend on this parameter and keep on increasing as 1/L_z. The slope of the corresponding linear fit significantly differs for various k-point grids (Fig. 2); extrapolations to zero 1/L_z yield quasiparticle gaps of 8.08 eV, 7.37 eV, 7.24 eV, and 7.19 eV for 6 × 6 × 1, 12 × 12 × 1, 18 × 18 × 1, and 24 × 24 × 1 k-point grid, respectively. Finally, the quasiparticle bandgap for a full k-point grid (N_k → ∞) was obtained using extrapolation of quasiparticle gap values obtained in Fig. 2 (1/L_z → 0; Fig. 3). The final production value of the quasiparticle gap obtained by the GW method is $ΔfGW$ = 7.14 eV. The agreement of the above-mentioned unconverged value of 7.01 eV⁷ vs our converged result may be attributed to the fortuitous bias cancellation.

For completeness, we calculated the optical gap (position of the first excitonic peak) by the BSE method. Taking all extrapolations of the relevant parameters into account, we obtain $ΔoptBSE$ = 5.21 eV. The binding energy of the first exciton in turn amounts to $EbBSE=$ 1.92 eV.

Finally, we provide an exciton wave function analysis showing contributions of the respective orbitals around VBM maximum and CBM minimum in k-space. The wave function of the λ-th exciton is expressed in an electron–hole product basis as $∑cvkAcvkλϕckϕvk$⁠. The eigenstate corresponding to the first bright (doubly degenerate) exciton in FG is visualized in Fig. 4 by plotting circles with $|Acvk1|$ radius into the band structures.⁸⁶ Pairs of large circles visible in Fig. 4 represent electron–hole pairs, which contribute the most to the wave function of the first excitonic peak $Eexc1$⁠. It turns out that region close to the Γ point (for the highest occupied band and lowest unoccupied band) dominates, whereas other regions of the Brillouin zone provide only negligible contributions to the studied first excitonic state. This salient feature of the FG electronic structure indicates that the nodal surface of Ψ^T of its first excitonic state will likely be well approximated by the multideterminant expansions within restricted active-space consisting of one-particle states in the vicinity of doubly degenerate VBM and CBM in k-space (two electrons in three orbitals).

The FNDMC production results from both considered codes, QWalk and QMCPACK, are reported in Tables I and II, respectively. Note that the presented FNDMC total energy differences (⁠$ΔfFN$⁠) approximate fundamental gap Δ_f although they were obtained by promotion approach based on two total energies [Eq. (4), elaborated in Sec. II C] instead of the conventional approach based on three total energies [Eq. (1)]. The results agree very well for the overlapping system sizes (3 × 3 and 4 × 4), ruling out biases possibly coming from electron interaction modification (Coulomb vs MPC interaction), Jastrow (two-particle vs three-particle), or slab vacuum size (20 Å vs 40 Å).

Linear extrapolations of the QWalk data for sizes ranging from 2 × 2–4 × 4 (Table I) to the thermodynamic limit lead to the gaps 6.5 ± 0.1 eV (three-point) and 6.9 ± 0.1 eV (using two latter points), respectively, whereas the latter value is more reliable, since, the first point corresponding to 2 × 2 supercell appears to be overly influenced by the pronounced size effects.

Linear extrapolations of the FNDMC data obtained from QMCPACK (Table II) to the 1/κ → 0 limit lead to the gaps 6.9 ± 0.1 eV (two-point and sizes 3 × 3 and 4 × 4, which well agrees vs the corresponding QWalk value) and 7.1 ± 0.1 eV (2-point and sizes 4 × 4 and 5 × 5), respectively. Linear fit for the full set leads to $ΔfFN=7.0±0.1$ eV, which is our best linear-fit FNDMC estimate of the desired quantity for FG. The quadratic three-point fit of the same results leads to the value of 7.1 ± 0.1 eV that is, within the statistical uncertainty, same as for the QWalk data.

Finally, we have evaluated promotion energy from the ground state to the triplet state (Table II; T₁) in order to verify the magnitude of K₁₂ [Eq. (22)] contribution to E_b (expected to vanish in the thermodynamic limit). It turns out that already for a finite supercell of size 4 × 4, the exchange contribution is negligible (0.03 ± 0.1 eV, well below our target accuracy standard).

The results from the two presented GW and FNDMC methods provide independent, state-of-the-art result in a strikingly similar fundamental gap estimates,

$ΔfGW=7.14$ eV,and

$ΔfFN=7.1±0.1$ eV.

We conclude that such an excellent agreement of the two complementary benchmark methods leaves only very little room for further significant corrections at this accuracy level, and we believe that our results settle an ultimate reference value for an ideal freestanding FG sheet

$Δftheor≈7.1±0.1$ eV.

Additional BSE computations enable us to ascertain the best available theoretical estimate for an optical gap to

$Δopttheor≈5.2±0.1$ eV.

It is interesting to observe that FG belongs to the class of materials exhibiting universal scaling of exciton binding energy and fundamental gap,⁸⁷ which predicts E_b ≈ Δ_f/4 (1.78 eV vs our value of 1.92 eV). On the other hand, this value is at the upper bound of experiments that report a lower optical gap value ∼3 eV to 5 eV. We must conclude that additional significant effects, beyond the scope of our study of an ideal FG at zero temperature, are likely to be responsible for the remaining discrepancy. The most probable candidates that would cause midgap states and lower optical absorption energies include structural imperfections such as vacancies^5,6 and contaminants.³ In addition, multilayer effects¹⁰ and termination group effects^1,12,88 cannot be ruled out either. Zero point energy and finite-temperature effects may also cause significant bandgap renormalization, as demonstrated for covalent systems such as diamond⁸⁹ and hexagonal BN,⁹⁰ and thus leave room for further FG studies.

The reference value of the fundamental gap of fluorographene was estimated by the state-of-the-art quantum ab initio methods: GW and FNDMC. Both approaches have been pushed to their limits independently and have arrived at the benchmark value of $Δftheor≈7.1±0.1$ eV for a freestanding material sheet free of defects at zero temperature, making little room for further corrections at the 0.1 eV accuracy level. The Bethe–Salpeter equation determined the first exciton stabilization of $EbBSE=1.92$ eV so that our best optical gap estimate amounts to $Δopttheor≈5.2±0.1$ eV. In addition, we pointed out a possibility to extract approximations to Δ_f from neutral single-reference Bloch-orbital FNDMC computations of promotion gaps in finite small-to-medium supercell sizes (vs the actual exciton size) extrapolated to the thermodynamic limit. Such an approach is appreciable since only two total energy FNDMC computations of neutral systems are required to obtain reasonable estimates of the fundamental gap. The alternative involves the usual three total energies that involve charged states and that can typically exhibit further complications such as more profound finite size effects.

See the supplementary material for details on the derivation of Coulomb and exchange integrals in 1D and 2D with Bloch orbitals in the thermodynamic limit [Eqs. (28) and (29)].

The authors are grateful to Rene Derian, Matej Ditte, and Jan Brndiar for fruitful discussions. Financial support by the Czech Science Foundation (Grant Nos. 18-25128S and 18-24321Y), the Institution Development Program of the University of Ostrava (Grant No. IRP201826), the Slovak Research and Development Agency (Grant No. APVV-18-0161), and the European Regional Development Fund (Grant No. ITMS2014+:313011W085) is gratefully acknowledged. The computations were performed at the IT4Innovations National Supercomputing Center (Grant No. LM2018140). L.M. acknowledges support by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, Theoretical Condensed Matter Physics, under the Award No. de-sc0012314.

The data that support the findings of this study are available within the article and its supplementary material.