The HAL QCD potential in $I=1$ $\pi \pi$ system with the $\rho$ meson bound state

In this paper, we investigate the HAL QCD potential in the $I=1$ $\pi \pi$ scattering using the hybrid method for all-to-all propagators, in which a propagator is approximated by low-eigenmodes and the remaining high-eigenmode part is stochastically estimated. To verify the applicability of the hybrid method to systems containing quark creation$/$annihilation contributions such as the $\rho$ meson, we calculate the $I=1$ $\pi\pi$ potential 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},m_{\rho}) \approx (870, 1230)$ MeV, in which the $\rho$ meson appears as a deeply-bound state. While we find that the naive stochastic evaluations for quark creation$/$annihilation contributions lead to extremely large statistical fluctuations, additional noise reduction methods enable us to obtain a sufficiently precise potential, which shows a strong attractive force. We also confirm that the binding energy and $k^3 \cot \delta$ obtained from our potential are roughly consistent with an existing $\rho$ meson bound state, within a large systematic error associated with our calculation, whose possible origin is also discussed.


Introduction
One of the most challenging issues in particle and nuclear physics is to understand hadronic resonances in terms of the fundamental theory of quarks and gluons, Quantum Chromodynamics(QCD). To achieve this goal, two methods to study hadron-hadron interactions non-perturbatively in lattice QCD have been employed so far: the Lüscher's finite volume method [1][2][3] and the HAL QCD method [4][5][6][7]. The Lüscher's finite volume method enables us to calculate scattering phase shifts directly from finite-volume energy spectra. Pole structures of bound states and resonances can be addressed by the analytic continuation of the S-matrix into the complex energy region, which however would require some ansatz for the structure of the S-matrix, in particular for coupled channel systems. Until now, several mesonic resonances, such as the ρ meson, have been studied in lattice QCD by this method [8][9][10].
In the HAL QCD method, on the other hand, an energy-independent but non-local potentials of hadron interactions are constructed from the Nambu-Bethe-Salpeter(NBS) wave function calculated in lattice QCD, from which physical observables are extracted afterward. This method has a unique advantage for the understanding of hadronic resonances from the first-principle. In this method, once the potential is obtained, one can directly address the pole structure of the S-matrix without any additional model-dependent ansatz. The extension to coupled channel systems, which are often essential to understand resonances, can be achieved in a straightforward manner [11]. Another strengh of this method is that the signal of the potential can be extracted not only from the ground state but also from excited states, which is crucial for reliable calculations for baryon-baryon systems [7,12]. Various interesting results have been reported in this method, for example, the identification of the Z c (3900) as the threshold cusp effect [13,14] and predictions on the existence of ΩΩ and N Ω di-baryons at the physical point [15,16].
At present, however, studies of resonances with the HAL QCD method are restricted due to the difficulty to treat all-to-all propagators within reasonable numerical costs and sufficient precisions. In our previous attempts [17,18], we utilized the LapH method [19] to treat all-to-all propagators, and it was revealed that the LapH smearing on the sink operator enhances the non-locality of the potential, so that the leading order approximation in the derivative expansion for the potential become insufficient. To establish a more suitable way for all-to-all propagators, we have recently applied the hybrid method [20], which treats all-to-all propagators by the low-eigenmode approximation plus the stochastic estimation for the remaining high modes, to the HAL QCD method [21]. In contrast to the LapH method, the hybrid method can keep the locality of quark operators since it contains full information of eigenmodes of the Dirac operator. In Ref. [21], we have studied the I = 2 ππ S-wave scattering with the hybrid method, and we confirmed that the combination of the HAL QCD method and the hybrid method gave us reliable results with better convergence of the derivative expansion, as long as appropriate choices of parameters for the hybrid method have been made.
In this paper, we apply the hybrid method to I = 1 ππ system and study the ρ meson. Since all-to-all propagators are mandatory to calculate quark creation/annihilation contributions, this is a best benchmark system to verify the applicability of the hybrid method. We calculate the potential on the gauge configurations at (m π , m ρ ) ≈ (870, 1230) MeV, in which the ρ meson is not a resonance but a deeply-bound state. It is revealed that stochastic estimations in the hybrid method for quark creation/annihilation contributions extremely enhance statistical fluctuations of the HAL QCD potential, and therefore we have to take some additional noise reductions to obtain a sufficiently precise potential. As a consistency check, we calculate the binding energy and k 3 cot δ from the resultant potential, and confirm that a deeply-bound ρ state is reproduced within a somewhat large systematic error. This paper is organized as follows. In Sec. 2, we briefly explain the HAL QCD method and the hybrid method. Simulation details in this study are given in Sec. 3. Our main result, the potential of the I = 1 ππ system, is presented in Sec. 4. We also discuss physical observables computed by the potential and the origin of their systematic uncertainty here. Our conclusion and outlook are given in Sec. 5.

HAL QCD method
The fundamental quantity in the HAL QCD method is the Nambu-Bethe-Salpeter (NBS) wave function, which is defined for the I = 1 two-pion system as ψ W (r, ∆t) = 0|(ππ) I=1,Iz=0 (r, 0, ∆t)|ππ; I = 1, where |ππ; I = 1, I z = 0, k is an asymptotic state for the elastic I = 1 ππ system in the center-of-mass frame with a relative momentum k, the total energy W = 2 m 2 π + k 2 and k = |k|. The operator (ππ) I=1,Iz=0 (r, t, ∆t) is a local two-pion operator projected to the I = 1, I z = 0 channel, explicitly given by where π + (x, t)(π − (x, t)) is the positively (negatively) charged pion operator defined as ) with up and down quark fields u(x, t) and d(x, t).
The above definition of the NBS wave function is more general than the equal time (∆t = 0) NBS wave function, conventionally employed in the HAL QCD method, where two sink hadron operators are put on the same time slice. In general, the HAL QCD potential depends on the choice of hadron operators in the definition of the NBS wave function, and we call it "scheme"-dependence of the potential [17,22]. The potential derived from the NBS wave function with ∆t = 0 belongs to the same scheme as the conventional equal time scheme if we take ∆t → 0 in the continuum limit, while it belongs to a different scheme if we keep physical ∆t finite in the continuum limit. Since calculations in this study are performed only at one lattice spacing, we consider the results from ∆t = 0 and ∆t = 0 as those in two different schemes between which the discretization artifact appears differently.
While the potentials are scheme-dependent, physical quantities such as phase shifts and binding energies, of course, do not depend on the scheme (up to the discretization errors). On can even take advantage of this arbitrariness by choosing a better scheme so that statistical/systematic errors are minimized. As discussed later, the main reason why we introduce the scheme with non-zero ∆t is to reduce statistical fluctuations of the potential for the I = 1 ππ system, which are caused by stochastic estimations for all-to-all quark propagators in the hybrid method.
As discussed in Ref. [5,23] for the case of the ∆t = 0 scheme, we can show the radial part of the l-th partial component in the NBS wave function with the non-zero ∆t scheme behaves at large r = |r| as where A l (∆t, k) is an overall factor and δ l (k) is the scattering phase shift, which is equal to the phase of the S-matrix implied by its unitarity. By using this behavior, we can construct an energy-independent but non-local potential through the Schrödinger-type equation as where µ = m π /2 is a reduced mass of two-pions, and a subscript ∆t of U represents the scheme for the potential. In practice, the non-locality of the potential is treated by the derivative expansion as The normalized ππ correlation function, numerically calculable in lattice QCD, is related to the NBS wave functions as where W n and B n are the energy and overlap factor of the n-th excited elastic state, and an ellipsis indicates inelastic contributions. Here C(t) and F (r, t, ∆t) are π and ππ correlation functions defined by where ρ 0 i is the neutral ρ meson operator, ρ 0 i =ūγ i u −dγ i d. Since this source operator strongly overlaps with the ρ meson state, we expect that the truncation error of the derivative expansion in the effective leading-order analysis is suppressed around the mass of the ρ meson.
The normalized correlation function R(r, t, ∆t) satisfies [7] at a sufficiently large t where inelastic contributions in R(r, t, ∆t) becomes negligible. From eq.(10), the effective leading-order(LO) potential is obtained as Using the rotational invariance of the potential, we can rewrite the above definition to improve signals as [24] V LO ∆t (r) = where the O h is the cubic rotation group. We also note that we employ the 4th order difference approximation for ∇ 2 to reduce discretization errors at short distances, since it turns out that physical observables in the deeply-bound system are sensitive to the potential at short distances.

All-to-all propagator: the hybrid method
In this subsection, we briefly explain the hybrid method, a technique for the all-to-all propagator calculation employed in this study. Let us consider the spectral decomposition of the quark propagator as where v (i) (x) and λ i are eigenvectors and eigenvalues of the Hermitian Dirac operator H = γ 5 D, respectively, with N being the total number of eigenmodes, and color and spinor indices are implicit for simplicity. We here assume |λ i | ≤ |λ j | for i < j. The low-eigenmode approximation for the propagator with the spectral decomposition is introduced as while the remaining high-eigenmode part is estimated by using the Z 4 noise vector η (i) [r] , together with the variance reduction by dilution as where is a projection onto the remaining high-eigenmode part, N r (N dil ) is a number of noise vectors (dilutions), and ψ (i) [r] are solution vectors obtained by solving H · ψ (i) [r] = P 1 η (i) [r] . In this study, the temporal coordinate is diluted with the J-interlace as For spatial coordinates, we introduce not only s2 (even-odd) and s4 dilutions used in the previous study [21], but also a s8 dilution. In the s8 dilution, one noise vector is split into (n x , n y , n z ) = (odd,odd,odd) and n x + n y + n z = 1 mod 4 (n x , n y , n z ) = (even,even,even) and n x + n y + n z = 2 mod 4 (n x , n y , n z ) = (odd,odd,odd) and n x + n y + n z = 3 mod 4 (n x , n y , n z ) = (even,even,even) and n x + n y + n z = 0 mod 4 (n x , n y , n z ) = (odd,even,even) and n x + n y + n z = 1 mod 4 (n x , n y , n z ) = (even,odd,odd) and n x + n y + n z = 2 mod 4 (n x , n y , n z ) = (odd,even,even) and n x + n y + n z = 3 mod 4 (n x , n y , n z ) = (even,odd,odd) and n x + n y + n z = 0 mod 4 (n x , n y , n z ) = (even,odd,even) and n x + n y + n z = 1 mod 4 (n x , n y , n z ) = (odd,even,odd) and n x + n y + n z = 2 mod 4 (n x , n y , n z ) = (even,odd,even) and n x + n y + n z = 3 mod 4 (n x , n y , n z ) = (odd,even,odd) and n x + n y + n z = 0 mod 4 (n x , n y , n z ) = (even,even,odd) and n x + n y + n z = 1 mod 4 (n x , n y , n z ) = (odd,odd,even) and n x + n y + n z = 2 mod 4 (n x , n y , n z ) = (even,even,odd) and n x + n y + n z = 3 mod 4 (n x , n y , n z ) = (odd,odd,even) and n x + n y + n z = 0 mod 4 . (17) See Fig. 1 for a schematic figure of the s8 dilution. Color and spinor indices are fully diluted in this study.
Combining the low-eigenmode and the high-eigenmode parts, the all-to-all propagator is written as where the hybrid lists u with N hl = N eig + N dil .

Correlation function with the hybrid method
The correlation function F (r, t, ∆t) with the ρ-type source operator, is expressed in terms of the hybrid method (up to an overall sign) as where Note that equal-time quark propagators would appear due to contractions in the sink operator if we took ∆t = 0 in the calculation.

Simulation details
In this study, we employ 2+1 flavor full QCD configurations generated by JLQCD and CP-PACS Collaborations [25] on a 16 3 × 32 lattice with the Iwasaki gauge action [26] at β = 1.83 and a non-perturbatively improved Wilson-clover action [27] at c SW = 1.7610 and hopping parameters (κ ud , κ s ) = (0.1376, 0.1371). These parameters correspond to the lattice spacing a = 0.1214 fm, the pion mass m π ≈ 870 MeV, and the ρ meson mass m ρ ≈ 1230 MeV. Note that the ρ meson is not a resonance but a bound state of two pions in this calculation. The periodic boundary condition is employed for all spacetime directions.  Tab. 1 shows details of our numerical setup, whereas parameters for the hybrid method are summarized in Tab. 2. In case 0, the source operator in eq. (9) is constructed from the point quark source. In case 1, on the other hand, we employ the smeared quark source q s (x, t) = y f (x − y)q(y, t) with the Coulomb gauge fixing, so that inelastic contributions are reduced at earlier imaginary times. The smearing function f is given by [29] f with a = 1.0, b = 0.47 in lattice unit. As regards the setup for the random noise vectors, case 0 is calculated with three independent Z 4 noise vectors corresponding to r, s, p in eq. (22). In case 1, we generate four different sets of three Z 4 noise vectors, and take an average over 4 × 3! = 24 samples (3! = 6 samples for each set using the permutation of r, s, p) to reduce noise contamination. Statistical errors are estimated by the jackknife method with bin-size 1 (6) in case 0 (case 1).
In case 0, we employ the equal-time (∆t = 0) scheme. As will be shown in Sec. 4.1, however, the statistical errors of the potential are found to be too large to obtain physical results, probably due to the statistical fluctuations associated with the equal-time quark propagations in the sink operator. We therefore employ ∆t = 1 scheme in case 1 and avoid equal-time quark propagations. In addition, we make a spatial dilution finer in the sinkto-sink propagator to reduce noise contamination in the Laplacian part, whose increased numerical costs are partly compensated by decreasing the temporal dilution from the 16interlace to the 4-interlace. Since we found in the previous study [21] that propagations along the temporal direction from t 0 to t 0 + t with t < J/2 are not distorted much by the J-interlace dilution, the 4-interlace temporal dilution reduces the computational cost for sink-to-sink propagations without additional strong noise enhancements. tion [21]) obtained in case 1 for pion m π (t) and ρ meson m ρ (t), which are calculated from C(t) and F (t) ≡ rȲ l=1,m=0 (Ω r )F (r, t, ∆t), respectively. Note that we insert the spherical harmonics for the P-waveȲ l=1,m=0 in the summation to obtain F (t), which is relevant to the ρ meson. The fit to C(t) at t = 4 − 11 gives m π = 871(4) MeV, while the fit to F (t) at t = 6 − 11 gives m ρ = 1228(5) MeV. The ratio of m π and m ρ becomes m π /m ρ = 0.709(4), which is consistent with m π /m ρ = 0.7076 (18) reported in the previous study [25]. Fig. 2 also shows that the ground state saturations in C(t) and F (t) are achieved at least t = 4 and t = 6, respectively, in case 1. In case 0, while the ground state saturation in C(t) is achieved at later time than case 1 (see Fig. 2 (Left) in [21]), t = 6 is found to be sufficient since errors in the potential is dominated by the statistical fluctuations as will be shown in Sec. 4.1. In the following, we take results at t = 6 as our central values and use results at t = 5, 7 to estimate systematic errors associated with their time dependence.

Potential in case 0
We first consider the case 0 for the I = 1 ππ potential, whose setup for the hybrid method is the same as the case 3 for the I = 2 ππ potential in Ref. [21]. In the previous study, we have found that the I = 2 ππ potential is reasonably accurate at t < 8. Fig. 3(Left) shows the potential obtained at t = 6. As can be seen, the potential has extremely large statistical fluctuations in this setup. Since equal-time quark propagations at the sink were absent for the I = 2 ππ potential in the previous study, we suspect that extremely large statistical fluctuations for the I = 1 ππ potential are caused by noise contaminations from the hybrid method to evaluate such equal-time propagations at the sink.
To suppress such noise contaminations, we additionally employ three noise reduction  techniques, (1) the different-time scheme for the NBS wave function to avoid the equaltime propagation, (2) the finer space dilution in the quark annihilation part to reduce noise contamination in spatial indices, (3) the average over different noise vectors. In the following, we will show the result in case 1 with these three improvements, whose details were already explained in Sec. 3.

Potential with additional noise reductions in case 1
The potential in case 1 at t = 6 is shown in Fig. 3 (Right). Since the smeared quark sources are employed in case 1, t = 6 is large enough to suppress elastic contributions to the potential. Thanks to additional noise reduction techniques mentioned in Sec. 4.1, statistical fluctuations of the potential are drastically reduced. The potential shows a strong attraction without repulsive core, which is consistent with existence of the deeplybound ρ meson in this system. As shown in Fig. 4, the potentials is almost independent of time at t = 5, 6, 7, as expected from the effective energy shown in Fig. 2. Interestingly, we notice that statistical fluctuations of the potential increase as the distance r increases. We interpret this behavior qualitatively as follows. Two-pion scattering states give dominant contributions to the long-distance part of the potential, as the two-pion sink operator in the NBS wave function at large r strongly couple to them. The bound ρ meson state, on the other hand, give large contributions to the short distant part of the potential. Since the ρ-type operator we employ at the source hardly creates such two-pion scattering states, it is hard to determine the long-distance part of the potential precisely, and thus statistical fluctuations become large. We also observe that the short-distance part of the potential has non-smooth behaviors, which probably come from higher partial wave contaminations, for example, the l = 3 partial wave in our case, as similar behaviors have been sometimes observed for the HAL QCD potentials in previous studies and the rotational breaking by the discretization artifact is expected to be enhanced at short-distance.
To calculate physical observables such as binding energies and scattering phase shifts, we fit the potential on discrete lattice points by a sum of three Gauss functions given by Several issues for the fit of the potential are in order here. The first one is the finite volume effect. As seen in Fig. 3 (Right), the potential deviates from zero even at r = La/2 = 0.9712 fm due to the finite volume effect of the periodic boundary condition. We therefore partly include this finite volume effect into the fit as V (r) PBC = V (r) + n∈{(0,0,±1),(0,±1,0),(±1,0,0)} V (r + Ln).
The second issue is the non-smooth behavior of the potential at short distance, as mentioned before. To make the fit stable, we have to exclude two points of the potential at r = 0.2428 and 0.3642 fm, which largely deviate from other data points. We expect that the exclusion of these points partly reduces the systematic uncertainty associated with the contaminations from higher partial waves at short distances. We leave a more detailed analysis for future investigations with finer lattices and a new partial wave decomposition method [28]. Tab. 3 gives the result of the fit at t = 6 and Fig. 5 shows the original potential and the fitting result. Note that the χ 2 /d.o.f.=7.59 is much larger than 1 even with the exclusion of two data points at r = 0.2428 and 0.3642 fm in the fit, since remaining data points at a short distance still have scattered central values with small statistical errors.

Physical observables
Using the potential given by eq. (25), we calculate the ground state energy of the I = 1 ππ system in the infinite volume. We employ the Gaussian expansion method(GEM) [30] to evaluate the ground state energy, which is given by where the first error denotes the statistical error and the second error the systematic one estimated by the time dependence of the binding energy at t = 6 ± 1. Comparing with the binding energy E bind = |m ρ − 2m π | ≈ 515 MeV from m π and m ρ (See Sec. 3), the results are consistent with each other within a large systematic error in eq. (27). We also remark the systematic error associated with the fit of the potential. As mentioned in the previous subsection, some unreliable points have to be excluded in the fit, and data points at short distances are still scattered with small statistical errors, which leads to large χ 2 /d.o.f. In such a situation, the fit at the short-range part as well as the resultant binding energy may have additional large uncertainty, since the latter is rather sensitive to the structure of the potential at short distances.
While the corresponding systematic error is not fully quoted in eq. (27), part of such a systematics seems to be reflected in the systematic error estimated from the time dependence. In fact, we find that the time dependence of the results is substantial even though the potential is rather time independent as shown in Fig. 4. This indicates that the large time dependence is mostly originated from the uncertainty in the fit of the potential. To make systematic uncertainties fully under control, we need to employ calculations at finer lattice spacings to obtain more data points at short distances or to find a better scheme for the NBS wave function to have smoother behaviors at short distances. Having remarked the above open issue, we can still positively conclude that it is possible to calculate reasonably precise potentials in the systems including quark creation/annihilation processes by the combination of the hybrid method and the HAL QCD method.
We finally discuss a relation between k 3 cot δ 1 (k) and the bound state pole in detail, as the normality check proposed in Ref. [29]. In the P-wave scattering, k 3 cot δ 1 (k) is related to the scattering S-matrix S 1 (k) as Generally, the scattering S-matrix in P-wave near the bound state pole (k ≈ iκ b ) behaves as [31] where κ b is an absolute value of k of the pole and β 2 b is positive real constant related to the normalization factor of the wave function of the bound state. By using Eq. (28) and (29), the physical pole condition in P-wave becomes In Fig. 6, we show typical behaviors of k 3 cot δ 1 (k) calculated by the square well potential in several cases. We can see how k 3 cot δ 1 (k) evolves when the attraction becomes stronger from Fig. 6 (a) to Fig. 6 (c). As seen in Fig. 6 (b), the deeply-bound state appears as the intersection (blue solid star) between −k 2 √ −k 2 (the bound state condition, black dashed line) and a branch of k 3 cot δ 1 (k) (red solid line) disconnected from a branch at the origin (k 2 = 0). Moreover, k 3 cot δ 1 (k) satisfies the physical pole condition, Eq.(30) (See Fig. 6 (b)(lower right)).
These two typical behaviors of k 3 cot δ 1 in the presence of one deeply-bound state in the P-wave are indeed observed for our data obtained from the potential at t = 6: Fig 7  (Left) shows that an intersection between −k 2 √ −k 2 (black dashed line) and k 3 cot δ 1 (k) (red solid line) in the branch disconnected from the origin appears at k 2 /m 2 π ≈ −0.623, corresponding to the GEM result, E bind ≈ 668 MeV. Shown in Fig 7 (Right) is (k 3 cot δ 1 (k) − (−k 2 √ −k 2 ))/m 3 π , and one can explicitly see how the physical pole condition is satisfied.

Summary and outlook
In this paper, we calculate the HAL QCD potential of the I = 1 ππ system at (m π , m ρ ) ≈ (870, 1230) MeV, using the hybrid method for all-to-all propagators. While statistical fluctuations in the straightforward calculation are found to be extremely large due to the   Figure 7: (Left) k 3 cot δ 1 (k) (red solid lines) calculated with the potential at t = 6, together with the bound state condition (black dashed line). Note that red solid lines diverge to ±∞ around k 2 /m 2 π ≈ −0.6. (Right) A difference between k 3 cot δ 1 (k) and −k 2 √ −k 2 around the intersection, together with k 2 obtained by the Gaussian expansion method (blue solid star). Only the central values are used for the visibility. quark creation/annihilation process, we have successfully obtained the precise potential by developing various noise reduction techniques such as space dilutions and the non-equal time scheme for the potential. We have calculated physical quantities such as the binding energy and phase shifts from the potential. It is observed that our potential reproduces the characteristic features of the deeply-bound ρ meson, whose binding energy is consistent with that obtained from the temporal correlation within a large systematic error in the former. The large systematic error in the present calculations is caused by the uncertainty of the fit for the potential at short distances, whose origin is attributed to the contaminations from higher partial wave components.
Finally, we would like to comment on further improvements to our calculation in the future. This and previous studies [21] on the hybrid method reveals that one can obtain reasonably precise potentials as long as appropriate setups of calculations are introduced, but on the other hand, it is also found that the numerical cost for noise reductions seems too large to perform such calculations on larger lattice volumes. Therefore, we have to investigate possibilities to achieve both small noise contamination and small numerical costs. Fortunately, we find that the combination of some other techniques such as the all-mode-averaging [32], the one-end trick and sequential propagators [33] is promising to achieve above requirements. As a first step toward this direction, we are now working on the ρ resonance at m π ≈ 410 MeV with new improved methods, and results will be reported in near future.