Two-nucleon S-wave interactions at the SU(3) flavor-symmetric point with m ud ≃ m phys s : a first lattice QCD calculation with the stochastic Laplacian Heaviside method
Ben Hörz, Dean Howarth, Enrico Rinaldi, Andrew Hanlon, Chia Cheng Chang, Christopher Körber, Evan Berkowitz, John Bulava, M.A. Clark, Wayne Tai Lee, Colin Morningstar, Amy Nicholson, Pavlos Vranas, André Walker-Loud
LLLNL-JRNL-813871, RIKEN-iTHEMS-Report-20, MITP/20-055
Two-nucleon S-wave interactions at the SU (3) flavor-symmetric point with m ud (cid:39) m phys s :a first lattice QCD calculation with the stochastic Laplacian Heaviside method Ben H¨orz, Dean Howarth,
2, 1
Enrico Rinaldi,
3, 4
Andrew Hanlon, Chia Cheng Chang ( 張 家 丞 ),
4, 6, 1
Christopher K¨orber,
7, 6, 1
Evan Berkowitz, John Bulava, M.A. Clark, Wayne Tai Lee, Colin Morningstar, Amy Nicholson,
13, 1
Pavlos Vranas,
2, 1 and Andr´e Walker-Loud
1, 6, 2 Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA Physics Division, Lawrence Livermore National Laboratory, Livermore, CA 94550, USA Arithmer Inc., R&D Headquarters, Minato, Tokyo 106-6040, Japan RIKEN iTHEMS, Wako, Saitama 351-0198, Japan Helmholtz Institute Mainz, 55099 Mainz, Germany;GSI Helmholtzzentrum f¨ur Schwerionenforschung, 64291 Darmstadt, Germany Department of Physics, University of California, Berkeley, CA 94720, USA Institut f¨ur Theoretische Physik II, Ruhr-Universit¨at Bochum, D-44780 Bochum, Germany Maryland Center for Fundamental Physics, University of Maryland, College Park, MD 20742, USA CP -Origins & Dept. of Mathematics and Computer Science, Universityof Southern Denmark Campusvej 55, 5230 Odense M, Denmark NVIDIA Corporation, 2701 San Tomas Expressway, Santa Clara, CA 95050, USA Department of Statistics, Columbia University, New York, NY 10027, USA Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA Department of Physics and Astronomy, University of North Carolina, Chapel Hill, NC 27516-3255, USA (Dated: September 25, 2020 - 0:31)We report on the first application of the stochastic Laplacian Heaviside method for computingmulti-particle interactions with lattice QCD to the two-nucleon system. Like the Laplacian Heavi-side method, this method allows for the construction of interpolating operators which can be used toconstruct a set of positive definite two-nucleon correlation functions, unlike nearly all other applica-tions of lattice QCD to two nucleons in the literature. It also allows for a variational analysis in whichoptimal linear combinations of the interpolating operators are formed that couple predominantly tothe eigenstates of the system. Utilizing such methods has become of paramount importance in orderto help resolve the discrepancy in the literature on whether two nucleons in either isospin channelform a bound state at pion masses heavier than physical, with the discrepancy persisting even inthe SU (3)-flavor symmetric point with all quark masses near the physical strange quark mass. Thisis the first in a series of papers aimed at resolving this discrepancy. In the present work, we employthe stochastic Laplacian Heaviside method without a hexaquark operator in the basis at a latticespacing of a ∼ .
086 fm, lattice volume of L = 48 a (cid:39) . m π (cid:39)
714 MeV. Withthis setup, the observed spectrum of two-nucleon energy levels strongly disfavors the presence of abound state in either the deuteron or dineutron channel.
CONTENTS
I. Introduction 1II. Stochastic Laplacian Heaviside Method 3III. Lattice Calculation 4A. Correlation Functions 4B. Energy spectrum 61. The pion 62. Single nucleon analysis 63. Two-nucleon analysis 7C. Phase shift analysis 71. Deuteron channel 92. Dineutron channel 11IV. Discussion and Outlook 11A. Comparing with the literature 12Acknowledgments 13 A. Energy levels for the deuteron and dineutron 14B. Computational and algorithmic optimization 14C. Error propagation for q cot( δ ) 151. Unbiasedness in Weighted Regression 152. Assumptions for unbiasedness are not met 163. Our modified weighted least squares 16References 17 I. INTRODUCTION
Quantum Chromodynamics (QCD), the fundamentaltheory of nuclear strong interactions, encodes the inter-actions of nearly massless quarks and massless gluonswhich are confined into protons and neutrons, the nucle-ons, with a mass of O(1) GeV. These nucleons, whichform the basis of matter, have a residual strong inter-action that leads to the formation of nuclei with bind- a r X i v : . [ h e p - l a t ] S e p ing energies that are typically two orders of magnitudesmaller than this confinement scale: in the case of thedeuteron, the smallest nucleus made of one proton andone neutron, with a binding energy of B d ∼ . qq virtual pairs) 25 years ago [4]. The next calcula-tion was performed in 2006 using dynamical quarks withpion masses ranging from 350 (cid:46) m π (cid:46)
600 MeV [5].Since that time, there has been a measured growth inthe application of LQCD to systems with two or morebaryons, with the first calculation of a bound two-baryonsystem appearing in 2010 [6, 7]. However, significantchallenges, most notably the exponentially bad signal-to-noise (StoN) ratio [8], have prevented substantiveprogress: some 15 years later, even with all the growthin computing power and algorithmic advances, there arestill no computations of two-nucleon systems utilizing theL¨uscher method [9, 10] with pion masses lighter than m π ∼
300 MeV [11].On the other hand, we have seen the emergence ofLQCD calculations of light nuclei at m π ∼
800 MeV(up to A = 4) [12, 13] which have been used to matchto a pion-less effective field theory of few-nucleon inter-actions and used to calibrate and predict nuclei up to A = 6 [14]. We have also seen the development of anew method which first constructs a two-nucleon poten-tial, known as the HAL QCD potential [15–19], fromwhich a Schr¨odinger equation is solved and can then beused to predict the scattering phase shifts. If the HALQCD method can be demonstrated to be equivalent tothe L¨uscher method, it offers a promising alternative forcomputing the interactions of baryons from LQCD.However, there is controversy in the literature concern-ing the aforementioned equivalence and, in turn, thereis a discrepancy on whether or not two nucleons forma bound state at medium and heavy pion masses. Inshort, most calculations of two-nucleons that utilize theL¨uscher method observe the presence of (deeply) bounddeuteron and dineutron systems at pion masses largerthan ∼
300 MeV [13, 20–24]—with the exception of the Mainz group which found a bound dineutron to be un-likely at m π ∼
960 MeV [25]—while the HAL QCD Col-laboration, utilizing their potential method, concludesthat there are no bound states in either channel [17, 26].For a more detailed discussion of the controversy, see therecent review in Ref. [1].Some have found it tempting to think this disagree-ment is a demonstration that the HAL QCD method hasuncontrolled systematic uncertainties. However, whilethe L¨uscher method provides a rigorous mapping be-tween the finite-volume energy spectrum and the infinite-volume scattering amplitudes, there are potential unre-solved systematic uncertainties in the application of themethod, particularly in properly identifying the multi-particle energy spectrum. All applications that observethe existence of bound two-nucleon systems rely upona local hexaquark creation operator at the source anddilute, two-nucleon momentum-space annihilation opera-tors at the sink. The HAL QCD Collaboration has shownthat the extracted spectrum in many of these cases doesnot pass basic consistency checks, demonstrating thatthere are larger systematic uncertainties than have beenreported [27]. Combined with the StoN challenges andthe very small elastic scattering energy gaps, this has ledHAL QCD to speculate the calculations which observebound states have been misled by “false plateaux” in theeffective masses of the system, which can arise with non-positive-definite correlation functions [28]. These non-positive definite correlation functions re-quire the assumption that the overlap onto the eigen-states of the system are dominated by a single interpo-lating field constructed from the projection of each nu-cleon individually onto a state of definite momentum atthe sink side. Consider the center of mass (CoM) forsimplicity, such that (cid:104)
N N ( q ) | ∼ (cid:88) x , y c ( p ) (cid:104) | e i p · x e − i p · y N ( x ) N ( y ) , (1.1)where q is the relative interacting momentum, which isdetermined through the L¨uscher quantization condition; p is given by a non-interacting plane-wave momentummode allowed in the finite periodic volume, p = πL n with n vector of integers; and c ( p ) is a weight whichcan be chosen arbitrarily if a single momentum mode p dominates the overlap (in practice, the existing calcu-lations have chosen c ( p ) = 1). For weakly interactingsystems, such as I = 2 ππ scattering, this type of inter-polating fields works reasonably well as demonstrated bythe consistency between the results from NPLQCD [32]and HadSpec [33] which utilized this simplistic operatorand a full variational basis, respectively. For strongly in-teracting systems, such as the two-nucleon system, theresults in the literature are insufficient to draw a conclu-sion one way or the other as to how well this simplistic This has been challenged, but not conclusively demonstrated tobe wrong [29–31]. basis of interpolating fields couples cleanly to the spec-trum.In contrast, with a variational basis of interpolatingfields, one is not restricted to this assumption and insteadutilizes a linear combination of creation/annihilation op-erators (cid:104)
N N ( q ) | ∼ (cid:88) p (cid:88) x , y c ( p ) (cid:104) | e i p · x e − i p · y N ( x ) N ( y ) , (1.2)where now the c ( p ) coefficients are determined througha diagonalization of the set of interpolating fields usedand constrained by the numerical values of the correla-tion function. Even with the variational basis, experienceshows that it is still necessary to have a large basis of op-erators which provide sufficient overlap onto the variousstates of the system. For example, in the I = 1 ππ sys-tem, one must include operators that look both like local ρ operators (¯ qγ µ q ) as well as displaced two-pion oper-ators to obtain a spectrum that is consistent with theexpected ρ resonance [34]. A similar study of the nega-tive parity nucleon found that a non-local N π operatoris required [35]. In the case of the two-nucleon system,it could be that the hexaquark operator is importantfor coupling to a deeply bound state, as speculated inRef. [22]. In the present work we take a first step towards try-ing to resolve this discrepancy by performing the firstLQCD calculation of the two-nucleon systems in bothisospin channels using a positive-definite correlation ma-trix with a variational basis of operators. The first appli-cation of a variational basis to two-baryon systems wasapplied recently by the Mainz group to the h-dibaryonand dineutron systems [25, 36] in which significant ten-sion with the local hexaquark results from NPLQCD wasobserved [6, 13], although some tension with the HALQCD potential method exists as well [7, 26]. One possi-ble explanation for this is the use of only two dynamicalquark flavors (with a “quenched” strange quark) fromMainz, giving rise to potentially large systematic effectsin the determination of the binding energy as comparedto HAL QCD and NPLQCD. Furthermore, all calcula-tions were performed with a single lattice spacing withdifferent lattice actions.For the present work, we focus on the two two-nucleonchannels with total isospin I = 0 and 1, which we re-fer to as the deuteron and dineutron channel, respec-tively. We utilize the stochastic Laplacian Heavisidemethod [37], which is a stochastic variant of the distil-lation method [38]. We will summarize the method inSec. II. As discussed in Sec. III, our calculation strongly A recent study showed the use of hexaquark operators, at boththe source and sink, gives effective energies above threshold (andeven above effective energies utilizing two-baryon interpolatorsat the sink) in the dineutron channel, suggesting a hexaquarkoperator may not be necessary to accurately extract the spec-trum [25]. disfavors a bound state in either the deuteron or dineu-tron channel. We discuss the implications of this work aswell as the limitations and provide an outlook in Sec. IV.In order for the broader community to have confidencein the application of LQCD to nuclear physics, it is ofparamount importance to resolve the issue underlyingthe contradictory results in the literature.
II. STOCHASTIC LAPLACIAN HEAVISIDEMETHOD
A successful computation of two-nucleon energies reliesheavily on the construction of optimal operators. Unfor-tunately, this leads to several sources of computationaldifficulty. First, the six valence quarks present in two-nucleon correlation functions give rise to a large numberof Wick contractions. Next, in order to maximize over-lap onto the individual finite-volume two-nucleon states,both nucleon interpolating operators should be projectedonto a definite spatial momentum at the source and sink.Finally, two-nucleon interpolators which transform irre-ducibly under the finite-volume remnant of rotationalsymmetry require a summation over two-nucleon momen-tum combinations which transform among themselvesunder the little group.The considerations above necessitate a flexible and ef-ficient treatment of all-to-all quark propagation in whichthe quark propagator is determined between all spa-tial lattice points. The stochastic Laplacian Heaviside(LapH) method, which is employed here, has been suc-cessful for two-meson [39, 40] and meson-baryon [41] cor-relators. This method enables a particular choice ofquark smearing in which the quark fields are projectedonto the space spanned by the lowest N ev eigenmodesof the gauge-covariant three-dimensional laplace opera-tor [38]. Stochastic estimators with N r noise sources arethen introduced in this N ev × N spin dimensional LapHsubspace rather than the entire spatial lattice, signifi-cantly improving the variance [37].The stochastic estimators are improved by ‘dilu-tion’ [42], in which each stochastic field is partitionedinto N dil fields, each of which has support on a uniquesubset of the LapH subspace. For this work we employ N ev = 384, full spin dilution, and 12 LapH eigenvectorprojectors, so N dil = 4 ×
12 = 48. In order to ensureunbiased estimates of products of quark propagators, in-dependent stochastic fields are required for each valencequark line, so that estimates for the quark propagatorsare given by Q aα,bβ ( x, y ) = lim N r →∞ N r (cid:88) r,d φ ( r,d ) aα ( x ) ρ bβ ( y ) ( r,d ) ∗ (2.1) In practice only the upper two spin components are used for thecomputation of states propagating forward in time, reducing theeffective N dil by a factor two in the correlator construction. where ( r, d ) denote the noise and dilution indices respec-tively, ρ ( y ) is a stochastic combination of LapH eigenvec-tors and φ = Qρ is the result of a linear system solve,which are solved efficiently in GPU accelerated nodeswith QUDA [43, 44].The computation of two-nucleon correlation functionsis also simplified with stochastic LapH. Each two-nucleoninterpolator trasnforming irreducibly is given by O II Λ λ(cid:96) ( P ) = (cid:88) p p c II Λ λ(cid:96) (cid:96) N (cid:96) N (cid:96) (2.2)with definite isospin ( I, I ), little group irrep Λ, irrep row λ , and total momentum P . The additional identifier (cid:96) distinguishes multiple linearly independent operators ofthis type, while the labels (cid:96) , for the single nucleon op-erators denote the individual I and momenta p , . Notethat all terms have p + p = P , with the p , in differ-ent terms related via little group transformations. Thecoefficients c II Λ λ(cid:96) (cid:96) are determined according to Ref. [45]and are available upon request.When using stochastic LapH estimates for quark prop-agators, temporal correlators factorize into ‘source’ and‘sink’ functions which depend on the fields in Eq. (2.1) ata given Euclidean time separation. For single nucleons,these fields areΦ ( i ,i ,i ) (cid:96) ( p , t ) = c (Λ ,λ ) αβγ (cid:15) abc (cid:88) x e i p · x φ ( i ) aα ( x ) φ ( i ) bβ ( x ) φ ( i ) cγ ( x ) (2.3)and Ω ( i ,i ,i ) (cid:96) ( p , t ) (in which the φ ( x ) are replaced with ρ ( x )), where we have used the shorthand i k = ( r k , d k ) tocombine the noise and dilution indices.The rank-three tensors of Eq. (2.3) are contracted overthe i k to project onto definite ( I, I ) and treat all Wickcontractions for each of the terms in Eq. (2.2). In order toproduce an unbiased stochastic estimate, each of the sixvalence quark lines in a two nucleon correlation functionrequire a different r k . However, for a given set of stochas-tic sources, each permutation and combination of six r k produces a new (in principle correlated) estimate. How-ever, even using the minimal number of noise sources butmoderately increasing the number of permutations re-sults in a scaling of the statistical errors consistent withindependent measurements [46]. We use the maximalnumber of permutations of 6 noise sources, accountingfor nucleon-level symmetries, giving a total of 180 per-mutations. Further details of algorithmic improvementsand the optimization of our developed code are given inAppendix B. III. LATTICE CALCULATION
We employ an isotropic clover-Wilson action with N f = 2 + 1 dynamical fermions that matches the setupbeing used by the CLS Collaboration [47]. We have TABLE I. Bare parameters for the lattice action of the C103ensemble.ens. β V c κ u,d = κ s c sw C103 3.4 96 × generated the new C103 ensemble with m u = m d = m s (cid:39) m phys s , using the openQCD code [48] on the Blue-Gene/Q machine at LLNL (Vulcan). The lattice spacingis a (cid:39) .
086 fm [49] with a lattice extent V = 48 ×
96, pe-riodic boundary conditions in space and thermal bound-ary conditions in Euclidean time. The C103 ensemble has4 thermalized replicas (streams) of ∼
400 configurations,and each replica is started from different thermalized con-figurations and with different random seeds. Each config-uration is saved after 2 HMC trajectories of length τ = 2in molecular dynamic time units. The bare parametersof the lattice action are provided in Table I. The presentcomputation of two-nucleon correlation functions uses 4time-sources (cfr. t in Eq. (2.3)) on 802 configurationsspanning two of the replicas, for a total of 3208 time-sources. A. Correlation Functions
At low temperatures the spectral decomposition of atwo-point correlation function is given by C ij ( t ) = (cid:88) n z i,n ˜ z † j,n e − E n t , (3.1)where z i,n = (cid:104) Ω | O i | n (cid:105) is the overlap of the n th en-egy eigenstate onto the vacuum through the annihila-tion operator O i . If the creation and annihilation oper-ators come from a Hermitian-conjugate basis, this corre-lation function is positive definite such that all z i,n ˜ z † i,n = | z i,n | ≥
0. This simple fact greatly simplifies the anal-ysis of excited state contamination to the ground statecontribution in Eq. (3.1). Specifically it eliminates thepossibility of having a false plateau which could be gen-erated by opposite sign contributions to Eq. (3.1) fromthe lowest lying states in the spectrum.For single-hadron correlation functions, a calculationwhich uses local point or gaussian-smeared (Wupper-tal [50, 51]) quark sources for the hadron creation opera-tor while using momentum-space annihilation operatorsis still positive definite since translation invariance en-sures that, up to a multiplicative constant arising fromthe Fourier transform, the creation and annihilation op-erators are still Hermitian conjugate to each other. Ifthe annihilation operator of a two-hadron correlationfunction was constructed with a single total-momentumFourier transform, then it would also be positive definitefor the same reason, but it is well known that such opera-tors do not provide enough control over the eigenstates ofthe system to reliably extract the multitude of energy lev-els corresponding to the two hadrons interacting at differ-ent values of relative momentum. Therefore, two-hadroncorrelation functions are typically computed with eachof the two final-state hadrons separately Fourier trans-formed to a particular final-state momentum. Unfortu-nately, such annihilation operators are no longer Hermi-tian conjugates of the spatially local creation operators,and thus the correlation functions lose their positive def-inite quality.The sLapH (and LapH) methods allow for the con-struction of Hermitian-conjugate pairs of creation andannihilation operators in which each hadron at the sourceand sink can be separately Fourier transformed. Theadvantage is twofold: a volume-averaging effect at thesource as well as the sink, improving the stochastic pre-cision and a positive definite matrix of correlation func-tions.Another well-known feature of two-nucleon calcula-tions is that a ratio of correlation functions constructedas R ( t ) = C NN ( t ) C N ( t ) C N ( t ) , (3.2)provides the best way to estimate the interaction en-ergy. The stochastic correlation between the two-nucleonand single-nucleon correlation functions, C NN and C N ,is very strong, and the ratio R benefits from a large can-cellation of the single-hadron inelastic excited states: theeffective mass of this ratio correlation function R yields aprecise estimate of the interaction energy. However, priorto a time separation when the single-hadron correlationfunction has relaxed to the ground state, this ratio corre-lation function can be susceptible to false plateaus: theTaylor expansion of the single-hadron correlators in thedenominator leads to opposite sign contributions to theratio correlation function which are precisely the kind ofcorrections that can lead to false plateaus.To describe this feature more precisely, suppose thesingle-nucleon correlation function was described by justthe ground state and a single excited state C N ( t ) = A e − E t + A e − E t . (3.3)In this simplistic model, the two-nucleon correlationfunction would be given by C NN ( t ) = (cid:88) q B ,q e − (2 E +∆ E ( q )) t + (cid:88) ˜ q B , ˜ q e − ( E + E +∆ E (˜ q )) t + (cid:88) q (cid:48) B ,q (cid:48) e − (2 E +∆ E ( q (cid:48) )) t (3.4)where the sums over q , ˜ q and q (cid:48) run over the elasticscattering modes between two ground state nucleons, aground and excited state, and between two excited statesrespectively, as allowed by the L¨uscher quantization con-dition. The interaction energies ∆ E , ∆ E and ∆ E depend upon the relative momentum between the states( q ) and are typically much smaller than the inelastic ex-cited state energy E − E as these elastic scattering ener-gies must vanish as L → ∞ except in the case of a boundstate. The large-time behavior of the ratio correlationfunction is then approximated by R ( t ) (cid:39) b , e − ∆ E ( q ) t + b , e − ∆ E ( q ) t + b , e − ( E − E +∆ E (˜ q )) t − a e − ( E − E ) t + · · · (3.5)where the b ,n , b , and a are ratios of overlap factors.The observed near-exact cancellation of inelastic excitedstates in the ratio manifests as near-exact cancellationbetween the b , and − a terms on the second line ofEq. (3.5). Such cancellations with opposite signs canlead to false plateaus early in Euclidean time before thetime-separation in which the single-nucleon correlationfunction is saturated by the ground state. To avoid this problem, the NPLQCD Collaborationhas long advocated that a sufficient amount of statisticsshould be used such that the interaction energies can beprecisely determined without the need of relying uponthe ratio correlation function, but rather the two-nucleonand single-nucleon correlation functions can be fit inde-pendently and ∆ E ( q ) can be extracted under jackknifeor bootstrap resampling of the ground state energies suchdetermined [6, 52, 53].For many calculations, including the present one, thestatistical precision is insufficient to achieve a multi-sigma determination of the interaction energy from fitsto the two-nucleon and single-nucleon correlation func-tions separately. A simple measure of the feasibility ofsuch a strategy is whether one can use the ratio corre-lation function R ( t ) only at sufficiently late times thatthe single-nucleon has plateaued and still achieve a con-vincing energy extraction of the interaction energy [24].This is almost the case in the present calculation, but werequire the use of a few time slices (O(0 . − .
35) fm)prior to the ground state saturation of the single-nucleoncorrelation functions.The desire to leverage the positive definite nature ofthe two-nucleon correlation functions with sLapH, andthat, with the present stochastic precision, we must relyupon values of the correlation function prior to the single-nucleon being saturated by just the ground state, moti-vates the following set of correlation functions and theirparameterizations. First, we factorize the spectral de-composition by pulling out the ground state contributionas a prefactor. For a single nucleon of momentum q , we The false plateaus HAL QCD has speculated occur for the localsource, momentum-space sink correlation functions are not fromthis early time interference, but rather from non-positive definitecontributions from various elastic scattering states which pollutethe correlation function at late time, up until O(4) fm [28]. parameterize the correlation function as C N q ( t ) = z q, e − E q t (cid:16) z q,n e − ∆ E qn, t (cid:17) , (3.6)with an implicit sum over all excited states n >
0. Theground state energy is E q and∆ E qn, ≡ E qn − E q . (3.7)The ground state overlap factor is given by z q, , and the z q,n are the ratio of overlap factors of the n th state to theground state which all satistfy the bound z q,n > R ( t ) = r e − ∆ E NN t (cid:16) r l e − ∆ E NNl, t (cid:17)(cid:16) z q,n e − ∆ E qn, t (cid:17) (cid:16) z p,m e − ∆ E pm, t (cid:17) , (3.8)with implicit summations over the l , n and m excitedstates. The various new terms in this expression are • ∆ E NN = E NN − E q − E p , the ground state interactionenergy of interest for total momentum P = p + q ; • ∆ E NNl, = E NNl − E NN , the energy gap between the l th two-nucleon excited state and the two-nucleon groundstate energy. The l th energy gap can arise from eitheran elastic scattering state of the two ground state nu-cleons or when one or both nucleons are in an inelasticexcited state; • r = ( z NN ) / ( z q, z p, ), the ratio of the ground statetwo-nucleon overlap factor to the product of single-nucleon overlap factors; • r l = ( z NNl /z NN ) >
0, the ratio of the l th two-nucleonoverlap factors to the ground state two-nucleon overlapfactor, which are all positive.With this fit function, Eq. (3.8), if an equal number of“inelastic” excited states are included in the numeratoras in the denominator, as well as possibly extra “elastic”excited states, the fit function can naturally capture thecancellation of the inelastic excited states from the singlenucleon that are observed to also pollute the two-nucleoncorrelation functions at early times, without forcing thiscancellation to be exact. In the next section, we willdemonstrate the stability of the analysis with respect tothe time-range and number of states used in the analysis. B. Energy spectrum
1. The pion
We first look at the pion correlation function to esti-mate m p i . A single operator was used to construct this . . . E t min Q FIG. 1. A t min stability plot of the pion ground state masswith different N states in the fit function. See text for de-scription. correlation function which is fit to the cosh version ofEq. (3.1) to take into account wrap-around effects C π ( t ) = (cid:88) n z n z † n (cid:16) e − E n t + e − E n ( T − t ) (cid:17) . (3.9)Fig. 1 shows an N -state stability plot of the ground statepion mass versus t min with the chosen fit (given by thefilled symbol) coming from N = 3 states and t min =3. The fits were performed with a Bayesian constrainedanalysis [54], resulting in a determination of the pionmass in lattice units of m π = 0 . . (3.10)
2. Single nucleon analysis
We then move on to study the mass of single nucle-ons. The single-nucleon correlation functions were fitwith Eq. (3.6), also using a Bayesian constrained anal-ysis [54]. The ground state energy prior is estimatedfrom the long-time behavior of the effective mass and theground state overlap factor is estimated from an effectiveoverlap construction: m eff ( t ) = ln (cid:18) C N q ( t ) C N q ( t + 1) (cid:19) ,z eff N q ( t ) = (cid:104) e m eff ( t ) t C N q ( t ) (cid:105) / . (3.11)The prior central values are taken from the mean of theseat a late reference time of t = 10 and the prior widths aretaken to be 10 times the uncertanties on the relative ef-fective quantity at this time. For the excited state energysplittings, we use a log-normal distributed prior such thatthe total energies are ordered. The mean values of thepriors are estimated at twice the pion mass with a widththat comes down a little lower than the first N π p-wave . . . . . . E n = 1 n = 2 n = 3 n = 4 n = 5 . . . Q . . . . . . . . . t min . . . w FIG. 2. Stability of the single nucleon ground state at zeromomentum. The filled (black) circles are the effective massdata from the correlation function. The open squares arethe resulting ground state mass as a function of t min and thenumber of states n used in the analysis. The filled square isthe chosen fit. The vertical gray bands indicate time-regionsexcluded from this fit and the gray filled curve is the effectivemass reconstructed from its posteriors and the red horizontalband is the value of m N . scattering state. The central value of the l th state energyand the excited state overlap factors are then priored as E ql = E q + l × ∆ E q , ln(∆ E q ) = (ln(2 m π ) , . ,z q,l = (1 . , .
0) (3.12)where ∆ E q is the mean value and l = 1 is the first excitedstate. We use the notation ( p c , p w ) to represent a priorwith central value p m and width p w assuming that itsdistribution is Normal unless the prior name is ln( · ), inwhich case a log-normal distribution is assumed.In Fig. 2, we show the resulting ground state en-ergy of the nucleon at rest versus t min and the num-ber of excited states. It is sufficient to chose n = 3states (2 excited states) to fit the single nucleon asearly as t min = 2 to achieve an answer that is con-sistent with the general stability displayed. We ob-serve very similar stability of the ground state mass forall of the boosted single-nucleon correlation functionswhich are shown in the github repository accompanyingthis publication https://laphnn.github.io/nn_c103_qcotd_swave_only/ [55]. In all cases, we observe an n = 3 fit from t min = 2 is in excellent agreement withthe general stability of the ground state as well as n = 2with t min = 5. We use these two choices for our analysisand to explore systematics associated with the choice ofthe number of states and fit range. We find in latticeunits m N = 0 .
3. Two-nucleon analysis
In order to determine the two-nucleon eigenstates, acorrelation matrix, C ij ≡ (cid:104) ˆ O i ( t ) ˆ O † j ( t ) (cid:105) , is formed fromthe set of operators, ˆ O i ( t ), which have been projectedonto a given ( P , Λ). Solutions to the following General-ized Eigenvalue Problem (GEVP), C ( t d ) v n ( t d , t ) = λ n C ( t ) v n ( t d , t ) , (3.14)for given reference times t d , t , may then be used to rotatethe correlation matrices,ˆ C n ( t ) = ( v n ( t , t d ) , C ( t ) v n ( t , t d )) (3.15)to a basis consisting of linear combinations of operatorshaving optimal overlap (for the given basis) onto theeigenstates of the system.From this set of correlation functions we form the ratio R using Eq. (3.2) with single-nucleon correlators corre-sponding to momenta p n of the nearest non-interactingenergy level for a given state, n . This ratio is then fitto the functional form of Eq. (3.8) with similar Bayesianmethods as the single-nucleon case. Priors for the variousparameters are chosen as follows: • ∆ E NN : similar to the single nucleon, these are esti-mated from the effective mass of the ratio correlationfunction at a reference time of t = 10 with the priormean estimated from the mean of the effective massand a prior-width that is 10 times larger than the un-certainty of the effective mass; • ∆ E NNl, : We add two towers of excited states, one cor-responding to inelastic excited states with prior meansand widths estimated as with the single nucleon inelas-tic excited states and a second tower with energy gapsestimated to arise from elastic scattering states. Sincewe expect the GEVP to remove the low-lying elasticscattering excited states, the gap to the first excitedstate is estimated to be several levels above the groundstate with a prior width that allows it to be as small asthe first anticipated excited scattering state or as largeas an inelastic single nucleon excited state. • r : As with the single nucleon, the ground state ratiooverlap factor is estimated through an effective overlapfactor of the ratio correlation funcion, Eq. (3.11); • r l = (1 . , . C. Phase shift analysis
The L¨uscher finite-volume formalism [9, 10], and itsextension to moving frames [56, 57] and various gener-alizations [58–64], allows one to faithfully connect thefinite-volume two-particle spectrum to the correspondinginfinite-volume scattering phase shifts at the momentaassociated with those energies. The reduced hypercu-bic symmetry of the lattices, however, mixes the par-tial waves associated with spherical symmetry in infinitevolume. Thus, a system which has been projected ontothe given irreps of the hypercubic group will, in general,have non-zero overlap with an infinite number of partialwaves, and therefore a truncation in the partial wavesconsidered is required. Fortunately, at low energies, weexpect contributions from a partial wave, l , to fall off as E − l , justifying a truncation to the lowest partial wavesthat couple to a given finite-volume irrep.In this first look at NN interactions with sLapH (thepresent work), we ignore all partial-wave mixing inducedby the finite, periodic volume and restrict ourselves toconsidering the s -wave interactions (a standard choicein the field so far for two-baryons). We also restrictthe energies considered to those below the t -channel cut, q ∗ ≤ m π / q ∗ is the magnitude of the momentumof each nucleon in the center-of-mass (CoM) frame. Wewill relax these restrictions/assumptions in a forthcom-ing paper where we explore the partial wave mixing andenergies up to the inelastic pion-production threshold.For the low energies considered here, the K-matrix foreach partial wave is expected to be well described by asmooth polynomial in q ∗ , known as the effective rangeexpansion (ERE) q ∗ cot δ ( q ∗ ) = − a + 12 r q ∗ + 16 r q ∗ + · · · , (3.16)where δ ( q ∗ ) is the scattering phase shift, a is the scat-tering length, r is the effective range, and r n , n > qR (cid:28)
1, where R is the range of the potential.Under the assumption that partial wave mixingis negligible, the L¨uscher quantization condition pro-vides a one-to-one mapping between the spectrum and q ∗ cot δ ( q ∗ ), which for the s -wave is q ∗ cot δ ( q ∗ ) = 2 γL √ π Z d (cid:18) , q ∗ L π (cid:19) , (3.17)where γ is the ratio of the energy to the CoM energy, γ = E/E ∗ , and Z d is a generalized zeta function definedin Ref. [56], characterized by the boost vector d ≡ L π P . (3.18)The input values q ∗ are derived starting from the latticeextracted energies, E = (cid:112) E ∗ + | P | , where E ∗ is thenrelated to q ∗ via E ∗ = 2 (cid:113) q ∗ + m N . (3.19)There are several ways to proceed in fitting the nu-merical results to extract the ERE parameters; here, wediscuss three. In the first method, referred to as the determinant residual method, the L¨uscher quantizationcondition (truncated to some maximum partial wave andparameterized appropriately) is used directly to form theresiduals of the χ function, which is subsequently min-imized [65]. Using the quantization condition directly inthe fitting procedure is a natural way to include multiplepartial waves. A convenient feature of this method, asopposed to methods that directly solve the quantizationcondition, is that the generalized zeta functions can allbe computed once before the minimization process starts.However, one cannot avoid recomputing the covariancematrix each time the parameters are adjusted during thefit, since the model cannot be separated from the data. Insubsequent papers, we will explore this method in moredetail when we consider the partial wave-mixing inducedboth by the finite volume as well as the physical mixingof the S – D waves in the deuteron channel.The second method, which we refer to as the q cot δ method, has been common in the application to two-baryon systems under the truncated partial-wave expan-sion (also considered here). First, one converts the en-ergy levels, which are determined typically with Gaus-sian distributed noise, to values of the CoM momentumEq. (3.19) which are used to determine the phase shiftvalues through Eq. (3.17). These values of the q ∗ cot δ ( q ∗ )are then fit with the ERE Eq. (3.16), to determine thevalues of a , r and other shape parameters that describethe low-energy interactions.There is a subtlety in the q cot δ method, that to thebest of our knowledge has not been handled correctly inthe literature. As is well known, the zeta functions ap-pearing in Eq. (3.17) have non-linear dependence upon q ∗ in the typical range over which the momenta can bedetermined. This transforms the roughly Gaussian dis-tributed determination of E (and hence q ∗ ) into a highlyasymmetric distribution of q ∗ cot δ ( q ∗ ). Moreover, it iscommon to perform the ERE fit by treating q ∗ cot δ ( q ∗ )data points as having uncertainties in both the x and y directions. However, under the assumption of no par-tial wave mixing, the L¨uscher quantization condition,Eq. (3.17) provides a one-to-one mapping between the x ( q ∗ ) and y ( q ∗ cot δ ( q ∗ )) values, such that there isreally only a single variable with uncertainty. For suf-ficiently precise determinations of q ∗ values such that alinear approximation to Eq. (3.17) describes the results,treating the pairs of ( q ∗ , q ∗ cot δ ( q ∗ )) points with cor-related uncertainties is expected to faithfully reproducethe true uncertainty with the standard linear transforma-tions for handling x and y uncertainty. However, whenthe non-linearity of Eq. (3.17) is important, this methodcan produce biased results.We propose an alternative method that properly han-dles this non-linear relationship, which we refer to asthe spline/gradient method. Consider a bootstrap (BS)resampling of the values of ( X i , Y i ) = ( q ∗ i , q ∗ i cot δ ( q ∗ i ))pairs on irrep i . For a given BS sample, one can define thesquared distance between this point and the intersectionof the ERE function with the i th irrep as the distancealong the curve defined through Eq. (3.17) which we de-note Y i = f ( X i ) for convenience, with the distance givenby s i ( f (cid:48) i , z i, ˆ β , z i,bs ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) z i,bs z i, ˆ β dx (cid:113) f (cid:48) i ( X ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (3.20)where f (cid:48) i is the derivative of f along the curve, z i is ageneralized coordinate along the curve, z i, ˆ β is the loca-tion of intersection of the ERE function with the i th irrepand z i,bs is the coordinate of the bs th sample of irrep i .These BS distances can be used to construct the objec-tive function that penalizes data discrepancy between allirreps with the intersection of the ERE parameterization.Then, an uncorrelated, unweighted least square penaltyfor BS sample bs would be given by (cid:80) i s i ( f (cid:48) i , z i, ˆ β , z i,bs ) .To estimate the appropriate covariance, we leverage theDelta Method [66] that scales the covariance from X us-ing the gradient of f ( X ). For a normally distributed setof X variables of mean µ X ( ¯ X − µ X ) ⇒ N (0 , Σ X ) , with Σ X is the covariance of the X variables over theirreps, the Delta Method states for a differentiable f , thedistribution of f ( X ) follows f ( ¯ X ) − f ( µ X ) ⇒ N (cid:0) , ∇ f ( µ X ) T Σ X ∇ f ( µ X ) (cid:1) , where ∇ f ( µ X ) is the vector of f (cid:48) i over the irreps. Andthus, the correlated objective function we can minimizeto estimate the ERE parameters ( ˆ β ), for a given BS sam-ple is given by˜ χ bs = (cid:88) ij s i ( f (cid:48) i , z i, ˆ β , z i,bs ) W bs,ij s j ( f (cid:48) j , z j, ˆ β , z j,bs ) , (3.21)where the inverse “covariance matrix” for sample bs is W bs = (cid:104) ∇ f ( X bs ) T ˆΣ X ∇ f ( X bs ) (cid:105) − . (3.22)While the variance of X is fixed for each BS sample, thegradients are evaluated sample by sample. In order to es-timate the distance and gradient along the L¨uscher curve,a cubic spline is fit to each pair of values using the BSsamples. The ERE parameters are then estimated withthe resulting BS distribution of ˆ β BS . For more detail, seeAppendix C and Ref. [67].A third method we consider is essentially the “spec-trum method” described in Ref. [65], which directly min-imizes spectrum residuals, thus avoiding skewed q cot δ distributions. This method is composed into two steps:An outer step, effectively computing the spectrum as afunction of ERE parameters and an inner root-findingstep, solving for the value of q ∗ which satisfies theequation (3 .
16) = (3 .
17) for given q cot δ parameteriza-tion. This inner step performs a least-squares minimiza-tion for fixed ERE parameters that minimize the resid-ual of the predicted q ∗ values with those determined from the spectrum. The outer step computes a func-tion f ( a, r , ... ; irrep , P , n ) that returns the value of q ∗ for a given irrep at boost P of the n th principle cor-relator, which are compared against the spectrum by q ∗ , P ,n = E ∗ / − m N which is the numerical valueof q ∗ for the same state. We then minimize the χ withrespect to the ERE parameters χ = (cid:88) i,j ( f ( β ; i ) − q ∗ i )Cov − q ∗ ,ij ( f ( β ; j ) − q ∗ j ) , (3.23)where β = { a, r , . . . } and i, j are master indices runningover the combinations of irrep , P , n . The covariance isconstructed from the bootstrap distributions of q ∗ i withrespect to the bootstrap means q ∗ i ,Cov q ∗ ,ij = 1 N bs (cid:88) bs ( q ∗ i,bs − q ∗ i )( q ∗ j,bs − q ∗ j ) . (3.24)There are many other variants of extracting the physi-cal parameters from the two-particle spectrum which arediscussed in the literature, see for example Refs. [68–76].As we will show in Sec. III C 1 and III C 2, of the twomethods used in this work, the spectrum method is lesssusceptible to outliers, since the q ∗ values determinedare bound to finite intervals and have a near Gaussiandistribution following from their parent E distributions,and therefore, the resulting uncertainty on the extractedERE parameters is smaller. This is in contrast to the val-ues of q cot δ determined with Eq. (3.17) as these distri-butions become highly non-symmetric and heavy tailed.Nevertheless, the spline/gradient method we introducereproduces the same values of the ERE parameters andis less susceptible to the heavy-tailed fluctuations thanthe more standard analysis of q cot δ values one finds inthe literature.
1. Deuteron channel
To extract results for the deuteron channel, we con-sider all irreps whose lowest partial-wave contributioncorresponds to s -wave scattering of nucleons with isospin I = 0 and spin s = 1. In order to determine the spec-trum, we first perform a stability analysis of the two-nucleon correlation function as a function of t min and thenumber of “elastic” excited states used above and beyondthe n = 2 states used for the single nucleon. In Fig. 3,we show sample stability plots for fits to the N N ratiocorrelation functions in two different irreps. In all irreps,we find that the choices • N, Eq. (3.1): N states = 2, t = [5 , • NN, Eq. (3.8): N states = 2, n el = 0, t = [5 , • Good quality of fit, Q ;0 − . − . − . − . − . . . . ∆ E T n el = 0 n el = 1 n el = 2 . . . Q t min . . . w − . − . − . − . − . − . . ∆ E E n el = 0 n el = 1 n el = 2 . . . Q t min . . . w FIG. 3. Stability plot of the ground state energy in the T1girrep with d = 0 (top) and the second principal correlator inthe E irrep with d = 1 (bottom). The filled (black) circlespoints are the effective mass of the ratio correlation function,Eq. (3.8). The open squares are the resulting ∆ E g.s. energy asa function of t min and the number of “elastic” excited statesused, see the text for more detail. The filled square is thechosen fit. The vertical gray bands indicate time-regions ex-cluded from this fit and the gray filled curve is the effectivemass reconstructed from its posteriors and the red horizontalband is the value of ∆ E g.s. . • For a given t min , the highest weight w = e logGBF asmeasured by the relative Bayes Factor, see Ref. [77] forfurther discussion on this point where fits with differentamounts of data are also considered; • Consistency with the long time values of the effectivemass of the ratio correlation function.We opted to select the values of t min = 5 and n el = 0to be the same for all irreps analyzed to minimize thechance of accidentally biassing the result through a morefine-grained optimization.Stability plots for all irreps can be found with thegit repository accompanying this publication https://laphnn.github.io/nn_c103_qcotd_swave_only/ [55].In Appendix A in Table II, we list the irreps and theresulting ground state energies of the two-nucleon sys-tem and corresponding boosted single nucleons as wellas the processed values of q ∗ and q ∗ cot δ in m π unitsused in this analysis.In Fig. 4 we show the resulting values of ( q ∗ , q ∗ cot δ )from all irreps along with an ERE fit using the spline/gradient method (left) and the spectrum method(right). To cleanly display the correlated distributions of( q ∗ , q ∗ cot δ ) pairs, we bootstrap our energy results andshow the resulting 68% confidence intervals in the data.We note that the results from different irreps agree nicelywithin their respective energy ranges. The purely s -wavecontributions from each irrep are expected to be consis-tent with each other, with any discrepancies arising frommixing of higher partial waves. The smooth q cot δ be-havior taken from multiple irreps thus gives some confi-dence that mixing from higher partial waves is negligiblewithin our errors.We find that the fits to the ERE to q ∗ , next-to-leading order (NLO) and q ∗ , next-to-next-to leading or-der (NNLO) give consistent results for the phase shiftwithin our energy range at our given uncertainties. This,coupled with the smooth behavior of the data, stronglyindicates a convergence of the expansion within the en-ergies considered. Our results for the effective range pa-rameters are: method order m π a m π r m π r q cot δ NLO − . +3 . − . ) 5 . +1 . − . ) − spec NLO − . .
6) 5 . − q cot δ NNLO − . +3 . − . ) 5 . +2 . − . ) 2( +56 − )spec NNLO − . .
7) 4 . .
3) 29(37) (3.25)
Using the NLO ERE expansion, one can solve aquadratic equation for solutions of q cot δ = iq , (3.26)resulting in the two solutions q ± m π = im π r (cid:18) ± (cid:114) − r a (cid:19) . (3.27)Taking the results from the more stable spectrum analy-sis, the plus solution is found to be q + m π = i . . (3.28)In principle, this could correspond to a bound state so-lution. However, this solution lies well outside the rangewhere our results are constraining the amplitude (it isthe crossing of our q cot δ and − (cid:112) − q at large, negativevalue of q ). However, this can not be a physical boundstate as the slope of the q cot δ curve is larger than thetangent of the − iq curve at this crossing, as discussed indetail in Ref. [27]. The negative solution q deuteron − m π = − i . − . − .
05 0 .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . q /m π − . − . . . . . . . . q c o t δ / m π T1g( P = 0)A2( P = 1) E( P = 1)E( P = 3) E( P = 4)A2( P = 2) A2( P = 3)A2( P = 4) B1( P = 2)B2( P = 2) − . − .
05 0 .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . q /m π − . − . . . . . . . . q c o t δ / m π T1g( P = 0)A2( P = 1) E( P = 1)E( P = 3) E( P = 4)A2( P = 2) A2( P = 3)A2( P = 4) B1( P = 2)B2( P = 2) ( q cot δ analysis) (spectrum analysis)FIG. 4. Phase shift analysis, ignoring all partial wave mixing, in the irreps that overlap with the S-wave deuteron. Thespline/gradient method is shown on the left and the spectrum method on the right. The smaller (magenta) band is the 1-sigmaresult of the NLO ( q ∗ ) order fit while the larger (gray) band is the 1-sigma result of the NNLO ( q ∗ ) order analysis. The cyanline is the solution of q cot δ = iq where a bound state would occur if it were in the spectrum. solution would move towards zero and become a positiveimaginary solution which is the bound state. There havebeen few previous identifications of virtual states withlattice QCD in the two-meson sector [78, 79].Such a bound state solution would have a positive scat-tering length, such that the intercept of q cot δ at thresh-old ( q = 0) would be negative. This implies one shouldfind negative values of q cot δ for small, positive q , whichwe do not find with our results. Thus, the rsults of thiscomputation strongly disfavor the existence of a bounddeuteron at this pion mass, and with this particular ac-tion at finite lattice spacing.These results are not sufficient to rule out a boundstate in the system. For example, the operator basis wehave chosen, which does not include a hexaquark opera-tor, may not have sufficient overlap with a bound state tocorrectly extract energy levels. If this were the case, allof our results would have to systematically shift down-wards by several sigma with the inclusion of this other-wise missing operator. In a forthcoming publication, wewill investigate the impact of including such a hexaquarkoperator in the basis, which has yet to be included dueto its numerical cost.
2. Dineutron channel
For the dineutron channel, we similarly chose all irrepscorresponding to I = 1, s = 0, having overlap onto the s -wave as the leading contribution at low energies for valuesof q ∗ < m π /
2. After performing a stability analysis, wealso observe the same choice of t min = 5 and n el . = 2provides an optimal or near-optimal fit for all irreps. Theirreps and resulting energies and processed values of q ∗ and q ∗ cot δ are given in Table III in Appendix A.The resulting values of q ∗ cot δ also suggest minimalpartial wave-mixing and a smooth q ∗ dependence. TheERE analysis with the two methods described above is displayed in Fig. 5 and yields the parameters method order m π a m π r m π r q cot δ NLO − . +3 . − . ) 8 . +4 . − . ) − spec NLO − . .
0) 8 . . − q cot δ NNLO − . +3 . . ) 7 . +5 . − . ) 14( +117 − )spec NNLO − . .
0) 8 . . − Similar to the deuteron, the results are consistent withno bound state and a virtual bound state at q dineutron − m π = − i . . (3.31)Taken together, our results, while not conclusive,strongly disfavor the presence of a bound state in eitherthe deuteron or dineutron channel. IV. DISCUSSION AND OUTLOOK
We have presented the first lattice QCD calculationof two-nucleon systems using the stochastic LaplacianHeaviside method [37]. There are only two such two-baryon calculations using a variational operator basis inthe literature, the other being an application to the h-dibaryon system and the dineutron system [25, 36] us-ing the more common “distillation” method [38]. In thiswork and Refs. [25, 36], the pion mass is rather heavy.In our case it is set approximately equal to the physicalstrange quark mass resulting in m π ∼
714 MeV, and inRef. [25] it corresponds to m π ∼
960 MeV in an N f = 2calculation. Ref. [25] also performed calculations at pion masses as low as m π ∼
436 Mev, but reliable fits to the phase shift were unattain-able. − . − .
05 0 .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . q /m π − . − . . . . . . . . q c o t δ / m π A1g( P = 0) A1( P = 1) A1( P = 2) A1( P = 3) A1( P = 4) − . − .
05 0 .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . q /m π − . − . . . . . . . . q c o t δ / m π A1g( P = 0) A1( P = 1) A1( P = 2) A1( P = 3) A1( P = 4) ( q cot δ analysis) (spectrum analysis)FIG. 5. Same as Fig. 4 except for the s -wave dineutron channel. Even at the SU (3) flavor symmetric point with a heavypion mass, the pion is the lightest propagating degree offreedom emerging from QCD. It is therefore natural tomeasure other length scales with respect to the pion massand possibly natural to expect that the range R of thepotential would approximately be given by m − π . Whilethe scattering length can take on any value (with theunitary limit, a → ∞ , being the crossover between BECand BCS like systems), the effective range is typically thesize of the potential.In the present work, we have found that in both thedineutron and deuteron channels, the effective rangeis r m π ∼ −
9, which is an unusually large value.The delta-nucleon mass splitting is another small en-ergy scale, and it is found that the splitting decreaseswith increasing pion mass at a mild rate such that it is m ∆ − m N (cid:39)
200 MeV at the SU (3) flavor symmetricpoint near the physical strange quark mass [80]. Whilethis is a small energy scale compared to m π , it is notclear how this translates into a range of the two-nucleonpotential, but one should keep in mind that QCD doesnaturally produce such an energy scale. NPLQCD sim-ilarly found large values of the effective range in theircalculations at m π (cid:39)
800 MeV [21, 24], though not quiteas large.Causality and unitarity can be used to place a boundon the size of the effective range in terms of the rangeof the potential [81] with corrections arising from a finitescattering length [82] r ≤ (cid:20) R − R a + R a (cid:21) . (4.1)Since we do not know the range of the potential, R , asit is dynamically generated by QCD (and it is not anobservable), we can invert this relation and use our de-termination of the effective range to place a lower boundon R . Using the NLO ERE parameters from the spec-trum fit of the deuteron, the real solution of Eq. (4.1)provides the limit m π R (cid:38) . , R (cid:38) .
55 fm , which is roughly the same or larger than the size of thenucleon: as the pion mass increases, the pion cloud ofthe nucleon shrinks till the size of the nucleon roughlycorresponds to a size r N ∼ Λ − , similar to this value.Perhaps the range of the potential is set by the nucleonscoming “into contact” with each other. A. Comparing with the literature
Several groups [13, 20–25] have used the L¨uschermethod to compute the scattering phase shifts of thetwo-nucleon systems, deuteron and dineutron, at pionmasses larger than 300 MeV. In all cases, except theMainz group, they have found (deeply) bound states witha reasonable degree of certainty. On the other hand, theHAL QCD Collaboration [17, 26] has used their poten-tial method to conclude that there are no bound states(again at higher than physical pion masses). Below wediscuss possible sources of discrepancy.We will focus our comparison with the results fromNPLQCD at m π ∼
800 MeV [21, 24] as their results arethe most similar to ours also being at the SU (3) flavorsymmetric point near the physical strange quark mass. They have found that both the deuteron and dineutronchannels form bound states with a relatively large bind-ing energy of B (cid:39)
20 MeV at the SU (3) flavor symmet-ric point. In Fig. 6 we show our present determinationof q cot δ in the deuteron channel along with the valuesfrom Ref. [24].As is clearly visible from the figure, the results fromNPLQCD and the present work are not compatible witheach other: In order to have a bound state, there mustbe negative values of q cot δ at positive values of q inorder for the ERE to cross the − (cid:112) − q line with a slope NPLQCD also has results at m π ∼
450 MeV with bound states,however, these results are self-inconsistent as pointed about byHAL QCD [27] as well as with a low-energy scattering analy-sis [83], and so we do not compare with them. − . − .
05 0 .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . q /m π − . − . . . . . . . . q c o t δ / m π m π ∼ : single stout m π ∼ : CLS FIG. 6. Values of q cot δ in the deuteron channel in thepresent work compared with NPLQCD results at m π ∼
800 MeV. While there is some communal expectation thatdiscretization effects should be sub-dominant, one should becautious to note that these computations have been per-formed with different lattice actions at only a single latticespacing each: in the present case, the CLS clover-Wilson ac-tion [47] and with NPLQCD, a single-stout smeared [84], tad-pole improved [85] clover-Wilson action [13]. Assuming thediscretizaton effects are relatively small, and that the phaseshift does not have a strong pion mass dependence, these re-sults are in conflict. smaller than the tangent to this line [27]. Further, wehave no evidence of such large negative values of q asdoes NPLQCD (the clustering of (green) points around q /m π (cid:39) − . • There are larger systematic uncertainties in the HALQCD method and/or the local two-nucleon creation op-erators typically used by NPLQCD and Yamazaki etal. than are currently understood. HAL QCD hasspeculated that these calculations suffer from “falseplateaus” that arise through an unfortunate linearcombination of elastic scattering states [28]; • Our variational calculation has not utilized a local hex-aquark creation/annihilation operator. It is possiblethat such a local operator may couple to a deep boundstate with a significantly larger overlap such that, with-out it, the operator basis is not sufficient to identifythe state. If this were the case, the addition of thehexaquark operator would have to shift the resultingspectrum in all the irreps presented in this work down in a coordinated way that does not spoil the otherwisevery smooth q dependence observed. • None of the calculations have been performed withmore than a single lattice spacing and so there could belarger-than-expected discretization effects which pro-hibit the rigorous identification (or exclusion) of boundstates.All of these sources of potential systematic uncertaintymay need to be explored in more detail to resolve thediscrepancy. A good first start would be the use of allmethods in the literature on the same set of gauge config-urations such that one could eliminate all the systematicuncertainties aside from the method in the determina-tion of the spectrum. Provided the dispersion relation iscontinuum like in the range of momentum considered, theL¨uscher method remains valid. HAL QCD has performedthe computation of the ΞΞ spectrum and interactions us-ing their potential method and the L¨uscher method witha local source [87], which led them to conclude the localsource method has elastic excited state pollution leadingto a false plateau.In a forthcoming publication, we will compare and con-trast our present work (with more statistics) to the localsource method used by NPLQCD and Yamazaki et al.,as well as the displaced nucleons used by CalLat. Withboth methods, we will implement the HALQCD poten-tial approach such that we can isolate possible sourcesof systematic uncertainty arising in the method. We alsoplan to implement a hexaquark operator into the basis tosee how much it shifts the spectrum, if at all. Resolvingthis discrepancy is critical if we are to have confidencein the application of LQCD to multi-nucleon systems,and more importantly for the NP and HEP long-rangescience goals, to be able to compute the response of few-nucleon systems to SM and BSM currents. NPLQCDhas invested significant effort in computing such matrixelements, see for example the recent review [88], but ifthe spectrum has been misidentified, it is not clear howmuch these systematic uncertainties would modify theirresults and conclusions.
ACKNOWLEDGMENTS
We would like to thank Ra´ul Brice˜no, Evgeny Epel-baum, Daniel Phillips and Bira van Kolck for helpful cor-respondence regarding the parameterization of the two-nucleon amplitudes.Computing time for this work was provided throughthe Innovative and Novel Computational Impact on The-ory and Experiment (INCITE) program and the LLNLMultiprogrammatic and Institutional Computing pro-gram for Grand Challenge allocations on the LLNL su-percomputers. This research utilized the NVIDIA GPU-accelerated Summit supercomputer at Oak Ridge Lead-ership Computing Facility at the Oak Ridge NationalLaboratory, which is supported by the Office of Science4
TABLE II. Energy levels and phase shifts for the deuteron channel. These are determined with a 2-exponential fit to boththe single-nucleon ( t = [5 , t = [5 , d , Eq. (3.18). The state indicates the resulting principle correlation functionafter performing the GEVP. d , are the squared boost vectors of the individual nucleons used in the ratio correlation functionwhich is fit to determine the interacting energy ∆ E NN and total energy E NN which is converted to the CoM frame E ∗ NN andprocessed to get q ∗ cot δ . d irrep state d E d E ∆ E NN E NN E ∗ NN q ∗ /m π q ∗ cot δ/m π of the U.S. Department of Energy under Contract No.DE-AC05-00OR22725 as well as the Lassen (NVIDIA-GPU), and Vulcan (BG/Q) supercomputers at LawrenceLivermore National Laboratory.The computations were performed using the chroma laph and last laph software suites. chroma laph uses the USQCD chroma [89] library andthe QDP++ library. The propagator solves were efficientlyperformed with
QUDA [43, 44] and much of the sLapHworkflow has been ported to
QUDA as well. The contrac-tions were optimized with contraction optimizer [90].The computations were managed with
METAQ [91, 92].The correlation function analysis was performed with lsqfit [93] and gvar [94] and Sigmond. The resultingcorrelation functions will be released in with a futurepublication that includes a larger number of correlation functions with a full partial-wave analysis. The phaseshift analysis code and resulting bootstrap results ofthe data presented here are included with the githubrepository [55].This work was supported by the NVIDIA Corpo-ration (MAC), the Alexander von Humboldt Founda-tion through a Feodor Lynen Research Fellowship (CK),the RIKEN Special Postdoctoral Researcher Program(ER), the U.S. Department of Energy, Office of Sci-ence, Office of Nuclear Physics under Award Num-bers DE-AC02-05CH11231 (CCC, BH, AWL), DE-AC52-07NA27344 (DH, PV), DE-FG02-93ER-40762 (EB); theNuclear Physics Double Beta Decay Topical Collabora-tion (AN, AWL); and the DOE Early Career AwardProgram (AWL). CJM acknowledges support from theU.S. NSF under award PHY-1913158.
Appendix A: Energy levels for the deuteron and dineutron
The deuteron and dineutron irreps, extracted energies and processed values of q and q cot δ are provided in Table IIand Table III respectively. Appendix B: Computational and algorithmic optimization
Several of the kernels required in the stochastic LapH workflow have been implemented to run on NVIDIA V100GPUs using the QUDA library [43]. We constructed new routines that compute the cross product and contractionof color vectors, as well as specialized routines that compute time-slice reductions. Due to the reduction strategywe employ, the bulk of these contractions is expressed in terms of BLAS3 (matrix-matrix) operations, automaticallyimproving the arithmetic intensity of the computation. These operations take the form C i = A B i (B1)5 TABLE III. Energy levels and phase shifts for the dineutron channel. These are determined with a 2-exponential fit to boththe single-nucleon ( t = [5 , P irrep state N E N E ∆ E NN E NN E cm NN q /m π q cm cot δ/m π for dense matrices A, B, C , and batch index i . The matrix A is constant with respect to the batch index. As such,we wrote interfaces in QUDA for the cuBLAS function stridedBatchZGEMM which minimizes data-transfer latencybetween the host and device by caching the A matrix. We found that by using the 4 V100 accelerators on a singleLassen (LLNL) node we gained speed-up factors of ∼
30x over using the host IBM Power-9 CPUs. All the code weconstructed specifically for this computation is publicly available in the QUDA GitHub repository. The contraction,time-slice reductions, and cuBLAS interface components are in mainline QUDA.Overall, correlation function construction requires a large number of tensor contractions, and therefore dedicatedcomputational optimizations are vital. A total of 32,960 correlation functions is computed on each gauge configuration.A strategy to minimize the amount of computational work by optimizing the contraction order as well as re-usingcommon subexpressions is described in Ref. [95]. Following the nomenclature of that reference, 200,370,960 diagramsare left to evaluate after consolidating duplicates in the initial set of 2,052,792,360 diagrams. Eliminating commonsubexpressions in the set of remaining diagrams reduces the number of computationally dominant contractions with N scaling from 344,163,600 to 9,969,360 for a combined speedup by roughly a factor 350x compared to the naiveevaluation of all tensor contractions. Appendix C: Error propagation for q cot( δ )
1. Unbiasedness in Weighted Regression
In this section we offer a short refresher on the guarantees and assumptions behind classic regression. The statisticsregression model specifies the relationship between data (
X, Y ) and the associated parameters of interest β via anadditive error: Y = Xβ + (cid:15) (C1)where Y is a n × X is a n × p matrix, β is a p × (cid:15) is a n × n is thenumber of data points and p is the number of coefficients in the regression. The usual inference task is to infer β usingthe noisy observed values of Y . Any estimate for β is often denoted as ˆ β and an estimate is unbiased if E ( ˆ β ) = β .The weighted regression estimate minimizes the squared loss between the Y values and the inferred line X ˆ β ,ˆ β reg = arg min (cid:107) W / ( Y − X ˆ β ) (cid:107) = ( X T W X ) − X T W Y . We will show that ˆ β reg is unbiased if Eq. (C1) and thefollowing assumption holds: E ( (cid:15) | X ) = 0 (C2)The addition of mean 0 error in the Y dimension intuitively justifies why minimizing the squared loss in Y alonewould yield the unbiased results. Mathematically the derivation further highlights the assumption on treating X asgiven implied in both Eq. (C1) and Eq. (C2) E ( ˆ β reg | X ) = E (( X T W X ) − X T W Y | X ) (C3)= E (( X T W X ) − X T W ( Xβ + (cid:15) ) | X ) (C4)= β + ( X T W X ) − X T W E ( (cid:15) | X ) (C5)= β (C6)(C7)6 FIG. 7. Plotting the bootstrap samples within a single irrep, demonstrating the difference between the modified distance vsthe usual regression
The most common choice for W = Σ − where Cov ( (cid:15) | X ) = Σ but its choice affects the variance for ˆ β instead of itsunbiasedness. It is also worth pointing out that Normality nor symmetry in (cid:15) are necessary for the unbiasedness tohold.
2. Assumptions for unbiasedness are not met
In this section we explain the possible errors in fitting the classic weighted least squares for our curve fitting exerciseat hand. The first violation is the existence of error in the horizontal axis for each data point. Given the derivationsabove, if X is measured with error as well, our objective would incorporates both error in X and Y instead of solelyminimizing errors in the vertical axis.The second deviation from the classic setting is the strong non-linear relationship between the errors in X and Y within each irrep. Although we have errors in both axes, knowing the error in one dimension allows us to infer theerror in the other dimension. The seemingly 2 dimensional error is therefore more appropriately modeled as havinga single source of variability.
3. Our modified weighted least squares
Our modification essentially converts the problem at hand into the classical settings by re-defining the error interms of squared distance along the curve rather than the vertical axis. Fig. 7 demonstrates this modified distancebetween a data point and any candidate regression line, implied by ˆ β is the distance along the blue curve between thered points X ∗ i, ˆ β to X i,j . Classic regression on the other hand computes the vertical distance between Y i,j and X i,j ˆ β which forces the regression line towards implausible values.Let f i denote the functional relationship between X i, · and Y i, · for irrep i so that Y i, · = f i ( X i, · ). The generalequation for the length along a differentiable curve between 2 points is s ( f (cid:48) i , a, b ) = | (cid:82) ba (cid:112) f (cid:48) i ( X ) dx | . Let thepoint of intersection between f i and the regression line will be denoted as ( X ∗ i, ˆ β , X ∗ i, ˆ β ˆ β ). An unweighted least squarepenalty would then be (cid:80) i s ( f (cid:48) i , X i, ˆ β , X i,j ) , where j denotes the index for different bootstrap values.To weigh the different irreps, we estimate the covariance using a similar approach as in the Delta Method [66].The Delta Method states that if ( ¯ X − µ X ) ⇒ N (0 , Σ X ), then for a differentiable f we have f ( ¯ X ) − f ( µ X ) ⇒ N (0 , ∇ f ( µ X ) T Σ X ∇ f ( µ X )). We recycle the same f (cid:48) from calculating the lengths before, and we estimate the covari-7ance of X treating each irrep as a different dimension, specifically:ˆΣ X = 1 n − p (cid:88) j X ,j − ¯ X ... X k,j − ¯ X k (cid:2) X ,j − ¯ X . . . X k,j − ¯ X k (cid:3) (C8) ∇ f ( X · ,j ) = f (cid:48) ( X ,j )) . . . . . . f (cid:48) k ( X k,j ) (C9) W j = (cid:104) ∇ f ( X · ,j ) T ˆΣ X ∇ f ( X · ,j ) (cid:105) − (C10)where ¯ X i = n (cid:80) j X i,j is the average over the bootstraps, k is the total number of irreps, and W is the weights we willuse to scale the modified errors. It is worth noting that the Delta Method is more practical than directly estimatingthe covariance empirically because the the distance along the curve between several bootstrap samples on certainirreps are infinite if they are across the critical point of the cot function. These points make empirically estimatingthe covariance of Y infeasible.Our optimization for each bootstrap j , is to solveˆ β j = arg min β L ( f (cid:48)· , ˆ β, X · ,j ) T W j L ( f (cid:48)· , ˆ β, X · ,j ) (C11)where L ( f (cid:48)· , ˆ β, X · ,j ) = s ( f (cid:48) , X , ˆ β , X ,j )... s ( f (cid:48) k , X k, ˆ β , X k,j ) .The estimation of f i and f (cid:48) i were done using splines which are piece-wise cubic polynomials that ensure the derivativeis continuous. [1] Christian Drischler, Wick Haxton, Kenneth McElvain,Emanuele Mereghetti, Amy Nicholson, Pavlos Vranas,and Andr´e Walker-Loud, “Towards grounding nuclearphysics in QCD,” (2019), arXiv:1910.07961 [nucl-th].[2] Ingo Tews, Zohreh Davoudi, Andreas Ekstr¨om, Jason D.Holt, and Joel E. Lynn, “New Ideas in ConstrainingNuclear Forces,” (2020), arXiv:2001.03334 [nucl-th].[3] Vincenzo Cirigliano, William Detmold, Amy Nichol-son, and Phiala Shanahan, “Lattice QCD In-puts for Nuclear Double Beta Decay,” (2020),10.1016/j.ppnp.2020.103771, arXiv:2003.08493 [nucl-th].[4] M. Fukugita, Y. Kuramashi, M. Okawa, H. Mino,and A. Ukawa, “Hadron scattering lengths in latticeQCD,” Phys. Rev. D , 3003–3023 (1995), arXiv:hep-lat/9501024.[5] S.R. Beane, P.F. Bedaque, K. Orginos, and M.J. Savage,“Nucleon-nucleon scattering from fully-dynamical latticeQCD,” Phys. Rev. Lett. , 012001 (2006), arXiv:hep-lat/0602010.[6] S.R. Beane et al. (NPLQCD), “Evidence for a BoundH-dibaryon from Lattice QCD,” Phys. Rev. Lett. ,162001 (2011), arXiv:1012.3812 [hep-lat].[7] Takashi Inoue, Noriyoshi Ishii, Sinya Aoki, Takumi Doi,Tetsuo Hatsuda, Yoichi Ikeda, Keiko Murano, HidekatsuNemura, and Kenji Sasaki (HAL QCD), “Bound H-dibaryon in Flavor SU(3) Limit of Lattice QCD,” Phys.Rev. Lett. , 162002 (2011), arXiv:1012.5928 [hep-lat].[8] G. Peter Lepage, “The Analysis of Algorithms for LatticeField Theory,” in Boulder ASI 1989:97-120 (1989) pp. 97–120.[9] M. Luscher, “Volume Dependence of the Energy Spec-trum in Massive Quantum Field Theories. 2. ScatteringStates,” Commun. Math. Phys. , 153–188 (1986).[10] Martin Luscher, “Two particle states on a torus and theirrelation to the scattering matrix,” Nucl. Phys. B ,531–578 (1991).[11] Takeshi Yamazaki, Ken-ichi Ishikawa, Yoshinobu Kura-mashi, and Akira Ukawa, “Study of quark mass de-pendence of binding energy for light nuclei in 2+1 fla-vor lattice QCD,” Phys. Rev. D , 014501 (2015),arXiv:1502.04182 [hep-lat].[12] Silas R. Beane, William Detmold, Thomas C Luu,Kostas Orginos, Assumpta Parreno, Martin J. Savage,Aaron Torok, and Andre Walker-Loud, “High StatisticsAnalysis using Anisotropic Clover Lattices. II. Three-Baryon Systems,” Phys. Rev. D , 074501 (2009),arXiv:0905.0466 [hep-lat].[13] S.R. Beane, E. Chang, S.D. Cohen, William Detmold,H.W. Lin, T.C. Luu, K. Orginos, A. Parreno, M.J. Sav-age, and A. Walker-Loud (NPLQCD), “Light Nucleiand Hypernuclei from Quantum Chromodynamics in theLimit of SU(3) Flavor Symmetry,” Phys. Rev. D ,034506 (2013), arXiv:1206.5219 [hep-lat].[14] N. Barnea, L. Contessi, D. Gazit, F. Pederiva, and U. vanKolck, “Effective Field Theory for Lattice Nuclei,” Phys.Rev. Lett. , 052501 (2015), arXiv:1311.4966 [nucl-th].[15] N. Ishii, S. Aoki, and T. Hatsuda, “The Nuclear Forcefrom Lattice QCD,” Phys. Rev. Lett. , 022001 (2007), arXiv:nucl-th/0611096.[16] Sinya Aoki, Tetsuo Hatsuda, and Noriyoshi Ishii, “The-oretical Foundation of the Nuclear Force in QCD and itsapplications to Central and Tensor Forces in QuenchedLattice QCD Simulations,” Prog. Theor. Phys. , 89–128 (2010), arXiv:0909.5585 [hep-lat].[17] Noriyoshi Ishii, Sinya Aoki, Takumi Doi, Tetsuo Hatsuda,Yoichi Ikeda, Takashi Inoue, Keiko Murano, HidekatsuNemura, and Kenji Sasaki (HAL QCD), “Hadron-hadron interactions from imaginary-time Nambu-Bethe-Salpeter wave function on the lattice,” Phys. Lett. B ,437–441 (2012), arXiv:1203.3642 [hep-lat].[18] Sinya Aoki, Takumi Doi, Tetsuo Hatsuda, Yoichi Ikeda,Takashi Inoue, Noriyoshi Ishii, Keiko Murano, HidekatsuNemura, and Kenji Sasaki (HAL QCD), “Lattice QCDapproach to Nuclear Physics,” PTEP , 01A105(2012), arXiv:1206.5088 [hep-lat].[19] Sinya Aoki, Bruno Charron, Takumi Doi, Tetsuo Hat-suda, Takashi Inoue, and Noriyoshi Ishii, “Constructionof energy-independent potentials above inelastic thresh-olds in quantum field theories,” Phys. Rev. D , 034512(2013), arXiv:1212.4896 [hep-lat].[20] Takeshi Yamazaki, Ken-ichi Ishikawa, Yoshinobu Kura-mashi, and Akira Ukawa, “Helium nuclei, deuteron anddineutron in 2+1 flavor lattice QCD,” Phys. Rev. D ,074514 (2012), arXiv:1207.4277 [hep-lat].[21] S.R. Beane et al. (NPLQCD), “Nucleon-Nucleon Scatter-ing Parameters in the Limit of SU(3) Flavor Symmetry,”Phys. Rev. C , 024003 (2013), arXiv:1301.5790 [hep-lat].[22] Evan Berkowitz, Thorsten Kurth, Amy Nicholson, BalintJoo, Enrico Rinaldi, Mark Strother, Pavlos M. Vranas,and Andre Walker-Loud, “Two-Nucleon Higher Partial-Wave Scattering from Lattice QCD,” Phys. Lett. B ,285–292 (2017), arXiv:1508.00886 [hep-lat].[23] Kostas Orginos, Assumpta Parreno, Martin J. Sav-age, Silas R. Beane, Emmanuel Chang, and WilliamDetmold, “Two nucleon systems at m π ∼
450 MeVfrom lattice QCD,” Phys. Rev. D , 114512 (2015),arXiv:1508.07583 [hep-lat].[24] Michael L. Wagman, Frank Winter, Emmanuel Chang,Zohreh Davoudi, William Detmold, Kostas Orginos,Martin J. Savage, and Phiala E. Shanahan, “Baryon-Baryon Interactions and Spin-Flavor Symmetry fromLattice Quantum Chromodynamics,” Phys. Rev. D ,114510 (2017), arXiv:1706.06550 [hep-lat].[25] A. Francis, J. R. Green, P. M. Junnarkar, Ch. Miao,T. D. Rae, and H. Wittig, “Lattice QCD study of the H dibaryon using hexaquark and two-baryon interpola-tors,” Phys. Rev. D99 , 074505 (2019), arXiv:1805.03966[hep-lat].[26] Takashi Inoue, Sinya Aoki, Takumi Doi, Tetsuo Hatsuda,Yoichi Ikeda, Noriyoshi Ishii, Keiko Murano, HidekatsuNemura, and Kanji Sasaki (HAL QCD), “Two-BaryonPotentials and H-Dibaryon from 3-flavor Lattice QCDSimulations,”
Progress in strangeness nuclear physics.Proceedings, ECT Workshop on Strange Hadronic Mat-ter, Trento, Italy, September 26-30, 2011 , Nucl. Phys.
A881 , 28–43 (2012), arXiv:1112.5926 [hep-lat].[27] Takumi Iritani, Sinya Aoki, Takumi Doi, Testuo Hat-suda, Yoichi Ikeda, Takashi Inoue, Noriyoshi Ishii,Hidekatsu Nemura, and Kenji Sasaki, “Are two nucleonsbound in lattice QCD for heavy quark masses? Consis-tency check with L¨uscher’s finite volume formula,” Phys. Rev. D , 034521 (2017), arXiv:1703.07210 [hep-lat].[28] Takumi Iritani et al. , “Mirage in Temporal Correla-tion functions for Baryon-Baryon Interactions in LatticeQCD,” JHEP , 101 (2016), arXiv:1607.06371 [hep-lat].[29] Takeshi Yamazaki, Ken-Ichi Ishikawa, Yoshinobu Ku-ramashi, and Akira Ukawa (PACS), “Systematicstudy of operator dependence in nucleus calculation atlarge quark mass,” PoS LATTICE2016 , 108 (2017),arXiv:1702.00541 [hep-lat].[30] Takeshi Yamazaki, Ken-ichi Ishikawa, and YoshinobuKuramashi (PACS), “Comparison of different source cal-culations in two-nucleon channel at large quark mass,”EPJ Web Conf. , 05019 (2018), arXiv:1710.08066[hep-lat].[31] Silas R. Beane et al. , “Comment on ”Are two nucleonsbound in lattice QCD for heavy quark masses? - Sanitycheck with L¨uscher’s finite volume formula -”,” (2017),arXiv:1705.09239 [hep-lat].[32] S.R. Beane, E. Chang, W. Detmold, H.W. Lin, T.C. Luu,K. Orginos, A. Parreno, M.J. Savage, A. Torok, andA. Walker-Loud (NPLQCD), “The I=2 pipi S-wave Scat-tering Phase Shift from Lattice QCD,” Phys. Rev. D ,034505 (2012), arXiv:1107.5023 [hep-lat].[33] Jozef J. Dudek, Robert G. Edwards, Michael J. Pear-don, David G. Richards, and Christopher E. Thomas,“The phase-shift of isospin-2 pi-pi scattering from latticeQCD,” Phys. Rev. D , 071504 (2011), arXiv:1011.6352[hep-ph].[34] Jozef J. Dudek, Robert G. Edwards, and Christopher E.Thomas (Hadron Spectrum), “Energy dependence of the ρ resonance in ππ elastic scattering from lattice QCD,”Phys. Rev. D , 034505 (2013), [Erratum: Phys.Rev.D90, 099902 (2014)], arXiv:1212.0830 [hep-ph].[35] C. B. Lang and V. Verduci, “Scattering in the π N neg-ative parity channel in lattice QCD,” Phys. Rev.
D87 ,054502 (2013), arXiv:1212.5055 [hep-lat].[36] Andrew Hanlon, Anthony Francis, Jeremy Green, Parik-shit Junnarkar, and Hartmut Wittig, “The H dibaryonfrom lattice QCD with SU(3) flavor symmetry,” PoS LATTICE2018 , 081 (2018), arXiv:1810.13282 [hep-lat].[37] Colin Morningstar, John Bulava, Justin Foley, Keisuke J.Juge, David Lenkner, Mike Peardon, and Chik HimWong, “Improved stochastic estimation of quark prop-agation with Laplacian Heaviside smearing in latticeQCD,” Phys. Rev. D , 114505 (2011), arXiv:1104.3870[hep-lat].[38] Michael Peardon, John Bulava, Justin Foley, Colin Morn-ingstar, Jozef Dudek, Robert G. Edwards, Balint Joo,Huey-Wen Lin, David G. Richards, and Keisuke JimmyJuge (Hadron Spectrum), “A Novel quark-field creationoperator construction for hadronic physics in latticeQCD,” Phys. Rev. D , 054506 (2009), arXiv:0905.2160[hep-lat].[39] John Bulava, Brendan Fahy, Ben H¨orz, Keisuke J. Juge,Colin Morningstar, and Chik Him Wong, “ I = 1 and I = 2 π − π scattering phase shifts from N f = 2 +1 lattice QCD,” Nucl. Phys. B , 842–867 (2016),arXiv:1604.05593 [hep-lat].[40] Ruair´ı Brett, John Bulava, Jacob Fallica, Andrew Han-lon, Ben H¨orz, and Colin Morningstar, “Determinationof s - and p -wave I = 1 / Kπ scattering amplitudes in N f = 2 + 1 lattice QCD,” Nucl. Phys. B , 29–51(2018), arXiv:1802.03100 [hep-lat].[41] Christian Walther Andersen, John Bulava, Ben H¨orz, and Colin Morningstar, “Elastic I = 3 / p -wave nucleon-pion scattering amplitude and the ∆(1232) resonancefrom N f =2+1 lattice QCD,” Phys. Rev. D , 014506(2018), arXiv:1710.01557 [hep-lat].[42] Justin Foley, K. Jimmy Juge, Alan O’Cais, Mike Pear-don, Sinead M. Ryan, and Jon-Ivar Skullerud, “Practicalall-to-all propagators for lattice QCD,” Comput. Phys.Commun. , 145–162 (2005), arXiv:hep-lat/0505023.[43] M.A. Clark, R. Babich, K. Barros, R.C. Brower, andC. Rebbi, “Solving Lattice QCD systems of equa-tions using mixed precision solvers on GPUs,” Com-put.Phys.Commun. , 1517–1528 (2010), https://github.com/lattice/quda , arXiv:0911.3191 [hep-lat].[44] R. Babich, M.A. Clark, B. Joo, G. Shi, R.C. Brower, et al. , “Scaling Lattice QCD beyond 100 GPUs,” (2011),arXiv:1109.2935 [hep-lat].[45] C. Morningstar, J. Bulava, B. Fahy, J. Foley, Y.C. Jhang,K.J. Juge, D. Lenkner, and C.H. Wong, “Extendedhadron and two-hadron operators of definite momentumfor spectrum calculations in lattice QCD,” Phys. Rev. D , 014511 (2013), arXiv:1303.6816 [hep-lat].[46] John Bulava, Ben H¨orz, and Colin Morningstar, “Multi-hadron spectroscopy in a large physical volume,” EPJWeb Conf. , 05026 (2018), arXiv:1710.04545 [hep-lat].[47] Mattia Bruno et al. , “Simulation of QCD with N f =2 + 1 flavors of non-perturbatively improved Wilsonfermions,” JHEP , 043 (2015), arXiv:1411.3982 [hep-lat].[48] M. L¨uscher, “ openQCD : simulation programs for lat-tice qcd,” http://luscher.web.cern.ch/luscher/openQCD/ .[49] Mattia Bruno, Tomasz Korzec, and Stefan Schaefer,“Setting the scale for the CLS 2 + 1 flavor ensembles,”Phys. Rev. D , 074504 (2017), arXiv:1608.08900 [hep-lat].[50] S. Gusken, U. Low, K.H. Mutter, R. Sommer, A. Patel,and K. Schilling, “Nonsinglet Axial Vector Couplings ofthe Baryon Octet in Lattice QCD,” Phys. Lett. B ,266–269 (1989).[51] C. Alexandrou, F. Jegerlehner, S. Gusken, K. Schilling,and R. Sommer, “B meson properties from lattice QCD,”Phys. Lett. B , 60–67 (1991).[52] Silas R. Beane, William Detmold, Huey-Wen Lin,Thomas C. Luu, Kostas Orginos, Martin J. Savage,Aaron Torok, and Andre Walker-Loud (NPLQCD),“High Statistics Analysis using Anisotropic Clover Lat-tices: (III) Baryon-Baryon Interactions,” Phys. Rev. D , 054505 (2010), arXiv:0912.4243 [hep-lat].[53] S.R. Beane, E. Chang, W. Detmold, H.W. Lin, T.C. Luu,K. Orginos, A. Parreno, M.J. Savage, A. Torok, andA. Walker-Loud (NPLQCD), “The Deuteron and ExoticTwo-Body Bound States from Lattice QCD,” Phys. Rev.D , 054511 (2012), arXiv:1109.2889 [hep-lat].[54] G.P. Lepage, B. Clark, C.T.H. Davies, K. Hornbostel,P.B. Mackenzie, C. Morningstar, and H. Trottier, “Con-strained curve fitting,” Nucl. Phys. B Proc. Suppl. ,12–20 (2002), arXiv:hep-lat/0110175.[55] sLapHnn, “Phase shift analysis code for C103 s -wave nn ,” (2020), https://github.com/laphnn/nn_c103_qcotd_swave_only .[56] K. Rummukainen and Steven A. Gottlieb, “Resonancescattering phase shifts on a nonrest frame lattice,” Nucl.Phys. B , 397–436 (1995), arXiv:hep-lat/9503028.[57] C.h. Kim, C.T. Sachrajda, and Stephen R. Sharpe, “Finite-volume effects for two-hadron states in movingframes,” Nucl. Phys. B , 218–243 (2005), arXiv:hep-lat/0507006.[58] Thomas Luu and Martin J. Savage, “Extracting scatter-ing phase shifts in higher partial waves from lattice qcdcalculations,” Phys. Rev. D , 114508 (2011).[59] Ziwen Fu, “Rummukainen-Gottlieb’s formula on two-particle system with different mass,” Phys. Rev. D ,014506 (2012), arXiv:1110.0319 [hep-lat].[60] Luka Leskovec and Sasa Prelovsek, “Scattering phaseshifts for two particles of different mass and non-zero to-tal momentum in lattice QCD,” Phys. Rev. D , 114507(2012), arXiv:1202.2145 [hep-lat].[61] Maxwell T. Hansen and Stephen R. Sharpe, “Multiple-channel generalization of Lellouch-Luscher formula,”Phys. Rev. D , 016007 (2012), arXiv:1204.0826 [hep-lat].[62] M. Gockeler, R. Horsley, M. Lage, U.-G. Meissner, P.E.L.Rakow, A. Rusetsky, G. Schierholz, and J.M. Zanotti,“Scattering phases for meson and baryon resonances ongeneral moving-frame lattices,” Phys. Rev. D , 094513(2012), arXiv:1206.4141 [hep-lat].[63] Ra´ul A. Brice˜no, Zohreh Davoudi, and Thomas C. Luu,“Two-Nucleon Systems in a Finite Volume: (I) Quan-tization Conditions,” Phys. Rev. D , 034502 (2013),arXiv:1305.4903 [hep-lat].[64] Raul A. Briceno, “Two-particle multichannel systems ina finite volume with arbitrary spin,” Phys. Rev. D ,074507 (2014), arXiv:1401.3312 [hep-lat].[65] Colin Morningstar, John Bulava, Bijit Singha, Ru-air´ı Brett, Jacob Fallica, Andrew Hanlon, and BenH¨orz, “Estimating the two-particle K -matrix for mul-tiple partial waves and decay channels from finite-volume energies,” Nucl. Phys. B , 477–507 (2017),arXiv:1707.05817 [hep-lat].[66] P.J. Bickel and K.A. Doksum, Mathematical Statistics:Basic Ideas and Selected Topics , Mathematical Statistics:Basic Ideas and Selected Topics No. v. 1 (Prentice Hall,2001).[67] Andrew W. Steiner, “Two- and multi-dimensionalcurve fitting using bayesian inference,” (2018),arXiv:1802.05339 [physics.data-an].[68] Veronique Bernard, Ulf-G. Meissner, and Akaki Ruset-sky, “The Delta-resonance in a finite volume,” Nucl.Phys. B , 1–20 (2008), arXiv:hep-lat/0702012.[69] S.R. Beane, E. Chang, S.D. Cohen, W. Detmold, H.-W. Lin, T.C. Luu, K. Orginos, A. Parreno, M.J. Savage,and A. Walker-Loud, “Hyperon-Nucleon Interactions andthe Composition of Dense Nuclear Matter from QuantumChromodynamics,” Phys. Rev. Lett. , 172001 (2012),arXiv:1204.3606 [hep-lat].[70] J.M.M. Hall, A. C. P. Hsu, D.B. Leinweber, A.W.Thomas, and R.D. Young, “Finite-volume matrix Hamil-tonian model for a ∆ → Nπ system,” Phys. Rev. D ,094510 (2013), arXiv:1303.4157 [hep-lat].[71] Jia-Jun Wu, T.-S.H. Lee, A.W. Thomas, and R.D.Young, “Finite-volume Hamiltonian method for coupled-channels interactions in lattice QCD,” Phys. Rev. C ,055206 (2014), arXiv:1402.4868 [hep-lat].[72] David J. Wilson, Jozef J. Dudek, Robert G. Edwards,and Christopher E. Thomas, “Resonances in coupled πK, ηK scattering from lattice QCD,” Phys. Rev. D ,054008 (2015), arXiv:1411.2004 [hep-ph].[73] David J. Wilson, Raul A. Briceno, Jozef J. Dudek, Robert G. Edwards, and Christopher E. Thomas, “Cou-pled ππ, K ¯ K scattering in P -wave and the ρ resonancefrom lattice QCD,” Phys. Rev. D , 094502 (2015),arXiv:1507.02599 [hep-ph].[74] Jozef J. Dudek, Robert G. Edwards, and David J. Wilson(Hadron Spectrum), “An a resonance in strongly cou-pled πη , KK scattering from lattice QCD,” Phys. Rev.D , 094506 (2016), arXiv:1602.05122 [hep-ph].[75] Dehua Guo, Andrei Alexandru, Raquel Molina, andMichael D¨oring, “Rho resonance parameters fromlattice QCD,” Phys. Rev. D , 034501 (2016),arXiv:1605.03993 [hep-lat].[76] Antoni J. Woss, David J. Wilson, and Jozef J. Dudek(Hadron Spectrum), “Efficient solution of the multi-channel L¨uscher determinant condition through eigen-value decomposition,” Phys. Rev. D , 114505 (2020),arXiv:2001.08474 [hep-lat].[77] William I. Jay and Ethan T. Neil, “Bayesian model aver-aging for analysis of lattice field theory results,” (2020),arXiv:2008.01069 [stat.ME].[78] Jozef J. Dudek, Robert G. Edwards, Christopher E.Thomas, and David J. Wilson (Hadron Spectrum), “Res-onances in coupled πK − ηK scattering from quantumchromodynamics,” Phys. Rev. Lett. , 182001 (2014),arXiv:1406.4158 [hep-ph].[79] David J. Wilson, Raul A. Briceno, Jozef J. Dudek,Robert G. Edwards, and Christopher E. Thomas,“The quark-mass dependence of elastic πK scatter-ing from QCD,” Phys. Rev. Lett. , 042002 (2019),arXiv:1904.03188 [hep-lat].[80] A. Walker-Loud et al. , “Light hadron spectroscopy usingdomain wall valence quarks on an Asqtad sea,” Phys.Rev. D , 054502 (2009), arXiv:0806.4549 [hep-lat].[81] Eugene P. Wigner, “Lower Limit for the Energy Deriva-tive of the Scattering Phase Shift,” Phys. Rev. , 145–147 (1955).[82] Daniel R. Phillips and Thomas D. Cohen, “How short istoo short? Constraining contact interactions in nucleon-nucleon scattering,” Phys. Lett. B , 7–12 (1997),arXiv:nucl-th/9607048.[83] V. Baru, E. Epelbaum, and A.A. Filin, “Low-energytheorems for nucleon-nucleon scattering at M π = 450MeV,” Phys. Rev. C , 014001 (2016), arXiv:1604.02551[nucl-th].[84] Colin Morningstar and Mike J. Peardon, “Analyticsmearing of SU(3) link variables in lattice QCD,” Phys.Rev. D , 054501 (2004), arXiv:hep-lat/0311018.[85] G.Peter Lepage and Paul B. Mackenzie, “On the viabilityof lattice perturbation theory,” Phys. Rev. D , 2250–2264 (1993), arXiv:hep-lat/9209022.[86] Takashi Inoue, Sinya Aoki, Takumi Doi, Tetsuo Hatsuda,Yoichi Ikeda, Noriyoshi Ishii, Keiko Murano, HidekatsuNemura, and Kenji Sasaki, “Two-baryon potentials andh-dibaryon from 3-flavor lattice qcd simulations,” Nu-clear Physics A , 28?43 (2012).[87] Takumi Iritani, Sinya Aoki, Takumi Doi, Shinya Gongyo,Tetsuo Hatsuda, Yoichi Ikeda, Takashi Inoue, NoriyoshiIshii, Hidekatsu Nemura, and Kenji Sasaki (HAL QCD),“Systematics of the HAL QCD Potential at Low Ener-gies in Lattice QCD,” Phys. Rev. D , 014514 (2019),arXiv:1805.02365 [hep-lat].[88] Zohreh Davoudi, William Detmold, Kostas Orginos, As-sumpta Parre˜no, Martin J. Savage, Phiala Shanahan,and Michael L. Wagman, “Nuclear matrix elements from lattice QCD for electroweak and beyond-Standard-Modelprocesses,” (2020), arXiv:2008.11160 [hep-lat].[89] Robert G. Edwards and Balint Joo (SciDAC, LHPC,UKQCD), “The Chroma software system for latticeQCD,” Nucl. Phys. B Proc. Suppl. , 832 (2005),arXiv:hep-lat/0409003.[90] Ben H¨orz, “Contraction optimizer,” (2019), https://github.com/laphnn/contraction_optimizer .[91] Evan Berkowitz, “METAQ: Bundle SupercomputingTasks,” (2017), https://github.com/evanberkowitz/metaq , arXiv:1702.06122 [physics.comp-ph].[92] Evan Berkowitz, Gustav R. Jansen, Kenneth McEl-vain, and Andr´e Walker-Loud, “Job Management andTask Bundling,” EPJ Web Conf. , 09007 (2018),arXiv:1710.01986 [hep-lat].[93] Peter Lepage, “gplepage/lsqfit: lsqfit version 11.5.1,”(2020), https://github.com/gplepage/lsqfit .[94] Peter Lepage, “gplepage/gvar: gvar version 11.2,”(2020), https://github.com/gplepage/gvar .[95] Ben H¨orz and Andrew Hanlon, “Two- and three-pion finite-volume spectra at maximal isospin fromlattice QCD,” Phys. Rev. Lett.123