Systematics of the HAL QCD Potential at Low Energies in Lattice QCD

The $\Xi\Xi$ interaction in the $^1$S$_0$ channel is studied to examine the convergence of the derivative expansion of the non-local HAL QCD potential at the next-to-next-to-leading order (N$^2$LO). We find that (i) the leading order potential from the N$^2$LO analysis gives the scattering phase shifts accurately at low energies, (ii) the full N$^2$LO potential gives only small correction to the phase shifts even at higher energies below the inelastic threshold, and (iii) the potential determined from the wall quark source at the leading order analysis agrees with the one at the N$^2$LO analysis except at short distances, and thus, it gives correct phase shifts at low energies. We also study the possible systematic uncertainties in the HAL QCD potential such as the inelastic state contaminations and the finite volume artifact for the potential and find that they are well under control for this particular system.


I. INTRODUCTION
In lattice QCD, two methods have been proposed so far to study the baryon-baryon interactions. One is the direct method [1][2][3], where the energy spectrum on finite volume(s) is extracted from the temporal correlation of two baryons and is converted to the scattering phase shift and/or the binding energy in the infinite volume through the Lüscher's finite volume formula [4,5]. The other is the HAL QCD method [6][7][8][9][10], where the potential between baryons is first derived from the spatial correlations of two baryons, and it is used to calculate the observables through the Schrödinger-type equation in the infinite volume.
While both methods are supposed to give the same results in principle, previous numerical studies for two-nucleon (N N ) systems show clear discrepancy: The direct method indicates that both dineutron ( 1 S 0 ) and deuteron ( 3 S 1 ) are bound for heavy pion masses (m π ≥ 300 MeV), while the HAL QCD method does not provide such bound states in both channels for heavy pion masses. This discrepancy was recently discussed in a series of papers [11][12][13][14], where it was pointed out that the effective two-particle energy as a function of the Euclidean time may significantly suffer from elastic scattering states of two nucleons. To elucidate such uncertainties, certain "normality checks" for the finite-volume spectrums were introduced [11][12][13][14].
The advantage of the time-dependent HAL QCD method [8] over the direct method is that the former is free from the ground state saturation problem in principle, since the energy-independent potential controls both ground state and the elastic excited states simultaneously as long as the inelastic scatterings in the small Euclidean time are properly suppressed. 1 In practice, there appear systematic uncertainties associated with the truncation of the derivative expansion for the non-local potential. Therefore, the main purpose of the present paper is to study the convergence of the derivative expansion, as well as other sources of systematic uncertainties such as the inelastic state contaminations and the distortion of the interaction under finite volume. We consider the ΞΞ system in the 1 S 0 channel and perform the (2+1)-flavor lattice QCD calculation at m π = 0.51 GeV and m K = 0.62 GeV.
Because of the large quark masses, the statistical errors in this case become relatively small, so that one can focus on the detailed analysis of the systematic errors. Also, this channel and the N N system in the 1 S 0 channel belong to the same multiplet in the flavor SU (3) limit.
This paper is organized as follows. In Sec. II, we review the time-dependent HAL QCD method. In Sec. III, we present the lattice QCD results for the ΞΞ interaction in the 1 S 0 channel at the next-to-next-to-leading order (N 2 LO) in the derivative expansion. The N 2 LO potential is extracted from a specific combination of the ΞΞ correlations with different source operators. The systematic errors associated with the inelastic state contaminations and the distortion in the finite volume are also examined. In Sec. IV, we calculate the scattering phase shifts in this channel, and check the convergence of the derivative expansion in the HAL QCD method. In Sec. V, we demonstrate the self-consistency between the phase shifts obtained from the HAL QCD potential and those obtained from the energy spectra obtained from the HAL QCD potential combined with the Lüscher's formula. Sec. VI is devoted to the conclusion. In Appendix A, we discuss the relation between the energy-independent non-local potential and the energy-dependent local one.

II. FORMALISM
The key quantity in the HAL QCD method [6][7][8][9][10] is the Nambu-Bethe-Salpeter (NBS) wave function, defined by where |0 is the vacuum state of QCD, |2B, W is the QCD eigenstate for two baryons with eigenenergy W , and B( x, t) is a single baryon operator with spin indices omitted for simplicity. We then define a non-local and energy-independent potential U ( r, r ) so as to below inelastic threshold, W < W th = 2m B + m π , with m B the baryon mass, m π the pion mass, and W = 2 m 2 B + k 2 . Here we define E k = k 2 /(2µ) and H 0 = −∇ 2 /(2µ) with a reduced mass µ = m B /2. We note that U ( r, r ) depends on the specific choice of the interpolating operator B(x) used in Eq. (1). Nevertheless, the S-matrix is free from the choice of B(x) as long as it is an "almost-local operator field" [15] (Nishijima-Zimmermann-Haag theorem).
To extract the NBS wave function in lattice QCD, we start with the two-baryon correlation function, where J 2B (t 0 ) is a source operator for two-baryon. By inserting the complete set, we obtain where W n = 2 m 2 B + k 2 n is the n-th energy eigenvalue, A n ≡ 2B, W n |J 2B (0)|0 corresponds to the overlap with each elastic eigenstate, and the ellipses represent the inelastic contributions. In principle, one can extract A 0 ψ W 0 ( r) for the lowest energy W 0 from the large t behavior of C 2B ( r, t).
In practice, however, since C 2B ( r, t) becomes too noisy at large t, we need to employ the time-dependent HAL QCD method [8]. Let us define the ratio of correlation functions, which we call the R-correlator, as with ∆W n = W n − 2m B , ∆W th = W th − 2m B and A n = A n /C 2 , where C B (t) and C are a single baryon correlation function and the corresponding overlap factor, respectively. They are given by where J B (t 0 ) is a single baryon source operator and ellipses represent the inelastic states contributions.
Since the non-local potential U ( r, r ) is defined to be energy-independent [7], all elastic scattering states below the threshold share the same U ( r, r ). Therefore, Eq. (2) with an where the effect of the inelastic channel of O(e −∆W th t ) is neglected in the right hand side, while there is no term beyond ∂ 2 /∂t 2 in the left hand side of Eq. (7), i.e., Eq. (7) is derived without non-relativistic approximation.
Note that the ground state saturation is no more required in this time-dependent HAL QCD method. Instead, the required condition is that R( r, t) is saturated by the contributions from elastic states ("the elastic state saturation"), which can be achieved by a moderate value with m NG being the mass of the lightest Nambu-Goldstone boson). 2 This is the fundamental difference between the HAL QCD method and the direct method.
As discussed in [9,10], U ( r, r ) in Eq. (7) is not determined uniquely by R( r, t), though different U ( r, r )s give same observables below the inelastic threshold. In the HAL QCD method, the derivative expansion scheme enables one to extract one of the possible U ( r, r )s in a unique manner. Let us consider the two-baryon system in the spin-singlet channel.
Then the leading order (LO) analysis neglecting the higher orders leads to with In order to examine the convergence of the derivative expansion, we consider the N 2 LO analysis in this paper, The relation between the potential from the LO analysis, V LO 0 (r), and those from the N 2 LO analysis, V N 2 LO 0 (r) and V N 2 LO 2 (r) is given by which shows that the N 2 LO correction in V LO 0 (r) depends on both V N 2 LO 2 (r) and the spatial profile of the R-correlator, the latter of which depends not only on the spatial profile of the NBS wave functions ψ Wn (r) but also on their magnitude A n in the R-correlator. The (r) are t-independent as long as the elastic state saturation is achieved and the higher order contributions in the derivative expansion can be neglected. One may also estimate the magnitude of systematic errors from the truncation of the derivative expansion and from the inelastic state contaminations by studying the t-dependence of the potentials. 2 There is a possibility that the inelastic contributions cancel partially between the numerator and the denominator of R( r, t), so that the elastic state saturation in R( r, t) may appear for smaller t than those in C 2B ( r, t) and C B (t).

A. Lattice Setup
Throughout this paper, we use 2+1 flavor QCD ensembles [16], generated by using the We employ the wall source q wall (t) = y q( y, t), which has been mainly used in the previous studies by the HAL QCD method, and the smeared source q smear ( [16]. For the smeared source, the same x is taken as the center of the smeared source for all six quarks in two baryons as has been done in Ref. [16]. For both sources, the point-sink operator for each baryon ("point-sink scheme" in the HAL QCD method [17]) is exclusively employed in this study. The correlation functions are calculated by the unified contraction algorithm (UCA) [18]. A number of configurations and other parameters are summarized in Table I. Statistical errors are evaluated by the jack-knife method. For more details on the simulation setup, see Ref. [11].
In the present study, we focus on the ΞΞ system in the 1 S 0 channel: This is one of the most convenient choices to obtain the insights of N N systems, since it belongs to the same 27 representation as the N N system in the 1 S 0 channel in the flavor SU(3) limit but has much better signal to noise ratio than the N N ( 1 S 0 ) case. We use the relativistic interpolating operators [11] for Ξ, which are given by where C = γ 4 γ 2 is the charge conjugation matrix, α and (a, b, c) are the indices for the spinor and color, respectively.

B. The R-correlator
We first consider the behaviors of the R-correlator defined in Eq. (5). Shown in Fig. 1  smeared source (Right). The results show strong quark-source dependence: The R-correlator from the wall source (R wall ( r, t)) is delocalized with a weak t-dependence, while that from the smeared source (R smear ( r, t)) is localized and has a strong t-dependence. If the R-correlator is saturated by the ground state, its spatial profile should be independent of the source and its temporal profile should be simply dictated by an overall factor, exp(−∆W n=0 t).
To see more closely the t-dependence of the spatial profile of the R-correlator, we plot R( r, t) normalized to be unity at r = 3.5 fm for the wall source and at r = 1.0 fm for the smeared source in Fig. 2. The shape of the R-correlator from the wall source has a weak t-dependence, which indicates that the contribution from the elastic scattering states other than the ground state in R wall ( r, t) are relatively small. On the other hand, the shape of R smear ( r, t) show a sizable t-dependence, which indicates that it has a substantial admixture from the several elastic scattering states. Although the parameters (A, B) of the smeared source shown in Table I are tuned to suppress the excited states of a single baryon, the same parameters are not guaranteed to suppress the elastic scattering states for two baryons. Indeed, one of the most relevant parameters which control the magnitudes of elastic state contributions is the relative distance r between two baryons at the source, as The smeared source operator in all previous works in the direct method (except for [3]) essentially corresponds to r = 0 and could be coupled to all elastic scattering states with almost an equal magnitude. 3 See Ref. [14] for more detailed studies on this point.

C. HAL QCD potential at the leading order
Let us now study the potential in the HAL QCD method at the leading order, V LO 0 (r). Fig. 3 shows the one for ΞΞ( 1 S 0 ) and its breakups (H 0 , ∂/∂t and ∂ 2 /∂t 2 terms in Eq. (9)) on L = 64 at t = 13 from the wall source (Left) and the smeared source (Right). For the wall source, the H 0 term is dominant with sizable contributions from the ∂/∂t term, while the ∂ 2 /∂t 2 term is negligible. The ∂/∂t term is not constant as a function of r, which indicates that there exist small but non-negligible contributions from the excited states in R wall ( r, t).
For the smeared source, on the other hand, all terms are important. In particular, the ∂/∂t term (green triangles) shows substantial r-dependence indicating large contributions from the excited states in the smeared source. However, such dependence is cancelled by the H 0 term (blue squares) and is further corrected by the ∂ 2 /∂t 2 term (black diamonds). The final results (red circles) with the smeared source and the wall source show qualitatively similar behaviors, i.e., the repulsive core at the short distance and the attractive pocket at the intermediate distance. This illustrates that the time-dependent HAL QCD method works well for extracting the ΞΞ potential irrespective of the source structures. Fig. 4 is a comparison among the LO potentials (V LO 0 (r)) for different t in each source. For the wall source, the potentials at t = 10−16 are consistent with each other within statistical errors, while those from the smeared source show the detectable t-dependence. Fig. 5 is a comparison of V LO 0 (r) between two sources at t = 10, 12, 14, 16. As t increases, the LO potential from the smeared source gradually converges to that from the wall source. The relatively large t-dependence of the potentials from the smeared source as well as the remaining small discrepancy of potentials between two sources even at t = 16

Shown in
indicate that the N 2 LO analysis in the derivative expansion is necessary to understand the data from the smeared source. This is a natural consequence of the fact that the N 2 LO (11), is much more significant in the smeared source than the wall source as shown in Fig. 3.
D. HAL QCD potential at the next-to-next-to-leading order We next apply the N 2 LO analysis in the derivative expansion to R-correlators for both sources. The potential at the LO analysis, V LO 0 (r), and those at the N 2 LO analysis, (r), satisfy the linear equations given by where source = wall or smear.
To extract V N 2 LO 0,2 (r), we first consider the following relation derived from Eq. (14), (r) can be determined from Eq. (11). (r) by m 2 π to make its mass dimension +1 for a comparison to V 0 (r)'s. We find that V N 2 LO 0 (r) agrees well with the V LO(wall) 0 (r) except at short distances. We also find that V N 2 LO 2 (r) is localized within the range of 1 fm, which is much shorter than the range of V LO(wall),N 2 LO 0 (r). We note here that the negative sign of V N 2 LO 2 (r) does not necessarily imply attraction, since the N 2 LO potential is given by As already mentioned, ∇ 2 R( r, t)/R( r, t) from the smeared source is much larger than that of the wall source (see Fig. 3). Intuitively, this is because R smear (r, t) (R wall (r, t)) contains larger (smaller) contributions from excited states and thus is more (less) sensitive to higher order terms in the derivative expansion of the potential. Therefore, the N 2 LO analysis is mandatory for the smeared source, while the LO analysis for the wall source leads to the potential which is almost identical to V N 2 LO 0 (r). Fig. 7 is the t-dependence of V N 2 LO 0,2 (r) in the range of t = 13 − 16. Since appreciable t-dependence is not seen within the error bars, the N 4 LO contribution is expected to be small.   is stable for t much less than 16, suggesting that the systematic error originating from the inelastic contributions of the single-baryon cancels largely between the numerator and the denominator of the R-correlator for the wall source.

F. Effect of the finite volume
In Fig. 9, we show the volume dependence of the potential at the LO analysis for the wall source at t = 13 with L = 40, 48 and 64. All the potentials are consistent with each other within statistical errors. This indicates that the artifact due to finite volume is negligible for the potential, mainly because the potential is short ranged.

IV. SCATTERING PHASE SHIFTS
In the previous section, we examine systematic uncertainties on the HAL QCD potential.
In this section, we examine how these systematic uncertainties affect the physical observables such as the scattering phase shifts, in particular the effect of the derivative expansion. To calculate the scattering phase shifts, δ 0 (k), we first fit the potentials by a sum of Gaussians,  Table II. In Fig. 10, we show the comparison of the scattering phase shifts from V LO(wall) 0 (r), been also observed for the N N systems in the 1 S 0 and 3 S 1 channels in quenched QCD with m π 530 MeV [21] and the I = 2 ππ system in (2+1)-flavor QCD with m π 870 MeV [17].
The scattering length a 0 obtained through lim k→0 k cot δ 0 (k) = 1/a 0 from V LO(wall) 0 (r), (r)∇ 2 at t = 13−16 is shown in Fig. 11. The result indicates that the scattering length is almost insensitive to the degrees of the approximation but has a small variation in t, which is, however, within statistical errors. We thus conclude that the systematic errors from the derivative expansion and the inelastic state contaminations are well under control for this observable. Numerical values for the scattering length are summarized in Table III, where the central value and statistical errors are evaluated at t = 13 and the systematic errors are estimated from the t-dependence among t = 13 − 16.
We have checked that alternative fitting functions of the potential such as the combination of two Gaussians + (Yukawa) 2 form as employed in [9,22] give results consistent with those from the present fitting function within errors.

V. FINITE VOLUME FORMULA AND EFFECTIVE RANGE EXPANSION
Before closing the paper, we discuss the relation among the energy spectrum, the Lüscher's finite volume formula and the effective range expansion (ERE). Once the energy shift of the two-body system on a finite volume is measured, the scattering phase shift (r) (black diamonds) ,  is obtained by the Lüscher's formula as where k 2 is related to the energy shift on a finite volume as ∆E L = 2 m 2 B + k 2 − 2m B . For the attractive interaction, k 2 can be negative on a finite volume. Note that the poles of the S-matrix with k cot δ 0 (k) = − √ −k 2 in the infinite volume correspond to the bound states.
For the unbound two-body system, the asymptotic behavior of ∆E L for large L reads with the reduced mass µ, the scattering length a 0 , c 1 = −2.837297, and c 2 = 6.375183 [4,5].
Let us now calculate k 2 from eigenvalue spectra  Table II. Fig. 12 (Left) shows the volume dependence of the lowest eigenvalues: The data are found to be well described by Eq. (17), which indicates that the system does not have a bound state. By fitting the data with Eq. (17), we obtain the scattering length as (a 0 m π ) −1 = 0.402 (14) consistent with the values in Table III, (a 0 m π ) −1 = 0.352(36)( +80 −0 ). As extensively discussed in Ref. [12], the ERE, k cot δ 0 (k) = 1/a 0 + (1/2)r eff k 2 + · · ·, provides a systematic and reliable way to relate the volume dependence of ∆E L , the scattering phase shifts and the bound state pole around k 2 = 0. 6 In Fig. 12 (Right), we plot the finite volume spectra on the (k 2 , k cot δ 0 (k)) plane, using the lowest eigenvalues of H on L = 40, 48, and 64, and the eigenvalue of the first excited state on L = 64. Note that the data (triangle, square and diamonds) and their errors are plotted together with the Lüscher's formula (dotted lines). The blue band corresponds to the results obtained by solving the Schrödinger equation in the infinite volume. We find that the finite volume energy spectra at k 2 < 0 and k 2 > 0 are smoothly connected around k 2 = 0 along with the blue band, as is expected from the analytic properties of S-matrix and the ERE. In fact, the ERE at the NLO determined from these 4 data (pink band) is consistent with the blue band at |(k/m π ) 2 | < ∼ 0.2 within errors. One also observes that the positive intercept at k 2 = 0 (1/a 0 ) supports the conclusion from Fig. 12 (Left) that the system has no bound state.

VI. SUMMARY
In this paper, we have made critical investigations on the systematic uncertainties in the HAL QCD method. While the time-dependent HAL QCD method is free from the issue associated with the ground state saturation, the approximation of the energy-independent non-local potential by the derivative expansion introduces systematic uncertainties, so that it is necessary to check the errors introduced by the expansion. and N 2 LO potentials in the derivative expansion. Scattering phase shifts calculated from these potentials reveal that the LO potential is sufficient to reproduce observables at low energies (k 2 /m 2 π < 0.1), while the N 2 LO correction becomes non-negligible but remains small even at high energies (k 2 /m 2 π 0.5), confirming the good convergence of the derivative expansion below the inelastic threshold for this particular system.
We have also found that the potential at the LO analysis for the wall source agrees with the LO potential at the N 2 LO analysis except at short distances and can reproduce the scattering phase shifts precisely at low energies. Other systematic uncertainties such as the inelastic state contaminations and the finite volume effect to the potential are investigated and are found to be well under control.
After establishing the reliability of the HAL QCD potential, we have calculated the eigenvalues of the Hamiltonian in finite boxes with the potential. The volume dependence of the lowest eigenvalues is well described by 1/L-expansion for scattering states obtained from the Lüscher's finite volume formula. We have also discussed the relation among the energy spectrum, phase shifts and the effective range expansion.
In a forthcoming paper [14], we will perform the spectral decomposition of the correlation function based on the eigenmodes of the Hamiltonian in a finite box with the HAL QCD potential, which will enable us to better understand the requirements for reliably extracting finite-volume energies.
Appendix A: Non-locality vs. Energy dependence Here we examine the relation between the energy-independent non-local potential with the derivative expansion, U ( r, r ) = {V 0 (r) + V 2 (r)∇ 2 + · · ·}δ( r − r ), and the energydependent local potential, V eff (r; E). For simplicity, in this appendix, we restrict ourselves to the N 2 LO analysis. In other words, we assume as if the non-local potential were given exactly by U ( r, r ) = {V 0 (r) + V 2 (r)∇ 2 }δ( r − r ).
In this case, it is easy to show that the Schrödinger equation with this non-local potential, given by can be written in terms of the energy-dependent local potential as where which gives an exact relation between the energy-independent non-local potential and the energy-dependent local potential (within the N 2 LO analysis). Although both descriptions for the potential are theoretically equivalent as shown above, we stress that the HAL QCD method is based on the energy-independent non-local potential, which can be extracted from arbitrary linear combinations of the NBS wave function ψ W ( r) thanks to the time dependent method, while the energy-dependent local potential requires the eigenstate saturation, which is difficult to achieve in practice, particularly for excited states. Also V eff (r; E) gives the correct scattering phase shift at each E (one potential per energy), while V 0 (r) + V 2 (r)∇ 2 gives the correct scattering phase shifts (within the N 2 LO analysis) at all E ≤ E th (one potential for all). In these figures, we use V N 2 LO 0 (r) and V N 2 LO 2 (r) obtained at t = 13 for V 0 (r) and V 2 (r), respectively. The energy dependent correction is small at low energies, while it is no longer negligible at higher energies. As the energy increases, the attractive pocket at an intermediate distance becomes shallower and the radius of the repulsive core becomes larger.
This result also demonstrates how the non-locality of the energy-independent potential, U (r, r ), (Note that V 0 (r), V 2 (r), · · · are energy-independent by definition 7 ), is related to the energy dependence of the local potential, V eff (r; E). 7 In the literature, there appears a confusion on the relation between the energy-independent non-local potential and the energy-dependent local potential. See Ref. [30], which clarifies the relation between the two in detail.