Investigation of Ionization Potential in Quantum Dots Using the Stratified Stochastic Enumeration of Molecular Orbitals Method.

The overarching goal of this work is to investigate the size-dependent characteristics of the ionization potential of PbS and CdS quantum dots. The ionization potentials of quantum dots provide critical information about the energies of occupied states, which can then be used to quantify the electron-removal characteristics of quantum dots. The energy of the highest-occupied molecular orbital is used to understand electron-transfer processes when invesigating the energy-level alignment between quantum dots and electron-accepting ligands. Ionization potential is also important for investigating and interpreting electron-detachment processes induced by light (photoelectron spectra), external voltage (chemiresistance), and collision with other electrons (impact ionization). Accurate first-principles calculations of ionization potential continue to be challenging because of the computational cost associated with the construction of the frequency-dependent self-energy operator and the numerical solution of the associated Dyson equation. The computational cost becomes prohibitive as the system size increases because of the large number of 2particle-1hole (2p1h) and 1particle-2hole (1p2h) terms needed for the calculation. In this work, we present the Stratified Stochastic Enumeration of Molecular Orbitals (SSE-MO) method for efficient construction of the self-energy operator. The SSE-MO method is a real-space method and the central strategy of this method is to use stochastically enumerated sampling of molecular orbitals and molecular-orbital indices for the construction of the 2p1h and 1p2h terms. This is achieved by first constructing a composite MO-index Cartesian coordinate space followed by transformation of the frequency-dependent self-energy operator to this composite space. The evaluation of both the real and imaginary components of the self-energy operator is performed using a stratified Monte Carlo technique. The SSE-MO method was used to calculate the ionization potentials and the frequency-dependent spectral functions for a series of PbS and CdS quantum dots by solving the Dyson equation using both single-shot and iterative procedures. The ionization potentials for both PbS and CdS quantum dots were found to decrease with increasing dot size. Analysis of the frequency-dependent spectral functions revealed that for PbS quantum dots the intermediate dot size exhibited a longer relative lifetime whereas in CdS the smallest dot size had the longest relative lifetime. The results from these calculations demonstrate the efficacy of the SSE-MO method for calculating accurate ionization potentials and spectral functions of chemical systems.

Introduction

Ionization potential (IP) ω (or ionization energy) is defined as the energy needed to remove an electron from a chemical system.The ejection of the electron can be facilitated using incident photons, scattering by high-energy electrons, or applying a strong electric field. As one of the fundamental properties of a material, IP has relevance to the areas of photovoltaics, mass spectrometry, photoelectron spectroscopy, electrochemistry, photocatalysis, and light-induced electron-transfer processes. The IP of atoms, molecules, and various chemical compounds is a quantity of interest when performing X-ray spectroscopy. Recently, it has been found that small gas-phase polyatomic molecules with a heavy atom, such as iodomethane, when bombarded with hard X-ray pulses display surprisingly enhanced ionization relative to an individual heavy atom with the same absorption cross-section.[1] Following the excitation of an electron from an inner orbital of one atom, another electron from a higher energy orbital of the same atom can occupy this inner orbital. Instead of undergoing de-excitation, the newly excited electron can transfer energy via photon emission by ionizing an electron from an outer orbital of a neighboring atom. This process is called interatomic Coulombic decay (ICD).[2,3] The IPs of both the inner and outer orbitals of molecules, dimers, and clusters influence which de-excitation mechanism occurs in a highly excited neutral or a highly excited ionized state of these types of systems. Knowledge of the IPs of these systems’ inner and outer orbitals can assist in the prediction of which de-excitation mechanism is likely to occur in a system of interest.[4−9] High precision IP measurements in atoms and molecules provide important information about electron–electron correlation and serves as a benchmark for the development and testing of theories. Ionization potential also is a quantity of interest in biological systems. For example, ionizing radiation causes permanent heritable DNA damage,[10] and the IPs of nucleic acid tautomers are quantities of interest, due to the role of these molecules in cancer.[11]

In quantum dots (QDs) and nanomaterials, knowledge of the IP for ionization from HOMO serves as an important metric for quantifying electron-transfer rates.[12−18] Photoejection of electrons by X-ray and UV radiation has been used to study valence-band states in QDs.[19] Cyclic voltametry has also been used to calculate HOMO energies in QDs and IPs.[20] Transient photoemission in two-photon experiments has provided information on the energy levels of unoccupied orbitals.[21] Knowledge of the relative positions of the HOMO and lowest unoccupied molecular orbital (LUMO) levels of a QD with respect to the surface ligands is an important factor in extraction of a hot carrier from the QD.[19]

In molecular quantum chemistry, the simplest approximation of IP (denoted as ω0) is given by the Koopmans’ approximation, where the exact IP is approximated as the negative of the orbital energies.Koopmans’ treatment utilizes the Hartree–Fock (HF) approximation to obtain a single N-electron Slater determinant from which an electron from one of the occupied states is annihilated, as shown in eq . Although Koopmans’ theorem is limited by the use of single Slater determinants and does not account for orbital relaxation and electron–electron correlation effects, it still provides an acceptable first approximation for the IP of a system of interest.[22] Going beyond Koopmans’ approximation by including the effect of electron−electron correlation can be achieved in a variety of different ways such as with electron-propagator methods,[23−31] algebraic diagrammatic construction (ADC),[32,33] equation-of-motion coupled-cluster (IP-EOM-CCSD),[34−39] many-body perturbation theory (MBPT),[40,41] GW method,[42] correlated-orbital theory (COT),[43] and time-dependent density functional theory (TDDFT).[44−46] The IP of the HOMO has a special significance in DFT because of Janak’s theorem and plays a prominent role in the development and testing of DFT functionals.[47] Without loss of generality, the many-body correction to the IP can always be written as,where Δω accounts for all of the correction terms missing from the Koopmans’ approximation; post-HF methods mentioned earlier offer different approximations and formulations for calculating Δω. However, efficient first-principles calculation of Δω for large chemical systems continues to be challenging and is an active field of research. Using the electron-propagator method, Ortiz and co-workers developed a series of approximations that offer an order-by-order treatment of electron correlation to the many-body correction for IPs.[23−25,48,49] Open-shell systems possess additional complexities compared to their closed-shell counterparts, and the spin-flip EOM-IP approach has been used to treat open-shell systems.[50] In the GW formulation, the projective eigen decomposition of the dielectric screening (PDEP) algorithm has been used for calculating the quasiparticle gap of QDs.[42] In recent works, methods using stochastic techniques have been demonstrated to achieve the low-scaling needed for applications to large chemical systems. Specifically, the use of stochastic orbitals in the stochastic Green’s function method developed by R. Baer can be used for the calculation of IPs.[51] A different strategy of combining the Laplace-transformed expression of the self-energy operator with a real-space Monte Carlo integration scheme developed by Hirata and co-workers has been used for the calculation of IPs at the second-order MP2 level.[52] The approach was recently extended for ground-state MP4 level.[53] In a related work by Li et al., the Laplace-transformed MP2 approach has been combined with the density-of-states approach to reduce the overall computational cost of the MP2 calculation.[54] This approach demonstrated the effectiveness of using the intrinsic degeneracies present in chemical systems to reduce the overall cost of MP2 calculations.

In this article, we present the stratified stochastic enumeration of molecular orbitals (SSE-MO) method for the efficient computation of the IP through the iterative solution of the Dyson equation. The SSE-MO method was originally inspired by the 2013 paper, “Stochastic Enumeration Method for Counting NP-Hard Problems” by Rubinstein.[55] The original stochastic enumeration by Rubinstein was based on the importance sampling scheme. In the field of computer science, the stochastic enumeration technique has been applied to traversing deep tree structures and implementing backtracking algorithms.[55] For the SSE-MO method, we have combined stochastic enumeration with stratified sampling to perform the necessary summations over a direct-product space of molecular orbital indices and 6D Cartesian coordinate space. The reduced computational cost of the SSE-MO method allowed us to investigate the full frequency-dependent pole structure of the 1-particle Green’s function by an iterative solution of the Dyson equation. In this work, the SSE-MO method was used to investigate IPs of PbS and CdS QDs. The motivation for the SSE-MO method comes from the fact that not all terms contribute equally to the overall self-energy operator. A similar observation has also been made for the density-of-state MP2 method.[54] For example, for the Pb4S4 system, the contribution for each term as a function of the hole-index, (i), is presented in Figure . Without loss of generality, the self-energy operator can be written as the sum of the 2-particle 1-hole (2p1h) and 1-particle 2-hole (1p2h) terms as shown in eq .It is seen that some terms contribute more than other terms. The SSE-MO method aims to distribute the computation effort used to calculate of the overall self-energy operator in proportion to the contribution of each of the terms and their associated errors. Here, we provide a general formulation of stratified stochastic tensor contraction (section ), details about the construction of the self-energy operator (section ), formulation of stochastic stratified sampling in MO space (section ), performing low-discrepancy sampling (section ), construction of the control-variate functions for the real-space integrals (section ), performing joint MO-space real-space stratified sampling (section ), the iterative solution of the Dyson equation (section ), and the calculation of pole-strengths of the Green’s function and calculation of higher-order derivatives of the self-energy operator (section ). Benchmark calculations on well-studied systems and new results on PbS and CdS QDs are presented in section . Comparisons to existing methods and future directions are discussed in section .

Figure 1

Percent contribution of terms A2p1h and A1p2h (defined in eq ) to the total self-energy as a function of hole index (i) for Pb4S4.

Theory

Stratified Stochastic Enumerated Tensor Contraction

As an introduction to the application of SSE for the calculation of the self-energy operator, we present the SSE approach for performing a general N-index tensor contraction. We start by considering the following general tensor contraction S = Tr{ABCD}.One situation for which this type of tensor contraction is encountered is when the integration of a 4-point kernel in real space ⟨A(r1)B(r2)C(r3)D(r4) G(r1, r2, r3, r4)⟩ on a spatial grid with N point per dimension is being performed. This tensor contraction has N12 terms and a simple sequential evaluation will require N12 terms. For the stratified stochastic enumeration approach, we will first define a composite index K such that K = 1, ..., N12. The composite index, K, uniquely maps each ordered set of indices (i1, i2, i3, ..., i12) to an integer in 1, ..., N12.Using the composite index K, we can define the summation as displayed in the following equation.Next, we divide the entire range of K into N nonoverlapping segments Nseg = N. The number of terms in each segment is N = N11. The summation over K can be written in terms of the segmented summation.The partial averages are defined as follows.The total sum can be written as displayed in the equation below.In SSE, the sequential segment average, X̅, is approximated using the stochastic average,where the subscript in denotes “sampling-without-replacement”. For any segment “p”, the SSE average approaches the sequential average as M → N.The SSE estimate of the total summation is presented in the equation below.The allocation of the sampling points for each segment is proportional to the variance in the SSE segment average X̅SSE.The SSE approach is based on stratified sampling, which has been used extensively for reducing sampling error in Monte Carlo calculations[56−58] and a brief description stratified sampling is presented in the Supporting Information (SI). The SSE method is not restricted to square tensors and can be applied to rectangular tensors as well. We recommend a row-major composite indexing scheme. For an index vector (i1, i2, i3, ..., i), where D is the dimension of the tensor and each index (i, d = 1, D) is in the range (i = 1, ..., N), the row-major composite index K can be calculated using the following expression.

Second-Order Dyson Equation

In this work we are interested in calculating the IPs of chemical systems. The component of the 1-particle Green’s function, G(ω), that contains information about the IPs can be expressed in the Lehman representation as follows.[59]where, G is the matrix element of the matrix representation of the operator in canonical Hartree–Fock (HF) orbital basis {χ},{ a†, a,} are creation and annihilation operators defined with respect to the HF orbitals, {E0, Ψ0} are the exact ground state energies and wave function, and H is the electronic Hamiltonian. By inserting a complete set of projectors, 1 = ∑ |Ψ⟩ ⟨Ψ|, we obtain the following equation.As shown in eq , the Green’s function’s poles correspond to the vertical IPs of a many-electron system.The quantity |A|2 is the residue of the pole and is known as the pole strength.The limit η → 0+ in eq is traditionally associated with this expression because of its use in performing the Fourier transform from the time-domain to frequency-domain and will be suppressed in the rest of the derivation. Analogous to the many-body Green’s function, the uncorrelated HF Green’s function G0 is given by the following expression.[59]which immediately simplifies to following diagonal representation.[59]where i, j = 1, Nocc are indices for occupied orbitals and ϵ is the orbital energy. The above equation recovers the Koopmans’ approximation to the IPs, which is defined as EIPKoopmans = −ϵ.

In the frequency representation, the relationship between the correlated 1-particle Green’s function, G, and the uncorrelated Green’s function G0 is given by the well-known Dyson equation.[60]Here, Σ is the self-energy operator. The relationship between the correlated and uncorrelated Green’s function can be derived using various techniques including time-dependent perturbation theory, time-independent perturbation theory, coupled-cluster theory, the configuration interaction method, and electron-propagator methods.[26,34,59,61,62] The above operator equation can also be presented in various representations such as plane-waves, real-space grids, and canonical HF orbitals.[42,63,64] In this work, we use canonical HF orbitals to represent the Dyson equation.[60] To facilitate the calculation of the poles, it is useful to express eq in terms of inverse operators by multiplying G–1(ω) from left and G0–1(ω) from the right.Once G(ω) is determined for an appropriate set of values of ω, the poles can be observed by constructing a plot of G(ω) versus ω. Approximating the total self-energy operator by diagonal representation,allows for analytical inversion of the Dyson equation into the following simplified expression,[60]where ω0 is the orbital energy of occupied orbital (i) and ω = −EIP. We have used the second-order approximation to the self-energy operators, which in the canonical MO basis is defined as,[60]where i, j, and k indicate occupied spin orbitals and a and b indicate virtual spin orbitals.

Using the restricted Hartree–Fock (RHF) formulation, the correction to the orbital energies using the second-order self-energy expression can be written as,where the RHF expressions for the self-energy terms are defined as,Here, we have used the following compact notation for the energy denominators,and the r12–1 matrix elements are defined using the chemist’s notation for the indices.Next, we will develop the SSE-MO approach for evaluating the self-energy operator.

Stratified Stochastic Enumeration of Self-Energy

We begin by defining a set of ordered integers (i, a, b),which contains all the possible combinations of indices that occur in 2p1h self-energy expression. We will use the composite index K = (i, a, b) to enumerate this ordered set of integers. The size of set is given as,Using this notation, we can define a general form of the 2p1h self-energy term as follows.where,In eq the summation is performed sequentially for all terms. In the stochastic enumeration (SE) approach the sequential sum is replaced by a stochastic summation. We define a new operator Σ̃ which is defined as follows,and where K is sampled from the set . This sampling is performed without replacement and the notation is used to emphasize this procedure (sample-without-replacement). The NsampleMO is sample size and bounded from above by Kmax2p1h. In the limit when the sample size approaches Kmax2p1h the following limiting condition is satisfied.Simple stochastic enumeration will involve performing the sampling over multiple runs and averaging the final results.The variance is defined as follows.We expect the variance to disappear when NsampleMO approaches Kmax2p1h,To reduce the variance of the overall calculations, we introduce stratification in the sampling procedure. This is achieved in two steps. First, the set is decomposed into a union of nonintersecting subsets,where the subsets are nonoverlapping.The number of elements in subset M is denoted as Kmax2p1h,M.In the second step, X̃2 is calculated using summation over all the subsets.The stratified stochastic enumeration of the MO indices (SSE-MO), which uses stochastic enumeration for segment sampling, is described by the following equation.To write the expressions in compact notation, we introduce the following for stochastic summation.Using this notation, we can write the following expression.A similar treatment is performed for the 1p2h terms. The combined result for the total self-energy operator is given as,

Calculation of Optimal Sampling Points for MO-Space Stratified Sampling

To calculate optimal sampling points, we define the segment average as,which allows us to write the following expression.It is important to note that Ỹ2p1h,M is a stochastic variable for which the average value, Ỹ2p1h,M,avg, can be obtained by sampling over multiple runs.The variance is defined as,The variance of Ỹ goes to zero as NsampleMO,M approaches Kmax2p1h,M,To distribute the sampling points optimally, we define the following weight factor,In addition to that we also define another weight factor that depends on the magnitude of the terms,The sampling is performed in batches, and the variance is updated after completion of a batch. The number of sampling points for each segment has the general form of,where all segments get Nbase sample points per batch irrespective of the segment average and variance. Nopt indicates the number of additional sample points that are distributed in a manner that is proportional to the normalized weights, which depend on the segment average and variance.

Low-Discrepancy Sampling without Replacement Using Quasi-Monte Carlo Method

It is important to note that obtaining is a correlated sampling process. It is intrinsically non-Markovian and depends on the entire history of the string of previously generated indices. One way to achieve this in discrete integer space is by performing self-avoiding random walks. However, sampling in the self-avoiding random walker is local in nature, and therefore is not ideal for variance reduction in each segment. Here we use quasi-Monte Carlo sampling and a low-discrepancy integer sequence to perform sampling within each segment.[56−58] The linear congruent generator for low-discrepancy quasi-random numbers is modified for generation of integer sequences.[56] The sampling index for a segment M is defined as K( and can have values in the range [1, ..., Kmax(]. The exact value of Kmax( for each segment is known at the start of the calculation and is a consequence of the stratification procedure described in subsection . Associated with each segment are two integer random numbers which we define as q( and r(. The variable q( impacts the discrepancy of the points and is an integer random number chosen randomly from the interval q( ∼ [10, 50]. Using q(, we define the following sequence from which r( is selected randomly.Using q( and r(, the low-discrepancy sequence is defined as,where is the floor of the ratio.

Control Variate for Monte Carlo Evaluation of Two-Electron Integrals

In this section, we extend the stochastic procedure developed in the previous section for numerical evaluation of the two-electron integrals. Monte Carlo evaluation of the two-electron integrals is not a requirement for implementing the SSE-MO method, and the SSE-MO procedure described in section can be used whenever the two-electron repulsion integrals, V, are available in MO representation. However, for large systems, it is computationally advantageous to avoid the AO-to-MO two-electron transformation and instead, to numerically integrate directly in the MO representation using the Monte Carlo scheme.

Associated with each MO pair function, ψ(r)ψ(r), we define a control-variate function, ψ(r). The control variate function must satisfy two important features. First, ψ(r) must be nonfactorizable as a product of functions that depend only p and q indices.Second, the two-electron integrals, [ψcv(r1)|r12–1|ψcv(r2)], must be known analytically. Adding and subtracting the control variate function, we express the MO product function as follows.Here, α is the control variate and d(r) is the difference function. Using the above expression, the two-electron integral can be expressed as,where,and,The control variate, α, is defined as the quantity that minimizes the following weighted-variance function.Here, Sspace is a set of sampling points in 3D space from which r is drawn at random. The control-variate functions, ψcv, are represented by Gaussian functions. For p = q, a single Gaussian function is used and for p ≠ q a linear combination of two Gaussian functions is used.This form of the control variate function guarantees the orthonormality conditions for the MOs.The widths and the centers of the Gaussian functions are determined using a moment-matching condition. The weighted moments for any pair of molecular orbitals are calculated as,and the moments for the control-variate functions are obtained analytically.The coefficients for the control-variate function are obtained by performing a steepest-descent search on the following loss-function.The equations eq and eq completely define the control-variate function and are used to calculate the Vcv term.

Monte Carlo Evaluation of Two-Electron Integrals Using Real-Space Stratified Sampling

Since Vcv is analytical, only the D terms are calculated numerically using the stratified Monte Carlo procedure. We use a combination of ratio estimator, control variate, and stratified sampling techniques to efficiently and accurately evaluate the MO integrals. For calculating D, we define the following two-electron kernel function,where,Associated with each D integral, we define a control variate fcv as,where the integral of the control variate is 1 for all values of p, q, s, t.The calculation of D requires evaluation of the six-dimensional integral over all space. Traditionally, Monte Carlo integration is performed over an N-dimensional unit cube by transforming the integral range from [−∞, + ∞] to [0, 1]. One approach to achieve this is by using the following transformation,However, this procedure introduces singularity in the form of the Jacobian J(t) in the integration kernel. In this work, we use a finite-grid approximation to evaluate the integral over a finite volume,where the limits selected are large enough so that ϵ ≤ 10–5 for all MOs. Using fcv, we define the following ratio estimator for Monte Carlo evaluation of the D integral.The averages ⟨Tr12–1⟩ and ⟨fcv⟩ are defined as follows.We introduce stratification in the sampling of points in real-space by dividing the entire space into a set of Nsegspace nonoverlapping regions with identical volumes.The stratified sampling estimate of the averages is then defined as,Nsamplespace,M,pqst is the number of sampling points associated with the spatial segment, M for indices p, q, s, t. Similar to the stratification strategy in eq , the sampling points for each segment are proportional to the variance of the integral kernel in that segment.We use common-random-number (CRN) sampling for sampling within a segment. The CRN method has been used extensively for reducing variance[56−58] and a brief summary of the method is presented in ref (65). In the evaluation of the integral, this means that at any point in time, if a random number η1 is used in the evaluation of T, then that same random number is used for evaluation of fcv in the same segment. To emphasize this usage, we use superscript “CRN” (common random number) in the following expression.

Iterative Solution of the Dyson Equation

Combining the results from section and section , we can express the full self-energy operator as the sum of two terms,The first term Σcv(ω) depends only on control-variate functions and is evaluated analytically.The second term contains the different terms and is defined as follows.

The single-shot determination of the self-energy operator is performed by evaluating the self-energy at the HOMO energy.The full iterative solution of the Dyson equation is obtained by evaluating Σ(ω) for a range of ω and then finding the point where,The SSE-MO method allows for a third approximation for ω. We can solve the Dyson equation iteratively using Σcv(ω)and then include the correction from ΔΣ.

Calculation of Derivatives

The first derivative of the self-energy, with respect to ω, is useful for locating the poles of G(ω), and is defined as follows.The higher-order derivatives of the self-energy operator can then be obtained from the higher-order powers of the denominator.In the SSE-MO method, the derivative of Σ is obtained by replacing the energy denominator with the higher powers of the denominator,The contributions from Σcv are obtained analytically, and the contributions from the ΔΣ are obtained from the stochastic enumeration procedure described above. Because the derivative does not impact r12–1, which is present in expression for the self-energy, the calculation of the derivatives can be performed concurrently while the calculation of ΔΣ is being performed. This approach was used for the construction of the spectral function and is presented in section .

Computational Details

The SSE-MO method was applied to investigate the ionzation potentials of PbS and CdS QDs. In addition, benchmark calculations for Ne, H2O, and CH4 were also performed. The single-particle states and energies were obtained from HF calculations using the 6-31G* basis for Ne, H2O, and CH4 and the LANL-2DZ ECP basis for the quantum dots. These HF calculations were performed using the TERACHEM electronic structure package. The SSE-MO calculations were performed by dividing the MO-index space into Nocc number of segments. The 6D Cartesian space was divided into 100 nonoverlapping regions. A total of Nsample ∼ 109 sampling points were used for calculating the self-energy at each value of ω and the sampling points were distributed using the stratification strategy described earlier. We used the relative standard deviation, σrel, also known as coefficient of variance for defining the convergence criteria for the calculated IPs (eq ).In this work, we enforced σrel < 10–2 to be the criteria for convergence for each segment.

Results

10-Electron System

For benchmarking and testing, the SSE-MO method was used to calculate the IPs of Ne, H2O, and CH4. The results for these chemical systems are presented in Table . The IPs were calculated using both single-shot and iterative solution of the Dyson equation and the results between the two approaches were found to be very similar to a maximum difference of 0.17 eV. In all cases, the SSE-MO results were found to be in good agreement with the previously reported results.

Table 1

Ionization Potentials (eV) of Ten Electron Systems: Comparison with Benchmark Literature Values

system	Koopmans’	single-shot solution	iterative solution	lit. value[66,67]	IP-EOM-CCSD(T)[68]
CH4	14.86	13.94 ± 0.03	13.95 ± 0.03	13.91	12.76
Ne	22.59	21.47 ± 0.06	21.38 ± 0.07	21.13	20.98
H2O	13.56	10.44 ± 0.13	10.61 ± 0.08	10.74	11.37

Ionization Potential of PbS and CdS Quantum Dots

The SSE-MO method was applied to Pb4S4, Pb44S44, Pb140S140, Cd6S6, Cd24S24, and Cd45S45 and the ionization potentials from the single-shot and iterative solution of the Dyson equation are presented in Table and Table , respectively. We note that the calculated IPs are vertical ionization potentials and do not include contributions from the quantum mechanical treatment of nuclear degrees of freedom. Figure illustrates the graphical verification of the self-consistency of the iterative procedure for Pb140S140. We observe that the curve for Σ(ω) + ω0 versus ω intersects with the curve for ω versus ω at that the value of ω for which the diagonal approximation to the Dyson equation converges. The frequency dependence of the 1-particle Green’s function was evaluated near the poles and is presented in Figures , 4, and 5.

Table 2

Self-Energy and Ionization Potentials (eV) of PbS and CdS Quantum Dots from Single-Shot Solution

system	Koopmans’	self-energy from single-shot solution	IP from single-shot solution
Pb4S4	8.28	0.65	7.63 ± 0.05
Pb44S44	7.13	0.22	6.91 ± 0.04
Pb140S140	6.91	0.09	6.82 ± 0.05
Cd6S6	5.25	0.42	4.837 ± 0.04
Cd24S24	6.25	0.09	6.16 ± ϵ < 0.01
Cd45S45	6.09	0.23	5.86 ± 0.02

Table 3

Self-Energy and Ionization Potentials (eV) of PbS and CdS Quantum Dots from Iterative Solution

system	Koopmans’	self-energy from iterative solution	IP from iterative solution
Pb4S4	8.28	0.61	7.66 ± 0.01
Pb44S44	7.13	0.17	6.96 ± 0.01
Pb140S140	6.91	0.28	6.71 ± 0.08
Cd6S6	5.25	0.41	4.84 ± 0.01
Cd24S24	6.25	0.08	6.16 ± ϵ < 0.01
Cd45S45	6.09	0.22	5.87 ± 0.01

Figure 2

ω (ordinate) versus ω (abscissa) displayed as the curve labeled ω. The curve labeled ω0+Σ displays the relationship between the HOMO energy + the self-energy (ordinate) and ω. The value of ω at which these two curves intersect is equivalent to the value of ω for which the diagonal approximation to the Dyson equation converges.

Figure 3

Poles of G0(ω) and G(ω) for the Pb4S4 system. G(ω) (ordinate) versus ω (abscissa) is labeled as G(ω) in the legend. The curve labeled G0(ω) displays the relationship between G0(ω) (ordinate) and ω (abscissa).

Figure 4

Poles of G0(ω) and G(ω) for the Pb44S44 system. G(ω) (ordinate) versus ω (abscissa) is labeled as G(ω) in the legend. The curve labeled G0(ω) displays the relationship between G0(ω) (ordinate) and ω (abscissa).

Figure 5

Poles of G0(ω) and G(ω) for the Pb140S140 system. G(ω) (ordinate) versus ω (abscissa) is labeled as G(ω) in the legend. The curve labeled G0(ω) displays the relationship between G0(ω) (ordinate) and ω (abscissa).

When compared to G0(ω), the poles of G(ω) were found to have higher values of ω indicating that for these systems, inclusion of electron correlation effects resulted in a lower IP than Koopmans’ IP values. Comparison between G and G0 shows that inclusion of electron correlation in IPs becomes more important for larger dots. Comparison between the single-shot versus iterative solution of Dyson equation also exhibits similar trends, where the need for iterative solutions become more important for larger dots.

Single-Pole Approximation to the Spectral Function

We define the single-pole approximation to the spectral function as,where the subscript “sp” in Asp denotes that we are looking at the form of the spectral function near the pole (ω= ϵ). The imaginary part of the self-energy operator can be approximated from the first derivative of the self-energy, with respect to ω, as described in section . For example, the imaginary part of the following quantity I(ω),is given by,In Figure , the ratio Awp/Awpmax is plotted as a function of ω/ωopt for the three PbS QDs.

Figure 6

Ratio Awp/Awpmax (ordinate) plotted as a function of ω/ωopt (abscissa) for a series of PbS quantum dots. Awp/Awpmax is the ratio of Asp(ω) and the maximum value of Asp(ω). ωopt is the value of ω for which convergence of the diagonal approximation to the Dyson equation is achieved.

The line width of the plot was found to be narrowest for the Pb4S4 and broadest for Pb44S44. This feature indicates that the relative lifetime of the quasi-hole in the intermediate dot size (Pb44S44) is longer than the other dots in the series. Similar analysis for the CdS QDs in Figure revealed that the line width decreases with increasing dot size. The results from the spectral analysis highlight the importance of including frequency dependency in the self-energy operator. The plots also demonstrate the impact of many-body correlations in the these systems. Specifically, in the absence of electron–electron correlation, the limit σ → 0 will reduce the plots to a Dirac delta function.

Figure 7

Ratio Awp/Awpmax (ordinate) plotted as a function of ω/ωopt (abscissa) for a series of CdS quantum dots. Awp/Awpmax is the ratio of Asp(ω) and the maximum value of Asp(ω). ωopt is the value of ω for which convergence of the diagonal approximation to the Dyson equation is achieved.

Discussion

Correlated Sampling in the Combined Cartesian and Molecular Orbital Index Space

The main philosophy of the SSE-MO method is to perform correlated sampling in a joint real-space and occupation-number space (Table ). Assuming a discretization of 100 points per Cartesian coordinate, the total number of points needed for exhaustive sampling is in the range of 1016 to 1020 as shown in Table . However, not all spatial components of all the molecular orbitals contribute equally and uniformly to the calculation of the self-energy. There are certain combinations of MOs whose form in specific regions of the Cartesian space correlate strongly with the error in the self-energy calculations. Through the use of a two-step stratified sampling scheme in both Cartesian and MO space, the SSE-MO method provides a systematic and adaptive procedure to identify the important contributors. We have used a combination of ratio estimator, control-variate, and stratified sampling techniques for the efficient and accurate evaluation of the MO integrals. The key quantity that implements and controls this concept is the Nsamplespace,M,pqst term. This term represents the number of spatial sampling points for the Mth spatial segment for the correction term D associated with indices p, q, r, s and depends on both the spatial and MO indices. The total number of sampling points is given by the following expression.As shown in eq , this number was directly obtained from the variance of the integral kernel, which also includes the contribution from the r12–1 operator. Note that these sampling points were not used to evaluate the full r12-integral kernel, but instead were used to evaluate only the component of the full r12-integral kernel not included in the control-variate expression. The Cartesian space sampling for each spatial segment was performed using simple Monte Carlo sampling. This process can be enhanced by using low-discrepancy random numbers, which is a quasi-Monte Carlo approach. We expect that using the quasi-Monte Carlo approach will accelerate the overall calculation process.

Table 4

Total Number of Sampling Points in the Combined MO-Cartesian Space Assuming 100 Points Per Cartesian Coordinate

system	N2p1h	N1p2h	NspaceMO	NspaceMO × Nspace6D
Pb4S4	3.87 × 104	1.76 × 104	5.63 × 104	5.63 × 1016
Pb44S44	5.15 × 107	2.34 × 107	7.50 × 107	7.50 × 1019
Pb140S140	1.18 × 109	6.37 × 108	1.82 × 109	1.82 × 1021
Cd6S6	7.02 × 105	3.32 × 105	1.03 × 106	1.03 × 1018
Cd24S24	4.50 × 107	2.13 × 107	6.62 × 107	6.62 × 1019
Cd45S45	2.96 × 108	1.40 × 108	4.36 × 108	4.36 × 1020

Segment-Based Analysis of Sampling Error

The error in the calculated IP using the SSE-MO method originates from the sampling error associated with sampling the integral kernel in the combined Cartesian-MO space. However, not all segments contribute equally to the numerical error. The goal of SSE-MO is to distribute the computational effort in proportion to the contributions from each segment. One insight generated from the SSE-MO calculation is information about the contribution of each segment to the total self-energy operator. We define the cumulative percent contribution for the segments as,where S(i) is the contribution to the self-energy for each segment. The cumulative percent contribution of the segments to the total self-energy operator is denoted as C(M) and is presented for the PbS and CdS quantum dots in Figures and 9, respectively.

Figure 8

Cumulative sum of the percent contributions of the segments composing the sample space for the 2p1h and 1p2h terms of the self-energy (ordinate) versus the segment index (abscissa) is displayed. Parts A, B, and C, are for QDs Pb4S4, Pb44S44, and Pb140S140, respectively.

Figure 9

Cumulative sum of the percent contributions of the segments composing the sample space for the 2p1h and 1p2h terms of the self-energy (ordinate) versus the segment index (abscissa) is displayed. Parts A, B, C are for QDs Cd6S6, Cd24S24, and Cd45S45, respectively.

Analysis of the results revealed that the 2p1h and 1p2h terms show very different behavior. In all cases it was found that only few segments, typically ≤50, had significant contributions to the 1p2h component of the self-energy operator. In contrast, for the 2p1h component, the cumulative sum of the percent contribution increased in a much more gradual manner. The distributions of the standard deviations associated with the segments for the 2p1h and 1p2h terms of the self-energy for the two largest quantum dots, (Pb140S140 and Cd45S45), are presented in Figures and 11, respectively.

Figure 10

Frequency distributions of the standard deviations (in eV) for the segments that compose the sample space of the 2p1h and 1p2h terms of the self-energy computed for Pb140S140, are displayed.

Figure 11

Frequency distributions of the standard deviations (in eV) for the segments that compose the sample space of the 2p1h and 1p2h terms of the self-energy computed for Cd45S45, are displayed.

Analysis of the distributions reveals that the sampling error in the 2p1h term is significantly smaller than the sampling error in the 1p2h term for the two largest quantum dots. These plots also show that the overall sampling error in the calculated IP is dominated by the sampling error in the 1p2h term. The advantage of the SSE-MO method is that, by construction, the SSE-MO scheme is able to extract this information dynamically during the course of the calculation and then allocate more sampling points to segments that have high sampling errors. Because SSE-MO is based on stratified sampling, the conventional stratified sampling error analysis[69] is applicable for the sampling error in the IP calculations. In addition to this segment-based analysis, the overall sampling error in the calculated IPs as a function of the number of sampling points used to construct the self-energy for Pb140S140 quantum dot is presented in Figure .

Figure 12

Sample standard deviation (σ ) in the ionization potential (eV) (ordinate) versus the number of sampling points (abscissa) used to construct the self-energy for the Pb140S140 dot is displayed.

Connection with Diagrammatic Monte Carlo

The SSE-MO method is conceptually similar to diagrammatic MC (diagMC), where terms are evaluated stochastically. However, there are key differences between the two methods. SSE-MO is not diagram-based and the relative importance of the terms are not evaluated using topological connectivity of the vertices. Also, the SSE-MO method uses stratified sampling as opposed to importance sampling, where emphasis is placed on reducing numerical error through variance minimization and numerical effort is predominantly spent on computing the correction term to the self-energy operator. As an intrinsically adaptive approach, the calculation puts more points where they are needed to achieve reduction of numerical error.

Comparison with Laplace-Transformed Approach

The SSE-MO method does not perform Laplace-transformation, but instead relies on stochastic enumeration to reduce the computational cost. Consequently, only 3D and 6D spatial integrals are solved numerically. As a consequence, higher-order derivatives of the self-energy operator (dΣ/dω) can be obtained with relative ease and with very little additional computational cost during the self-energy calculation. This not only allows for calculation of the imaginary component of the self-energy operator, but also open doors for iterative solution of the Dyson equation by Taylor-series expansion of the self-energy operator.Because we are not using a Laplace transformation, the SSE-MO method is well-suited to extending the self-energy calculation to Σ(3) using the P3 correction developed by Ortiz and co-workers.[67] For example, the Laplace transformation of the following term in the P3 expressionwill involve the Monte Carlo numerical integration of a 20-dimensional integral. In the SSE-MO implementation, the dimensionality of the spatial integral will still be six and the MO index will be sampled from the 3p-2h space.Equation can be viewed as the stochastic tensor contraction over the MO indices and can potentially be applied to other branches of quantum mechanics.

Selection of Control-Variate Functions

The use of moment-based fitting ensures that the integral of the Mth-order multinomial comes out to be exact. For this work, a maximum of two Gaussian functions were used and was found to be adequate. For more challenging systems, the number of Gaussian functions can be systematically increased. In addition, metrics other than the moments can be used as criteria for the selection of the Gaussian functions. The choice of the control-variate functions is not restricted to Gaussian functions. For QDs, it is possible to take advantage of the approximate spherical symmetry of the system and construct the control-variate functions from the hydrogenic wave functions with effective hole and particle masses.

Although both density fitting[70] and the control variate schemes use Gaussian functions, their purpose and implementation are very different. When using the control variate scheme, the goal is to reduce numerical error. When using density fitting, the goal is to approximate it. Specifically in the control variate, the integral of f is expressed as

There are two main differences between control variate and density fitting:

When using the control variate scheme, the error in fitting the integral is always calculated. The error in the estimation of the integral comes from the numerical approximation to the analytical fitting error. If we were to replace the numerical integral, ⟨(f–f0)⟩, by an analytical integral, we would recover the exact integral. The origin of error in density fitting comes from the finite expansion of the auxiliary basis. While in the control variate the error is from the numerical approximation to the residue-error integral, ⟨(f–f0)⟩.

Unlike density-fitting’s attribute of fit-once-use-everywhere, the control variate approach is kernel dependent. This means that the integrals ⟨fK⟩ and ⟨fK⟩ will have different control variate parameters α and α, respectively. These parameters are obtained by minimizing the variance as shown below:One approach to do the above integrals efficiently is to first expand the square term and then perform the α-independent integrals separately as shown below.

Conclusions

This work presents the development and implementation of the stratified stochastic enumeration of molecular orbitals (SSE-MO) method for construction of the self-energy operator. The central idea of this method is to express the self-energy operator in a composite space, which is generated by combining the 3D Cartesian space of molecular orbitals with the discrete integer space of the molecular orbital indices. In conjunction, a stratified sampling Monte Carlo scheme was also developed for the efficient evaluation of the complex self-energy operator and its frequency derivatives. The SSE-MO method was applied to a series of CdS and PbS QDs, and the IPs of these QDs were obtained from both single-shot and iterative solution of the second order diagonal approximation to the Dyson equation. The results from these calculations showed that the IPs decreased with increasing dot size. The imaginary component of the self-energy operator was used to construct the single-pole frequency-dependent spectral functions of the quantum dots. The quantum dots with the longest relative lifetimes of the quasi-hole state were identified. The strategy of stochastic enumeration used in the SSE-MO method can also be interpreted in the broader context of stochastic tensor contraction methods and can be applied to other areas of quantum mechanics, where the sequential enumeration of summations is computationally prohibitive.