$I=2$ $\pi\pi$ potential in the HAL QCD method with all-to-all propagators

In this paper, we perform the first application of the hybrid method (exact low modes plus stochastically estimated high modes) for all-to-all propagators to the HAL QCD method. We calculate the HAL QCD potentials in the $I=2$ $\pi\pi$ scattering in order to see how statistical fluctuations of the potential behave under the hybrid method. All of the calculations are performed with the 2+1 flavor gauge configurations on $16^3 \times 32$ lattice at the lattice spacing $a \approx 0.12$ fm and $m_{\pi} \approx 870$ MeV. It is revealed that statistical errors for the potential are enhanced by stochastic noises introduced by the hybrid method, which, however, are shown to be reduced by increasing the level of dilutions, in particular, that of space dilutions. From systematic studies, we obtain a guiding principle for a choice of dilution types/levels and a number of eigenvectors to reduce noise contaminations to the potential while keeping numerical costs reasonable. We also confirm that we can obtain the scattering phase shifts for the $I=2$ $\pi\pi$ system by the hybrid method within a reasonable numerical cost, which are consistent with the result obtained with the conventional method. The knowledge we obtain in this study will become useful to investigate hadron resonances which require quark annihilation diagrams such as the $\rho$ meson by the HAL QCD potential with the hybrid method.


Introduction
Understanding hadronic resonances in terms of quantum chromodynamics (QCD) is one of the most challenging subjects in both particle and nuclear physics. In order to study hadronic resonances in QCD, we have to analyze hadron-hadron scatterings using nonperturbative methods such as lattice QCD, where two methods have been employed so far. One is the Lüscher's finite volume method (and the extensions thereof) [1][2][3], which has been mainly applied to meson-meson systems [4] including resonant states such as ρ [5,6] and σ [7]. The other one is the HAL QCD method [8][9][10][11], in which non-local but energy-independent potentials are constructed from the Nambu-Bethe-Salpeter (NBS) wave functions in lattice QCD, and then physical observables such as scattering phase shifts are extracted by solving the Schrödinger equations with the potentials. The HAL QCD method has been applied for a wide range of hadronic systems [12][13][14][15][16][17][18][19][20][21][22] including, in particular, the candidate of the exotic state, Z c (3900) [23]. The consistency between the Lüscher's method and the HAL QCD method for two-baryon systems is extensively studied in Refs. [24][25][26][27]. Recently, the first realistic calculation of baryon interactions with nearly physical quark masses is performed in the HAL QCD method [28] and ΩΩ and N Ω systems are found to form di-baryons located near the unitarity [29,30].
Although the HAL QCD method is a powerful method to study hadron-hadron scatterings, we still do not have a mature way to treat scattering processes containing quark annihilation diagrams, due to the need for all elements of propagators (so-called all-to-all propagators). Such quark annihilation diagrams typically appear in resonant channels, thus establishing an efficient technique to treat them is an urgent issue for deeper understanding of hadronic resonances in the HAL QCD method. A previous study [31], which utilized the LapH method [32] for the all-to-all propagator in the HAL QCD method, revealed that non-locality in the definition of Nambu-Bethe-Salpeter (NBS) wave function introduced by the LapH smearing enhances higher order contributions in the derivative expansion of the potential, so that the leading order potential suffers large systematic uncertainties even at low energies.
In this study, we employ another technique to obtain all-to-all propagator, namely the hybrid method [33], to the HAL QCD method for the first time. This technique combines a low-mode spectral decomposition of the propagator together with stochastic estimations for remaining high modes. Contrary to the LapH method, this method keeps the locality of quark fields and thus the locality of hadron operators in the HAL QCD potential. In order to investigate how statistical errors of potentials and physical observables behave under the hybrid method, we calculate potentials for the I = 2 ππ S-wave scattering on gauge configurations at m π ≈ 870 MeV with all-to-all propagators. Since the calculation of quark annihilation diagrams is not necessary, this system is suitable to benchmark the hybrid method. As a result, we find that the potential by the hybrid method gives us reliable results as long as statistical errors caused by stochastic noises are kept sufficiently small by noise dilution techniques. Our study also provides an optimal choice of parameters in the hybrid method to calculate the potential, which will be useful for future studies. This paper is organized as follows. In Sec. 2, after briefly explaining the HAL QCD method, we introduce the hybrid method and discuss the application to the HAL QCD method. Simulation details in this study are shown in Sec. 3. As main results of our study, we present systematic studies on behaviors of the HAL QCD potentials with the hybrid method in Sec. 4, and a comparison with previous results without all-to-all propagators in Sec. 5. Our conclusion is given in Sec. 6.

HAL QCD method
A basic quantity in the HAL QCD method is the Nambu-Bethe-Salpeter (NBS) wave function, which is given for the I = 2 two-pion system as where |π + π + ; k is an asymptotic state of the elastic π + π + S-wave system in the centerof-mass frame with a relative momentum k and an energy W = 2 m 2 π + k 2 , k = |k| and the positively charged pion operator is given as π + (x, t) =d(x, t)γ 5 u(x, t) with up and down quark fields u(x, t) and d(x, t).
For r = |r| > R with an interaction range R, the NBS wave function satisfies the Helmholtz equation as and its radial part with angular momentum l behaves as where δ l (k) is the scattering phase shift corresponding to the phase of the S-matrix constrained by the unitarity [9,34] and C is a constant. An energy-independent and non-local potential U (r, r ) is defined from where µ = m π /2 is the reduced mass of the two-pion. Thanks to Eq. (3), this potential is faithful to the scattering phase shifts by construction, but the potential depends on the scheme such as the choice of operators in the definition of the NBS wave function [9,10].
In practice, we deal with the non-locality of the potential in the derivative expansion as whose convergence property depends on the scheme of the potential. In the previous study [31], it is found that locality of the operator in the NBS wave function is a major factor which governs the non-locality of the potential: The scheme with non-local (smeared) operators leads to the potential with large non-locality and inclusion of V N 2 LO (r) is required to reproduce the I = 2 ππ scattering phase shifts even at low energies, while the scheme with local operators leads to the potential with small non-locality and V LO (r) alone is found to be sufficient. In this paper, we employ the scheme with local operators and calculate the potential at the leading order (LO), V LO (r). We define a 4-pt correlation function as where a source operator J π + π + (t 0 ) creates I = 2 π + π + scattering states in the A + 1 representation. In this paper, we employ a bi-local operator in which each pion operator is projected to zero-momentum, We then normalized it as which can be written as R(r, t) = n B n ψ Wn (r)e −∆Wnt + (inelastic contributions), where ∆W n = W n − 2m π and W n = 2 k 2 n + m 2 π is the energy of the n-th elastic state.
the time-dependent HAL QCD method reads [11] as long as t is large enough to suppress inelastic contributions in R(r, t). Thus the LO potential is given by

The hybrid method for all-to-all propagators
In this paper, we employ the hybrid method for all-to-all propagators [33], where the dominant part of the quark propagator is represented by low eigenmodes of the Dirac operators while remaining contributions from high eigenmodes are estimated stochastically. The quark propagator can be expressed in terms of eigenmodes of the hermitian Dirac operator H = γ 5 D as where λ i and v (i) are the i-th eigenvalue and eigenvector of H, respectively, and N is a total number of eigenmodes. Color and spinor indices are implicit here. We approximate the propagator using the low-lying N eig eigenmodes as which is expected to be a reasonable approximation for pions at low energies. The rest, can be estimated stochastically as where a noise vector η satisfies |η aα (x)| 2 = 1 (no summation) (18) for color indices a, b and spinor indices α, β, and ψ is the solution of H · ψ = P 1 η. The symbol indicates an expectation value over probability distribution of noise vectors. In practice, the exact noise average is approximated by the finite sum, together with the variance reduction by diluted noises as where N r is the total number of noises, N dil is the total number of dilution, and ψ [r] . In our study, color and spinor are fully diluted, while several types of dilutions are employed for time and space. The time coordinate is either fully or J-interlaced diluted by the noises as where i runs from 0 to J − 1, thus J = N t corresponds to the full dilution. On the other hand, no dilution, s2 (even/odd) dilution or s4 dilution is employed for spatial coordinates. Noise vectors for the s2 dilution are given by while those for s4 dilution read n y , n z ) = (even,even,even) or (odd,odd,odd) η (1) = 0 if (n x , n y , n z ) = (odd,even,even) or (even,odd,odd) η (2) = 0 if (n x , n y , n z ) = (even,odd,even) or (odd,even,odd) η (3) = 0 if (n x , n y , n z ) =(even,even,odd) or (odd,odd,even) .
In total, the all-to-all propagator in the hybrid method is given by where the hybrid lists u [r] are defined as [r] , ..., η [r] , ..., ψ with N hl = N eig +N dil . The HAL QCD potential with the local operator can be constructed by using this all-to-all propagator.

Correlation functions in the hybrid method
The 2-pt correlation function for the charged pion, is expressed as where and the · indicates an inner product in color and spinor indices. The 4-pt correlation function defined by leads to

Simulation details
We employ the same 2+1 flavor QCD ensemble in the previous study [31], generated by JLQCD and CP-PACS collaborations [35,36] on a 16 3 × 32 lattice with an Iwasaki gauge action [37] at β = 1.83 and a non-perturbatively improved Wilson-clover action [38] at c SW = 1.7610 and hopping parameters (κ ud , κ s ) = (0.1376, 0.1371), which corresponds to the lattice spacing a = 0.1214 fm and the pion mass m π ≈ 870 MeV. In addition to the local quark source, we use the smeared quark source [24], q s (x, t) = y f (x − y)q(y, t) with the Coulomb gauge fixing, where with a = 1.0, b = 0.47 in lattice unit. Note that, regardless of the type of quark sources, all calculations are made with wall sources at the hadron level (see Eq. (7)). The periodic boundary condition is used in all directions. The setups for the hybrid method in this study are presented in Tab. 1. In this study, we use a single noise vector for each propagator (N r = 1), and the noise vectors are generated by Z 4 random noises. Statistical errors are estimated by the jackknife method with binsize 1 except for the case5a, where binsize is 6. Fig.1 (Left) shows the effective masses of the single pion, where effective masses are calculated by solving the following equation, Note that we use a half-integer time convention for m eff , whose convention is also used for the effective energy shift. We find that the result from the smeared source reaches the plateau at earlier time slices than the one from the point source. The fit to smeared data at t = 4 − 11 gives m π = 870(4) MeV, while the fit to point data at t = 8 − 11 leads to m π = 874(8) MeV, and both agree with 870 MeV in the previous study [31]. Effective energy shifts for two pions from smeared and point sources are plotted in Fig. 1 (Right), where the energy shift is defined by ∆E ππ = E ππ − 2m π with the two pion energy E ππ . In our setup, an energy gap between the ground and the 1st excited states is estimated to be ∆E 1 ∼ 420 MeV in the non-interacting case, thus the ground state saturation is expected to be achieved at roughly t ∼ 1 ∆E 1 ∼ 3.8 in lattice unit. Fits give ∆E ππ = 9.7(0.7) MeV (t = 3 − 10) from the smeared source, and ∆E ππ = 11(1) MeV (t = 5 − 10) from the point source, which are consistent with ∆E ππ = 14(5) MeV in the previous study [31]. These results suggest that R(r, t) is dominated by the ground state actually at t ≥ 5, so that the leading order potential obtained by the time-dependent HAL QCD method becomes reliable at low energies.

The hybrid method and the HAL QCD potential
In this section, we show how statistical errors of the HAL QCD potential for the I = 2 ππ system depend on various setups of the hybrid method. We mainly discuss data at t = 6, which is sufficiently large for the elastic state domination as discussed in the previous section. In the followings, we use a quartet (time dilution, space dilution, N eig , source type) to specify calculation setups. Fig.2 shows the I = 2 ππ potential at t = 6 from the point source with the full timedilution (case1), together with its decomposition into the first (Laplacian), the second (the 1st time derivative) and the third (the 2nd time derivative) terms of the right-hand side in Eq. (12). Although a bulk behavior of the potential agrees with the previous result [31], the potential obtained by the hybrid method suffers from much larger statistical errors, which mainly come from the Laplacian term.  The potential from the hybrid method with case1 (full, none, 100, point) at t = 6 (blue circles), together with its breakdown to three contributions, the Laplacian (red triangles), the 1st time derivative (green squares) and the 2nd time derivative (yellow diamonds).

Dilution in spatial directions
In order to reduce noise contaminations in the Laplacian term, we introduce the s2 dilution for spatial coordinates (case2) in addition to the full time dilution. Shown in Fig.3 (Left) is the corresponding result together with that from the full time dilution only (case1) for comparison. As can be seen from the figure, statistical errors of the potential are much more reduced by the s2 dilution for spatial coordinates. Although a numerical cost in case2 becomes approximately two times larger than that in case1, the statistical errors decrease by a factor of ∼ 3, and therefore, we can conclude that the space dilution actually reduces the statistical noise. On the other hand, as seen in Fig.3 (Right), which compares the effective energy shifts ∆E eff between case1 and case2, their difference is moderate in contrast to the potentials, while the magnitude of errors becomes smaller by space dilutions. This observation suggests that the summation over spatial coordinates for the calculation of the energy shift reduces noise contamination, thanks to cancellation among different spatial points. We thus conclude that noise reduction by dilutions in spatial directions is more important for the HAL QCD potential than for the energy shift calculation, as the potential is extracted from spatial as well as temporal dependences of correlation functions.

Dilution in a temporal direction
In order to compensate the increased numerical cost by the introduction of the space dilution, we investigate a possibility to reduce the cost by employing less dilutions in a temporal direction. In Fig. 4, we compare the 16-interlace time dilution (case3) with the full (32-interlaced) time dilution (case2) at various t(= 4, 6, 8, 10) together with the s2 space dilution for both cases. We observe that statistical errors are comparable between two cases at small t (t = 4, 6), while errors in the 16-interlace time dilution (case3) becomes much larger than those in the full time dilution (case2) at larger t (t = 8, 10).  These behaviors may be qualitatively explained as follows. If a quark propagator from t 0 to t + t 0 is estimated by the hybrid method with the 16-interlaced time dilution, signals propagated from t 0 to t + t 0 behave as exp[−Et], while noises contaminated from t 0 + 16 to t + t 0 decrease as exp[−E|t − 16|]. Therefore, signals we need are larger than this type of noise contaminations at t < 8. On the other hand, the signals are largely contaminated by the noises at t ≥ 8, so that the potential cannot be reliably extracted.
This observation suggests that one can reduce numerical costs by J-interlace time dilutions, which reduce N dil to N dil × J/N t , while keeping the quality of the potential, as long as the potential is calculated at t < J/2. Therefore, the most efficient way to calculate the potential would be to combine J-interlace time dilution with as large J as possible together with a lattice setup which enables us to extract the potential at smaller t, such as the smeared source instead of the point source.

Effects of the smeared source
In this subsection, we investigate the behavior of statistical fluctuations of the potential in the case of the smeared source.
We first compare the potential from the smeared source (case4) with that from the point source (case3), keeping 16-interlace time dilution and s2 space dilution in both cases. Fig. 5 (Left) shows that the source smearing makes the potential noisier. This enhancement of noises by the smeared source may be explained by the fact that spatial summations in the source smearing accumulate fluctuations associated with noise vectors. In addition, the gauge fixing may make possible cancellations among gauge variant noises less effective. To reduce noise contaminations due to the source smearing, we introduce finer space dilution, s4. We compare case5 (16-interlace, s4, 100, smear) with case3 (16-interlace, s2, 100, point) in Fig. 5 (Right), which shows a finer space dilution makes statistical errors in the potential with the smeared source comparable to those in the point source calculation.
In Fig. 6, we compare time dependence of the potential with the smeared source and that with the point source. As seen in Fig. 6, while the potential from the point source (Left) has a significant t dependence at small t, the potential from the smeared source (Right) is almost t independent even at t = 2. Therefore, the smeared source actually enables us to extract a reliable potential at an earlier time than the point source. While the use of the smeared source may not be mandatory in the present case, the smeared source becomes more useful when the potential from the point source show slower convergences in time.

Dependence on N eig
We finally investigate how noises of the potential depend on a number of low modes N eig . (16int, s4, 100, smear) t = 1 (16int, s4, 100, smear) t = 2 (16int, s4, 100, smear) t = 3 (16int, s4, 100, smear) t = 4 (16int, s4, 100, smear) t = 5 (16int, s4, 100, smear) t = 6 Naively, we expect that statistical errors become smaller for larger N eig , since the relative segment of the propagator estimated exactly, D −1 0 , becomes larger. In order to confirm this point explicitly, we compare potentials at t = 6 between N eig = 100 (case3) and 200 (case6) with (16-interlace, s2, point) in Fig. 7 (Left), which shows that the potential with N eig = 200 (red) has smaller noise contamination than the one with N eig = 100 (blue). Therefore, we can confirm that our expectation is indeed the case. Next, in order to see which is more important for noise reductions, finer space dilution or larger N eig , we compare (16-interlace, s2, N eig = 484, smear) (case7) with (16-interlace, s4, N eig = 100, smear) (case5), keeping N hl = N eig + N dil = 868 same for both cases. Note that, in our setup, N hl is effectively a good measure for numerical costs because the most time-consuming part of our calculations is the contraction part, whose numerical costs are controlled by N hl . Fig. 7 (Right) indicates that larger N eig (case7, red) is a little better for the potential to have smaller noise contaminations than smaller N eig (case5, blue).
From these studies, taking larger N eig is slightly advantageous over finer space dilutions in our case. However, this relative efficiency could depend on the actual value of N eig and the lattice setup such as the size of the lattice volume. In particular, we cannot freely make N eig larger and larger, since numerical costs for the calculation of eigenmodes become non-negligible at some point. Therefore, increasing N eig as long as the numerical cost for eigenmodes remains to be subdominant will be the first guiding principle before performing the detailed optimization on N eig .

Lessons in this section
In order to extract the HAL QCD potential with the hybrid method for all-to-all propagators, lessons learned from investigations in this section are summarized as follows.
(1) A finer space dilution should be used to reduce noise contaminations to the potential, in addition to full color and spinor dilutions.
(2) The smeared source accumulates noise contamination, thus additional dilutions in spatial directions are mandatory. We should take the smeared source to extract the potential as small t as possible if potentials with the point source become reliable only at larger t.
(3) A J-interlace time dilution can be used for the potential to be extracted at t < J/2, to reduce N dil by a factor N t /J. (4) It is better to increase N eig until the costs for the calculation of eigenmodes become significant. A total number of N hl becomes N hl = 12 × J × 2 s/2 + N eig , where s corresponds to a level of the space dilution, s = 0 (no dilution), s = 2 (s2), s = 4 (s4).

Comparison with the result without all-to-all propagators
In this section, we compare the I = 2 ππ potential and the corresponding scattering phase shifts obtained in the hybrid method for all-to-all propagators with those in the standard method without all-to-all propagators. Fig. 8 (Left) compares the LO potential for the I = 2 ππ system obtained by the hybrid method (16-interace, s4, N eig = 100, smear) with N conf = 60 (case5a) at t = 6 with the one with the conventional setup in the HAL QCD method (32 wall quark sources/conf with N conf = 700) at t = 10. Both results agree with each other within statistical errors.
Then the potential is fitted by from which the I = 2 ππ scattering phase shifts are extracted. Fig. 8 (Right) shows the fit line, and Tab. 2 summarizes fit parameters.  In Fig. 9, we present the I = 2 ππ scattering phase shifts δ 0 (k) (Left) and k cot δ 0 (k) (Right) as a function of k 2 , together with results from the HAL QCD method with wall quark source and those from the Lüscher's finite volume method [31]. As expected from an agreement of the potential, we confirm that results by the hybrid method agree with the ones without all-to-all propagators (wall quark source) and the results from the Lüscher's method. Luscher's method Hybrid method t = 6 Wall source t = 10 Figure 9: (Left) Scattering phase shifts δ 0 (k) as a function of k 2 . (Right) k cot δ 0 (k) as a function of k 2 . Blue (red) bands correspond to the results from the HAL QCD method with the hybrid method (with the wall quark source). Black bands in the right figure correspond to the results from the Lüscher's method [31].

Conclusion
In this paper, we employ the hybrid method of all-to-all propagators [33] for the HAL QCD method and study the interaction of the I = 2 ππ system at m π ≈ 870 MeV. Even though the hybrid method brings extra statistical fluctuations for the results, we obtain the reasonably accurate potential by increasing dilution levels, which gives the I = 2 ππ scattering phase shifts consistent with the result using the conventional method. Our findings for appropriate choices of parameters in the hybrid method to calculate the potentials are summarized in Sec. 4.5.
As the hybrid method works in the HAL QCD method, we will calculate the I = 1 ππ potential using all-to-all propagators. It is interesting to see whether the ρ resonance is correctly reproduced by the HAL QCD potential. A preparatory study has already been made, and the results will be published in near future.