Emergence of the $\rho$ resonance from the HAL QCD potential in lattice QCD

We investigate the $I=1$ $\pi \pi$ interaction using the HAL QCD method in lattice QCD. We employ the (2+1)-flavor gauge configurations on $32^3 \times 64$ lattice at the lattice spacing $a \approx 0.0907$ fm and $m_{\pi} \approx 411$ MeV, in which the $\rho$ meson appears as a resonance state. We find that all-to-all propagators necessary in this calculation can be obtained with reasonable precision by a combination of three techniques, the one-end trick, the sequential propagator, and the covariant approximation averaging (CAA). The non-local $I=1$ $\pi \pi$ potential is determined at the next-to-next-to-leading order (N$^2$LO) of the derivative expansion for the first time, and the resonance parameters of the $\rho$ meson are extracted. The obtained $\rho$ meson mass is found to be consistent with the value in the literature, while the value of the coupling $g_{\rho \pi \pi}$ turns out to be somewhat larger. The latter observation is most likely attributed to the lack of low-energy information in our lattice setup with the center-of-mass frame. Such a limitation may appear in other P-wave resonant systems and we discuss possible improvement in future. With this caution in mind, we positively conclude that we can reasonably extract the N$^2$LO potential and resonance parameters even in the system requiring the all-to-all propagators in the HAL QCD method, which opens up new possibilities for the study of resonances in lattice QCD.


I. INTRODUCTION
Understanding the hadronic resonances from the first-principle lattice QCD simulation is one of the most important subjects in particle and nuclear physics.At present, the finite-volume method [1][2][3] and the HAL QCD method [4][5][6][7] are employed to extract hadron interactions.The finite-volume method extracts scattering phase shifts through finite-volume energy spectra obtained from temporal correlation functions.It has been successfully applied for various two-meson interactions and related mesonic resonances [8].
The HAL QCD method directly constructs inter-hadron potentials from spatial and temporal correlation functions calculated in lattice QCD.In this method, the potentials can be extracted even without ground state saturations for correlation functions [7].Scattering parameters are then obtained by solving the Schrödinger equation in infinite-volume without any model-dependent ansatz.These features make this method particularly useful to study baryonic systems [18] and coupled channel systems [19].The HAL QCD method has been successfully applied to many hadronic systems (see Ref. [20] and references therein for the recent status), including the detailed coupled channel studies for the tetraquark candidate Z c (3900) [21,22] and H-dibaryon [23].
There exists, however, a practical challenge in the HAL QCD method when expanding the scope to many other resonances, since the expensive computations of all-to-all quark propagators are necessary in most cases.To overcome this difficulty, we have previously explored two different all-to-all techniques, the LapH method [24] and the hybrid method [25].
It turned out [26,27] that the LapH smearing on sink operators with a small number of LapH vectors enhances non-locality of the HAL QCD potential and thus the systematic errors associated with the truncation of the derivative expansion.Increasing the number of LapH vectors to minimize such non-locality is practically impossible for larger volumes.On the other hand, the hybrid method with local sink operators is free from the enhancement of the non-locality and is more suitable for the HAL QCD method.A series of studies with the hybrid method [28,29], however, has revealed that it requires too much numerical cost for a reduction of stochastic noises to perform large-scale simulations for hadronic resonances.k = |k|.The operator (ππ) I=1,Iz=0 (r, t) is a two-pion operator projected to the I = 1, I z = 0 channel given by (ππ) I=1,Iz=0 (r, t) = 1 √ 2 {π + s (r + x, t)π − s (x, t) − π − s (r + x, t)π + s (x, t)}, π + s (x, t) = ds (x, t)γ 5 u s (x, t), π − s (x, t) = ūs (x, t)γ 5 d s (x, t), where u s , d s are smeared up and down quark fields.A detail of the quark smearing is given in Sect.III.A radial part of the l-th partial component in the NBS wave function behaves at large r = |r| as [5,34] ψ l W (r) ≈ A l (k)e iδ l (k) sin(kr where A l (k) is an overall factor and δ l (k) is a scattering phase shift, which is equal to a phase of the S-matrix implied by its unitarity.By using this property, we can construct an energy-independent non-local potential U (r, r ) as with µ = m π /2 a reduced mass of two pions.In general, the HAL QCD potential depends on a choice of hadron operators in the NBS wave function (the sink operator (ππ) I=1,Iz=0 (r, t) in our case), and it is referred to as the scheme dependence of the potential [26,35].Physical observables extracted from potentials in different schemes, of course, agree with each other by construction.Therefore, we can utilize this scheme dependence to reduce statistical and/or systematic uncertainties of observables.As discussed in Appendix C, a comparison of different schemes shows that the I = 1 ππ potential has much smoother r dependences in an "equal-time smeared-sink scheme" (Eq.( 2)), where sink quark fields are slightly smeared and two pion operators are put on the same time slice.We employ this scheme for the whole analysis in this study.
To extract the potential in lattice QCD simulations, we begin with a normalized correlation function defined as where F π and F ππ are a single-pion and a two-pion correlation function, respectively, x,y,t 0 π − (x, t + t 0 )π + (y, t 0 ) , Here J T − 1 I=1,Iz=0 (t 0 ) is a source operator which creates ππ scattering states with (I, I z ) = (1, 0) in an irreducible representation T − 1 of the cubic group.Thus R(r, t) is related to the NBS wave function as where W n and B n are energy and overlap factor of the n-th elastic state, and the ellipses indicate inelastic contributions.Using the energy independence of the potential, we can at a sufficiently large t where inelastic contributions in R(r, t) becomes negligible.In actual calculations, we introduce a derivative expansion to treat the non-local potential as and the effective LO potential is given by where invariance of the potential under the cubic rotation group O h is utilized to improve signals [36].In this study, we further determine the effective N 2 LO potential in order to extract resonance parameters more accurately.The effective N 2 LO potential U N 2 LO (r, r ) = ) is determined by solving the following linear equations [37]: where R i (i = A, B) are the normalized correlation functions with different source operators ).Note that coefficients V 0 (r), V 2 (r) in the full derivative expansion (Eq.( 11)) are independent of source operators, while effective depend on a choice of source operators.In other words, effective potentials implicitly depend on discrete energy levels included in their determination due to the truncation of the derivative expansion [20].Therefore, systematic errors in the derivative expansion for physical observables depend on the magnitude of non-locality in the true potential as well as on the difference between the energy region relevant for physical observables and that employed to determine the effective potentials [20].For the source operators, we choose ρ-type J where p 3 = (0, 0, 2π/L).(ππ) I=1,Iz=0 (p, t) and ρ 0 3 are given as where we use local quark fields for source operators.
Calculations of correlation functions with momentum projected sources generally need all-to-all propagators, which requires too much numerical cost to calculate exactly.Therefore, we evaluate all-to-all propagators by the combination of the one-end trick, the sequential propagator, and the CAA.We give a brief introduction of the one-end trick in Appendix A, and details of diagram calculations are presented in Appendix B.

III. SIMULATION DETAILS
We employ (2+1)-flavor full QCD configurations generated by the PACS-CS Collaborations [38] on a 32 3 × 64 lattice with the Iwasaki gauge action [39] at β = 1.90 and a non-perturbatively improved Wilson-clover action [40] at c SW = 1.715 and hopping parameters (κ ud , κ s ) = (0.13754, 0.13640).These parameters correspond to a lattice spacing a = 0.0907 fm, and a pion mass m π ≈ 411 MeV, where the ρ meson appears as a resonance with m ρ ≈ 892 MeV [9].The calculations are performed in the center-of-mass frame with the periodic boundary condition for all spacetime directions.In this report, dimensionful quantities without the corresponding unit are written in lattice unit unless otherwise stated.Table I and II show general setups and parameters of the one-end trick and the CAA, respectively.We employ smeared quark operators q s (x, t) = y f (x − y)q(y, t) at the sink with the Coulomb gauge fixing, in order to improve signals of potentials at short distance.
A smearing function f is given by with A = 1.0,B = 1.0,R = 3.5.As discussed in Appendix C, these parameters make potential smoother without worsening the convergence of the derivative expansion.For the one-end trick, we generate a single Z 4 noise vector for each insertion.To suppress the corresponding stochastic noises, we employ a dilution technique [25] in color, spinor and space indices.Color and spinor indices are fully diluted, and for the space dilution, we take s2 (even-odd) dilution and s4 dilution [29] in the ππ-type source and the ρtype source, respectively.In the CAA, we exactly estimate a low-mode part with 300 eigenmodes, and a high-mode part is estimated by an average over loosely solved solutions on 64 different spatial points x = (x 0 +8l, y 0 +8m, z 0 +8n) mod 32, with l, n, m ∈ {0, 1, 2, 3}.
Finer and looser solutions are obtained with 1.0 × 10 −24 and 9.0 × 10 −9 for the squared residue, respectively.We randomly choose the reference point x 0 = (x 0 , y 0 , z 0 ) for each configuration.
Figure 1 (left) show an effective mass of a pion obtained by an average over 200 configurations (×64 time slice average).A fit to the pion propagator F π (t) at t = [t min , t max ] = [14,29] with a cosh function gives m π = 413.5 (1.4) MeV.We also check a t min dependence of the effective mass, and the dependence is negligible compared with statistical errors as far as t min ≥ 13.Therefore we confirm that a ground state saturation in F π (t) is achieved at t = 13.A possible leading inelastic contribution for two pions in this setup comes from a P-wave KK state with energy W KK = 2 m 2 K + (2π/L) 2 ≈ 1530 MeV in non-interacting case, while the two-pion ground state energy is reported as E 0 = 914 (11) MeV in Ref. [9].We therefore expect inelastic contributions in F ππ (r, t) are suppressed at t ≈ 1/[W KK − E 0 ] ≈ 3.5.These considerations suggest that inelastic contributions in R(r, t) become negligible at t ≥ 13, so that potentials can be reliably extracted at t ≥ 13.Hereafter, we show results at t = 14 and 18 for ρ-type source and ππ-type source, respectively.
In lattice QCD, the rotational symmetry is broken to the cubic symmetry, and there exist higher partial wave components in the irreducible representation of the cubic group (l = 3, 5, • • • partial waves in this study).This leads to systematic uncertainties in the HAL QCD potential, which exhibit as multi-valued structures of potentials as a function of r.
We address this issue by performing the approximated partial wave decomposition recently introduced to lattice QCD [41].In practice, we remove the dominant contaminations, the l = 3 partial wave component, when we evaluate the potential at r = [2, 14.8].Tunable parameters of the decomposition [41], a number of radial bases n max , a number of partial waves considered l max and a width of the shell ∆, are taken as (n max , l max , ∆) = (4, 5, 1.2) at 2 ≤ r ≤ 10 or (4, 5, 1.5) at r > 10, where we use larger ∆ at larger r to avoid artificial oscillation of decomposed data due to too small ∆.  represents potentials after the partial wave decomposition with the P-wave centrifugal term added, which become much smoother as multi-valued structures are eliminated.The potentials with the centrifugal term reveal characteristic features for an existence of a resonance state such as an attractive pocket at short distances and a potential barrier around r = 0.5 fm.We also notice that potentials obtained from different source operators are different from each other, which fact suggests a presence of non-negligible higher-order contributions in the derivative expansion.We fit LO potentials with a sum of Gaussian terms given by For the fit, we utilize data projected to the l = 1 component by the partial wave decomposition at r = [2, 14.8] as already discussed, combined with the original lattice data at r ≤ 2, to which the partial wave decomposition cannot be reliably applied.We also remove data at very short distances (r = 0, 1) since they suffer from large discretization errors.Remaining systematic uncertainties caused by non-smoothness at short distances  are estimated by differences among three different fit results: a result using all allowed data (Fit), a result removing data at r ≤ 0.32 fm which are significantly larger than Fit (Fit − ), and a result removing data at r ≤ 0.32 fm which are significantly smaller than Fit (Fit + ).The fit results using all allowed data (Fit) are shown in Figure 3, and comparisons of three fit results at short distances are given in Figure 4. Resultant fit parameters and χ 2 /dof are given in Table .III and IV.As seen in Fig. 4, the non-smooth behavior of the potential at short distances, which is probably caused by contaminations from higher partial waves, affects the fit result at r 0.25 fm.Since the removal of such contaminations at short distances is impractical, we estimate systematic errors for physical observables by differences among fit results, taking the result using all allowed data as a central value.the derivative expansion for the LO potential.Since the ρ-type source strongly overlaps the ρ resonance state, which corresponds to the ground state in this setup, the resultant phase shift with the ρ source reproduces the ρ resonance structure relatively well.On the other hand, since the ππ-type source mainly overlaps P-wave ππ scattering states, which appear in the energy region far above the ρ resonance in this lattice setup, it is difficult for the phase shift with the ππ-type source to capture the resonance structure correctly.

B. The N 2 LO analysis
As we have two LO potentials, we can proceed to the N 2 LO analysis.The effective N 2 LO potentials are obtained through eq.( 13) as In Fig. 6 (upper left), we show V N 2 LO 2 obtained from raw data (blue points), and l = 1 data with the partial wave decomposition (red points).Thanks to the removal of higher partial wave contaminations, we can significantly reduce fluctuations of , as seen in the figure .A somewhat singular behavior at r ≈ 0.5 fm is caused by a vanishing denominator of in Eq.( 20).We however expect that this singular behavior is canceled by a vanishing numerator at the same point.As discussed in Appendix D, this expectation is shown to be true as long as the N 4 LO (and higher order) terms in the derivative expansion are negligible.Furthermore, we assume that 1 − 2µV N 2 LO 2 > 0, which is also shown to be true in Appendix D if the N 4 LO or higher order terms are negligible.We thus fit

2
(red point) by a smooth function, a 3-Gaussian function in Eq (19), where data with 1 − 2µV Let us consider a determination of V N 2 LO 0 (r) next.We first fit the Laplacian term  To obtain the N 2 LO phase shifts, we solve the radial Schrödinger equation with the N 2 LO potential, rewritten as The N 2 LO phase shifts and the corresponding k 3 cot δ 1 / √ s are shown in Figure 7, together with the LO phase shifts and the previous finite-volume result for comparisons.We have checked that the N 2 LO phase shifts do not vary beyond the magnitude of statistical errors even if we choose a different timeslice for the ππ-type source in the N 2 LO analysis.As can be seen in Fig. 7, except for the region s < 0.75 GeV 2 ( √ s < 870 MeV), the N 2 LO phase shifts and k 3 cot δ 1 / √ s are roughly consistent with the finite-volume results.
The deviation observed in the low-energy region can be understood from the truncation error of the derivative expansion as discussed in Sect.II.In this study, the calculations are performed only in the center-of-mass energy frame, where the corresponding energy levels on the current lattice volume do not cover the low-energy region near the ππ threshold.
Therefore, the N 2 LO approximation in this study could suffer from the large truncation error of the derivative expansion in such a low-energy region.This discrepancy actually affects a determination of some resonance parameters as will be discussed later.The detailed investigation is left for future studies since it needs much higher precision with possibly an additional technical development of the laboratory-frame calculation [33].

C. Resonance parameters
In this subsection, we extract resonance parameters for the ρ meson in the N 2 LO analysis using two different methods.

Breit-Wigner fit
We first extract resonance parameters in the conventional way, by fitting the scattering phase shifts with the Breit-Wigner form as where m ρ and g ρππ are fit parameters corresponding to a resonance mass and a ρ → ππ effective coupling, respectively.We show the fit result in Fig. 8, which gives g ρππ = 13.4(2.6)(+0.8 −0.0 ), with χ 2 /dof = 0.18, where the first errors are statistical and the second ones are systematic errors associated with the short-range behavior of V N 2 LO 0 .
We have checked that resonance parameters remain unchanged within statistical errors even if we add a centrifugal barrier modification as a higher order term in k 2 [42] to the standard Breit-Wigner form in Eq. ( 23).

Direct pole search
Theoretically, a resonance state is defined as a pole of the S-matrix on the second Riemann sheet, which provides us the second method to extract resonance parameters in the HAL QCD method.To access the S-matrix in complex energy region, we solve the Schrödinger equation with arguments rotated by r → re iθ , k → ke −iθ [43][44][45], which reads The regular solution φ to this equation behaves at long distances as where ĥ± l (z) = nl (z) ± i ĵl (z) are the Riccati-Hankel functions and J l is the Jost function for the angular momentum l.The S-matrix on the ray of ke −iθ can therefore be obtained as from which we can search a pole position k pole = |k pole |e −iθ pole by changing an input θ and k.The resonance mass and the decay width are extracted from the pole position √ s pole as where the decay width is related to the coupling constant g ρππ as The direct pole search gives g ρππ = 12.7(2.9)(+0.7 −0.0 ), (33) where the first errors are statistical while the second ones are systematic errors associated with the short-range behavior of V

Comparison to the previous result
Let us compare our N 2 LO results with the previous PACS-CS (2011) result using the finite-volume method [9], both of which employ the same gauge configurations.We plot m ρ and g ρππ in Fig. 9.While m ρ 's are consistent with each other in all three cases, the Breit-Wigner fit, the pole search and the PACS-CS (2011) result, coupling constants in our both results are about twice as large as the previous one.This discrepancy can be clearly seen as a difference in slopes of k 3 cot δ 1 / √ s data at s < 0.9 GeV 2 in Fig. 8, which directly correspond to the coupling as −6π/g 2 ρππ .In particular, a significant disagreement between the lowest energy level in PACS-CS and our data at s ≈ 0.75 GeV 2 is a main source of the discrepancy for the slope.We note that the lowest energy level in PACS-CS (2011) is obtained in the laboratory frame with P = (0, 0, 2π/L).Since such a low-energy region cannot be covered by the center-of-mass frame employed in this study and the non-locality of the potential in I = 1 ππ systems turns out to be large, our N 2 LO approximation is likely to suffer from large truncation errors in the derivative expansion at low energies.
This observation gives us a useful lesson for the study of P-wave (or higher partial wave) resonances by the HAL QCD method with the center-of-mass frame.If the nonlocality of the potential happens to be large, the truncation errors could be large at lowenergies near the threshold.While a resonance mass is likely to be well reproduced as long as the resonance appears in the energy region accessible in the center-of-mass frame, the decay width (the effective coupling) may suffer from larger systematics since it is sensitive to energy dependence on a much wider range around the resonance.As a possible option to control this systematics, if the resonance mass can be roughly guessed, one may tune lattice parameters such as a box size carefully so as to cover a wide energy range even in the center-of-mass frame.This procedure, however, is difficult in practice and applicability for searches of unknown resonances is also limited.The second option is to establish the existence of a resonance and to estimate its mass by the HAL QCD method in the center-of-mass frame, which is supplemented by the finite-volume method in the laboratory frame to estimate its width reliably.The third option is to perform the HAL QCD method with a combination of both the center-of-mass and laboratory frames.In fact, a theoretical framework has been already proposed for the HAL QCD method in the laboratory frame [33].While extraction of HAL QCD potentials from NBS wave functions in the laboratory frame is indeed a numerical challenge, the first numerical trial is now ongoing.

V. SUMMARY
We study the I = 1 ππ interaction at m π ≈ 411 MeV, where the ρ meson emerges as the resonance with m ρ ≈ 892 MeV.We calculate all-to-all propagators by a combination of the one-end trick, the sequential propagator, and the covariant approximation averaging.
Thanks to those techniques, we successfully determine the potential in this channel at the estimator.Let us consider a combination of quark propagators given by where D −1 f is a quark propagator with a flavor f , Γ is some product of gamma matrices, and x i = (x i , t i ) are arbitrary.We abbreviate color and spin indices for simplicity.Such a structure typically appears at the source side of correlation functions including meson operators.For example, in the separated diagram in Fig. 10, it appears twice as (+) The calculation of those structures naively needs two stochastic estimations for each, since each of them contains two all-to-all propagators.The one-end trick, however, utilize the γ 5 -Hermiticity of the Dirac operator to estimate that structure with a single noise insertion as follows.
where we insert the stochastic estimator δ y,z ≈ 1

Nr
Nr−1 r=0 η [r] (z)η † [r] (y) in the second line and use the γ 5 -Hermiticity in the last line.We define "one-end vectors" as then the final expression becomes The one-end vectors ξ and χ are obtained by solving the linear equation Dξ = ηe ip•y and Dχ = γ 5 Γ † η, respectively.The dilution technique for noise reduction can be combined as well.This trick is particularly suitable for the HAL QCD method since it does not introduce any stochastic estimations at the sink side, which otherwise strongly affects spatial dependences of the NBS wave function.Moreover, a numerical cost and stochastic noises are also reduced in accordance with a decrease in the number of noise vectors.
FIG. 10.Representative diagrams appeared in this study.Blue solid, orange dashed and green dotted lines are calculated with the one-end trick, sequential propagator and point-to-all propagator, respectively.Statistical improvement by the CAA is also employed for green dotted lines.

Separated diagram
As already seen, a separated diagram in Fig. 10 is written as (B1) By using the one-end trick twice, we obtain where i, j are indices for dilutions and r, s distinguish independent noise vectors.In practice, the center-of-mass coordinate x is averaged over the whole spacetime to improve the statistical errors, .

Box diagrams
A box diagram shown in Fig. 10 is written as (B5) We next exactly calculate another all-to-all propagator D −1 (x + r, t; y 1 , t 0 ) by the sequential propagator technique [31], where we consider a linear equation with a sequential source vector e ipz•y 1 γ 5 ξ To increase statistics of the box diagrams, instead of an average over all x with an additional noisy estimation, we employ the covariant approximation averaging (CAA) for x, which is given by G box,imp x 0 ;t 0 (r, t) = G  where λ n and v (n) are the n-th eigenvalue and eigenvector of H, respectively, N low is the number of low-eigenmodes used in the CAA, while H −1 high,exact/relaxed is an inverse of H projected onto a space spanned by remaining high-eigenmodes solved with a tight/relaxed stopping condition.Since χ and ζ are already solved with high precision, we only relax a precision of the sink-to-sink propagator (green dotted line in Fig. 10).Furthermore, we averaged over all x in the low-eigenmode part to maximize statistics.Using the one-end trick for a summation over y, we obtain G tri x;t 0 (r, t) = (−) i χ (i) † γ 3 ,t 0 [r] (x, t)H −1 (x, t; r + x, t)ξ (i) † γ 5 ,t 0 [r] (x, t)v (n) (x, t)v (n) † (x + r, t)ξ (i) 0,t 0 [r] (r + x, t) + χ (i) † γ 5 ,t 0 [r] (x 0 , t)H −1 high,exact/relaxed (x 0 , t; r + x 0 , t)ξ the derivative expansion less reliable.Therefore we would like to check whether our sinksmearing scheme given in eq.( 2) is free from such a problem.For this purpose, we calculate I = 2 ππ potential in both point-sink and smeared-sink schemes and compare LO phase shifts between the two schemes.
Calculations of NBS wave functions in both schemes are performed by using the one-end trick with full color/spin dilution and s2 space dilution for a single Z 4 noise.A number of configuration is N conf = 10 (×64 timeslice average), and statistical errors are estimated by the jackknife method with bin-size 1.
Figure 13(left) shows effective LO potentials at t = 14.Potentials between two schemes show significantly different behaviors only at short distances, which however do not affect

FIG. 1 .
FIG.1.Effective mass of a pion (blue circles) and the fit result by a cosh function at t = [t min , t max ] =[14,29] (cyan solid line with bands).
FIG. 2. (Left) Effective LO potentials.Blue and red points show the results from the ρ-type source and the ππ-type source, respectively.Inset shows an enlarged view of potentials.(Right) Improved potentials obtained by the partial wave decomposition with the P-wave centrifugal term, V c (r) = 1 2µ

Figure 2 (
Figure 2 (Left) show the results for effective LO potentials without the partial wave decomposition.We observed that the potentials are attractive at all distances.Fig. 2 (Right)

FIG. 3 .FIG. 4 .
FIG. 3. (Left) Fit result with the ρ-type source.Inset shows an enlarged view of them.(Right)The same plot with the ππ-type source.Both results are obtained with all allowed data points.

FIG. 5 .
FIG.5.Phase shifts at the LO analysis.Blue (orange) band shows the ρ-type (ππ-type) source result.Statistical errors are given by dark color bands, whereas systematic errors estimated by three different fits at short distances are represented by light color bands.The previous finitevolume results by the PACS-CS Collaboration[9] are also given by navy stars for comparison.

Figure 5 4 r
Figure 5 shows phase shifts obtained from fitted potentials, where systematic errors associated with removals of data at short distances are shown by light color bands on top of statistical errors by dark color bands.Shown together is the previous finite-volume result reported in Ref. [9], which employs the same gauge configurations.The phase shift obtained with the ρ-type source crosses 90 degrees around √ s ≈ 870 MeV, while it only reaches around 130 degrees as the energy increases.On the other hand, the phase shift obtained with the ππ-type source crosses 90 degrees at much higher energy, around √ s ≈ 1050 MeV, with much broader width.These behaviors are probably caused by truncation errors of

2
determined from raw data (blue circles) and data obtained with the partial wave decomposition (red triangles).(Upper right) The fit result (cyan band) using the decomposed data (red triangles).(Lower) V N 2 LO 0 obtained by three fit results , Fit(cyan), Fit − (magenta) and Fit + (green).Shown together are the effective LO potentials for a comparison.
r) by the same 3-Gaussian function, and resultant parameters are given in TableVI.We then obtain V N 2 LO 0 (r) by combining all the fit results in Eq. (21).We estimate systematic errors of V N 2 LO 0 (r) at short distances through those of V LO ρ and ∇ 2 R ρ (r)/R ρ (r).

FIG. 7 .
FIG. 7. The N 2 LO phase shifts (left) and k 3 cot δ 1 / √ s (right), together with LO results (left figure) and previous finite-volume result by the PACS-CS Collaboration [9] (both figure) for comparisons.Large statistical errors at s > 0.9 GeV 2 in k 3 cot δ 1 / √ s (right) are mainly caused by a divergent behavior of cot δ at the phase shift around 180 degrees.

)FIG. 8 .
FIG.8.The Breit-Wigner fit for the N 2 LO phase shifts k 3 cot δ 1 (k)/ √ s (blue points).The green band represents the fit with statistical errors and a range of the energy used in the fit.We also show data and Breit-Wigner fit of PACS-CS (2011)[9] by the black points and the dashed line, respectively, for a comparison.

FIG. 9 .
FIG. 9. A comparison of N 2 LO resonance parameters from the direct pole search (blue circle) and the Breit-Wigner fit (red triangle).Vertical and horizontal axes represent the coupling g ρππ and the mass m ρ , respectively.Statistical (systematic) errors are represented by solid (dashed) bars.Shown together is the result from PACS-CS (2011) [9] (black star).Errors for PACS-CS (2011) are statistical only.

N 2
LO of the derivative expansion for the first time and calculate the resonance parameters of the ρ meson.The mass and decay width of the ρ resonance are directly extracted from the pole position of the S-matrix, √ s pole = 886(17) − i22(9), whose real part agrees with the ρ resonance mass in the previous study but whose imaginary part leads to the coupling constant g ρππ twice as large as the previous one.Larger coupling g ρππ originates from the discrepancy in phase shifts at s 0.75 GeV 2 , whose energy region cannot be covered in the center-of-mass frame of our lattice setup.This observation provides a useful lesson for studies of P-wave resonances by the HAL QCD method and future direction for the improvement is discussed.Although the issue above remain to be verified explicitly, the result in this study shows that hadronic resonances which require all-to-all calculations can be investigated with reasonable precisions even at the N 2 LO level in the HAL QCD method.This study opens new doors toward understanding hadronic resonances by the HAL QCD method, including more challenging systems such as I = 0 ππ (σ resonance), I = 1/2 Kπ (κ resonance) and exotic resonances.

Appendix B :
Numerical evaluation for each diagramHere we outline details of numerical evaluation for each diagram calculations in this study.Figure10gives representative quark contraction diagrams appearing in the two-pion correlation functions, and the techniques utilized in the evaluations of quark propagator are shown by different colors and symbols.Other diagrams similar to these representatives are calculated similarly.In the following, we assume to employ a single noise vector for each insertion.Color and spin indices are implicit for simplicity as well.

t 0 < l a t e x i t s h a 1 _ b a s e 6 4 =
" q w d 6 + O Y u T 4 X 5 o A l n e 1 m 3 i G M w T O c = " > A A A C Z n i c h V G 7 S g N B F D 1 Z X z F q j I o o 2 A R D x C r c S E C x E m 0 s 1 R g T i C H s r p O 4 u N l d d i c B D f 6 A Y G s K K w U R 8 T N s / A E L / 0 C x V L C x8 G a z I B r U O 8 z M m T P 3 3 D k z o z m m 4 U m i x 5 D S 0 9 v X P x A e j A w N j 0 R H Y 2 P j O 5 5 d d 3 W R 0 2 3 T d g u a 6 g n T s 7 2 j b x y 1 X r P L W 8 n m H F 3 S C / u / o E e 6 4 x t Y j T f 9 a l N s n S P C H 5 D + + d z d Y G c h l c 6 k M p u Z x M p q 8 B V h z G A W 8 / z e i 1 j B O j a Q 4 3 O r O M U Z W q E n J a p M K l O d V C U U a C b w L Z T 4 J 1 q 5 i p 4 = < / l a t e x i t > t < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 o s D 1 s n B E / w 4 0 / 4 K I / I I j L C m 5 c e J s G R I t 6 h 5 k 5 c + a e O 2 d m N M c 0 P E n U C C l d 3 T 2 9 f f 0 D 4 c G h 4 Z H R y N h 4 3 r O r r i 5 y u m 3 a b k FT P W E a l s h J Q 5 q i 4 L h C r W i m 2 N I O 1 l v 7 W z X h e o Z t b c p D R 5 Q q 6 p 5 l l A 1 d l U y l 5 U 4 k S n H y Y 7 Y T J A I Q R R A p O 3 K L b e z C h o 4 q K h C w I B m b U O F x K y I B g s N c C X X m X E a G vy 9 w j D B r q 5 w l O E N l 9 o D H P V 4 V A 9 b i d a u m 5 6 t 1 P s X k 7 r J y F j F 6 o j t q 0 i P d 0 w t 9 / F q r 7 t d o e T n k W W t r h b M z e j K V f f 9 X V e F Z Y v 9 L 9 a d n i T K W f a 8 G e 3 d 8 p n U L v a 2 v H V 0 0 s y u Z W H 2 e r u m V / V 9 R g x 7 4 B l b t T b 9 J i 8 w l w v w B i Z / P 3 Q n y i / F E M p 5 M J 6 O r a 8 F X 9 G M a c 1 j g 9 1 7 C K j a Q Q o 7 P F T j F G c 5 D z 8 q Q M q F M t l O V U K C Z w L d Q Z j 4 B 8 J K J + w = = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 / Z P S I j c y q o R W I o 3 b 4 4 c u Y j s j o Q 7 2 j b x y 1 X r P L W 8 n m H F 3 S C / u / o E e 6 4 x t Y j T f 9 a l N s n S P C H 5 D + + d z d Y G c h l c 6 k M p u Z x M p q 8 B V h z G A W 8 / z e i 1 j B O j a Q 4 3 O r O M U Z W q E n J a p M K l O d V C U U a C b w L Z T 4 J 1 q 5 i p 4 = < / l a t e x i t > t < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 o s D 1 s nE N q + l b F W C J U L N p w J S C K A = " > A A A C Z H i c h V H L S s N A F D 2 N 7 / q q i i A I I p a K q 3 I r B c W V 6 M Z l H 7 Y W a p E k T j W Y J i G Z F r T 4 A 7 p V X L h S EB E / w 4 0 / 4 K I / I I j L C m 5 c e J s G R I t 6 h 5 k 5 c + a e O 2 d m N M c 0 P E n U C C l d 3 T 2 9 f f 0 D 4 c G h 4 Z H R y N h 4 3 r O r r i 5 y u m 3 a b k FT P W E a l s h J Q 5 q i 4 L h C r W i m 2 N I O 1 l v 7 W z X h e o Z t b c p D R 5 Q q 6 p 5 l l A 1 d l U y l 5 U 4 k S n H y Y 7 Y T J A I Q R R A p O 3 K L b e z C h o 4 q K h C w I B m b U O F x K y I B g s N c C X X m X E a G vy 9 w j D B r q 5 w l O E N l 9 o D H P V 4 V A 9 b i d a u m 5 6 t 1 P s X k 7 r J y F j F 6 o j t q 0 i P d 0 w t 9 / F q r 7 t d o e T n k W W t r h b M z e j K V f f 9 X V e F Z Y v 9 L 9 a d n i T K W f a 8 G e 3 d 8 p n U L v a 2 v H V 0 0 s y u Z W H 2 e r u m V / V 9 R g x 7 4 B l b t T b 9 J i 8 w l w v w B i Z / P 3 Q n y i / F E M p 5 M J 6 O r a 8 F X 9 G M a c 1 j g 9 1 7 C K j a Q Q o 7 P F T j F G c 5 D z 8 q Q M q F M t l O V U K C Z w L d Q Z j 4 B 8 J K J + w = = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 / Z P S I j c y q o R WI o 3 b 4 4 c u Y j s j o Q = " > A A A C a n i c h V H L S g M x F D 0 d X 7 W + a t 0 o b o p V c V V S E X y s i m 5 c t t V a Q U V m x l S D 8 2 I m L e r Q H 3 D n S r A r B R H x M 9 z 4 A y 7 6 C a K 7 C m 5 c e D s d E B X 1 h i Q n J / f c n C S a Y w h P M t a I K B 2 d X d 0 9 0 d 5 Y X / / A 4 F B 8 O L H h 2 R V X 5 0 X d N m x 3 U 1 M 9 b g i L F 6 W Q B t 9 0 X K 6 a m s F L 2 u F K a 7 9 U 5 a 4 n b G t d H j t 8 x 1 T 3 L V E W u i q J K v n b W j l 5 V N u N p 1 i a B ZH 8 C T I h S C G M n B 2 / w T b 2 Y E N H B S Y 4 L E j C B l R 4 1 L a Q A Y N D 3 A 5 8 4 l x C I t j n q C F G 2 g p l c c p Q i T 2 k c Z 9 W W y F r 0 b p V 0 w v U O p 1 i U H d J m c Q U e 2 S 3 r M k e 2 B 1 7 Y u + / 1 v K D G i 0 v x z R r b S 1 3 d o d O R 9 f e / l W Z N E s c f K r + 9 C x R x k L g V Z B 3 J 2 B a t 9 D b + u r J e X N t q T D l T 7 M r 9 k z + L 1 m D 3 d M N r O q r f p 3 n h T p i 9 A G Z 7 8 / 9 E 2 z M p j N z 6 c X 8 X C q 7 H H 5 F F O O Y w A y 9 9 z y y W E U O x c D d G S 5 Q j 7 w o C W V M G W + n K p F Q M 4 I v o U x + A F s E j H w = < / l a t e x i t > x + r < l a t e x i t s h a 1 _ b a s e 6 4 = " C s c H Y T w 0 4 k 2

t 0 < l a t e x i t s h a 1 _ b a s e 6 4 =
" q w d 6 + O Y u T 4 X 5 o A l n e 1 m 3 i G M w T O c = " > A A A C Z n i c h V G 7 S g N B F D 1 Z X z F q j I o o 2 A R D x C r c S E C x E m 0 s 1 R g T i C H s r p O 4 u N l d d i c B D f 6 A Y G s K K w U R 8 T N s / A E L / 0 C x V L C x 8 G a z I B r U O 8 z M m T P 3 3 D k z o z m m 4 U m i x 5 D S 0 9 v X P x A e j A w N j 0 R H Y 2 P j O 5 5 d d 3 W R 0 2 3 T d g u a 6 g n T sE R O G t I U B c c V a k 0 z R V 4 7 W G v v 5 x v C 9 Q z b 2 p a H j i j V 1 K p l V A x d l U x l Z Z n K s Q S l y I 9 4 N 0 g H I I E g N u z Y N X a x B x s 6 6 q h B w I J k b E K F x 6 2 I N A g O c y U 0 m X M Z G f 6 + w D E i r K 1 z l u A M l d k D H q u 8K g a s x e t 2 T c 9 X 6 3 y K y d 1 l Z R x J e q A b e q V 7 u q V n + v i 1 V t O v 0 f Z y y L P W 0 Q q n P H o y n X 3 / V 1 X j W W L / S / W n Z 4 k K l n y v B n t 3 f K Z 9 C 7 2 j b x y 1 X r P L W 8 n m H F 3 S C / u / o E e 6 4 x t Y j T f 9 a l N s n S P C H5 D + + d z d Y G c h l c 6 k M p u Z x M p q 8 B V h z G A W 8 / z e i 1 j B O j a Q 4 3 O r O M U Z W q E n J a p M K l O d V C U U a C b w L Z T 4 J 1 q 5 i p 4 = < / l a t e x i t > t < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 o s D 1 s n E N q + l b F W C J U L N p w J S C K A = " > A A A C Z H i c h V H L S s N A F D 2 N 7 / q q i i A I I p a K q 3 I r B c W V 6 M Z l H 7 Y W a p E k T j W Y J i G Z F r T 4 A 7 p V X L h S EB E / w 4 0 / 4 K I / I I j L C m 5 c e J s G R I t 6 h 5 k 5 c + a e O 2 d m N M c 0 P E n U C C l d 3 T 2 9 f f 0 D 4 c G h 4 Z H R y N h 4 3 r O r r i 5 y u m 3 a b k F T P W E a l s h J Q 5 q i 4 L h C r W i m 2 N I O 1 l v 7 W z X h e o Z t b c p D R 5 Q q 6 p 5 l l A 1 d l U y l 5 U 4 k S n H y Y 7 Y T J A I Q R R A p O 3 K L b e z C h o 4 q K h C w I B m b U O F x K y I B g s N c C X X m X E a G v y 9 w j D B r q 5 w l O E N l 9 o D H P V 4 V A 9 b i d a u m 5 6 t 1 P s X k 7 r J y F j F 6 o j t q 0 i P d 0 w t 9 / F q r 7 t d o e T n k W W t r h b M z e j K V f f 9 X V e F Z Y v 9 L 9 a d n i T K W f a 8 G e 3 d 8 p n U L v a 2 v H V 0 0 s y u Z W H 2 e r u m V / V 9 R g x 7 4 B l b t T b 9 J i 8 w l w v w B i Z / P 3 Q n y i / F E M p 5 M J 6 O r a 8 F X 9 G M a c 1 j g 9 1 7 C K j a Q Q o 7 P F T j F G c 5 D z 8 q Q M q F M t l O V U K C Z w L d Q Z j 4 B 8 J K J + w = = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 / Z P S I j c y q o R W I o 3 b 4 4 c u Y j s j o Q = " > A A A C a n i c h V H L S g M x F D 0 d X 7 W + a t 0 o b o p V c V V S E X y s i m 5 c t t V a Q U V m x l S D 8 2 I m L e r Q H 3 D n S r A r B R H x M 9 z 4 A y 7 6 C a K 7 C m 5 c e D s d E B X 1 h i Q n J / f c n C S a Y w h P M t a I K B 2 d X d 0 9 0 d 5 Y X / / A 4 F B 8 O L H h 2 R V X 5 0 X d N m x 3 U 1 M 9 b g i L F 6 W Q B t 9 0 X K 6 a m s F L 2 u F K a 7 9 U 5 a 4 n b G t d H j t 8 x 1 T 3 L V E W u i q J K v n b W j l 5 V N u N p 1 i a B Z H 8 C T I h S C G M n B 2 / w T b 2 Y E N H B S Y 4 L E j C B l R 4 1 L a Q A Y N D 3 A 5 8 4 l x C I t j n q C F G 2 g p l c c p Q i T 2 k c Z 9 W W y F r 0 b p V 0 w v U O p 1 i U H d J m c Q U e 2 S 3 r M k e 2 B 1 7 Y u + / 1 v K D G i 0 v x z R r b S 1 3 d o d O R 9 f e / l W Z N E s c f K r + 9 C x R x k L g V Z B 3 J 2 B a t 9 D b + u r J e X N t q T D l T 7 M r 9 k z + L 1 m D 3 d M N r O q r f p 3 n h T p i 9 A G Z 7 8 / 9 E 2 z M p j N z 6 c X 8 X C q 7 H H 5 F F O O Y w A y 9 9 z y y W E U O x c D d G S 5 Q j 7 w o C W V M G W + n K p F Q M 4 I v o U x + A F s E j H w = < / l a t e x i t > x + r < l a t e x i t s h a 1 _ b a s e 6 4 = " C s c H Y T w 0 4 k 2

(
0[r] (r + x, t).(B12)As in the case of the box diagram, we employ the CAA for x, which gives an improved triangle diagram asG tri,imp x 0 ;t 0 (r, t) = G tri,exact x 0 ;t 0 (r, t) − G tri,relaxedx 0 ;t 0 FIG. 11.A comparison in I = 1 ππ potential between two schemes at t = 14.(Left) The effective LO potentials from the ππ-type source operator.Blue (red) points show data in the point-sink (smeared-sink) scheme.(Right) Those from the NBS wave function without box diagrams.

TABLE I .
Numerical setup for the calculation.

TABLE II .
Setups for the one-end trick and the CAA in this study.N eig is the number of low eigenmodes.Color and spinor dilutions are always used.

TABLE III .
Fit parameters for effective LO potential with the ππ-type source.

TABLE IV .
Fit parameters for effective LO potential with the ρ-type source.

TABLE V .
Fit parameters of the V N 2 LO N 2 LO 2 ≤ 0 are excluded in the fit.Fit parameters for V N 2 LO 2 are summarized in Table V and the fit result is shown by a cyan band in Fig. 6 (upper right).Since significant non-smooth behavior is not observed for V N 2 LO 2 at short distances, systematic errors associated with removals of data mentioned before are not included in the analysis for V N 2 LO 2 .

TABLE VI .
Fit parameters of the Laplacian term ∇ 2 R ρ /R ρ .