Neutral band gap of carbon by quantum Monte Carlo methods

We present a method of calculating the energy gap of a charge-neutral excitation using only ground-state calculations. We report Quantum Monte Carlo calculations of $\Gamma\rightarrow\Gamma$ and $\Gamma\rightarrow X$ particle-hole excitation energies in diamond carbon. We analyze the finite-size effect and find the same $1/L$ decay rate as that in a charged excitation, where $L$ is the linear extension of the supercell. This slow decay is attributed to the delocalized nature of the excitation in supercells too small to accommodate excitonic binding effects. At larger system sizes, the apparent $1/L$ decay crosses over to a $1/L^3$ behavior. Estimation of the scale of exciton binding can be used to correct finite-size effects of neutral gaps.


Introduction
Quantitative predictions of electronic excitation spectra in solids has been a central topic from the very beginning of condensed matter theory. A first insight is usually provided by the band structure obtained by density functional theory (DFT). Many-body perturbation theory has successfully extended DFT or self-consistent field approaches to provide electronic excitation spectra including perturbative electronic interaction effects. However, the precise values of spectral quantities like the band gap are still sensitive to the underlying functional and approximation scheme [1]. Their accuracy is frequently established by comparison with experiment involving further modelling. Based on the variational principle, quantum Monte Carlo (QMC) methods [2][3][4][5] provide an alternative way to obtain ground and excited state properties whose quality may be judged independently of experimental measurements.
Particularly adapted to study electronic correlation effects in liquids and solids, QMC calculations for extended systems have become more and more affordable. Similar to many classical simulation methods, QMC calculations are done in a finite (periodic) simulation cell of linear extension , containing electrons. Within projector Monte Carlo methods, e.g., fixed-node diffusion Monte Carlo (DMC) [6], the systematic bias is entirely determined by the nodes of the underlying many-body wave function. Backflow [7], −body and iterated backflow [8,9] and related neural network wave functions [10,11] have offered a systematic way of improving the accuracy judged by the variational principle. Results in the thermodynamic limit are then obtained by numerical extrapolation of different system sizes, possibly accelerated via analytical correction terms. Dominant sources of finite-size errrors in electronic systems are fermionic shell [12] and Coulomb charge effects [13]. Since the studies of a fixed-node error are in general limited to small system sizes, the control of both of these systematic biases, fixed-node and finite-size error, will eventually determine the overall accuracy of QMC calculations.
Recently, we have shown that QMC calculations of the fundamental (charged) gap of an electron system, frequently called quasi-particle (QP) gap, defined by the electron addition and removal energies, converge slowly, inversely proportional to the linear extension ∼ −1/3 of the simulation cell [14]. Heuristically, the 1/ dependance results from the interaction of the additional charges across the periodic boundaries [15,16] which should be absent in neutral systems. Then, finite-size effects would be drastically reduced in calculations of the neutral gap, obtained by exciting the system at a fixed number of electrons. Therefore, QMC calculations of charge neutral, particle hole excitations [17][18][19][20][21][22][23][24][25][26][27][28] commonly assume a faster 1/ 3 behavior to extrapolate finitesize errors. However, QMC results of charged and neutral gaps of diamond Si [18,29] yields comparable values at fixed finite system size, not supporting the expected qualitative change of finite-size effects.
In this paper, we study finite-size effects of QMC calculations of the neutral gap. Similar to [14], we explain that leading-order finite-size effects are encoded in the asymptotic behavior of the static structure factor, (k), which is encoded in the many-body wave function underlying the calculation. The size effects are not caused by the charged nature of the excitation, but by their extended or localized character. In particular, promoting a single Bloch orbital from the ground state to an excited state in the determinantal part of the correlated wave function, leads to an excitation of extended character, thus it suffers from the same slow 1/ decay as their charged counterparts.
In section 2 we briefly introduce our theoretical methods which allows us to obtain the neutral gap as a ground-state calculation, discuss a possible use of twisted boundary conditions, and present the general form of the trial wave functions. In section 3 we present our QMC results for diamond carbon compared with quasiparticle-gap from [14]. In the next section 4 we discuss finite-size effects and we finally recollect our conclusions and perspectives in the final section 5.

Neutral gap
The neutral gap as defined above, equation (1.2), explicitly involves the energy of the first excited state. In practice, QMC methods of excited state are either based on exact symmetry constraints [30,31], or are approximately implemented in a minimization procedure [19,32]. However, any contamination of the excited state by the ground state will violate the variational principle, and may lead to a bias for the gap. Here, we show how to restate the neutral gap within a pure ground-state formalism.
For ideal solid structures, the Hamiltonian is invariant under the simultaneous translation of all electrons by any crystal lattice vector, t. Therefore, electronic energies, k , and eigenstates, Ψ k , can be labelled by a wave vector k in the first Brillouin zone in the reciprocal crystal lattice. Wave functions of a given crystal momentum k should satisfy Ψ k (r 1 + t, r 2 + t, . . . , r + t) = e ik·t Ψ k (r 1 , r 2 , . . . , r ). (2.1) In general, we can assume that the ground state is in the trivial representation, k = 0, with energy 0 ( ). In an insulator, all excitations of the ground state are gapped, k ( ) − 0 ( ) ⩾ Δ for k ≠ 0. Approaching the thermodynamic limit, k becomes a continuous parameter and we can write (2.2)

33701-2
Therefore, similar to the fundamental gap [14], also the neutral gap can in principle be obtained by ground state calculations imposing a non-vanishing crystal momentum. By continuity, assuming a normal band insulator, lim k→0 min k ( ) − 0 ( ) continuously connects to the optical or vertical gap of the bulk in the zero momentum sector (k = 0), so that optical excitation energies at vanishing crystal momentum can be obtained by extrapolation.
Although the conservation of crystal momentum followed from the assumption of a perfect crystalline lattice, this requirement can be relaxed to include also ionic zero point motion and temperature effects [33,34].

Twisted boundary conditions
Twist averaged boundary conditions [12] can be generally used to reduce finite-size errors due to fermionic shell effects. For a given twist , the many-body wave function is then chosen to obey whenever any electron moves across the simulation cell, e.g., in the -direction. Electronic eigenstates in the solid can then be labelled by and k.
The ground state energy per volume of a solid in the thermodynamic limit can then be approximated by averaging over twisted boundary conditions assuming a uniform grid containing twist angles. Here, k ( , ) denotes the ground state energy of an electron system at a given twist and crystal momentum k, is the volume of the simulation cell; and k denote the number of electrons and crystal momentum which minimize the ground state energy at fixed .
For a normal band insulator with non-vanishing fundamental gap, the number of electrons in the ground state at fixed does not depend on ; canonical and grand canonical twist-averaging give the exact same ground state energy in this case [14]. Although the crystal momentum, k , of the ground state wave function at a given does not vanish in general, the total crystal momentum of the ground state in the thermodynamic limit will be zero, k = 0, due to time-reversal symmetry. A non-vanishing total crystal momentum, q can be imposed by considering a crystal momentum k + q for the wave function at a single twist . Within canonical twist averaging we can then obtain the neutral gap by However, in contrast to grand canonical twist averaging for the fundamental gap [14], only a finite set of different crystal momenta, q, are accessible in a finite simulation cell. In a grand-canonical version of the neutral gap, in addition to changes of the crystal momentum by q at fixed number of electrons , the number of electrons can vary at different twists with possibly different changes of the crystal momentum. For example, we may increase/decrease the number of electrons at some 1/2 compared to the ground state, ( 1/2 ) = ± 1. Now, all values of q become accessible, and the neutral gap will be the minimum of the canonical one and the fundamental gap, at least to first order, since energy calculations at different twists are essentially independent.

Trial wave functions
Many-body trial wave functions are conveniently written as where k (r|R) may depend symmetrically on all other electron coordinates, R, as it is the case in generalized backflow or neural network orbitals [8-11, 35, 36]. If both, k (r|R) as well as the symmetric factor (R), obey periodic boundary conditions and are invariant against the translation of all electron positions by crystal lattice vectors, we have The crystal momentum of the wave function is then given by the sum over all occupied Bloch vectors, which must be chosen in the grid of wave vectors compatible with . Generalizations to include spin and multideterminant wave functions are straightforward.

QMC results for diamond carbon
We have performed electronic QMC calculations on carbon in the diamond structure at ambient pressure = 1.318 (lattice constant 3.567 Å). We used a Slater-Jastrow trial wave function, which was fully optimized within variational Monte Carlo, including the long-range (reciprocal lattice) contributions. The orbitals in the Slater determinant were taken from DFT calculations with QUANTUM ESPRESSO [37,38] using the LDA functional with a cutoff of 200 Ry for the kinetic energy and of 400 Ry for the electron density and with 12×12×12 k-point grid. QMC calculations have been performed with the QMCPACK code [39]. Burkatzki-Filippi-Dolg pseudopotential [40] was used to remove the core electrons. Diffusion Monte Carlo calculations were performed with 0.01 Hartree time step as in [14]. We used two system sizes: the cubic cell containing eight atoms and a 2 × 2 × 2 supercell containing 64 atoms. All calculations are done such that they can be directly compared with those of the quasiparticle gap of [14] (see also Supplementary material of [14]). Figure 1 shows the DFT-LDA band structure of carbon diamond and the chosen band path in the Brillouin zone. The indirect gap amounts to 3.9 eV. The direct Γ → Γ gap is 5.4 eV and the Γ → , which is slightly larger that the indirect gap, is 4.4 eV. We will further limit our QMC calculations to only two neutral electron-hole excitations: Γ → Γ and Γ → . Table 1 summarizes our DMC results for the neutral gap at a fixed wave vector imposed by promoting an electron from Γ to , and from Γ to Γ. We will mainly focus on Γ to in what follows. The  For Γ → , we can compare our neutral gap with the quasiparticle gap given in the supplementary material of [14]. We find both gaps to coincide within our statistical error. Such a coincidence of neutral and quasiparticle gap for supercells of different sizes has already been found in diamond silicon [18,29].

33701-4
De facto, these results exclude the use of different finite-size extrapolation laws for quasiparticle and neutral gap. Previously, we derived the asymptotic 1/ law for the fundamental quasi-particle gap as an exact result based on a microscopic analysis with the underlying trial wave function [14]. Below we will extend this result to the neutral gap and show that the same 1/ asymptotic behavior is recovered as long as the excitations are extended, as it is the case for our single-determinant wave function.

Finite size effects for neutral gaps: extended vs localized excitations
For the quasiparticle gap, the 1/ behavior was directly related to the asymptotic → 0 value [14] of the structure factor (k) = ⟨ k −k ⟩/ where k = exp[−ik · r ]. In figure 2, we compare the difference of the structure factor of the excited state and the ground state, ( ), from neutral excitation to the ones obtained from electron addition and removal. Although affected by considerable noise, and sensitive to the optimisation, especially of the long-range part of the wave function, our results for the neutral gap extrapolate to a non-vanishing value. The variance of these values at fixed is large, since different angular directions seem to follow the envelope given by the electron addition and removal  curves (green) with opposite slopes. The behavior of the structure factor further supports the same 1/ extrapolation law of the quasiparticle gap also for neutral excitations. Below we analyze the behavior of the structure factor for single and multideterminant wave functions and relate their asymptotic limits to the extended character of the excitation.

Single determinant wave function
Let us consider a single determinant wave function of the form equation (2.6) where the Slater determinant is formed out of Bloch orbitals k labelled by band index and crystal momentum k. In the context of a many-body wave function, two orbitals of crystal momentum k and k+q , → 0, are defined to belong to the same band if they are adiabatically connected in the sense that e iq·r k (r|R) ≃ k+q (r|R), (4.1) and k (r) = k+G (r) is a periodic function in the reciprocal space of the crystal and we have dropped the dependence on R to simplify the notation. This property is of course verified in the case the orbitals are taken from any effectively independent electron theory, e.g., Kohn-Sham DFT orbitals as used in our calculations above. In a normal band insulator with non-vanishing fundamental gap [14], we can assume that the occupied orbitals (those entering the Slater determinant 0 ) consist of completely filled bands, e.g., for each filled band all of the crystal momentum states within the first Brillouin zone are occupied and adiabatically connected with each other. Applying a density fluctuation q with long wavelength, → 0, we then have where we have used in the second line that k+q is already occupied, together with the antisymmetry of the determinant. Density fluctuations are thus suppressed in the insulating state by the gap in the density of states, and we have 0 ( → 0) = 0 in the ground state.
Let us now consider a particle-hole excitation by promoting an electron from the valence band = at k ℎ to the conduction band = with Bloch momentum k , described by a Slater determinant Applying again a density fluctuation, we obtain in the long wavelength limit

33701-6
We then get for the structure factor . (4.5) The Slater determinants entering the numerator are orthogonal, so that we have k ℎ → k ( → 0) = [ + ( → 0) + − ( → 0)] = 2 in the noninteracting limit, (R) = 0, where ± ( ) is the structure factor of a quasiparticle excitation with a single electron added or removed [14]. In general, although (R) will violate orthogonality, the limiting value will not vanish and remains on the order of + ( → 0) + − ( → 0). For simplicity, so far, we have assumed that the limiting value → 0 can be taken at fixed number of electrons, neglecting that for finite systems the structure factor is only known on a finite grid of q. Extending the above calculations to the next order, one can see that the coefficient of order 2 is extensive ∼ . However, almost all of the contributions occur equally for ground and excited states, so that the prefactor of the 2 behavior cancels to a large extent in the difference, ( ) = [ k ℎ → k ( ) − 0 ( )] [and similar for ± ( )]. Therefore, as in the case of adding or removing electrons [14], the difference between the excited and ground state structure factor will remain finite when approaching = 0. Assuming that the asymptotic value of k ℎ → k coincides with + + − from the quasiparticle gap, we can follow exactly the argument given in [14]. Then, the leading-order finite-size corrections will be of order 1/ and proportional to the inverse dielectric constant −1

Localized excitations: multideterminant wave functions
Let us now discuss the case of localized excitations, the simplest possibility is via coherent particlehole superpositions, e.g., with coefficients p . We now get . (4.8) We now see that the asymptotic value of the structure factor will vanish if the coefficients p become a smooth, differentiable function of p in the thermodynamic limit. In this case, ( ) will vanish as 2 , and the finite-size error of the neutral gap will be of order 1/ .

Excitonic effects
As already mentioned above, we have not observed any improvement by the use of multi-determinant wave functions in the case of diamond in the 8 atom cell. Still, due to particle-hole attraction, localization of the excitation will eventually occur giving rise to excitonic effects.
Let us denote by the length scale of localization effects such that | | decays exponentially with for increasing . For simulation cells of extension ≲ 2 , the decay of is not resolvable due to the finite resolution 2π/ in momentum space. The wave function is then indistinguishable from a single 33701-7 determinant one, and we have (π/ ≲ ≲ 2π/ ) ≃ + + − in this case. Only for sizes ≳ 2 the particle hole correlations become effective, ( ≲ π/ ) ∼ 2 .
Assuming a simulation cell with ≪ , we see that finite-size effects will still be dominated by the approximate flatness of the structure factor until very small wave vectors ≲ π/ are reached. Until sizes of order of the localization length are reached, an apparent 1/ will remain. A simple estimate for finite-size effects including excitonic effects is given by where we basically subtract the overshooting of the correction given in equation (4.6) once system sizes of order 2 are reached and neglect any further corrections from the ultimate 2 asymptotics. Let us stress, that size effects described by equation (4.9), correct for excitonic effects of neutral excitations at the onset of the continuum. They do not address size effects of calculations aiming directly at the binding energy of excitons.

Discussion and conclusions
We have performed QMC calculations of the neutral gap in diamond carbon for two supercell sizes. Our values of the neutral gap was found to practically coincide with those of quasi-particle gaps of diamond carbon for both supercells. Such a quantitative agreement of neutral and quasiparticle gap has already been observed in diamond silicon [29], indicating that finite-size effects are not necessarily different in calculations of neutral or charged gaps.
We have given further analytical arguments for the observed finite-size effects. Single determinant wave functions imposing exact crystal momentum will necessarily result in extended, delocalized excitations. The resulting values for the neutral gap will then obey the same leading order 1/ in the thermodynamic limit extrapolation as in the case of quasiparticle gap calculations.
We have shown how multi-determinant wave functions can describe localized excitations for a fixed crystal momentum. The localization ultimately restores a 2 behavior of the structure factor at small wave vectors. This implies the expected faster 1/ convergence of the neutral gap, once the simulation cell exceeds the scale needed to observe localization.
The natural scale of localization is given by the exciton's Bohr radius, ∼ ℏ 2 / 2 where is the effective band mass of particle-hole excitations, e.g., the curvature of the energies involved in the transition. This sets also the scale of simulation cells needed to observe the cross-over from 1/ to 1/ 3 asymptotics in extrapolating neutral gaps. We have shown how finite-size corrections can be adapted to include the expected change of the behavior.
In diamond carbon, a rough estimation, based on the measured particle-hole effective mass ∼ 0.2 0 , 0 being a free electron mass, and the dielectric constant = 5.7 [41], gives an excitonic length scale of ∼ 28.5 Bohr, necessitating a simulation cell roughly twice as large. Using equation (4.9) indicating a possible lowering of ∼ | (2 )|/ = 0.24 eV of the neutral gap due to excitonic effects with respect to the quasi-particle gaps extrapolated with equation (4.6). Such a reduction would also bring the neutral QMC gap closer to experimental values, which however need to be corrected by electron-phonon coupling effects and the use of pseudopotentials for the core electrons, see table II of reference [14].
Discussing the methodology of QMC calculations for the neutral gap, we have included the possibility of twisted boundary conditions. In our QMC calculations on diamond carbon, we have not made use of them, in contrast to previous calculations of the fundamental gap [14]. Changes of the neutral gap with respect to small twists in the boundary conditions can be used to probe the effective band mass of particle-hole excitations, such that estimates for can be obtained without involving knowledge external to the QMC calculations.