Quark mass dependence of on-shell and half off-shell scattering amplitudes from Bethe-Salpeter wave function inside the interaction range

We evaluate scattering amplitudes at on-shell and half off-shell for $I=2$ S-wave two-pion system using the Bethe-Salpeter wave function inside the interaction range in the quenched QCD. The scattering length and effective range are extracted from these scattering amplitudes. Quark mass dependence of them is investigated with the pion mass ranged in $0.52$--$0.86$~GeV. We examine consistency between a result by the conventional finite volume method and our estimate, as well as the phenomenological value.


I. INTRODUCTION
Lattice QCD has contributed to quantitative understanding of hadrons from the first principle of the strong interaction. Hadron scattering represented by the scattering length a 0 , the effective range r eff , and the scattering phase shift δ(k) itself, can be obtained by lattice QCD. The recent works are summarized in Ref. [1].
The scattering information is extracted by the finite volume method and the extensions, formulated by M.Lüscher [2,3]. An analytic function describes the relation between two-hadron energies in a finite volume and the scattering phase shift in the infinite volume. This relation is derived by use of the Bethe-Salpeter (BS) wave function outside the interaction range R of two hadrons in quantum field theory [4,5]. A method to define a potential between hadrons from the BS wave function is also proposed [6].
An associated issue between the on-shell scattering amplitude and the BS wave function inside R is argued in the infinite volume [4,5,7]. The S-wave two-pion scattering amplitude is defined by the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula under an assumption that inelastic scattering effects are negligible. In the center of mass frame, the half off-shell amplitude H(p; k) is related to the pion four-point function, = d 4 z d 4 y 1 d 4 y 2 K( p, z)K(− k 1 , y 1 )K(− k 2 , y 2 ) × 0|π 1 ( z)π 2 ( x)π 1 ( y 1 )π 2 ( y 2 )|0 , where the vectors are four-dimensional. The two-pion energy is E k = 2 m 2 π + k 2 , where k = |k| is a magnitude of three-dimensional momentum. The kernel is K( p, z) = ie i p· z (− p 2 + m 2 )/ √ Z with Z being the renormalization factor of the pion operator π i . The three momenta p, k 1 , and k 2 are on-shell. Contrarily, q is in general off-shell. The on-shell scattering amplitude H(k; k) is associated with the scattering phase shift δ(k) through H(k; k) = 4π k e iδ(k) sin δ(k).
Some ratio is needed to cancel out e iδ(k) , which can not be measured directly on the lattice. Based on this issue, we accomplished a simulation to obtain scattering amplitudes from the BS wave function inside R on a finite volume lattice [8]. The simulation utilized the isospin I = 2 S-wave two-pion in the quenched QCD with the lattice spacing of a −1 = 1.207 GeV at the pion mass m π = 0.86 GeV. Using the on-shell amplitude, we numerically confirmed agreement between a result by the finite volume method and that by our approach. We also presented lattice QCD can successfully give the half off-shell scattering amplitude. Although the half off-shell scattering amplitude is not an observable in experiments, it gives constraints on parameters of effective theories and models of the hadron interaction. Lattice QCD results of the half off-shell scattering amplitude can be useful as supplementary data to effective theories and models.
In this paper, we extend our previous simulation [8] to investigate the scattering amplitudes with several pion masses m π = 0.52-0.86 GeV at the same lattice spacing. The results are extrapolated to the physical point to examine consistency among the previous result by the conventional finite size method, the phenomenological estimate, and our result. This paper is organized as follows. Sec. II devotes the formulation in the lattice theory. In Sec. III, details of our simulation set up are presented. Sec. IV explains our results at each m π , as well as those at the physical point. Comparison with the previous lattice QCD results and the phenomenological values is also given in this section. Sec. V summarizes this paper. Appendix contains the operator dependence of the scattering amplitude, formulation of the scattering amplitude using the BS wave function in the momentum space, and discussion on time dependence of the scattering amplitude.

II. FORMULATION ON THE LATTICE
Formulation of the scattering amplitude is explained. We use the same notation in our previous work [8], following Refs. [4,5,7]. We restrict ourselves to the S-wave scattering of I = 2 two-pion in the center of mass frame, assuming inelastic scattering effects are negligible. The two-pion BS wave function on the lattice φ(x; k) is obtained by where |π + π + , E k is a ground state of two pions, and E k = 2 m 2 π + k 2 is the two-pion energy. Φ(x, t) is an operator of two pions, is performed to obtain the S-wave scattering state with an assumption that higher angular momentum contributions of l ≥ 4 scattering are negligible. The pion interpolating operator is defined by where Ω(t src ) is an operator of two-pion at the source time slice t src , and t = |t sink − t src |. The dot term corresponds to excited state contributions. C k is an overall constant. An essential quantity to calculate the scattering amplitude from φ(x; k) is the reduced BS wave function h(x; k) [5,7]. It is defined through φ(x; k) as where h(x; k) possesses an important property that h(x; k) equals to zero outside the interaction range of two pions R, except for the exponential tail, h(x; k) defines the half off-shell amplitude on the lattice H L (p; k) with an off-shell momentum p, where j 0 (px) is the spherical Bessel function. If R is less than half of the lattice extent L, R < L/2, and the exponential tail is negligible, then the range of the summation can be changed from L to ∞ due to Eq. (10), as discussed in Ref. [8]. H L (p; k) becomes the half off-shell amplitude in the infinite volume H(p; k) in Eq. (1) as except for an overall finite volume correction, C 00 = C k /F (k, L). F (k, L) is a finite volume correction of the two-pion state, called the Lellouch and Lüscher factor [9]. As explained above, R < L/2 with a negligible exponential tail of h(x; k) is the sufficient condition for the scattering amplitude calculation on the lattice. We numerically confirm this condition is satisfied, in the later section.

III. SET UP
Our simulation set up is presented. We perform a quenched QCD simulation following Refs. [5,10]. We generate 200 gauge configurations on 24 3 × 96 lattice by Hybrid Monte Carlo algorithm, stored at every 100 trajectories. Our gauge action is Iwasaki-type [11]. The lattice spacing is a −1 = 1.207 GeV at the bare coupling is β = 2.334. Our valence quark action is Clovertype [12]. The Clover coefficient is mean field improved, C SW = 1.398 [10]. Our pion masses are in range of m π = 0.86-0.52 GeV with the valence quark hopping parameters of κ val = 0.1340-0.1369. Our simulation parameters are collected in Table I. We also generate gauge configurations on 24 3 ×64 with the same setup to investigate the source operator dependence of the scattering amplitude on the lattice. Details of the investigation are explained in Appendix A.
The two-pion four-point function C ππ (x, t) in Eq. (6) is calculated using complex random Z2 sources to avoid Fierz contamination, where Ω(t) is the source operator in Eq. (6). π + (t, η) is defined by (14) N r is the number of the random sources η i (x), satisfying the following condition, N r = 4 is employed in our simulation. The source resides at a time t src in Eq. (6) and all spatial points, as well as all colors and spins. The latter reduces the simulation cost [13]. In our set up, we found the gain is roughly a factor of three. We also calculate a single pion correlator C π (t) with the same source to measure m π in the large time region, t ≫ 1, We perform the measurements with every four time slice per configuration, i.e., the total number of t src is 24. We adopted the periodic boundary condition in space and the Dirichlet boundary condition in time. Distance between the Dirichlet boundary and t src is kept to be 12.
We employ two methods to determine the interaction momentum k. One is a momentum from E k , denoted by k t , E k is obtained from the temporal correlator of two pions using C ππ (x, t) in Eq. (6), The other is a momentum from the BS wave function φ(x; k) outside the interaction range, denoted by k s , The condition Eq. (10) is convinced by definition.

IV. RESULT
A. Effective mass and energy Figure 1 represents our results of effective masses of a single pion, defined by A plateau of the effective mass starts from t = 14 in all κ val cases. We determine m π from a single exponential fit to C π (t) in the range of [t min , t max ] = [14,74]. The values of m π are listed in Table II.
We employ the same temporal fitting procedure as that of a single pion to determine an interaction energy of two pions. The fitted values of E k is summarized in Table II.
Using the results for m π and E k obtained from the fits, k 2 t in Eq. (17) is evaluated. The values of k 2 t are tabulated in Table III.

B. BS wave function
The BS wave function φ(x; k) in Eq. (3) is calculated as follows. An effective BS wave function is defined by We choose the reference position of x ref = (12,7,2). φ eff (t, x) and k 2 t tabulated in Table III. Figure 4 illustrates our results of h(x; k). We employ a ratio of h(x; k) over φ(x ref ; k) at a reference point x ref to cancel out the overall factor. For x 10, h(x; k t ) = 0 is found to be satisfied in our statistical precision. This result shows the interaction range R ∼ 10 < L/2, and the exponential tail of h(x; k) is negligible compared to our statistical error. Our data guarantees the sufficient condition of Eq. (11) in our quark mass region.
Using data in the outside region of R, an alternative interacting momentum k 2 s in Eq. (19) can be determined, which is more precise than k 2 t [5]. We obtain k 2 s from a constant fit to −∆φ Table III collects our results of k 2 s , as well as k 2 t . k 2 s is consistent with k 2 t with a smaller error than that of k 2 t by a factor of two.

D. Scattering amplitude
Once the sufficient condition R ∼ 10 < L/2 is satisfied, the scattering amplitude can be computed using Eq. (11). We choose k = k s in the following analyses, unless explicitly stated. Figure 5 displays off-shell momentum dependence of the half off-shell amplitude H(p; k) in Eq. (12). The overall factor of H L (p; k) in Eq. (11) are eliminated by taking a ratio of H L (p; k) over its on-shell value H L (k; k), A clean signal of the ratio is observed throughout our p 2 range. The validity of H(p; k) is ensured below the threshold drawn in the figure at p 2 = 3m 2 π i.e. E k = 4m π , though the quenched approximation prohibits dynamical inelasticity.
The operator dependence of H L (p; k) is examined for both source and sink operators. The dependence is under control in our simulations. The details are explained in Appendix A.
We discuss the lattice artifacts in our result of H L (p; k). The rotational symmetry breaking at the finite lattice spacing causes deviation between on-axis and off-axis h(x; k) values. The influence to H L (p; k) is evaluated to be 3 % in our simulation at a −1 = 1.207 GeV. The size of the error is comparable to our statistical error. Another issue is the finite lattice artifact in the short distance. It highly affects the data especially around x = 0. They are suppressed in H L (p; k), however, due to the Jacobian factor r 2 in the integration, Contribution from h(x; k) near x = 0 is not significant. Nevertheless, the continuum limit is required to remove these lattice artifacts. H L (p; k) can be also calculated with the BS wave function in the momentum space, φ(p; k) = x φ(x; k)e −ip·x (see for example Ref. [14]). We numerically confirmed both approaches give a consistent value. Appendix B explains details of the formulation in the momentum space.
Two supplemental confirmations of the validity of R ∼ 10 are possible by use of the scattering amplitude. One is the range of the summation for H L (p; k) in Eq. (11). H L (p; k) using a summation over all spatial volume gives a consistent value with that using a summation up to x = 10 ∼ R, which implies correctness of the estimate of R ∼ 10. The other utilizes the on-shell scattering amplitude H L (k; k) and an analytic solution of (∆ + k 2 )φ(x; k) = 0 in x > R. In the S-wave case, the analytic solution φ x>R (x; k) can be expressed using the Green function on the lattice G(x; k), where v 00 is a constant. φ x>R (x; k) can be also expressed by the phase shift δ(k), where (l ≥ 4) contains only the spherical Bessel functions j l (px) of l ≥ 4. Comparing Eq. (26) with Eq. (29) using the expansion by j 0 (kx) and l = 0 spherical Neumann function n 0 (kx) leads to two simple equations [5]. The coefficient of n 0 (kx) gives In parallel, the coefficient of j 0 (kx) provides where Using Eqs. (30) and (31), one obtains the finite volume formula [3], Therefore, Eq. (30) must be satisfied in the finite volume method. Based on this argument, we define an indicator R(x; k) to test the equality in Eq. (30), Outside the interaction range R, R(x; k) becomes unity, if Eq. (30) is satisfied. Figure 6 represents our result of R(x; k). It increases monotonically and approaches to unity. R(x; k) is consistent with unity in x 10, as expected. It validates our estimate of R ∼ 10.
We also discuss t dependence of the on-shell amplitude H L (k; k). H L (k; k) is calculated at each t for the t dependence analysis, in contrast to the above calculation using the fit result in Sec. IV B. Figure 7 illustrates t dependence of H L (k; k)/(C k φ(x ref ; k)). The ratio is almost flat in t. .
We used the expansion of C k φ(x ref ; k) in Eq. (29) and assumed l ≥ 4 terms are negligible. The phase factor as well as the overall constants are canceled in the ratio. Inversely, δ(k) is given by .
The reference point is chosen to be x ref = (12,7, 2) by the following procedure. Evaluation of tan δ(k) through x > R, is an alternative expression of Eq. (29). Assuming A l (k), B l (k) = O(1), the size of the l = 4 contribution at each x is estimated by using a ratio Y (x; k), Y (x; k) with k = k t at κ val = 0.1340 is presented in Fig. 9. The values of Y (x; k) at some positions in x > 10 are found to be close to zero, satisfying |Y (x; k)| < 10 −6 . Similar results of Y (x; k) are obtained in other κ val . From the estimation we choose a reference point as x ref = (12, 7, 2). The effective range expansion defines the scattering length a 0 and the effective range r eff , We estimate a 0 by The tiny value of k 2 presented in Table III justifies the estimation. Similarly, r eff can be evaluated by where we assume

Scattering length
We evaluate a 0 /m π through Eq. (40) using tan δ(k) obtained from Eq. (36). In the evaluation, not only k = k s but also k = k t is employed. A smaller error of a 0 /m π is obtained by k = k s . The results are tabulated in Ta x ref = (12,12,12) gives the largest value of |Y (x; k)| as presented in Fig. 9. The result clearly indicates a sizable l = 4 contribution at x ref = (12,12,12). Our analysis of x ref dependence of a 0 /m π suggests the l = 4 contribution in a 0 /m π is suppressed well by our choice of x ref = (12, 7, 2). a 0 /m π is extrapolated to the physical m π using a formula motivated by chiral perturbation theory [16], where A a0 , B a0 , C a0 are fitting parameters. Our results at the physical point are listed in Table IV. Figure 11 summaries chiral extrapolations of a 0 /m π . Our value at the physical point is consistent with the previous result by lattice QCD using the conventional finite volume method based on φ(x; k) outside the interaction range [5] and the phenomenological estimate [15]. The agreement ensures our approach with φ(x; k) inside the interaction range.

Effective range
The effective range r eff is evaluated by the slope of H L (p; k) with respective to p 2 and δ(k) shown in Eq. (41). Our results with k = k s and k = k t are compiled in Table V. The result using k t has a larger error than that using k s , as in the case of a 0 . Figure 12 plots m π dependence of m π r eff . Our result of m π r eff agrees with the  value by the finite volume method, calculated using data in Ref. [5]. Our data are more accurate than those of the finite volume method due to the explicit p 2 dependence of H(p; k). r eff is extrapolated to the physical point using a formula based on chiral perturbation theory [18].
where A r eff , B r eff are fitting parameters. Our result at the physical point, whose value is summarized in Table V, underestimates the phenomenological value [15]. The reason seems to be the chiral extrapolation of m π r eff . It rapidly grows toward the physical point. We also need to validate the two assumptions in Eq. (41), though consistency between our result and that of the conventional finite volume method is confirmed at each simulation point. Another possibility is the quenching effect. In fact, N f = 2 + 1 lattice QCD using m π = 390 MeV successfully reproduces the phenomenological estimate [17].
More realistic N f = 2 + 1 data around the physical point are required to draw a definite conclusion.

V. SUMMARY
We have successfully calculated the on-shell and half off-shell scattering amplitudes from the BS wave function inside the interaction range, as reported in our previous paper [8]. Our approach utilizes the BS wave function inside the interaction range, while the conventional finite volume method is based on the BS wave function outside the interaction range. The on-shell scattering amplitude gives direct scattering information through the phase shift. The half off-shell amplitude is not an observable in experiments, on the other hand, but it could be an important input of theoretical effective theories and models to constrain their parameters. Furthermore, the half off-shell amplitude gives the effective range under two assumptions.
In this article, we extended our study to investigate quark mass dependence of I = 2 S-wave two-pion onshell and half off-shell scattering amplitudes at the center of mass in the quenched QCD. Our simulation was performed at the lattice spacing of a −1 = 1.207 GeV using pion masses of m π = 0.52 − 0.86 GeV on 24 3 × 96 lattice. We first checked the interaction range is within half of our spatial lattice, which satisfies a sufficient condition of our method as well as the finite volume method. It allows us to evaluate on-shell and half off-shell scattering amplitudes. We obtained clean signals of them. We also discussed the source and sink operator independence of our scattering amplitudes.   11. Quark mass dependence of the ratio of the scattering length over the pion mass a0/mπ. Open circles are lattice QCD results with the conventional finite volume method [5]. Open square is a phenomenological estimate of Chiral Perturbation Theory(ChPT) [15]. Open symbols are shifted to some extent for clarification of data.
We then extracted the scattering length from the onshell scattering amplitude. Our results at each m π and the physical point agree with those obtained by the finite volume method. It proves our approach is an alternative to the conventional finite volume method. Open circles are N f = 0 lattice QCD results with the conventional finite volume method calculated using data in Ref. [5].
Open triangle is N f = 2 + 1 lattice QCD result with the finite volume method [17]. Open square is a phenomenological estimate of Chiral Perturbation Theory(ChPT) [15]. Open symbols are shifted to some extent for clarification of data.
We also extracted the effective range from the slope of the half off-shell amplitude at the on-shell momentum under two assumptions to be validated. Our result agrees with that by the conventional finite volume method at

ACKNOWLEDGMENTS
We thank J. Carbonell and V. A. Karmnov for pointing out the momentum space formulation, and members of PACS collaboration for useful discussion. We also appreciate INT at University of Washington for the kind hospitality and partial financial support during an early stage of this work. Our simulation was performed on COMA under Interdisciplinary Computational Science Program of Center for Computational Sciences, University of Tsukuba. This work is based on Bridge++ code (http://bridge.kek.jp/Lattice-code/) [19]. This work is supported in part by JSPS KAKENHI Grant Numbers 16H06002 and 18K03638.
The source operator dependence is simply explained by an overall factor of φ(x; k). When the ground state dominates C ππ (t) in large t region, the source opera- In order to check the source operator independence of h(x; k)/φ(x ref ; k), we compare results using the Z2 random and wall sources on the 24 3 × 64 ensemble with κ val = 0.1340, which was used partly in our previous paper [8]. The same simulation set up is adopted as that in Sec. III. We use 400(200) configurations with 32(16) measurements per configuration using the random Z2(wall) source. The simulation parameters are listed in Table VI. The wall source operator Ω wall (t) in C ππ (x, t) of Eq. (6) is defined by Ω wall (t) = π + wall (t src )π + wall (t src + 1), The pion operators in Ω wall (t) are placed at t src and t src + 1 to prevent from Fierz rearrangement [20,21]. Figure 13 exhibits the effective mass m eff π and energy E eff k , defined in Eqs. (20) and (21) respectively. Our results from the two sources are consistent with each other, as well as that on 24 3 ×96 lattice. We fit data in the range of [t min , t max ] = [14,44] to extract m π and E k . The results are summarized in Table VII together with those for k 2 t and k 2 s . In contrast to the case on 24 3 × 96 lattice in Fig. 3, the data of φ(x; k) on 24 3 × 64 lattice available for analysis of the BS wave function are limited to those at a time slice of t = 44. We use φ(x; k) at t = 44 for our analysis on 24 3 × 64 lattice. Figure 14 presents the two source results of h(x; k)/φ(x ref ; k). They agree with each other, and with the result obtained on 24 3 × 96 lattice in Fig. 4. Consistency of the two results proves the source independence of h(x; k s )/φ(x ref ; k). We also confirmed the source operator independence of the physical quantities a 0 /m π and r eff extracted from H L (p; k). The results are in Table VII. Our analysis establishes the source operator dependence is well under control.

Sink operator dependence
The sink smearing of the pion operator in C ππ (t) produces an extra overall factor. It can not be removed by the same ratio as that in the source operator case. Using a smearing function of the pion f (x), the sink smeared BS wave functionφ(x; k) is obtained bȳ where only one of the pion operators is smeared.
is not the same as φ(x; k)/φ(x ref ; k) in contrast to the source operator case [7]. Nevertheless, the extra sink smearing factor can be analytically removed [22]. The sink smeared scattering amplitudeH L (p; k) defined byφ(x; k) relates to H L (p; k) such that where The local operator corresponds to f (x) = δ (3) (x) with C f (p) = 1 in all p. Once f (x) is given, C f (p) can be analytically calculated and removed fromH L (p; k). We numerically check the sink smearing independence of our result.φ(x; k) andH L (p; k) are calculated by replacing the integration and e −ip·x in Eqs. (A3), (A4), and (A5) to the summation and j 0 (px). We employ an exponential sink smearing function f (x) = e −Ax with a constant A. We confirm the values ofH L (p; k)/C f (p) with different f (x) are consistent, unless f (x) is too broad comparing to L/2 − R. A broad sink smearing function is found to lift up h(x; k) in |x| > R, which violates the sufficient condition Eq. (10). same on-shell scattering amplitude. In the zero momentum limit, Eq. (B3) leads to the scattering length a 0 , as in Ref. [14]. Equation (B3) is not suitable, however, for the lattice calculation on a finite volume. φ ∞ (x; k) at all x in the infinite volume is demanded.
In the following, we consider a formulation with the finite integration range [8]. If h ∞ (x; k) = 0 outside the interaction range R, then the integration range of H(p; k) can be changed from ∞ to R, The partial integration gives where the second term is the surface term. At on-shell p = k, the first term vanishes. H(p; k) is expressed by the surface term only.
In the spherical coordinate, H(p; k) can be simplified.

Finite volume
The lattice version of Eq. (B1) is where C k is an overall constant in Eq. (11). h(x; k) is the reduced wave function on the lattice. A choice of L max = L/2 and L min = L/2 − 1 corresponds to the summation over the entire spatial volume with its extent L. L max and L min can be decreased, as long as L max , L min > R is satisfied.
The partial integration on the lattice leads to The surface term on the lattice surf(p; k) is given by where X(a) = x except for X i (a) = a. The surface term is not zero, in general. If L max = L/2 and L min = L/2 − 1 are chosen and p i L max = n i π, n i ∈ Z is satisfied under the periodic boundary condition, the surface term becomes zero.
In the case of the S-wave scattering on the lattice, H L (p; k) becomes where ∆ is the symmetric lattice laplacian defined in Eq. (9). The surface term in this case is given by where X(a) = |X(a)|. If L max = L/2 and L min = L/2−1 are chosen, the periodicity and isotropy of φ(x; k) leads to a simpler form of surf(p; k), where X ′ 1,2 (a) = x 1,2 and X ′ 3 (a) = a.
The t dependence of H L (k; k)/(C k φ(x ref ; k)) in Fig. 7 is discussed under several assumptions.
We define H L (t, k; k) to study the scattering amplitude H L (k; k) at each t, We evaluate a ratio of H L (t, k; k)/C ππ (x, t). It can be split into the ground state and the excited state parts, where the excited state contributions δH L (t, k; k) and δC ππ (x ref , t) for the numerator and denominator, respectively. Figure 7 illustrates H L (t, k; k)/C ππ (x ref , t) is almost flat against t. On the other hand, Fig. 8 clearly reveals non-negligible contributions of δH L (t, k; k) and δC ππ (x ref , t) in the small t region. It suggests a possibility of cancellation between the numerator and the denominator, This is a sufficient condition of flat t dependence of H L (t, k; k)/C ππ (x ref , t).
We demonstrate the sufficient condition is realized, if C ππ (x ref , t) is dominated by scattering states with almost zero momentum. In the small t region, C ππ (x, t) includes large contributions from not only the ground state of two pions, but also scattering states with the first radial excited state of π, denoting π ′ . We restrict our consideration below the energy of π ′ π ′ → π ′ π ′ scattering, neglect-ing inelasticities. Then, C ππ (x, t) is expressed as where A q (t) = C q e −Eqt , A ′ q ′ (t) = C ′ q ′ e −E ′ q ′ t and A ′′ q ′′ (t) = C ′′ q ′′ e −E ′′ q ′′ t with C q , C ′ q ′ and C ′′ q ′′ are overall constants for each contribution, and E ′′ q = 2 m 2 π ′ + q 2 .
The excited state contributions δH L (t, k; k) and δC ππ (x ref , t) for the numerator and denominator are expressed as The cancellation condition of Eq. (C3) implies the coefficients of each state contribution in δH L (t, k; k) and δC ππ (x ref , t) coincide. At some momentum p, the relation between the coefficients is We investigate p 2 dependence of R surf (k; p) in the range of 0.1k 2 ≤ p 2 ≤ 10k 2 , supposing unmeasured q ′ 0 and q ′′ 0 are in this range. The result is plotted in Fig. 15. Difference of R surf (k; p) is less than 3 % even at p 2 = 10k 2 . Our data support validity of Eq. (C12).
In summary, the excited states contamination in the ratio of H L (k; k)/(C k φ(x ref ; k)) is understood under the following conditions.
• Contribution from higher momentum states is small, due to our choice of the source operator.
Then, contamination between H L (t, k; k) and C ππ (x ref , t) can be explained. For further analysis, we need the variational method to distinguish the excited states.