Lattice QCD and baryon-baryon interactions: HAL QCD method
LLattice QCD and baryon-baryon interactions:HAL QCD method
Sinya Aoki , , ∗ and Takumi Doi , Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, KyotoUniversity, Kyoto 606-8502 , Japan Quantum Hadron Physics Laboratory, RIKEN Nishina Center, Wako, 351-0198,Japan Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS),RIKEN, Wako, 351-0198, Japan
Correspondence*:Sinya [email protected]
ABSTRACT
In this article, we review the HAL QCD method to investigate baryon-baryon interactions such asnuclear forces in lattice QCD. We first explain our strategy in detail to investigate baryon-baryoninteractions by defining potentials in field theories such as QCD. We introduce the Nambu-Bethe-Salpeter (NBS) wave functions in QCD for two baryons below the inelastic threshold. We thendefine the potential from NBS wave functions in terms of the derivative expansion, which isshown to reproduce the scattering phase shifts correctly below the inelastic threshold. Using thisdefinition, we formulate a method to extract the potential in lattice QCD. Secondly, we discusspros and cons of the HAL QCD method, by comparing it with the conventional method, whereone directly extracts the scattering phase shifts from the finite volume energies through theL ¨uscher’s formula. We give several theoretical and numerical evidences that the conventionalmethod combined with the naive plateau fitting for the finite volume energies in the literature sofar fails to work on baryon-baryon interactions due to contaminations of elastic excited states. Onthe other hand, we show that such a serious problem can be avoided in the HAL QCD method bydefining the potential in an energy-independent way. We also discuss systematics of the HALQCD 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 (central) nuclearforce, tensor force, spin-orbital force, and three nucleon force. We finally show the latest resultscalculated at the nearly physical pion mass, m π (cid:39) MeV, including hyperon forces which leadto form ΩΩ and N Ω dibaryons. Keywords: lattice QCD, nuclear forces, baryon-baryon interactions, dibaryons, equation of state, neutron stars
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 fundamentaltheory of strong interaction, quantum chromodynamics (QCD). In particular, nuclear forces serve asone of the most basic constituents in nuclear physics, which are yet to be understood from QCD. While a r X i v : . [ h e p - l a t ] A ug inya Aoki et al. Lattice QCD and baryon-baryon interactions so-called realistic nuclear forces [1, 2, 3] have been established with a good precision, they are constructedphenomenologically based on scattering data experimentally obtained. Recent development in effectivefield theory (EFT) provides a more systematic approach for nuclear forces from a viewpoint of chiralsymmetry in QCD [4, 5, 6, 7, 8], whose unknown low-energy constants, however, cannot be determinedwithin 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 first-principles 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-bodytechniques developed in nuclear physics. The outcome is expected to make a significant impact on ourunderstanding of nuclear astrophysical phenomena such as supernovae, binary neutron star merges andnucleosynthesis.In this paper, we review the HAL QCD method to determine baryon-baryon interactions in latticeQCD. In this method, integral kernels, or so-called “potentials”, are first extracted from lattice QCD, andphysical observables such as scattering phase shifts and binding energies are calculated by solving theSchr¨odinger equation with obtained potentials in the infinite volume. We show that the notion of potentialcan be rigorously introduced as a representation of the S-matrix in quantum field theories as QCD. Theessential 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-localand energy-independent potential where the non-locality is defined through the derivative expansion. Inparticular, energy-independence of the potential is useful since one can extract the potential from the groundstate as well as elastic excited states simultaneously. This enables us to avoid the notorious signal-to-noiseissue for multi-baryon systems in lattice QCD (or the ground state saturation problem), and to make areliable determination of baryon-baryon interactions.In lattice QCD, there also exists a conventional method, in which phase shifts are obtained from finitevolume energies through the L ¨uscher’s formula. For meson-meson systems, a number of works have beenperformed based on the L ¨uscher’s formula [9], where finite volume energies are extracted utilizing thevariational method [10]. The L ¨uscher’s formula has been generalized for various systems, such as boostedsystems [11], arbitrary spin/partial waves [12, 13] and three-particle systems [14, 15]. While theoreticalbases are well established for both conventional method and HAL QCD method, numerical results forbaryon-baryon systems 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 systematic errors.In particular, we examine whether the systematic errors associated with excited state contaminationsare controlled or not in the procedure of the conventional method in the literature (“the direct method”),namely, simple plateau fitting for the ground state at early Euclidean times. We also examine systematicerrors in the HAL QCD method, in particular, the truncation error of the derivative expansion. We showtheoretical and numerical evidences that the inconsistency between two methods originates from excitedstate contaminations in the direct method. We also demonstrate that the inconsistency can be actuallyresolved if and only if finite energy spectra are properly obtained with an improved method rather than thenaive plateau fitting in the conventional method.After establishing the reliability of the HAL QCD method, we present the numerical results of nuclearforces from the HAL QCD method at various lattice QCD setups. The results at heavy pion masses forcentral and tensor forces are shown and their quark mass dependence as well as physical implicationsare discussed. The calculations of spin-orbit forces and three-nucleon forces are also given. Once nuclearforces are obtained, one can solve nuclear many-body systems with the obtained potentials. We study finite inya Aoki et al. Lattice QCD and baryon-baryon interactions nuclei, nuclear equation of state and structure of neutron stars based on lattice nuclear forces at heavy pionmasses. Finally, the latest results of nuclear forces near the physical pion mass are presented, as well ashyperon 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 interactionsfrom 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 demonstrationfor the newly developed method, the HAL QCD method. In Sec. 3, we discuss pros and cons of thesetwo methods, and compare the numerical results at heavy pion masses. We present evidences that theresults from the direct method suffer from uncontrolled systematic errors associated with the excited statecontaminations. In Sec. 4, we summarize results on nuclear potentials in the HAL QCD method. Afterreviewing the results obtained at heavy pion masses for central and tensor forces in the parity-even channelas well as spin-orbit forces and three-nucleon forces, we present nuclear many-body calculations based onlattice nuclear forces for double-magic nuclei, equation of state and the structure of neutron stars. Latestresults for nuclear forces near the physical pion mass are also given. In Sec. 5, we present hyperon forcesnear the physical pion mass, which lead to ΩΩ and N Ω dibaryons. Sec. 6 is devoted to the summary andconcluding remarks. In lattice QCD, the 2-pt function for a hadron H , created by O † H and annihilated by O H , is expressed as (cid:104) | O H ( (cid:126)p, t ) O † H ( (cid:126)p, | (cid:105) = ∞ (cid:88) n =0 Z n ( (cid:126)p ) e − E n ( (cid:126)p ) t + · · · , Z n ( (cid:126)p ) = |(cid:104) | O H ( (cid:126)p, | n, E n ( (cid:126)p ) (cid:105)| , (1)where | n, E n ( (cid:126)p ) (cid:105) is the n -th one-particle state with a mass m n , a momentum (cid:126)p and an energy E n ( (cid:126)p ) = (cid:112) m n + (cid:126)p , and ellipses represent contributions from multi particle states. We here assume m < m n> ,so that m is the hadron mass for the ground state, which can be extracted from the asymptotic behavior ofthe 2-pt function in the large t as (cid:104) | O H ( (cid:126)p, t ) O † H ( (cid:126)p, | (cid:105) (cid:39) Z ( (cid:126)p ) e − E ( (cid:126)p ) t + O (cid:16) e − E n> ( (cid:126)p ) t (cid:17) , t → ∞ , (2)where finite volume artifact is exponentially suppressed and can be eliminated by an infinite volumeextrapolation.So far, this method in lattice QCD (and the extension to lattice QCD+QED) has successfully reproducedlight hadron spectra [16] including the proton–neutron mass splitting [17]. A simple application of themethod, however, does not work well for an investigation of hadron interactions. For example, the 2-ptfunction of two baryons in the center of mass system behaves in the large t as (cid:104) | O BB ( (cid:126) , t ) O BB ( (cid:126) , † | (cid:105) (cid:39) Z BB e − E BB t + · · · , (3)where we obtain the lowest energy E BB . In the infinite volume limit, E BB behaves as E BB = 2 m B or E BB = 2 m B − ∆ E depending on an absence or presence of bound state. Here m B is the correspondingbaryon mass and ∆ E > is the binding energy of the lowest bound state. Only the binding energy ofthe bound state can be extracted by this simple method and thus more sophisticated methods are required. Frontiers 3 inya Aoki et al.
Lattice QCD and baryon-baryon interactions
Currently there are two methods to investigate hadron interactions in lattice QCD, the direct method (orfinite volume method) and the HAL QCD method, which are explained in this section.
The method most widely used to investigate hadron interactions in lattice QCD is to extract scatteringphase shifts from energy eigenvalues in 3-dimensional finite boxes through the L¨uscher’s finite volumeformula [18]. For example, in the case of the S -wave scattering phase shift, δ ( k ) , the formula reads k cot δ ( k ) = 1 πL (cid:88) (cid:126)n ∈ Z (cid:126)n − q , q = kL π , (4)where k is determined through E BB ( L ) = 2 (cid:113) k + m B with E BB ( L ) being the energy of the two baryonmeasured in lattice QCD on a finite box with the spatial extension L as in eq. (3). We here neglect thepartial wave mixing in the cubic group and spin degrees of freedom, for simplicity. Only the discrete setsof point ( k , k cot δ ( k )) , which satisfies eq. (4), are realized on a given volume L . Thus the scatteringphase shift δ ( k ) at the corresponding k can be extracted in lattice QCD, simply by measuring the finitevolume energy, E BB ( L ) . Note that the formula assumes that the hadron interaction is accommodatedwithin the lattice box and is not distorted by the finite volume artifact, which condition should be examinednumerically to be satisfied in actual calculations.In Fig. 1, we illustrate how scattering phase shifts and the bound state energy can be extracted by thismethod in the case of the N N scatterings. In the figure, the red solid line represents the effective rangeexpansion (ERE) for k cot δ ( k ) /m π at the Next-to-Leading order (NLO) as km π cot δ ( k ) = 1 a m π + r m π k m π (5)where the scattering length a and the effective range r are taken to be a m π = 16 . , r m π = 1 . for N N ( S ) (Left) or a m π = − . , r m π = 1 . for N N ( S ) (Right) with m π = 140 MeV, while coloreddashed lines represent the L¨uscher’s finite volume formula, eq. (4) on L = 10 , , , fm. Discretepoints which satisfy both the L¨uscher’s finite volume formula and the ERE are realized on each volume, asshown by the open squares, up/down triangles and diamonds.A distribution of the allowed k for k > becomes denser as the volume increases, so as to be continuousin the infinite volume limit, while a sequence of discrete points for k < leads to an accumulation point,which corresponds to the scattering state at k = 0 in the left figure or the bound state pole, denoted by theblack solid circle in the right figure. It is noted here that the bound state pole appears as the intersectionbetween the ERE and the bound state condition, − (cid:112) − ( k/m π ) (black solid line). To see this, we firstwrite k cot δ ( k ) = ik · S ( k ) + 1 S ( k ) − , S ( k ) = e iδ ( k ) , (6)where S ( k ) is the S-matrix for the N N elastic scattering. The bound state energy κ b can be extracted fromthe pole of this S-matrix as S ( k ∼ iκ b ) (cid:39) − iβ b k − iκ b , (7) inya Aoki et al. Lattice QCD and baryon-baryon interactions ( k/m π ) k c o t δ / m π NN( S ) finite volume ERE − q − ( k/m π ) L = 10 fm L = 12 fm L = 14 fm L = 18 fm ( k/m π ) k c o t δ / m π NN( S ) finite volume ERE − q − ( k/m π ) pole L = 10 fm L = 12 fm L = 14 fm L = 18 fm Figure 1.
A determination of k cot δ ( k ) /m π from energies of the two nucleon state in the finite volume.Taken from [19].where β b is real and positive for physical poles [20]. Thus at k (cid:39) − κ b , we have k cot δ ( k ) | k = iκ b = − κ b = − (cid:112) − k (cid:12)(cid:12)(cid:12) k = iκ b , (8)which means that the binding momentum k = iκ b is given by an intersection between k cot δ ( k ) and −√− k . Moreover, since ddk (cid:104) k cot δ ( k ) − ( − (cid:112) − k ) (cid:105)(cid:12)(cid:12)(cid:12)(cid:12) k = − κ b = − β b < , (9)the slope of k cot δ ( k ) must be smaller than that of −√− k as a function of k at the bound state pole,as in the case of Fig. 1 (Right). The finite volume analysis thus provides not only an infinite volumeextrapolation of the binding energy but also a novel way to examine the normality of the result in the directmethod [19]. The HAL QCD method, another method to investigate hadron interactions in lattice QCD, employs theequal time Nambu-Bethe-Salpeter (NBS) wave function, defined by φ k ( r ) e − W k t ≡ (cid:104) | N ( x + r , t ) N ( x , t ) | N N, W k (cid:105) , (10)where | N N, W k (cid:105) is the N N eigenstate in QCD with the center of mass energy W k = 2 (cid:113) k + m N andthe nucleon mass m N , and N ( x , t ) is a nucleon (annihilation) operator, made of quarks. Other quantumnumbers such as spin/isospin of two nucleons are suppressed for simplicity. We mainly use N α ( x ) = ε abc (cid:16) u a T ( x ) Cγ d b ( x ) (cid:17) q cα ( x ) , x ≡ ( x , t ) , (11) Frontiers 5 inya Aoki et al.
Lattice QCD and baryon-baryon interactions where C = γ γ is the charge conjugation matrix, q = u ( d ) for proton (neutron). Other choices such assmeared quarks are possible here, and such arbitrariness is considered to be a choice of the scheme forthe definition of the NBS wave function or the potential. (See [21] for such an example.) Throughout thispaper, we consider the N N elastic scattering, so that W k < W th ≡ m N + m π , where m π is the pionmass. 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 , beyondwhich ( i.e. r ≡ | r | > R ) the NBS wave function satisfies the Helmholtz equation as ( k + ∇ ) φ k ( r ) (cid:39) , k = | k | . (12)Furthermore, it behaves for large r > R as φ k ( r ) (cid:39) (cid:88) l,m Z l,m sin( kr − lπ/ δ l ( k )) kr Y lm (Ω r ) , (13)where Y lm is the spherical harmonic function for the solid angle Ω r of r , and we ignore spins of nucleonfor simplicity . 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 theS-matrix in QCD [24, 25].In the HAL QCD method, the non-local but energy-independent potential is defined from the NBS wavefunction through the following equation, ( E k − H ) φ k ( r ) = (cid:90) d r (cid:48) U ( r , r (cid:48) ) φ k ( r (cid:48) ) , E k = k m , H = − ∇ m , m = m N , (14)for W k < W th , and eq. (12) implies U ( r , r (cid:48) ) = 0 for r > R . While an existence of U ( r , r (cid:48) ) has been shownin [26, 23, 27], the non-local potential which satisfies eq. (14) is not unique. Thus we have to define thepotential uniquely, by specifying how to extract it. For this purpose, we introduce the derivative expansion, U ( r , r (cid:48) ) = V ( r , ∇ ) δ (3) ( r − r (cid:48) ) , whose lowest few orders for the N N with a given isospin channel arewritten as V ( r , ∇ ) = V ( r ) + V σ ( r )( σ · σ ) + V T ( r ) S (cid:124) (cid:123)(cid:122) (cid:125) LO + V LS ( r ) L · S (cid:124) (cid:123)(cid:122) (cid:125) NLO + O ( ∇ ) , (15)where V ( r ) is the central potential, V σ ( r ) is the spin dependent potential with σ i being the Pauli matrixacting on the spinor index of the i -th nucleon, V T ( r ) is the tensor potential with the tensor operator S = 3(ˆ r · σ )(ˆ r · σ ) − ( σ · σ ) ( ˆ r ≡ r /r ), and V LS ( r ) is the spin-orbit (LS) potential with the angularmomentum L = r × p and the total spin S = ( σ + σ ) / . It is noted that an expansion of the non-localpotential is not unique. For example, we may improve the convergence of the expansion by modifying the ∇ operator [28].Once we obtain the approximated potential at lowest few orders, we can calculate the scattering phaseshifts or the binding energies of possible bound states by solving the Schr¨odinger equation with thispotential in the infinite volume. As is the case for the finite volume method, it is necessary that the potentialis not distorted by the finite volume artifact, but this can be checked easily since the potential itself is The formula becomes more complicated if the nucleon spins are considered [22, 23]. inya Aoki et al. Lattice QCD and baryon-baryon interactions explicitly obtained. We can also check how good the approximated potential is, by increasing the order ofthe expansion. Needless to say, the approximated potential depends on momenta of input wave functions.As pointed out in [29], these dependences of the approximated potentials have been misidentified withthose of the non-local potential in the literature [30]. In the next subsection, we will explicitly demonstratehow this procedure works.
In order to see how the scattering phase shifts can be obtained by the HAL QCD method, we consider thequantum mechanics for a spinless system with a separable potential, defined by U ( r , r (cid:48) ) = ωv ( r ) v ( r (cid:48) ) , v ( r ) ≡ e − µr . (16)The S -wave solution of the Schr¨odinger equation with this potential is given exactly by φ k ( r ) = e iδ ( k ) kr (cid:20) sin { kr + δ ( k ) } − sin δ ( k ) e − µr (cid:18) r ( µ + k )2 µ (cid:19)(cid:21) , (17)where k cot δ ( k ) = − µ (cid:20) µ ( µ − k ) − µ + k µ ( µ + k ) + ( µ + k ) πmω (cid:21) , (18)which is the 4-th order polynomials in k . In order to make the scattering phase shift a more complicatedfunction of k , we artificially modify the wave function from φ k ( r ) to φ k ( r ) which is defined by φ k ( r ) = φ k ( r ) ( r ≤ R ) C ( k ) e iδ R ( k ) kr sin { kr + δ R ( k ) } ( r > R ) , (19)where R is an infrared cutoff, and it is understood that the potential is modified accordingly. The continuityof φ k ( r ) and φ (cid:48) k ( r ) at r = R gives k cot δ R ( k ) = k Y cot( kR ) + XX cot( kR ) − Y , X = φ k ( R ) , Y = ddr [ rφ k ( r )] (cid:12)(cid:12)(cid:12)(cid:12) r = R , (20)as well as C ( k ) = X/ sin( kR + δ R ( k )) . Hereafter, we study how the scattering phase shifts are obtainedin the HAL QCD method.The derivative expansion for the S -wave scatterings leads to V ( r , ∇ ) = V ( r ) + V ( r ) ∇ + O ( ∇ ) , (21)and we consider to determine the potential in each order from φ k ( r ) .The leading order (LO) potential is given by V LO ( r , ∇ ) = V LO0 ( r ; k ) = ( E k − H ) φ k ( r ) φ k ( r ) , (22) Frontiers 7 inya Aoki et al.
Lattice QCD and baryon-baryon interactions while the next-to-leading order (NLO) potential is extracted as V NLO ( r , ∇ ) = V NLO0 ( r ; k , k ) + V NLO1 ( r ; k , k ) ∇ , (23)where (cid:18) V NLO0 ( r ; k , k ) V NLO1 ( r ; k , k ) (cid:19) = 1 D ( r ; k , k ) (cid:18) m (cid:2) V LO0 ( r ; k ) E k − V LO0 ( r ; k ) E k (cid:3) V LO0 ( r ; k ) − V LO0 ( r ; k ) (cid:19) ,D ( r ; k , k ) = 2 m (cid:104) V LO0 ( r ; k ) − V LO0 ( r ; k ) − ( E k − E k ) (cid:105) . (24)Note that the potential in each order in the derivative expansion { V ( r ) , V ( r ) , · · · } are defined to be k -independent, while the potentials approximately obtained in each LO/NLO analysis, { V LO0 ( r ; k ) } and { V NLO0 ( r ; k , k ) , V NLO1 ( r ; k , k ) } , have implicit k -dependence due to the truncation error in thederivative expansion [29].We calculate S -wave scattering phase shifts corresponding to these approximated potentials, and comparethem 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 , where all (dimensionful) quantities are measured in unitsof µ . In this example, we take ω = − . µ , m = 3 . µ and R = 2 . /µ . In the figures, the exact phaseshift δ 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, theLO approximation at k = 0 (orange), exact at k = 0 by construction, gives a reasonable approximation atlow energies ( k (cid:39) ) but deviates from the exact one at high energies near k (cid:39) µ . On the other hand, theLO approximation at k = µ (green) becomes accurate at higher energies near k (cid:39) µ but inaccurate at lowenergies near k (cid:39) . Combining two NBS wave functions, φ k =0 ( r ) and φ k = µ ( r ) , one can determine theapproximated potential at the NLO, V NLO ( r , ∇ ) , whose scattering phase shifts are represented by the redsolid lines in the figures. The phase shifts at the NLO (red lines) gives reasonable approximations of theexact results (blue solid lines) in the whole range ( ≤ k ≤ µ ), as they are exact at k = 0 and k = µ by construction. If we increase the order of the expansion more and more, the approximation becomesbetter and better. Using this model, let us compare the direct method and the HAL QCD method. At the LO, the directmethod gives either k cot δ ( k ) at k = 0 or k = µ 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 theLO 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 = 0 and k = µ as k cot δ ( k ) = 1 a + r eff k , a = lim k → k cot δ ( k ) , r eff δ ( µ ) µ − µ a , (25)which is given by a red dashed line in the right figure. By comparing the HAL QCD method with potentialsat NLO (the red solid line) and the direct method with NLO ERE (the red dashed line), the former leadsto a better approximation of the exact result than the latter, since higher order effects in ERE in terms of k are included in the former. Note, however, that sufficiently precise data in the direct method can alsoevaluate higher order ERE terms than NLO, in principle. A similar attempt to represent an arbitrary potential in terms of a separable potential is given in [31, 32]. inya Aoki et al. Lattice QCD and baryon-baryon interactions
LO at k = 0 NLOexact
LO at k = µ NLO
LO at k = 0 exactNLO ERE LO at k = µ LO ERE ar k = 0 LO ERE at k = µ Figure 2.
The scattering phase shifts δ ( k ) and k cot δ ( k ) as a function of k . See the main text for moredetails. -2 -1 0 1 2 -2 -1 0 1 2 1 1.5 φ (x,y,z=0) PBCx[fm] y[fm] -2 -1 0 1 2 -2 -1 0 1 2 0 0.5 1 1.5 φ (x,y,z=0) APBCx[fm] y[fm] Figure 3.
The NBS wave function for
N N ( S ) at E k (cid:39) MeV with the PBC (Left) and at E k (cid:39) MeV with the APBC (Right). Both are normalized to unity at r = 1 fm. Taken from [33]. 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 inQCD.Fig. 3 shows three dimensional plots of the NBS wave functions φ k ( x, y, z = 0) for N N ( S ) with theperiodic boundary condition (PBC) at E k (cid:39) MeV (Left) and with the anti-periodic boundary condition(APBC) at E k (cid:39) MeV (Right), in quenched lattice QCD at a (cid:39) . fm on L (cid:39) . fm with m π (cid:39) MeV [33]. As seen from the figure, two NBS wave functions look very different from eachother. In particular, the right one vanishes on the boundary due to the APBC constraint.
Frontiers 9 inya Aoki et al.
Lattice QCD and baryon-baryon interactions V c , s ( L O ) (r) [ M e V ] r[fm] 45[MeV]0[MeV]-50 0 50 0 1 2 V c , s ( L O ) (r) [ M e V ] [fm] D S -50 0 50 100 0 0.5 1 1.5 2 Figure 4. (Left) The LO potential for
N N ( S ) as a function of r at E k (cid:39) MeV (red solid circles)and at E k (cid:39) MeV (blue open circles). (Right) The LO potential as a function of r at E k (cid:39) MeV for
N N ( S ) (red open circles) and for N N ( D ) (cyan solid circles). Taken from [33].Fig. 4 (Left) compares the LO potentials for N N ( S ) obtained from the corresponding NBS wavefunctions in Fig. 3. While the NBS wave functions at different energies have different spatial structures, thepotentials look very similar. This suggests that the higher order terms in the derivative expansion of thepotential have negligible contributions at this energy interval, ≤ E k ≤ MeV.Fig. 4 (Right) compares the LO potential for
N N ( S ) (red open circles) with the one for N N ( D ) (cyan solid circles) at E k (cid:39) MeV. Although statistical fluctuations are larger for the latter, they looksimilar, suggesting that L dependence of the potential is also small in this setup. If more accurate datashow a difference of potentials between N N ( S ) and N N ( D ) , one may determine the L dependentterm of the potential in the spin-singlet channel. In order to extract the NBS wave functions on the finite volume in lattice QCD, we consider the 4-ptfunction given by F J ( r , t − t ) = (cid:104) | N ( x + r , t ) N ( x , t ) ¯ J NN ( t ) | (cid:105) = (cid:88) n A Jn φ k n ( r ) e − W kn ( t − t ) + · · · , (26)where ¯ J NN ( t ) is an operator which creates two nucleon states at time t , A Jn ≡ (cid:104) N N, W k n | ¯ J NN (0) | (cid:105) ,and ellipses represent inelastic contributions, which become negligible at W th ( t − t ) (cid:29) . Like the directmethod, one can extract the NBS wave function for the ground state from the above 4-pt function as F J ( r , t ) (cid:39) A J φ k ( r ) e − W k t (27)for ( W k − W k ) t (cid:29) , where W k ( W k ) is the lowest (second-lowest) energy on the finite volume. TheLO potential from the NBS wave function for the ground state is then extracted from F J ( r , t ) at large t . Aswill be discussed in the next section, however, it is numerically very difficult to determine F J ( r , t ) for twonucleons at such large t due to the bad signal-to-noise (S/N) ratio. inya Aoki et al. Lattice QCD and baryon-baryon interactions
Fortunately, an alternative extraction is available for the HAL QCD method [34]. Let us consider theratio of 4-pt function to the 2-pt function squared as R J ( r , t ) ≡ F J ( r , t ) G N ( t ) , G N ( t ) = (cid:88) x (cid:104) | N ( x , t ) N ( , | (cid:105) (cid:39) Z N e − m N t + · · · , (28)which behaves R J ( r , t ) = (cid:88) n ˜ A Jn φ k n ( r ) e − ∆ W kn t , ˜ A Jn ≡ A Jn Z N , ∆ W k ≡ W k − m N , (29)for W th t (cid:29) , where inelastic contributions can be neglected. Noticing that ∆ W k = k m N − (∆ W k ) m N , (cid:18) k m N − H (cid:19) φ k ( r ) = V ( r , ∇ ) φ k ( r ) , (30)we obtain (cid:26) − H − ∂∂t + 14 m N ∂ ∂t (cid:27) R J ( r , t ) = V ( r , ∇ ) R J ( r , t ) . (31)We can approximately extract V ( r , ∇ ) from R J ( r , t ) for (different) J ’s, as long as t satisfies the conditionthat W th t (cid:29) (elastic state saturation), which is much easier than to achieve ( W k − W k ) t (cid:29) (groundstate saturation). We call this alternative extraction the time-dependent HAL QCD method. It is interesting to ask whether the attractions of the nuclear forces at low energies would become weakeror stronger if the pion mass were larger than the value in Nature. In principle, such a question can beanswered 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 summarizedin Table 1, the direct method tends to indicate that attractions between two nucleons become stronger asthe pion mass increases, so that both deuteron and di-neutron form bound states, while the HAL QCDmethod suggests that the attractions become weaker and the bound deuteron does not exist at heavierpion 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 [19]. 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 [40] gives 75(5)MeVat m π = 810 MeV, which is much larger than 19(10) MeV at m π = 960 MeV by Mainz2018 [43]. Onthe other hand, HAL2012 [44] gives 38(5) MeV at m π = 837 MeV from the HAL QCD method. Thesedeviations seem to be too large to be explained by lattice artifacts.In order to understand origins of these discrepancies, we have performed extensive investigations, whoseresults have been published in a series of papers [46, 19, 47, 48], which will be explained in the followingsubsections.
Frontiers 11 inya Aoki et al.
Lattice QCD and baryon-baryon interactions
Collaboration Ref. N f m π − ∆ E ( S ) − ∆ E ( S ) − ∆ E ( H ) The direct methodYKU2011 [35] 0 800 4.4(1.2) 7.5(1.0) —YIKU2012 [36] 2+1 510 7.4(1.4) 11.5(1.3) —NPL2015 [37] 2+1 450 12.5( +3 . − . ) 14.4( +3 . − . ) —NPL2012 [38] 2+1 390 7.1(9.0) 11(13) 13.2(4.4)YIKU2015 [39] 2+1 300 8.5( +1 . − . ) 14.5( +2 . − . ) —NPL2013 [40] 3 810 15.9(3.8) 19.5(4.8) 74.6(4.7)NPL2017 [41] 3 810 20.6( +3 . − . ) 27.9( +3 . − . ) —CalLat2017 [42] 3 810 21.8( +3 . − . ) 30.7( +2 . − . ) —3 8.35(1.1)* 3.3( +1 . − . ) —Mainz2018 [43] 3 †
960 0 — 19(10)2+1 †
440 — — 18.8(5.5)*The HAL QCD methodIAH2007 [26] 0 530 0 0 —AHI2009 [23] 0 380, 530, 730 0 0 —HAL2012 [44] 3 1171 0 0 49.1(6.5)3 1015 0 0 37.2(4.4)3 837 0 0 37.8(5.2)3 672 0 0 33.6(5.9)3 469 0 0 26.0(6.5)HAL2012a [34] 2+1 701 0 — —HAL2013 [45] 2+1 411, 570, 701 0 — —
Table 1.
Summary of binding energies [MeV] for
N N ( S ) , N N ( S ) and H -dibaryon in lattice QCD.NPL2013, NPL2017 and CalLat2017 employed the same set of gauge configurations. CalLat2017 foundtwo states in each channel. In Mainz2018, dynamical 2-flavor with quenched strange quark configurationsare employed and N f in the table (with † symbol) denotes the information in the valence quark sector. Allvalues of ∆ E correspond to those in the infinite volume limit except ones with ∗ , which are values on thefinite volumes. The number in ∆ E indicates the system is unbound in this channel. In the direct method, reliable extractions of the two nucleon ground state energies are crucially important.As long as ( W k − W k ) t (cid:29) , the two nucleon correlation function is dominated by the ground state as G NN ( t ) = (cid:104) | J NN ( t ) ¯ J (cid:48) NN (0) | (cid:105) (cid:39) Z Jk ¯ Z J (cid:48) k e − W k t , Z J ( J (cid:48) ) k ≡ (cid:104) | J NN ( J (cid:48) NN ) | N N, W k (cid:105) , (32)so that the extracted ground state energy W k depends neither the source operator ¯ J (cid:48) NN nor the sink operator J NN , while magnitudes of contaminations from excited states are affected by the choices of these operators.Since W k − W k (cid:39) (2 π/L ) /m N on the finite box with the spacial extension L , t (cid:29) fm is required,for example, for L (cid:39) m N (cid:39) 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 ∼ fm,by tuning the source operators ¯ J (cid:48) NN in order to achieve a plateau of the effective energy shift ∆ E eff NN ( t ) atsuch a small t , where ∆ E eff NN ( t ) = − a log R NN ( t + a ) R NN ( t ) , R NN ( t ) ≡ G NN ( t ) G N ( t ) , (33) inya Aoki et al. Lattice QCD and baryon-baryon interactions t [ a ]201510505 ∆ E e ff NN ( t ) [ M e V ] L = PACS-CSsmeared src. NR NN ( S )wall src. NR NN ( S )
10 11 12 13 14 15 16 t [ a ] ∆ E e ff ΞΞ ( t ) [ M e V ] smeared src. ΞΞ ( S ) ( g ( r ) = 1) smeared src. g ( r ) = 1 + 0 . − . smeared src. g ( r ) = 1 − . − . smeared src. g ( r ) = 1 − . − . Figure 5. (Left) The effective energy shift ∆ E eff NN ( t ) for N N ( S ) from the wall source (red circles)and the smeared source (blue squares) on L = 48 a (cid:39) . fm at m π = 0 . GeV, m N = 1 . GeV and m Ξ = 1 . GeV [46]. (Right) The effective energy shift ∆ E effΞΞ ( t ) for ΞΞ( S ) from the smeared sourcewith different sink operators on the same gauge configurations [46]. Figure 6. ∆ E eff NN ( t ) − ∆ E NN from the mockup data R mockup NN ( t ) with fluctuations and errors as a functionof t . (Left) b = 0 . , ± . , . and c = 0 . . (Right) c = 0 . , . , . and b = − . .Unfortunately, such a naive plateau fitting at earlier t may not be reliable due to contaminations fromnearby excited states, which may easily produce (incorrect) plateau-like behaviors in effective energies. Itwas indeed demonstrated that plateau-like behaviors in effective energy shifts at small t can depend notonly on the source operator but also on the sink operator: Plateaux disagree between the wall source (redcircle) and the smeared source (blue square) in the left of Fig. 5, while plateaux depend on sink operatorsfor the same smeared source in the right figure.In order to see how easily contaminations from elastic-excited states can produce plateau-like behaviorsat earlier t , let us consider the effective energy shift from the mockup data for R NN ( t ) , given by R mockup NN ( t ) = e − ∆ E NN t (cid:16) b e − δE el . t + c e − δE inel . t (cid:17) , (34)where we take δE el . = 50 MeV for the typical lowest elastic excitation energy on L (cid:39) fm at m N (cid:39) . GeV, and δE inel . (cid:39) m π (cid:39) MeV for the lowest inelastic energy. Naively, it is expected that the correct
Frontiers 13 inya Aoki et al.
Lattice QCD and baryon-baryon interactions plateau at ∆ E NN for the ground state appears at t (cid:29) /δE el . (cid:39) fm, which however is too large to havegood signals for two baryons such as N N . By tuning the source operator, one may reduce coefficients b and c . Since the N N operator does not strongly couple to
N N π state, we expect small c and take c = 0 . . On the other hand, N N operators easily couple to both ground and 1st elastic excited statesas they become almost identical to each other in the infinite volume limit. We therefore take b = 0 . (the highly tuned operator) , b = ± . (the tuned ones) as well as b = 0 . (the untuned one). Fig. 6(Left) shows ∆ E eff NN ( t ) for these 4 examples with c = 0 . , where random fluctuations and errors whosemagnitude increase exponentially in t are assigned to R mockup NN ( t ) . All examples show plateau-like behaviorsat t (cid:39) fm, but these four plateaux disagree with each other. As | b | increases, the deviation between thevalues of these “pseudo plateaux” and the true value becomes larger. Contaminations of the elastic excitedstates can easily produce the plateau-like behavior at earlier t , and the t dependence of data alone cannottell 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 shownin Fig. 6 (Right), where the effective energy shift for c = 0 . , . , . with b = − . is plotted. Allcases converge to almost the same pseudo plateau, while a pseudo plateau starts at later t for larger c . It isnoted that the multi-exponential fit does not work in this case at t (cid:39) . fm, which is much smaller thanthe necessary t (cid:29) /δE el . . The multi-exponential fit at such small t only separates the pseudo plateau fromthe inelastic contributions but is difficult to distinguish the ground state and the 1st excited state for theelastic states. ( k/m π ) k c o t δ / m π YIKU2012 NN( S ) − q − ( k/m π ) L/a = 32
L/a = 40
L/a = 48
L/a = 64
L/a = ∞ ERE ( k < ) ( k/m π ) k c o t δ / m π NPL2015 NN( S ) − q − ( k/m π ) L/a = 24
L/a = 32
L/a = 48
L/a = ∞ ERE (NPL2015)ERE ( k < ) Figure 7. (Left) k cot δ ( k ) /m π in YIKU2012 [36] for N N ( S ) as a function of ( k/m π ) . The solidred line and light red band represent the ERE fit and the corresponding error (statistical and systematicadded in quadrature), respectively. The dashed lines are the finite volume formula for the correspondingvolume. (Right) k cot δ ( k ) /m π in NPL2015 [37] for N N ( S ) as a function of ( k/m π ) . Two ERE fitsare performed depending on the lattice data to be used for the fit. The red line with the band represents thefit made by the authors in [19], while the blue line with the band is plotted by the authors in [19] using thefit result of NPL2015. Both figures are taken from [19].While the check through operator dependence is useful, it requires extra calculations. We find that thefinite volume formula in eq. (4) provides a simpler test, which tells us whether the ground state energiesextracted by the plateau fitting give a reasonable ERE or not without extra calculations. We call this test a inya Aoki et al. Lattice QCD and baryon-baryon interactions normality check [19]. Fig. 7 (Left) shows k cot δ ( k ) /m π in YIKU2012 [36] as a function of k /m π for N N ( S ) , where the solid red line represents the NLO ERE fit in eq. (5), and the light red bands showsstatistical and systematic errors added in quadrature [19]. Contrary to a naive expectation from non-singularERE behaviors, data align almost vertically, since ∆ E NN is almost independent of the volume. In otherwords, according to the finite volume formula, the claimed “binding energy” (open circle) is too shallowto have such volume independent ∆ E . Not only the central value of the NLO ERE fit gives singularparameters as (( a m π ) − , r m π ) = (5 . , . 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 operatordependence of these data, strongly indicate that the naive plateau fitting employed in the direct method isunreliable. Another example is shown in Fig. 7 (Right) for N N ( S ) from NPL2015 [37]. In this case, twodifferent NLO ERE fits (red line/band and blue line/band) are performed depending on the lattice data tobe used for the fit. It turns out that two ERE are inconsistent with each other, indicating that their latticedata themselves are “self-inconsistent”. In addition, one of ERE (blue line/band) is found to violate thephysical pole condition, eq. (9), at the crossing point (open circle). Similar symptoms are observed for allother data in the direct method claiming the existence of N N bound states at heavy quark masses [19]. The source operator dependence of the HAL QCD potential has been investigated in [47]. Fig. 8 (Left)compares the LO potentials, V LO0 ( r ) , for ΞΞ( S ) between the wall source (red open circles) and thesmeared source (blue open squares). We observe a small difference at short distances, from which onecan determine the N LO potential, V N LO ( r , ∇ ) = V N LO0 ( r ) + V N LO2 ( r ) ∇ . Note that the NLO term, V N LO1 ( r ) ∇ = V N LOLS ( r ) L · S is absent in the S channel. Fig. 8 (Right) shows V N LO2 ( r ) , which isnonzero only at r < . fm, where two LO potentials differ. We then extract the scattering phase shifts,using this N LO potential.The N LO corrections turn out to be negligible at low energies, as shown in Fig. 9 (Left), where k cot δ ( k ) is almost identical between V N LO ( r , ∇ ) (red solid circles) and V N LO0 ( r ) (blue solid squares).Furthermore, even the LO analysis for the wall source, V LO(wall)0 ( r ) (black open diamond), is sufficientlygood at low energies. As energy increases, the N LO corrections become visible as seen in Fig. 9 (Right),where ( k/m π ) = 0 . corresponds to ∆ E (cid:39) MeV for the energy shift from the threshold. It is notedthat V N LO0 ( r ) (blue solid squares) gives a little closer results to N LO results (red solid circles) than V LO(wall)0 ( r ) (black open diamond) does. In this subsection, we explain why the wall source and the smeared source give inconsistent plateaubehaviors, in the case of ΞΞ correlation functions as an example.To this end, we consider the Hamiltonian H = H + V LO(wall)0 , where we employ V LO(wall)0 ( r ) , theLO potential from the wall source, since it works rather well at low energies as shown in the previoussubsection. We first decompose R J ΞΞ ( r , t ) for J = wall/smear in terms of finite volume eigenfunctions of H as R J ΞΞ ( r , t ) = (cid:88) n a Jn ( t )Ψ n ( r ) e − ∆ E n t , a Jn ( t ) = (cid:88) r Ψ † n ( r ) R J ΞΞ ( r , t ) e ∆ E n t . (35) After these problems were pointed out in [19], revised data of NPL2013 have been presented in [41], whose EREs are still marginal to satisfy/violate thephysical pole condition.
Frontiers 15 inya Aoki et al.
Lattice QCD and baryon-baryon interactions
Figure 8. (Left) The LO potential, V LO0 ( r ) , for ΞΞ( S ) from the wall source (red open circles) and thesmeared source (blue open square). (Right) The second order term, V N LO2 ( r ) (blue solid squares), in theN LO potential V N LO ( r , ∇ ) = V N LO0 ( r ) + V N LO2 ( r ) ∇ for ΞΞ( S ) . Both are taken from [47]. ( k / m ) k c o t ( k ) / m , t = V LO(wall)0 V N LO0 V N LO0 + V N LO2 2 ( k / m ) ( k ) [ deg .], t = V LO(wall)0 V N LO0 V N LO0 + V N LO2 2
Figure 9. (Left) k cot δ ( k ) /m π as a function of ( k/m π ) at low energies, where δ ( k ) is the scatteringphase shift for ΞΞ( S ) , calculated from V N LO ( r , ∇ ) (red solid circles), V N LO0 ( r ) (blue solid squares)and V LO(wall)0 ( r ) (black open diamond). (Right) The corresponding δ ( k ) . Both are taken from [47].where Ψ n ( r ) and ∆ E n are normalized-eigenfunction and eigenenergy in the finite volume, respectively,and a Jn ( t ) is the overlapping coefficient extracted at t .Then the correlation function for the source J in the direct method is given by R J ΞΞ ( t ) = (cid:88) r R J ΞΞ ( r , t ) = (cid:88) n b Jn ( t ) e − ∆ E n t b Jn ( t ) = a Jn ( t ) (cid:88) r Ψ n ( r ) . (36) inya Aoki et al. Lattice QCD and baryon-baryon interactions t [ a ] ( S ) E e ff ( t ) [ M e V ] L = E E effwall E effsmear wall src.smeared src. t [ a ] ( S ) E e ff ( t ) [ M e V ] L = E E effwall E effsmear wall src.smeared src. Figure 10.
The reconstructed effective energy shift ∆ E J eff ( t, t = 13 a ) for the wall source (red bands)and the smeared source (blue bands) on L = 48 a , while the effective energy shifts directly from R J ΞΞ ( t ) are shown for J = wall (red open circles) and J = smear (blue open squares). The black dashed lines arethe energy shifts for the ground state of H in the finite box. (Left) ≤ t/a ≤ . (Right) ≤ t/a ≤ .Taken from [48].Finally, approximating a sum over n by the lowest few orders, we reconstruct the behavior of the effectiveenergy shift as a function of t as ∆ E J eff ( t, t ) = 1 a log (cid:18) R J ( t, t ) R J ( t + a, t ) (cid:19) , R J ( t, t ) = n max (cid:88) n =0 b Jn ( t ) e − ∆ E n t , (37)where we fix the overlapping coefficient b Jn ( t ) at t = t , and n max is a number of excited states used inthe approximation.In Fig. 10, we show reconstructed effective energy shift ∆ E J eff ( t, t = 13 a ) on L = 48 a with n max = 4 ,together with the effective energy shifts from R J ΞΞ ( t ) , for the wall source (red bands and red open circles)and the smeared source (blue bands and blue open squares). The black dashed line represents the energyshift for the ground state of H = H + V LO(wall)0 on L = 48 a .We find that the plateau-like structures in the direct method around t/a = 15 are well reproduced by ∆ E J eff ( t, t = 13 a ) for both sources in Fig. 10 (Left). This indicates that the plateau-like structures in thedirect 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 theground state. The fate of these structures is shown in Fig. 10 (Right), where we plot ∆ E J eff ( t, t =13 a ) at asymptotically large t . While the plateau-like structure for the wall source is almost unchanged, ∆ E J eff ( t, t = 13 a ) for the smeared source gradually increases and reaches to the true value at t/a ∼ .The above results clearly reveal that the plateau-like structures at t/a ∼ for the smeared source arepseudo-plateaux caused by the contaminations of the excited states. Large contaminations from excitedstates in the case of the smeared source are not caused by the smearing, but are indeed implied by putting Frontiers 17 inya Aoki et al.
Lattice QCD and baryon-baryon interactions
10 11 12 13 14 15 t [ a ] E e ff ( t ) [ M e V ] L = E g.s. proj. wall src.g.s. proj. smeared src. noninteractingwall src.smeared src.
10 11 12 13 14 15 t [ a ] E e ff ( t ) [ M e V ] L = E noninteracting 1st proj. wall src.1st proj. smeared src. Figure 11.
The effective energy shift ∆ E J,n eff ( t ) from R J,n ΞΞ ( t ) , the correlation function projected to the n -th eigenstate at the sink on L = 48 a , for J = wall (black open up-triangles) and J = smear (purpleopen down-triangle). Red bands represent the energy shifts from the eigenvalues of H in the finite box,while black lines denote those of a free Hamiltonian H . (Left) The projection to the ground state ( n = 0 ),together with the effective energy shift in the direct method without projection for the wall source (redopen circles) and the smeared source (blue open squares). (Right) The projection to the 1st excited state( n = 1 ). Taken from [48].two baryon operators on the same space-time point as L (cid:88) x B ( x , t ) B ( x , t ) = (cid:88) p ˜ B ( p , t ) ˜ B ( − p , t ) , ˜ B ( p , t ) ≡ (cid:88) x B ( x , t ) e − i p · x , (38)where the above source operator couples to all momentum modes with almost equal weight. Since almostall previous studies on N N interactions in the direct method employed this type of the source operator, theirconclusions on the existences of both deuteron and di-neutron are not valid due to large contaminations. Once eigenmodes of H in the finite box are obtained, we can construct an improved sink operator for aparticular eigenstate, whose correlation function with the J source is given by R J,nBB ( t ) = (cid:88) r Ψ † n ( r ) R JBB ( r , t ) . (39)Fig. 11 shows the effective energy shift ∆ E J,n eff ( t ) calculated from R J,n ΞΞ ( t ) on L = 48 a with J = wall(black open up-triangles) and J = smear (purple open down-triangle), for the ground state (Left) and the1st excited state (Right), together with ∆ E or ∆ E , eigenvalues of H in the finite box (red bands) as wellas those of H (black lines). For the ground state in Fig. 11 (Left), the effective energy shift in the directmethod without projection are also plotted for the wall source (red open circles) and the smeared source(blue open squares). Note that Mainz2018 employed a source operator as ˜ B ( p = , t ) ˜ B ( − p = , t ) and they reported that “In the 27-plet (dineutron) sector, the finite volumeanalysis suggests that the existence of a bound state is unlikely.” . inya Aoki et al. Lattice QCD and baryon-baryon interactions
After the sink projection, the effective energy shifts agree well between wall and smeared sources around t/a ∼ , not only for the ground state but also for the 1st excited state. while the effective energy shiftsfor the ground state in the direct method without projection disagree between two sources. In particular, anagreement between two sources with sink projection for the 1st excited state is rather remarkable, sincevariational 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 , of H (redbands). Note that the effective energy shift for the 1st excited state, ∆ E wall , ( t ) , has larger errors since thecontribution 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 constructeigenfunctions, agreements in the effective energy shifts for the ground state as well as the 1st excited stateprovide a non-trivial consistency check between the HAL QCD method and the L¨uscher’s finite volumeformula (with proper projections to extract the finite volume spectra). We thus conclude from Fig. 11 notonly that the HAL QCD potential correctly describes the energy shifts of two baryons in the finite box forboth ground and excited states but also that these energy shifts can be extracted even for baryon-baryonsystems if and only if the sink/source operators are highly improved. We emphasize that improvement ofoperators has to be performed not by the tuning of the plateau-like structures but by a sophisticated methodsuch as the variational method [10] (or a method presented here). See [43] for a recent study toward sucha direction. In this section, we summarize results on nuclear potentials in the HAL QCD method.
We first show the results of nuclear forces in the parity-even channel ( S and S - D channels) atheavy quark masses obtained by the LO analysis for the derivative expansion of the potential. Since thestatistical fluctuations are smaller at heavier quark masses in lattice QCD, this study is a good starting pointto grasp the nature of lattice QCD nuclear forces. In addition, quark mass dependence of nuclear forces isof fundamental importance from a point of view of, e.g., anthropic principle, which cannot be studied byexperiments.In the case of S channel, we obtain the LO central force following Eq. (31). In the case of S - D channel, the LO potentials consist of the central and tensor forces, which can be obtained from the coupledchannel analysis between the S - and D -wave components as (cid:26) − H − ∂∂t + 14 m N ∂ ∂t (cid:27) R J ( r , t ) = [ V C ( r ) + V T ( r ) S + · · · ] R J ( r , t ) , (40)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 ) canbe obtained [23]. Since the A +1 representation couples to the angular momentum l = 0 , , , · · · , theseprojections are expected to serve as the relevant partial wave decomposition at low energies. We find In lattice QCD studies for the meson-meson scatterings [9], serious systematics from the excited state contaminations in the simple plateau fitting have beenwidely recognized and the variational method has been utilized to obtain the finite volume spectra rather reliably, which can be combined with the L¨uscher’sfinite volume formula to extract phase shifts.
Frontiers 19 inya Aoki et al.
Lattice QCD and baryon-baryon interactions -80 -40 0 40 80 1200.0 0.5 1.0 1.5 2.0 2.5 V (r) [ M e V ] r [fm] V (27) M PS = 1171 [MeV] M PS = 1015 [MeV] M PS = 837 [MeV] M PS = 672 [MeV] M PS = 469 [MeV] -80 -40 0 40 80 1200.0 0.5 1.0 1.5 2.0 2.5 V (r) [ M e V ] r [fm] V C (10 * ) M PS = 1171 [MeV] M PS = 1015 [MeV] M PS = 837 [MeV] M PS = 672 [MeV] M PS = 469 [MeV] -100 -80 -60 -40 -20 0 200.0 0.5 1.0 1.5 2.0 V (r) [ M e V ] r [fm] V T (10 * ) M PS = 1171 [MeV] M PS = 1015 [MeV] M PS = 837 [MeV] M PS = 672 [MeV] M PS = 469 [MeV] -10 0102030405060 0 50 100 150 200 250 300 350 δ [ D eg ] T lab [MeV] NN S M PS = 1171 [MeV] M PS = 1015 [MeV] M PS = 837 [MeV] M PS = 672 [MeV] M PS = 469 [MeV] Partial-Wave-Analysis -10 0102030405060 0 50 100 150 200 250 300 350 δ – , ε – [ D eg ] T lab [MeV] NN S - D M PS = 1171 [MeV] δ – S M PS = 1015 [MeV] δ – S M PS = 837 [MeV] δ – S M PS = 672 [MeV] δ – S M PS = 469 [MeV] δ – S δ – D ε – Figure 12. (Upper) Nuclear forces obtained from 3-flavor lattice QCD at M ps = S channel (27-plet in SU(3) f representation). (Middle) Central force in the S - D channel (10 ∗ -plet in SU(3) f representation). (Right) Tensor force in the S - D channel. (Lower) N N scattering phase shifts as a function of energy in the laboratory frame (colored solid lines), obtainedfrom 3-flavor lattice QCD at M ps = S channel. (Right) Results in the S - D channel (with Stapp’s convention).Figures are taken from [44].that the NBS correlation functions after P A +1 and (1 − P A +1 ) are dominated by S -wave and D -wavecomponents, respectively, indicating that the contaminations from l ≥ components are indeed small. Fora more advanced partial wave decomposition, see [49].We perform the calculations in quenched [26, 23], dynamical 2-flavor [50], dynamical 3-flavor [51, 52, 44] and dynamical (2+1)-flavor [53, 34, 45, 47] lattice QCD with variousquark masses. We here present the results obtained in 3-flavor lattice QCD at ( M ps , M oct ) =(1171 , , (1015 , , (837 , , (672 , , (469 , MeV [51, 52, 44]. In the case of ( M ps , M oct ) = (837 , , the value of quark masses m u = m d = m s nearly correspond to the physicalstrange quark mass. We generate gauge configurations with the RG-improved Iwasaki gauge action andnon-perturbatively O ( a ) -improved Wilson quark action on a L × T = 32 × lattice. The lattice spacingis a = 0 . 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 isemployed for each baryon at the sink.Shown in Fig. 12 (Upper) are the lattice QCD results for the potentials. We find that the results areinsensitive to the Euclidean time t , at which the NBS correlation function is evaluated, indicating thatthe derivative expansion is well converged. The obtained potentials are found to reproduce the qualitative M ps = m π = m K and M oct = m N = m Λ = m Σ = m Ξ in 3-flavor QCD. inya Aoki et al. Lattice QCD and baryon-baryon interactions 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 intriguingfeatures in the quark mass dependence of the potentials. At long distances, it is observed that the ranges ofthe tail structures in the central and tensor forces become longer at lighter quark masses. Such a behaviorcan be understood from the viewpoint of one-boson-exchange potential. At short distances, the repulsivecores in the central forces are found to be enhanced at lighter quark masses. This could be explained by theshort-range repulsion due to the one-gluon-exchange in the quark model, whose strength is proportionalto the inverse of the (constituent) quark mass. In fact, our systematic studies including hyperon forceswith the same lattice setup revealed that the nature of repulsive core is well described by the quark Pauliblocking effect together with the one-gluon-exchange effect [51, 44, 54].As noted before, the potentials themselves are not physical observables and quantitative lattice QCDpredictions shall be given in terms of scattering observables. Shown in Fig. 12 (Lower) are the scatteringphase shifts (and mixing angles) obtained from lattice nuclear forces. We find that
N N systems do notbound at these pseudoscalar masses as discussed in Sec. 3. Behaviors of phase shifts are qualitativelysimilar to the experimental ones, while the strength of the attraction is weaker due to the heavy quarkmasses 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-rangeattraction becomes stronger and the short-range repulsive core also becomes stronger. We also note thatlighter 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 ) (cid:39) (701, 1584), (570, 1412) and (411,1215) MeV [45]. Note that only up and down quark massesare varied with a strange quark mass being fixed to the physical value in this study. We employ the gaugeconfigurations generated by the PACS-CS Collaboration with the RG-improved Iwasaki gauge action andnon-perturbatively O ( a ) -improved Wilson quark action on a L × T = 32 × lattice. The lattice spacingis a (cid:39) . fm ( a − = 2 . GeV), which leads to the spatial extension L (cid:39) . fm.In Fig. 13, we show the lattice QCD results for the potentials in the S and S - D channels, togetherwith the corresponding phase shifts in the S channel. Qualitative features are similar to those in 3-flavorcase: (i) the central forces have repulsive cores at short distance and attractive wells at long and mediumdistances, both of which are enhanced at lighter quark masses (ii) the tensor force is strong with a negativesign, which increases at lighter quark masses. If we consider an interaction at higher order terms in the derivative expansion, there appear morestructures in the potentials. In particular, the extension from LO analysis to NLO analysis enables us todetermine the spin-orbit (LS) force. The LS force is known to play an important role in the LS-splittings ofnuclear spectra and the nuclear magic numbers. In addition, the LS force in the P - F channel attractsgreat interest in nuclear astrophysics, since it could lead to the P -wave superfluidity in the neutron starsand affect the cooling process of neutron stars.We here present the calculation in parity-odd channels ( P , P , P , P - F channels) at heavy quarkmasses and show the results of LS forces as well as central/tensor forces [50]. In order to construct thesource operator which couples to parity-odd states, we employ the two nucleon operators as J αβ ( f i ) ≡ N α ( f ( i ) ) N β ( f ( i ) ∗ ) for i = ± , ± , ± (41) Frontiers 21 inya Aoki et al.
Lattice QCD and baryon-baryon interactions -500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 0.5 1 1.5 2 2.5 V C (r ; S ) [ M e V ] r [fm]-40-20 0 20 40 60 80 0 0.5 1 1.5 2 2.5m π =411 MeVm π =570 MeVm π =700 MeV -10-5 0 5 10 15 0 50 100 150 200 250 300 δ ( S ) [ deg r ee ] E [MeV]m pion =700 MeVm pion =570 MeVm pion =411 MeV -500 0 500 1000 1500 2000 2500 3000 3500 0 0.5 1 1.5 2 2.5 V C (r ; S - D ) [ M e V ] r [fm]-40-20 0 20 40 60 80 0 0.5 1 1.5 2 2.5m π =411 MeVm π =570 MeVm π =700 MeV -140-120-100-80-60-40-20 0 20 0 0.5 1 1.5 2 2.5 V T (r ; S - D ) [ M e V ] r [fm] m π =411 MeVm π =570 MeVm π =700 MeV Figure 13.
Nuclear forces obtained from (2+1)-flavor lattice QCD at m π (cid:39)
411 (red), 570 (green), 701(blue) MeV: (Upper-Left) Central forces in the S channel (Lower) Central forces (left) and tensor forces(right) in the S - D channel. (Upper-Right) The scattering phase shifts in the S channel at m π (cid:39) N denotes a nucleon operator with a momentum, N α ( f ( i ) ) = (cid:88) (cid:126)x ,(cid:126)x ,(cid:126)x (cid:15) abc (cid:16) u Ta ( (cid:126)x ) Cγ d b ( (cid:126)x ) (cid:17) q c,α ( (cid:126)x ) f ( i ) ( (cid:126)x ) (42)with f ( ± j ) ( (cid:126)x ) ≡ exp ( ± πix j /L ) . A cubic group analysis shows that this source operator contains theorbital contribution T − ⊕ A +1 ⊕ E + , whose dominant components have l = 1 , , , respectively, andthus covers all the two-nucleon channels with J ≤ . Combined with the spin degrees of freedom, weconsider the T − representation in the spin singlet channel and the A − , T − , ( E − ⊕ T − ) representationsin the spin triplet channel. At low energies, these representations correspond to the P channel and the P , P and P - F channels, respectively, from which we extract the central force in the spin singletchannel ( 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 spintriplet channel.Calculations are performed in 2-flavor lattice QCD at quark masses corresponding to ( m π , m N ) (cid:39) (1133,2158) MeV [50]. We employ the gauge configurations generated by the CP-PACS Collaboration with theRG-improved Iwasaki gauge action and a mean field O ( a ) -improved Wilson quark action on a × lattice. The lattice spacing a = 0 . fm leads to the spatial extension L (cid:39) . fm.Shown in Fig. 14 (Upper-Left) are the lattice QCD results for the potential, V I =0 C,S =0 ( r ) , V I =1 C,S =1 ( r ) , V I =1 T ( r ) , V I =1 LS ( r ) . We find that (i) the central forces V I =0 C,S =0 ( r ) and V I =1C; S =1 ( r ) are repulsive, inya Aoki et al. Lattice QCD and baryon-baryon interactions -400-200 0 200 400 0 1 2 V (r) [ M e V ] r [fm] V I=0C;S=0 V I=1C;S=1 V I=1T V I=1LS -50 0 50 100 150 200 250 300 350 400 450 0 0.5 1 1.5 2 V (r) [ M e V ] r [fm] V(r; P )V(r; P )V(r, P )V(r; P ) -35-30-25-20-15-10-5 0 5 10 15 0 50 100 150 200 250 300 350 pha s e s h i ft ( deg r ee ) E lab = 2 k /m N [MeV] (a) -4-2 0 2 4 6 8 10 12 14 16 18 0 50 100 150 200 250 300 350 pha s e s h i ft ( deg r ee ) E lab = 2 k /m N [MeV] (b) P F mixing parameter P (EXP) F (EXP)mixing parameter (EXP) Figure 14. (Upper-Left) Central ( S = 0 and ), tensor and spin-orbit potentials in parity-odd channelsobtained by 2-flavor lattice QCD at m π (cid:39) MeV. (Upper-Right) The potentials for the P , P , P and P channels. (Lower-Left) Phase shifts in the P , P and P channels, together with theexperimental ones for comparisons. (Lower-Right) Phase shifts and mixing parameter (with Stapp’sconvention) in the P – F channel, together with the experimental ones. Figures are taken from [50].(ii) the tensor force V I =1T ( r ) is positive and weak compared to V I =1C; S =1 ( r ) and V I =1LS ( r ) , and (iii) theLS force V I =1LS ( r ) is negative and strong. These features are qualitatively in line well with those ofthe phenomenological potential. One can also see these properties in terms of the potential in eachchannel. In Fig. 14 (Upper-Right), we plot the potentials in the P , P , P and P channels,which are defined by V ( r ; P ) = V I =0C ,S =0 ( r ) , V ( r ; P ) = V I =1C ,S =1 ( r ) − V I =1T ( r ) − V I =1LS ( r ) , V ( r ; P ) = V I =1C ,S =1 ( r ) + 2 V I =1T ( r ) − V I =1LS ( r ) , V ( r ; P ) = V I =1C ,S =1 ( r ) − V I =1T ( r ) + V I =1LS ( r ) .To obtain the scattering observables, we fit the potentials and solve the Schr¨odinger equation in theinfinite volume. In Fig. 14 (Lower), we show the results for the scattering phase shifts. Compared with theexperimental phase shifts, we find that behaviors of phase shifts are generally well reproduced, while themagnitudes are smaller due to the heavier pion mass in lattice QCD calculations. In the P channel, weobserve that the attraction is missing compared with the experimental one, which however is also likely dueto the weak tensor force V T caused by the heavier pion mass. Among others, the most interesting featureis the attraction in the P channel as shown in Fig. 14 (Lower-Right), originated from the strong (andnegative) LS forces. As noted before, it is this interaction which is relevant to the paring correlation of theneutrons and possible P -wave superfluidity in the neutron stars. Frontiers 23 inya Aoki et al.
Lattice QCD and baryon-baryon interactions
We now turn to the study of three-nucleon forces. Determination of three-nucleon forces is one of themost challenging problems in nuclear physics: Three-nucleon forces are known to play important rolein nuclear spectra/structures such as the binding energies of (light) nuclei and properties of neutron-richnuclei. They are also essential ingredients to understand properties of nuclear matters such as the equationof state (EoS) at high density, which is relevant to the structures of neutron stars and nucleosynthesis atthe binary neutron star mergers. While there have been many studies to construct three-nucleon forces byphenomenological approaches [55, 56] or by chiral EFT approaches [57, 6, 7, 8], it is most desirable tocarry out the direct determination from QCD.To study three-nucleon forces in lattice QCD, we consider the NBS wave function for a n ( ≥ -particlesystem, | α (cid:105) , Ψ nα ([ x ]) e − W α t = (cid:104) | N ( x , t ) N ( x , t ) · · · N ( x n , t ) | α (cid:105) , [ x ] = x , x , · · · , x n (43)where W α is the center of mass energy of the system and we ignore the spins of nucleon for simplicity.In [58, 59, 60], we show that the asymptotic behavior of the NBS wave function with the non-relativisticapproximation can be written as Ψ n [ L ] , [ K ] ( R, Q ) ∝ (cid:88) [ N ] U [ L ][ N ] ( Q ) e iδ [ N ] ( Q ) sin (cid:0) QR − ∆ L + δ [ N ] ( Q ) (cid:1) ( QR ) D − U † [ N ][ K ] ( Q ) (44)where D = 3( n − is the dimension of a n -particle system, ∆ L = (2 L + D − π/ , Ψ n [ L ] , [ K ] ( R, Q ) isthe radial component of the NBS wave function in D -dimension with R and Q being the hyper radiusand 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 unitarymatrix, which parameterize the T -matrix as T [ L ][ K ] ( Q, Q ) = (cid:88) [ N ] U [ L ][ N ] ( Q ) T [ N ] ( Q ) U † [ N ][ K ] ( Q ) , (45) T [ N ] ( Q ) = − n / m N Q n − e iδ [ N ] ( Q ) sin δ [ N ] ( Q ) . (46)Therefore, as in the case of n = 2 system (See Sec. 2.2.1), the information of T -matrix is encoded inthe asymptotic behavior of the NBS wave function. Based on this property, we can define the energy-independent 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, R N ( (cid:126)r, (cid:126)ρ, t − t ) ≡ G N ( (cid:126)r, (cid:126)ρ, t − t ) / { G N ( t − t ) } (47) G N ( (cid:126)r, (cid:126)ρ, t − t ) ≡ L (cid:88) (cid:126)R (cid:104) | ( N ( (cid:126)x ) N ( (cid:126)x ) N ( (cid:126)x ))( t ) ( N (cid:48) N (cid:48) N (cid:48) )( t ) | (cid:105) (48)where (cid:126)R ≡ ( (cid:126)x + (cid:126)x + (cid:126)x ) / , (cid:126)r ≡ (cid:126)x − (cid:126)x , (cid:126)ρ ≡ (cid:126)x − ( (cid:126)x + (cid:126)x ) / are the Jacobi coordinates. In the time-dependent HAL QCD method at the LO analysis for the derivative expansion and with the non-relativisticapproximation, we can extract the three-nucleon forces V NF ( (cid:126)r, (cid:126)ρ ) through the following Schr¨odinger inya Aoki et al. Lattice QCD and baryon-baryon interactions -100-50 0 50 100 0 0.5 1 V N F ( r ) [ M e V ] r [fm]m π = 0.76 GeVm π = 0.93 GeVm π = 1.13 GeV -100-50 0 50 100 0 0.5 1 V N F ( r ) [ M e V ] r [fm] t=12t=11t=10-100 0 100 1 1.5 2 2.5 Figure 15.
Three-nucleon forces in the triton channel with the linear setup. (Left) Results from 2-flavorlattice QCD at m π = m π = 0 . GeV.equation, (cid:20) − µ r ∇ r − µ ρ ∇ ρ + (cid:88) i Lattice QCD and baryon-baryon interactions other hand, three nucleon forces are found to be suppressed at long distances. This is in accordance withthe suppression of two-pion-exchange due to the heavier pion masses.Shown in 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 [36] with the RG-improved Iwasakigauge action and non-perturbatively O ( a ) -improved Wilson quark action on a L × T = 64 × lattice(work in progress). The lattice spacing is a = 0 . fm and L = 5 . fm. Avoiding the very short-distanceregion where lattice discretization error could affect the results, we again find the short-range repulsivethree-nucleon forces at r (cid:39) m π = m π = 0 . GeV. It is important to pursue the study at lighter pion masses towards the physicalpion mass. Once nuclear potentials are obtained by lattice QCD, we can use them to study various phenomena innuclear physics and astrophysics. We here present the study of nuclear spectra/structures and Equation ofState (EoS) of dense matter relevant to neutron star physics. Potentials used in this subsection are of theleading order only, and therefore are all hermitian. We can make non-hermitian higher order potentials inthe HAL QCD method hermitian in the derivative expansion[63], which may be used for future applicationsin nuclear many body calculations.In [64], binding energies and structures of doubly magic nuclei, He, O, Ca, are studied by anab initio nuclear many-body calculation based on lattice nuclear forces. We employ the nuclear forcesobtained in 3-flavor lattice QCD at M ps = 469 MeV (See Fig. 12). We consider two-body nuclear forcesin S , S and D channels, while nuclear forces in other channels and spin-orbit forces as well asthree-nucleon forces are neglected. For simplicity, the Coulomb force between protons is not taken intoaccount, 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 solvedself-consistently in a nonperturbative manner. In a practical calculation, the self-energy is calculated byAlgebraic 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 [64] for details.)In Tab. 2, we summarize the results for the ground state energies, together with the results from BruecknerHartree-Fock (BHF) calculation [65] and exact stochastic variational calculation [66] using the same latticenuclear 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 He as a benchmark. Since the SCGF result deviates from the exact solution by less than 10% for He,and the SCGF approach is size extensive, we take a conservative estimate of 10% error for O and Ca,which are quoted in the second parentheses. The BHF results are sensibly more bound than the SCGFresults, and we interpret this as a limitation of BHF theory. For the results shown in Tab. 2, there existadditional errors associated with the statistical fluctuations in the input lattice nuclear forces, which areestimated to be ∼ 10% [65]. Note that statistical fluctuations are correlated among nuclei, so we expect ourobservations described below are rather robust against statistical errors.We find that at M ps =469 MeV in the SU(3) limit of QCD, both He and Ca have bound ground stateswhile the deuteron is unbound. O is likely to decay into four separate alpha particles, though it is alreadyclose to become bound. Moreover, we find that asymmetric isotopes, like O, are strongly unboundsystems. These results suggest that, when lowering the pion mass toward its physical value, closed shell inya Aoki et al. Lattice QCD and baryon-baryon interactions E A [MeV] He O CaSCGF -4.80(0.03) -17.9 (0.3) (1.8) -75.4 (6.7) (7.5)BHF -8.2 -34.7 -112.7Exact calc. -5.09 – –Experiment -28.3 -127.7 -342.0Separation into He clusters: -2.46 (0.3) (1.8) 24.5 (6.7) (7.5) Table 2. Ground state energies of He, O and 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 QCDmethod. Comparison is given with those obtained with BHF [65] and the exact calculation [66]. The lastline is the breakup energy for splitting the system in He clusters (of total energy A/ × He O Ca r pt − matter [fm]: SCGF 1.67 2.64 2.97BHF 2.09 2.35 2.78HF 1.62 2.39 2.78 r charge [fm]: SCGF 1.89 2.79 3.10Experiment 1.67 2.73 3.48 Table 3. Matter and charge radii of He, O and Ca at M P S =469 MeV computed by the SCGF method,which are compared with those by BHF [65], by Hartree-Fock (HF) and by experiments [67, 68]. Forcharge radii, we assumed the physical charge distributions of the nucleons. Taken from [64].isotopes are created at first around the traditional magic numbers and the region of M ps ∼ 500 MeV marksa 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 centralvalues. Although the total binding energies are 15-20% of the experimental value (Tab. 2), the computedcharge radii are about the same as the experiment. We also find that the calculated one-nucleon spectraldistributions are qualitatively close to those of real nuclei even for M ps =469 MeV considered here. Thisis due to the fact that the heavy nucleon mass ( m N =1161.1 MeV) used here reduces the motion of thenucleons inside the nuclei and counterbalances the effect of weak attraction of the lattice nuclear forces atthis 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 S , S and D channels obtained in 3-flavor lattice QCD. Tostudy 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 theBrueckner-Hartree-Fock (BHF) theory [70], which is known to be quantitative enough to give the essentialunderlying physics for infinite matter.In Fig. 16 (Upper), we show the results of the ground state energy per nucleon ( E/A ) as a function ofthe Fermi momentum k F for the symmetric nuclear matter and the pure neutron matter. Shown togetherare the so-called APR EoS [69], which are obtained by the variational chain summation method fromphenomenological nuclear forces with (APR(Full)) and without (APR(AV18)) three-nucleon forces. InFig. 16 (Upper-Left), we find that the symmetric nuclear matter becomes a self-bound system with asaturation point ( k F , E/A ) (cid:39) (1 . − , − . at the lightest quark mass ( M ps = 469 Frontiers 27 inya Aoki et al. Lattice QCD and baryon-baryon interactions -20 0 20 40 60 800.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 E / A [ M e V ] k F [fm −1 ] Weizsackermass formula M PS = 1171 [MeV] M PS = 1015 [MeV] M PS = 837 [MeV] M PS = 672 [MeV] M PS = 469 [MeV] APR(AV18)APR(Full) E / A [ M e V ] k F [fm −1 ] M PS = 1171 [MeV] M PS = 1015 [MeV] M PS = 837 [MeV] M PS = 672 [MeV] M PS = 469 [MeV] APR(AV18)APR(Full) M [ M Sun ] R [km] M PS = 1171 [MeV] M PS = 1015 [MeV] M PS = 837 [MeV] M PS = 672 [MeV] M PS = 469 [MeV] Figure 16. (Upper) Ground state energy per nucleon ( E/A ) as a function of the Fermi momentum k F bythe BHF theory with nuclear forces from 3-flavor lattice QCD at M ps = 469-1171 MeV, together with thatfrom APR [69] with and without phenomenological three-nucleon forces. (Left) Results for the symmetricnuclear matter. filled square indicates the empirical saturation point. (Right) Results for the pure neutronmatter. (Lower) Mass-radius relation of the neutron star based on EoS obtained by the BHF theory withnuclear forces from 3-flavor lattice QCD at M ps = 469-1171 MeV. Figures are taken from [70].MeV). This is the first time that the saturation in the symmetric nuclear matter is obtained through first-principles lattice QCD simulations. The saturation point, however, deviates from the empirical pointprimarily due to heavy pion (pseudo-scalar meson) mass in lattice simulation and the lack of three-nucleonforces in BHF calculation.We also find a nontrivial M ps dependence of the EoS: the saturation disappears at intermediate pion masses( M ps = 672 , MeV) and possibly appears again at the heavy pion mass region ( M ps = 1015 , MeV). This implies that the saturation originates from a subtle balance between short-range repulsionand the intermediate attraction of the nuclear force, which have different m q dependence [44]. A similarnontrivial M ps dependence originated from the balance between repulsion and attraction is also observedfor 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 wedecrease the pion mass, EoS is found to become stiffer. To further study the impact on phenomena in nuclearastrophysics, we calculate the mass ( M ) versus the radius ( R ) relation of neutron stars at each pion mass. inya Aoki et al. Lattice QCD and baryon-baryon interactions Figure 17. An illustration of the complementary role of lattice QCD and experiments in the determinationof baryon forces.Here, we solve the Tolman-Oppenheimer-Volkoff (TOV) equation by using the EoS of the neutron-starmatter with neutron, proton, electron and muon under the charge neutrality and beta equilibrium, where weuse the standard parabolic approximation for asymmetric nuclear matters.Shown in 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 themaximum mass of the neutron star ( M max ) in this calculation is much smaller than the recent observations, M max (cid:39) M (cid:12) , the deviation is most likely due to the heavy pion masses and lack of interactions asthree-nucleon forces. A naive extrapolation of M max and the corresponding radius to M ps = 137 MeVwould give M max ∼ . M (cid:12) and R ∼ km, which are encouraging for more quantitative studies infuture. Another hottest topic in the context of neutron star physics is the effect of hyperon on the EoSat high density (so-called “hyperon puzzle”). Lattice QCD can play an unique role to study this effectby determining the hyperon forces which suffer from large uncertainties in experiments to date. For theon-going study in this direction, see [71]. While the results of nuclear forces at heavy pion masses are very intriguing and useful to extract thephysical picture of nuclear forces, the quantitative results require the study at the physical pion mass. Notethat the pion mass dependence of nuclear forces is quite non-trivial as discussed in Secs. 4.1 and 4.3, so thedirect calculation near the physical point is desirable.To this end, we have recently performed the first calculation of nuclear forces near the physical up, downand 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 − in parity-even channels ( S - and D -waves). As mentioned before, the statistical fluctuations in latticeQCD are smaller (larger) for larger (smaller) quark masses, and thus the results have better precision insectors involving more number of strange quarks (larger strangeness | S | ). On the other hand, experimentsin such larger | S | sectors are more difficult due to the short life time of hyperons. Therefore, lattice QCDstudies and experiments are complementary with each other in the determination of baryon forces (SeeFig. 17).(2+1)-flavor gauge configurations are generated on a L × T = 96 × lattice with the RG-improvedIwasaki gauge action and non-perturbatively O ( a ) -improved Wilson quark action and APE stout smearing.The lattice spacing is a (cid:39) . fm ( a − (cid:39) . GeV), so that spatial extent, L = 8 . fm, is sufficientlylarge 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 ) (cid:39) (146 , , MeV. NBS correlation functionsfor two-baryon systems are calculated for 55 channels in total and we extract the central and tensorforces in parity-even channel at the LO analysis for the derivative expansion (work in progress, and see Frontiers 29 inya Aoki et al. Lattice QCD and baryon-baryon interactions Figure 18. Nuclear forces from (2+1)-flavor lattice QCD near the physical point, m π = 146 MeV. Thecentral force in the S channel (Left). The central force (Middle) and the tensor force (Right) in the S - D channel.also [72]). In order to make this first calculation a reality, “trinity” of state-of-the-art developments wascrucial: (a) time-dependent HAL QCD method (theory), (b) unified contraction algorithm (software) and(c) K-computer, HOKUSAI and HA-PACS supercomputers (hardware).Shown in Fig. 18 are the results for the central force in the S channel (Left), and the central force(Middle) and tensor force (Right) in the S - D channel. As noted above, nuclear forces are the mostchallenging interactions in lattice QCD calculation, and one can see that the results suffer from largestatistical fluctuations. Nevertheless, the obtained results exhibit several interesting properties.First of all, the repulsive core at short-range is clearly observed in central forces. In order to clarify thephysical picture for the repulsive core, it is useful to compare them with hyperon forces obtained in thesame lattice setup. We find that the strength of repulsive core (or attractive core) highly depends on theflavor SU(3) (SU(3) f ) classification, in a consistent way with the quark Pauli blocking effect. In addition, ifwe compare interactions which belong to the same SU(3) f classification, such as N N ( S ) and ΞΞ( S ) both of which belong to 27-plet, we find that the strength differs in a way which can be understood fromthe viewpoint of one-gluon-exchange (e.g., repulsive core in N N ( S ) is stronger than that in ΞΞ( S ) ).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 quarkmasses. See also [73] 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 isattractive, resembling the phenomenological potential as one-pion-exchange potential (OPEP). The tensorforce has relatively smaller statistical errors than the central forces, showing that the tensor force becomesstronger (with a negative sign) and has a longer tail, as compared with the tensor forces at heavier pionmasses (Sec. 4.1). This property can be understood by the picture of OPEP. These results are encouragingand serve as the first step to establish a direct connection between QCD and nuclear physics. At the sametime, statistical errors remain to be large and there also exist systematic errors associated with inelastic statecontaminations. The studies to resolve these issues are in progress, and the second generation calculation isplanned on the forthcoming Exascale computer, “Fugaku’ (See https://postk-web.r-ccs.riken.jp/). Before closing this review, we present our latest results on dibaryon searches in lattice QCD near thephysical pion mass [72]. A dibaryon, a bound-state (or a resonance) with a baryon number B = 2 in QCD, inya Aoki et al. Lattice QCD and baryon-baryon interactions can be classified in the SU(3) f as ⊗ = ⊕ s ⊕ ⊕ ⊕ ⊕ a (50)for the octet-octet system, where the deuteron, the only stable dibaryon observed in nature so far, appearsin the representation while H dibaryon has been predicted in the representation [74] and activelyinvestigated in lattice QCD [51, 52, 44, 75, 43]. For the decuplet-octet system, the classification leads to ⊗ = ⊕ ⊕ ⊕ (51)and N Ω ( N ∆ ) dibaryon has been predicted in the ( ) representation [76, 77, 78], and ⊗ = ⊕ ⊕ ⊕ (52)for the decuplet-decuplet system, where ΩΩ dibaryon has been predicted in the representation [79]while ∆∆ has been predicted in the [78, 80] and the corresponding d ∗ (2380) was indeed observed [81].Note that among decuplet baryons, only Ω is stable against strong decays. We first consider the ΩΩ system in the representation of SU(3) f in the S channel [82].Fig. 19 (Upper-Left) shows ΩΩ potentials at t/a = 16 , , , which has qualitative features similar tothe central potentials for N N but whose repulsion is weaker and attraction is shorter-ranged. This potentialpredicts 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 searchedexperimentally by two-particle correlations in relativistic heavy-ion collisions [84]. N Ω dibaryon We next consider the N Ω system with S = − in the representation in the S channel [83]. Near thephysical point, N Ω ( S ) 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 – , which is attractive at all distanceswithout repulsive core, so that one bound state appears in this channel. In Fig. 19 (Lower-Right), thebinding energy (vertical) and the the root-mean-square distance (horizontal) are plotted for n Ω − with noCoulomb interaction (red) and p Ω − with Coulomb attraction (blue). These binding energies are muchsmaller than B = 18 . . +12 . − . ) MeV at heavy pion mass m π = 875 MeV [85]. Such a N Ω statecan be searched through two-particle correlations in relativistic nucleus-nucleus collisions [84] and anexperimental indication was also reported [86]. Let us consider the scattering length a and the effective range r eff for ΩΩ ( S ) and N Ω ( S ). In Fig. 20,the ratio r eff /a as a function of r eff are plotted for ΩΩ ( S ) and N Ω ( S ) obtained in lattice QCD nearthe physical pion mass, together with the experimental values for N N ( S ) (deuteron) and N N ( S )(di-neutron). Small values of | r eff /a | in all cases indicate that these systems are located close to the unitarylimit. Frontiers 31 inya Aoki et al. Lattice QCD and baryon-baryon interactions ������� � � � �� � ��� � � � ��� � � � ��� � � � ��� � � � ��� � � � � � ����� ������������������ ������������������ � � � ��� � � � � � � � � � � � � � ���� � � � � � � � ��� � � � � � � � � � ���������������� � �������� � ���� �������������� r [fm] B [ M e V ] p n Figure 19. (Upper) The ΩΩ system in the S channel in flavor QCD at m π (cid:39) MeV and a (cid:39) . fm on a (8.1 fm) box. (Left) The ΩΩ potential V ( r ) at t/a = 16 , , . (Right) The bindingenergy of the ΩΩ system and the root-mean-square distance between two Ω ’s are shown by blue soliddiamond (red solid triangle), calculated from the ΩΩ potential V ( r ) at t/a = 17 without (with) theCoulomb repulsion. Taken from [82]. (Lower) The N Ω system in the S channel with the same latticesetup for ΩΩ . (Left) The N Ω potential V C ( r ) at t/a = 11 , , , . (Right) The binding energy and theroot-mean-square distance for the n Ω − (red open circle) and p Ω − (blue open square). Taken from [83]. In this paper, we have reviewed the recent progress in lattice QCD study of baryon-baryon interactions bythe HAL QCD method. We first presented the detailed account on how to define the potentials in quantumfield theories such as QCD. The key observation is that the Nambu-Bethe-Salpeter (NBS) wave functionscontain the information of scattering phase shifts below inelastic threshold in their asymptotic behaviorsoutside the range of the interactions. The potentials at the interaction region can then be defined throughthe NBS wave functions so as to be faithful to the phase shifts by construction, where the non-localityof the potential is defined by the derivative expansion. In addition, by constructing the potentials inenergy-independent way, the potentials can be extracted from two-baryon correlation functions without therequirement of the ground state saturation.We then made a detailed comparison between the HAL QCD method and the conventional method, inwhich phase shifts are obtained from the finite volume energies through the L ¨uscher’s formula. We pointedout that, while the validity of the latter method relies on the ground state saturation of the correlation inya Aoki et al. Lattice QCD and baryon-baryon interactions r eff [fm] r e ff / a ( S ) N ( S ) NN ( S ) NN ( S ) Figure 20. The ratio of the effective range and the scattering length r eff /a as a function of r eff for ΩΩ( S ) (blue open diamond) and N Ω( S ) (red open circle) obtained in lattice QCD, as well as for N N ( S ) (purple open up-triangle) and N N ( S ) (green open down-triangle) in experiments. Taken from[83]. The sign convention for the scattering length is opposite to Eq. (5) in this figure.function, its practical procedure for multi-baryon systems (“direct method”) so far has utilized only theplateau-like structures of the effective energies at Euclidean times much earlier than the inverse of thelowest excitation energy. We showed theoretical and numerical evidences that such a procedure generallyleads to unreliable results due to the contaminations from the elastic excited states: For instance, the resultswere found to be dependent on the operators and unphysical behaviors were exposed by the normalitycheck. This invalidates the claim of the literature in the direct method that N N bound states exist for pionmasses heavier than 300 MeV.On the other hand, HAL QCD method is free from such a serious problem since the signal of potentialscan be extracted from not only the ground state but also elastic excited states. While there instead exists thetruncation error of the derivative expansion of the potential, the calculation of the higher order term in thederivative 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 statecontaminations in the direct method were explicitly quantified. It turns out that the plateau-like structuresof effective energies at early time slices are indeed pseudo-plateaux contaminated by elastic excited statesand that the plateau for the ground state is realized only at a much larger time corresponding to the inverseof the lowest excitation energy gap. We also showed that, by employing an optimized operator utilizingthe finite volume eigenmodes, the effective energies from the correlation functions give consistent resultswith those from the HAL QCD potential. Thus the long-standing issue on the consistency between theconventional method based on the L¨uscher’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 nuclearforces from the HAL QCD method at various lattice QCD setups. At heavy pion masses, where goodsignal-to-noise ratio can be achieved in lattice QCD, we observed that the obtained N N potentials in theparity-even channel ( S , S - D ) reproduce the qualitative features of the phenomenological potentials,namely, attractive wells at long and medium distances, accompanied with repulsive cores at short distancein the central potentials and the strong tensor force. The net interactions were found to be attractive, whichhowever are not strong enough to form a bound N N state, probably due to the heavy pion masses. We Frontiers 33 inya Aoki et al. Lattice QCD and baryon-baryon interactions observed that the tail structures are enhanced at lighter pion masses, which can be understood from theviewpoint of one-pion exchange contributions. We also found the repulsive cores are enhanced at lighterpion masses. Combined with our systematic studies including hyperon forces, the nature of repulsive coreswas found to be well described by the quark Pauli blocking effect together with the one-gluon-exchangecontribution.The HAL QCD method can be extended to determine more complicated nuclear forces, such as spin-orbitforces and three-nucleon forces. In this paper, we considered two-nucleon systems in the parity-oddchannels ( P , P , P , P - F channels) and calculated spin-orbit forces as well as central and tensorforces. We found that qualitative features of experimental results are generally well reproduced, whilethe 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 P channel. Three-nucleon forces were studied inthe triton channel, ( I, J P ) = (1 / , / + ) , thank to the unified contraction algorithm (UCA), which canenormously speed up calculations of multi-baryon correlation functions. It was found that there exists arepulsive three-nucleon forces at short distances. These observations are of interest in the context of notonly the structures of nuclei but also those of neutron stars, e.g., P -wave superfluidity and the maximummass of neutron stars.We carried out the applications to nuclei, nuclear equation of state (EoS) and structure of neutron starsbased on lattice nuclear forces at heavy quark masses. We performed ab initio self-consistent Green’sfunction (SCGF) calculations for closed shell nuclei with nuclear forces at M ps =469 MeV in the SU(3) limitof QCD. We found that He, Ca nuclei are bound and O is close to become bound, while asymmetricisotopes are strongly unbound. The results suggest that, when lowering the pion mass toward its physicalvalue, islands of stability appear at first around the traditional doubly magic numbers. The nuclear EoS wasalso studied by the BHF theory with nuclear forces in the flavor SU(3) limit. We found that the saturationproperty appears in the symmetric nuclear matter at M ps =469 MeV. A mass-radius relation of the neutronstar was also studied based on the EoS obtained from lattice nuclear forces and we observed a tendencythat 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 theparity-even channel. The computation is quite challenging particularly for nuclear forces due to badsignal-to-noise ratio near the physical point. Nevertheless, we observed prominent characteristics of nuclearforces, such as repulsive cores at short distances as well as attractive interactions at mid and long distancesin central forces, and a strong (and negative) tensor force. We also presented the results for the hyperonforces obtained near the physical point. We found that both ΩΩ( S ) and N Ω( S ) systems have strongattractions, and (quasi) bound dibaryons are formed near the unitary limit. These systems could be searchedexperimentally through two-particle correlations in relativistic nucleus-nucleus collisions.Present results shown in this paper already indicate a clear pathway which connects nuclear physics withits underlying theory of the strong interaction, QCD. While there remain many challenges to accomplishresearches in this direction, there is no doubt that successive theoretical developments together withnext-generation supercomputers will further deepen the connection between the two. The outcome is alsoexpected to play a crucial role to understand the nuclear astrophysical phenomena such as supernovaexplosions and mergers of binary neutron stars, as well as the nucleosynthesis associated with theseexplosive phenomena. inya Aoki et al. Lattice QCD and baryon-baryon interactions ACKNOWLEDGEMENTS This work is supported in part by the Grant-in-Aid of the Japanese Ministry of Education, Sciencesand Technology, Sports and Culture (MEXT) for Scientific Research (Nos. JP16H03978, JP18H05236,JP18H05407, JP19K03879), by HPCI programs (hp120281, hp130023, hp140209, hp150223, hp150262,hp150085, hp160211, hp160093, hp170230, hp170170, hp180179, hp180117, hp190160, hp190103),by a priority issue (Elucidation of the fundamental laws and evolution of the universe) to be tackled byusing Post “K” Computer, and by Joint Institute for Computational Fundamental Science (JICFuS). Theauthors thank members of the HAL QCD Collaboration for providing lattice QCD results and for fruitfulcollaborations based on which this paper is prepared. Figure 12 is reprinted from [44] with permissionfrom Elsevier. Figure 14 is taken from [50], and Figures 19 and 20 are taken from [83], under the term ofthe https://creativecommons.org/licenses/by/3.0/ . REFERENCES [ 1] Stoks VGJ, Klomp RAM, Terheggen CPF, de Swart JJ. Construction of high quality N N potentialmodels. Phys. Rev. C49 (1994) 2950–2962. doi:10.1103/PhysRevC.49.2950. [ 2] Wiringa RB, Stoks VGJ, Schiavilla R. An Accurate nucleon-nucleon potential with chargeindependence breaking. Phys. Rev. C51 (1995) 38–51. doi:10.1103/PhysRevC.51.38. [ 3] Machleidt R. The High precision, charge dependent Bonn nucleon-nucleon potential (CD-Bonn).Phys. Rev. C63 (2001) 024001. doi:10.1103/PhysRevC.63.024001. [ 4] Weinberg S. Nuclear forces from chiral Lagrangians. Phys. Lett. B251 (1990) 288–292. doi:10.1016/0370-2693(90)90938-3. [ 5] Weinberg S. Effective chiral Lagrangians for nucleon - pion interactions and nuclear forces. Nucl.Phys. B363 (1991) 3–18. doi:10.1016/0550-3213(91)90231-L. [ 6] Epelbaum E, Hammer HW, Meissner UG. Modern Theory of Nuclear Forces. Rev. Mod. Phys. (2009) 1773–1825. doi:10.1103/RevModPhys.81.1773. [ 7] Machleidt R, Entem DR. Chiral effective field theory and nuclear forces. Phys. Rept. (2011)1–75. doi:10.1016/j.physrep.2011.02.001. [ 8] Hammer HW, Koenig S, van Kolck U. Nuclear effective field theory: status and perspectives (2019). [ 9] Briceno RA, Dudek JJ, Young RD. Scattering processes and resonances from lattice QCD. Rev. Mod.Phys. (2018) 025001. doi:10.1103/RevModPhys.90.025001. [ 10] Luscher M, Wolff U. How to Calculate the Elastic Scattering Matrix in Two-dimensional QuantumField Theories by Numerical Simulation. Nucl. Phys. B339 (1990) 222–252. doi:10.1016/0550-3213(90)90540-T. [ 11] Rummukainen K, Gottlieb SA. Resonance scattering phase shifts on a nonrest frame lattice. Nucl.Phys. B (1995) 397–436. doi:10.1016/0550-3213(95)00313-H. [ 12] Briceno RA, Davoudi Z, Luu TC. Two-Nucleon Systems in a Finite Volume: (I) QuantizationConditions. Phys. Rev. D (2013) 034502. doi:10.1103/PhysRevD.88.034502. [ 13] Briceno RA. Two-particle multichannel systems in a finite volume with arbitrary spin. Phys. Rev. D (2014) 074507. doi:10.1103/PhysRevD.89.074507. [ 14] Hansen MT, Sharpe SR. Relativistic, model-independent, three-particle quantization condition. Phys.Rev. D (2014) 116003. doi:10.1103/PhysRevD.90.116003. [ 15] Hansen MT, Sharpe SR. Lattice QCD and Three-particle Decays of Resonances. Ann. Rev. Nucl.Part. Sci. (2019) 65–107. doi:10.1146/annurev-nucl-101918-023723. Frontiers 35 inya Aoki et al. Lattice QCD and baryon-baryon interactions [ 16] Durr S, et al. Ab-Initio Determination of Light Hadron Masses. Science (2008) 1224–1227.doi:10.1126/science.1163233. [ 17] Borsanyi S, et al. Ab initio calculation of the neutron-proton mass difference. Science (2015)1452–1455. doi:10.1126/science.1257050. [ 18] Luscher M. Two particle states on a torus and their relation to the scattering matrix. Nucl. Phys. B354 (1991) 531–578. doi:10.1016/0550-3213(91)90366-6. [ 19] Iritani T, Aoki S, Doi T, Hatsuda T, Ikeda Y, Inoue T, et al. Are two nucleons bound in lattice QCDfor heavy quark masses? Consistency check with L ¨uscher’s finite volume formula. Phys. Rev. D96 (2017) 034521. doi:10.1103/PhysRevD.96.034521. [ 20] Sitenko AG. Lectures in scattering theory (New York: Pergamon Press Oxford) (1971). [ 21] Kawai D, Aoki S, Doi T, Ikeda Y, Inoue T, Iritani T, et al. I = 2 ππ scattering phase shift from theHAL QCD method with the LapH smearing. PTEP (2018) 043B04. doi:10.1093/ptep/pty032. [ 22] Ishizuka N. Derivation of Luscher’s finite size formula for N pi and NN system. PoS LAT2009 (2009)119. doi:10.22323/1.091.0119. [ 23] Aoki S, Hatsuda T, Ishii N. Theoretical Foundation of the Nuclear Force in QCD and its applicationsto Central and Tensor Forces in Quenched Lattice QCD Simulations. Prog. Theor. Phys. (2010)89–128. doi:10.1143/PTP.123.89. [ 24] Lin CJD, Martinelli G, Sachrajda CT, Testa M. K → ππ decays in a finite volume. Nucl. Phys. B619 (2001) 467–498. doi:10.1016/S0550-3213(01)00495-3. [ 25] Aoki S, et al. I=2 pion scattering length from two-pion wave functions. Phys. Rev. D71 (2005)094504. doi:10.1103/PhysRevD.71.094504. [ 26] Ishii N, Aoki S, Hatsuda T. The Nuclear Force from Lattice QCD. Phys. Rev. Lett. (2007) 022001.doi:10.1103/PhysRevLett.99.022001. [ 27] Aoki S, Doi T, Hatsuda T, Ikeda Y, Inoue T, Ishii N, et al. Lattice QCD approach to Nuclear Physics.PTEP (2012) 01A105. doi:10.1093/ptep/pts010. [ 28] Sugiura T, Ishii N, Oka M. Derivative Expansion of Wave Function Equivalent Potentials. Phys. Rev. D95 (2017) 074514. doi:10.1103/PhysRevD.95.074514. [ 29] Aoki S, Doi T, Hatsuda T, Ishii N. Comment on “Relation between scattering amplitude and Bethe-Salpeter wave function in quantum field theory”. Phys. Rev. D98 (2018) 038501. doi:10.1103/PhysRevD.98.038501. [ 30] Yamazaki T, Kuramashi Y. Relation between scattering amplitude and Bethe-Salpeter wave functionin quantum field theory. Phys. Rev. D96 (2017) 114511. doi:10.1103/PhysRevD.96.114511. [ 31] Ernst DJ, Shakin CM, Thaler RM. Separable Representations of Two-Body Interactions. Phys. Rev. C8 (1973) 46–52. doi:10.1103/PhysRevC.8.46. [ 32] Ernst DJ, Shakin CM, Thaler RM. Separable representations of T matrices valid in the vicinity ofoff-shell points. Phys. Rev. C9 (1974) 1780–1783. doi:10.1103/PhysRevC.9.1780. [ 33] Murano K, Ishii N, Aoki S, Hatsuda T. Nucleon-Nucleon Potential and its Non-locality in LatticeQCD. Prog. Theor. Phys. (2011) 1225–1240. doi:10.1143/PTP.125.1225. [ 34] Ishii N, Aoki S, Doi T, Hatsuda T, Ikeda Y, Inoue T, et al. Hadron-hadron interactions from imaginary-time Nambu-Bethe-Salpeter wave function on the lattice. Phys. Lett. B712 (2012) 437–441. doi:10.1016/j.physletb.2012.04.076. [ 35] Yamazaki T, Kuramashi Y, Ukawa A. Two-Nucleon Bound States in Quenched Lattice QCD. Phys.Rev. D84 (2011) 054506. doi:10.1103/PhysRevD.84.054506. [ 36] Yamazaki T, Ishikawa Ki, Kuramashi Y, Ukawa A. Helium nuclei, deuteron and dineutron in 2+1flavor lattice QCD. Phys. Rev. D86 (2012) 074514. doi:10.1103/PhysRevD.86.074514. inya Aoki et al. Lattice QCD and baryon-baryon interactions [ 37] Orginos K, Parreno A, Savage MJ, Beane SR, Chang E, Detmold W. Two nucleon systems at m π ∼ 450 MeV from lattice QCD. Phys. Rev. D92 (2015) 114512. doi:10.1103/PhysRevD.92.114512. [ 38] Beane SR, Chang E, Detmold W, Lin HW, Luu TC, Orginos K, et al. The Deuteron and ExoticTwo-Body Bound States from Lattice QCD. Phys. Rev. D85 (2012) 054511. doi:10.1103/PhysRevD.85.054511. [ 39] Yamazaki T, Ishikawa Ki, Kuramashi Y, Ukawa A. Study of quark mass dependence of binding energyfor light nuclei in 2+1 flavor lattice QCD. Phys. Rev. D92 (2015) 014501. doi:10.1103/PhysRevD.92.014501. [ 40] Beane SR, Chang E, Cohen SD, Detmold W, Lin HW, Luu TC, et al. Light Nuclei and Hypernucleifrom Quantum Chromodynamics in the Limit of SU(3) Flavor Symmetry. Phys. Rev. D87 (2013)034506. doi:10.1103/PhysRevD.87.034506. [ 41] Wagman ML, Winter F, Chang E, Davoudi Z, Detmold W, Orginos K, et al. Baryon-Baryon Interactionsand Spin-Flavor Symmetry from Lattice Quantum Chromodynamics. Phys. Rev. D96 (2017) 114510.doi:10.1103/PhysRevD.96.114510. [ 42] Berkowitz E, Kurth T, Nicholson A, Joo B, Rinaldi E, Strother M, et al. Two-Nucleon Higher Partial-Wave Scattering from Lattice QCD. Phys. Lett. B765 (2017) 285–292. doi:10.1016/j.physletb.2016.12.024. [ 43] Francis A, Green JR, Junnarkar PM, Miao C, Rae TD, Wittig H. Lattice QCD study of the H dibaryon using hexaquark and two-baryon interpolators. Phys. Rev. D99 (2019) 074505. doi:10.1103/PhysRevD.99.074505. [ 44] Inoue T, Aoki S, Doi T, Hatsuda T, Ikeda Y, Ishii N, et al. Two-Baryon Potentials and H-Dibaryonfrom 3-flavor Lattice QCD Simulations. Nucl. Phys. A881 (2012) 28–43. doi:10.1016/j.nuclphysa.2012.02.008. [ 45] Ishii N. Baryon-baryon Interactions from Lattice QCD. PoS CD12 (2013) 025. doi:10.22323/1.172.0025. [ 46] Iritani T, et al. Mirage in Temporal Correlation functions for Baryon-Baryon Interactions in LatticeQCD. JHEP (2016) 101. doi:10.1007/JHEP10(2016)101. [ 47] Iritani T, Aoki S, Doi T, Gongyo S, Hatsuda T, Ikeda Y, et al. Systematics of the HAL QCD Potentialat Low Energies in Lattice QCD. Phys. Rev. D99 (2019) 014514. doi:10.1103/PhysRevD.99.014514. [ 48] Iritani T, Aoki S, Doi T, Hatsuda T, Ikeda Y, Inoue T, et al. Consistency between L¨uscher’s finitevolume method and HAL QCD method for two-baryon systems in lattice QCD. JHEP (2019) 007.doi:10.1007/JHEP03(2019)007. [ 49] Miyamoto T, Akahoshi Y, Aoki S, Aoyama T, Doi T, Gongyo S, et al. Partial wave decomposition onthe lattice and its applications to the HAL QCD method (2019). [ 50] Murano K, Ishii N, Aoki S, Doi T, Hatsuda T, Ikeda Y, et al. Spin–orbit force from lattice QCD. Phys.Lett. B735 (2014) 19–24. doi:10.1016/j.physletb.2014.05.061. [ 51] Inoue T, Ishii N, Aoki S, Doi T, Hatsuda T, Ikeda Y, et al. Baryon-Baryon Interactions in the FlavorSU(3) Limit from Full QCD Simulations on the Lattice. Prog. Theor. Phys. (2010) 591–603.doi:10.1143/PTP.124.591. [ 52] Inoue T, Ishii N, Aoki S, Doi T, Hatsuda T, Ikeda Y, et al. Bound H-dibaryon in Flavor SU(3) Limit ofLattice QCD. Phys. Rev. Lett. (2011) 162002. doi:10.1103/PhysRevLett.106.162002. [ 53] Ishii N, Aoki S, Hatsuda T. Nuclear forces from quenched and 2+1 flavor lattice QCD using thePACS-CS gauge configurations. PoS LATTICE2008 (2008) 155. doi:10.22323/1.066.0155. [ 54] Oka M, Shimizu K, Yazaki K. Quark cluster model of baryon baryon interaction. Prog. Theor. Phys.Suppl. (2000) 1–20. doi:10.1143/PTPS.137.1. Frontiers 37 inya Aoki et al. Lattice QCD and baryon-baryon interactions [ 55] Coon SA, Han HK. Reworking the Tucson-Melbourne three nucleon potential. Few Body Syst. (2001) 131–141. doi:10.1007/s006010170022. [ 56] Pieper SC, Pandharipande VR, Wiringa RB, Carlson J. Realistic models of pion exchange threenucleon interactions. Phys. Rev. C64 (2001) 014001. doi:10.1103/PhysRevC.64.014001. [ 57] Weinberg S. Three body interactions among nucleons and pions. Phys. Lett. B295 (1992) 114–121.doi:10.1016/0370-2693(92)90099-P. [ 58] Aoki S, Charron B, Doi T, Hatsuda T, Inoue T, Ishii N. Construction of energy-independent potentialsabove inelastic thresholds in quantum field theories. Phys. Rev. D87 (2013) 034512. doi:10.1103/PhysRevD.87.034512. [ 59] Aoki S, Ishii N, Doi T, Ikeda Y, Inoue T. Asymptotic behavior of Nambu-Bethe-Salpeter wavefunctions for multiparticles in quantum field theories. Phys. Rev. D88 (2013) 014036. doi:10.1103/PhysRevD.88.014036. [ 60] Gongyo S, Aoki S. Asymptotic behavior of Nambu-Bethe-Salpeter wave functions for scalar systemswith a bound state. PTEP (2018) 093B03. doi:10.1093/ptep/pty097. [ 61] Doi T, Endres MG. Unified contraction algorithm for multi-baryon correlators on the lattice. Comput.Phys. Commun. (2013) 117. doi:10.1016/j.cpc.2012.09.004. [ 62] Doi T, Aoki S, Hatsuda T, Ikeda Y, Inoue T, Ishii N, et al. Exploring Three-Nucleon Forces in LatticeQCD. Prog. Theor. Phys. (2012) 723–738. doi:10.1143/PTP.127.723. [ 63] Aoki S, Iritani T, Yazaki K. Hermitizing the HAL QCD potential in the derivative expansion. PTEP (2020) 023. doi:10.1093/ptep/ptz166. [ 64] McIlroy C, Barbieri C, Inoue T, Doi T, Hatsuda T. Doubly magic nuclei from lattice QCD forces at M P S = 469 MeV/c . Phys. Rev. C97 (2018) 021303. doi:10.1103/PhysRevC.97.021303. [ 65] Inoue T, Aoki S, Charron B, Doi T, Hatsuda T, Ikeda Y, et al. Medium-heavy nuclei from nucleon-nucleon interactions in lattice QCD. Phys. Rev. C91 (2015) 011001. doi:10.1103/PhysRevC.91.011001. [ 66] Nemura H. Recent developments on LQCD studies of nuclear force. Int. J. Mod. Phys. E23 (2014)1461006. doi:10.1142/S0218301314610060. [ 67] De Vries H, De Jager CW, De Vries C. Nuclear charge and magnetization density distributionparameters from elastic electron scattering. Atom. Data Nucl. Data Tabl. (1987) 495–536. doi:10.1016/0092-640X(87)90013-1. [ 68] Angeli I, Marinova KP. Table of experimental nuclear ground state charge radii: An update. Atom.Data Nucl. Data Tabl. (2013) 69–95. doi:10.1016/j.adt.2011.12.006. [ 69] Akmal A, Pandharipande VR, Ravenhall DG. The Equation of state of nucleon matter and neutronstar structure. Phys. Rev. C58 (1998) 1804–1828. doi:10.1103/PhysRevC.58.1804. [ 70] Inoue T, Aoki S, Doi T, Hatsuda T, Ikeda Y, Ishii N, et al. Equation of State for Nucleonic Matterand its Quark Mass Dependence from the Nuclear Force in Lattice QCD. Phys. Rev. Lett. (2013)112503. doi:10.1103/PhysRevLett.111.112503. [ 71] Inoue T. Strange Nuclear Physics from QCD on Lattice. AIP Conf. Proc. (2019) 020002.doi:10.1063/1.5118370. [ 72] Doi T, et al. Baryon interactions from lattice QCD with physical quark masses – Nuclear forces and ΞΞ forces –. EPJ Web Conf. (2018) 05009. doi:10.1051/epjconf/201817505009. [ 73] Park A, Lee SH, Inoue T, Hatsuda T. Baryon-baryon interactions at short distances – constituent quarkmodel meets lattice QCD (2019). [ 74] Jaffe RL. Perhaps a Stable Dihyperon. Phys. Rev. Lett. (1977) 195–198. doi:10.1103/PhysRevLett.38.617,10.1103/PhysRevLett.38.195. [Erratum: Phys. Rev. Lett.38,617(1977)]. inya Aoki et al. Lattice QCD and baryon-baryon interactions [ 75] Beane SR, et al. Evidence for a Bound H-dibaryon from Lattice QCD. Phys. Rev. Lett. (2011)162001. doi:10.1103/PhysRevLett.106.162001. [ 76] Goldman JT, Maltman K, Stephenson GJ Jr, Schmidt KE, Wang F. STRANGENESS -3 DIBARYONS.Phys. Rev. Lett. (1987) 627. doi:10.1103/PhysRevLett.59.627. [ 77] Oka M. Flavor Octet Dibaryons in the Quark Model. Phys. Rev. D38 (1988) 298. doi:10.1103/PhysRevD.38.298. [ 78] Dyson F, Xuong NH. Y=2 States in Su(6) Theory. Phys. Rev. Lett. (1964) 815–817. doi:10.1103/PhysRevLett.13.815. [ 79] Zhang ZY, Yu YW, Shen PN, Dai LR, Faessler A, Straub U. Hyperon nucleon interactions in a chiralSU(3) quark model. Nucl. Phys. A625 (1997) 59–70. doi:10.1016/S0375-9474(97)00033-X. [ 80] Kamae T, Fujita T. Possible Existence of a Deeply Bound Delta-Delta System. Phys. Rev. Lett. (1977) 471–475. doi:10.1103/PhysRevLett.38.471. [ 81] Adlarson P, et al. ABC Effect in Basic Double-Pionic Fusion — Observation of a new resonance?Phys. Rev. Lett. (2011) 242302. doi:10.1103/PhysRevLett.106.242302. [ 82] Gongyo S, et al. Most Strange Dibaryon from Lattice QCD. Phys. Rev. Lett. (2018) 212001.doi:10.1103/PhysRevLett.120.212001. [ 83] Iritani T, et al. N Ω dibaryon from lattice QCD near the physical point. Phys. Lett. B792 (2019)284–289. doi:10.1016/j.physletb.2019.03.050. [ 84] Morita K, Gongyo S, Hatsuda T, Hyodo T, Kamiya Y, Ohnishi A. Probing ΩΩ and p Ω dibaryonswith femtoscopic correlations in relativistic heavy-ion collisions. Phys. Rev. C101 (2020) 015201.doi:10.1103/PhysRevC.101.015201. [ 85] Etminan F, Nemura H, Aoki S, Doi T, Hatsuda T, Ikeda Y, et al. Spin-2 N Ω dibaryon from LatticeQCD. Nucl. Phys. A928 (2014) 89–98. doi:10.1016/j.nuclphysa.2014.05.014. [ 86] Adam J, et al. The Proton- Ω correlation function in Au+Au collisions at √ s NN =200 GeV. Phys. Lett. B790 (2019) 490–497. doi:10.1016/j.physletb.2019.01.055.(2019) 490–497. doi:10.1016/j.physletb.2019.01.055.