Lattice QCD and baryon-baryon interactions: HAL QCD method

In this article, we review the HAL QCD method to investigate baryon-baryon interactions such as the nuclear forces in lattice QCD. We first explain our strategy to investigate baryon-baryon interactions by defining potentials in field theories such as QCD in detail. We introduce the Nambu-Bethe-Salpeter (NBS) wave functions in QCD for two baryons below the inelastic threshold. We then define the potential from NBS wave functions in terms of the derivative expansion, which is shown to reproduce the scattering phase shifts correctly below the inelastic threshold. Using this definition, we propose a method to extract the potential in lattice QCD. Secondly, we discuss pros and cons of the HAL QCD method, by comparing it with the conventional method, where one directly extracts the scattering phase shifts from the finite volume energies through the L\"uscher's formula. We give several theoretical and numerical evidences that the conventional method so far fails to work on baryon-baryon interactions due to contaminations of elastic excited states. On the other hand, we show that such a serious problem can be avoided in the HAL QCD method by defining the potential in an energy-independent way. We also discuss systematics of the HAL QCD method, in particular errors associated with a truncation of the derivative expansion. Thirdly, we present several results obtained from the HAL QCD method, which include the (central) nuclear force, the tensor force, the spin-orbital force, and three nucleon force. We finally show the latest results calculated at the nearly physical pion mass, $m_\pi \simeq 146$ MeV, including hyperon forces which lead to form $\Omega\Omega$ and $N\Omega$ dibaryons.


INTRODUCTION
How do nuclear many-body systems emerge from the fundamental degrees of freedom, quarks and gluons? It has been a long-standing problem to establish a connection between nuclear physics and the fundamental theory of strong interaction, quantum chromodynamics (QCD). In particular, nuclear forces serve as one of the most basic constituents in nuclear physics, which are yet to be understood from QCD. While so-called realistic nuclear forces [1,2,3] have been established with a good precision, they are constructed phenomenologically based on scattering data experimentally obtained. Recent development in effective field theory (EFT) provides a more systematic approach for nuclear forces from a viewpoint of chiral symmetry in QCD [4,5,6,7,8], whose unknown low-energy constants, however, cannot be determined within its framework but are obtained only by the fit to the experimental data. Under these circumstances, it is most desirable to determine nuclear forces as well as general baryon-baryon interactions from firstprinciples calculations of QCD, the lattice QCD method. Once baryon forces are extracted from QCD, we can solve finite nuclei, hypernuclei and nuclear/hyperonic matter by employing various many-body techniques developed in nuclear physics. The outcome is expected to make a significant impact on our understanding of nuclear astrophysical phenomena such as supernovae, binary neutron star merges and nucleosynthesis.
In this paper, we review the HAL QCD method to determine baryon-baryon interactions in lattice QCD. In this method, integral kernels, or so-called "potentials", are first extracted from lattice QCD, and physical observables such as scattering phase shifts and binding energies are calculated by solving the Schrödinger equation with obtained potentials in the infinite volume. We show that the notion of potential can be rigorously introduced as a representation of the S-matrix in quantum field theories as QCD. The essential point is that the potentials are defined through the Nambu-Bethe-Salpeter (NBS) wave functions, in which the information of phase shifts are encoded in their asymptotic behaviors. We employ a non-local and energy-independent potential where the non-locality is defined through the derivative expansion. In particular, energy-independence of the potential is useful since one can extract the potential from the ground state as well as elastic excited states simultaneously. This enables us to avoid the notorious signal-to-noise issue for multi-baryon systems in lattice QCD (or the ground state saturation problem), and to make a reliable determination of baryon-baryon interactions.
In lattice QCD, there also exists a conventional method, in which phase shifts are obtained from the finite volume energies through the Lüscher's formula. While theoretical bases are well established for both conventional method and HAL QCD method, the numerical results in some cases at heavy pion masses have shown inconsistencies with each other. In this paper, we make a detailed comparison between two methods, scrutinizing possible sources of the systematic errors. In particular, we examine whether the systematic errors associated with the excited state contaminations are controlled or not in the procedure of the conventional method in the literature ("the direct method"), namely, plateau fitting for the ground state at early Euclidean times. We also examine the systematic errors in the HAL QCD method, the truncation error of the derivative expansion. We show theoretical and numerical evidences that the inconsistency between two methods is originated from excited state contaminations in the direct method. We also present how the inconsistency can be actually resolved if and only if finite energy spectra are properly obtained in the conventional method, with an improved method rather than the naive plateau fitting.
After establishing the reliability of the HAL QCD method, we present the numerical results of nuclear forces from the HAL QCD method at various lattice QCD setups. The results at heavy pion masses for central and tensor forces are shown and their quark mass dependence as well as physical implications are discussed. The calculations of spin-orbit forces and three-nucleon forces are also given. Once nuclear forces are obtained, one can solve nuclear many-body systems with the obtained potentials. We study finite nuclei, nuclear equation of state and structure of neutron stars based on lattice nuclear forces at heavy pion masses. Finally, the latest results of nuclear forces near the physical pion mass are presented, as well as hyperon forces, which are shown to generate ΩΩ and N Ω dibaryons. This paper is organized as follows. In Sec. 2, we discuss methods to study baryon-baryon interactions from lattice QCD. After briefly introducing the conventional method and its actual practice, called the "direct method", we describe the detailed theoretical formulation as well as its practical demonstration for the newly developed method, the HAL QCD method. In Sec. 3, we discuss pros and cons of these two methods, and compare the numerical results at heavy pion masses. We present evidences that the results from the direct method suffer from uncontrolled systematic errors associated with the excited state contaminations. In Sec. 4, we summarize results on nuclear potentials in the HAL QCD method. After reviewing the results obtained at heavy pion masses for central and tensor forces in the parity-even channel as well as spin-orbit forces and three-nucleon forces, we present nuclear many-body calculations based on lattice nuclear forces for double-magic nuclei, equation of state and the structure of neutron stars. Latest results for nuclear forces near the physical pion mass are also given. In Sec. 5, we present hyperon forces near the physical pion mass, which lead to ΩΩ and N Ω dibaryons. Sec. 6 is devoted to the summary and concluding remarks.

TWO BARYON SYSTEMS IN LATTICE QCD
In lattice QCD, the 2-pt function for a hadron H, created by O † H and annihilated by O H , is expressed as where |n, E n ( p) is the n-th one-particle state with a mass m n , a momentum p and an energy E n ( p) = m 2 n + p 2 , and ellipses represent contributions from multi particle states. We here assume m 0 < m n>0 , so that m 0 is the hadron mass for the ground state, which can be extracted from the asymptotic behavior of the 2-pt function in the large t as So far, this method in lattice QCD (and QED) has successfully reproduced light hadron spectra [9] including the proton-neutron mass splitting [10]. A simple extension of the method, however, does not work well for an investigation of hadron interactions. For example, the 2-pt function of two baryons in the center of mass system behaves in the large t as where we obtain the lowest energy E BB , which behaves as E BB = 2m B in the absence of bound states or E BB = 2m B − ∆E in the presence of bound states. Here m B is the corresponding baryon mass and ∆E > 0 is the binding energy of the lowest bound state. Only the binding energy of the bound state can be extracted by this simple method and thus more sophisticated methods are required. Currently there are two methods to investigate hadron interactions in lattice QCD, the direct method (or finite volume method) and the HAL QCD method, which are explained in this section.

Direct method
The most widely used method to investigate hadron interactions in lattice QCD is to extract scattering phase shifts from energy eigenvalues in 3-dimensional finite boxes through the Lüscher's finite volume  Figure 1. A determination of k cot δ 0 (k)/m π from energies of the two nucleon state in the finite volume. Taken from [12].
formula [11]. For example, in the case of the S-wave scattering phase shift, δ 0 (k), the formula reads where k is determined through E BB (L) = 2 k 2 + m 2 B with E BB (L) being the energy of the two baryon measured in lattice QCD on a finite box with the spatial extension L as in eq. (3). Only the discrete sets of point (k 2 , k cot δ 0 (k)), which satisfies eq. (4), are realized on a given volume L 3 . Thus the scattering phase shift δ 0 (k) at the corresponding k can be extracted in lattice QCD, simply by measuring the finite volume energy, E BB (L).
In Fig. 1, we illustrate how scattering phase shifts and the bound state energy can be extracted by this method in the case of the N N scatterings. In the figure, the red solid line represents the effective range expansion (ERE) for k cot δ 0 (k)/m π at the Next-to-Leading order (NLO) as where the scattering length a 0 and the effective range r 0 are taken to be a 0 m π = 16.8, r 0 m π = 1.9 for A distribution of the allowed k 2 for k 2 > 0 becomes denser as the volume increases, so as to be continuous in the infinite volume limit, while a sequence of discrete points for k 2 < 0 leads to an accumulation point, which corresponds to the scattering state at k 2 = 0 in the left figure or the bound state pole, denoted by the black solid circle in the right figure. It is noted here that the bound state pole appears as the intersection between the ERE and the bound state condition, − −(k/m π ) 2 (black solid line). To see this, we first where S(k) is the S-matrix for the N N elastic scattering. The bound state energy κ b can be extracted from the pole of this S-matrix as where β 2 b is real and positive for physical poles [13]. Thus at k 2 −κ 2 b , we have which means that the binding momentum k = iκ b is given by an intersection between k cot δ 0 (k) and − √ −k 2 . Moreover, since the slope of k cot δ 0 (k) must be smaller than that of − √ −k 2 as a function of k 2 at the bound state pole, as in the case of Fig. 1 (Right). The finite volume analysis thus provides not only an infinite volume extrapolation of the binding energy but also a novel way to examine the normality of the result in the direct method [12].

Formulation
The HAL QCD method, another method to investigate hadron interactions in lattice QCD, employs the equal time Nambu-Bethe-Salpeter (NBS) wave function, defined by where |N N, W k is the N N eigenstate in QCD with the center of mass energy W k = 2 k 2 + m 2 N and the nucleon mass m N , and N (x, t) is a nucleon (annihilation) operator, made of quarks. Other quantum numbers such as spin/isospin of two nucleons are suppressed for simplicity. We mainly use where C = γ 2 γ 4 is the charge conjugation matrix, q = u(d) for proton (neutron). Other choices such as smeared quarks are possible here, and such arbitrariness is considered to be a choice of the scheme for the definition of the NBS wave function or the potential. (See [14] for such an example.) Throughout this paper, we consider the N N elastic scattering, so that W k < W th ≡ 2m N + m π , where m π is the pion mass. Note that this condition is also necessary for the finite volume method in the previous subsection.
Since interactions among hadrons are all short-ranged in QCD, there exists some length scale R, beyond which (i.e. r ≡ |r| > R ) the NBS wave function satisfies the Helmholtz equation as Furthermore, it behaves for large r > R as where Y lm is the spherical harmonic function for the solid angle Ω r of r, and we ignore spins of nucleon for simplicity 1 . Here it is important to note that the NBS wave function contains information of the phase δ l (k) of the S-matrix for the orbital angular momentum l, which is a consequence of the unitarity of the S-matrix in QCD [17,18].
In the HAL QCD method, the non-local but energy-independent potential is defined from the NBS wave function through the following equation, for W k < W th , and eq. (12) implies U (r, r ) = 0 for r > R. While an existence of U (r, r ) has been shown in [19,16,20], the non-local potential which satisfies eq. (14) is not unique. Thus we have to define the potential uniquely, by specifying how to extract it. For this purpose, we introduce the derivative expansion, , whose lowest few orders for the N N with a given isospin channel are written as where V 0 (r) is the central potential, V σ (r) is the spin dependent potential with σ i being the Pauli matrix acting on the spinor index of the i-th nucleon, V T (r) is the tensor potential with the tensor operator S 12 = 3(r · σ 1 )(r · σ 2 ) − (σ 1 · σ 2 ) (r ≡ r/r), and V LS (r) is the spin-orbit (LS) potential with the angular momentum L = r × p and the total spin S = (σ 1 + σ 2 )/2. It is noted that an expansion of the non-local potential is not unique. For example, we may improve the convergence of the expansion by modifying the ∇ operator [21].
Once we obtain the approximated potential at lowest few orders, we can calculate the scattering phase shifts or the binding energies of possible bound states by solving the Schrödinger equation with this potential in the infinite volume. We can also check how good the approximated potential is, by increasing the order of the expansion. Needless to say, the approximated potential depends on momenta of input wave functions. As pointed out in [22], these dependences of the approximated potentials have been misidentified with those of the non-local potential in the literature [23]. In the next subsection, we will explicitly demonstrate how this procedure works.

Demonstration
In order to see how the scattering phase shifts can be obtained by the HAL QCD method, we consider the quantum mechanics for a spinless system with a separable potential, defined by The S-wave solution of the Schrödinger equation with this potential is given exactly by where which is the 4-th order polynomials in k 2 . In order to make the scattering phase shift a more complicated function of k 2 , we artificially modify the wave function from φ 0 k (r) to φ k (r) which is defined by where R is an infrared cutoff, and it is understood that the potential is modified accordingly. The continuity of φ k (r) and φ k (r) at r = R gives as well as C(k) = X/ sin(kR + δ R (k)). Hereafter, we study how the scattering phase shifts are obtained in the HAL QCD method.
The derivative expansion for the S-wave scatterings leads to and we consider to determine the potential in each order from φ k (r).
The leading order (LO) potential is given by while the next-to-leading order (NLO) potential is extracted as where Note that the potential in each order in the derivative expansion {V 0 (r), V 1 (r), · · · } are defined to be k-independent, while the potentials approximately obtained in each LO/NLO analysis, {V LO 0 (r; k)} and {V NLO 0 (r; k 1 , k 2 ), V NLO 1 (r; k 1 , k 2 )}, have implicit k-dependence due to the truncation error in the derivative expansion [22].
We calculate S-wave scattering phase shifts corresponding to these approximated potentials, and compare them with the exact phase shifts, δ R (k). Considering µ as a typical inelastic threshold energy in this model, we take k = 0 and/or k = µ for the following analysis. Fig. 2 shows the S-wave scattering phase shift δ(k) (Left) and k cot δ(k) (Right) as a function of k 2 , where all (dimensionful) quantities are measured in units of µ. In this example, we take ω = −0.017µ 4 , m = 3.30µ and R = 2.5/µ. In the figures, the exact phase shift δ R (k)(Left) or k cot δ R (k) (Right) is given by the blue solid line, while the LO approximations at k = 0 or k = µ are represented by orange and green solid lines, respectively. As seen from the figures, the LO approximation at k = 0 (orange), exact at k 2 = 0 by construction, gives a reasonable approximation at low energies (k 2 0) but deviates from the exact one at high energies near k 2 µ 2 . On the other hand, the LO approximation at k = µ (green) becomes accurate at higher energies near k 2 µ 2 but inaccurate at low energies near k 2 0. Combining two NBS wave functions, φ k 1 =0 (r) and φ k 2 =µ (r), one can determine the approximated potential at the NLO, V NLO (r, ∇), whose scattering phase shifts are represented by the red solid lines in the figures. The phase shifts at the NLO (red lines) gives reasonable approximations of the exact results (blue solid lines) in the whole range (0 ≤ k 2 ≤ µ 2 ), as they are exact at k 2 = 0 and k 2 = µ 2 by construction. If we increase the order of the expansion more and more, the approximation becomes better and better. 2 Using this model, let us compare the direct method and the HAL QCD method. At the LO, the direct method gives either k cot δ(k) at k 2 = 0 or k 2 = µ 2 without any information about the effective range, which only gives the LO ERE (an orange dashed line or a green dashed line in the right figure.) Thus the LO potentials approximate the exact k cot δ(k) much better (the orange solid line or the green solid line). In the direct method, the ERE at NLO is obtained by combining the data at k 2 = 0 and k 2 = µ 2 as which is given by a red dashed line in the right figure. By comparing the HAL QCD method with potentials at NLO (the red solid line) and the direct method with NLO ERE (the red dashed line), the former leads to a better approximation of the exact result than the latter, since higher order effects in ERE in terms of k 2 are included in the former. Note, however, that sufficiently precise data in the direct method can also evaluate higher order ERE terms than NLO, in principle. Figure 2. The scattering phase shifts δ(k) and k cot δ(k) as a function of k 2 . See the main text for more details.

Dependence of the LO N N potential on energy and partial waves
In this subsection, we consider effects of higher order terms in the derivative expansion for the N N in QCD. 530 MeV [26]. As seen from the figure, two NBS wave functions look very different from each other. In particular, the right one vanishes on the boundary due to the APBC constraint.    (cyan solid circles) at E k 45 MeV. Although statistical fluctuations are larger for the latter, they look similar, suggesting that L 2 dependence of the potential is also small in this setup. If more accurate data show a difference of potentials between N N ( 1 S 0 ) and N N ( 1 D 2 ), one may determine the L 2 dependent term of the potential in the spin-singlet channel.

Time-dependent HAL QCD method
In order to extract the NBS wave functions on the finite volume in lattice QCD, we consider the 4-pt function given by (26) whereJ N N (t 0 ) is an operator which creates two nucleon states at time t 0 , A J n ≡ N N, W kn |J N N (0)|0 , and ellipses represent inelastic contributions, which become negligible at W th (t − t 0 ) 1. Like the direct method, one can extract the NBS wave function for the ground state from the above 4-pt function as for is the lowest (second-lowest) energy on the finite volume. The LO potential from the NBS wave function for the ground state is then extracted from F J (r, t) at large t. As will be discussed in the next section, however, it is numerically very difficult to determine F J (r, t) for two nucleons at such large t due to the bad signal-to-noise (S/N) ratio.
Fortunately, an alternative extraction is available for the HAL QCD method [27]. Let us consider the ratio of 4-pt function to the 2-pt function squared as which behaves for W th t 1, where inelastic contributions can be neglected. Noticing that we obtain We can approximately extract V (r, ∇) from R J (r, t) for (different) J's, as long as t satisfies the condition that W th t 1 (elastic state saturation), which is much easier than to achieve (W k 1 − W k 0 )t 1 (ground state saturation). We call this alternative extraction the time-dependent HAL QCD method.

A COMPARISON OF THE TWO METHODS AT HEAVIER PION MASSES
It is interesting to ask whether the attractions of the nuclear forces at low energies would become weaker or stronger if the pion mass were larger than the value in Nature. In principle, such a question can be answered by employing either the direct method or the HAL QCD method in lattice QCD. There exists, however, a qualitative discrepancy between the two methods on the answer to this question. As summarized in Table 1, the direct method tends to indicate that attractions between two nucleons become stronger as the pion mass increases, so that both deuteron and di-neutron form bound states, while the HAL QCD method suggests that the attractions become weaker and the bound deuteron does not exist at heavier pion masses. Note that the results from the direct method in the flavor SU(3) limit (N f = 3 in the table), NPL2013/NPL2017, CalLat2017 and Mainz2018, exhibit discrepancies with each other [12]. In addition, while both methods lead to the bound H-dibaryon at heavier pion masses, in particular, in the flavor SU(3) limit, the predicted binding energies differ even within the direct method: NPL2013 [33] gives 75(5)MeV at m π = 810 MeV, which is much larger than 19(10) MeV at m π = 960 MeV by Mainz2018 [36]. On the other hand, HAL2012 [37] gives 38(5) MeV at m π = 837 MeV from the HAL QCD method. These deviations seem to be too large to be explained by lattice artifacts.
In order to understand origins of these discrepancies, we have performed extensive investigations, whose results have been published in a series of papers [39,12,40,41], which will be explained in the following subsections.
and H-dibaryon in lattice QCD. NPL2013, NPL2017 and CalLat2017 employed the same set of gauge configurations. CalLat2017 found two states in each channel. In Mainz2018, dynamical 2-flavor with quenched strange quark configurations are employed and N f in the table (with † symbol) denotes the information in the valence quark sector. All values of ∆E correspond to those in the infinite volume limit except ones with * , which are values on the finite volumes. The number 0 in ∆E indicates the system is unbound in this channel.

Operator dependence in the direct method
In the direct method, reliable extractions of the two nucleon ground state energies are crucially important. As long as (W k 1 − W k 0 )t 1, the two nucleon correlation function is dominated by the ground state as so that the extracted ground state energy W k 0 depends neither the source operatorJ N N nor the sink operator J N N , while magnitudes of contaminations from excited states are affected by the choices of these operators. Since W k 1 − W k 0 (2π/L) 2 /m N on the finite box with the spacial extension L, t 4 fm is required, for example, for L 4 fm and m N 2 GeV at heavier pion masses. Due to the bad S/N ratio at such large t, however, authors in previous literature extracted the ground state energies at much smaller t, t ∼ 1 fm, by tuning the source operatorsJ N N in order to achieve a plateau of the effective energy shift ∆E eff N N (t) at such a small t, where [MeV]  Unfortunately, such a naive plateau fitting at earlier t may not be reliable due to contaminations from nearby excited states, which may easily produce (incorrect) plateau-like behaviors in effective energies. It was indeed demonstrated that plateau-like behaviors in effective energy shifts at small t can depend not only on the source operator but also on the sink operator: Plateaux disagree between the wall source (red circle) and the smeared source (blue square) in the left of Fig. 5, while plateaux depend on sink operators for the same smeared source in the right figure.
In order to see how easily contaminations from elastic-excited states can produce plateau-like behaviors at earlier t, let us consider the effective energy shift from the mockup data for R N N (t), given by where we take δE el. = 50 MeV for the typical lowest elastic excitation energy on L 4 fm at m N 1.5 GeV, and δE inel. m π 500 MeV for the lowest inelastic energy. Naively, it is expected that the correct plateau at ∆E N N for the ground state appears at t 1/δE el. 4 fm, which however is too large to have good signals for two baryons such as N N . By tuning the source operator, one may reduce coefficients b 1 and c 0 . Since the N N operator does not strongly couple to N N π state, we expect small c 0 and take c 0 = 0.01. On the other hand, N N operators easily couple to both ground and 1st elastic excited states as they become almost identical to each other in the infinite volume limit. We therefore take b 1 = 0.01 (the highly tuned operator) , b 1 = ±0.1 (the tuned ones) as well as b 1 = 0.5 (the untuned one). Fig. 6 (Left) shows ∆E eff N N (t) for these 4 examples with c 0 = 0.01, where random fluctuations and errors whose magnitude increase exponentially in t are assigned to R mockup N N (t). All examples show plateau-like behaviors at t 1 fm, but these four plateaux disagree with each other. As |b 1 | increases, the deviation between the values of these "pseudo plateaux" and the true value becomes larger. Contaminations of the elastic excited states can easily produce the plateau-like behavior at earlier t, and the t dependence of data alone cannot tell us which plateau is correct, or in other words, cannot tell which tuning is good.
Contaminations from inelastic states seem unimportant to produce the plateau-like behavior, as shown in Fig. 6 (Right), where the effective energy shift for c 0 = 0.01, 0.05, 0.1 with b 1 = −0.1 is plotted. All cases converge to almost the same pseudo plateau, while a pseudo plateau starts at later t for larger c 0 . It is noted that the multi-exponential fit does not work in this case at t 1.0 fm, which is much smaller than the necessary t 1/δE el. . The multi-exponential fit at such small t only separates the pseudo plateau from the inelastic contributions but is difficult to distinguish the ground state and the 1st excited state for the elastic states.  [29] for N N ( 1 S 0 ) as a function of (k/m π ) 2 . The solid red line and light red band represent the ERE fit and the corresponding error (statistical and systematic added in quadrature), respectively. The dashed lines are the finite volume formula for the corresponding volume. (Right) k cot δ 0 (k)/m π in NPL2015 [30] for N N ( 1 S 0 ) as a function of (k/m π ) 2 . Two ERE fits are performed depending on the lattice data to be used for the fit. The red line with the band represents the fit made by the authors in [12], while the blue line with the band is plotted by the authors in [12] using the fit result of NPL2015. Both figures are taken from [12].

Normality check in the direct method
While the check through operator dependence is useful, it requires extra calculations. We find that the finite volume formula in eq. (4) provides a simpler test, which tells us whether the ground state energies extracted by the plateau fitting give a reasonable ERE or not without extra calculations. We call this test a normality check [12]. Fig. 7 (Left) shows k cot δ 0 (k)/m π in YIKU2012 [29] as a function of k 2 /m 2 π for N N ( 1 S 0 ), where the solid red line represents the NLO ERE fit in eq. (5), and the light red bands shows statistical and systematic errors added in quadrature [12]. Contrary to a naive expectation from non-singular ERE behaviors, data align almost vertically, since ∆E N N is almost independent of the volume. In other words, according to the finite volume formula, the claimed "binding energy" (open circle) is too shallow to have such volume independent ∆E. Not only the central value of the NLO ERE fit gives singular parameters as ((a 0 m π ) −1 , r 0 m π ) = (5.27, 303.6) but also it violates the physical pole condition, eq. (9), at the crossing point (open circle). The singular and unphysical behaviors, in addition to the operator dependence of these data, strongly indicate that the naive plateau fitting employed in the direct method is unreliable. Another example is shown in Fig. 7 [30]. In this case, two different NLO ERE fits (red line/band and blue line/band) are performed depending on the lattice data to be used for the fit. It turns out that two ERE are inconsistent with each other, indicating that their lattice data themselves are "self-inconsistent". In addition, one of ERE (blue line/band) is found to violate the physical pole condition, eq. (9), at the crossing point (open circle). Similar symptoms are observed for all other data in the direct method claiming the existence of N N bound states at heavy quark masses [12]. 3

The source dependence and the derivative expansion in the HAL QCD method
The source operator dependence of the HAL QCD potential has been investigated in [40]. We observe a small difference at short distances, from which one can determine the N 2 LO potential, , which is nonzero only at r < 1.0 fm, where two LO potentials differ. We then extract the scattering phase shifts, using this N 2 LO potential.
The N 2 LO corrections turn out to be negligible at low energies, as shown in Fig. 9 (Left), where k cot δ 0 (k) is almost identical between V N 2 LO (r, ∇) (red solid circles) and V N 2 LO 0 (r) (blue solid squares). Furthermore, even the LO analysis for the wall source, V LO(wall) 0 (r) (black open diamond), is sufficiently good at low energies. As energy increases, the N 2 LO corrections become visible as seen in Fig. 9 (Right), where (k/m π ) 2 = 0.5 corresponds to ∆E 90 MeV for the energy shift from the threshold. It is noted that V N 2 LO 0 (r) (blue solid squares) gives a little closer results to N 2 LO results (red solid circles) than V LO(wall) 0 (r) (black open diamond) does.

Understanding pseudo plateaux
In this subsection, we explain why the wall source and the smeared source give inconsistent plateau behaviors, in the case of ΞΞ correlation functions as an example.
To this end, we consider the Hamiltonian (r), the LO potential from the wall source, since it works rather well at low energies as shown in the previous subsection. We first decompose R J ΞΞ (r, t) for J = wall/smear in terms of finite volume eigenfunctions of H as . Both are taken from [40].  Figure 9. (Left) k cot δ 0 (k)/m π as a function of (k/m π ) 2 at low energies, where δ 0 (k) is the scattering The corresponding δ 0 (k). Both are taken from [40].
where Ψ n (r) and ∆E n are normalized-eigenfunction and eigenenergy in the finite volume, respectively, and a J n (t) is the the overlapping coefficient extracted at t.
Then the correlation function for the source J in the direct method is given by ( 1 wall src. smeared src. Taken from [41].
Finally, approximating a sum over n by the lowest few orders, we reconstruct the behavior of the effective energy shift as a function of t as where we fix the overlapping coefficient b J n (t 0 ) at t = t 0 , and n max is a number of excited states used in the approximation.
In Fig. 10, we show reconstructed effective energy shift ∆E We find that the plateau-like structures in the direct method around t/a = 15 are well reproduced by ∆E J eff (t, t 0 = 13a) for both sources in Fig. 10 (Left). This indicates that the plateau-like structures in the direct method at this time interval are explained by the contributions from several low-lying states.
These plateau-like structures of course do not necessarily correspond to the true energy shift of the ground state. The fate of these structures is shown in Fig. 10 (Right), where we plot ∆E J eff (t, t 0 = 13a) at asymptotically large t. While the plateau-like structure for the wall source is almost unchanged, ∆E J eff (t, t 0 = 13a) for the smeared source gradually increases and reaches to the true value at t/a ∼ 100.
The above results clearly reveal that the plateau-like structures at t/a ∼ 15 for the smeared source are pseudo-plateaux caused by the contaminations of the excited states. Large contaminations from excited states in the case of the smeared source are not caused by the smearing, but are indeed implied by putting  two baryon operators on the same space-time point as where the above source operator couples to all momentum modes with almost equal weight. Since almost all previous studies on N N interactions in the direct method employed this type of the source operator, their conclusions on the existences of both deuteron and di-neutron are not valid due to large contaminations. 4

Consistency between the two methods
Once eigenmodes of H in the finite box are obtained, we can construct an improved sink operator for a particular eigenstate, whose correlation function with the J source is given by After the sink projection, the effective energy shifts agree well between wall and smeared sources around t/a ∼ 13, not only for the ground state but also for the 1st excited state. while the effective energy shifts for the ground state in the direct method without projection disagree between two sources. In particular, an agreement between two sources with sink projection for the 1st excited state is rather remarkable, since variational methods, usually mandatory for excited states in lattice QCD, are not used here. Furthermore, the plateaux of the effective energy shifts after the sink projection also agree with ∆E 0,1 of H (red bands). Note that the effective energy shift for the 1st excited state, ∆E wall,1 eff (t), has larger errors since the contribution of the 1st excited state in R wall ΞΞ (t) is much smaller. Although the sink operator projection utilizes the information of the HAL QCD potential to construct eigenfunctions, agreements in the effective energy shifts for the ground state as well as the 1st excited state provide a non-trivial consistency check between the HAL QCD method and the Lüscher's finite volume formula (with proper projections to extract the finite volume spectra). We thus conclude from Fig. 11 not only that the HAL QCD potential correctly describes the energy shifts of two baryons in the finite box for both ground and excited states but also that these energy shifts can be extracted through the Lüscher's finite volume formula if and only if the sink/source operators are highly improved. We emphasize that improvement of operators has to be performed not by the tuning of the plateau-like structures but by a sophisticated method such as the variational method [42] (or a method presented here). See [36] for a recent study toward such a direction.

NUCLEAR POTENTIAL
In this section, we summarize results on nuclear potentials in the HAL QCD method.

Parity-even channel with LO analysis at heavy pion masses
We first show the results of nuclear forces in the parity-even channel ( 1 S 0 and 3 S 1 -3 D 1 channels) at heavy quark masses obtained by the LO analysis for the derivative expansion of the potential. Since the statistical fluctuations are smaller at heavier quark masses in lattice QCD, this study is a good starting point to grasp the nature of lattice QCD nuclear forces. In addition, quark mass dependence of nuclear forces is of fundamental importance from a point of view of, e.g., anthropic principle, which cannot be studied by experiments.
In the case of 1 S 0 channel, we obtain the LO central force following Eq. (31). In the case of 3 S 1 -3 D 1 channel, the LO potentials consist of the central and tensor forces, which can be obtained from the coupled channel analysis between the Sand D-wave components as where ellipses represent higher order terms in the derivative expansion. Using the projection to the A + 1 representation of the cubic group (S-wave projection), P A + 1 , and the orthogonal one (D wave projection), (1 − P A + 1 ), the above equation reduces to two independent equations, from which V C (r) and V T (r) can be obtained [16]. Since the A + 1 representation couples to the angular momentum l = 0, 4, 6, · · · , these projections are expected to serve as the relevant partial wave decomposition at low energies. We find that the NBS correlation functions after P A + 1 and (1 − P A + 1 ) are dominated by S-wave and D-wave components, respectively, indicating that the contaminations from l ≥ 4 components are indeed small. For a more advanced partial wave decomposition, see [43]. We perform the calculations in quenched [19,16], dynamical 2-flavor [44], dynamical 3flavor [45,46,37] and dynamical (2+1)-flavor [47,27,38,40] lattice QCD with various quark masses. We here present the results obtained in 3-flavor lattice QCD at (M ps , M oct ) = (1171, 2274), (1015, 2031), (837, 1749), (672, 1484), (469, 1161) MeV [45,46,37]. 5 In the case of (M ps , M oct ) = (837, 1749), the value of quark masses m u = m d = m s nearly correspond to the physical strange quark mass. We generate gauge configurations with the RG-improved Iwasaki gauge action and non-perturbatively O(a)-improved Wilson quark action on a L 3 × T = 32 3 × 32 lattice. The lattice spacing is a = 0.121(2) fm and hence lattice size L is 3.87 fm. In the calculation of the NBS correlation function, parity-even states are created by a two-baryon operator with a wall quark source, while a point operator is employed for each baryon at the sink. Fig. 12 (Upper) are the lattice QCD results for the potentials. We find that the results are insensitive to the Euclidean time t, at which the NBS correlation function is evaluated, indicating that the derivative expansion is well converged. The obtained potentials are found to reproduce the qualitative features of the phenomenological N N potentials, namely, attractive wells at long and medium distances, central repulsive cores at short distance and strong tensor force with a negative sign. We also find intriguing features in the quark mass dependence of the potentials. At long distances, it is observed that the ranges of the tail structures in the central and tensor forces become longer at lighter quark masses. Such a behavior can be understood from the viewpoint of one-boson-exchange potential. At short distances, the repulsive cores in the central forces are found to be enhanced at lighter quark masses. This could be explained by the short-range repulsion due to the one-gluon-exchange in the quark model, whose strength is proportional to the inverse of the (constituent) quark mass. In fact, our systematic studies including hyperon forces with the same lattice setup revealed that the nature of repulsive core is well described by the quark Pauli blocking effect together with the one-gluon-exchange effect [45,37,48].

Shown in
As noted before, the potentials themselves are not physical observables and quantitative lattice QCD predictions shall be given in terms of scattering observables. Shown in Fig. 12 (Lower) are the scattering phase shifts (and mixing angles) obtained from lattice nuclear forces. We find that N N systems do not bound at these pseudoscalar masses as discussed in Sec. 3. Behaviors of phase shifts are qualitatively similar to the experimental ones, while the strength of the attraction is weaker due to the heavy quark masses in this calculation. It is also observed that quark mass dependence of phase shifts is quite non-trivial. In fact, if we decrease the quark masses, there appear competing effects in the interaction: the long-range attraction becomes stronger and the short-range repulsive core also becomes stronger. We also note that lighter quark masses correspond to lighter nucleon mass, which leads to larger kinetic energies.
We also present the results obtained in (2+1)-flavor lattice QCD at quark masses corresponding to (m π , m N ) (701, 1584), (570, 1412) and (411,1215) MeV [38]. Note that only up and down quark masses are varied with a strange quark mass being fixed to the physical value in this study. We employ the gauge configurations generated by the PACS-CS Collaboration with the RG-improved Iwasaki gauge action and non-perturbatively O(a)-improved Wilson quark action on a L 3 × T = 32 3 × 64 lattice. The lattice spacing is a 0.091 fm (a −1 = 2.16(31)GeV), which leads to the spatial extension L 2.9 fm.
In Fig. 13, we show the lattice QCD results for the potentials in the 1 S 0 and 3 S 1 -3 D 1 channels, together with the corresponding phase shifts in the 1 S 0 channel. Qualitative features are similar to those in 3-flavor case: (i) the central forces have repulsive cores at short distance and attractive wells at long and medium distances, both of which are enhanced at lighter quark masses (ii) the tensor force is strong with a negative sign, which increases at lighter quark masses.

More structures: spin-orbit forces in the parity-odd channel and three nucleon forces
If we consider an interaction at higher order terms in the derivative expansion, there appear more structures in the potentials. In particular, the extension from LO analysis to NLO analysis enables us to determine the spin-orbit (LS) force. The LS force is known to play an important role in the LS-splittings of nuclear spectra and the nuclear magic numbers. In addition, the LS force in the 3 P 2 -3 F 2 channel attracts great interest in nuclear astrophysics, since it could lead to the P -wave superfluidity in the neutron stars and affect the cooling process of neutron stars.
We here present the calculation in parity-odd channels ( 1 P 1 , 3 P 0 , 3 P 1 , 3 P 2 -3 F 2 channels) at heavy quark masses and show the results of LS forces as well as central/tensor forces [44]. In order to construct the source operator which couples to parity-odd states, we employ the two nucleon operators as  [38].
where N denotes a nucleon operator with a momentum, with f (±j) ( x) ≡ exp (±2πix j /L). A cubic group analysis shows that this source operator contains the orbital contribution T − 1 ⊕ A + 1 ⊕ E + , whose dominant components have l = 1, 0, 2, respectively, and thus covers all the two-nucleon channels with J ≤ 2. Combined with the spin degrees of freedom, we consider the T − 1 representation in the spin singlet channel and the A − representations in the spin triplet channel. At low energies, these representations correspond to the 1 P 1 channel and the 3 P 0 , 3 P 1 and 3 P 2 -3 F 2 channels, respectively, from which we extract the central force in the spin singlet channel (V I=0 C,S=0 (r)), and the central, tensor and LS forces (V I=1 C,S=1 (r), V I=1 T (r), V I=1 LS (r)) in the spin triplet channel.
Calculations are performed in 2-flavor lattice QCD at quark masses corresponding to (m π , m N ) (1133, 2158) MeV [44]. We employ the gauge configurations generated by the CP-PACS Collaboration with the RG-improved Iwasaki gauge action and a mean field O(a)-improved Wilson quark action on a 16 3 × 32 lattice. The lattice spacing a = 0.156(2) fm leads to the spatial extension L 2.5 fm. Fig. 14  (ii) the tensor force V I=1 T (r) is positive and weak compared to V I=1 C;S=1 (r) and V I=1 LS (r), and (iii) the LS force V I=1 LS (r) is negative and strong. These features are qualitatively in line well with those of the phenomenological potential. One can also see these properties in terms of the potential in each channel. In Fig. 14 (Upper-Right), we plot the potentials in the 1 P 1 , 3 P 0 , 3 P 1 and 3 P 2 channels, which are defined by V (r;

Shown in
To obtain the scattering observables, we fit the potentials and solve the Schrödinger equation in the infinite volume. In Fig. 14 (Lower), we show the results for the scattering phase shifts. Compared with the experimental phase shifts, we find that behaviors of phase shifts are generally well reproduced, while the magnitudes are smaller due to the heavier pion mass in lattice QCD calculations. In the 3 P 0 channel, we observe that the attraction is missing compared with the experimental one, which however is also likely due to the weak tensor force V T caused by the heavier pion mass. Among others, the most interesting feature is the attraction in the 3 P 2 channel as shown in Fig.14 (Lower-Right), originated from the strong (and negative) LS forces. As noted before, it is this interaction which is relevant to the paring correlation of the neutrons and possible P -wave superfluidity in the neutron stars.
We now turn to the study of three-nucleon forces. Determination of three-nucleon forces is one of the most challenging problems in nuclear physics: Three-nucleon forces are known to play important role in nuclear spectra/structures such as the binding energies of (light) nuclei and properties of neutron-rich nuclei. They are also essential ingredients to understand properties of nuclear matters such as the equation of state (EoS) at high density, which is relevant to the structures of neutron stars and nucleosynthesis at the binary neutron star mergers. While there have been many studies to construct three-nucleon forces by phenomenological approaches [49,50] or by chiral EFT approaches [51,6,7,8], it is most desirable to carry out the direct determination from QCD.
To study three-nucleon forces in lattice QCD, we consider the NBS wave function for a n(≥ 3)-particle system, |α , where W α is the center of mass energy of the system and we ignore the spins of nucleon for simplicity.
In [52,53,54], we show that the asymptotic behavior of the NBS wave function with the non-relativistic approximation can be written as where D = 3(n − 1) is the dimension of a n-particle system, is the radial component of the NBS wave function in D-dimension with R and Q being the hyper radius and momentum, respectively, and [L], [K] denotes the quantum numbers of the angular momentum in D-dimension. δ [N ] (Q) is the generalized "phase shift" for a n-particle system and U [L][N ] (Q) is a unitary matrix, which parameterize the T -matrix as Therefore, as in the case of n = 2 system (See Sec. 2.2.1), the information of T -matrix is encoded in the asymptotic behavior of the NBS wave function. Based on this property, we can define the energyindependent non-local potential for a n-particle system, which can be extracted from the (time-dependent) HAL QCD method.
We calculate the six-point correlation function divided by two-point correlation function cubed, where R ≡ ( are the Jacobi coordinates. In the timedependent HAL QCD method at the LO analysis for the derivative expansion and with the non-relativistic approximation, we can extract the three-nucleon forces V 3N F ( r, ρ) through the following Schrödinger where V 2N ( r ij ) with r ij ≡ x i − x j denotes two-nucleon forces between (i, j)-pair, µ r = m N /2, µ ρ = 2m N /3 the reduced masses.
In our first study of three-nucleon forces, we consider the total 3N quantum numbers of (I, J P ) = (1/2, 1/2 + ), the triton channel. We also consider a particular spacial geometry of the 3N, i.e., the "linear setup" ( ρ = 0), where 3N are aligned linearly with equal spacing of r 2 ≡ | r|/2. This setup makes the analysis much simpler. In addition, we consider the following channel, ψ S ≡ 1 √ 6 − p ↑ n ↑ n ↓ + p ↑ n ↓ n ↑ − n ↑ n ↓ p ↑ + n ↓ n ↑ p ↑ + n ↑ p ↑ n ↓ − n ↓ p ↑ n ↑ , and calculate the corresponding matrix element of V 3N F , so that we can suppress the statistical fluctuations in subtracting the contribution from V 2N .
One of the biggest challenges in the lattice QCD study of three-nucleon forces is the enormous computational cost required for the calculation of correlation functions. In fact, in terms of a mass number A, the cost grows with the multiplication of two factors, one of which scales factorially in A due to the Wick contraction (permutation of quarks), and the other of which scales exponentially in A due to the color/spinor contractions. On this point, we have developed a novel computational algorithm, called the unified contraction algorithm (UCA), in which two contractions are unified and redundant calculations are eliminated systematically [55]. In particular, the computation becomes faster by a factor of 192 for a calculation of three-nucleon forces.
We perform the calculation in 2-flavor lattice QCD at (m π , m N ) = (0.76, 1.81), (0.93, 1.85), (1.13, 2.15) GeV [56]. We employ the gauge configurations generated by CP-PACS Collaboration with mean field O(a)-improved Wilson fermion and RG-improved Iwasaki gauge action on a L 3 × T = 16 3 × 32 lattice. The lattice spacing is a = 0.1555 (17) fm and thus L = 2.5 fm. Shown in Fig. 15 (Left) are the lattice QCD results for the three-nucleon forces. We find a repulsive interaction at short-distances, r 2 0.2-0.4 fm. (Results at r 2 0.2 fm would suffer from lattice discretization error.) Note that a repulsive short-range three-nucleon force is phenomenologically required to explain the properties of high density matter. On the other hand, three nucleon forces are found to be suppressed at long distances. This is in accordance with the suppression of two-pion-exchange due to the heavier pion masses. Fig. 15 (Right) is the latest preliminary result obtained at m π = 510 MeV. In this calculation, we employ (2+1)-flavor lattice QCD gauge configurations generated in [29] with the RG-improved Iwasaki gauge action and non-perturbatively O(a)-improved Wilson quark action on a L 3 × T = 64 3 × 64 lattice (work in progress). The lattice spacing is a = 0.090 fm and L = 5.8 fm. Avoiding the very short-distance region where lattice discretization error could affect the results, we again find the short-range repulsive three-nucleon forces at r 2 0.2-0.7 fm. We find that, while the pion mass dependence of three-nucleon forces is not significant at m π = 0.76-1.13 GeV, the range of repulsive three-nucleon forces tend to be enlarged at m π = 0.51 GeV. It is important to pursue the study at lighter pion masses towards the physical pion mass.

Applications to nuclei, nuclear equation of state and structure of neutron stars
Once nuclear potentials are obtained by lattice QCD, we can use them to study various phenomena in nuclear physics and astrophysics. We here present the study of nuclear spectra/structures and Equation of State (EoS) of dense matter relevant to neutron star physics. Potentials used in this subsection are of the leading order only, and therefore are all hermitian. We can make non-hermitian higher order potentials in the HAL QCD method hermitian in the derivative expansion [57], which may be used for future applications in nuclear many body calculations.
In [58], binding energies and structures of doubly magic nuclei, 4 He, 16 O, 40 Ca, are studied by an ab initio nuclear many-body calculation based on lattice nuclear forces. We employ the nuclear forces obtained in 3-flavor lattice QCD at M ps = 469 MeV (See Fig. 12). We consider two-body nuclear forces in 1 S 0 , 3 S 1 and 3 D 1 channels, while nuclear forces in other channels and spin-orbit forces as well as three-nucleon forces are neglected. For simplicity, the Coulomb force between protons is not taken into account, either. As the ab initio many-body calculation, we employ self-consistent Green's function (SCGF) method, in which the single-particle propagator (Green's function) and the self-energy is solved self-consistently in a nonperturbative manner. In a practical calculation, the self-energy is calculated by Algebraic Diagrammatic Construction (ADC) formalism at third order for the so-called (low-momentum) P -space and Bethe-Goldstone equation (BGE) for the Q = 1 − P space. (see [58] for details.) In Tab. 2, we summarize the results for the ground state energies, together with the results from Brueckner Hartree-Fock (BHF) calculation [59] and exact stochastic variational calculation [60] using the same lattice nuclear forces. For the results from SCGF, the first parentheses show the errors associated with the infrared (IR) extrapolation in the SCGF calculation. We also estimate the errors from many-body truncations using 4 He as a benchmark. Since the SCGF result deviates from the exact solution by less than 10% for 4 He, and the SCGF approach is size extensive, we take a conservative estimate of 10% error for 16 O and 40 Ca, which are quoted in the second parentheses. The BHF results are sensibly more bound than the SCGF results, and we interpret this as a limitation of BHF theory. For the results shown in Tab. 2, there exist additional errors associated with the statistical fluctuations in the input lattice nuclear forces, which are estimated to be ∼ 10% [59]. Note that statistical fluctuations are correlated among nuclei, so we expect our observations described below are rather robust against statistical errors.
We find that at M ps =469 MeV in the SU(3) limit of QCD, both 4 He and 40 Ca have bound ground states while the deuteron is unbound. 16 O is likely to decay into four separate alpha particles, though it is already close to become bound. Moreover, we find that asymmetric isotopes, like 28 O, are strongly unbound systems. These results suggest that, when lowering the pion mass toward its physical value, closed shell  Table 2. Ground state energies of 4 He, 16 O and 40 Ca calculated by self-consistent Green's function (SCGF) method using nuclear forces at M P S =469 MeV obtained from 3-flavor lattice QCD with the HAL QCD method. Comparison is given with those obtained with BHF [59] and the exact calculation [60]. The last line is the breakup energy for splitting the system in 4 He clusters (of total energy A/4×5.09 MeV). Taken from [58].   [61,62]. For charge radii, we assumed the physical charge distributions of the nucleons. Taken from [58].
isotopes are created at first around the traditional magic numbers and the region of M ps ∼ 500 MeV marks a transition between an unbound nuclear chart and the emergence of bound isotopes.
We calculate the root mean square radii, which are given in Tab. 3, where we show only the central values. Although the total binding energies are 15-20% of the experimental value (Tab. 2), the computed charge radii are about the same as the experiment. We also find that the calculated one-nucleon spectral distributions are qualitatively close to those of real nuclei even for M ps =469 MeV considered here. This is due to the fact that the heavy nucleon mass (m N =1161.1 MeV) used here reduces the motion of the nucleons inside the nuclei and counterbalances the effect of weak attraction of the lattice nuclear forces at this pion mass.
We next present the study of properties of dense matter, namely, Equation of State (EoS) of nuclear matter. We again employ the nuclear forces in 1 S 0 , 3 S 1 and 3 D 1 channels obtained in 3-flavor lattice QCD. To study the quark mass dependence, we use lattice results for all five quark masses, at M ps = 469, 672, 837, 1015, 1171 MeV, which are shown in Fig. 12. As a method for a many-body calculation, we employ the Brueckner-Hartree-Fock (BHF) theory [64], which is known to be quantitative enough to give the essential underlying physics for infinite matter.
In Fig. 16 (Upper), we show the results of the ground state energy per nucleon (E/A) as a function of the Fermi momentum k F for the symmetric nuclear matter and the pure neutron matter. Shown together are the so-called APR EoS [63], which are obtained by the variational chain summation method from phenomenological nuclear forces with (APR(Full)) and without (APR(AV18)) three-nucleon forces. In Fig. 16 (Upper-Left), we find that the symmetric nuclear matter becomes a self-bound system with a saturation point (k F , E/A) MeV). This is the first time that the saturation in the symmetric nuclear matter is obtained through firstprinciples lattice QCD simulations. The saturation point, however, deviates from the empirical point primarily due to heavy pion (pseudo-scalar meson) mass in lattice simulation and the lack of three-nucleon forces in BHF calculation.
We also find a nontrivial M ps dependence of the EoS: the saturation disappears at intermediate pion masses (M ps = 672, 837 MeV) and possibly appears again at the heavy pion mass region (M ps = 1015, 1171 MeV). This implies that the saturation originates from a subtle balance between short-range repulsion and the intermediate attraction of the nuclear force, which have different m q dependence [37]. A similar nontrivial M ps dependence originated from the balance between repulsion and attraction is also observed for N N scattering phase shifts, as was discussed in Sec. 4.1.
In Fig. 16 (Upper-Right), we find that neutron matter is not self-bound due to large Fermi energy. If we decrease the pion mass, EoS is found to become stiffer. To further study the impact on phenomena in nuclear astrophysics, we calculate the mass (M ) versus the radius (R) relation of neutron stars at each pion mass. Here, we solve the Tolman-Oppenheimer-Volkoff (TOV) equation by using the EoS of the neutron-star matter with neutron, proton, electron and muon under the charge neutrality and beta equilibrium, where we use the standard parabolic approximation for asymmetric nuclear matters. Fig. 16 (Lower) is the M -R relation of the neutron star for different pion masses. As M ps decreases, the M -R curve shifts to the upper right direction, due to the stiffening of the EoS. While the maximum mass of the neutron star (M max ) in this calculation is much smaller than the recent observations, M max 2M , the deviation is most likely due to the heavy pion masses and lack of interactions as three-nucleon forces. A naive extrapolation of M max and the corresponding radius to M ps = 137 MeV would give M max ∼ 2.2M and R ∼ 12 km, which are encouraging for more quantitative studies in future. Another hottest topic in the context of neutron star physics is the effect of hyperon on the EoS at high density (so-called "hyperon puzzle"). Lattice QCD can play an unique role to study this effect by determining the hyperon forces which suffer from large uncertainties in experiments to date. For the on-going study in this direction, see [65].

Challenge: nuclear forces near the physical pion mass
While the results of nuclear forces at heavy pion masses are very intriguing and useful to extract the physical picture of nuclear forces, the quantitative results require the study at the physical pion mass. Note that the pion mass dependence of nuclear forces is quite non-trivial as discussed in Secs. 4.1 and 4.3, so the direct calculation near the physical point is desirable.
To this end, we have recently performed the first calculation of nuclear forces near the physical up, down and strange quark masses. Actually, our aim is to calculate not only nucleon forces but also hyperon forces, hereby achieve the comprehensive determination of two-baryon interactions from the strangeness S = 0 to −6 in parity-even channels (S-and D-waves). As mentioned before, the statistical fluctuations in lattice QCD are smaller (larger) for larger (smaller) quark masses, and thus the results have better precision in sectors involving more number of strange quarks (larger strangeness |S|). On the other hand, experiments in such larger |S| sectors are more difficult due to the short life time of hyperons. Therefore, lattice QCD studies and experiments are complementary with each other in the determination of baryon forces (See Fig. 17).
(2+1)-flavor gauge configurations are generated on a L 3 × T = 96 3 × 96 lattice with the RG-improved Iwasaki gauge action and non-perturbatively O(a)-improved Wilson quark action and APE stout smearing. The lattice spacing is a 0.0846 fm (a −1 2.333 GeV), so that spatial extent, L = 8.1 fm, is sufficiently large to accommodate two baryons in a box. Quark masses are tuned so as to be near the physical point, and the hadron masses are found to be (m π , m K , m N ) (146, 525, 955) MeV. NBS correlation functions for two-baryon systems are calculated for 56 channels in total and we extract the central and tensor forces in parity-even channel at the LO analysis for the derivative expansion (work in progress, and see also [66]). In order to make this first calculation a reality, "trinity" of state-of-the-art developments was crucial: (a) time-dependent HAL QCD method (theory), (b) unified contraction algorithm (software) and (c) K-computer, HOKUSAI and HA-PACS supercomputers (hardware). Fig. 18 are the results for the central force in the 1 S 0 channel (Left), and the central force (Middle) and tensor force (Right) in the 3 S 1 -3 D 1 channel. As noted above, nuclear forces are the most challenging interactions in lattice QCD calculation, and one can see that the results suffer from large statistical fluctuations. Nevertheless, the obtained results exhibit several interesting properties.

Shown in
First of all, the repulsive core at short-range is clearly observed in central forces. In order to clarify the physical picture for the repulsive core, it is useful to compare them with hyperon forces obtained in the same lattice setup. We find that the strength of repulsive core (or attractive core) highly depends on the flavor SU(3) (SU(3) f ) classification, in a consistent way with the quark Pauli blocking effect. In addition, if we compare interactions which belong to the same SU(3) f classification, such as N N ( 1 S 0 ) and ΞΞ( 1 S 0 ) both of which belong to 27-plet, we find that the strength differs in a way which can be understood from the viewpoint of one-gluon-exchange (e.g., repulsive core in N N ( 1 S 0 ) is stronger than that in ΞΞ( 1 S 0 )). These observations confirm the physical picture for the repulsive core obtained in the 3-flavor lattice QCD (Sec. 4.1), the quark Pauli blocking effect and the one-gluon-exchange, is relevant even at physical quark masses. See also [67] for a more detailed study on this point.
At middle and long distances, while statistical errors are quite large, we observe that the central force is attractive, resembling the phenomenological potential as one-pion-exchange potential (OPEP). The tensor force has relatively smaller statistical errors than the central forces, showing that the tensor force becomes stronger (with a negative sign) and has a longer tail, as compared with the tensor forces at heavier pion masses (Sec. 4.1). This property can be understood by the picture of OPEP. These results are encouraging and serve as the first step to establish a direct connection between QCD and nuclear physics. At the same time, statistical errors remain to be large and there also exist systematic errors associated with inelastic state contaminations. The studies to resolve these issues are in progress, and the second generation calculation is planned on the forthcoming Exascale computer, "Fugaku' (See https://postk-web.r-ccs.riken.jp/).

DIBARYONS
Before closing this review, we present our latest results on dibaryon searches in lattice QCD near the physical pion mass [66]. A dibaryon, a bound-state (or a resonance) with a baryon number B = 2 in QCD, can be classified in the SU(3) f as for the octet-octet system, where the deuteron, the only stable dibaryon observed in nature so far, appears in the 10 representation while H dibaryon has been predicted in the 1 representation [68] and actively investigated in lattice QCD [45,46,37,69,36]. For the decuplet-octet system, the classification leads to and N Ω ( N ∆) dibaryon has been predicted in the 8 (27) representation [70,71,72], and for the decuplet-decuplet system, where ΩΩ dibaryon has been predicted in the 28 representation [73] while ∆∆ has been predicted in the 10 [72,74] and the corresponding d * (2380) was indeed observed [75]. Note that among decuplet baryons, only Ω is stable against strong decays.

The most strange dibaryon
We first consider the ΩΩ system in the 28 representation of SU(3) f in the 1 S 0 channel [76]. Fig. 19 (Upper-Left) shows ΩΩ potentials at t/a = 16, 17, 18, which has qualitative features similar to the central potentials for N N but whose repulsion is weaker and attraction is shorter-ranged. This potential predicts an existence of one shallow bound state, whose binding energy is plotted in Fig. 19 (Upper-Right) as a function of the root-mean-square distance, with (red) and without (blue) Coulomb repulsion between ΩΩ. We may call this ΩΩ bound state "the most strange dibaryon". Such a system can be best searched experimentally by two-particle correlations in relativistic heavy-ion collisions [78].

N Ω dibaryon
We next consider the N Ω system with S = −3 in the 8 representation in the 5 S 2 channel [77]. Near the physical point, N Ω( 5 S 2 ) may couple to D-wave octet-octet channels below the N Ω threshold (ΛΞ and ΣΞ), but such couplings are assumed to be small in this calculation. Fig. 19 (Lower-Left) shows the N Ω potential at t/a = 11-14, which is attractive at all distances without repulsive core, so that one bound state appears in this channel. In Fig. 19 (Lower-Right), the binding energy (vertical) and the the root-mean-square distance (horizontal) are plotted for nΩ − with no Coulomb interaction (red) and pΩ − with Coulomb attraction (blue). These binding energies are much smaller than B = 18.9(5.0)( +12.1 −1.8 ) MeV at heavy pion mass m π = 875 MeV [79]. Such a N Ω state can be searched through two-particle correlations in relativistic nucleus-nucleus collisions [78] and an experimental indication was also reported [80].

Comparison among dibaryons
Let us consider the scattering length a 0 and the effective range r eff for ΩΩ( 1 S 0 ) and N Ω( 5 S 2 ). In Fig. 20, the ratio r eff /a 0 as a function of r eff are plotted for ΩΩ( 1 S 0 ) and N Ω( 5 S 2 ) obtained in lattice QCD near the physical pion mass, together with the experimental values for N N ( 3 S 1 ) (deuteron) and N N ( 1 S 0 ) (di-neutron). Small values of |r eff /a 0 | in all cases indicate that these systems are located close to the unitary limit.

CONCLUSIONS
In this paper, we have reviewed the recent progress in lattice QCD study of baryon-baryon interactions by the HAL QCD method. We first presented the detailed account on how to define the potentials in quantum field theories such as QCD. The key observation is that the Nambu-Bethe-Salpeter (NBS) wave functions contain the information of scattering phase shifts below inelastic threshold in their asymptotic behaviors outside the range of the interactions. The potentials at the interaction region can then be defined through the NBS wave functions so as to be faithful to the phase shifts by construction, where the non-locality of the potential is defined by the derivative expansion. In addition, by constructing the potentials in energy-independent way, the potentials can be extracted from two-baryon correlation functions without the requirement of the ground state saturation.
We then made a detailed comparison between the HAL QCD method and the conventional method, in which phase shifts are obtained from the finite volume energies through the Lüscher's formula. We pointed out that, while the validity of the latter method relies on the ground state saturation of the correlation function, its practical procedure for multi-baryon systems ("direct method") so far has utilized only the plateau-like structures of the effective energies at Euclidean times much earlier than the inverse of the lowest excitation energy. We showed theoretical and numerical evidences that such a procedure generally leads to unreliable results due to the contaminations from the elastic excited states: For instance, the results were found to be dependent on the operators and unphysical behaviors were exposed by the normality check. This invalidates the claim of the literature in the direct method that N N bound states exist for pion masses heavier than 300 MeV.
On the other hand, HAL QCD method is free from such a serious problem since the signal of potentials can be extracted from not only the ground state but also elastic excited states. While there instead exists the truncation error of the derivative expansion of the potential, the calculation of the higher order term in the derivative expansion showed that the convergence of the expansion is sufficiently good at low energies. Furthermore, utilizing the finite volume eigenmodes of the HAL QCD Hamiltonian, the excited state contaminations in the direct method were explicitly quantified. It turns out that the plateau-like structures of effective energies at early time slices are indeed pseudo-plateaux contaminated by elastic excited states and that the plateau for the ground state is realized only at a much larger time corresponding to the inverse of the lowest excitation energy gap. We also showed that, by employing an optimized operator utilizing the finite volume eigenmodes, the effective energies from the correlation functions give consistent results with those from the HAL QCD potential. Thus the long-standing issue on the consistency between the conventional method based on the Lüscher's formula and the HAL QCD method was positively resolved.
After establishing the reliability of the HAL QCD method, we presented the numerical results of nuclear forces from the HAL QCD method at various lattice QCD setups. At heavy pion masses, where good signal-to-noise ratio can be achieved in lattice QCD, we observed that the obtained N N potentials in the parity-even channel ( 1 S 0 , 3 S 1 -3 D 1 ) reproduce the qualitative features of the phenomenological potentials, namely, attractive wells at long and medium distances, accompanied with repulsive cores at short distance in the central potentials and the strong tensor force. The net interactions were found to be attractive, which however are not strong enough to form a bound N N state, probably due to the heavy pion masses. We observed that the tail structures are enhanced at lighter pion masses, which can be understood from the viewpoint of one-pion exchange contributions. We also found the repulsive cores are enhanced at lighter pion masses. Combined with our systematic studies including hyperon forces, the nature of repulsive cores was found to be well described by the quark Pauli blocking effect together with the one-gluon-exchange contribution.
The HAL QCD method can be extended to determine more complicated nuclear forces, such as spin-orbit forces and three-nucleon forces. In this paper, we considered two-nucleon systems in the parity-odd channels ( 1 P 1 , 3 P 0 , 3 P 1 , 3 P 2 -3 F 2 channels) and calculated spin-orbit forces as well as central and tensor forces. We found that qualitative features of experimental results are generally well reproduced, while the magnitudes differ due to the heavy pion mass. In particular, we observed the strong (and negative) spin-orbit forces, which lead to the attraction in the 3 P 2 channel. Three-nucleon forces were studied in the triton channel, (I, J P ) = (1/2, 1/2 + ), thank to the unified contraction algorithm (UCA), which can enormously speed up calculations of multi-baryon correlation functions. It was found that there exists a repulsive three-nucleon forces at short distances. These observations are of interest in the context of not only the structures of nuclei but also those of neutron stars, e.g., P -wave superfluidity and the maximum mass of neutron stars.
We carried out the applications to nuclei, nuclear equation of state (EoS) and structure of neutron stars based on lattice nuclear forces at heavy quark masses. We performed ab initio self-consistent Green's function (SCGF) calculations for closed shell nuclei with nuclear forces at M ps =469 MeV in the SU(3) limit of QCD. We found that 4 He, 40 Ca nuclei are bound and 16 O is close to become bound, while asymmetric isotopes are strongly unbound. The results suggest that, when lowering the pion mass toward its physical value, islands of stability appear at first around the traditional doubly magic numbers. The nuclear EoS was also studied by the BHF theory with nuclear forces in the flavor SU(3) limit. We found that the saturation property appears in the symmetric nuclear matter at M ps =469 MeV. A mass-radius relation of the neutron star was also studied based on the EoS obtained from lattice nuclear forces and we observed a tendency that the maximum mass of a neutron star increases as the pion mass decreases.
Finally, we presented the first lattice QCD study of baryon forces near the physical pion mass in the parity-even channel. The computation is quite challenging particularly for nuclear forces due to bad signal-to-noise ratio near the physical point. Nevertheless, we observed prominent characteristics of nuclear forces, such as repulsive cores at short distances as well as attractive interactions at mid and long distances in central forces, and a strong (and negative) tensor force. We also presented the results for the hyperon forces obtained near the physical point. We found that both ΩΩ( 1 S 0 ) and N Ω( 5 S 2 ) systems have strong attractions, and (quasi) bound dibaryons are formed near the unitary limit. These systems could be searched experimentally through two-particle correlations in relativistic nucleus-nucleus collisions.
Present results shown in this paper already indicate a clear pathway which connects nuclear physics with its underlying theory of the strong interaction, QCD. While there remain many challenges to accomplish researches in this direction, there is no doubt that successive theoretical developments together with next-generation supercomputers will further deepen the connection between the two. The outcome is also expected to play a crucial role to understand the nuclear astrophysical phenomena such as supernova explosions and mergers of binary neutron stars, as well as the nucleosynthesis associated with these explosive phenomena.