Mirage in Temporal Correlation functions for Baryon-Baryon Interactions in Lattice QCD

Single state saturation of the temporal correlation function is a key condition to extract physical observables such as energies and matrix elements of hadrons from lattice QCD simulations. A method commonly employed to check the saturation is to seek for a plateau of the observables for large Euclidean time. Identifying the plateau in the cases having nearby states, however, is non-trivial and one may even be misled by a fake plateau. Such a situation takes place typically for the system with two or more baryons. In this study, we demonstrate explicitly the danger from a possible fake plateau in the temporal correlation functions mainly for two baryons ($\Xi\Xi$ and $NN$), and three and four baryons ($^3{\rm He}$ and $^4{\rm He})$ as well, employing (2+1)-flavor lattice QCD at $m_{\pi}=0.51$ GeV on four lattice volumes with $L=$ 2.9, 3.6, 4.3 and 5.8 fm. Caution is given for drawing conclusion on the bound $NN$, $3N$ and $4N$ systems only based on the temporal correlation functions.


Introduction
In lattice QCD, observables such as the energies and the matrix elements of hadrons are commonly extracted from temporal correlation functions at large Euclidean time where ground state saturation is expected to be realized. For example, a two point correlation function C(t) for the operator O 1,2 (t, x) is related to physical quantities as where |n is the n-th one-particle (zero momentum) eigenstate of QCD with mass m n which couples to the operator O 1,2 , and Z n is the corresponding pole residue, Z n = 0|O 1 (0, 0)|n n|O 2 (0, 0)|0 . The ellipsis represents contributions from two or more particle states. Assuming the ordering that 0 < m 1 < m 2 < m 3 · · · , we can extract the mass and the matrix element for the lowest energy state from the large t behavior of C(t) as where contributions from two or more particle states are suppressed for t → ∞.
In practice, we take large but finite t, so that e −(m 2 −m 1 )t becomes negligibly small. If m 2 − m 1 = O(Λ QCD ), which is generally true for single hadron states in QCD, it requires t ≥ O(1) fm. Therefore, one can safely extract single-hadron masses as long as C(t) is accurate enough at t ∼ O(1) fm. To check whether C(t) is dominated by the ground state, the effective mass, defined by is often employed, where a is the lattice spacing. If m eff (t) becomes almost independent of t at t ≥ t min ("the plateau"), C(t) is considered to be dominated by the ground state and the mass is extracted from C(t) by using the data at t ≥ t min .
For multi-hadrons, the energy shift of the whole system on the lattice relative to the threshold defined by the sum of each hadron masses is of interest, since it has information on the binding energy and the scattering phase shift [1]. For the energy shift of the twobaryon system, ∆E BB ≡ E BB − 2m B , where E BB is the lowest energy of the two-baryon system and m B is the baryon mass, one introduces the effective energy shift defined by where R BB is the two-baryon propagator C BB (t) divided by the one-baryon propagator C B (t) squared as (1.6) and the effective energy of two-baryon system E eff BB (t), which is defined by (1.7) In actual numerical simulations, it is often observed that the statistical error for ∆E eff BB (t) is substantially reduced from the individual errors for E eff BB (t) and 2m eff B (t) due to their mutual correlations. In addition, ∆E eff BB (t) shows a plateau-like behavior at relatively earlier time t than it is supposed to be, so that one may be tempted to extract physical information from such a behavior.
In this paper, we address the issue whether the plateau-like behavior observed for the effective energy shift of the multi-baryons system is reliable or not. Indeed, it was previously claimed, by fitting the plateau-like behavior of the effective energy shifts, that dineutron, deuteron, 3 He and 4 He are all bound for heavy pion masses, m π 510 MeV [2] and m π 300 MeV [3]. For making detailed comparisons with such previous results, we employ the same lattice setup as Ref. [2]. We perform more measurements of baryon correlation functions than those in previous studies to investigate the reliability of the plateau-like behavior from the point of view of statistics, while we take two different source operators (the smeared source used in [2] and the wall source 1 ) as well as two different single-baryon operators (the non-relativistic type used in [2] and the relativistic one) to study the reliability of the plateau-like behavior from the point of view of systematics. This paper is organized as follows. In Sec. 2, we give general considerations on the plateau identification in multi-baryon system, and explicitly demonstrate the danger of the fake plateau using the mock-up data. In Sec. 3, lattice simulation parameters used in this paper are summarized. In this paper, we consider the effective energy shift for ΞΞ as well as N N, 3N, 4N systems. In Sec. 4, we study the ΞΞ systems in detail, since signal to noise ratio (S/N ) in lattice QCD is better for Ξ than N . This is due to the fact that Ξ contains two heavier strange quarks, while N consists of lighter up and down quarks only. As demonstrated in Sec. 2, we observe plateau-like behaviors in the effective energy shift around t ∼ 1 fm, which however disagree between the smeared source and the wall source. We then discuss that it is difficult to judge which plateau (or neither) is true only from the information of time correlation functions. In Sec. 5 and 6, we analyze the N N, 3N, 4N systems in a similar manner. Although statistical errors are larger, we observe similar disagreements between two sources as in the case of ΞΞ systems. In Sec. 7, conclusions in this paper are given and some discussions follow. In appendix A, we present the study on the sink operator dependence for effective energy shifts. Disagreements are observed among plateau values from different sink operators for the smeared source, but not for the wall source. In appendix B, we collect the figures for effective energy shifts on various volumes.

Difficulties in multi-baryon systems
Even though the plateau method works in principle, and indeed in practice for, e.g., the ground state meson masses, the method sometimes suffers from difficulties, in particular, in the case of multi-baryon systems.
First of all, we note that the requirement of the ground state dominance encounters a fundamentally new challenge when one studies multi-hadron systems instead of singlehadron systems. In fact, t min required for the ground state dominance becomes much larger for multi-hadron states, since δE ≡ E 2 − E 1 is much smaller where E 1 is the ground state energy while E 2 is the lowest excited state energy. For example, in the case of bound states, δE is a few MeV for deuteron and a few tens of MeV for 4 He. With the absence of bound states as is the case for dineutron or diproton, there exist only continuum states and thus 1 The wall source has been adopted in the HAL QCD method [4][5][6][7][8], which utilizes the space-time correlation functions instead of just the temporal correlation to study multi-hadrons. In this method, bound states for dineutron and deuteron are not found at similar values of pion masses [9][10][11]. Detailed comparison between the HAL QCD approach and the approach discussed in this paper by using the same lattice data will be given in independent publications under preparation and will not be discussed in the present paper. no intrinsic energy gap exists. In lattice calculations, the energy spectrum is discretized in a finite box with the spatial extension L, leading to δE (2π) 2 /(L 2 m N ), which is also small as, for instance, δE 25 MeV at large enough L 8 fm for two baryons at the physical quark masses. These small splittings are in sharp contrast to the single-hadron systems, where δE ∼ O(Λ QCD ).
The requirement of taking large t causes a serious difficulty in lattice QCD, since the data at larger t are in general accompanied with much worse S/N . The situation is severe in particular for multi-baryon systems at large t, for which we have [12,13] where m B and m M are the ground state baryon mass and the meson mass coupled to the BB annihilation channel, respectively. The signal S A (t) is given by a propagator for an A-baryon system, schematically denoted as with (zero momentum) baryon creation and annihilation operatorsB(t) and B(t), while the noise N A (t) is given by The asymptotic formula Eq. (2.1) says that S/N becomes worse for bigger t as well as larger numbers of baryons and/or smaller quark mass (i.e. lighter meson). This may prevent us from taking sufficiently large t to guarantee the t independence of E eff A (t), so that we can not reliably control systematic errors from excited state contaminations.
In order to demonstrate the danger of such excited state contaminations, we consider the mock-up data given by where ∆E BB = E BB − 2m B with the ground state energy E BB , while δE el = E * BB − E BB and δE inel = E inel − E BB with the first excited elastic state energy E * BB and the lowest inelastic state energy E inel , respectively. Thus the effective energy shift becomes so that ∆E eff BB (t) − ∆E BB corresponds to the deviation in ∆E eff BB (t) from its true value. Note that both b 2 /b 1 and c 1 /b 1 can be negative if source and sink operators are different. As an example, we consider δE el = 50 MeV, which is the typical lowest excitation energy of elastic two-baryon scattering states in our numerical setup with La = 4.3 fm lattice (see Sec. 3), while we take δE inel = 500 MeV, which is roughly the order of m π in our simulations. In lattice QCD, one often tries to tune the interpolating operator so that excited state contaminations are suppressed. Since the difference between inelastic states and the ground state is expected to be intrinsic in QCD, one may ideally take a good operator for baryons which have small overlaps with inelastic states. We therefore adopt a small value c 1 /b 1 = 0.01 as the contamination from the inelastic state. On the other hand, it is much more difficult to separate the ground state from the elastic excited state by tuning the operator, since the difference between these states do not originate from QCD, but from the use of a finite lattice box. Accordingly, we take b 2 /b 1 = ±0.1 as well as b 2 /b 1 = 0 for a comparison, as the contamination of the excited elastic state.
In Fig. 1 (Left), we plot ∆E eff BB (t) − ∆E BB as a function of t for the above choice of parameters. Let us consider the case with b 2 /b 1 = 0 (black line) first. In the absence of the excited elastic state, the effective energy shift ∆E eff BB (t) smoothly approaches to the plateau (from above for the positive c 1 /b 1 ) and t 1 fm is sufficient to reduce the systematic error from the contamination to the level of accuracy we need for ∆E BB . Unfortunately, this ideal situation cannot be realized in practice, and for a more realistic cases with ±10 % contamination of the 1st excited elastic state at t = 0, we need t 8 − 10 fm to achieve the level of accuracy we need, as shown by red and blue lines. In practice, however, t min 8 − 10 fm is too large to have a good signal due to the exponentially decreasing signal to noise ratio for multi-baryons as mentioned before. Shown in Fig. 1 (Right) are ∆E eff BB (t)−∆E BB as a function of the discrete time (integer t/a with lattice spacing a = 0.1 fm) for t ≤ 2.5 fm, which would appear in typical numerical simulations. To obtain the data in this demonstration, we assign random fluctuations to R(t) whose magnitude increases exponentially in time and is comparable to that of our lattice data, and then calculate the central value and statistical error of ∆E eff BB (t) at each t. This figure clearly demonstrates that it is almost impossible to have data with enough accuracy at t 8 − 10 fm in current simulations.
Another point which is noteworthy in Fig. 1 (Right) is that the plateau-like behaviors show up at t 1 − 2 fm. Provided that t 1 − 2 fm is the region where statistical errors for two-baryon system may be controlled in present-day lattice simulations, one may easily misidentify this plateau-like behaviors as a real plateau. The estimate for ∆E BB then contains the systematic error of ±4 MeV (b 2 /b 1 = ±0.1), which is significant to the typical value of ∆E BB , 10 MeV or less, for the two-baryon system.
The behaviors demonstrated in Fig. 1 certainly depend on parameters such as δE el , δE inel , c 1 /b 1 , b 2 /b 1 , and a fake plateau may or may not appear in a specific lattice QCD simulation. There exists a potential danger, however, that a fake plateau appears during a search of a plateau at accessible t, by tuning, for example, the interpolating operators. Thus it is always necessary to find the explicit evidence that the obtained plateau-like structure is not fake. Due to the exponentially increasing noise in time, this task is extremely difficult, and becomes even impossible practically at physical quark masses with a larger lattice box, since δE el becomes much smaller as discussed before.

Fitting range for temporal correlations
In the following sections, we will analyze the lattice data to show explicitly the problem raised in Sec. 2.1. In the time correlation analysis of the actual lattice data for two baryons, one can only utilize the data at the moment up to t = t max ∼ 2 fm in the temporal direction due to the exponential decrease of S/N in E eff BB (t) for large t. Also, the lower limit of t = t min is constrained by the ground state saturation by a single hadron in m eff B (t). Therefore, a practical procedure adopted in many of the previous works are to look for the plateau of ∆E eff BB (t) = E eff BB (t) − 2m eff B (t) under the expectation that some cancellation of systematic as well as statistical errors. We will adopt the same practical procedure below for choosing the fitting window in the temporal direction, and show that the procedure leads to inconsistent results as expected.
For measurements of multi-baryon correlation functions, we employ two different sources, one is the smeared quark source, the other is the wall quark source, to check whether plateau-like behaviors agree between two sources. For the smeared source, we take exactly the same smearing function and parameters used in Ref. [2]: Quark propagators are solved using the exponentially smeared source of the form that 2 after the Coulomb gauge fixing is applied to gauge configurations. For the wall source, we take q w (t) = y q( y, t). (3.2) Relativistic interpolating operators for proton, neutron and Ξ are given by where C = γ 4 γ 2 is the charge conjugation matrix, α and a, b, c are the spinor index and color indices, respectively. Non-relativistic operator exclusively used in Ref. [2] is obtained by replacing Cγ 5 in Eq. (3.3) by Cγ 5 (1 + γ 4 )/2. We employ both non-relativistic and relativistic operators in this paper to estimate the systematic errors from the different choices.
For the source operators, we insert q s or q w in each flavor of Eq. (3.3) or its nonrelativistic variant. In the case of the smeared source, we take the same x for all quarks in Eq. (3.3), as is done in Ref. [2]. For sink operators, on the other hand, each baryon operator is composed of point quark fields, and projected to zero spatial momentum by averaging over the spatial position. For the choice of relativistic and non-relativistic baryon operator, we consider the same choice at both source and sink in this study. In the case of 4 He, however, the non-relativistic nucleon operator is used for the source regardless of the choice for the sink operator, in order to reduce the numerical cost. Altogether, we consider four different combinations for each multi-baryon system, two from wall and smeared quark sources times two from relativistic and non-relativistic baryon operators.
Quark propagators are solved with the periodic boundary condition in all directions using the quark source described above. Correlation functions (with relativistic and nonrelativistic baryon operators) are then calculated accordingly, where we use the unified contraction algorithm (UCA) [15]. UCA significantly reduces the computational cost of correlation functions, in particular for those of 3 He and 4 He. (See also related works [16][17][18][19].) On each gauge configuration, we repeat the measurement of correlation functions for a number of smeared sources at different spatial point and time slices and a number of wall sources at different time slices. For the 48 3  are calculated not only in one direction but also in other three as the time direction on each configuration using the rotational symmetry. In order to reduce the computational cost for the quark solver, the following stopping conditions |r crit | for the residual error are employed: |r crit | = 10 −4 (10 −12 ) for smeared (wall) source on the 32 3 × 48 lattice, |r crit | = 10 −6 (10 −4 ) for smeared (wall) source on the 40 3 × 48 lattice, |r crit | = 10 −6 (10 −4 ) for smeared (wall) source on the 48 3 × 48 lattice, |r crit | = 10 −6 for smeared source on the 64 3 × 64 lattice and |r crit | = 10 −6 (for half of the total statistics) or 10 −12 (for the other half) for wall source on the 64 3 × 64 lattice. In all cases, we check that systematic errors associated with the choice of the stopping condition is much smaller than the statistical fluctuations in this study. Nonetheless, we correct these errors by using the all-modeaveraging (AMA) technique [20,21] with the translational invariance. 3 Here, the AMA corrections for relaxed stopping conditions of |r crit | = 10 −4 or 10 −6 data are obtained by the corresponding computations with "exact" solver (|r crit | = 10 −12 ) with the following measurements: 1 source for smeared source on the 32 3 ×48 lattice, 1 (2) sources for smeared (wall) source on the 40 3 × 48 lattice, 4 × 1 (4 × 2) sources for smeared (wall) source on the 48 3 × 48 lattice and 1 × 1 (4 × 1) sources for smeared (wall) source on the 64 3 × 64 lattice, where the factor of 4 or 1 for the 48 3 × 48 and 64 3 × 64 lattices denotes the enhancement factor in statistics by the rotational symmetry. The lattice parameters, and the number of configurations as well as the number of smeared sources and wall sources are tabulated in Table 1. As noted above, the number of measurements for the 48 3 × 48 and 64 3 × 64 lattices can be increased by exploiting the rotational symmetry, and the factor of 4 in Table 1 represents this enhancement. In addition, for each measurement on any lattice volumes, we calculate correlation functions in forward and backward propagations (t > 0 and t < 0, respectively) and take an average to improve the signal. The corresponding factor of 2 is not included in Table 1. We note that the numbers of configurations and measurements for the smeared source in this work are much larger than those in [2].
on a lattice volume with La = 2.9, 3.6, 4.3 and 5.8 fm, respectively, the ratio of the number of measurements in this work to Ref. [2] amounts to be about 4.0, 2.8, 8.0 and 1.7 for each volume.
In our analyses, statistical errors are estimated by the jackknife method. We find that the auto-correlation in terms of configuration trajectory is small by observing that the statistical errors are almost independent among the choices of bin-size of 2, 5, 10, 20 configurations. Hereafter, we show the results obtained with the bin-size of 10 configurations (100 trajectories), unless otherwise stated.
Let us first consider the ΞΞ system in the spin-singlet channel with zero orbital angular momentum, ΞΞ( 1 S 0 ), where the interpolating operator is given by The reasons to choose this channel are twofold: Firstly, the signal to noise ratio for strange baryons is better than non-strange baryons. Secondly, in the flavor SU(3) limit, it belongs to the same 27 multiplet as the N N ( 1 S 0 ), so that one may obtain some insights into the bound dineutron suggested in previous works. To make a firm connection to previous works, we start our analyses with the smeared source Eq. (3.1) and later consider the case with the wall source.
One finds a plateau-like behavior in Fig. 2 (Middle left) for 10 ≤ t/a ≤ 18 before the explosion of the noise over the signal for larger t. As we have argued in Sec. 2 by the mock data, such a plateau is likely to be fake due to the contamination of the higher scattering states. Nevertheless, following the practical procedure taken by the previous works, let us try an exponential fit of R ΞΞ (t) in this "plateau" region and to take a large volume extrapolation. The fitted result is shown by the horizontal bars (the thick line and the thin lines are the central value and the 1σ statistical errors, respectively). We perform similar analyses for other volumes and also for relativistic interpolating operators, and results for ∆E ΞΞ ( 1 S 0 ) are summarized in Table 2. The numbers in the first parenthesis denote the statistical error, while the numbers in the second parenthesis denote systematic errors from the fit. Taking the same criterion adopted in Ref. [2], we estimate the systematic errors by variations among the fitting window as [t min ± 1, t max ± 1] 4 .
In . They indicate the existence of a ΞΞ bound state in the 1 S 0 channel, which is qualitatively "consistent" with previous studies finding dineutron bound state at this quark mass. Obviously the big question is whether such conclusion is reliable or not as we have discussed in Sec. 2.  To answer the above question solely in terms of the lattice data, let us move on to analyze ΞΞ( 3 S 1 ) with the same fitting procedure. In this case, the interpolating operator is given by Ξ 0 In the flavor SU(3) limit, this channel is in the 10 multiplet where no N N channels belong to.
One finds again that a "plateau"-like behavior in 11 ≤ t/a ≤ 18 before the explosion of the noise over the signal in larger t as shown in Fig. 2 (Middle right) in the case of nonrelativistic operator on the 48 3 × 48 lattice. We perform the same analysis in other volumes and also for relativistic interpolating operators. In Fig. 2 (Lower right), ∆E ΞΞ ( 3 S 1 ), with statistical and systematic errors added in quadrature, is plotted as a function of 1/L 3 , together with the values at infinite volume obtained by linear extrapolation. The results are summarized in Table 2 with the infinite volume limit, ∆E ΞΞ ( 3 S 1 ) = 6.81(1.04)( +0.52 −0.48 ) MeV (non-relativistic operator) and 12.20(94)( +0.02 −0.12 ) MeV (relativistic operator). These results in the 3 S 1 channel clearly indicate that the procedure to analyze the data was wrong as expected. If one could correctly identify the ground state energy of a two particle system on the finite lattice, its infinite volume extrapolation must be either zero (for the scattering state) or negative (for the bound state): Positive definite ∆E ΞΞ ( 3 S 1 ) as seen in Fig. 2 (Lower right) cannot be allowed. Therefore, we conclude that the plateaux seen in ∆E eff ΞΞ (t) for spin-singlet and spin-triplet channels are fake and are likely to be the mirages of true plateaux located in much larger t as we have discussed in Sec. 2.
To backup the conclusion obtained with the smeared source, let us now analyze the lattice data with the wall source. For the ∆E ΞΞ ( 3 S 1 ) channel, one finds again that a "plateau"-like behavior in 14 ≤ t/a ≤ 18 before the explosion of the noise over the signal in larger t as shown in Fig. 3 (Middle right) for the non-relativistic interpolating operators. We perform the same analysis in other volumes and also for the relativistic operators. In Fig. 3 (Lower right), ∆E ΞΞ ( 3 S 1 ) with the wall source is plotted as a function of 1/L 3 , together with the values at infinite volume obtained by linear extrapolation. The results are summarized in Table 2 Table 3, we summarize the results ∆E ΞΞ in all four cases which we have studied in this section. The positive ∆E ΞΞ ( 3 S 1 ) for the smeared source is not allowed physically, and there are apparent inconsistencies between the results of the smeared source and those  of the wall source. These are convincing enough that the previous works on temporal correlations have been looking at the fake plateaux just before the explosion of the noise over the signal as discussed in Sec.  Table 2. Summary of ∆E ΞΞ for both 1 S 0 (upper) and 3 S 1 (lower) channels from smeared and wall sources with range of an exponential fit, together with infinite volume extrapolations. On each volume, results from both relativistic and non-relativistic baryon operators are given. smeared source wall source ∆E ΞΞ ( 1 S 0 ) < 0 (bound state) 0 (no bound state) ∆E ΞΞ ( 3 S 1 ) > 0 (physically not allowed) 0 (no bound state) Table 3. Comparison of ∆E ΞΞ for different channels and different sources at infinite volume. Those are obtained by fitting "plateau"-like structure of ∆E eff ΞΞ (t) in the region of t just before explosion of the signal to noise ratio.
After clarifying the problem of fitting fake plateaux in ΞΞ systems, let us now turn our discussions to the N N systems in the N N ( 1 S 0 ) and N N ( 3 S 1 ) channels to show that the same problem takes place. Interpolating operators for these channels are given by N 1 N 2 − N 2 N 1 with N = p, n for N N ( 1 S 0 ), which belongs to 27 in the irreducible representation of the flavor SU(3), and by p α n α − n α p α with α = 1, 2 for N N ( 3 S 1 ), which belongs to 10 * representation. Note that the N N ( 1 S 0 ) is in the same flavor-SU(3) multiplet with ΞΞ( 1 S 0 ), while N N ( 3 S 1 ) belongs to the different flavor-SU(3) multiplet with ΞΞ( 3 S 1 ), so that we do not expect qualitative resemblance between N N ( 3 S 1 ) and ΞΞ( 3 S 1 ).
The upper two panels of Fig. 4 shows 2m eff N (t) and the effective energy E eff N N (t) for the smeared source with the non-relativistic nucleon operator on the 48 3 × 48 lattice in the N N ( 1 S 0 ) channel (Left) and in the N N ( 3 S 1 ) channel (Right).
The effective energy shifts from the smeared source on the 48 3 × 48 lattice are shown in the middle two panels in Fig. 4: Left (Right) panel for the 1 S 0 ( 3 S 1 ) channel with the non-relativistic nucleon operator. The explosion in the noise to signal ratio takes place for smaller t/a than that for ∆E eff ΞΞ (t) due to larger statistical errors in the N N case. We try to fit the plateau-like structure just before the explosion typically in the range 12 ≤ t/a ≤ 16. Obviously, we already knew from the discussions in the previous sections that such a plateau-like structure is fake. Our aim here (as in the case of the ΞΞ) is to show that the results of such fitting procedure adopted in previous literature do not make much sense.
In Table 4, results of ∆E N N on four volumes for the smeared source and for nonrelativistic and relativistic operators are summarized in the middle column. The fitting range for N N is relatively earlier than that for ΞΞ due to larger statistical errors. Systematic errors are estimated by changing the upper and lower limit of the fitting window by one unit of t/a as we have done in the case of ΞΞ.
The lower panels of Fig. 4 shows ∆E N N in the 1 S 0 channel (Left) and in the 3 S 1 channel (Right) as a function of 1/L 3 , together with the linear extrapolation in 1/L 3 to the infinite volume. In each figure, results of the non-relativistic and relativistic operators are plotted with the numerical data given in Table 4 . This agreement simply implies that the our present analysis and the previous analysis are consistent with each other and does not necessarily imply that there is indeed a bound state in these channels. This can be seen explicitly by the results of the wall source as shown below.
We now repeat the same analyses by changing the smeared source to the wall source. The results are summarized in Fig. 5 with the data in the right column in Table 4. The lower panels of Fig. 5 indicate that (i) the numbers are significantly different from those with the smeared source and (ii) there is no strong evidence of the bound states in both 1 S 0 and 3 S 1 channels. In fact, we obtain, for the wall source with non-relativistic operator,    Ref. [2] (non-rela.) −11.5(1.1)(0.6) − Table 4. A summary of ∆E N N for smeared and wall sources with both relativistic and nonrelativistic operators on four volumes and corresponding exponential fit ranges, together with infinite volume extrapolations. The result of the previous work with the same lattice setup is shown in the column [2] (non-rela.).

3 He and 4 He systems
We now consider 3 He (2 protons and 1 neutron) and 4 He (2 protons and 2 neutrons). Since m u = m d in 2+1 flavor QCD, 3 He is identical to triton, 3 H (1 proton and 2 neutrons), as far as its mass is concerned. Upper left (right) panel of Fig. 6 shows the effective energy shift ∆E eff 3 ) on the 48 3 × 48 lattice for both smeared and wall sources with the non-relativistic operator. The explosion in the noise to signal ratio from even smaller t/a than that of the two-nucleon case. We try to fit the plateau-like structure just before the explosion typically in the range 10 ≤ t/a ≤ 14. In Table 5, results of ∆E3 He and ∆E4 He on four volumes for smeared as well as wall sources and for non-relativistic as well as relativistic operators are summarized. Systematic errors are estimated by changing the upper and lower limit of the fitting window by one unit of t/a.
Middle left (right) panel of Fig. 6 shows ∆E3 He (∆E4 He ) from the smeared source as Ref. [2] (non-rela.) −43(12)(8) − Table 5. A summary of ∆E3 He and ∆E4 He for smeared and wall sources with both relativistic and non-relativistic operators on four volumes and corresponding exponential fit ranges, together with infinite volume extrapolations. a function of 1/L 3 , together with the linear extrapolation in 1/L 3 to the infinite volume, for both non-relativistic and relativistic operators. Lower left (right) panel of Fig. 6 shows ∆E3 He (∆E4 He ) from the wall source as a function of 1/L 3 , together with the linear extrapolation in 1/L 3 to the infinite volume, for both non-relativistic and relativistic operators. As in the case of N N , the result of the smeared source and that of the wall source do not agree: The former indicates the bound states for both 3 He and 4 He as suggested in [2], while the latter shows no strong evidence of such bound states.

Conclusions
In this paper, we have addressed the issue of the single state saturation of the temporal correlation function for the multi-baryons by employing (2+1)-flavor lattice QCD at m π = 0.51 GeV on four lattice volumes with L = 2.9, 3.6, 4.3 and 5.8 fm. A major difference between the single baryon and multi-baryons on the lattice is that there appears energy levels corresponding to the elastic scattering for the latter. The level spacings become smaller as L becomes larger, since they correspond to the continuum states for L → ∞. Therefore, it is required to take large temporal distance t between the source and sink operators to isolate the ground state of multi-baryons. This is, however, very difficult due to the exponential growth of the noise over the signal which has been known to be a characteristic feature of the multi-baryon correlations. In such a situation, one may be misled by a fake plateau of the effective energy shift ∆E eff (t) at intermediate values of t before the explosion of the noise takes place. We have demonstrated, by using the mock data, that the above situation can easily happen with a slight contamination of the elastic scattering state. Then we analyzed the lattice data in ΞΞ( 1 S 0 ) and ΞΞ( 3 S 1 ) channels to show explicitly that the same situation indeed takes place for the real data. By adopting the smeared source operator used in [2] and the wall source operator, we fit the plateau-like structure around t/a ∼ 15 and find that the results of ∆E ΞΞ at each L as well as those extrapolated to L → ∞ turn out to disagree with each other between two sources. This implies that the ground state saturation is not achieved in such intermediate values of t. Moreover, we found that ∆E ΞΞ ( 3 S 1 ) > 0 for the smeared source at L → ∞, which is not physically acceptable.
One may suspect that the above disagreement originates from slower temporal convergence of single baryon for the wall source than the smeared source. However, this is not necessarily the case, since the plateau of the effective energy shift shows much stronger dependence on the change of the two-baryon sink operators for the smeared source as shown in Appendix A. In fact, one can explicitly show, by using the HAL QCD method, that the smeared source has significantly larger contamination from the two-baryon excited states. The details will be reported in a forthcoming paper [22].
We have applied the same analysis also to N N ( 1 S 0 ), N N ( 3 S 1 ), 3 He and 4 He, although the statistical errors become lager for non-strange baryons than those for ΞΞ. Again, the results of the two sources do not agree with each other: The smeared source indicates that there are bound states in all these channels, while no definite signatures of the bound states are found for the wall source.
By combining the general theoretical considerations and the numerical evidences, we conclude that the plateau-like structure seen at the moderate values of t in the temporal correlation for multi-baryons should be considered as a "mirage" in the sense that the true signals are located in much larger t with different values of ∆E eff (t). This also casts strong doubt on the recent works on the basis of the plateau fitting of the temporal correlations [2,3,16,[23][24][25][26][27][28][29][30][31][32][33][34][35][36], almost all of which claim the existence of bound multi-baryons (such as dineutron, deuteron, 3 He, 4 He, and other strange multi-baryons). At least, one should use more sophisticated approaches than the plateau fitting, such as the Bayesian fitting, Black box, or variational methods to extract the ground state energy from the temporal correlators (see e.g. Ref. [37] for the review of these methods and the applications to singlehadron spectroscopy.) A trustable way to examine the validity of these results is to study the L-dependence of ∆Eà la Lüscher's finite volume formula. Detailed analysis along this line will be reported in another forthcoming paper [38].
It should also be noted that the use of the full space-time correlations (HAL QCD method) instead of only the temporal correlations has been shown to solve the single-state saturation problem discussed in this paper [9]. Detailed examination between the results from the temporal correlation alone and those from the space-time correlation by using the same lattice data as those in the present paper will be also reported in a forthcoming paper [22].

A Sink operator dependence
In the main text, we investigated the reliability of the plateau-like behavior using different source operators. In this Appendix, we make similar analysis by using different sink operators. We consider the ΞΞ system in the 1 S 0 channel as a representative case and start with the following temporal correlation function, The interpolating operator for Ξ( x, t) is given by Eq. (3.3) and we consider only the relativistic operator in this Appendix. The source operator, J ΞΞ , is taken to be the same as those used in Sec. 3, with either of the smeared source or of the wall source. The sink operator is a combination of the two local Ξ operators with a smearing function g(r) [45,46]. The temporal correlation C ΞΞ (t) in Sec. 3 corresponds to the case g(r) = 1. The effective energy, E eff ΞΞ (t), and the effective energy shift, ∆E eff In the following analysis, we adopt g(r) with the following form, In Fig. 7 (Left), we plot ∆E eff ΞΞ (t) for four different sink operators in the case of the smeared source. Although we find a plateau-like behavior for each sink operator, the values of ∆E eff ΞΞ (t) do not agree among different sink operators. Such a large sink-operator dependence indicates that the contamination from the elastic scattering states in the finite volume causes fake plateaux as demonstrated in Sec. 2. The true plateau may be identified at much larger values of t, but the explosion of the noise prohibits to extract sensible signal at large t as we discussed in the text. Shown in Fig. 7   sink operators in the case of the wall source. For each sink operator, we find a plateau-like behavior: In this case, it happens that the values of ∆E eff ΞΞ (t) agree among different sink operators within statistical errors.