Quantum Monte Carlo calculations of light nuclei with local chiral two- and three-nucleon interactions

TL;DR

This paper evaluates local chiral effective field theory interactions in quantum Monte Carlo calculations, analyzing their impact on light nuclei properties, including binding energies and distribution functions, with focus on nucleon interactions and methodological details.

Contribution

It provides a detailed analysis of local chiral nucleon-nucleon and three-nucleon interactions and their effects on light nuclei properties using quantum Monte Carlo methods.

Findings

Local chiral interactions affect deuteron properties.

Operator ambiguities influence three-nucleon interactions.

Results on binding energies and distributions for light nuclei.

Abstract

Local chiral effective field theory interactions have recently been developed and used in the context of quantum Monte Carlo few- and many-body methods for nuclear physics. In this work, we go over detailed features of local chiral nucleon-nucleon interactions and examine their effect on properties of the deuteron, paying special attention to the perturbativeness of the expansion. We then turn to three-nucleon interactions, focusing on operator ambiguities and their interplay with regulator effects. We then discuss the nuclear Green's function Monte Carlo method, going over both wave-function correlations and approximations for the two- and three-body propagators. Following this, we present a range of results on light nuclei: Binding energies and distribution functions are contrasted and compared, starting from several different microscopic interactions.

Tables

Table 1. Table 1: Deuteron properties including the binding energy E b subscript 𝐸 𝑏 E_{b} , asymptotic D / S 𝐷 𝑆 D/S ratio η d subscript 𝜂 𝑑 \eta_{d} , quadrupole moment Q d subscript 𝑄 𝑑 Q_{d} (impulse approximation), and root-mean-square (rms) matter radius ⟨ r d 2 ⟩ expectation-value subscript superscript 𝑟 2 𝑑 \sqrt{\expectationvalue{r^{2}_{d}}} . Electromagnetic interaction effects are neglected here (when included they change the values below only within the uncertainties). The uncertainties for the local chiral interactions represent the discussed truncation error estimate. See text for more details. Experimental values are from Refs.  Leun and Alderliesten ( 1982 ); Rodning and Knutson ( 1990 ); Bishop and Cheung ( 1979 ); Simon  et al. ( 1981 ) .

	$R_0=1.0$ fm	$R_0=1.2$ fm	Exp
$E_b$ (MeV)	$2.21\left(2\right)$	$2.20\left(3\right)$	$2.224575\left(9\right)$
${\eta }_d$	$0.0263\left(3\right)$	$0.0267\left(6\right)$	$0.0256\left(4\right)$
$Q_d\left({\text{fm}}^2\right)$	$0.286\left(5\right)$	$0.289\left(6\right)$	$0.2859\left(3\right)$
$\sqrt{⟨r_d^2⟩}$ (fm)	$1.97\left(2\right)$	$1.97\left(3\right)$	$1.9660\left(68\right)$

Table 2. Table 2: Kinetic and potential energy contributions to the GFMC energy at LO , NLO , and N 2 LO for both cutoffs and for a particular choice of 3 ​ N 3 𝑁 3N D 𝐷 D and E 𝐸 E operators [Eqs. ( 65b ) and ( 66a )]. For R 0 = 1.0 subscript 𝑅 0 1.0 R_{0}=1.0  fm, c D = 0.0 subscript 𝑐 𝐷 0.0 c_{D}=0.0 , c E = − 0.63 subscript 𝑐 𝐸 0.63 c_{E}=-0.63 , while for R 0 = 1.2 subscript 𝑅 0 1.2 R_{0}=1.2  fm, c D = 3.5 subscript 𝑐 𝐷 3.5 c_{D}=3.5 , c E = 0.085 subscript 𝑐 𝐸 0.085 c_{E}=0.085   Lynn  et al. ( 2016 ) . For comparison, we also show results for the Argonne v 18 subscript 𝑣 18 v_{18} N ​ N 𝑁 𝑁 NN interaction supplemented with the UIX 3 ​ N 3 𝑁 3N interaction. V 3 ​ N subscript 𝑉 3 𝑁 V_{3N{}} stands for the sum of all 3 ​ N 3 𝑁 3N contributions. All energies are in MeV.

		$R_0$ (fm)	$K$	$V_{NN}$	$V_{3N}$	$V_{C,c_1}$	$V_{C,c_3}$	$V_{C,c_4}$	$V_{D2}$	$V_{E\tau }$
\isotope[3]H	LO	1.0	$60.2\left(2\right)$	$-74.0\left(2\right)$
		1.2	$55.5\left(1\right)$	$-67.8\left(1\right)$
	NLO	1.0	$46.3\left(2\right)$	$-54.4\left(2\right)$
		1.2	$36.9\left(2\right)$	$-45.0\left(2\right)$
	N2LO	1.0	$42.7\left(2\right)$	$-50.6\left(2\right)$	$-1.32\left(2\right)$	$-0.08\left(1\right)$	$-1.22\left(2\right)$	$-0.53\left(7\right)$	$0.0$	$0.51\left(1\right)$
		1.2	$37.6\left(1\right)$	$-45.9\left(1\right)$	$-0.87\left(1\right)$	$-0.06\left(1\right)$	$-0.27\left(1\right)$	$-0.35\left(3\right)$	$-0.09\left(1\right)$	$-0.10\left(1\right)$
	AV18+UIX		$51.4\left(2\right)$	$-59.4\left(2\right)$	$-1.23\left(1\right)$
\isotope[3]He	LO	1.0	$60.0\left(1\right)$	$-73.0\left(1\right)$
		1.2	$55.0\left(1\right)$	$-67.4\left(1\right)$
	NLO	1.0	$43.9\left(3\right)$	$-51.5\left(3\right)$
		1.2	$36.4\left(2\right)$	$-44.3\left(2\right)$
	N2LO	1.0	$41.3\left(3\right)$	$-50.5\left(2\right)$	$-1.27\left(2\right)$	$-0.08\left(1\right)$	$-1.16\left(2\right)$	$-0.53\left(9\right)$	$0.0$	$0.49\left(1\right)$
		1.2	$36.8\left(1\right)$	$-45.1\left(1\right)$	$-0.83\left(1\right)$	$-0.05\left(1\right)$	$-0.26\left(1\right)$	$-0.34\left(3\right)$	$-0.08\left(1\right)$	$-0.09\left(1\right)$
	AV18+UIX		$50.4\left(1\right)$	$-58.4\left(1\right)$	$-1.20\left(1\right)$
\isotope[4]He	LO	1.0	$142.0\left(2\right)$	$-193.4\left(2\right)$
		1.2	$132.1\left(2\right)$	$-183.5\left(2\right)$
	NLO	1.0	$90.2\left(3\right)$	$-115.9\left(3\right)$
		1.2	$73.0\left(3\right)$	$-99.4\left(2\right)$
	N2LO	1.0	$90.9\left(2\right)$	$-116.1\left(2\right)$	$-7.46\left(4\right)$	$-0.41\left(1\right)$	$-6.74\left(5\right)$	$-2.6\left(2\right)$	$0.0$	$2.34\left(2\right)$
		1.2	$79.9\left(2\right)$	$-106.3\left(2\right)$	$-5.56\left(4\right)$	$-0.30\left(1\right)$	$-1.78\left(3\right)$	$-1.7\left(2\right)$	$-1.24\left(4\right)$	$-0.51\left(1\right)$
	AV18+UIX		$115.8\left(1\right)$	$-140.4\left(1\right)$	$-6.73\left(2\right)$

Table 3. Table 3: Point-proton radii as calculated in Eq. ( 59 ) at LO , NLO , and N 2 LO for both cutoffs for the A = 3 , 4 𝐴 3 4 A=3,4 nuclei. The theoretical uncertainties are from both the GFMC statistical uncertainties as well as the theoretical uncertainty coming from the truncation of the chiral expansion as described in Sec.  II . Experimental values are from Refs.  Patrignani  et al. ( 2016 ); Friar  et al. ( 1997 ); Sick ( 2014 ); Amroun  et al. ( 1994 ) with uncertainties calculated using standard propagation of uncertainty methods. All radii are in fm.

	\isotope[3]H		\isotope[3]He		\isotope[4]He
$R_0$	$1.0\text{fm}$	$1.2\text{fm}$	$1.0\text{fm}$	$1.2\text{fm}$	$1.0\text{fm}$	$1.2\text{fm}$
LO	$1.27\left(35\right)$	$1.27\left(37\right)$	$1.36\left(56\right)$	$1.36\left(52\right)$	$1.02\left(55\right)$	$1.00\left(53\right)$
NLO	$1.62\left(10\right)$	$1.64\left(13\right)$	$1.92\left(16\right)$	$1.88\left(18\right)$	$1.57\left(15\right)$	$1.53\left(18\right)$
N2LO	$1.55\left(4\right)$	$1.55\left(6\right)$	$1.77\left(5\right)$	$1.77\left(6\right)$	$1.43\left(5\right)$	$1.42\left(7\right)$
Exp	1.59(10)		1.78(2)		1.46(1)

Equations

$$\underset{\tau \to \infty }{\lim }{\mathrm{e}}^{-H\tau }|{\Psi }_T⟩\to |{\Psi }_0⟩,$$

$$V_{NN}=V_{NN}^{(0)}+V_{NN}^{(2)}+V_{NN}^{(3)}+\dots ,$$

$$\mathbf{p}$$

$$\mathbf{q}$$

$$\mathbf{k}$$

$$f(p^2)=\exp \left[-(p^2/{\Lambda }_{NN}^2{)}^n\right],$$

$$V(\mathbf{p},{\mathbf{p}}^{\prime })\to V(\mathbf{p},{\mathbf{p}}^{\prime })f(p^2)f(p^{\prime 2}).$$

$${\delta }^{(3)}(\mathbf{r})\to {\delta }_{R_0}(r)=\frac{{\mathrm{e}}^{-(r/R_0{)}^n}}{\frac{4\pi }{n}\Gamma (\frac{3}{n})R_0^3}.$$

$$\int \text{differential}[3]r{\delta }_{R_0}(r)=1.$$

$$f_{\text{long}}(r)=1-{\mathrm{e}}^{-(r/R_0{)}^4}.$$

$$\begin{array}{rl}\Delta X^{{\text{N}}^2\text{LO}}=\max & (Q^4\left|X^{\text{LO}}\right|,Q^2\left|X^{\text{NLO}}-X^{\text{LO}}\right|,\\ & Q\left|X^{{\text{N}}^2\text{LO}}-X^{\text{NLO}}\right|),\end{array}$$

$${\psi }_d^{(M_J)}(\mathbf{r})=\left[\frac{u(r)}{r}+\frac{S_{12}(\hat{\mathbf{r}})}{\sqrt{8}}\frac{w(r)}{r}\right]\frac{{\chi }_{M_J}}{\sqrt{4\pi }},$$

$${\int }_0^{\infty }\text{differential}r{r}^2\left[{\left(\frac{u(r)}{r}\right)}^2+{\left(\frac{w(r)}{r}\right)}^2\right]=1.$$

$$\frac{\tilde{u}(q)}{q}$$

$$\frac{\tilde{w}(q)}{q}$$

$${\int }_0^{\infty }\frac{\text{differential}q{q}^2}{(2\pi {)}^3}\left[{\left(\frac{\tilde{u}(q)}{q}\right)}^2+{\left(\frac{\tilde{w}(q)}{q}\right)}^2\right]=1.$$

$$n(q)=\frac{1}{4\pi }\left[{\left(\frac{\tilde{u}(q)}{q}\right)}^2+{\left(\frac{\tilde{w}(q)}{q}\right)}^2\right],$$

$$\int \frac{\text{differential}[3]q}{(2\pi {)}^3}n(q)=1.$$

$$F_{C,M_J}(q)=\frac{1}{2}\int \text{differential}[3]r^{\prime }{\rho }_d^{(M_J)}({\mathbf{r}}^{\prime }){\mathrm{e}}^{\mathrm{i}\mathbf{q}\text{dotproduct}{\mathbf{r}}^{\prime }},$$

$${\rho }_d^{(0)}({\mathbf{r}}^{\prime })=$$

$${\rho }_d^{(\pm 1)}({\mathbf{r}}^{\prime })=$$

$$C_0(r)$$

$$C_2(r)$$

$$T_{20}(q)\approx -\sqrt{2}\frac{F_{C,0}^2(q)-F_{C,1}^2(q)}{F_{C,0}^2(q)+2F_{C,1}^2(q)}.$$

$$V_{\text{pert}}\equiv V_{{\text{N}}^2\text{LO}}-V_{\text{NLO}}$$

$$\begin{gathered} \text{expectationvalue}{H}_{\text{NLO}}+V_{\text{pert}}{\psi }_d^{(\text{NLO})} \\ =E_{\text{NLO}}+\text{expectationvalue}{V}_{\text{pert}}{\psi }_d^{(\text{NLO})}. \end{gathered}$$

$$V_C$$

$$V_D$$

$$V_E$$

$${\mathcal{X}}_{ij}(\mathbf{r})\equiv X_{ij}(\mathbf{r})-\frac{4\pi }{m_{\pi }^2}{\delta }_{R_{3N}}(r){\sigma }_i\text{dotproduct}{\sigma }_j.$$

$$V_{D1}$$

$$V_{D2}$$

$$\begin{array}{rl}\{\text{mathbbm}1,& {\sigma }_i\text{dotproduct}{\sigma }_j,{\tau }_i\text{dotproduct}{\tau }_j,{\sigma }_i\text{dotproduct}{\sigma }_j{\tau }_i\text{dotproduct}{\tau }_j,\\ & {\sigma }_i\text{dotproduct}{\sigma }_j{\tau }_i\text{dotproduct}{\tau }_k,[({\sigma }_i\text{crossproduct}{\sigma }_j)\text{dotproduct}{\sigma }_k][({\tau }_i\text{crossproduct}{\tau }_j)\text{dotproduct}{\tau }_k]\}.\end{array}$$

$$V_{E\tau }$$

Full text

I Introduction

Theoretical nuclear physics has undergone a renaissance in recent decades because of two main developments: The increasing reach and precision of nuclear many-body methods, and the formulation of systematic nuclear interactions based on chiral effective field theory (EFT).

Ab initio many-body methods in nuclear physics include the no-core shell model Barrett et al. (2013), nuclear lattice simulations Epelbaum et al. (2011), the coupled-cluster method Hagen et al. (2014, 2016), the in-medium similarity renormalization group (SRG) method Hergert et al. (2016), self-consistent Green’s function methods Carbone et al. (2013); Somà et al. (2014), and quantum Monte Carlo (QMC) methods Carlson et al. (2015). Among these, QMC methods, which are based on the imaginary-time evolution of a trial wave function and include the Green’s function Monte Carlo (GFMC) method and the auxiliary-field diffusion Monte Carlo (AFDMC) method, are notable for their high accuracy across various physical systems.

In a typical calculation, QMC methods reach uncertainties of $\sim 1\%$. By design, QMC methods introduce only a limited number of approximations that can be controlled and accounted for systematically. Both the GFMC method and the AFDMC method rely on the diffusion equation

[TABLE]

where $H$ is the Hamiltonian of the system, $\tau$ is imaginary time, and $\ket{\Psi_{T}}$ is a trial state for the system not orthogonal to the ground state $\ket{\Psi_{0}}$. These “diffusion” methods solve Eq. (1) stochastically by casting it as a path integral and sampling the paths using Monte Carlo methods. This allows one to extract ground- and low-lying excited-state properties of nuclear systems with high accuracy.

Furthermore, QMC methods are notable because they approach the many-body problem with a correlated wave-function–oriented framework. For certain nuclear systems, e.g., the Hoyle state of \isotope[12]C, many-body methods that rely on basis-set expansions can experience difficulties in capturing physics that requires a large number of basis states to describe, such as clustering effects. For QMC methods, which rely on a trial wave function to describe the state of interest, these effects are more straightforward to incorporate. While the GFMC method has an unfavorable scaling behavior with respect to the nucleon number $A$, the above-mentioned strengths make QMC calculations of smaller systems an ideal benchmark for other methods.

Besides the exciting advancements for nuclear many-body methods, the development of chiral EFT as a tool for the derivation of systematic nuclear interactions connected to the underlying theory of strong interactions, quantum chromodynamics (QCD), represents a major step forward in nuclear theory. The idea, first presented by Weinberg in the 1990s Weinberg (1990, 1991, 1992), is to write down the most general Lagrangian consistent with all the symmetries of the underlying theory, including the chiral symmetry of low-energy QCD, in terms of the relevant degrees of freedom at low energies, i.e., nucleons and pions. Together with a power counting scheme to order the resulting contributions according to their importance, the result is a low-energy effective field theory for nuclear forces. The idea was further developed by van Kolck et al. in early pioneering work Ordóñez et al. (1994); van Kolck (1994); Ordóñez et al. (1996). The first “modern” chiral EFT interactions with a $\chi^{2}/\text{datum}$ around 1 in a fit to $NN$ scattering data were introduced in the early 2000s by Entem and Machleidt Entem and Machleidt (2003) and by Epelbaum et al. Epelbaum et al. (2005).

The advantages of the chiral EFT approach to nuclear interactions over commonly used phenomenological approaches include the ability to systematically determine consistent many-body interactions and electroweak currents, as well as to estimate theoretical uncertainties. The chiral EFT approach, however, is not without some open problems. These include, e.g., power counting schemes, residual cutoff dependences, and associated regulator artifacts. In the past few years, various groups have investigated several aspects involved in constructing nuclear forces from chiral EFT, e.g., the fitting protocol Ekström et al. (2013, 2015), regulators Tews et al. (2016); Dyhdalo et al. (2016), or uncertainty estimates Epelbaum et al. (2015a), with the goal of improving predictions based on chiral interactions.

For many years, chiral EFT interactions could not be implemented in QMC methods because these interactions are derived in momentum space and are typically nonlocal while QMC methods rely on local interactions. In spite of some work to remedy this shortcoming Lynn and Schmidt (2012), it remains technically challenging to develop QMC methods that both can use nonlocal interactions and lead to results without large statistical uncertainties; see also Ref. Roggero et al. (2014) for an alternative approach.

In recent years, however, it was realized that all sources of nonlocality can be removed up to next-to-next-to-leading order (N2LO) in the standard Weinberg power counting. This led to the development of local chiral interactions and their implementation in QMC methods Gezerlis et al. (2013, 2014); Tews et al. (2016); Piarulli et al. (2015, 2016) and has allowed for the first QMC studies of light nuclei, neutron matter, and other light neutron systems with chiral EFT interactions at N2LO including $3N$ interactions Lynn et al. (2014, 2016); Klos et al. (2016); Chen et al. ; Gandolfi et al. (2017). In this paper, we provide details for the calculations of light nuclei and present additional results.

The structure of this paper is as follows. In Sec. II, we discuss how local chiral EFT interactions have been derived, highlight some interesting features of these local interactions, and discuss open questions. In Sec. III, we describe the GFMC and AFDMC methods in more detail and discuss the necessary changes in order to accommodate local chiral EFT interactions. In Sec. IV, we provide a summary of results for light nuclei obtained with QMC methods and chiral EFT interactions. Finally, we give a summary in Sec. V.

II Local Chiral Interactions

As stated in the introduction, chiral EFT is a systematic way of organizing nuclear interactions. Based on the most general Lagrangian consistent with the symmetries of QCD, and combined with a power counting scheme, it is possible to expand nuclear interactions in a series with the expansion parameter $p/\Lambda_{b}$, where $p$ is a typical low-momentum scale in nuclear systems of the order of the pion mass $m_{\pi}$, and $\Lambda_{b}\sim 500$ MeV is the breakdown scale that determines the range of applicability of the EFT. Then, nuclear interactions can be arranged as

[TABLE]

where the superscript denotes the chiral order (the power of $Q\sim p/\Lambda_{b}$ in the corresponding contributions). At leading order (LO), $Q^{0}$, two contributions add to the nuclear interaction: the one-pion exchange (OPE) and momentum-independent short-range contact interactions. At higher orders, two-pion–exchange interactions (TPE) and momentum-dependent (derivative) contact interactions appear. For more details on chiral EFT, see Refs. Epelbaum et al. (2009); Machleidt and Entem (2011).

Because chiral EFT is naturally formulated in momentum space, it can contain nonlocal parts by construction. In this section, we review the strategy to remove all sources of nonlocality, present selected results for the deuteron, show details of the inclusion of $3N$ interactions at N2LO, and discuss several open questions regarding locality and regularization.

II.1 Locality in chiral EFT

Chiral EFT interactions, with the exception of early pioneering work Ordóñez et al. (1994), have been developed in momentum space. We define the incoming (outgoing) single-particle momenta in the $NN$ sector as $\mathbf{p}_{1}$, $\mathbf{p}_{2}$ ($\mathbf{p}_{1}^{\prime}$, $\mathbf{p}_{2}^{\prime}$). Then the incoming (outgoing) relative momentum $\mathbf{p}$ ($\mathbf{p}^{\prime}$), the momentum transfer $\mathbf{q}$, and momentum transfer in the exchange channel $\mathbf{k}$ are defined as

[TABLE]

The Fourier transformation of a function of $\mathbf{q}$ leads to a local function in coordinate space that depends only on the two-particle distance $\mathbf{r}$, whereas a function of $\mathbf{k}$ does not.

Chiral EFT $NN$ interactions depend on two linearly independent momenta out of the four possible momenta ($\mathbf{p}$ and $\mathbf{p}^{\prime}$ or $\mathbf{q}$ and $\mathbf{k}$). There are two possible sources of nonlocality ($\mathbf{k}$ dependence):

The momentum-space regulator functions used to regulate high-momentum contributions to the interaction and 2. 2.

momentum-dependent higher order contact operators.

We review the method to remove these sources of nonlocality, which was first discussed in Ref. Freunek and later employed in practice in Refs. Gezerlis et al. (2013, 2014).

II.1.1 Local regulators

When employing chiral EFT interactions in few- and many-body calculations, momentum-dependent regulator functions need to be introduced to cutoff divergences from high-momentum modes. The typical functional form employed to regulate both the short-range contact interactions and long-range pion exchanges in nonlocally regulated chiral EFT interactions is

[TABLE]

where $\Lambda_{NN{}}$ is the momentum-space cutoff for the $NN$ sector of the interaction and $n$ is an integer. Then, the interaction $V(\mathbf{p},\mathbf{p}^{\prime})$ is regulated as

[TABLE]

Even when these regulators are applied to a local interaction $V(\mathbf{p},\mathbf{p}^{\prime})=V(\mathbf{q})$, e.g., a momentum-independent contact interaction or the local one-pion-exchange interaction, the regularized interaction becomes nonlocal due to the explicit $\mathbf{k}$ dependence of the regulator functions.

A possible solution is to introduce local short- and long-range regulators. In our case, we regulate the chiral interactions directly in coordinate space. Short-range contact interactions, which Fourier transform to $\delta$ functions in coordinate space, are regulated by “smearing them out,” i.e.,

[TABLE]

In this work, we choose $n=4$. The constant $R_{0}$ serves as a coordinate-space cutoff parameter. The normalization is chosen such that

[TABLE]

For the long-range parts of the interaction, we use a similar functional form:

[TABLE]

In Fig. 1, we compare the short-range regulator used in the local chiral interactions for two values of the cutoff parameter $R_{0}$ with the short-range part used in the Argonne $v_{18}$ interaction Wiringa et al. (1995). Specifically, the short-range part of the Argonne $v_{18}$ interaction is given by $[P^{i}_{ST,{NN{}}}+\mu rQ^{i}_{ST,{NN{}}}+(\mu r)^{2}R^{i}_{ST,{NN{}}}]W(r)$, where $\mu$ is the average pion mass; $P$, $Q$, and $R$ are a set of parameters; and $W(r)$ is a Woods-Saxon potential. We display this short-range part of the Argonne $v_{18}$ interaction in the central channel for deuteron-like pairs, $i=c$ (central), $ST=01$, and $NN{}=np$, and normalize as in Eq. (7); see Ref. Wiringa et al. (1995) for details on the values of the parameters $P$, $Q$, $R$, and $\mu$ and the Woods-Saxon potential $W(r)$.

Regarding the range of cutoff parameters, one would like to take $R_{0}$ as small as possible in coordinate space to minimize regulator artifacts. However, as has been argued in Ref. Baru et al. (2012) in the context of the multiple-scattering series, the chiral expansion for the pion-exchange potentials breaks down for distances of $r\sim 0.8$ fm. For $r\gtrsim 1.0$ fm, the convergence of the multiple-scattering series, however, is found to be rather fast. Taking $R_{0}$ to be arbitrarily large, on the other hand, cuts off long-range pion physics that is resolved. We therefore adopt the range 1.0–1.2 fm for the cutoff $R_{0}$.

Although we stress that there is no direct correspondence between coordinate- and momentum-space cutoffs, a possibility of comparing the coordinate-space cutoff $R_{0}$ with typical momentum-space cutoff parameters $\Lambda_{NN{}}$ can be obtained by Fourier transforming the coordinate-space regulator function Eq. (6), integrating over all momenta, and identifying the result with a sharp cutoff. This gives $\Lambda_{NN{}}=\hbar c[6\pi^{2}\delta_{R_{0}}\!(0)]^{1/3}$, and thus we identify the corresponding momentum scales $\sim 500$ MeV with $R_{0}=1.0$ fm, and $\sim 400$ MeV with $R_{0}=1.2$ fm. While a clear translation between coordinate-space and momentum-space cutoffs can only be obtained when looking at a particular system or channel, we note that the estimated range encompassed by our cutoff choice is typical of other nonlocal chiral EFT interactions; see also Hoppe et al.

Regarding the long-range regulator, there are additional advantages in choosing a local regulator function. As has been argued recently Epelbaum et al. (2015a), the standard regulator choice Eq. (4) distorts the analytic structure of the partial-wave amplitude near threshold. Since the long-range interactions in chiral EFT are local [with the exception of relativistic corrections entering at next-to-next-to-next-to-leading order (N3LO)], it is logical to employ a local regulator in coordinate space, which cuts off the short-range part of the pion-exchange interactions but leaves the long-range part undisturbed. For this reason, a (different) local long-range regulator function is also chosen in the semilocal interactions of Epelbaum et al. Epelbaum et al. (2015a, b).

To regularize pion loops in the TPE contributions at NLO and higher orders, we use the framework of spectral function regularization (SFR). In SFR, the integrals over loop momenta in the spectral representation of the TPE contributions are cut off at $\tilde{\Lambda}$. In the following, we use the SFR cutoff $\tilde{\Lambda}=1000$ MeV since only a negligible dependence on its choice was found Gezerlis et al. (2014); Lynn et al. (2014). In particular, increasing the SFR cutoff from 1 to 1.4 GeV lowered the \isotope[4]He binding energy and the energy per particle of pure neutron matter (with only $NN$ interactions in both cases) by less than $\sim 2\%$, which is well within the $\sim 5\%$ truncation uncertainty at this order.

II.1.2 Local contact operators

Choosing local regulators removes the first source of nonlocality in chiral interactions. The second source of nonlocality originates in the momentum dependence of higher order contact interactions. Up to N2LO, these can be eliminated by exploiting Fierz ambiguities. At next-to-leading order (NLO), i.e., $Q^{2}$ in the chiral expansion, the general set of contact operators consistent with all the symmetries contains 14 different operators. In addition to spin-isospin dependences, these operators contain momentum dependences of the form $q^{2}$ and $k^{2}$ or $\mathbf{q}\times\mathbf{k}$, where the $k^{2}$ dependences are undesirable for local interactions. One can show using the Pauli principle that between antisymmetric states only 7 out of the 14 operators are linearly independent. Six linearly independent operators can be chosen to be local ($q^{2}$ dependent) while the 7th operator can be chosen to be the spin-orbit interaction; see Ref. Gezerlis et al. (2014) for more details.

At N3LO, there are an additional 15 linearly independent contact operators. Only 8 of these are local, while the other 7 operators contain $\mathbf{k}$ dependences that cannot be removed. Nevertheless, it is possible to construct maximally local N3LO interactions that contain, at most, nonlocalities of second order in momentum; see Ref. Piarulli et al. (2015) for initial work in this direction. To summarize, by choosing an appropriate set of contact operators and local regulator functions, all sources of nonlocality in chiral EFT can be removed up to N2LO.

II.1.3 Uncertainty estimates

To estimate the truncation uncertainty of the chiral expansion, we follow Ref. Epelbaum et al. (2015a) and estimate the uncertainty of an observable $X$ at N2LO as

[TABLE]

and correspondingly at lower orders. Furthermore, we require the uncertainties to be at least the size of the actual higher order corrections. We define the scale $Q$ as $Q=\max(p/\Lambda_{b},m_{\pi}/\Lambda_{b})$ with $p$ being a typical momentum scale of the system. For the work we present below, for nuclei, we choose $Q=m_{\pi}/\Lambda_{b}$, whereas for our neutron matter results, we take $Q$ from the average momentum in a Fermi gas $Q=\sqrt{3/5}k_{F}/\Lambda_{b}$, with Fermi momentum $k_{F}$; see Ref. Lynn et al. (2016). This choice is conservative, because typical binding momenta in nuclei are smaller than the pion mass. These uncertainty estimates provide a quantitative estimate of the effect of truncating the chiral expansion at some order $\nu$. A careful statistical analysis using Bayesian procedures has been undertaken in Ref. Furnstahl et al. (2015), where it was shown that the prescription we use, first introduced in Ref. Epelbaum et al. (2015a), results in $\nu/(\nu+1)\times 100\%$ degree-of-belief (DOB) intervals. That is, our NLO and N2LO uncertainty estimates are equivalent to $50\%$ and $\sim 67\%$ DOB intervals.

Further details of the $NN$ interaction, e.g., on the inclusion of charge-independence and charge-symmetry breaking terms, the values of the fitted low-energy constants (LECs), and phase shifts, are given in Ref. Gezerlis et al. (2014).

II.2 Deuteron properties

The deuteron is the lightest nucleus with $A>1$ in nature and provides a natural testing ground for the $NN$ interaction. In this section, we present some properties of this simple system using chiral interactions at N2LO. The deuteron wave function can be written in terms of its $S$- [$u(r)$] and $D$-wave [$w(r)$] components as

[TABLE]

where $\chi_{M_{J}}$ is the spin wave function for the total angular momentum projection $M_{J}$, and $S_{ik}(\mathbf{r})=3{\bf\it\sigma}_{i}\dotproduct\mathbf{\hat{r}}{\bf\it\sigma}_{k}\dotproduct\mathbf{\hat{r}}-{\bf\it\sigma}_{i}\dotproduct{\bf\it\sigma}_{k}$ is the tensor operator. The $S$- and $D$-wave components are normalized such that

[TABLE]

The $S$- and $D$-wave components in momentum space are then related by Fourier-Bessel transforms

[TABLE]

(where $j_{l}(x)$ is a spherical Bessel function) so that the normalization is

[TABLE]

We show the $S$- and $D$-wave components of the deuteron wave function in Fig. 2 for chiral interactions at N2LO with two different cutoff scales along with the deuteron wave function for the Argonne $v_{18}$ interaction. Compared to the hard Argonne $v_{18}$ interaction, the $S$-wave components of the local chiral interactions are softer, reflected in the larger value at vanishing pair separation $r$. As a result, the $D$-wave component is pushed away from $r=0$. In addition, the $D$-wave component at N2LO with cutoff $R_{0}=1.0$ fm ($R_{0}=1.2$ fm) has a node at $\sim 0.2$ fm ($\sim 0.02$ fm). This node has no physical consequences for the deuteron structure and for both cutoffs occurs at very short distances, where the uncertainty coming from the truncation of the chiral expansion is largest.

In Table 1, we collect a number of properties of the deuteron at N2LO and compare with experiment. The deuteron binding energy is not used in fits of the LECs and can be used as a check for the local potentials. At N2LO, the deuteron binding energy is consistent with experiment taking into account the uncertainties.

II.2.1 Momentum distribution

The deuteron momentum distribution can be written in terms of the $S$- and $D$-wave components as

[TABLE]

so that the normalization is

[TABLE]

In Fig. 3, we show the deuteron momentum distribution for our two cutoff choices along with the momentum distribution obtained for the Argonne $v_{18}$ interaction. It is interesting to note that the three momentum distributions display very similar behavior up to the respective cutoffs of the two chiral interactions. For $R_{0}=1.0\text{ fm}\sim 500\text{ MeV}\approx 2.5\text{ fm}^{-1}$, the blue curve begins to deviate significantly from the Argonne $v_{18}$ result at momenta $\sim 2.5\text{ fm}^{-1}$, while for $R_{0}=1.2\text{ fm}\sim 400\text{ MeV}\approx 2.0\text{ fm}^{-1}$, the red curve begins to deviate significantly from the Argonne $v_{18}$ result at $\sim 2.0\text{ fm}^{-1}$. However, we also emphasize that momentum distributions are necessarily renormalization scale and scheme dependent and are thus not observable Furnstahl and Hammer (2002).

II.2.2 Tensor polarization

Since momentum distributions are scheme and scale dependent, we now consider the tensor polarization. The charge form factors for different $M_{J}$ states are given by

[TABLE]

with the deuteron two-body density $\rho_{d}^{(M_{J})}(\mathbf{r}^{\prime})$ in state $M_{J}$ in terms of the distance $\mathbf{r}^{\prime}$ from the center of mass:

[TABLE]

The functions $C_{0}$ and $C_{2}$ are in turn written in terms of the $S$- and $D$-wave components of the deuteron wave function:

[TABLE]

The tensor polarization $T_{20}(q)$ is defined (in the impulse approximation) by Forest et al. (1996)

[TABLE]

We compare the tensor polarization for both cutoffs and for the Argonne $v_{18}$ interaction with experimental data Abbott et al. (2000) in Fig. 4. The first minimum of $T_{20}(q)$ is experimentally known at $q~{}\approx~{}3.5(5)\text{ fm}^{-1}$ Forest et al. (1996); Abbott et al. (2000); Zhang et al. (2011), in agreement with the predictions of all three cases displayed. At higher values of $q$, we expect some disagreement between our calculations and experiment given that we work in the impulse approximation.

II.2.3 Perturbation-theory calculations

The chiral expansion is meant to be a perturbative expansion in powers of a small parameter $Q\sim p/\Lambda_{b}$. One may well ask if the expected perturbative expansion is evident in the interactions themselves. To investigate this, we treat the difference of the N2LO and the NLO interactions as a perturbation

[TABLE]

and perform first-, second-, and third-order perturbation-theory calculations for the deuteron binding energy. For example, at first order,

[TABLE]

These results at second order and above are displayed in Fig. 5. As is evident from the figure, the softer interaction with $R_{0}=1.2$ fm is more perturbative than the harder interaction with $R_{0}=1.0$ fm. In both cases, the perturbative series appears to be converging to the value at N2LO, but the convergence is faster for the $R_{0}=1.2$ fm cutoff.

II.3 Three-nucleon interactions at

N2LO

Phenomenological models for $3N$ interactions, including the Urbana Pudliner et al. (1997), Illinois Pieper et al. (2001), and Tucson-Melbourne Coon and Han (2001) models, have been very successfully used in QMC calculations of nuclear systems. These models are based on the $3N$ TPE interaction that was first proposed by Fujita and Miyazawa nearly 60 years ago Fujita and Miyazawa (1957). Despite their undeniable success, they suffer from several shortcomings: They do not emerge naturally from the phenomenological $NN$ interactions and they are not systematically improvable.

In chiral EFT, however, $3N$ interactions naturally emerge in the expansion and are consistent with the $NN$ interactions. Furthermore, they are systematically improvable. The leading $3N$ interactions appear at N2LO in Weinberg power counting and can be visualized in terms of the diagrams in Fig. 6. The first diagram, proportional to the pion-nucleon LECs $c_{1}$, $c_{3}$, and $c_{4}$, corresponds to the long-range $S$- and $P$-wave TPE interactions by Fujita and Miyazawa. The LECs $c_{i}$ already appear in the subleading TPE interactions at the $NN$ level at the same chiral order, which highlights the consistency of the $NN$ and $3N$ interactions in chiral EFT. The second diagram, proportional to the LEC $c_{D}$, is an intermediate-range one-pion-exchange–contact interaction, and the third diagram, proportional to the LEC $c_{E}$, is a $3N$ contact interaction.

The diagrams in Fig. 6 give rise to the following momentum-space $3N$ interactions:

[TABLE]

where Roman indices refer to nucleon number, Greek indices refer to Cartesian coordinates, $\pi(ijk)$ gives all permutations of the indices, $g_{A}$ is the axial-vector coupling constant, $F_{\pi}$ is the pion decay constant, $\Lambda_{\chi}$ is taken to be a heavy meson scale, and $m_{\pi}$ is the pion mass. The function $F_{ijk}^{\alpha\beta}$ is defined in Ref. Tews et al. (2016) and depends on the LECs, $c_{1}$, $c_{3}$, and $c_{4}$. The two LECs $c_{D}$ and $c_{E}$ first appear in the $3N$ sector at N2LO and have to be fitted to $A\geq 3$ experimental data. We discuss our fitting procedure further below.

II.3.1 Local $\mathit{3N{}}$

interactions

The Fourier transformations of Eqs. (22a) to (22c) can be found in Ref. Tews et al. (2016). Here, we briefly review some important details from that work and point out additional details that arose in the implementation of the coordinate-space interactions in finite nuclei and neutron matter.

In commonly used phenomenological models, any short-range structures which arise in the Fourier transformation of long-range parts of the $3N$ forces are typically absorbed by other short-range structures (e.g., the scalar short-range structure in the Urbana IX (UIX) $3N$ interaction): However, we retain these additional structures explicitly. Our regularization scheme for the $3N$ interactions is consistent with that used in the $NN$ sector, i.e., $\delta$ functions denoting contact interactions are replaced with Eq. (6), long-range pion-exchange interactions are regulated by applying Eq. (8), and the $3N$ cutoff parameter is taken in the same range as the $NN$ cutoff parameter (in the following, we choose $R_{3N{}}=R_{0}=1.0-1.2$ fm). The full Fourier transformations of Eq. (22a) are available in Ref. Tews et al. (2016), but we note that a compact form of $V_{C,c_{3}}^{ijk}$ and $V_{C,c_{4}}^{ijk}$ can be obtained by writing them in the form of an anticommutator and a commutator of a modified coordinate-space pion propagator

[TABLE]

See the Appendix for details.

II.3.2 Regulator artifacts

As was discussed in Refs. Tews et al. (2016); Lynn et al. (2016); Dyhdalo et al. (2016), the use of local regulators in the $3N$ sector leads to two kinds of observable regulator artifacts. The first kind of regulator artifact affects the short-range parts of the interactions in Eqs. (22b) and (22c). These parts retain additional ambiguities at finite cutoff $R_{3N{}}\neq 0$. The first ambiguity concerns the choice of momentum variables in the Fourier transformation. Depending on how this choice is made, Eq. (22b) Fourier transforms to one of the following two equations:

[TABLE]

where $X_{ik}(\mathbf{r})=[S_{ik}(\mathbf{r})T(r)+{\bf\it\sigma}_{i}\dotproduct{\bf\it\sigma}_{k}]Y_{ik}(r)$ is the coordinate-space pion propagator, and the tensor and Yukawa functions are defined as $T(r)=1+3/(m_{\pi}r)+3/(m_{\pi}r)^{2}$ and $Y(r)=\mathrm{e}^{-m_{\pi}r}/r$. The sum with $i<j<k$ runs over all particles 1 to $A$, and the cyclic sum runs over the cyclic permutations of a given triple. It is clear that in the limit $R_{3N{}}\to 0$ the two possible $V_{D}$ structures are identical, because then the $\delta$ functions enforce $i=j$ $(k=j)$ in the first (second) term. The interaction $V_{D}$ does not distinguish which of the two nucleons in the contact interaction participates in the pion exchange. The term $V_{D2}$ can also be obtained by imagining a heavy fictitious scalar particle being exchanged between the two nucleons in the contact; see Fig. 7. This ambiguity was already pointed out in Ref. Navrátil (2007).

The second ambiguity in the $3N$ short-range interactions relates to the choice of the contact operator in Eq. (22c). The same Fierz-rearrangement freedom that allows for a selection of local contact operators entering in the $NN$ sector at NLO also allows for the selection of one out of the following six operators in the $3N$ sector Epelbaum et al. (2002):

[TABLE]

The usual choice is ${\bf\it\tau}_{i}\dotproduct{\bf\it\tau}_{j}$. This Fierz-rearrangement freedom holds as long as the regulator is symmetric under individual nucleon permutations. However, in the presence of local regulators, the Fierz-rearrangement freedom is violated, and different operator choices can lead to different results. Corrections to the violated Fierz rearrangement freedom are of higher order in chiral EFT. A systematic study of these effects in the $NN$ sector is in preparation Huth et al. . In the following, we have explored three different choices for the contact operator:

[TABLE]

The first two operator structures are chosen because $\mathbbm{1}$ and ${\bf\it\tau}_{i}\dotproduct{\bf\it\tau}_{j}$ have opposite signs in light nuclei but the same sign in neutron matter and thus give an estimate of the uncertainty due to this ambiguity. The last choice contains the projection operator $\mathcal{P}$ that projects on to triples with $S=\tfrac{1}{2}$ and $T=\tfrac{1}{2}$,

[TABLE]

where the sums are over pairs in a given triple. These are the triples that survive in the limit $\delta_{R_{3N{}}}\!(r)\to\delta^{(3)}(\mathbf{r})$, that is, the limit $R_{3N{}}\to 0$ (or $\Lambda\to\infty$ in momentum space).

The second regulator artifact affects the long-range $3N$ TPE interaction. It has been found that the effective $3N$ cutoff for a local regulator is lower (in momentum space) than for a typical nonlocal regulator Tews et al. (2016); Dyhdalo et al. (2016). As a consequence, one finds less repulsion from a local $3N$ TPE interaction than for the standard nonlocal formulation. This, again, is a regulator artifact that vanishes when $R_{3N{}}\to 0$. Lowering the $3N$ cutoff well below the $NN$ cutoff, however, leads to collapses because the increasing $3N$ attraction cannot be counteracted by additional $NN$ repulsion; see Ref. Tews et al. (2016).

II.3.3 Fitting procedure

We now turn to the fitting procedure for the LECs $c_{D}$ and $c_{E}$. This procedure was presented and discussed in Ref. Lynn et al. (2016), but we review it here for completeness. In the past, the binding energies of \isotope[3]H and \isotope[4]He or the binding energy of \isotope[3]H and the $nd$ doublet scattering length ${}^{2}a_{nd}$ have been used to fix $c_{D}$ and $c_{E}$. However, these observables are correlated and thus underconstrain the two LECs Gazit et al. (2009). The $3N$ couplings have also been fit to the \isotope[3]H binding energy and the \isotope[4]He radius Hebeler et al. (2011). Arguments can be made that $3N$ interactions should be fit in $A\leq 3$ systems only Gazit et al. (2009) or that reproducing observables over a wider range in the nuclear chart is more appropriate Pieper et al. (2001); Ekström et al. (2015). We take a middle-ground approach and have two goals with our fitting strategy: (1) to probe properties of light nuclei and (2) to probe $T=3/2$ physics. With these in mind, we take as observables the \isotope[4]He binding energy and $n$-$\alpha$ scattering $P$-wave phase shifts. The $n$-$\alpha$ system is the lightest nuclear system for which three neutrons can be found interacting and thus provides an indirect constraint on $T=3/2$ physics.

We first find contours for $c_{D}$ and $c_{E}$ that reproduce the \isotope[4]He binding energy. We further constrain $c_{D}$ and $c_{E}$ by calculating the $P\ 3/2^{-}$ and $P\ 1/2^{-}$ phase shifts for the $n$-$\alpha$ system as described in Ref. Nollett et al. (2007) and demanding a good reproduction of the splitting between these two $P$-wave phase shifts. See Ref. Lynn et al. (2016) for more details.

In Ref. Lynn et al. (2016), we explored the various combinations of $V_{D}$ [Eqs. (24a) and (24b)] and $V_{E}$ [Eqs. (26a)–(26c)] and found some dependence on these choices. In particular, no fit to both observables (the \isotope[4]He binding energy and the $n$-$\alpha$ $P$-wave scattering phase shifts) was obtained for the case with $V_{D1}$ and the softer cutoff $R_{0}=1.2$ fm. For all other combinations, results for light nuclei with $A=3,4$ were similar. Below we take a representative choice ($V_{D2},V_{E\tau}$) for the results we display.

III Quantum Monte Carlo methods

In this section, we provide details on QMC methods including the variational Monte Carlo (VMC) method, which is used as a starting point for both GFMC and AFDMC calculations.

III.1 Variational Monte Carlo

The variational Monte Carlo (VMC) method relies on the Rayleigh-Ritz variational principle:

[TABLE]

where $\ket{\Psi_{T}}=\ket{\Psi_{T}(\{c_{i}\})}$ is a trial wave function with a set of adjustable parameters $\{c_{i}\}$, and $E_{0}$ is the energy of the ground state of $H$. The equality above only holds if $\ket{\Psi_{T}}=\ket{\Psi_{0}}$, the ground state of $H$.

For few-body nuclei with $A=3,4$ the form of the variational trial wave function is given as

[TABLE]

The two-body part of the wave function consists of a symmetrized product of correlation operators acting on a Jastrow wave function,

[TABLE]

where the Jastrow wave function is

[TABLE]

The Jastrow factor is a product of central two-body correlations $f_{c}$ acting on an appropriate antisymmetric single-particle state. For few-body nuclei with $A=3,4$, $\ket{\Phi}$ can be taken as an appropriate antisymmetric linear combination of spin-isospin states. For example, for \isotope[4]He, one can take

[TABLE]

The correlation operators are defined as

[TABLE]

where the $\{O_{ij}^{(p)}\}$ are the operators

[TABLE]

taken from the two-body interaction. We use the short-hand notation $p=\{\sigma,\tau,\sigma\tau,t,t\tau,b\}$ for the operators as in Ref. Wiringa (1991). The symmetrizer in Eqs. (29) and (30) is necessary to maintain the overall antisymmetry of the wave function, since in general the $U_{ij}$ do not commute with each other. In Fig. 8, we display the two-body correlations $f_{c}$, and $\{u_{p}\}$ obtained in the simulation of \isotope[4]He with the N2LO interactions with both cutoffs as well as those obtained for the Argonne $v_{18}$ $NN$ interaction supplemented by the UIX $3N$ interaction. What can be seen from these correlations, most particularly in the case of the central correlation $f_{c}$, is the softening of the interaction as we take the cutoff from $R_{0}=1.0$ fm to $R_{0}=1.2$ fm. We find that the spin-orbit correlation has only a minimal effect on the variational energies we obtain and a relatively high computational cost, and therefore we set $u_{b}(r_{ij})=0$ in our calculations.

The three-body correlation operator takes the following form:

[TABLE]

where $\bar{r}$ is a scaled relative separation and $\epsilon$ is a small negative constant. This form is suggested by perturbation theory Pudliner et al. (1997). In addition to the explicit three-body correlations of Eq. (35), a central, geometric three-body correlation is wrapped into the two-body correlations,

[TABLE]

with

[TABLE]

where $R_{ijk}=r_{ij}+r_{jk}+r_{ik}$ and the $\{t_{i}\}$ are variational parameters. These correlations serve to reduce the repulsion which arises from the product of certain spin-isospin correlation operators when any two nucleons come close together. Reducing this repulsion was found to improve variational energies with wave functions of the form of Eq. (30) Lomnitz-Adler et al. (1981).

Equation (28) is evaluated by means of Monte Carlo integration,

[TABLE]

where $a$ and $b$ stand for a given order of operators in the product Eq. (30), a complete sum over all spin and isospin states is assumed, and the integrals are performed as a Monte Carlo integration over the coordinate-space configurations $\mathbf{R}=\{\mathbf{r}_{1},\mathbf{r}_{2},\ldots,\mathbf{r}_{A}\}$. The sums over the orders $a$ and $b$ are also performed via a Monte Carlo sampling as discussed below. The weight function can be taken as

[TABLE]

for example. In practice, because of the different orders $a$ and $b$ in the left and right wave functions, Eq. (39) is not guaranteed to be nonzero, and so we add to it an additional term proportional to $\sum_{s,t}|\Psi_{a}(\mathbf{R};s,t)^{\dagger}\Psi_{b}(\mathbf{R};s,t)|$. That is, we add a term proportional to the sum of the absolute value of the overlaps of the individual spin-isospin components of the wave functions.

The symmetrizer in Eqs. (29) and (30) requires, in principle, the evaluation of all $[A(A-1)/2]!$ possible orderings of the operators. To save computational cost, the order of operators is instead sampled. This approximation does not contribute much to the statistical variance since all orderings share the same linear (dominant) contributions and the differences between different orderings are proportional to $\{u_{p}^{2}\}$.

The Metropolis algorithm is employed and the result is, after sufficient equilibration, a set of configurations labeled by the $3A$ coordinates and the orderings of the operators, $\{\mathbf{R},a,b\}$, which are distributed according to the square of the trial wave function. As the integration and sum over all orderings is done stochastically, there is an error associated with the expectation value of any operator $\expectationvalue{O}$, given by

[TABLE]

where $N$ is the number of statistically independent evaluations. For more details, see Refs. Wiringa (1991); Carlson et al. (2015).

With the algorithm described above, the variational parameters $\{c_{i}\}$ are adjusted to minimize the expectation value of the Hamiltonian in Eq. (28). Wave functions so obtained can be used as reasonable approximations to the exact ground state (especially in few-body nuclei) and are a necessary starting point for the GFMC method.

III.2 Diffusion Monte Carlo

Even with the sophisticated wave functions described in Sec. III.1, it is not possible to construct by hand exact solutions to the many-body Schrödinger equation. Diffusion Monte Carlo methods including the AFDMC and GFMC methods rely on the fact that, given a nuclear system specified by the Hamiltonian $H$ with ground state $\ket{\Psi_{0}}$ and a trial wave function for that system $\ket{\Psi_{T}}$ with nonvanishing overlap with the ground state,

[TABLE]

The object $\mathrm{e}^{-H\tau}$ is the many-body imaginary-time Green’s function (or imaginary-time propagator) for the system with the imaginary time $\tau$. This “sifting” property of the imaginary-time propagator is easy to understand if the trial wave function is expanded in a complete set of eigenstates of $H$, $\{\ket{\phi_{n}}\}$, with energies $\{E_{n}\}$,

[TABLE]

where we have introduced the trial energy $E_{T}$ and $a_{n}=\langle\phi_{n}|\Psi_{T}\rangle$. In principle, $E_{T}$ can take any value, but it is often adjusted to be the ground-state energy (or the energy of the low-lying excited state sought). Then, since $E_{i}>E_{T}$ for all $i>0$, in the large-imaginary-time limit all of the excited-state components of the trial state are exponentially damped and one is left with the exact many-body ground state. In this language, we can say that with the VMC method alone it is not possible to avoid some contamination in nuclear wave functions from excited states. That is, while we can make $a_{0}$ of Eq. (42) the dominant contribution through the adjustment of the variational parameters $\{c_{i}\}$, it is not possible with the VMC method alone to guarantee that $a_{n>0}=0$.

In the remainder of this section, we discuss diffusion Monte Carlo methods, paying particular attention to the GFMC method, which we use to calculate properties of light nuclei. For more details, we refer to Ref. Carlson et al. (2015) and references therein. We begin with a discussion of the calculation of the imaginary-time propagator, which plays a central role in diffusion Monte Carlo methods.

In general, it is difficult to compute the exact many-body imaginary-time propagator for arbitrary imaginary times. Instead, the properties of the exponential are exploited to rewrite the propagation to large imaginary time as a product of small propagations,

[TABLE]

with $\Delta\tau=\tau/N$, and $N$ large enough ($\Delta\tau$ small enough) such that one of several approximations can be used to calculate the short-imaginary-time propagator. In the case of the AFDMC method, for example, a Trotter breakup is used Carlson et al. (2015),

[TABLE]

where $T$ is the kinetic energy operator and $V_{ij}$ is a local two-body interaction. In Eq. (44), the order in the product on the left is taken in the opposite order of the product on the right. This keeps the propagator unitary (in real time) and eliminates terms of $\mathcal{O}(\Delta\tau^{2})$.

In the GFMC method, the exact two-body propagator is used to construct the many-body propagator, as suggested by studies of condensed helium systems Ceperley (1995):

[TABLE]

Here, $\alpha$ and $\beta$ stand for the appropriate spin-isospin states for a given nucleus, $\mathbf{R}=\{\mathbf{r}_{1},\mathbf{r}_{2},\ldots,\mathbf{r}_{A}\}$ and $\mathbf{R}^{\prime}=\{\mathbf{r}^{\prime}_{1},\mathbf{r}_{2}^{\prime},\ldots,\mathbf{r}^{\prime}_{A}\}$ are the collections of $3A$ coordinates before and after the propagation step, $G_{0}(\mathbf{R},\mathbf{R}^{\prime};\Delta\tau)$ is the many-body free-particle imaginary-time propagator

[TABLE]

$g_{ij}$ is the exact two-body interacting imaginary-time propagator,

[TABLE]

which can be computed to high accuracy ($\sim$8- to 10-digit accuracy or better than half machine precision), and $g_{0,ij}$ is the two-body free-particle analog of $g_{ij}$. This construction allows for taking much larger time steps than in the Trotter breakup in Eq. (44). The trade-off is that the calculation of the exact two-body propagator of Eq. (47) is too costly to compute “on the fly” and must be carried out in advance and stored on a grid of points to be interpolated on during the GFMC propagation.

The complete two-body propagator depends on initial and final relative coordinates, the initial and final spin states of the pair, and the isospin of the pair,

[TABLE]

where the indices $ij$ as in Eq. (47) are suppressed here and in what follows for simplicity unless they are explicitly needed for clarity. Reference Schmidt and Lee (1995) originally proposed using fast Fourier transforms (FFT) and the Trotter break up for scalar interactions, and this idea was adapted to realistic nuclear interactions in Ref. Pudliner et al. (1997). In this method, interactions are first decomposed into partial waves. The nuclear Hamiltonian commutes with the operators $J^{2}$, $J_{z}$, $S^{2}$, $T^{2}$, and $T_{z}$, and, thus, sets them as good channel quantum numbers: $\mathbf{S}=\mathbf{S}_{1}+\mathbf{S}_{2}$ is the total spin, $\mathbf{J}=\mathbf{L}+\mathbf{S}$ is the total angular momentum, and $\mathbf{T}=\mathbf{T}_{1}+\mathbf{T}_{2}$ is the total isospin. Then, the channel propagators $\matrixelement{r^{\prime}JM_{J}L^{\prime}STM_{T}}{\mathrm{e}^{-H\Delta\tau}}{rJM_{J}LSTM_{T}}$ are computed and resummed to obtain the two-body propagator:

[TABLE]

Here, $\gamma$ stands for the set of quantum numbers $\{JMLL^{\prime}M_{L}M_{L}^{\prime}\}$, $C$ is a Clebsch-Gordan coefficient, $Y$ is a spherical harmonic, and $\Omega$ ($\Omega^{\prime}$) are the angular coordinates of $\mathbf{r}$ ($\mathbf{r}^{\prime}$).

Each of the channel propagators is calculated by breaking up the (already-small) time step $\Delta\tau$ into smaller steps $\delta\tau=\Delta\tau/N_{\tau}$, with $N_{\tau}$ large, using the symmetrized Trotter breakup, and FFT:

[TABLE]

In Eq. (50b), the right-most exponential acts upon an array of initial relative separations, the result is transformed to momentum space using FFT, the exponential of the kinetic energy acts upon that result, which is then transformed back to coordinate space using FFT, whereupon the left-most exponential acts upon the array. This method introduces errors of $\mathcal{O}(\delta\tau^{3})$, is fast, and is easy to implement.

An alternative method is to diagonalize the channel Hamiltonians in momentum space Lynn and Schmidt (2012). When the interaction is nonlocal (no longer diagonal in coordinate space), then the advantages of the Trotter breakup vanish. That is, it is just as difficult to calculate the matrix elements $\matrixelement{\mathbf{r}^{\prime}}{\mathrm{e}^{-V\Delta\tau}}{\mathbf{r}}$ as it is to calculate the original matrix elements $\matrixelement{\mathbf{r}^{\prime}}{\mathrm{e}^{-H\Delta\tau}}{\mathbf{r}}$. In order to diagonalize the channel Hamiltonians, we take as an orthonormal basis the set of spherical Bessel functions which solve the free radial Schrödinger equation with a Dirichlet boundary condition at some radius $R$ much beyond the range of the interaction,

[TABLE]

where $\{k_{n}\}$ is the set of discrete momenta for a given $L$ and $R$. In this basis, the kinetic energy is diagonal, and the potential-energy matrix elements can be obtained with simple matrix multiplications, which perform the necessary numerical integrals. While this method was originally developed to calculate two-body propagators for nonlocal interactions, it works equally well for local interactions, providing equal accuracy and speed when compared with the symmetrized Trotter break up with FFT.

So far, we have discussed only the contribution to the many-body propagator coming from $NN$ interactions. We include $3N$ interactions in the propagator as a symmetric linear approximation to $\mathrm{e}^{-V_{3N{}}\Delta\tau}$:

[TABLE]

where the sums $\sum_{p}V^{(p)}_{3N{}}$ are over all $3N$ operators of Eqs. (64a)–(64c), one of Eqs. (65a) and (65b), and one of Eqs. (66a)–(66c). As before, $\alpha$, $\beta$, $\gamma$, and $\delta$ are appropriate spin-isospin states and $\gamma$ and $\delta$ are summed over. This linear approximation is a controlled approximation that becomes more exact with smaller $\Delta\tau$. There are improvements to this linear approximation possible. For example, replacing $[\mathbbm{1}-\tfrac{\Delta\tau}{2}\sum_{p}V^{(p)}_{3N}(\mathbf{R})]$ with $\prod_{p}[\mathbbm{1}-\tfrac{\Delta\tau}{2}V^{(p)}_{3N}(\mathbf{R})]$ would capture at least some $\mathcal{O}(\Delta\tau^{2})$ effects. Another possibility is to include all parts of the $3N$ interaction that can be rewritten effectively as two-body operators into the two-body propagator as suggested in Ref. Wiringa et al. (2000) [these include the TPE $P$-wave anticommutator contribution Eq. (64b), the TPE $S$-wave contribution Eq. (64a), the $V_{D}$ contributions Eqs. (65a) and (65b), and two of the three $V_{E}$ contributions Eqs. (66a) and (66b)]. However, we have found that with the time step we typically use, $\Delta\tau=0.0005\text{ MeV}^{-1}$, the time-step error introduced by this linear approximation is negligible.

With the imaginary-time propagator so obtained, one would ideally like to calculate expectation values such as $\expectationvalue{O(\tau)}=\tfrac{\matrixelement{\Psi(\tau)}{O}{\Psi(\tau)}}{\innerproduct{\Psi(\tau)}{\Psi(\tau)}}$, with $\Psi(\tau)$ defined as

[TABLE]

However, in practice, one does not have direct access to the propagated wave function, and an evaluation of that expectation value is cumbersome for spin- and isospin-dependent operators, and especially for momentum-dependent operators. Thus, what is more commonly used is the mixed expectation value of a given operator (suppressing the spin-isospin indices), defined as

[TABLE]

with the paths $\differential{\bf\it\mathcal{R}}\equiv\prod_{i=0}^{N-1}\differential\mathbf{R}_{i}$, and the total imaginary time $\tau=N\Delta\tau$. The paths are Monte Carlo sampled to perform the integrals. Note that the operator $O$ must act on the trial wave function (to the left).

The mixed estimate introduces an explicit dependence on the trial wave function. However, if the trial wave function is a good approximation, we can write

[TABLE]

where $\delta\Psi(\tau)$ is the (small) correction to the trial wave function introduced by the imaginary-time propagation, and keep terms only of $\mathcal{O}\left(\delta\Psi(\tau)\right)$. Then we have

[TABLE]

where $\expectationvalue{O}_{T}$ is the variational estimate. Thus, if $\Psi_{T}$ is a good approximation to the exact wave function obtained through imaginary-time propagation (as measured by the relative smallness of the difference $\expectationvalue{O(\tau)}_{\text{mixed}}-\expectationvalue{O}_{T}$) when compared with $\expectationvalue{O(\tau)}_{\text{mixed}}$, then the mixed estimate introduces only a small systematic uncertainty. Typically we aim for the difference between the mixed and variational estimates to be no larger than $\sim 5\%$ of the mixed estimate. There are other ways to avoid the use of mixed estimates, such as computing the observable in the midpoint of the path Kalos (1966), but this requires a propagation time twice as long as in the mixed-estimate case. Note that in the case of the energy expectation value, $\expectationvalue{H}$, the Hamiltonian and the imaginary-time propagator commute. In this case,

[TABLE]

such that $\lim_{\tau\to\infty}\expectationvalue{H(\tau)}_{\text{mixed}}=E_{0}$. In short, for the Hamiltonian, the mixed estimate is identical to the normal estimate.

When performing the propagation, one has to employ another approximation. Nucleons are fermions and their many-body wave functions contain nodal surfaces. As a consequence, a configuration that crosses a nodal surface introduces a sign change in the matrix elements in Eq. (54). At large $\tau$, these sign changes contribute to a decreasing denominator, causing large statistical fluctuations (large variance). This is the famous fermion sign problem. One way to circumvent this problem is the so-called constrained path algorithm; for a detailed description, see Ref. Wiringa et al. (2000). In short, the idea is to discard configurations that in future propagations would only contribute to the variance. If one knew the exact wave function, then the overlap of these discarded configurations with the ground-state wave function would be zero $\innerproduct{\Psi_{\text{discarded}}}{\Psi_{0}}=0$. However, since we do not in general know the exact ground-state wave function, the constraint is imposed on the overlap with the trial wave function so that the average of the overlaps $\innerproduct{\Psi_{\text{discarded}}}{\Psi_{T}}$ over the random walk is approximately zero.This approximation was inspired by the fixed-node approximation used in condensed matter systems.

For scalar wave functions (no spin or isospin dependence) the fixed-node approximation provides both a way to tame the sign problem, and results in an upper bound to the ground-state energy. However, because of the spin and isospin dependence of the nuclear case, the constrained-path algorithm no longer supplies a strict upper bound, as has been discussed and demonstrated in Ref. Wiringa et al. (2000). To overcome this additional difficulty, in cases where the constrained-path algorithm is used, we take a number $n_{u}$ of unconstrained steps after convergence of the constrained-path calculation. We take $n_{u}$ as large as possible. Typically, $n_{u}\sim 20$ before the fermion sign problem overwhelms the signal. This “transient estimation” results in significantly improved estimates, introducing an error, for example, in \isotope[6]Li of just $\sim 0.5\%$; see Ref. Wiringa et al. (2000). Figure 9 gives an example of a constrained-path calculation of the ground-state energy of \isotope[4]He and the subsequent transient estimation. Note that the constrained-path propagation overbinds the system, demonstrating that for some trial wave functions the constrained-path estimate is not an upper bound.

IV Energies and other results for

$\boldsymbol{A=3,4}$

The light nuclei with $A=3,4$ are a minimal testing ground for any nuclear Hamiltonian: The reasonable reproduction of binding energies and radii in these nuclei is a basic yardstick against which our local chiral interactions can be measured.

In this section, we present the main results for light nuclei that we have obtained with our local chiral EFT $NN$ and $3N$ interactions at N2LO Lynn et al. (2016). We emphasize the order-by-order convergence of observables in the $A=3,4$ nuclei and present a detailed breakdown of the contributions to the energies from various components of the $NN$ and $3N$ interactions. We also show several one-body distributions and the related longitudinal charge form factor.

IV.1 Energies of light nuclei at

LO, NLO, and N2LO

At LO, the $NN$ interaction consists simply of the one-pion exchange potential and two contact interactions with LECs fit to $NN$ scattering phase shifts. With only the basic pion physics present and little freedom to fit the phase shifts, essentially only the large scattering length plus OPE physics can be reproduced, and the resulting potential is excessively attractive in low partial waves. This can be seen in Fig. 10, where at LO, the ground-state energies of $A=3,4$ nuclei are significantly lower than experiment. In fact, the LO $NN$ interaction overbinds by as much as $\sim 50$% ($\sim 30$%) for $A=4$ ($A=3$). At NLO, the $NN$ interaction is too repulsive and leads to underbinding. However, the deviation from experiment decreases to $\sim 25$% ($\sim 15$%) for $A=4$ ($A=3$). Finally, at N2LO, the $3N$ interaction with two free LECs enters. We fit $c_{D}$ and $c_{E}$ directly to the binding energy of \isotope[4]He, (see Sec. II.3) and, since the binding energies of the $A=3$ systems are highly correlated with the binding energy of \isotope[4]He (i.e., the Tjon line Tjon (1976); Platter et al. (2005)), the $A=3$ binding energies are also well reproduced.

The uncertainties in Fig. 10 contain contributions from the GFMC statistical uncertainties as well as from an estimate for the theoretical uncertainty coming from the truncation of the chiral expansion (as discussed in Sec. II). The theoretical uncertainties display at least three desirable features: (1) They encompass, order by order, the cutoff variation in the energy, (2) order by order, the experimental energy is within the uncertainty bands, and (3) as the chiral order increases, the uncertainty coming from the truncation of the chiral expansion decreases rapidly. Thus, at N3LO, we can expect that while the energies of these systems will not change dramatically, the uncertainties will continue to reduce.

In addition to Fig. 10, we show more details of various contributions to the $A=3,4$ energies in GFMC calculations in Table 2, where the softer nature of the interaction with $R_{0}=1.2$ fm is evident from the lower kinetic energies compared to the case with $R_{0}=1.0$ fm. Note, however, that the kinetic energy by itself is not an observable.

The trend represented in Fig. 10 is also present in the radii of the system: See Fig. 11 and Table 3. Here we compute the so-called point-proton radii of $A=3,4$ systems:

[TABLE]

where $(1+\tau_{z})/2$ is a projection operator onto protons and $Z$ is the number of protons. However, the measured charge radius includes effects from the charge densities of the finite-sized nucleons themselves. The relationship between the point-proton radius and the observable charge radius $r_{\text{c}}$ is given by

[TABLE]

where $\sqrt{\expectationvalue{r_{p}^{2}}}=0.8751(61)$ fm is the root-mean-square (rms) charge radius of the proton Patrignani et al. (2016), $N$ is the number of neutrons, $\expectationvalue{r_{n}^{2}}=-0.1161(22)\text{ fm}^{2}$ is the mean-square charge radius of the neutron Patrignani et al. (2016), and $m_{p}$ is the proton mass. The last term of Eq. (59) is the so-called Darwin-Foldy correction to the proton charge radius Friar et al. (1997). For larger $A$, there are also spin-orbit corrections to the charge radius Ong et al. (2010). The experimental charge radii are from Ref. Sick (2014) (\isotope[4]He and \isotope[3]He) and Ref. Amroun et al. (1994) (\isotope[3]H).

The correlation between the energies and radii of the nuclei are evident in Fig. 11 and Table 3. At LO, as the nuclei are significantly overbound, the point-proton radii are significantly smaller than the values extracted from experiment. At NLO, with the nuclei underbound, the point-proton radii are too large. At N2LO, with reasonable reproduction of the nuclear binding energies for the $A=3,4$ systems, the calculated point-proton radii are in good agreement within both the experimental and theoretical uncertainties. Note that the relatively large uncertainty quoted in the point-proton radius for \isotope[3]H extracted from experiment is due to the relatively large uncertainty in the charge radius for this nucleus: Compare $r_{\text{c}}(\isotope[3]{He})=1.973(14)$ fm with $r_{\text{c}}(\isotope[3]{H})=1.755(86)$ fm. The experimental uncertainty is roughly a factor of six larger for \isotope[3]H than for \isotope[3]He.

IV.2 More details on distributions

In addition to energies and radii, we have also calculated one-body distributions. The one-body point distributions are defined as

[TABLE]

with $N=p$ (taking the positive sign in the projector $\tfrac{1+\tau_{z}}{2}$) giving the point-proton distribution and $N=n$ (taking the negative sign in the projector $\tfrac{1-\tau_{z}}{2}$) giving the point-neutron distribution. When folded with the spatial proton charge distribution, the point proton distribution is promoted to the charge distribution, which is the Fourier transform of the charge form factor measured in electron scattering experiments. The short-distance behavior of the presented point-nucleon distributions are not as well constrained, because the high momentum-exchange charge form factor is challenging to measure and to calculate accurately. Nevertheless, the charge radius (or point-proton radius) as an integrated quantity is well constrained by experiment, and our results reproduce within uncertainties the point-proton radii extracted from experiment.

In Fig. 12, we show the point-proton distribution in \isotope[4]He for both cutoffs $R_{0}=1.0,1.2$ fm at N2LO with and without the $3N$ interaction. The corresponding point-proton radius is shown in a color-coded way on the right-hand side of the figure. Though it is not consistent from the EFT point of view to show the N2LO results without the $3N$ interaction, it is nevertheless instructive to see the effects of the $3N$ interaction in this way. One can see that its effect is to increase the density of protons at intermediate distances from the center of mass ($r\sim 1.0$ fm) while decreasing their density at short distances, yielding a peak at about $r~{}\sim 0.6$ fm. The effect of this shift is to bring the overall point-proton radius into better agreement with the number extracted from the experimental charge radius.

In Fig. 13 we show the one-body point-proton and neutron distributions for \isotope[3]He at N2LO for both cutoffs. At short distances from the center of mass, the distributions for the cutoff $R_{0}=1.2$ fm demonstrate a softer character: There is a higher probability of finding either a neutron or a proton at short distances from the center of mass than is the case for the distributions calculated with the $R_{0}=1.0$ fm cutoff. As is the case for \isotope[4]He, only the large-$r$ part of the distributions can be well constrained, and in this region, both cutoffs agree. We also show the corresponding point-proton radii with statistical GFMC uncertainties only, to demonstrate that integrated quantities such as the charge radius are essentially cutoff independent for these light systems at this order of the chiral expansion. Finally, in Fig. 13, one can see that the proton distribution is qualitatively twice the neutron distribution, but there are quantitative differences due to the presence of isospin-symmetry-breaking terms in the Hamiltonian.

The point-proton and point-neutron distributions we calculate are related to the experimentally observable electric charge form factor. In particular, the longitudinal electric charge form factor is given by

[TABLE]

where $\tilde{\rho}$ are the Fourier transforms of the point-nucleon distributions defined in Eq. (60), $G_{E}^{n,p}$ are the single nucleon electric charge form factors for the neutron $n$ and proton $p$ and $Q_{\text{el}}^{2}$ is the four-momentum squared:

[TABLE]

with

[TABLE]

Above, $m_{N}$ and $m_{A}$ are the average nucleon mass and the mass of the target nucleus, respectively. For the single-nucleon charge form factors $G_{E}^{n,p}$, we use the parametrizations of Kelly Kelly (2004), which enforce the correct asymptotic behavior as $Q^{2}_{\text{el}}\to 0$ and $Q^{2}_{\text{el}}\to\infty$.

In Fig. 14, we present the longitudinal electric charge form factor for \isotope[4]He compared with an unpublished compilation by Sick Sick of the data from Refs. Frosch et al. (1967); Erich et al. (1968); McCarthy et al. (1977); Arnold et al. (1978); Ottermann et al. (1985). The figure is log scaled as charge form factors are often plotted, but this scaling artificially enhances the apparent size of the uncertainties. However, the figure should be read as simply that at N2LO the uncertainty in the location of the first minimum in the \isotope[4]He charge form factor is roughly 0.6–0.8 fm*-1*. Note that calculations are performed without two-body currents, and thus the poorer comparison with data at higher $q$ is somewhat expected Carlson et al. (2015).

V Summary

In this paper, we presented additional details on and results for QMC calculations of light nuclei with local chiral $NN$ and $3N$ interactions. We discussed deuteron properties in detail, employing a soft and a hard local chiral interaction. We found that local chiral interactions give a reasonable description of the deuteron binding energy, rms radius, asymptotic $D/S$ ratio, and quadrupole moment. Furthermore, local chiral interactions reproduce the experimentally known first minimum of the deuteron tensor polarization.

We then performed perturbative calculations for both interactions in the deuteron, using the difference of the N2LO and NLO interactions as a perturbation around the NLO result. While both perturbative series seem to converge to the N2LO result, we found the softer interaction to be more perturbative, as expected.

We then presented additional details on our calculations of radii and binding energies of the light $A=3,4$ nuclei \isotope[3]H, \isotope[3]He, and \isotope[4]He. For each binding energy and radius and for both local chiral interactions, we observed an order-by-order convergence toward the experimental value. Finally, we discussed proton and neutron distributions for \isotope[3]He and \isotope[4]He.

Together with the results of Ref. Lynn et al. (2016), we have established QMC methods with local chiral interactions as a versatile tool to study properties of light nuclei and neutron matter.

Acknowledgements.

We thank S. Bacca, P. Klos, and D. Lonardoni for useful discussions. This work was supported by the ERC Grant No. 307986 STRONGINT, the National Science Foundation Grants No. PHY-1430152 (JINA-CEE) and No. PHY-1404405, the U.S. DOE under Grants No. DE-AC52-06NA25396 and No. DE-FG02-00ER41132, the NUCLEI SciDAC program, the Natural Sciences and Engineering Research Council (NSERC) of Canada, and the LANL LDRD program. Computational resources have been provided by the Jülich Supercomputing Center, the Lichtenberg high performance computer of the TU Darmstadt, and Los Alamos Open Supercomputing. We also used resources provided by NERSC, which is supported by the U.S. DOE under Contract No. DE-AC02-05CH11231.

Appendix A COMPLETE COORDINATE-SPACE EXPRESSIONS FOR THE

$\boldsymbol{3N{}}$ INTERACTION AT N2LO

As noted in Sec. II.3, the TPE parts of the $3N$ interaction $V_{C}$ can be compactly written in terms of the standard coordinate-space pion propagator $X_{ij}(\mathbf{r})$ (defined in that section) and a modified coordinate-space pion propagator $\mathcal{X}_{ij}(\mathbf{r})\equiv X_{ij}(\mathbf{r})-\frac{4\pi}{m_{\pi}^{2}}\delta_{R_{3N{}}}\!(r){\bf\it\sigma}_{i}\dotproduct{\bf\it\sigma}_{j}$. We also define the function $U(r)=1+1/(m_{\pi}r)$. Then, the complete TPE part of the interaction can be written as

[TABLE]

The remaining parts of the interaction are written as

[TABLE]

[TABLE]

We remind the reader that the projection operator $\mathcal{P}$ is defined in Eq. (27). We note that some differences exist between these expressions compared with those in Ref. Tews et al. (2016). Under the change $\sum_{\pi(ijk)}\to\sum_{\text{cyc}}$, Eqs. (64a) and (66a) to (66c) pick up an additional factor of 2 and Eq. (65a) picks up an additional term with $i\leftrightarrow k$. In addition, Eqs. (64b) and (64c) pick up factors of $\tfrac{1}{2}$ and $\tfrac{1}{2\mathrm{i}{}}$, respectively, from the replacements ${\bf\it\tau}_{i}\dotproduct{\bf\it\tau}_{j}=\tfrac{1}{2}\{{\bf\it\tau}_{i}\dotproduct{\bf\it\tau}_{k},{\bf\it\tau}_{k}\dotproduct{\bf\it\tau}_{j}\}$ and ${\bf\it\tau}_{i}\dotproduct({\bf\it\tau}_{j}\crossproduct{\bf\it\tau}_{k})=\tfrac{1}{2\mathrm{i}{}}[{\bf\it\tau}_{i}\dotproduct{\bf\it\tau}_{k},{\bf\it\tau}_{k}\dotproduct{\bf\it\tau}_{j}]$.