Electronic structure of semiconductor nanoparticles from stochastic evaluation of imaginary-time path integral

In the Kohn-Sham orbital basis imaginary-time path integral for electrons in a semiconductor nanoparticle has a mild Fermion sign problem and is amenable to evaluation by the standard stochastic methods. This is evidenced by the simulations of silicon hydrogen-passivated nanocrystals, such as $Si_{35}H_{36},~Si_{87}H_{76},~Si_{147}H_{100}$ and $Si_{293}H_{172},$ which contain $176$ to $1344$ valence electrons and range in size $1.0 - 2.4~nm$, utilizing the output of density functional theory simulations. We find that approximating Fermion action with just the leading order polarization term results in a positive-definite integrand in the functional integral, and that it is a good approximation of the full action. We compute imaginary-time electron propagators in these nanocrystals and extract the energies of low-lying electron and hole levels. Our quasiparticle gap predictions agree with the results of high-precision calculations using $G_0W_0$ technique. This formalism can be extended to calculations of more complex excited states, such as excitons and trions.

Non-perturbative high-precision quantum Monte Carlo (QMC) techniques, such as fixed node, diffusion, and auxiliary field, exist (see, e.g., [11][12][13]). Importantly, in these electronic systems the Fermion sign problem [14,15] is mild enough to allow precise simulations. However, these MC methods are mostly suited for studying the ground state. Work to develop a QMC technique for excited states has started, but the studies so far have been based on the tight-binding approximation and applied to model systems (see, e.g., [16,17]) and to graphene and carbon nanotubes [18][19][20][21][22][23][24].
Here we present a DFT-based comprehensive nonperturbative QMC technique for a semiconductor nanoparticle, where excited states, such as electrons and holes, excitons and trions, can be obtained from the output of the same MC simulation. The system-specific Kohn-Sham (KS) orbitals are used as a basis in the electron action in the path integral representation of the statistical sum. Our results suggest that in this approach there is only a mild Fermion sign problem and evaluation by the standard stochastic importance sampling methods employed in, e.g., lattice Quantum Chromodynamics (QCD) (see [25]  Si 293 H 172 , including low-lying single-particle energies.
In the Born-Oppenheimer approximation the nonrelativistic Hamiltonian for valence electrons is Here ψ α (x) is the electron field operator, α is the spin index, e is the electron charge, V eN (x, R I ) is a pseudopotential, i.e., an effective potential of ions at positions R I felt by the valence electrons [26,27]; A 0 (x) is the scalar potential operator which mediates electron electrostatic interactions. Note that the A 0 terms can be integrated out leading to the standard two-body Coulomb interaction operator. The Kohn-Sham (KS) equation of the orbital-based DFT with a semi-local exchange-correlation functional, such as that by Perdew, Burke and Ernzerhof (PBE) [28], is where φ i (x), i are the KS orbitals and eigenvalues, respectively, and E xc [n] is the exchange-correlation functional, and n(x) is the ground state density of valence electrons [2,28,29]. In general, the state label i may include band number, lattice wave vector and spin label. But in this work we only consider spin-symmetric aperiodic systems so that φ i↑ (x) = φ i↓ (x) ≡ φ i (x). Extending this approach to the case of a periodic and/or spin-polarized system is straightforward.
In order to utilize electronic structure information from the DFT output we introduce a iα , which is the annihilation operator of a Fermion in the KS state |i, α , so that ψ α (x) = i φ i (x)a iα (see, e.g., [30,31]). In terms of a iα the imaginary-time (Euclidean) action corresponding to the Hamiltonian (1) is where β = 1/T e , T e is the electronic temperature in energy units used in the simulation, and The spatial integration is over the simulation box volume; KS indices i and j vary over the range of KS orbitals included in the simulation -the "active window". Fermion fields a iα (τ ) are anti-periodic a i (τ ) = −a i (τ + β), while the A 0 -fields are periodic in time. The grand canonical chemical potential µ is set at the mid-gap which ensures charge neutrality. In order to obtain (3) we added and subtracted eV KS (x) from (2) to the electron lagrangian, and then used ∇ 2 V H = −4πen(x) combined with the shift-invariance of the functional integration over A 0 . The statistical sum is given by The advantage of using KS orbital basis is that these states approximate binding of electrons to the ions and some of the electron interactions. Also, since KS states are labeled by their energy, it is straightforward to only include few KS states near the Fermi level that are relevant to the description of low-energy excitations, which is more efficient than using spatial grid covering the whole simulation cell. Thus, in this approach low-energy excitations of a nanoparticle are described by the KS quasiparticles subject to the static potential V xc ij and interacting via A ij (τ ) exchanges. Strictly speaking, one should proceed using the basis of i δ ij − V xc ij eigenstates. However, setting V xc ij = V xc ii δ ij is an approximation often made in the GW method [6,32]. Also, for the systems we have simulated and for the range of KS states included, V xc ij is strongly dominated by its diagonal entries. So, here we also approximate V xc ij = V xc ii δ ij . Next we perform Grassmann integration over the Fermion variables and expand the resulting action S(A 0 ) in powers of A 0 (τ, x). The terms linear in A 0 cancel, which reflects the system's overall charge neutrality. Retaining the leading non-vanishing term in the expansion yields the following action where is the inverse non-interacting propagator and A iτ,jτ = A ij (τ )δ τ τ . The second term in (7) includes the random phase approximation (RPA) polarization insertion in the KS orbital basis. The approximate action (7) has been used in this work to simulate the nanoparticle electrons. But, in order to evaluate the statistical sum (6) numerically we define the action on a discretized space-time grid. The Lagrangian corresponding to the Hamiltonian in (1) is invariant under time-dependent, spatially uniform U (1) gauge transformations where Λ(τ ) is a function of time [18]. Note that the gauge field action -the last term in (3) -is invariant under (9) [18] [43].
To generate A 0 -field configurations we have used the frequency representation and a spatial grid, where the polarization term was expressed as a function of x, y using φ j (x). In this case A 0 (ω, x) can be used instead of the link variables required on a Euclidean time lattice [33]. However, a τ -KS lattice, where {A 0 (ω, x)} have been converted to the link variables, has been used to compute the observables. While more computationally expensive than a τ -KS lattice, the ωx basis is used since 1) numerical cancellation of the "tadpole" terms in a simulation requires perfect representation of the time derivative operator; 2) the Laplacian in the gauge field action cannot be represented accurately with the few KS orbitals included in a reasonably sized active window. Then, the action is where ω k = 2πkT e is a bosonic Matsubara frequency, and where a l , l = x, y, z, are the lattice spacings and n i = (exp [β˜ i ] + 1) −1 . Gauge invariance requires that each included A(ω) is coupled to at least one pair of fermion modes a † (ω 1 ), a(ω 2 ) with ω = ω 2 − ω 1 . Therefore, in the simulations where a finite frequency cutoff was used the expression for S k (x, y) was modified accordingly.
Importantly, S k (x, y) from (11) is real and symmetric under x ↔ y, which is due to the basic properties of KS eigenfunctions (see, e.g., [34]). Then the action S 2 (A 0 ) from (10) is non-negative, which suggests that in this approach meaningful simulations are possible without resorting to sign-suppression techniques (see, e.g., [35,36]). The size of neglected higher order terms in the action expansion will be discussed later.
The action (10) is quadratic in A 0 . Therefore, importance sampling is done by diagonalizing S(ω k , x, y) for each ω k which results in where λ i (ω k ) > 0, v(0) i ≡ 0, N x is the number of spatial grid points. Then one generates random values for the u(ω k ) i , v(ω k ) i variables that are distributed normally according to the corresponding eigenvalue λ i (ω k ), and changes the basis to get the A 0 (τ, x) configurations. For all the nanocrystals simulated here N 10 3 configurations have been found to be sufficient to obtain statistically significant results.
The observables considered here are the electron and hole quasi-particle energies, which have been extracted from two-point propagators, which are the averages of the matrix elements of the propagator matrix corresponding to the state i and time slices τ 1 and τ 2 , in the long-time limit. For a particle state i > HO, where HO labels the highest occupied orbital, it is E g β E g τ 1, where E g is the quasiparticle gap. Then where δ = β/N t , N t is the number of time grid points, A ij (τ ) is defined in (5),˜ i is defined in (8). The last line in (13) shows the behavior of the correlator at large times on the fully interacting excitation energy E i (C is an irrelevant coefficient). We fit for this energy and estimate the statistical error of E i via the bootstrapping technique [37]. A full error analysis, which would include an assessment of systematic errors due to, e.g., fit-window sizes and fitting functions, finite lattice spacings, etc., is left to future work, though we do not expect the uncertainties quoted here to change significantly. We have used time grids with δ = 0.025 eV −1 for all systems after checking that KS energies can be accurately extracted from the propagators in the non-interacting cases. Frequency cutoffs were chosen so that ω max δ 1. The procedure to extract hole excitation energies is analogous.
DFT simulations of the atomistic models of the nanocrystals (such as Si 293 H 172 shown in Fig. 1), including geometry relaxation, have been done using Quantum Espresso DFT program with the PBE exchangecorrelation functional [38]. Norm-conserving pseudopotentials [39] have been used ensuring that the KS orbitals are orthonormal. Kinetic energy cutoff, which determines lattice spacings a i , i = x, y, z, has been set to 340.1 eV, which is the same as in [32], resulting in a i 0.05 nm. Models of the nanocrystals ranging in size from 1.0 to 2.4 nm were placed in the periodic cubic simulation boxes with about 1 nm of vacuum between the surfaces in order to prevent spurious interactions between periodic images. The number of KS orbitals included in the simulations have been chosen so that imax − HO HO+1 − imin ≥ 1.5E P BE g , where i max , i min are the highest and the lowest included KS orbital labels and E P BE g = HO+1 − HO is the noninteracting gap. The corresponding number of states above/below Fermi level included in the simulation varied from 36 to 96 as the system's size increased.
For the nanoparticles considered in this work V xc ii shifts are sizable and tend to shrink the non-interacting gap E P BE g . This would require lowering T e in order to maintain E P BE g T e required for the tadpole term cancellation, which would significantly increase the computational expense. So, here we have treated V xc ii as selfenergy corrections, i.e., we have simulated with˜ i = i − µ and then subtracted V xc ii from the resulting singleparticle energies.
The calculation results are shown in Tables I and II. In order to check the size of the terms neglected in the approximate action (10) we used the A-configurations gen-  (1) 1.0865 (7) 1.073 (3) Whiles 2 is only the leading non-vanishing term in the expansion of s, we have found that in all cases s 2 /|s| is close to one. (See Table I.) This suggests that the full Fermion action can be reasonably approximated by just the leading term. The average of the action's phase ϕ, where s = |s|e iϕ , is close to zero, which suggests that the sign problem in these systems is mild. Shown in Table II are the quasiparticle gap predictions in these nanocrystals. Our results agree with the results of highprecision calculations using G 0 W 0 for the same nanocrystals [32]. Low-energy particle and hole levels in Si 293 H 172 are shown in Fig. 2.
In conclusion, we have performed initial steps toward development of a first-principles high-precision Monte Carlo technique for the excited states of a semiconductor nanoparticle, which utilizes KS orbital basis in the imaginary-time functional integral for the electrons. We  find that approximating Fermion action with the leading order RPA polarization term in the expansion in powers of A 0 leads to a positive definite integrand in the statistical sum and that it is a reasonable approximation to the full action; s 2 /|s| − 1 can be viewed as a source of systematic error. So, our results suggest that in this approach these systems have only a mild Fermion sign problem. Obvious improvements to our approximate calculations would be to use the full action instead of our quadratic approximation in (7), which could be done via re-weighting [40,41] or more advanced sampling techniques (e.g., hybrid Monte Carlo [42]). Work to develop technique for other excited states, such as excitons and trions, is in progress.