Uncertainties of Euclidean time extrapolation in lattice effective field theory

Several ab initio methods are currently being applied to the study of nuclear structure. These include coupled-cluster expansions [1], the no-core shell model [2, 3], the in-medium similarity renormalization group [4], self-consistent Green's functions [5], Green's function Monte Carlo [6], and auxiliary-field diffusion Monte Carlo [7].

Much of the recent progress in ab initio nuclear structure calculations is due to 'soft' chiral nuclear effective field theory (EFT) interactions. The lattice formulation of chiral nuclear EFT is described in [8], and a review of lattice EFT methods can be found in [9]. A comprehensive overview of chiral nuclear EFT is available in [10, 11]. This framework, known as nuclear lattice effective field theory (NLEFT), has been used to calculate the ground states of alpha nuclei from ⁴He to ²⁸Si, as well as to describe the structure of the Hoyle state [12, 13] and the dependence of the triple-alpha process on the fundamental parameters of nature [14].

NLEFT is an ab initio method, where chiral nuclear EFT is combined with auxiliary-field quantum Monte Carlo (AFQMC) lattice calculations. This Monte Carlo approach differs from other ab initio methods in that it does not require truncated basis expansions, many-body perturbation theory, or any constraint on the nuclear wave function. While our NLEFT results are thus unbiased Monte Carlo calculations, the AFQMC approach nevertheless relies on Euclidean time projection of projection Monte Carlo (PMC) data in order to compute the properties of the low-lying states of light and medium-mass nuclei. One of the largest sources of computational uncertainty is then due to the appearance of complex sign oscillations or the 'sign problem' which limits the extent of Euclidean time available for direct PMC calculations. While the sign problem is greatly suppressed by the soft interaction employed in NLEFT, it still represents a significant obstacle to practical PMC calculations, especially in cases, where the number of protons is not equal to the number of neutrons.

In this paper, we focus on the question of the accuracy and reliability of the Euclidean time extrapolations. We start in section 2 by reviewing the PMC formalism, along with the methodology for extrapolating to infinite Euclidean time in section 3. Next, we provide in section 4 an overview of the NLEFT results for the alpha nuclei ranging from ⁴He to ²⁸Si, as obtained from Euclidean time extrapolations of PMC data corresponding to the lattice action described in [12, 13, 15]. In section 5, we consider the effects of statistical and systematical errors on the accuracy of the extrapolation method, and we conclude in section 6 by a discussion of future improvements and refinements to our extrapolation methods.

Our NLEFT calculations are, as in chiral nuclear EFT, organized in powers of a soft scale Q, which is associated with factors of momenta and the pion mass. The contributions of $\mathcal{O}({{Q}^{0}})$ to the nuclear Hamiltonian are referred to as leading order (LO), the $\mathcal{O}({{Q}^{2}})$ terms are of next-to-leading order (NLO), and the $\mathcal{O}({{Q}^{3}})$ terms are of next-to-next-to-leading order (NNLO), at which point our present calculations are truncated. The LO lattice Hamiltonian includes a significant part of the NLO and higher-order contributions, as smeared contact interactions [8, 16, 17] are used. It should also be noted that since we are using a low-momentum power counting scheme, no additional two-nucleon contributions arise at NNLO beyond the terms already appearing at NLO, as these can be absorbed into redefinitions of the NLO couplings. A full discussion of the interactions used in the reported results can be found in [15].

In the present NLEFT calculations, we have used a periodic L = 6 cube and a lattice spacing of a = 1.97 fm, which translates into a cube length of La = 11.82 fm. Our initial wave functions, $|\Psi _{A}^{{\rm init}}\rangle$, are Slater-determinant states composed of delocalized standing waves in the periodic cube with A nucleons. Localized alpha-cluster trial states have also been used for studies of ¹²C and ¹⁶O [12, 13, 18]. These provide not only a consistency check on the Euclidean time extrapolation, but also an opportunity to assess the spatial structure of the nuclei. For simplicity, we describe our calculations using the language of continuous time evolution, even though our AFQMC calculations use transfer matrices with a temporal lattice spacing of $a_{t}^{{}}=1.32$ fm [9]. The Euclidean projection time is given by $t=N_{t}^{{}}a_{t}^{{}}$, where $N_{t}^{{}}$ denotes the number of Euclidean time slices.

We start the Euclidean time projection by means of a 'low-energy filter' based upon Wigner's SU(4) symmetry, where the spin and isospin degrees of freedom of the nucleon are all equivalent as four components of an SU(4) multiplet. The SU(4) symmetric Hamiltonian is of the form

$$\begin{eqnarray}&&{{H}_{{\rm SU}(4)}}\equiv {{H}_{{\rm free}}}+\frac{1}{2}{{C}_{{\rm SU}(4)}}\mathop{\sum }\limits_{\vec{n},\vec{n}^{\prime} }:\rho (\vec{n})f(\vec{n}-\vec{n}^{\prime} )\rho (\vec{n}^{\prime} ):,\end{eqnarray} \tag{ 1 }$$

where $f(\vec{n}-\vec{n}^{\prime} )$ is a Gaussian smearing function with its width set by the average effective range of the two S-wave interaction channels, and ρ is the total nucleon density. Application of the exponential of ${{H}_{{\rm SU}(4)}}$ gives

$$\begin{eqnarray}&&|{{\Psi }_{A}}({{t}^{\prime }})\rangle \equiv {\rm exp} (-{{H}_{{\rm SU}(4)}}{{t}^{\prime }})|\Psi _{A}^{{\rm init}}\rangle ,\end{eqnarray} \tag{ 2 }$$

referred to as a 'trial state'. This part of the calculation is computationally straightforward, as it only requires a single auxiliary field. Most significantly, it does not generate any sign oscillations in the Monte Carlo calculation.

Next, we use the full LO Hamiltonian $H_{{\rm LO}}^{{}}$ to evolve the trial state for a time t, and construct the Euclidean-time projection amplitude

$$\begin{eqnarray}&&{{Z}_{A}}(t)\equiv \langle {{\Psi }_{A}}({{t}^{\prime }})|{\rm exp} (-{{H}_{{\rm LO}}}t)|\Psi _{A}^{{}}({{t}^{\prime }})\rangle ,\end{eqnarray} \tag{ 3 }$$

from which we compute the 'transient energy'

$$\begin{eqnarray}&&{{E}_{A}}(t)=-\partial [{\rm ln} {{Z}_{A}}(t)]/\partial t,\end{eqnarray} \tag{ 4 }$$

by means of a numerical finite difference. Hence, if the lowest eigenstate of ${{H}_{{\rm LO}}}$ that possesses a non-vanishing overlap with the trial state $|{{\Psi }_{A}}({{t}^{\prime }})\rangle$ is denoted $|{{\Psi }_{A,0}}\rangle$, the energy ${{E}_{A,\,0}}$ of $|{{\Psi }_{A,0}}\rangle$ is obtained as the $t\to \infty$ limit of ${{E}_{A}}(t)$. Sign oscillations in the Monte Carlo calculation set the main limitation on the number of Euclidean time steps for which equation (4) can be evaluated. The coupling ${{C}_{{\rm SU}(4)}}$ is a free parameter which can be used either to optimize the convergence of the Euclidean time evolution, or to provide additional constraints for the extrapolation $t\to \infty$.

Higher-order corrections to ${{E}_{A,\,0}}$ are evaluated using perturbation theory. We compute expectation values using

$$\begin{eqnarray}&&Z_{A}^{\mathcal{O}}(t)\equiv \langle \Psi _{A}^{{}}({{t}^{\prime }})|{\rm exp} (-{{H}_{{\rm LO}}}t/2)\mathcal{O}{\rm exp} (-{{H}_{{\rm LO}}}t/2)|{{\Psi }_{A}}({{t}^{\prime }})\rangle ,\end{eqnarray} \tag{ 5 }$$

for any operator $\mathcal{O}$. Given the ratio

$$\begin{eqnarray}&&X_{A}^{\mathcal{O}}(t)=Z_{A}^{\mathcal{O}}(t)/{{Z}_{A}}(t),\end{eqnarray} \tag{ 6 }$$

the expectation value of $\mathcal{O}$ for the desired state $|{{\Psi }_{A,0}}\rangle$ is again obtained in the $t\to \infty$ limit according to

$$\begin{eqnarray}&&X_{A,0}^{\mathcal{O}}\equiv \langle {{\Psi }_{A,0}}|\mathcal{O}|{{\Psi }_{A,0}}\rangle =\mathop{{\rm lim} }\limits_{t\to \infty }\;X_{A}^{\mathcal{O}}(t),\end{eqnarray} \tag{ 7 }$$

which gives the corrections to ${{E}_{A,\,0}}$ induced by the NLO and NNLO contributions, including the effects of strong and electromagnetic isospin symmetry breaking.

The closer the trial state $|\Psi _{A}^{{}}({{t}^{\prime }})\rangle$ is to $|{{\Psi }_{A,0}}\rangle$, the less the required projection time t. The trial state can be optimized by adjusting both the SU(4) projection time ${{t}^{\prime }}$ and the strength of the coupling ${{C}_{{\rm SU}(4)}}$ of ${{H}_{{\rm SU}(4)}}$. As shown in section 4, the accuracy and reliability of the extrapolation $t\to \infty$ is greatly improved by simultaneously incorporating data from multiple trial states that differ in the choice of ${{C}_{{\rm SU}(4)}}$. This approach enables a 'triangulation' of the asymptotic behavior as the common limit of several different functions of t.

In most cases, reaching the limit $t\to \infty$ requires an extrapolation from finite values of t. Given the limited extent of the data, there are uncertainties in this extrapolation. The behavior of $Z_{A}^{{}}(t)$ and $Z_{A}^{\mathcal{O}}(t)$ at large t is controlled by the low-energy spectrum of ${{H}_{{\rm LO}}}$. Let $|E\rangle$ label the eigenstates of ${{H}_{{\rm LO}}}$ with energy E, and let ${{\rho }_{A}}(E)$ denote the density of states for a system of A nucleons. For simplicity, we omit additional labels needed to distinguish degenerate states. We can then express $Z_{A}^{{}}(t)$ and $Z_{A}^{\mathcal{O}}(t)$ in terms of their spectral representations

$$\begin{eqnarray}&&{{Z}_{A}}(t)=\int {\rm d}E{{\rho }_{A}}(E){{\left| \langle E|{{\Psi }_{A}}({{t}^{\prime }})\rangle \right|}^{2}}{\rm exp} (-Et),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{rcl} Z_{A}^{\mathcal{O}}(t) & = & \int {\rm d}E{\rm d}{{E}^{\prime }}{{\rho }_{A}}(E){{\rho }_{A}}({{E}^{\prime }}){\rm exp} (-(E+{{E}^{\prime }})t/2) \\ {} & {} & \times \langle {{\Psi }_{A}}({{t}^{\prime }})|E\rangle \langle E|\mathcal{O}|{{E}^{\prime }}\rangle \langle {{E}^{\prime }}|\Psi _{A}^{{}}({{t}^{\prime }})\rangle . \\ \end{array}\end{eqnarray} \tag{ 9 }$$

The spectral representations of ${{E}_{A}}(t)$ and $X_{A}^{\mathcal{O}}(t)$ are then obtained by using equations (4) and (6), respectively. We can approximate these to arbitrary accuracy over any finite range of t by taking ${{\rho }_{A}}(E)$ to be a sum of energy delta functions

$$\begin{eqnarray}&&{{\rho }_{A}}(E)\approx \mathop{\sum }\limits_{k=0}^{{{k}_{{\rm max} }}}{{c}_{A,k}}\delta (E-{{E}_{A,k}}).\end{eqnarray} \tag{ 10 }$$

Let us now consider the extrapolation formulas for the LO energy and the higher order perturbative corrections. From equations (8) and (10), we find for the LO energy

$$\begin{eqnarray}&&E_{A}^{j}({{N}_{t}})={{E}_{A,0}}+\mathop{\sum }\limits_{k=1}^{{{k}_{{\rm max} }}}|{{c}_{A,j,k}}|{\rm exp} \left( -\frac{{{\Delta }_{A,k}}{{N}_{t}}}{{{\Lambda }_{t}}} \right),\end{eqnarray} \tag{ 11 }$$

where $t={{N}_{t}}/{{\Lambda }_{t}}$ with ${{\Lambda }_{t}}=150$ MeV, corresponding to ${{a}_{t}}=1.32$ fm. The energy gaps are defined as ${{\Delta }_{A,k}}\equiv {{E}_{A,k}}-{{E}_{A,0}}$, and the index j denotes a specific choice of ${{t}^{\prime }}$ and ${{C}_{{\rm SU}(4)}}$ in the trial wave function $|{{\Psi }_{A}}({{t}^{\prime }})\rangle$. We take $\Delta _{A}^{k+1}\gt \Delta _{A}^{k}$, and ${{k}_{{\rm max} }}=3$ for ⁴He (A = 4) and ${{k}_{{\rm max} }}=2$ for $A\geqslant 8$. For the operator matrix elements that make up the perturbative NLO and NNLO corrections, we find

$$\begin{eqnarray}&&X_{A}^{\mathcal{O},j}({{N}_{t}})=X_{A,0}^{\mathcal{O}}+\mathop{\sum }\limits_{k=1}^{{{k}_{{\rm max} }}}{{x}_{A,j,k}}{\rm exp} (-\frac{{{\Delta }_{A,k}}{{N}_{t}}}{2{{\Lambda }_{t}}}),\end{eqnarray} \tag{ 12 }$$

where the dominant contributions are taken to be due to transition amplitudes involving the ground state and excited states. In order for this to be a good approximation, it is necessary that the overlap between our trial state and the ground state not be small compared to the overlap with the low-lying excited states. It should be noted that the coefficients ${{x}_{A,j,k}}$ can be positive as well as negative, which gives us the possibility of 'triangulating' the asymptotic values $X_{A,0}^{\mathcal{O}}$ from above and below. For this purpose, the parameters ${{t}^{\prime }}=N_{t}^{\prime }/{{\Lambda }_{t}}$ and ${{C}_{{\rm SU}(4)}}$ should be optimally chosen for each value of A. For equation (11), the dependence on t is monotonically decreasing and thus no triangulation of the asymptotic value from above and below is possible. However, we are helped by the fact that the rate of convergence is twice that of equation (12). In order to determine ${{E}_{A,\,0}}$ and $X_{A,0}^{\mathcal{O}}$, a correlated χ² fit to the LO energy and all NLO and NNLO matrix elements in NLEFT is performed for each value of A. This procedure also determines the coefficients ${{c}_{A,j,k}},{{x}_{A,j,k}}$ and ${{\Delta }_{A,k}}$. We find that using 2–6 distinct trial states for each A allows for a significantly more accurate and stable determination of ${{E}_{A,\,0}}$ and $X_{A,0}^{\mathcal{O}}$ than would be possible with a single trial state. Note that the energy gaps ${{\Delta }_{A,k}}$ in the extrapolation functions are taken to be independent of the trial wave function j, which gives an additional consistency criterion. We find that a simultaneous description using equations (11) and (12) accounts for all of the PMC data we have obtained for different $|{{\Psi }_{A}}({{t}^{\prime }})\rangle$.

We shall next elaborate on how the extrapolation methods of section 3 perform when confronted with actual PMC data. First, we show our data for the light nuclei ⁴He and ⁸Be in figures 1 and 2, respectively. Further, our results for ¹²C are given in figure 3, for ¹⁶O in figure 4, for ²⁰Ne in figure 5, for ²⁴Mg in figure 6, and for ²⁸Si in figure 7. The curves show a correlated χ² fit for all trial states with a given A, using the same spectral density ${{\rho }_{A}}(E)$. The upper row in each figure shows the LO energy, the total isospin-symmetric 2NF correction (NLO), the electromagnetic and isospin-breaking corrections and the total 3NF correction. The remaining panels show the matrix elements $X_{A}^{\mathcal{O}}(t)$ that form part of the NLO and 3NF terms.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the second and third rows of figure 1 through figure 7, the operators $\partial {{E}_{A}}/\partial {{C}_{i}}$ give the contributions of the NLO contact interactions. The interactions that involve ${{C}_{0}}$, ${{C}_{{{I}^{2}}}}$, ${{C}_{{{Q}^{2}}}}$, ${{C}_{{{I}^{2}}{{Q}^{2}}}},$ ${{C}_{{{S}^{2}}{{Q}^{2}}}}$, ${{C}_{{{S}^{2}}{{I}^{2}}{{Q}^{2}}}}$, ${{C}_{{{(QS)}^{2}}}}$ and ${{C}_{{{I}^{2}}{{(QS)}^{2}}}}$ are defined in equations (18)–(23) of [17]. Similarly, the interactions involving ${{C}_{{{S}^{2}}{{Q}^{2}}{\rm SO}(3){\rm B}}}$ and ${{C}_{{{I}^{2}}{{S}^{2}}{{Q}^{2}}{\rm SO}(3){\rm B}}}$ are given in equations (55) and (56) of [17], and $\Delta {{E}_{A}}(\Delta {{x}_{\pi }})$ denotes the energy shift due the $\mathcal{O}({{a}^{2}})$-improved pion–nucleon coupling in equation (57) of [17]. The operators $\partial {{E}_{A}}/\partial {{D}_{i}}$ give the individual contributions to the 3NF correction, which are defined in equations (37)–(41) of [17].

The error estimates given in parentheses in the captions of figure 1 through figure 7 have been obtained from a χ² minimization using the PMC errors for each datapoint as weight factors. This procedure could be affected by autocorrelations in Monte Carlo time, and secondly by the fact that the observables (except for the LO energy) are formed out of a ratio of amplitudes for which the PMC error is known separately, according to equation (6). Instead of simply adding these PMC errors in quadrature, a more realistic error estimate could be obtained by means of a resampling algorithm such as the Jackknife or the Bootstrap (for a pedagogical introduction, see e.g. [20]). In order to decrease the effects of any residual autocorrelations in the PMC data, the Jackknife method can be combined with 'blocking' of the data, whereby adjacent (in Monte Carlo time) samples are combined into blocks of increasing size until the variance of the sample converges as a function of the block size (see e.g. [21]).

In our PMC production runs, we do not use ${{k}_{{\rm max} }}=1$ (corresponding to a single energy gap) as this would in most cases lead to large values of χ² as well as inconsistent results for different choices of trial states (see table 1). Extrapolations with ${{k}_{{\rm max} }}=2$ account very well for the medium-mass nuclei, which appear to be highly compact objects for which the contamination from low-lying excited states is small. For ⁴He, we find that ${{k}_{{\rm max} }}=3$ is required in order to account for all the PMC data for all observables and all trial states. However, while extrapolations with ${{k}_{{\rm max} }}=3$ may provide a more accurate description over a larger range of trial states, such fits are also much more difficult to constrain adequately, due to the much larger number of adjustable parameters involved. The extent and resolution of our PMC data in Euclidean time also limits, in most cases, the number of resolvable energy gaps to ${{k}_{{\rm max} }}=2$.

Table 1.  Uncertainty analysis of the Euclidean time extrapolation for ¹²C with ${{k}_{{\rm max} }}=2.$ The values of ${{C}_{{\rm SU}(4)}}$ (in MeV⁻²) for each trial state shown in figure 3 are '1'  =  $-8.0\times {{10}^{5}}$, '2'  =  $-1.2\times {{10}^{4}}$, '9'  =  $-7.5\times {{10}^{5}}$, '10'  =  $-6.5\times {{10}^{5}}$, and '11'  =  $-7.0\times {{10}^{5}}$. The quantities shown (in MeV) are: the LO non-perturbative 2NF result, followed by the perturbative higher-order and isospin-breaking corrections as described in the main text. The fit labeled '5' (shown in figure 3) is a correlated extrapolation using all trial states. The fits labeled '4a'–'4e' check the consistency of fit '5' under the removal of a single trial state from the full analysis. We also show the (poorly constrained) extrapolations '1a'–'1e', where each trial state is treated separately. Note that this only allows for an analysis with ${{k}_{{\rm max} }}=1$. The one-standard deviation error estimates (given in parentheses) are obtained from a Marquardt–Levenberg minimization of χ² (per d.o.f.) with the Monte Carlo error estimates used as weights.

Fit	Trial states	LO (2NF)	NLO (2NF)	EMIB (2NF)	NNLO (3NF)	χ2
5	$1,2,9,10,11$	$-96.9(2)$	$10.48(3)$	$7.76(1)$	$-14.80(6)$	0.68
$4a$	$1,2,9,10$	$-96.8(2)$	$10.46(4)$	$7.76(1)$	$-14.84(6)$	0.71
$4b$	$1,2,9,11$	$-97.0(2)$	$10.45(3)$	$7.76(1)$	$-14.85(6)$	0.72
$4c$	$1,2,10,11$	$-96.8(2)$	$10.46(4)$	$7.76(2)$	$-14.80(6)$	0.72
$4d$	$1,9,10,11$	$-97.3(2)$	$10.54(5)$	$7.72(2)$	$-14.61(6)$	0.57
$4e$	$2,9,10,11$	$-96.9(2)$	$10.44(4)$	$7.75(2)$	$-14.94(7)$	0.69
$1a$	1	$-95.0(2)$	$10.17(2)$	$7.79(1)$	$-13.93(4)$	1.83
$1b$	2	$-94.4(2)$	$10.55(2)$	$7.98(1)$	$-14.46(5)$	3.35
$1c$	9	$-93.5(2)$	$10.03(2)$	$7.72(1)$	$-13.60(4)$	1.14
$1d$	10	$-94.0(9)$	$9.10(9)$	$7.32(3)$	$-13.73(16)$	2.36
$1e$	11	$-92.1(2)$	$9.83(2)$	$7.58(1)$	$-13.34(4)$	1.28

Our extrapolation procedure is examined in detail for the case of ¹²C in table 1. Our main fit, labeled 'fit 5', is also shown in figure 3, and consists of a simultaneous fit to five trial states that differ in the value of $C_{{\rm SU}(4)}^{{}}$. We observe that ${{\chi }^{2}}\simeq 0.68$, which indicates that the error bars of the individual MC data points are likely to be overestimated. This could plausibly happen as the observables are formed from the ratio of two amplitudes, the errors of which are at present simply added in quadrature. A full-fledged jackknife error analysis may yield a more realistic result. The uncertainties shown in parentheses correspond to the variances reported by the χ² minimization. In order to assess the accuracy of these error estimates as well as the stability of the central values, we have generated a number of fits, where each one of the five trial states has been excluded in turn, in the spirit of the Jackknife method [20]. These fits are denoted '4a' through '4e' in table 1. Evidently, these agree closely with the full analysis, the largest discrepancy being due to the exclusion of trial state '2', which generates the largest shift in the extrapolated values as well as a significant reduction in χ². As is evident from figure 3, trial state '2' (denoted by blue open squares) is furthest away from the 'triangulation point', and may therefore not be completely described by an extrapolation with ${{k}_{{\rm max} }}=2$. An even more stable result might be obtained by replacing that trial state in the analysis with one which is closer to the triangulation point.

For comparison, we also show in table 1 the results of independent, uncorrelated fits to each of the five trial states for ¹²C. These extrapolations do not benefit from the consistency requirements of the multi-trial state extrapolations, and furthermore these can only be taken to ${{k}_{{\rm max} }}=1$, as the extent of the data in Euclidean time is too short to constrain more than one energy gap. We observe that such extrapolations are clearly much less reliable, and suffer from several pronounced issues. One is the clear tendency for 'spurious early convergence', which is due to the lack of enforced independence on the value of ${{C}_{{\rm SU}(4)}}$. We also observe that the extrapolated values as well as the χ² fluctuate significantly between extrapolations to different trial states. Again, trial state '2' appears to be the most problematic, although we also find that fits with a smaller χ² are no more reliable than those with a larger value. Neither do the variances produced by the χ² minimization properly describe the uncertainties.

Also, as shown in table 1, the total uncertainty is dominated by that of the LO contribution, which does not consist of a ratio of amplitudes according to equation (6). This relatively large extrapolation error is due to the appearance of the absolute values of the coefficients ${{c}_{A,j,k}}$ in equation (11), which prevents a triangulation of the LO contribution. Nevertheless, the Jackknife method could be used to evaluate the effect of autocorrelations between consecutive auxiliary-field configurations in Monte Carlo time on the LO result. We find that the elimination of such autocorrelations from the start, by allowing for sufficient decorrelation time between consecutive measurements, does not present any difficulties. For the NLO and NNLO operators, more consistent and reliable error estimates for the individual data points could clearly be obtained by Jackknife resampling of the Monte Carlo data. While such an analysis would not alter any of the conclusions concerning the stability and consistency of the Euclidean time extrapolations, we have investigated the likely outcome of a full Jackknife resampling of our data by performing multiple independent Monte Carlo runs for ¹²C with ${{N}_{t}}=12$. For such runs, we obtain LO energies for ${{N}_{t}}=11.0$, 11.5 and 12.5 by means of a numerical finite difference. Similarly, for the higher-order corrections we compute matrix elements for ${{N}_{t}}=11.5$ and 12.5. As expected, we find that the uncertainties of the LO energies are accurately given by the Monte Carlo errors, whereas those of the individual NLO matrix elements appear overestimated by a factor of $\simeq 2$ due to cancellations between the numerator and denominator of equation (6). This result is consistent with the values of ${{\chi }^{2}}\lt 1$ obtained from the full extrapolation.

As the extent of our PMC data in Euclidean time is relatively short, we discuss next the expected reliability of our extrapolated results. Fortunately, the number of Euclidean time steps ${{N}_{t}}$ available for the extrapolation does not decrease drastically with the number of nucleons A. At this time, our method has been successfully applied to the spectrum, structure and electromagnetic properties of ¹²C in [13, 14], and to those of ¹⁶O in [18], where consistency between delocalized plane-wave and alpha-cluster trial wave functions was established. In our 'triangulation' method, the extrapolation is very strongly constrained by the requirement that all observables, for all trial states, should be described by the same exponential dependence on the Euclidean projection time t. Rapid convergence in t then translates into a small sensitivity to ${{C}_{{\rm SU}(4)}}$ at large values of t, which helps to guard against 'spurious early convergence', where a smaller energy gap is overlooked.

Nevertheless, in the absence of consistency conditions on the extrapolations, we find that our method is accurate even for ${{k}_{{\rm max} }}=1$ when the leading energy gap is very large, such as for the deuteron in a periodic L = 3 cube with lattice spacing a = 1.97 fm (see [14]). In that case, the extrapolated results can be directly compared with Lanczos diagonalization, as shown in table 2 and figure 8. In spite of this impressive agreement, we still need to consider the possibility that a sufficiently small energy gap can be missed in the extrapolation due to the limited extent and resolution of the PMC data in Euclidean time. Of all the results presented here, those for ⁸Be are likely to be most affected by the limited Euclidean projection time, as the convergence to the ground state is clearly the slowest and the data do not allow for extrapolations beyond ${{k}_{{\rm max} }}=2$, unlike the case of ⁴He, where up to three energy gaps could be constrained by the PMC data.

Zoom In Zoom Out Reset image size

Table 2.  Euclidean time extrapolation for the deuteron with ${{k}_{{\rm max} }}=1$ in a periodic L = 3 cube with lattice spacing a = 1.97 fm, with ${{\tilde{g}}_{\pi N}}\equiv {{g}_{A}}/(2{{f}_{\pi }})$ and ${{m}_{{\rm N}}}$ the nucleon mass. ${{E}_{d}}({\rm LO})$ denotes the (non-perturbative) LO energy, and the remaining quantities are perturbative contributions which quantify the sensitivity of the of ${{E}_{d}}$ to small shifts in the pion mass (for full details and definitions, see [14]). The appropriate units are given for each quantity, with '[l.u.]' indicating units of the inverse (spatial) lattice spacing. The second column shows the extrapolated Monte Carlo results, with one-standard-deviation errors similar to those in table 1 given in parentheses. The third column shows the results obtained from a Lanczos diagonalization of the two-nucleon Hamiltonian.

Observable	2H (PMC+ex)	2H (Lanczos)
${{E}_{d}}({\rm LO})$ (MeV)	−9.070(12)	−9.078
$\Delta {{E}_{d}}(\Delta {{\tilde{M}}_{\pi }})$ (MeV)	−0.003 548(12)	−0.003 569
$\Delta E_{d}^{{}}(\Delta \tilde{M}_{\pi }^{{\rm IB}})$ (MeV)	−0.002 372(8)	−0.002 379
$\partial {{E}_{d}}/\partial {{m}_{{\rm N}}}$	−0.003 82(2)	−0.003 809
$\partial {{E}_{d}}/\partial {{\tilde{g}}_{\pi N}}$ (l.u.)	0.010 24(11)	0.010 17
$\partial {{E}_{d}}/\partial {{C}_{0}}$ (l.u.)	0.138 97(15)	0.138 867
$\partial {{E}_{d}}/\partial {{C}_{I}}$ (l.u.)	−0.4171(4)	−0.416 60

We have presented an overview of the techniques and the analysis used for Euclidean time extrapolations in NLEFT. The core issue is that, due to computational constraints, one must extract asymptotic values from a limited region in Euclidean projection time. In order to reduce the uncertainty of the extrapolation, we perform multi-exponential fits as prescribed by spectral decomposition for the asymptotic behavior of the projection amplitudes. We find that the fitted asymptotic values are greatly stabilized by using multiple initial states and observables. As examples of this analysis, we have shown energies at LO, NLO, and NNLO, as well as individual operator expectation values for the alpha nuclei up to A = 28.

While these results are promising, the methods we have presented leave much room for further improvement. In particular, the current analysis is not adequate for the case when our trial states have only a small overlap with the ground state. In such cases, the raw PMC data will be far from their asymptotic values and it is unlikely that combining several different sets of data will provide any significant improvement. Another unfavorable situation that can arise is when the sign oscillations are severe and the quality of the PMC data degrades very quickly with projection time.

In these more difficult scenarios, one can improve on the situation by performing coupled multi-channel projections that evolve more than one initial state simultaneously in Euclidean time. This is the same approach which is being used in a technique called the 'adiabatic projection method' used for describing scattering states of clusters [22, 23]. As shown in [23], when projecting n states the exponential convergence of the ground state energy in projection time is given by the excitation energy of the nth excited state rather than by that of the first excited state. This has the potential to help significantly in circumventing the limited extent of the Euclidean time propagation.

We are grateful for assistance in automated data collection by Thomas Luu. We acknowledge partial financial support from the Deutsche Forschungsgemeinschaft (Sino-German CRC 110), the Helmholtz Association (Contract No. VH-VI-417), BMBF (Grant No. 05P12PDFTE), the US Department of Energy (DE-FG02-03ER41260), and the US National Science Foundation (PHY-1307453). Further support was provided by the EU HadronPhysics3 project and the ERC Project No. 259218 NUCLEAREFT. The computational resources were provided by the Jülich Supercomputing Centre at Forschungszentrum Jülich and by RWTH Aachen. TL acknowledges a grant from the Magnus Ehrnrooth foundation of the Finnish Society of Sciences and Letters.