Consistency between L\"uscher's finite volume method and HAL QCD method for two-baryon systems in lattice QCD

There exist two methods to study two-baryon systems in lattice QCD: the direct method which extracts eigenenergies from the plateaux of the temporal correlator and the HAL QCD method which extracts observables from the non-local potential associated with the tempo-spatial correlator. Although the two methods should give the same results theoretically, qualitatively different results have been reported. Recently, we pointed out that the separation of the ground state from the excited states is crucial to obtain sensible results in the former, while both states provide useful signals in the latter. In this paper, we identify the contribution of each state in the direct method by decomposing the two-baryon correlators into the finite-volume eigenmodes obtained from the HAL QCD method. We consider the $\Xi\Xi$ system in the $^1$S$_0$ channel at $m_\pi = 0.51$ GeV in 2+1 flavor lattice QCD using the wall and smeared quark sources. We demonstrate that the"pseudo-plateau"at early time slices (t = 1~2 fm) from the smeared source in the direct method indeed originates from the contamination of the excited states, and the true plateau with the ground state saturation is realized only at t>5~15 fm corresponding to the inverse of the lowest excitation energy. We also demonstrate that the two-baryon operator can be optimized by utilizing the finite-volume eigenmodes, so that (i) the finite-volume energy spectra from the HAL QCD method agree with those from the optimized temporal correlator and (ii) the correct spectra would be accessed in the direct method only if highly optimized operators are employed. Thus we conclude that the long-standing issue on the consistency between the L\"uscher's finite volume method and the HAL QCD method for two baryons is now resolved: They are consistent with each other quantitatively only if the excited contamination is properly removed in the former.


Introduction
The interactions between two baryons have been studied by two methods in lattice QCD. The first one is the direct method [3][4][5], which extracts the eigenenergies of the ground and/or the excited states from the temporal correlations of two-baryon systems. The binding energies and scattering phase shifts are calculated from eigenenergies using Lüscher's finite volume formula [6,7]. The second one is the HAL QCD method [8][9][10][11], which derives the energy-independent non-local kernel (called the "potential" in the literature) from the tempo-spatial correlations of two baryons. Then the binding energies and phase shifts in the infinite volume are calculated by using the Schrödinger-type equation with the kernel as the potential, which has field theoretical derivation on the basis of the reduction formula for composite operators. Both methods rely on the asymptotic behavior of the Nambu-Bethe-Salpeter (NBS) wave function, and should in principle give the same results for observables [9,11,12]. In practice, however, the current numerical results for twonucleon (N N ) systems seem to be inconsistent with each other: For heavy pion masses (m π > 0.3 GeV), both dineutron ( 1 S 0 ) and deuteron ( 3 S 1 ) are claimed to be bound in the direct method, while those are unbound in the HAL QCD method. Also, the discrepancy is ubiquitous in two-baryon systems: Although both methods indicate a bound H-dibaryon in the SU(3) flavor limit at m π = m K 0.8 GeV, the binding energy is 74.6(4.7) MeV in the direct method [13], while it is 37.8(5.2) MeV in the HAL QCD method [14]. 1 In a series of recent papers [1,2,[16][17][18][19], we have carefully examined the systematic uncertainties in both methods. The difficulty of two-baryon systems compared to a single baryon originates from the existence of elastic scattering states. Their typical excitation energies δE are one to two orders of magnitude smaller than O(Λ QCD ), so that one needs to probe large Euclidean time t (δE) −1 to extract the genuine signal of the ground state in the direct method. However, the statistical fluctuation increases exponentially in t as well as the baryon number A for multi-baryon systems as proved in [20,21]. This practically prevents one to identify the signal of the ground state in the naive analysis of the temporal correlation of two baryons.
Moreover, our extensive studies [1,2] showed that a commonly employed procedure in the direct method to identify plateaux at early time slices, t (δE) −1 , suffers from uncontrolled systematic errors from the excited state contaminations, since pseudo-plateaux 2 easily appear at early time slices. The typical symptoms of such systematics in the previous studies were explicitly exposed by the normality check 3 based on Lüscher's finite volume formula [7] and the analytic properties of the S-matrix [2].
As far as the HAL QCD method is concerned, the time-dependent formalism [10] is free from the problem of the ground state saturation, since the energy-independent potential is extracted from the spatial and temporal correlations with the information of both the ground and excited states associated with the elastic scattering. While in practical calculations there appears a systematic uncertainty associated with the truncation of the derivative expansion for the non-locality of the potential, the derivative expansion is found to be well converged at low energies [19,23]. Other systematic uncertainties such as the contaminations from the inelastic states and the finite volume effect for the potential are also shown to be well under control [19].
In this paper, we reveal the origin of the inconsistency between the direct method and 1 Recently, there appears another study [15] using the direct method, which indicates the dineutron is unbound while the H-dibaryon is bound 19(10) MeV at mπ = mK = 0.96 GeV. 2 In Refs. [1,2], they are called "fake plateaux" or "mirages" of the plateau of the ground state. 3 In Ref. [2], it is called "sanity check", a common terminology in computer science for a simple/quick test [22].
the HAL QCD method, by explicitly evaluating the magnitude of the excited states in the temporal correlation functions. We focus on the ΞΞ system in the 1 S 0 channel, which is a most convenient channel to obtain insights into the N N systems, since it belongs to the same SU(3) flavor multiplet as N N ( 1 S 0 ) but has much better statistical signals. Detailed studies in this channel were already performed with the direct method [1] as well as the HAL QCD method [19] in (2+1) flavor lattice QCD at m π = 0.51 GeV and m K = 0.62 GeV, so that the main purpose of this paper is to present an in-depth analysis by combining both results: In particular, the excited state contaminations in the temporal correlation functions are quantitatively evaluated by decomposing them in terms of the finite-volume eigenmodes of Hamiltonian with the HAL QCD potential (the HAL QCD Hamiltonian). We show how the pseudo-plateau actually appears at early time slices and also predict the time slice at which the ground state saturation is achieved. Moreover, we establish a consistency between the direct method and the HAL QCD method, by demonstrating that temporal correlation functions constructed from the optimized two-baryon operators by the eigenmode of the HAL QCD Hamiltonian show the plateaux with the values consistent with the eigenenergies at early time slices.
This paper is organized as follows. In Sec. 2, we introduce the theoretical framework of the direct method and the HAL QCD method, and present the numerical setup of the lattice QCD calculation. In Sec. 3, we recapitulate the previous analysis on the direct method [1] as well as on the HAL QCD method [19]. In Sec. 4, we decompose the correlation functions into the eigenmodes of the HAL QCD Hamiltonian. The anatomy of the excited state contaminations in the direct method is presented. We also demonstrate that eigenfunctions can be used to optimize two baryon-operators. The consistency between the temporal correlations with the optimized operators and the HAL QCD method is established. Sec. 5 is devoted to the conclusion. In Appendix A, we check how the nextto-next-leading order (N 2 LO) analysis for the HAL QCD potential affects the finite volume spectra. In Appendix B, we collect eigenfunctions of the HAL QCD Hamiltonian on various volumes. In Appendix C, we study the reconstruction of the R-correlator from the elastic states. In Appendix D, we collect the results for the reconstruction of the effective energy shifts. In Appendix E, we show the effective energy shifts from the optimized operators on various volumes. We note that a preliminary account of this study was reported in Refs. [16,18].

Methods and Lattice setup
In this section, we briefly summarize the direct method and the HAL QCD method for two-baryon systems, together with the lattice setup used in this paper.

Direct method
In the direct method for two-baryon systems, the energy eigenvalues (on a finite volume) are measured from the temporal correlation of the two-baryon operator, J sink,src BB (t), as where W n is the energy of n-th two-baryon elastic state and the ellipsis denotes the inelastic contributions. In order to obtain the energy shifts ∆E n ≡ W n − 2m B with m B being the single baryon mass, one often uses the ratio of the temporal correlation function of two-(one-) baryon system C BB (t) (C B (t)) as to reduce the statistical uncertainties as well as some systematics thanks to the correlations between C BB (t) and C B (t). The energy shift of the ground state can be obtained from the plateau value of the effective energy shift defined by with a being the lattice spacing. Here t needs to be sufficiently large compared to the inverse of the excitation energy. Once the energy shift of the ground (or excited) state on a finite volume is obtained, one can calculate the scattering phase shift in the infinite volume at that energy, δ 0 (k), via Lüscher's finite volume formula [7], where we consider the S-wave scattering for simplicity, k is defined through W n = 2 m 2 B + k 2 and L is the number of the spatial site of the lattice box. If the energy shift ∆E n is negative, the analytic continuation of the above formula to k 2 < 0 is understood. The state with a negative energy shift in the infinite volume limit corresponds to a bound state.
As noted before, the origin of the difficulty of two-baryon systems is the existence of elastic scattering states. Since the typical excitation energy of such states is (2π) 2 /((La) 2 m B ), the ground state saturation requires extremely large t, e.g., t O(4) fm at La = 4 fm and m B = 2 GeV, where the bad signal-to-noise ratio makes it practically impossible to obtain signals. In the literature of the direct method [3][4][5], however, one extracted the energy shift for the ground state from the plateau-like behavior of the effective energy shift at early time slices, t ∼ O(1) fm instead, assuming that the ground state saturation is achieved.

HAL QCD method
In the HAL QCD method, the energy-independent non-local potential U ( r, r ) is defined from with the Nambu-Bethe-Salpeter (NBS) wave function [8,9], Here |2B, W is the QCD eigenstate for two baryons with the eigenenergy W = 2 m 2 B + k 2 in the center of mass system, B( x, t) is a single baryon operator, E k = k 2 /(2µ), and H 0 = −∇ 2 /(2µ) with µ = m B /2 being the reduced mass. Eq. (2.5) has field theoretical derivation on the basis of the Nishijima-Zimmermann-Haag reduction formula for composite operators [24]. Below the inelastic threshold W th , the potential U ( r, r ) is shown to be faithful to the phase shifts, which are encoded in the behaviors of the NBS wave functions at large r.
The four-point correlation function of the two-baryon system F ( r, t) is given by where A n ≡ 2B, W n |J src BB (0)|0 is the overlap factor and the ellipsis represents the inelastic contributions.
In the time-dependent HAL QCD method [10,11], the potential is extracted directly from the so-called R-correlator as with the ellipsis being the inelastic contributions. Eq. (2.9) requires neither the ground state saturation nor the determination of individual eigenenergy W n and eigenfunction ψ Wn ( r), as all elastic states can be used to extract the energy-independent potential. Therefore, compared with the direct method, the condition required for the reliable calculation is much more relaxed in the time-dependent HAL QCD method as t O(Λ −1 QCD ) ∼ O(1) fm, where R( r, t) is saturated by the contributions from elastic states ("the elastic state saturation").
In practice, we expand the non-local potential in terms of derivatives as U ( r, r ) = n V n ( r)∇ n δ( r − r ). The leading order (LO) approximation gives U ( r, r ) V LO 0 (r)δ( r − r ), which can be determined as We can also examine the effect of higher order contributions to observables such as the scattering phase shifts: In this paper, we present the study on the correction to the LO potential for the spin-singlet channel at the next-to-next-leading order (N 2 LO) as

Lattice setup
Numerical data in previous literature [1,19] and in this paper are obtained from the (2+1)-flavor lattice QCD ensembles, generated in Ref. [25] with the Iwasaki gauge action and nonperturbatively O(a)-improved Wilson quark action at the lattice spacing a = 0.08995(40) fm (a −1 = 2.194(10) GeV). In the present paper, we make use of gauge ensembles on three lattice volumes, L 3 × T = 40 3 × 48, 48 3 × 48, and 64 3 × 64, with heavy up and down quark masses and the physical strange quark masses, corresponding to m π = 0.51 GeV, m K = 0.62 GeV, m N = 1.32 GeV and m Ξ = 1.46 GeV. We employ two different quark sources with the Coulomb gauge fixing, the wall source, q wall (t) = y q( y, t), mainly used in the HAL QCD method, and the smeared source, respectively, as in Ref. [25], and the center of the smeared source is same for all six quarks (i.e., zero displacement between two baryons), as has been employed in all previous studies in the direct method claiming the existence of the N N bound states for heavy quark masses [3][4][5]. 4 For both sources, the point-sink operator for each baryon is employed in this study. A number of configurations and other parameters are summarized in Table 1. The correlation functions are calculated by the unified contraction algorithm (UCA) [26] and the statistical errors are evaluated by the jack-knife method. For more details on the simulation setup, see Ref. [1].
In this paper, we focus on ΞΞ ( 1 S 0 ) system, which belongs to the same 27 representation as N N ( 1 S 0 ) in the flavor SU(3) transformation, but has much better signal-to-noise ratio than N N as the system contains four strange quarks. We use the relativistic interpolating operator [1] for Ξ, given by where C = γ 4 γ 2 is the charge conjugation matrix, α and a, b, c are the spinor and color indices, respectively.
3 Summary of previous studies [MeV] The source operator dependence of the effective energy shift ∆E eff (t) for ΞΞ ( 1 S 0 ) using the wall source (red circles) and the smeared source (blue squares) for L = 48 [1].
(Right) The sink operator dependence of the same quantity with the smeared source [1].
where a baryon operator B is given in Eq. (2.13). While plateau-like structures appear around t/a ∼ 15 for both wall and smeared sources, the values disagree with each other. Similar inconsistencies are found on other volumes and the infinite volume extrapolation implies that the system is bound (unbound) for the smeared (wall) source. These discrepancies indicate some uncontrolled systematic errors. Indeed, such an early-time pseudoplateau can be shown to appear even with 10% contamination of the excited states as demonstrated by the mockup data [1]. Fig. 1 (Right) shows the sink operator dependence of ∆E eff (t) for ΞΞ ( 1 S 0 ) with the smeared source fixed, where the sink operator is generalized as The last one corresponds to the simple sink operator (g(r) = 1) in Eq. (3.1). Although a plateau-like structure is observed for each sink operator, the values disagree with each other. This implies that the plateau-like behaviors at t 1 fm with the smeared source are not the plateau of the ground state but are pseudo-plateaux caused by contaminations of elastic scattering states other than the ground state. We note that such sink-operator dependence is not observed in the case of the wall source [1]. Figure 2. k cot δ 0 (k)/m π as a function of (k/m π ) 2 for N N ( 1 S 0 ) on each volume and the infinite volume in the direct method from Ref. [27] (Left) and Ref. [25] (Right). Black dashed lines correspond to Lüscher's formula for each volume, while the black solid line represents the bound-state condition, − −(k/m π ) 2 . The red line (with an error band) corresponds to the ERE obtained from the data at k 2 < 0 on finite volumes. In the left figure, the ERE fit to the data at k 2 > 0 on finite volumes together with only the infinite volume limit at k 2 < 0 is also shown by the blue line. Both figures from Ref. [2].

Normality check
Since the information on operator dependence as shown in the previous subsection is not always available, we have introduced the "normality check" in Ref. [2], based on Lüscher's finite volume formula together with the analytic properties of the S-matrix. Some examples of the normality check are given in Fig. 2, where k cot δ 0 (k) is plotted as a function of k 2 for N N ( 1 S 0 ). Red and blue lines in Fig. 2 represent fits to data by the effective range expansion (ERE) at the next-to-leading order (NLO) as where a 0 and r 0 are the scattering length and the effective range, respectively. In Fig. 2 (Left), inconsistency in ERE parameters is observed: The NLO ERE fit obtained from the data at k 2 < 0 on finite volumes (red line) disagrees with the fit to the data at k 2 > 0 on finite volumes together with the infinite volume limit at k 2 < 0 (blue line). For the latter fit (the blue line), the physical condition of the bound state pole is also violated. In Fig. 2 (Right), the NLO ERE fit exhibits a singular behavior as the divergent effective range. See Ref. [2] for more detailed discussions. As in the case of operator dependence, the normality check in Ref. [2] indicates that the plateau fitting at t 1 fm suffers from large uncontrolled systematic errors probably due to contaminations from the excited states.

Source-operator dependence and derivative expansion of the potentials
In Ref. [19], we investigated source operator dependence of the potential as a tool to estimate the systematics associated with the derivative expansion of the HAL QCD potential, Figure 3. The potential at the leading order (LO) analysis V LO 0 (r) (red circles) from the wall source (Left) and the smeared source (Right) for ΞΞ( 1 S 0 ) at t/a = 13 for L = 64 [19]. The blue squares, green triangles and black diamonds denote 1st, 2nd and 3rd terms in Eq. (2.11), respectively.
using two sources, wall and smeared sources. Fig. 3 shows the LO ΞΞ( 1 S 0 ) potential and its breakup into 1st, 2nd, and 3rd terms in Eq. (2.11) of the time dependent HAL QCD method from the wall source (Left) and the smeared source (Right). For the wall source, the 1st term dominates with moderate (negligible) contributions from the 2nd (3rd) term. As the 2nd term is not constant as a function of r, there exist small but non-negligible contributions from the excited states. For the smeared source, on the other hand, all terms are important. The substantial r dependence of the 2nd term (green triangles), which indicates large contributions from the excited states in the smeared source, is canceled by the 1st term (blue squares) and further corrected by the 3rd term (black diamonds). The total potentials (red circles) from two sources, however, show qualitatively similar behaviors, which illustrates that the time-dependent HAL QCD method works well for extracting the potential irrespective of the source types.  (r), at t/a = 13. The potential approaches to zero within errors at larger r for both sources, indicating that contributions from inelastic states are suppressed. While the potentials from two sources are very similar, there exists the non-zero difference between them. We find that this difference becomes smaller as t increases, where the potential from the wall (smeared) source is independent of (dependent on) t [19]. This indicates that the effects of higher order terms in the derivative expansion exist in the data from the smeared source.
(r) on L = 64 at t/a = 13 are plotted. As we find that V N 2 LO 0 (r) agrees well with V LO(wall) 0 (r) except at short distances, we expect that V LO(wall) 0 (r) works well to reproduce physical observables at low energies. Indeed, as shown in Fig. 5 (Right), the N 2 LO correction to the S-wave scattering phase shift δ 0 is small at low energies, showing not only that the derivative expansion converges well but also that the LO analysis for the wall source is sufficiently good at low energies.

Anatomy: Excited state contaminations in the effective energy shifts
We now show our main results of this paper, where we analyze the behaviors of ΞΞ( 1 S 0 ) temporal correlation functions for both wall and smeared sources using the HAL QCD potential and demonstrate that contaminations of elastic excited states cause pseudo-plateau behaviors at early time slices.
The strategy of our analysis is as follows. Provided that the leading order HAL QCD potential from the wall source is found to be a reasonable approximation of the exact potential, we evaluate eigenfunctions of the Hamiltonian with this potential in the finite box whose eigenvalues are below the inelastic threshold. We then calculate overlap factors between these eigenmodes and the ΞΞ ( 1 S 0 ) correlation functions, in terms of which we reconstruct pseudo-plateau behaviors of the energy shifts. We also show that the plateaux of the temporal correlation functions projected to the lowest or the 2nd lowest eigenfunction agree with their eigenvalues. This fact demonstrates that both HAL QCD potential method and accurate extraction of energy shifts from the temporal correlation function give consistent results.

Eigenvalues and eigenfunctions in the finite box
We first evaluate eigenvalues and eigenfunctions of the leading-order HAL QCD Hamiltonian in a finite box given by is a reasonable approximation of the exact potential U . 5 Note that the volume dependence of V LO(wall) 0 is negligible, thanks to the short-ranged nature of the interaction.
In a finite lattice box, eigenvalues and eigenfunctions of the Hermitian matrix H LO can be easily obtained. The n-th eigenvalue of H LO in the A + 1 representation below the inelastic threshold 6 is tabulated in Table 2, where we show the energy shift compared to the threshold, ∆E n ≡ W n − 2m Ξ . The number of excited states below the inelastic threshold is 3, 4 and 6 on L = 40, 48, and 64, respectively. For larger volume, the energy gap becomes smaller and the number of elastic excited states increases. Fig. 6 shows the eigenfunctions Ψ n ( r) on L = 48, with r |Ψ n ( r)| 2 = 1 and Ψ n ( 0) > 0. 7 Up to a normalization, Ψ n ( r) corresponds to the NBS wave function ψ W =Wn ( r) in a finite volume. The lowest state Ψ 0 ( r) has the similar shape to the R-correlator for the wall source, which shows a weak peak structure around r 1 fm and becomes flat at large 5 In Appendix A, we employ the potential at the N 2 LO analysis instead, and find that the results are consistent with the LO results. 6 For the ΞΞ system in the 1 S0 channel at mπ = 0.51 GeV, the lowest inelastic threshold is either Ξ * Ξ or ΣΩ in 5 D0 channel, which corresponds to W th − 2mΞ 0.25 − 0.30 GeV on L = 64 − 40, using mΣ = 1.40 GeV, mΞ * = 1.68 GeV and mΩ = 1.74 GeV.
7 Ψn( r) is a multi-valued function of r due to the effect of the (periodic) finite box.
distances without any nodes, while the eigenfunctions for the excited states have nodes, whose number increases as the eigenvalue becomes larger. The short distance structures for Ψ n>0 ( r), which has a steeper peak around r < 1 fm than that of Ψ 0 ( r) resemble the R-correlator for the smeared source. Eigenfunctions on other volumes are collected in Appendix B. Figure 6. Eigenfunctions of elastic states in the A + 1 representation of the HAL QCD Hamiltonian H LO below the inelastic threshold Ψ n ( r) on L = 48 at t/a = 13 for n = 0, 1, 2 (Left) and n = 3, 4 (Right), where red up-pointing triangles, blue squares, green circles, orange diamonds and black down-pointing triangles correspond to n = 0, 1, 2, 3 and 4, respectively. The eigenfunction is normalized as r |Ψ n ( r)| 2 = 1 and Ψ n ( 0) > 0. Errors are statistical only.

Decomposition of the R-correlator via eigenfunctions
Since the R-correlator is dominated by elastic states at moderately large t, we can expand it in terms of eigenfunctions of H LO as where the overlap coefficient a n characterizes the magnitude of the contribution from the corresponding eigenfunction. Using the orthogonality of Ψ n ( r), a n can be determined by a wall/smear n = r Ψ † n ( r)R wall/smear ( r, t)e ∆Ent . (4. 3) The magnitude of the corresponding excited state contamination in the R-correlator is represented by the ratio a n /a 0 . In Fig. 7, we plot a n /a 0 obtained at t/a = 13 as a function of ∆E n for the wall source (Left) and the smeared source (Right). Calculations at t/a = 14, 15, 16 confirm that the results are almost independent of t within statistical errors, indicating that the decomposition is reliable. 8 In the case of the wall source, R wall (r, t) has a large overlap with the ground state, and |a n>0 /a 0 | is smaller than 0.1. In the case of 8 In Appendix C, we check how well the decomposition in Eq. This difference of the magnitude in a n /a 0 between two sources can be understood from the overlap between the R-correlator and the eigenfunctions. Fig. 8 shows the spatial profile of the overlap, Ψ † n ( r)R( r, t)e ∆Ent , for the wall source (Left) and the smeared source (Right) on L = 48, whose spatial summation corresponds to a n . The contribution from the first excited state (blue squares) is highly suppressed for the wall source, thanks to the cancellation between positive values at short distances and negative values at long distances of Ψ † n=1 ( r)R wall ( r, t)e ∆E n=1 t . For the smeared source, on the contrary, the contamination from the 1st excited state remains non-negligible, due to the absence of the negative parts in Ψ † n=1 ( r)R smear ( r, t)e ∆E n=1 t . Figure 8. The overlap between the R-correlator and the eigenfunction, Ψ † n ( r)R( r, t)e ∆Ent , as a function of r at t/a = 13 on L = 48 for the ground state (n = 0, red circles) and the first excited state (n = 1, blue squares) in the case of the wall source (Left) and the smeared source (Right).
In the literature, the smeared source is often employed in the direct method, in order to suppress contributions from (inelastic) excited states in a "single-baryon" correlation function. The same smeared source, however, does not necessarily reduce the contaminations from the elastic excited states in the two-baryon correlation function, as is explicitly shown in Fig. 7. Indeed, one of the most relevant parameters which control the magnitudes of elastic state contributions is the relative distance r between two baryons at the source, which appears as in the center of mass system. Almost all literature of the direct method have employed the smeared source essentially corresponding to | r| = 0, which implies that elastic states for all p are equally generated at the source. Thus the choice | r| = 0 (or | r| 1) is one of the possible reasons for large contaminations from elastic excited states in the case of the smeared source. As long as | r| is non-zero and large, however, modes with non-zero p may be suppressed due to the oscillating factor e i p· r . 9 The study in which the momentum projection is performed for each baryon in the source operator is recently reported in Ref. [15].
The temporal correlation function R(t) is reconstructed in terms of eigenfunctions as where b n ≡ a n r Ψ n ( r), whose ratio b n /b 0 gives the magnitude of the contamination to R(t) from the n-th elastic excited state. Fig. 9 shows |b n /b 0 | obtained at t/a = 13 as a function of ∆E n on three volumes for the wall source (Left) and the smeared source (Right). Solid (open) symbols correspond to positive (negative) values for b n /b 0 . For the wall source, the contamination from the first excited state is found to be smaller than 1%, and |b n /b 0 | is further suppressed exponentially for higher excited states. In the case of the smeared source, the contamination from the first excited state is as large as ∼ 10% with a negative sign and the contamination remains to be ∼ 1% even for the higher excited state with ∆E n ∼ 100 MeV.

Reconstruction of the effective energy shift
Let us now examine the energy shifts obtained from the reconstructed R-correlators; Ref. [5] reported the discrepancy in the effective energies between the zero displaced (| r| = 0) and the non-zero displaced (| r| = 0) source operators, which can be naturally understood from this viewpoint, rather than the existence of two bound states claimed in Ref. [5].    where we take n max = 3, 4, 6 for L = 40, 48, 64, respectively, corresponding to the number of elastic excited states below the inelastic threshold, and b n (t 0 ) and ∆E n (t 0 ) are extracted at fixed t 0 . In Fig. 10, we show the reconstructed effective energy shifts ∆E eff (t, t 0 = 13a), together with numerical data of the effective energy shifts ∆E eff (t) from the R-correlators, for the wall source (Left) and the smeared source (Right) on L = 48. The bands correspond to ∆E eff (t, t 0 = 13a) with statistical errors coming from those of b n and ∆E n at t 0 /a = 13, while red circles or blue squares correspond to ∆E eff (t) obtained directly from the Rcorrelator in Sec. 3.1. Here we do not consider ∆E eff (t, t 0 = 13a) for t/a < 13, where inelastic contributions are expected to be larger. Shown together by the black dashed line represents the energy shift ∆E 0 (t 0 = 13a) for the ground state of the HAL QCD  Hamiltonian H LO on L = 48. We find that the results of the direct method, most notably the plateau-like structures around t/a = 15, are well reproduced by ∆E eff (t, t 0 ) for both wall and smeared sources, indicating that the behaviors of ∆E eff (t) at this time interval in the direct method are explained by the contributions from the several low-lying states. These plateau-like structures, however, do not necessarily correspond to the plateau of the ground state. Indeed, in the case of the smeared source, there is a clear discrepancy between the value of the plateau-like structure and the eigenvalue ∆E 0 of the ground state. This is a consequence of large excited state contaminations in the correlation function for the smeared source. In the case of the wall source, on the other hand, since the overlap with the ground state is large, the value of the plateau-like structure is consistent with the value ∆E 0 .
The fate of the plateau-like structures is more clearly seen in Fig. 11, where we plot the behaviors at asymptotically large t of the reconstructed effective energy shifts ∆E eff (t, t 0 = 13a) for the wall source (red band) and the smeared source (blue band). While the plateaulike structure at t/a ∼ 15 for the wall source is almost unchanged at larger t, the value of ∆E eff (t, t 0 = 13a) in the case of the smeared source gradually changes as t increases until it reaches to the value of the ground state, ∆E 0 , at t/a ∼ 100.
In Appendix D.1, we perform the same analysis on other volumes and observe essentially the same behaviors as in the case of L = 48: For the wall source, the value of the plateau-like structure at t/a ∼ 15 remains almost unchanged at larger t and is consistent with ∆E 0 . For the smeared source, the value of the plateau-like structure at t/a ∼ 15 is inconsistent with ∆E 0 . The deviation is found to be larger on a larger volume, due to severer contaminations from the excited states on larger volumes (See Sec. 4.2). 10 The value of ∆E eff (t, t 0 = 13a) for the smeared source gradually changes at t increases and the ground state saturation is realized at t/a 50, 100 and 150 or t 5, 10 and 15 fm on L = 40, 48 and 64, respectively. These time scales for the ground state saturation are actually not surprising but rather natural, considering the fact that the lowest excitation energy is as small as δE ≡ ∆E 1 − ∆E 0 84, 55 and 30 MeV on L = 40, 48 and 64, respectively.
These results clearly reveal that the plateau-like structures at t/a ∼ 15 for the smeared source are pseudo-plateaux caused by contaminations of the elastic excited states. 11 While the effective energy shifts from the wall source happen to be saturated by the ground state even at t/a ∼ 15, it is generally difficult to confirm that a plateau-like structure corresponds to the correct energy shift of the ground state without the help of other inputs, such as the HAL QCD potential analysis in the present case. Since the calculation of the energy shift from the R-correlator at t ∼ (δE) −1 is impractical due to the exponentially growing noises, one cannot obtain the correct spectra from the plateau identification in the direct method unless sophisticated variational techniques [29] are employed. 12

Projection with improved sink operator
Once the finite-volume eigenmodes of H LO with the HAL QCD potential are known, an improved two-baryon sink operator for a designated eigenstate can be constructed as which is expected to have a large overlap to the n-th elastic state. 13 This is equivalent to considering the generalized temporal correlation function with the choice of g( r) = Ψ † n ( r) in Eq. (3.2), from which we define the effective energy shift for the n-th eigenfunction as . (4.9) 10 Values of pseudo-plateaux do not strongly depend on volumes, while the correct values ∆E0 do. This is a counter example against the argument in Refs. [4,28] in which it is claimed that the volume-independence of the plateaux guarantees their correctness. 11 In Appendix D.2, we show that the dominant contamination comes from the first excited state. 12 Application of the variational method to two-baryon systems has started lately but mostly with respect to the flavor space [15]. It will be interesting to perform the variational method with respect to relative coordinate space for two baryons in addition. In the lattice study of meson-meson scatterings [30], the danger of the excited state contaminations in the plateau fitting has been already recognized and the use of the variational method is known to be mandatory. 13 Use of such an improved operator as a "source" requires additional calculations and thus is left for future study. (r) at t/a = 13, as well as that for a non-interacting system (black lines). In the case of the ground state ( Fig. 12 (Left)), the effective energy shifts from the direct method for the wall source (red circles) and the smeared source (blue squares) are also plotted for comparisons.   First of all, after the sink projections, the results with the wall source and those with the smeared source agree well around t/a ∼ 13 not only for the ground state but also for the first excited state. This is in sharp contrast with the fact that the results in the direct method without projections disagree between two sources for the ground state. Although a small overlap with the first excited state causes relatively large statistical errors in the case of the wall source, an agreement between two sources for the first excited states is rather striking, and serves as a non-trivial check for the reliability of the effective energy shifts with the sink projection. Moreover, results after the sink projections also agree with those from the HAL QCD Hamiltonian. Although the sink projection utilizes the information of the HAL QCD potential through eigenfunctions, agreements in effective energy shifts within statistical errors for both ground and first excited states provide a non-trivial consistency check between the HAL QCD method and the direct method with proper projection. In other words, results in Fig. 12 establish that (i) the HAL QCD potential correctly describes the energy shifts of two baryons in the finite volume for both ground and excited states, and that (ii) these energy shifts can also be extracted in the direct method if and only if interpolating operators are highly improved. Since the origins of systematic uncertainties are generally quite different between the two methods, such a "projection check" would be useful in future lattice QCD studies for two-baryon systems.
In recent years, it has been argued that seemingly inconsistent results for the N N systems at heavy pion masses between Lüscher's finite volume method and the HAL QCD method may indicate some theoretical deficits in one of the two methods. It is now clear from our analysis that Lüscher's method and the HAL QCD method agree quantitatively with each other, as it should be so theoretically.

Summary
In our previous works [1,2], it has been shown that the plateau fitting of the eigenenergies at early Euclidean times t, employed in the direct method, is generally unreliable for multibaryon systems, due to the appearance of pseudo-plateaux caused by contaminations of the excited states with small gap corresponding to the elastic scattering states on the finite volume. In this paper, we quantified the degree of contaminations from such excited states by decomposing the two-baryon correlation functions in terms of the finite-volume eigenmodes of the HAL QCD Hamiltonian.
By taking ΞΞ ( 1 S 0 ) system at m π = 0.51 GeV in (2+1)-flavor lattice QCD with the wall and smeared quark sources for La = 3.6, 4.3, 5.8 fm, we showed that the excited state contaminations are suppressed for the wall source, while those for the smeared source are substantial and become severer on a larger spatial extent. For the smeared source, the plateau-like structures at t = 1 ∼ 2 fm are shown to be pseudo-plateaux and the plateau with the ground state saturation is realized only at t > 5 ∼ 15 fm corresponding to the inverse of the lowest excitation energy. We also demonstrated that one can optimize the two-baryon operator utilizing the finite-volume eigenmode of the HAL QCD Hamiltonian. The effective energies from the temporal correlation functions with the optimized operators are found to be consistent with the finite volume spectra obtained from the HAL QCD Hamiltonian. This result establishes not only that the correct finite-volume spectra can be accessed by employing highly optimized operators even in the direct method but also that the HAL QCD method and the direct method agree in the finite volume spectra for the two baryon systems. Thus the long-standing issue on the consistency between Lüscher's finite volume method and the HAL QCD method is positively resolved at least for the particular system considered here. The next step is to carry out comprehensive studies of baryon-baryon interactions around the physical quark masses in the HAL QCD method, which are partly underway (see, e.g. [31][32][33]). Those will reveal not only the nature of exotic dibaryons but also the equation of state of dense baryonic matter.
A The finite volume spectra from the N 2 LO potential In the main text (Sec. 4), we study the finite volume spectra to be used for the spectral decomposition of the R-correlator using the LO potential. In this appendix, we employ the potential at the N 2 LO analysis (Eq. (2.12)) and examine a stability of the finite volume spectra against the order of the derivative expansion. The HAL QCD Hamiltonian with the N 2 LO potential is given by Note that, while H is non-Hermitian, its eigenenergies are real since the eigen equation can be rewritten as the definite generalized Hermitian eigenvalue problem [19]. The eigenenergies at L = 64 using V Table 3. The results from the different potentials are consistent with each other at low energies within statistical errors and the N 2 LO correction remains to be small (at most ∼ 1% difference) even for higher energies. These results are in line with the observation that the LO analysis from the wall source gives the correct phase shifts at low energies and the N 2 LO analysis gives only small correction even at higher energies [19].

B Eigenfunctions on various volumes
In this appendix, we show the shapes of eigenfunctions on various volume.

D Reconstructed effective energy shifts
We collect the results on the reconstructed effective energy shifts, Eq. (4.6).

D.1 The results on various volumes
As shown in Fig. 15, we observe that the reconstruction generally works well. A small deviation at early time slices on L = 64 for the smeared source is most likely due to the statistical fluctuations and/or the systematics due to contaminations from the states above the threshold discussed in Appendix C. We study how each elastic excited state contributes to the effective energy shift for the smeared source, by changing the number of elastic excited states used in the reconstruction (n max ) in Eq. (4.6). Fig. 16 shows the n max dependence of the reconstructed effective energy shifts for the smeared source at L = 40, 48 and 64, which are compared with the lowest eigenenergies ∆E 0 (red bands). Due to the negative sign of b n /b 0 < 0 (see Fig. 9 (Right)), the reconstructed effective energy shifts are smaller than ∆E 0 . For L = 40 and 48, the dominant contribution comes from the first excited state (n max = 1) besides the ground state and higher modes give only minor corrections, while the second excited state (n max = 2) also gives significant contribution for L = 64. These results indicate that the pseudo-plateau structures around t/a ∼ 15 are originated mostly from the scattering states below ∼ 90 MeV (See Table 2).