Second-Order Perturbation Theory in Continuum Quantum Monte Carlo Calculations

We report on the first results for the second-order perturbation theory correction to the ground-state energy of a nuclear many-body system in a continuum quantum Monte Carlo calculation. Second-order (and higher) perturbative corrections are notoriously difficult to compute in most \textit{ab~initio} many-body methods, where the focus is usually on obtaining the ground-state energy. By mapping our calculation of the second-order energy correction to an evolution in imaginary time using the diffusion Monte Carlo method, we are able to calculate these corrections for the first time. After benchmarking our method in the few-body sector, we explore the effect of charge-independence breaking terms in the nuclear Hamiltonian. We then employ the new approach to investigate the many-body, perturbative, order-by-order convergence that is fundamental in modern theories of the nucleon-nucleon interaction derived from chiral effective field theory. Our approach is quite general and promises to be of wide applicability.

Nuclear many-body theory has been experiencing a renaissance over the last two decades, with several techniques now being routinely used to carry out firstprinciples studies of many strongly interacting nucleons [1].These include exact diagonalization approaches (when the basis size is manageable) [2], quantum Monte Carlo QMC methods (stochastic evolution, whether in the continuum [3] or on a lattice [4]), more or less controlled truncations such as many-body perturbation theory (limiting oneself to the first several orders), [5] selfconsistent Green's function approaches (implementing a specific class of diagrams up to infinite order), [6] Coupled Cluster (generating npnh excitations off a reference state) [7] and its younger cousin the in-medium Similarity Renormalization group [8].Several of these approaches are also very important in broader applications of the quantum many-body problem, including, for example, condensed matter physics and the study of ultracold atoms [9][10][11][12][13].The difference between nuclear and other many-body techniques is driven by the two-(or higher-) body interaction at play.Using nucleons as degrees of freedom, the interaction between particles has historically been written down as the combination of short-range/contact terms and terms corresponding to the exchange of mesons (most notably pions).The study of nuclear interactions itself has been rejuvenated by the aforementioned renaissance in nuclear many-body techniques (see, e.g., Refs.); another driving force has been the recasting of earlier phenomenological forces in the language of effective field theory [35].The contact terms arising in the context of pionless effective field theory (EFT), also of applicability to cold atoms, where the range of the interaction is much shorter than the interparticle spacing, are then supplemented by the exchange of the lightest mesons, most notably in the context of chiral EFT.Several open questions remain at the forefront of the study of nuclear forces, such as the few-body observables that should be used as constraints, the importance of higher orders and uncertainty estimates, as well as the determination of an appropriate (renormalization-group invariant) powercounting scheme [36].Once again, much of the recent progress (and hope for further developments) has been driven by the interplay of novel nuclear many-body techniques and designer nuclear forces.
As noted, an important open question, on which much work has been dedicated recently, relates to the relative importance of power counting and nonperturbative renormalization in chiral EFT.Given the significance of the interface between nuclear forces and many-body techniques, a newcomer might be forgiven for assuming that the problem of perturbative vs nonperturbative force vs many-body technique will have received considerable attention.However, until very recently this has not been the case; to see why, let us recall the quantummechanical expressions for 1st and 2nd order corrections to the ground-state energy for a given perturbation V : where ψ 0 , ψ k , E 0 , and E k are the ground and kth excited eigenstates and energies of the unperturbed Hamiltonian, respectively; crucially, all of these are many-body entities, i.e., in contradistinction to the textbook case, they cannot be tackled analytically.The problem becomes immediately apparent: While the first-order correction E (1) 0 is easily computable in a many-body context (e.g., in the diffusion Monte Carlo (DMC) method the propagator would involve the unperturbed Hamiltonian and the observables would involve V ), computing the second-order correction E (2) 0 is dramatically more complicated, as it requires knowledge of the complete energy spectrum E k .While progress has been made on tackling individual excited states in an ab initio setting, e.g., by imposing the arXiv:2302.07285v1[nucl-th] 14 Feb 2023 appropriate symmetries on a variational wave function [37], that is still a manual process (separately designing a new wave function for each new state) that cannot be easily generalized beyond a handful of cases.This significant task is woefully underdiscussed in the many-body literature; even the few exceptions cast the problem in quite different language (cumulant autocorrelation functions in [38] and auxiliary fields in [39]).In the present Letter, we show how to map the computation of the second-order energy correction E (2) 0 to an evolution in imaginary time; we do this (for the first time) in the context of a continuum QMC technique oriented toward nuclear forces, but the recasting we employ is much more general and could plausibly be implemented in other areas of physics or distinct many-body techniques.Both for reasons of clarity and in order to establish our approach's credentials, we start from a trivial (two-body) application as a benchmarking exercise; we then turn to a more challenging application, resulting from the introduction of a charge-independence breaking term in the nuclear interaction.Significantly, we then turn to the main application, namely using the new approach as a detailed probe of the perturbativeness (or lack thereof) of chiral EFT interactions at the many-body level.Crucially, the nonperturbative calculation does not correspond to an ad hoc interaction employed merely to make higherorder terms small [39,40] but to an actual chiral EFT interaction in widespread use (the local next-to-next-to leading order (N 2 LO) interaction of Gezerlis et al [18]).In other words, the new approach will allow us to study how perturbative a given interaction is at the many-body level, going beyond first-order perturbation theory.The fact that we address three distinct applications reflects the generality of the proposed methodology.
We compute Eq. ( 1) using the DMC method.We first start with many sets of particle positions (walkers) that are distributed according to a trial wave function ψ T by a preliminary variational Monte Carlo (VMC) calculation [41].The DMC method then projects out the lowest energy eigenstate ψ 0 by treating the Schrödinger equation as a diffusion equation in imaginary time τ and propagating the trial wave function up to large τ .This can be demonstrated by expanding the trial wave function in terms of the complete set of exact eigenstates, i , and then applying the imaginary time propagation operator, Since DMC projects out the ground-state wave function of the system of interest, we start our calculation of the second-order energy correction by assuming we have access to the ground state, and consider the quantity, By an insertion of the identity, and then a splitting of the sum between k = 0 and k = 0, Eq. ( 3) can be recast as In the limit of long imaginary time the exponential in the second term can be neglected and by comparing with Eq. ( 1), we see we have recovered both the first-and second-order corrections to the ground-state energy, Therefore, for a given perturbation V , the second-order correction to the ground-state energy can be computed in a many-body context by implementing Eq. ( 3).The question then, is how to evaluate Eq. ( 3) in a QMC context, where we are starting with a realistic trial wave function, and not the ground state.The plan will be to compute the integrand of Eq. ( 3), fit its behavior to the form of Eq. ( 4), and finally employ an extrapolated estimate.In a DMC calculation the full imaginary time τ is broken into many successive small-τ segments, which means our quantity of interest can be written as or as a Monte Carlo integral, where ψ T is our trial wave function and G(R , R; ∆τ ) is the importance-sampled short-imaginary-time propagator [42].We can then finally rewrite Eq. ( 7) as a twice importance sampled Monte Carlo integral, where the R i are now positions distributed according to G and the R i,0 are positions distributed according to the trial wave function.An example calculation of Eq. ( 3) in a DMC context by way of Eq. ( 8)is shown in Fig. 1 (the physical details will be introduced in due course).By comparing with Eq. ( 8) we can trace the behavior of the integral I.We should expect that at very small τ the integral should be equal to V 2 , since the exponential in Eq. ( 3) has only just started to decay.In addition, we expect the integral to decay to V 2 in the limit of large-τ based on Eq. ( 4).
Therefore, to compute the second-order energy correction using Eq.(3), we compute Eq. ( 8) at every step of the DMC process.To find an accurate estimate for the second-order correction, we perform two DMC runs.8) for perturbation from a neutron-neutron interaction to a neutron-proton interaction.By comparing with Eq.( 8) we see the expected behavior.At very small τ we essentially recover V 2 and in the limit of long τ we have the expected result that Istep decays down to V 2 .(See text for more details.) First we perform an initial run starting from walkers distributed according to the trial wavefunction, and then a second run starting from walkers distributed according to the projected ground-state wave function.We then compute an extrapolated estimate between this groundstate initialized run and the initial trial wavefunction run, which is necessary for the expectation value of any operator that does not commute with the Hamiltonian when forward walking techniques are not employed [43].
To benchmark this new method, the second-order energy corrections were calculated and compared against nonperturbative results for a variety of few-body systems in a harmonic oscillator potential with oscillator frequency ω.To perform these calculations we used a harmonic oscillator eigenstate basis trial wave function, with variational parameters differing slightly from the ground-state wave function of trapped particles, in the spirit of Ref. [44].
Our initial few body tests involve starting from noninteracting particles in the harmonic trap, and perturbatively including an interaction in the form of a simple Gaussian interaction, V = ae −q 2 (r2−r1) 2 , where a = 1.0 ω and q = 1.0 fm −1 .We first carry this out for simple two-particle systems in a variety of harmonic trap strengths.As was to be expected, we see that the quality of the answer provided by perturbation theory (whether 1st or 2nd order) is reasonable, as long as the perturbation is small.(See the first three rows of Table I.) Though we produce the answer using the DMC machinery, for few-body systems the results can also be produced quasi-analytically.However, because our interest lies in applying our method to nuclear systems, we now introduce our general E [2]  Non-Pt.E [n] refers to the sum of energy corrections up to the nth order.All energies are in units of ω.Upper: Results for a Gaussian perturbation between two noninteracting particles, Lower: Results for a perturbation from the neutron-neutron interaction to the neutron-proton interaction.
unperturbed Hamiltonian where N is the total number of particles and V (r ij ) is the interaction between spin up (primed) and spin down (unprimed) particles.In the case of the nucleon-nucleon interaction, this V (r ij ) contains contributions like the tensor force, spin-orbit, etc.Here, with a view to abstracting away the details of the interaction (allowing us to employ sophisticated wave functions in DMC) we limit ourselves to low density/s-wave interactions.Historically, nuclear many-body calculations have used phenomenological potentials that are fit to experimental data [3,45].For this initial nuclear investigation we employ a phenomenological potential of the Pöschl-Teller type as in Ref. [43], where v 0 and µ are parameters that can be tuned in order to reproduce the scattering length and effective range of a nuclear interaction, which describes the interaction completely at low energies [46].Having benchmarked the method at the two-body level, we now consider a range of particle numbers up to N = 6, and set our unperturbed interaction to correspond to a neutron-neutron interaction (still in a harmonic trap).Starting from this unperturbed neutron system, we apply a perturbation such that we move the system from a two-neutron interaction to a neutron-proton one as illustrated in the upper panel of Fig. 2. The nonperturbative ground-state energies of particles interacting through the nn or np potentials are compared against our perturbative calculations in the bottom three Upper panel: Perturbing potential given by the difference between the s-wave neutron-neutron interaction and the neutron-proton interaction.Lower panel: Perturbation given by the difference between the 1 S0 N 2 LO and NLO chiral EFT interactions for a range of coordinate space cutoff values.
rows of Table I.We can see that while the first-order correction does a reasonable job of approximating the neutron-proton case, the second-order correction is required in order to accurately reproduce the nonperturbative neutron-proton ground-state energy to within one standard deviation.
Having applied our method to the few-body sector where the perturbation is not too difficult to handle, we now turn to a realistic nuclear system.For this study we modify our Hamiltonian from Eq. ( 9) removing the external one-body oscillator potential.We also use N = 66 particles as this gives a very good approximation to the thermodynamic limit [47].In addition, since both the VMC and the DMC methods depend strongly on the choice of trial wave function, and we are interested in the physics of pure infinite neutron matter, where pairing is known to be important [43,48], we make use of periodic boundary conditions and the Jastrow-BCS wave function [49], As in Eq. ( 9), the primed indices correspond to spin up neutrons, and the unprimed indices to spin down neutrons.The pairing functions φ(r) = β(r) + n α n e ikn•r capture both long-and short-range physics as described in Ref. [49].
In recent decades, the ab initio many-body community has largely moved on from phenomenological potentials like the ones used for our preliminary tests, and modern interactions have been built using a chiral EFT framework.In addition to reproducing known two-body experimental results, these chiral EFT interactions also capture the important symmetries from the underlying theory of quantum chromodynamics [50].Chiral EFT also provides a systematic expansion for the nuclear force, with terms included at successive orders based on their importance following a power-counting scheme [51], where the superscript keeps track of the power of the expansion parameter Q/Λ b , which depends on the momenta of the nucleons, or the pion mass (Q), and the scale at which the chiral EFT expansion breaks down (Λ b ).This expansion structure informed the choice of interactions used for this work.To study the physics of neutron matter, we use the chiral EFT interactions of Ref. [15] that have been tuned to reproduce the dominant 1 S 0 channel interaction between opposite spin neutrons.We start from an unperturbed system with an interaction containing terms up to NLO and perturb to include terms at N 2 LO, for a range of coordinate space cutoffs, as seen in the bottom panel of Fig. 2.
The result of carrying out such a many-body computation up to second-order in perturbation theory is shown in Fig. 3 for a number of densities.It is worth emphasizing that when the core of the potential is soft (R 0 ≥ 1.1 fm) the second-order perturbation theory results for treating the difference between N 2 LO and NLO as a perturbation match the nonperturbative N 2 LO results to within 1% error.
This follows the general trend found via an exact diagonalization for the deuteron in Ref. [52].In each case we see that the first-order correction is positive, and the second-order correction is negative.It can also be seen that at the higher densities the agreement between the N 2 LO results and the second-order perturbation theory results is worsened.This is likely due to the fact that at large densities other interaction channels become important and so we are no longer in the regime where the neutron-neutron interaction is correctly described by a pure 1 S 0 interaction.
The results for R 0 = 1.0 fm in Fig. 3 suggest that the strong repulsive core coincides with the need to use at least third-order corrections, and raises the question as to whether the perturbative-based chiral EFT interactions are appropriately behaved [53,54].We have also attempted to perturb from NLO to LO, with a perturbation V LO − V NLO , (and vice versa) with limited success.This provides further evidence that the difference between LO and NLO is not perturbative.
In summary, we have developed a method for calculating the second-order perturbation theory correction to the ground-state energy in a continuum quantum Monte Carlo context.We have benchmarked this method against nonperturbative results in the few-body sector, and performed initial calculations for few-body neutron systems in a trap that could pave the way for more complicated four-species calculations in the future.Finally, we have completed the first calculations for pure infinite neutron matter that include the second-order correction to the ground-state energy.In addition, we have also tested the behavior of popular chiral EFT interactions and found indications of nonperturbativeness.In the future, this method has broad applicability and could become a staple in the toolbox of anyone making use of DMC methods or its extensions (e.g. the auxiliary-field diffusion Monte Carlo method or the Green's Function Monte Carlo method. The work of R.C. and A.G. was supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada and the Canada Foundation for Innovation (CFI).The work of J.E.L. was supported by the ERC Grant No. 307986 STRONGINT and the BMBF under Contract No. 05P15RDFN1.Computational resources were provided by SHARCNET and NERSC.

2 FIG. 1 .
FIG.1.A DMC calculation of Eq.(8) for perturbation from a neutron-neutron interaction to a neutron-proton interaction.By comparing with Eq.(8) we see the expected behavior.At very small τ we essentially recover V 2 and in the limit of long τ we have the expected result that Istep decays down to V 2 .(See text for more details.)

FIG. 3 .
FIG.3.Comparison between the nonperturbative results for NLO, N 2 LO and the perturbative results for perturbing from NLO to N 2 LO at both first and second order.