Dibaryon with highest charm number near unitarity from lattice QCD
Yan Lyu, Hui Tong, Takuya Sugiura, Sinya Aoki, Takumi Doi, Tetsuo Hatsuda, Jie Meng, Takaya Miyamoto
MMost charming dibaryon near unitarity
Yan Lyu,
1, 2, ∗ Hui Tong,
1, 3, † Takuya Sugiura, Sinya Aoki,
4, 2
Takumi Doi,
2, 3
Tetsuo Hatsuda, Jie Meng,
1, 5 and Takaya Miyamoto State Key Laboratory of Nuclear Physics and Technology,School of Physics, Peking University, Beijing 100871, China Quantum Hadron Physics Laboratory, RIKEN Nishina Center, Wako 351-0198, Japan Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS), RIKEN, Wako 351-0198, Japan Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan
A pair of triply charmed baryons, Ω ccc Ω ccc , is studied as an ideal dibaryon system by (2+1)-flavorlattice QCD with nearly physical light-quark masses and the relativistic heavy quark action withthe physical charm quark mass. The spatial baryon-baryon correlation is related to their scatteringparameters on the basis of the HAL QCD method. The Ω ccc Ω ccc in the S channel taking intoaccount the Coulomb repulsion with the charge form factor of Ω ccc leads to the scattering length a C0 (cid:39) −
19 fm and the effective range r Ceff (cid:39) .
45 fm. The ratio r Ceff /a C0 (cid:39) − . − . ccc Ω ccc is located inthe unitary regime. Introduction. − Quantum chromodynamics (QCD) is afundamental theory of strong interaction and governs notonly the interaction among quarks and gluons but alsothe interaction between color-neutral hadrons. In partic-ular, the nucleon-nucleon (
N N ) interaction, which showsa characteristic mid-range attraction and a short-rangerepulsion, as well as the baryon-baryon ( BB ) interac-tions are important for describing the nuclear structureand dense matter relevant to nuclear physics and astro-physics [1–5].Although the deuteron is the only stable bound statecomposed of two nucleons, there are possible bound orresonant dibaryons with and without strange quarks [6–8]. Among others, p Ω( uudsss ) [9] and ΩΩ( ssssss ) [10],which were predicted by lattice QCD (LQCD) simula-tions near the physical point [11], stimulate experimen-tal searches in high energy hadron-hadron and heavy-ioncollisions [8, 12–14].As originally pointed out by Bjorken [15], the triplycharmed baryon (the charm number C = 3) Ω ccc is sta-ble against the strong interaction and provides an idealground to study the perturbative and non-perturbativeaspects of QCD in the baryonic sector. Although it hasnot been observed yet experimentally there have beennumerous LQCD studies on its mass and electromagneticform factor (see [18] and references therein). Accordingly,it is timely to study the Ω ccc Ω ccc as the simplest pos-sible system to study heavy-baryon interactions. Its re-cent phenomenological study using the constituent quarkmodel can be found in Ref. [19].The purpose of this Letter is to study a system withthe charm number C = 6 system, Ω ccc Ω ccc in the S channel, for the first time from first principle LQCD ap- Recently, excited states of C = 1 baryon Ω c [16] and a C = 2baryon Ξ ++ cc [17] were discovered at CERN LHC. proach. The reason why we consider the S -wave andtotal spin s = 0 system is that the Pauli exclusion be-tween charm quarks at short distance does not operate inthis channel, so that the maximum attraction is expectedin comparison to other channels. It is of critical impor-tance to examine the scattering parameters such as thescattering length and the effective range to unravel theproperties of such heavy dibaryons near threshold. TheHAL QCD method [11, 22, 23], which treats the spa-tial correlation between two baryons on the lattice, pro-vides a powerful tool for such analysis: Indeed, we showbelow that Ω ++ ccc Ω ++ ccc ( S ) with both strong interactionand Coulomb repulsion is located near unitarity [24, 25]just above the threshold with a large negative scatteringlength. HAL QCD Method. − The crucial steps in the HALQCD method [11, 22, 23] are to obtain the equal-time Nambu-Bethe-Salpeter (NBS) wave function ψ ( r )whose asymptotic behavior at a large distance repro-duces the phase shifts, along with the corresponding two-baryon irreducible kernel U ( r , r (cid:48) ). Since the same kernel U ( r , r (cid:48) ) governs all the elastic scattering states, separat-ing the ground state and the excited states on the lattice,which is exponentially difficult for baryon-baryon inter-actions [26, 27], is not required to calculate the physicalobservables [23]. The normalized four-point function (the R -correlator) related to the NBS wave function is definedas R ( r , t >
0) = (cid:104) | Ω ccc ( r , t )Ω ccc ( , t ) J (0) | (cid:105) /e − m Ω ccc t = (cid:88) n A n ψ n ( r ) e − (∆ W n ) t + O ( e − (∆ E ∗ ) t ) , (1) In the charm number C = 3 sector, there exist a few recentstudies on heavy dibaryons in LQCD [20] and in the constituentquark model [21]. a r X i v : . [ h e p - l a t ] F e b where ∆ W n = 2 (cid:113) m ccc + k n − m Ω ccc with the baryonmass m Ω ccc and the relative momentum k n . O ( e − (∆ E ∗ ) t )denotes the contributions from the inelastic scatteringstates with ∆ E ∗ being the inelastic threshold, which areexponentially suppressed when t (cid:29) (∆ E ∗ ) − ∼ Λ − with Λ QCD ∼
300 MeV. J (0) is a source operator whichcreates two-baryon states with the charm number C = 6at Euclidean time t = 0 and A n = (cid:104) n |J (0) | (cid:105) with | n (cid:105) representing the QCD eigenstates in a finite volume with∆ W n < ∆ E ∗ . In this study, we take a local interpolatingoperator, Ω ccc ( x ) ≡ (cid:15) lmn [ c Tl ( x ) C γ k c m ( x )] c n,α ( x ), where l , m , and n stand for color indices, γ k being the Diracmatrix, α being the spinor index, and C ≡ γ γ being thecharge conjugation matrix.When contributions from the inelastic scattering statesare negligible ( t (cid:29) (∆ E ∗ ) − ), the R -correlator satis-fies [23] (cid:18) m Ω ccc ∂ ∂t − ∂∂t − H (cid:19) R ( r , t ) = (cid:90) d r (cid:48) U ( r , r (cid:48) ) R ( r (cid:48) , t ) , (2)where H = −∇ /m Ω ccc . By using the derivative ex-pansion at low energies, U ( r , r (cid:48) ) = V ( r ) δ ( r − r (cid:48) ) + (cid:80) n =1 V n ( r ) ∇ n δ ( r − r (cid:48) ), the central potential V ( r ) in theleading order (LO) is given as V ( r ) = R − ( r , t ) (cid:18) m Ω ccc ∂ ∂t − ∂∂t − H (cid:19) R ( r , t ) . (3)The spatial and temporal derivatives of R ( r , t ) on thelattice are calculated in central difference scheme byusing the nearest neighbor points. To extract thetotal spin s = 0, the following interpolating opera-tors for the Ω ccc Ω ccc system is adopted, [Ω ccc Ω ccc ] = (Ω / ccc Ω − / ccc − Ω / ccc Ω − / ccc + Ω − / ccc Ω / ccc − Ω − / ccc Ω / ccc ).Here the spin and its z component of the interpolat-ing operator Ω s z ccc are 3 / s z = ± / , ± /
2, re-spectively, and Ω s z ccc is constructed by spin projectionas shown in Ref. [28]. To obtain the orbital angularmomentum L = 0 on the lattice, the projection to A representation of the cubic group SO (3 , Z ) is employed; P A R ( r , t ) = (cid:80) R i ∈ SO (3 , Z ) R ( R i [ r ] , t ). Lattice setup. − (2+1)-flavor gauge configurations aregenerated on the L = 96 lattice with the Iwasaki gaugeaction at β = 1 .
82 and nonperturbatively O ( a )-improvedWilson quark action combined with stout smearing atnearly physical quark masses ( m π (cid:39)
146 MeV and m K (cid:39)
525 MeV) [29]. The lattice cutoff is a − (cid:39) .
333 GeV( a (cid:39) . La (cid:39) . S charmonium, i.e. a weighted average of the spin-singletstate η c and the spin-triplet state J/ Ψ.For the source operator J (0), we use the wall typewith the Coulomb gauge fixing. We employ the peri-odic (Drichlet) boundary condition for spatial (tempo-ral) direction. We use 112 gauge configurations whichare picked up one per ten trajectories. In order to re-duce statistical fluctuations, forward and backward prop-agations are averaged, and four times measurements areperformed by shifting source position along the temporaldirection for each configuration. Then, the total mea-surements amount to 896 for each set. The statisticalerrors are estimated by the jackknife method with a binsize of 14 configurations. A comparison with a bin sizeof 7 configurations shows that the bin size dependenceis small. The quark propagators are calculated by theBridge++ code [32], and the unified contraction algo-rithm is utilized to obtain the correlation functions [33]. TABLE I. Spin-averaged 1 S charmonium mass (( m η c +3 m J/ Ψ ) /
4) and the Ω ccc mass ( m Ω ccc ) calculated in set 1 andset 2 with the statistical errors. The third row shows the in-terpolated values obtained from set 1 and set 2. Experimentalvalue of ( m η c + 3 m J/ Ψ ) / m η c + 3 m J/ Ψ ) / m Ω ccc [MeV]set 1 3096 . .
3) 4837 . . . .
3) 4770 . . . .
3) 4795 . . . .
1) -
Masses for spin-averaged 1 S charmonium (( m η c +3 m J/ Ψ ) /
4) and Ω ccc baryon ( m Ω ccc ) calculated in set1 and set 2 by utilizing the single exponential fittingfrom the interval t/a = 25 −
35 are listed in Ta-ble. I, together with the values from linear interpolation(0 . × set 1 + 0 . × set 2) as well as the exper-imental value. Our result for m Ω ccc is consistent with4789(6)(21) MeV obtained by the (2+1)-flavor PACS-CSconfigurations [34]. We have checked that our results forhadron masses are unchanged within errors by the fittinginterval t/a = 30 − Numerical results. − The Ω ccc Ω ccc potential V ( r ) in the S channel from the interpolation between set 1 andset 2 is shown in Fig. 1 for t/a = 25 ,
26 and 27. Sincethe potentials from set 1 and set 2 are found to be con-sistent within statistical errors, the uncertainty in theinterpolation is negligible. Our choice t/a = 26 corre-sponds to t (cid:39) . − ∼ . t/a =25, 26 and 27 are consistent with each other within sta-tistical errors. This indicates that systematic errors dueto inelastic states and higher order terms of the deriva- r [fm] V ( r ) [ M e V ] t / a = 27 t / a = 26 t / a = 25 FIG. 1. (Color online). The Ω ccc Ω ccc potential V ( r ) in the S channel as a function of separation r at Euclidean time t/a = 25 (red square), 26 (blue diamond) and 27 (green cir-cle). tive expansion do not largely exceed the size of statisticalerrors [23] as we show below.We find that the potential V ( r ) is repulsive at short-range and attractive at mid-range, which has the samequalitative behaviors with the N N potential [35] andthe ΩΩ potential [10]. The magnitude of the poten-tial in the repulsive region r < .
25 fm (correspond-ing to dV ( r ) /dr <
0) for Ω ccc Ω ccc is an order of mag-nitude smaller than that of ΩΩ obtained by the samemethod [10]. This may be qualitatively explained by thephenomenological quark model [36] as the color-magneticinteraction between constituent quarks is proportional tothe square of reciprocal constituent quark mass. Qual-itatively, V cc cm /V ss cm = ( m ∗ s /m ∗ c ) ∼ (500 / ∼ . V ff (cid:48) cm is the color-magnetic interaction between thequarks with flavor f and f (cid:48) with m ∗ f being the constituentquark mass. On the other hand, the attraction in the re-gion r > .
25 fm (corresponding to dV ( r ) /dr >
0) mayoriginate from the exchange of charmed mesons or ratherbe attributed to the direct exchange of charm quarksand/or multiple gluons. As can be seen in Fig. 1, therange of the potential is much smaller than the size ofthe lattice volume, indicating that the finite volume ar-tifact is negligible.In order to convert the potential to physical ob-servables such as the scattering phase shifts and bind-ing energy, we perform the uncorrelated fit for V ( r )in Fig. 1 in the range r ≤ . V fit ( r ) = (cid:80) i =1 , , α i exp( − β i r ). Fit-ting parameters with t/a = 26 for example are( α , α , α ) = (239 . . , − . . , − . . β , β , β ) = (48 . . , . . , . . − with an accuracy of χ / d.o.f. ∼ . ccc Ω ccc scattering phase E CM [MeV] [ d e g ] t / a = 27 t / a = 26 t / a = 25 FIG. 2. (Color online). The Ω ccc Ω ccc scattering phase shifts δ in the S channel obtained from the potential V ( r ) at t/a = 25 ,
26, and 27 as a function of the center of mass kineticenergy E CM . shifts δ in the S channel calculated by solving theSchr¨odinger equation with the potential V ( r ) at t/a =25 ,
26, and 27. The relativistic kinetic energy is definedas E CM = 2 (cid:113) k + m ccc − m Ω ccc with a momentum k in the center of mass frame. The error bands reflectthe statistical uncertainty of V ( r ). In all three cases, thephase shifts start from 180 ◦ at E CM = 0, which indicatesthe existence of a bound state in Ω ccc Ω ccc system withoutCoulomb repulsion.The low-energy scattering parameters are extracted byusing the effective range expansion up to the next-to-leading order (NLO), k cot δ = − a + r eff k + O ( k ),where a and r eff are the scattering length and the effec-tive range, respectively. The results are a = 1 . . +0 . − . ) fm ,r eff = 0 . . +0 . − . ) fm . (4)The central values and the statistical errors in the firstparentheses are obtained at t/a = 26, while the system-atic errors in the last parentheses are estimated from thevalues at t/a = 25, 26 and 27, which originates from theinelastic states and the higher order terms of the deriva-tive expansion.The binding energy B and the root-mean-square dis-tance (cid:112) (cid:104) r (cid:105) of the bound Ω ccc Ω ccc state are obtainedfrom the potential V ( r ) as B = 5 . . +0 . − . ) MeV , (cid:112) (cid:104) r (cid:105) = 1 . . +0 . − . ) fm . (5)These results are consistent with the general formula forloosely bound states [24, 37] with scattering parameters a and r eff : B = 1 / ( m Ω ccc r )(1 − (cid:112) − (2 r eff /a )) (cid:39) . (cid:112) (cid:104) r (cid:105) = a / √ (cid:39) . e / phys. e / a C [ f m ] r d = 0.410 fmw/ sta. errorw/ sta. & sys. error r d = 0 fm FIG. 3. (Color online). The inverse of the scattering length1 /a C0 as a function of α e /α phys .e . The red solid line is thecentral values for r d = 0 .
410 fm. The statistical errors areshown by the inner band (red), while the outer band (gray)corresponds to the statistical and systematic errors added inquadrature. The blue dashed line corresponds to the centralvalues for r d = 0 fm. Since the binding energy and the size of the boundstate from the strong interaction are not large, we needto take into account the Coulomb repulsion V Coulomb ( r )between Ω ++ ccc s with finite spatial size. For this pur-pose, we consider the dipole form factor for Ω ++ ccc ac-cording to the LQCD study on the charge distributionof heavy baryons [18]: In the coordinate space, it cor-responds to an exponential charge distribution ρ ( r ) =12 √ / ( πr d ) e − √ r/r d , where the charge radius r d = (cid:112) |(cid:104) r (cid:105) charge | of Ω ++ ccc is taken to be r d = 0 . V Coulomb ( r ) = α e (cid:90) (cid:90) d r d r ρ ( r ) ρ ( | (cid:126)r − (cid:126)r | ) | (cid:126)r − (cid:126)r | = 4 α e r F ( x ) , (6)where x = 2 √ r/r d and F ( x ) = 1 − e − x (1 + x + x + x ). The effective range expansion with Coulomb re-pulsion is written as k (cid:2) C η cot δ C0 ( k ) + 2 ηh ( η ) (cid:3) = − a C0 + 12 r Ceff k + O ( k ) , (7)where δ C0 ( k ) is the phase shift in the presence of Coulombrepulsion, C η = πηe πη − , η = 2 α e m Ω ccc /k , h ( η ) =Re[Ψ( iη )] − ln( η ), and Ψ is the digamma function [38].To see the effect of the Coulomb repulsion, we vary α e from zero to the physical value α phys .e = 1 / .
036 be-low. Note that the systematic errors originated from theuncertainty in r d are found to be much smaller than thestatistical errors and are neglected.In Fig. 3, we show the inverse of scattering length1 /a C0 under the change of α e /α phys .e from 0 to 1. Due r eff [fm] r e ff / a w/ Coulombw/o Coulombw/ Coulombw/o Coulomb nnpp ccc ccc ( S ) ( S ) np ( S - D ) NN ( S ) FIG. 4. (Color online). The dimensionless ratio of the effec-tive range r eff and the scattering length a as a function of r eff .The red up(down)-pointing triangle and the blue right(left)-pointing triangle correspond to Ω ccc Ω ccc system and ΩΩ sys-tem in the S channel with(without) the Coulomb repulsionrespectively. The black circle represents NN system in the S - D channel. The green square ( nn ) and diamond ( pp )correspond to NN system in the S channel. The error barsfor Ω ccc Ω ccc are the quadrature of the statistical and system-atic errors in Eqs. (4) and (8). to the large cancellation between the attractive stronginteraction and the Coulomb repulsion, the result at α e /α phys .e = 1 is located very close to unitarity with alarge scattering length a C0 = − +7 − ) fm ,r Ceff = 0 . . +0 . − . ) fm . (8)The ratio r Ceff /a C0 = − . . +0 . − . ) is consider-ably smaller in magnitude than that of the dineutron( − . r eff /a as afunction of r eff for Ω ++ ccc Ω ++ ccc ( S ) and Ω − Ω − ( S ) with(without) Coulomb repulsion together with the experi-mental values for N N ( S - D ) [39] and N N ( S ) [40,41]. Note that we consider the Coulomb repulsionin Ω − Ω − ( S ) with the charge radius r d = 0 .
57 fmfor Ω − [18]. Among all those dibaryon systems,Ω ++ ccc Ω ++ ccc ( S ) is the closest to unitarity.Finally, we briefly discuss other possible systematicerrors in this work: (i) The finite cutoff effect is O ( α s a Λ QCD , ( a Λ QCD ) ) thanks to the RHQ action forthe charm quark and the non-perturbative O ( a ) improve-ment for light ( u, d, s ) quarks, and thus amounts to be In Ref.[10], Ω − was assumed to be point-like in charge distribu-tion, which overestimates the repulsion and increases the scat-tering length by 1 fm. O (1)%. (ii) In the vacuum polarization, light quarkmasses are slightly heavier than the physical ones andcharm quark loop is neglected. The former effect is ex-pected to be small since light quarks are rather irrelevantfor Ω ccc Ω ccc system. In fact, the range of the Ω ccc Ω ccc potential is found to be shorter than 1 fm. The lattereffect is suppressed due to the heavy charm quark mass,and is typically O (1)% [42]. These estimates for (i) and(ii) are also in line with the observation that our valueof m Ω ccc is consistent with that in the literature or hasdeviation of ∼
1% at most, where we refer to LQCDstudies by (2+1)-flavor at the physical point with finite a [34], (2+1)-flavor with chiral and continuum extrapola-tion [43] and (2+1+1)-flavor with chiral and continuumextrapolation [44, 45]. In the future, these systematicerrors will be evaluated explicitly. Summary and discussions. − In this Letter, we pre-sented a first investigation on the scattering properties ofthe Ω ccc Ω ccc on the basis of the (2+1)-flavor lattice QCDsimulations with physical charm mass and nearly physi-cal light quark masses. The potential for Ω ccc Ω ccc ( S )obtained by the time-dependent HAL QCD method with-out the Coulomb interaction shows a weak repulsion atshort distance surrounded by a relatively strong attrac-tive well, which leads to a most charming ( C = 6)dibaryon with the binding energy B (cid:39) . (cid:112) (cid:104) r (cid:105) (cid:39) . ++ ccc s with their charge formfactor obtained from LQCD, the Ω ++ ccc Ω ++ ccc ( S ) systemturns into the unitary region with r Ceff /a C0 (cid:39) − . − bbb Ω − bbb ( S ) for revealingthe quark mass dependence of the scattering parameters.Finally, our results may further stimulate the future ex-perimental activities to measure pair-momentum corre-lations of heavy baryons in high energy pp , pA and AA collisions [8, 14]. Acknowledgments. − We thank the members of HALQCD Collaboration for technical supports and stimulat-ing discussions. We thank Yusuke Namekawa for pro-viding the RHQ parameters. We thank members ofPACS Collaboration for the gauge configuration gener-ation conducted on the K computer at RIKEN. Thelattice QCD measurements have been performed onHOKUSAI supercomputers at RIKEN. This work waspartially supported by HPCI System Research Project(hp120281, hp130023, hp140209, hp150223, hp150262,hp160211, hp170230, hp170170, hp180117, hp190103).We thank ILDG/JLDG [46], which serves as an essen-tial infrastructure in this study. We thank the authorsof cuLGT code [47] for the gauge fixing. We thank Tat-sumi Aoyama, Haozhao Liang, Shuangquan Zhang andPengwei Zhao for helpful discussions. Y.L., H.T. andJ.M. were partially supported by the National Key R&D Program of China (Contracts No. 2017YFE0116700and No. 2018YFA0404400) and the National Natu-ral Science Foundation of China (NSFC) under GrantsNo. 11935003, No. 11975031, No. 11875075, and No.12070131001. This work was partially supported by JSPSGrant (No. JP18H05236, JP16H03978, JP19K03879,JP18H05407) and MOST-RIKEN Joint Project “Ab ini-tio investigation in nuclear physics”, “Priority Issue onPost-K computer” (Elucidation of the Fundamental Lawsand Evolution of the Universe), “Program for PromotingResearches on the Supercomputer Fugaku” (Simulationfor basic science: from fundamental laws of particles tocreation of nuclei) and Joint Institute for ComputationalFundamental Science (JICFuS). ∗ Y. L. and H. T. contributed equally to this Letter andshould be considered as co-first [email protected] † Corresponding [email protected][1] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Rev.Mod. Phys. , 1773 (2009).[2] J. Meng, Relativistic Density Functional for NuclearStructure (World Scientific, Singapore, 2016).[3] S. Shen, H. Liang, W. H. Long, J. Meng, and P. Ring,Progress in Particle and Nuclear Physics , 103713(2019).[4] C. Drischler, W. Haxton, K. McElvain, E. Mereghetti,A. Nicholson, P. Vranas, and A. Walker-Loud (2019)arXiv:1910.07961 [nucl-th].[5] H. Tong, P. Zhao, and J. Meng, Phys. Rev. C ,035802 (2020).[6] H. Clement, Progress in Particle and Nuclear Physics ,195 (2017).[7] A. Gal, Acta Phys. Polon. B , 471 (2016),arXiv:1511.06605 [nucl-th].[8] S. Cho, T. Hyodo, D. Jido, C. M. Ko, S. H. Lee, S. Maeda,K. Miyahara, K. Morita, M. Nielsen, A. Ohnishi, T. Sek-ihara, T. Song, S. Yasui, and K. Yazaki, Progress inParticle and Nuclear Physics , 279 (2017).[9] T. Iritani, S. Aoki, T. Doi, F. Etminan, S. Gongyo,T. Hatsuda, Y. Ikeda, T. Inoue, N. Ishii, T. Miyamoto,and K. Sasaki, Physics Letters B , 284 (2019).[10] S. Gongyo, K. Sasaki, S. Aoki, T. Doi, T. Hatsuda,Y. Ikeda, T. Inoue, T. Iritani, N. Ishii, T. Miyamoto,and H. Nemura (HAL QCD Collaboration), Phys. Rev.Lett. , 212001 (2018).[11] S. Aoki and T. Doi, Frontiers in Physics , 307 (2020).[12] K. Morita, S. Gongyo, T. Hatsuda, T. Hyodo, Y. Kamiya,and A. Ohnishi, Phys. Rev. C , 015201 (2020),arXiv:1908.05414 [nucl-th].[13] S. Acharya et al. (ALICE Collaboration), Nature ,232 (2020).[14] L. Fabbietti, V. M. Sarti, and O. V. Doce, (2020),arXiv:2012.09806 [nucl-ex].[15] J. Bjorken, AIP Conf. Proc. , 390 (1985).[16] R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett. , 182001 (2017). [17] R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett. , 112001 (2017).[18] K. U. Can, G. Erkol, M. Oka, and T. T. Takahashi,Phys. Rev. D , 114515 (2015).[19] H. Huang, J. Ping, X. Zhu, and F. Wang,arXiv:2011.00513 [hep-ph].[20] P. Junnarkar and N. Mathur, Phys. Rev. Lett. ,162003 (2019).[21] J.-M. Richard, A. Valcarce, and J. Vijande, Phys. Rev.Lett. , 212001 (2020).[22] N. Ishii, S. Aoki, and T. Hatsuda, Phys. Rev. Lett. ,022001 (2007).[23] N. Ishii, S. Aoki, T. Doi, T. Hatsuda, Y. Ikeda, T. In-oue, K. Murano, H. Nemura, and K. Sasaki (HAL QCDCollaboration), Physics Letters B , 437 (2012).[24] E. Braaten and H.-W. Hammer, Physics Reports ,259 (2006).[25] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev.Mod. Phys. , 1225 (2010).[26] G. P. Lepage, in From Actions to Answers: Proceedingsof the TASI 1989 , edited by T. Degrand and D. Toussaint(World Scientific, Singapore, 1990).[27] T. Iritani, S. Aoki, T. Doi, T. Hatsuda, Y. Ikeda, T. In-oue, N. Ishii, H. Nemura, and K. Sasaki (HAL QCD),JHEP , 007 (2019), arXiv:1812.08539 [hep-lat].[28] M. Yamada, K. Sasaki, S. Aoki, T. Doi, T. Hat-suda, Y. Ikeda, T. Inoue, N. Ishii, K. Murano, andH. Nemura (HAL QCD), PTEP , 071B01 (2015),arXiv:1503.03189 [hep-lat].[29] K.-I. Ishikawa, N. Ishizuka, Y. Kuramashi, Y. Nakamura,Y. Namekawa, Y. Taniguchi, N. Ukita, T. Yamazaki,and T. Yoshie (PACS), Proc. Sci.
LATTICE2015 , 075(2016), arXiv:1511.09222 [hep-lat].[30] S. Aoki, Y. Kuramashi, and S.-i. Tominaga, Prog. Theor. Phys. , 383 (2003), arXiv:hep-lat/0107009.[31] Y. Namekawa (PACS),
Proc. Sci.
LATTICE2016 , 125(2017).[32] http://bridge.kek.jp/Lattice-code/index_e.html .[33] T. Doi and M. G. Endres, Computer Physics Communi-cations , 117 (2013).[34] Y. Namekawa, S. Aoki, K.-I. Ishikawa, N. Ishizuka,K. Kanaya, Y. Kuramashi, M. Okawa, Y. Taniguchi,A. Ukawa, N. Ukita, and T. Yoshi´e (PACS-CS Collabo-ration), Phys. Rev. D , 094512 (2013).[35] T. Doi et al. , Proc. Sci.
LATTICE2016 , 110 (2017),arXiv:1702.01600 [hep-lat].[36] M. Oka, K. Shimizu, and K. Yazaki, Nuclear Physics A , 700 (1987).[37] P. Naidon and S. Endo, Reports on Progress in Physics , 056001 (2017).[38] P. G. Burke, R-Matrix Theory of Atomic Collisions (Springer-Verlag, Berlin, 2011).[39] R. W. Hackenburg, Phys. Rev. C , 044002 (2006).[40] J. R. Bergervoet, P. C. van Campen, W. A. van derSanden, and J. J. de Swart, Phys. Rev. C , 15 (1988).[41] I. Slaus, Y. Akaishi, and H. Tanaka, Physics Reports , 257 (1989).[42] S. Aoki et al. (Flavour Lattice Averaging Group), Eur.Phys. J. C , 113 (2020), arXiv:1902.08191 [hep-lat].[43] Z. S. Brown, W. Detmold, S. Meinel, and K. Orginos,Phys. Rev. D , 094507 (2014).[44] R. A. Brice˜no, H.-W. Lin, and D. R. Bolton, Phys. Rev.D , 094504 (2012).[45] C. Alexandrou, V. Drach, K. Jansen, C. Kallidonis, andG. Koutsou, Phys. Rev. D , 074501 (2014).[46] and .[47] M. Schrock and H. Vogt, Computer Physics Communi-cations184