# Baryon-Baryon Interactions and Spin-Flavor Symmetry from Lattice Quantum Chromodynamics

## Abstract

Lattice quantum chromodynamics is used to constrain the interactions of two octet baryons at the flavor-symmetric point, with quark masses that are heavier than those in nature (equal to that of the physical strange quark mass and corresponding to a pion mass of ). Specifically, the -wave scattering phase shifts of two-baryon systems at low energies are obtained with the application of Lüscher’s formalism, mapping the energy eigenvalues of two interacting baryons in a finite volume to the two-particle scattering amplitudes below the relevant inelastic thresholds. The values of the leading-order low-energy scattering parameters in the irreducible representations of are consistent with an approximate spin-flavor symmetry in the nuclear and hypernuclear forces that is predicted in the large- limit of QCD. The two distinct -invariant interactions between two baryons are constrained at this value of the quark masses, and their values indicate an approximate accidental symmetry. The irreps containing the , and channels unambiguously exhibit a single bound state, while the irrep containing the channel exhibits a state that is consistent with either a bound state or a scattering state close to threshold. These results are in agreement with the previous conclusions of the NPLQCD collaboration regarding the existence of two-nucleon bound states at this value of the quark masses.

###### pacs:

11.15.Ha, 12.38.Gc, 12.38.-t, 21.30.Fe, 13.75.Cs, 13.85.-t.^{1}

NPLQCD Collaboration

## I Introduction

It is speculated that hyperons, the counterparts of nucleons in which some of the valence quarks in the nucleon are replaced by strange quarks, play an important role in the composition of dense matter, such as that in the interior of neutron stars (for a comprehensive review, see Ref. Page and Reddy (2006)). The interactions between two nucleons are precisely constrained by experiment over a wide range of energies. However, those between a nucleon and a hyperon, or between two hyperons, are not well known Gal and Hungerford (2005); Hashimoto and Tamura (2006); Balewski et al. (1998); Sewerin et al. (1999); Kowina et al. (2004); Bilger et al. (1998); Abdel-Bary et al. (2004); Gasparyan et al. (2004); Batty et al. (1997); Ahn et al. (2005); Holzenkamp et al. (1989); Reuber et al. (1996); Rijken et al. (1999); Haidenbauer and Meissner (2005); Rijken and Yamamoto (2006); Polinder et al. (2006); Haidenbauer et al. (2013); Beane et al. (2012a), and are challenging to probe experimentally because of the short lifetime of hyperons and hypernuclei. Precise information on how hyperons interact, in particular in a nuclear medium, is essential to establish their effects on the equation of state of dense matter and other observables. On the theoretical side, the only reliable method with which to determine these interactions is to calculate them from the underlying strong interactions among quarks and gluons described by quantum chromodynamics (QCD). This can be achieved using the non-perturbative method of lattice QCD (LQCD), which involves numerically evaluating path integrals representing Euclidean correlation functions using Monte Carlo sampling methods. This approach is taken in this work to constrain the scattering amplitudes of several classes of nucleon-nucleon, hyperon-nucleon and hyperon-hyperon systems, albeit in world that exhibits an exact flavor symmetry, with degenerate light and strange quark masses tuned to produce pions and kaons with masses of . The calculations are performed in the absence of quantum electrodynamics (QED). This work extends our previous studies of such systems using the same ensembles of gauge-field configurations Beane et al. (2013a, b), and complements previous and ongoing studies of hyperon interactions using LQCD, see for example Refs. Beane et al. (2007); Nemura et al. (2009); Beane et al. (2010, 2011a, 2011b, 2012b); Inoue et al. (2011, 2012); Beane et al. (2012a); Etminan et al. (2014); Yamada et al. (2015); Nemura et al. (2017); Doi et al. (2017); Ishii et al. (2017); Sasaki et al. (2017).

Lüscher’s finite-volume (FV) methodology M. Lüscher (1986a, 1991); Rummukainen and Gottlieb (1995); Beane et al. (2004); Feng et al. (2004); Kim et al. (2005); Davoudi and Savage (2011); Leskovec and Prelovsek (2012); Gockeler et al. (2012); Briceno et al. (2013); He et al. (2005); Hansen and Sharpe (2012); Briceno and Davoudi (2013); Li and Liu (2013); Briceno et al. (2014); Briceno (2014) is used to constrain the scattering amplitudes of two-baryon systems below the relevant inelastic thresholds from the corresponding energy eigenvalues of two interacting baryons in a finite cubic volume with periodic boundary conditions. The extraction of energies relies on the identification of low-lying states in Euclidean correlation functions. Discussions of the methods of energy determination used in this work are presented in Sec. III.2, along with careful analyses of the scattering amplitudes that result from the extracted energies in Sec. III.3. In particular, it is shown that, in agreement with our previous conclusions in Refs. Beane et al. (2013a, b), there is clear evidence for the existence of a bound state in each of the -symmetric two-baryon channels containing the , and systems. The phase shifts in these channels, and in the channels containing the system, pass all of the so-called “sanity checks” introduced in Ref. Iritani et al. (2017), contradicting the claims in that reference.

flavor symmetry has important consequences for the interactions of two octet baryons. Despite there being 64 flavor states that can be constructed from two octet baryons, symmetry dictates that there are only six independent interactions between two octet baryons, namely those in the , , , , and irreducible representations (irreps). flavor symmetry is only an approximate symmetry in nature, given the different masses of the light quarks and the strange quark, but is exact within the present numerical study, enabling a simple classification of these interactions. In particular, at leading order (LO) in an effective field theory (EFT) expansion Savage and Wise (1996), only six coefficients, corresponding to the six irreps, need be determined. Because of the structure of the interpolating fields implemented in this work, scattering information has been obtained only for channels belonging to the first four irreps listed above.

The results of our calculations allow an exploration of the spin-flavor symmetries of nuclear and hypernuclear interactions that are predicted from QCD (with degenerate flavors) in the limit of a large number of colors, Kaplan and Savage (1996). In this limit, the six LO interactions of the -symmetric low-energy theory are defined by only two independent constants, reflecting a manifest spin-flavor symmetry. Corrections to the constraints imposed by the symmetry scale as . The calculations performed in this work provide an opportunity to examine the large- relations without contamination from breaking effects that are present in nature. Sec. III.3 includes the results of this investigation, which demonstrate for the first time the -symmetric nature of interactions in the two-baryon channels (even at ), and further point to an accidental symmetry. Assuming the spin-flavor symmetry, predictions are made for the leading LECs of the effective -symmetric baryon-baryon interactions. Future studies of the two-baryon channels belonging to the and irreps are needed to confirm these conclusions.

As the values of the parameters of QCD in this study differ from those in nature, there is no direct connection between the present results and phenomenology. However, our work presents an exploration of a non-Abelian gauge theory that is continuously connected to the strong-interaction sector of nature through the variation of the masses of the light quarks. These calculations establish and verify formal, numerical and algorithmic technologies that are needed for future explorations of multi-baryon systems at the physical values of the quark masses. Additionally, the possibility of changing the parameters of QCD in LQCD studies is itself a unique feature that has been shown to reveal insights into the structure of QCD that would be impossible to discover experimentally. Perhaps most importantly, refinements of the chiral nuclear forces requires calculations over a range of quark masses. Nucleon-nucleon interactions are speculated to be finely tuned in nature, and the calculations in this work test how robust this fine tuning is with regard to changes in the quark masses Beane and Savage (2003); Bedaque et al. (2011); Chen et al. (2012); Soto and Tarrus (2012), investigations that are only possible with LQCD. The first study of the unnaturalness of nucleon-nucleon interactions using the same ensembles of gauge-field configurations as in this work has already been conducted in Ref. Beane et al. (2013b). Here, this study is extended to scattering channels involving hyperons. In addition, important progress has been made recently in applying EFTs and nuclear many-body techniques to extend the range of predictions of QCD with to heavier nuclei Barnea et al. (2015); Kirscher et al. (2015, 2017), pointing to the ground state of O being likely unbound Contessi et al. (2017). Further investigations are needed to confirm this result and study its implications on the periodic table of nuclide at heavy quark masses. Such studies can be extended to hyperon systems with the aid of the LQCD results presented in this work, and will provide more insight into the robustness of the properties of nuclear and hypernuclear systems with respect to variations in the parameters of QCD.

The rest of this paper is organized as follows. The formalism required to analyze and interpret the numerical results of this study is presented in Sec. II. In particular, Sec. II.1 contains a summary of the method used to extract the scattering amplitudes below the relevant inelastic thresholds from LQCD energy eigenvalues, along with discussions of the volume dependence of bound-state energies. Sec. II.2 summarizes the expectations of flavor symmetry for two-baryon channels, and subsequent predictions for an extended symmetry present in the limit of large . Details of the numerical study and the results are presented in Sec. III. A summary and conclusion follow in Sec. IV. The paper includes four appendices: Appendix A presents more detail on the structure of baryon-baryon systems. Appendix B tabulates the values of LO scattering amplitudes in mixed flavor channels. Appendix C contains the full tables of the results for energies and phase shifts. Finally, Appendix D is devoted to examining the so-called “sanity checks” of Ref. Iritani et al. (2017), demonstrating the definitive presence of physical bound states in the two-nucleon channels at this value of the quark masses, and the source-independence of the results presented, contrary to the claims presented in Ref. Iritani et al. (2017).

## Ii Formalism

The goal of the numerical calculations presented here is to constrain scattering amplitudes in various baryon-baryon channels. Scattering information is obtained from the energy spectra of two baryons in a finite volume, and a summary of Lüscher’s methodology for mapping finite-volume energy eigenvalues to scattering amplitudes is presented in this section. Signatures of bound states in LQCD calculations of two-baryon spectra are further discussed. This section also contains theoretical background relevant for scattering processes with symmetry and the predictions of the large- limit of QCD.

### ii.1 Two-baryon systems in a finite volume and Lüscher’s methodology

Below all relevant inelastic thresholds, the interacting energies of two particles in a finite volume determine the scattering amplitudes through a direct mapping given by Lüscher’s quantization condition (QC) M. Lüscher (1986a, 1991). This mapping is valid as long as the interactions have a finite range that is contained inside the lattice volume. For typical hadronic systems, the range of interactions is set by the Compton wavelength of the pion. This gives rise to corrections to the QC that are suppressed as , where denotes the spatial extent of a cubic volume M. Lüscher (1986b).

Two octet baryons can be in either a spin-singlet () or a spin-triplet (coupled ) state. However, states in a finite cubic volume with periodic boundary conditions can not be characterized with well-defined angular-momentum quantum numbers. As a consequence, the FV QC mixes scattering amplitudes in all partial waves, preventing the extraction of scattering parameters. At low energies, however, only the lowest partial waves are expected to be significant, and the QC can be truncated to a finite space. Therefore, for two-baryon systems in a spin-singlet state, a simple algebraic relation enables the -wave scattering phase shift, , to be accessed from the FV energy eigenvalues at low energies,^{2}

(2) |

Here, is the relative momentum of each baryon in the center-of-mass (CM) frame and denotes the total CM momentum of the system in units of . is a kinematic function related to the three-dimensional zeta function, ,

(3) |

where is the relativistic gamma factor, with and denoting the total energy of two baryons in the lab and CM frames, respectively M. Lüscher (1986a, 1991); Rummukainen and Gottlieb (1995); Christ et al. (2005); Kim et al. (2005); Briceno et al. (2013); Briceño et al. (2013). Further,

(4) |

where, for two baryons with equal masses, . denotes a triplet of integers and acting on a vector rescales the component of the vector parallel to the boost vector by while leaving the perpendicular component intact. The zeta function in Eq. (4) can be numerically evaluated most efficiently using an equivalent exponential form M. Lüscher (1991); Beane et al. (2012c); Leskovec and Prelovsek (2012). For two-baryon systems at the energies considered below, the relativistic corrections due to the deviation of the factor from unity are at the sub-percent level. As a result, for boost vectors whose components (in units of ) are equal to each other modulo a factor of 2, the corresponding QCs in Eq. (2) are approximately the same. The other consequence of the proximity to the non-relativistic (NR) limit is that for boost vectors of the form , with each being an integer, the leading contamination to the -wave QC arises from nonvanishing -wave interactions, which are expected to be suppressed relative to the -wave interactions. With boost vectors that do not take the above form, the leading contamination arises from -wave interactions Briceno et al. (2013).

For two-baryon systems in a spin-triplet state at low energies, the physical mixing between and partial waves must be taken into account. A low-energy EFT of two-baryon interactions suggests that the - mixing parameter, , contributes to the low-energy expansion of the scattering amplitude at the same order as the effective range parameter, and may not be ignored Briceño et al. (2013). In the *Blatt-Biedenharn* parametrization of a coupled-channel scattering amplitude Blatt and Biedenharn (1952), the mixing parameter has an analytic expansion in energy near the bound-state pole, and the scattering amplitude exhibits a simple condition for the location of such a pole, . Here, is the counterpart of the -wave phase shift of the *barred* parametrization Stapp et al. (1957), which has a small -wave admixture as well. The Blatt-Biedenharn parametrization will be adopted in this work for scattering in the spin-triplet channels.

The mixing parameter, , adds an extra unknown to the QC in the spin-triplet channels. Constraining this parameter, as discussed in Ref. Briceño et al. (2013), requires knowledge of the spectra of two-baryon systems with the total spin aligned both parallel and perpendicular to the boost vector, and with boost momenta that have at least one component equal to unity modulo (in units of ). As not all distinct orientations of total spin with respect to the boost momenta are constructed in forming the correlation functions of spin-triplet systems in this work, the parameter cannot be constrained here for the spin-triplet channels. This also implies that for boost vector , the corrections to the QC from the s-d mixing might be significant at the order of low-energy EFT considered, and constraints on the -wave phase shift arising from a simple -wave QC may be contaminated. On the other hand, for boost vectors of the form , with each being an integer, in particular for the and boost vectors that are considered in this work, the -wave QC,

(5) |

is exact up to corrections from -wave interactions.^{3}

Once the phase shifts are determined at several CM energies, a low-energy parametrization of the scattering amplitude as a function of energy, with only a few unknown parameters, can be constrained over a given range of energies. In the baryon-baryon channels well below the t-channel cut, the most common parametrization is the effective range expansion (ERE). For -wave (-wave) interactions, the ERE is an expansion of the function in ,

(7) |

where , and are the scattering length, effective range and the leading shape parameter, respectively. The ellipsis denotes terms that are higher order in the momentum expansion. Lüscher’s QC condition provides (up to exponentially small volume corrections and discretization effects) an exact constraint on the amplitude at corresponding energies regardless of the complexities present in the analytic structure of the amplitude below the inelastic thresholds. It is the output of the QC that allows the efficacy of given parametrizations of the amplitude to be assessed. For example, although the ERE is guaranteed to have a nonzero radius of convergence around , the convergence rate is not known a priori, and fits with higher order terms in the ERE may be needed. With numerical calculations for a range of momenta, the appropriateness of a given truncation of the ERE must be carefully tested.

Lüscher’s QC contains information about possible bound states in the system through an analytic continuation of the condition to negative energies. In particular, it is straightforward to show that for , and for boost vectors of the type ,

(8) |

in the NR limit Koenig et al. (2011); Bour et al. (2011); Davoudi and Savage (2011); Briceño et al. (2013). Here, is the infinite-volume binding momentum of the state and is the residue of the scattering amplitude at the bound-state pole.
Note that the occurrence of negative values in a system in a finite volume is not necessarily an indication of a bound state, and the movement of the state on the real energy axis must be examined as function of volume, according to the above form, to ascertain that the energy (shift) remains in the negative region towards infinite volume. Here, this will be referred to as a *direct* method to obtain the binding energy. A crucial feature of calculations performed in this work is that two-baryon systems are studied at multiple volumes in order to provide unambiguous signatures for the existence of bound states once negative-valued energy shifts are observed. In particular, for the largest volume used, with a spatial extent of , the FV corrections to the infinite-volume binding momenta are very small for the bound states in the , and irreps, see Sec. III.3. Since the closed form of the FV corrections to the binding momenta are known Koenig et al. (2011); Bour et al. (2011); Davoudi and Savage (2011); Briceño et al. (2013), the significance of the terms that are dropped from the expansion in Eq. (8) can be evaluated order by order.

Another method of obtaining information about a bound state is to first constrain the scattering amplitude and its parametrization in terms of energy using Lüscher’s methodology. An analytic continuation to negative energies then allows the bound state energy to be obtained from the pole location(s) of the scattering amplitude,

(9) |

Since this method involves an intermediate step to obtain the binding energies, it is referred to here as an *indirect* method. The advantages of this method are that it makes no assumption about the suppression of higher-order exponentials in the extrapolation form as in Eq. (8), and that it provides information about the existence or absence of a bound state even near threshold. The disadvantage of this method is that it relies on a parametrization of the scattering amplitude. Often, including additional parameters to improve the goodness of the fit increases the uncertainty of constraints on the location of the pole. Bound state(s) extracted this way must be shown to be robust against changes in the parameterization, and the scattering amplitude at the bound state energy must be shown to satisfy certain physical conditions. These features will become more apparent in Sec. III.3, where the determinations of the binding energies in the various baryon-baryon channels are discussed.

### ii.2 Two-baryon scattering with flavor symmetry and large- predictions

The number of distinct FV spectra in the baryon-baryon systems is dictated by the flavor symmetry of the present calculations. The flavor representation of two octet baryons, each transforming in the irrep of , has a decomposition of the form:

(10) |

Flavor channels belonging to the totally symmetric irreps , and have a total spin equal to zero, while those belonging to the totally antisymmetric irreps , and have a total spin equal to one. The classification of the flavor channels is summarized in Appendix A for reference. The use of interpolating operators that transform under irreps of the decomposition of the product of two octet baryons allows for these distinct spectra to be determined in a LQCD calculation. The two-baryon interpolating operators used in this study, however, transform under the isospin subgroup of , with strangeness treated as a quantum number. As a result, the excited spectra corresponding to the and irreps cannot be rigorously determined unless multiple interpolating operators in flavor space are used to isolate the lowest-lying states of the systems. For example, to obtain the energy eigenvalues beyond the ground state in the irrep, a matrix of correlation functions in flavor space must be formed from interpolating operators corresponding to spin-singlet , and states, see Fig. 18. Since such a complete basis of operators was not used to form the correlation functions Beane et al. (2013a), direct constraints on scattering amplitudes in these two irreps could not be obtained.^{4}

At low energies, the leading -wave interactions of two octet baryons can be described by a Lagrange densitySavage and Wise (1996) in a pionless EFT Chen et al. (1999) of the form,

(11) | |||||

Here, is the octet baryon matrix,

(12) |

where Roman indices on the B fields denote spin components. The Savage-Wise (SW) coefficients can be matched to scattering amplitudes at LO in a momentum expansion. For natural interactions, i.e., when the scattering length is comparable to the range of interactions, the relationships between the scattering lengths and the SW coefficients are presented in Ref. Savage and Wise (1996) for various baryon-baryon channels. However, as is known in nature, and was deduced previously for the heavy quark masses of this work Beane et al. (2013b), the -wave interactions in both two-nucleon channels appear to be unnatural. The present investigation reconfirms the unnatural nature of interactions in the two-nucleon channels (belonging to the and irreps) and further points to the similar feature in channels belonging to the and irreps. For unnatural -wave interactions, the required power counting of the amplitude is produced in the Kaplan, Savage and Wise Kaplan et al. (1998a, b) and van Kolck van Kolck (1999) (KSW-vK) schemes. The relations of Ref. Savage and Wise (1996) for unnatural scattering lengths in terms of coefficients become

(13) |

where denotes the baryon mass, and the coefficients on the right-hand side are evaluated at the renormalization scale . For natural interactions, the renormalization scale is set equal to zero in the left-hand side of these equations, corresponding to a tree-level expansion of the scattering amplitude in these couplings.

The large- limit has interesting consequences and gives rise to further simplification of the interactions of two baryons Kaplan and Savage (1996). As argued in Ref. Kaplan and Savage (1996), in the limit of flavor symmetry, the interactions among two nucleons are invariant under a spin-flavor symmetry up to corrections that scale as . Including the strange quarks and in the limit of flavor symmetry, interactions are invariant under an symmetry up to corrections that scale as . Focussing on the latter case (which contains the former case as a subgroup), it can be shown that there are only two independent dimension-six -symmetric interactions of two octet baryons, with coefficients and .^{5}

(14) |

where the coefficients on the right-hand side are evaluated at the renormalization scale . For natural interactions, is set equal to zero in the left-hand side of these equations, corresponding to a tree-level expansion of the amplitudes in these couplings. Note that the scattering lengths in channels belonging to the and are the same up to corrections. Recalling that the and states belong to the and irreps, respectively, this equality is a manifestation of the accidental Wigner symmetry Wigner (1937), indicating that the spin-dependent -wave interaction vanishes in the large- limit Kaplan and Savage (1996); Kaplan and Manohar (1997). Additionally, a larger accidental symmetry of two-baryon interactions can be realized in the limit where the coefficient is of or larger while the coefficient is of or smaller. In this case, the contributions from the coefficient to the amplitudes is suppressed relative to those of the coefficient through a numerical suppression observed in the terms in Eqs. (14). This results in an symmetry of LO interactions, with only one coefficient, , to be constrained.

Observation of spin-flavor symmetry and an accidental symmetry of nuclear and hypernuclear forces, although at an unphysical value of the quark masses, will be the first confirmation of the large- QCD predictions in hyperon-nucleon and hyperon-hyperon systems. Such investigation are presented in Sec. III.3. Identifying these spin-flavor symmetries, along with the theoretical estimate of their violation, provides important constraints on the hyperon-nucleon interactions that can be included in calculations of finite density systems. It is important to note that the calculations presented in this work exhibit an exact flavor symmetry, making the large- predictions above free of the breaking contaminations that are present in nature.

## Iii Numerical calculations and results

This section contains the main results of this paper. These include the scattering phase shifts and the constraints on the ERE parametrization of two-baryon channels belonging to the , , and irreps of the decomposition of the product of two octet baryons, as well as an investigation of spin-flavor symmetries of interactions and a subsequent accidental symmetry predicted at large . The main inputs to the scattering amplitude determinations are energy eigenvalues obtained from LQCD calculations of correlation functions. Some of these energies have been previously presented in Refs. Beane et al. (2013a, b). Here, multiple analyses are performed to determine the ground states and first excited states of two-baryon channels. The extracted energies are found to be consistent with our previous determinations.

### iii.1 Details of LQCD computations

The ensembles of gauge-field configurations and the two-baryon correlation functions used in this work have previously been analyzed to obtain binding energies in two, three and four-baryon systems Beane et al. (2013a), as well as low-energy scattering phase shifts in two-nucleon systems Beane et al. (2013b). Additionally, the same gauge-field configurations have been used to study the magnetic structure of light nuclei Beane et al. (2014); Chang et al. (2015); Detmold et al. (2016) and some of the simplest reactions in the few-nucleon systems, such as the radiative capture process Beane et al. (2015) and single and double- decays Savage et al. (2016); Shanahan et al. (2017); Tiburzi et al. (2017). Details of the ensemble generation, as well as of the construction of the nuclear correlation functions, have been presented in those works, see for example Ref. Beane et al. (2013a). Here, some of the technical details are reviewed for completeness.

The gauge-field configurations were generated using a tadpole-improved Lüscher-Weisz gauge action M. Lüscher and P. Weisz (1985) and a clover action for fermions Sheikholeslami and Wohlert (1985). The choice of stout smearing and the tadpole-improved clover coefficient used in generating the gauge configurations alleviate discretization effects to , where denotes the lattice spacing. This spacing is determined, from spectroscopy on these ensembles, to be , see Ref. Beane et al. (2013a) and references therein. The physical spatial extents of these ensembles are approximately , and . Throughout this paper, these ensembles will be referred to as: , and , respectively. The first three dimensions refer to the spatial extent of the hypercubic volume, , while the last dimension refers to the temporal extent, , both in lattice units (l.u.). The configurations are separated by ten Hybrid Monte Carlo evolution trajectories to reduce autocorrelations, with the total number of configurations used for each ensemble, given in Table 1. An average of measurements are performed on each configuration. Various other properties of the ensembles are listed in Table 1. Given the large values of and relative to the inverse pion mass, both the thermal contamination and the exponential finite-volume contamination of single-hadron masses and two-baryon energies from pion propagation through the boundaries are strongly suppressed.

Sources are smeared with a gauge-invariant Gaussian profile with stout-smeared gauge links. The quark propagators at the sink are either not smeared (smeared-point combination, SP) or are smeared with the same smearing profile as that of the source operators (smeared-smeared combination, SS). The plateau regions of the effective mass plots (EMPs) formed out of the SP and SS correlation functions are found consistent in every case. Propagators are contracted at the sink in blocks of three quarks to assemble a baryon field with given quantum numbers at the sink. In particular, the baryon blocks are projected to a fixed three-momentum, enabling the two-baryon interpolators at the sink to have either zero or non-zero CM momentum, with various possibilities for the momentum of each baryon. As the next step, a fully-antisymmetrized quark-level wavefunction with overall quantum numbers of the two-baryon system of interest is formed at the location of the source. The contraction step is defined by the selection of the appropriate indices from the baryon blocks at the source, in a way that is dictated by the quark-level wavefunction. More details regarding the contraction algorithm for a general -nucleon system are presented in Ref. Detmold and Orginos (2013) (with a similar approach proposed in Refs. Doi and Endres (2013); Günther et al. (2013)). The final products of the contraction step are two-baryon correlation functions as a function of Euclidean time. These correspond to a definite total momentum resulting from several (nearly orthogonal) choices of baryon momentum at the sink.

### iii.2 Analysis of correlation functions

To maximize confidence in the energy determinations and their uncertainties, five different analysis procedures were used, and the results obtained from each method were found to be consistent. The statistical and fitting systematic uncertainties on the final results are taken from one analysis, with an additional systematic uncertainty added to account for the small variations between the five analyses.

The correlation function of a single or two-baryon system, projected to the total momentum , can be written as

(15) |

where () denotes the overlap of the interpolating operator () onto the corresponding eigenstates of the system, with subscripts “” and “” referring to the ground state and the first excited state, respectively. and denote the ground and excited-state energies, respectively, and the ellipsis denotes contributions from additional higher-energy states. In principle, these correlation functions can be used to obtain the tower of energy eigenvalues of the system in a finite volume. In practice, a reliable determination of even the few lowest-lying energies is challenging. Since only a single source operator and two different sink operators were used for any given momentum configuration, it is not possible to use a Hermitian variational approach here. Nonetheless, given the exponential form in Eq. (15), it is clear that a linear combination of the two correlation functions can be used to remove the excited-state contamination of the lowest lying state at earlier times. Various realizations of this approach are the Matrix Prony Beane et al. (2009a, b) and the GPoF Hua and Sarkar (1989); Sarkar and Pereira (1995); Aubin and Orginos (2011) methods. Alternatively, a correlated -squared function can be formed to fit directly to single or two-exponential forms, with the correlations both in time and between the different source and sink structures accounted for. As another alternative, the effective energy function defined as

(16) |

can be fit to a constant at late times to obtain the ground-state energy, . in Eq. (16) is a non-zero integer. A detailed account of various analysis techniques employed herein has been presented in Ref. Beane et al. (2009a). While the results and the plots corresponding to a single analysis are presented below, a systematic uncertainty accounting for the small variation among the energies from different analyses is incorporated in the numbers that are reported.

Statistical uncertainties were obtained for each analysis technique using bootstrap or jackknife procedures. The systematic uncertainty in each analysis includes a fitting uncertainty obtained by allowing the fit region to vary within an acceptable window. Representative fits obtained from the primary analysis are shown in the EMPs in Figs. 1-5. These correspond to the largest fit intervals with a in each case.

The central values and uncertainties in the masses of the octet baryon obtained in this manner for all ensembles are overlaid with the SP (SS) EMPs in the center (right) panels of Fig. 1. Both the SP and SS EMPs are additionally shown at a larger scale in the left panels. The extracted values of the baryon mass at each volume are given in Table 9 in Appendix C. The masses extracted for the three different ensembles agree within uncertainties, as expected from the large values of in this calculation, and the infinite-volume value of the mass is taken to be the mass extracted for the largest ensemble, in lattice units (l.u.).

The upper panels of each segment in Figs. 2–5 show the EMPs of the two-baryon systems for both the SP and SS correlation functions. The ground-state energy associated with each correlation function is determined as described above. However, the quantity that is of most interest in two-baryon channels is the shift in the energy of the system resulting from two-body interactions. The energy of two free baryons at rest, , can be subtracted from the two-baryon energies in a correlated manner to extract this small energy shift. Another approach that retains the correlations between the single and two-baryon correlations functions, thus reducing the statistical noise, is to form the ratio

(17) |

Here, and ( and ) are interpolating operators for the single(two)-baryon system and , and are known ratios of overlap factors of given states. At late times when the exponential factors in both the numerator and the denominator of the ratio on the right-hand side of Eq. (17) are negligible compared with unity, a fit to a single exponential can be performed at large times, following the analysis steps described above, to obtain the energy shift . The effective energy-shift function associated with the ratio in Eq. (17) can be defined as

(18) |

Given the form of , flat behavior of in time is not a sufficient indicator that the function is a single exponential. The values of overlap ratios in the numerator and the denominator in Eq. (17) may conspire to give rise to flat behavior, despite neither the single-baryon nor the two-baryon systems being in their respective ground states. As a result, in fitting the quantity , none of the fit intervals must begin earlier than the beginning of the single-exponential regions in the single-baryon and two-baryon EMPs.

In principle, two-baryon correlation functions contain spectral information beyond ground-state energies. Although this study did not use a large basis of operators, physical intuition regarding the differing nature of bound and scattering states of a two-baryon system suggested constructing not only the two-baryon operators that interpolate to two baryons at rest or in motion with equal velocity, but also those that interpolate to two baryons with relative back-to-back momenta. While the former can have significant overlap onto a compact state in a finite volume (corresponding to a bound state in infinite volume), they are not optimal interpolators for states corresponding to the scattering states of infinite volume. This results in correlation functions that are dominated by the ground state after a short time interval. On the other hand, interpolators with back-to-back momenta appear to predominantly overlap with states with positive energy shifts in the volume, and are almost orthogonal to the operators of the first type. The quality of plateaus in the EMPs with both types of interpolators was found to be comparable, suggesting that each set primarily overlaps onto one state and not the other. This allows the first excited states of the two-baryon systems to be extracted using the simplest back-to-back momentum configurations for baryons. The only exception is for the ensemble, where the splitting between the energy levels of the systems is small (being comparable to the uncertainties in the energies) and it can not be established that the first excited state is only minimally mixed into the nearby ground state. As a result, while for the smaller volumes two energy levels are extracted, for the ensemble only the ground-state energies are reported.

The upper panels in each segment in Figs. 2-5 include not only the lowest-lying state, but also the second lowest-lying state of the two-baryon systems obtained from the correlation functions with back-to-back momenta. The lower panels of each segment show EMPs corresponding to the quantity for the SS or SP (depending on the channel) correlation functions, as defined in Eq. (18). The same quantity can be constructed for the correlation functions that project to the first excited state, with . In the irrep, both the fit to the SP correlation function and a correlated fit to both the SP and SS correlation functions exhibited consistent plateaus, but the fit to the SS/SS correlation function ratio was found most precise. Similarly, in other irreps, fits to the SP/SS and SS/SS correlation function ratios, as well as a correlated fit to both of these ratios, were performed and the fit corresponding to the least uncertainty was selected, as is indicated in Figs. 3–5. The energy shifts and their uncertainties are denoted as horizontal bands in the plots, and are compiled for all two-baryon channels studied in this work in Fig. 6. The corresponding values are tabulated in Tables 10-13 of Appendix C for reference.

Recently, there have been comments by Iritani, et al. Iritani (2016a); Iritani et al. (2016); Iritani (2016b) questioning the extraction of energy eigenvalues from the late-time behavior of correlation functions, and methods for identification of energies such as those used here. These authors present an example of two-nucleon correlation functions that exhibit a considerable mismatch in the location of the naive plateaus in the EMPs when different source and sink operators are used (namely locally-smeared and wall sources). However, as is shown by the PACS-CS collaboration Yamazaki et al. (2017), such a mismatch disappears once both the single-nucleon and the two-nucleon systems are required to be in their ground states. The failure of wall sources to overlap well onto the ground state at early times is a well-known problem, and has no bearing on the results reported by other groups using more optimal sources, such as those used in this work. Indeed, the quality of plateaus in the two-baryon systems are comparable to those of the single-nucleon system in the present study, demonstrating that the ground state (and the first excited state) of these systems can be obtained efficiently, with the results from two different source-sink combinations being fully consistent.

Another argument to consider when assessing the claims by Iritani et al. regarding the occurrence of so-called “mirage plateaus” in two-baryon systems follows from observations of the volume dependence of the correlation functions. Fig. 7 shows the EMPs in each of the two-baryon channels studied in this work in the three different lattice volumes. Volume dependence is clearly visible in states identified as scattering states. No significant volume dependence is observed for the lowest-lying state, strongly supporting the hypothesis that the ground state in these channels is a bound state. If the plateaus observed for the lowest-lying state are to be identified as “mirages” (that is, if the systems do not exhibit bound ground states), such fake plateaus could only result from cancellations between the FV states above the two-baryon threshold that contribute to the correlation function with opposites signs, and whose contributions depend upon the source structures. The spectrum of these states changes rapidly with power-law scaling as the volume is increased, in contrast with an exponential scaling for a compact state. In order for the “mirage plateaus” to be nearly coincident over the large range of volumes considered here, – 300 fm, the linear combination of states would also have to change very rapidly, and in a finely-tuned manner, in order to keep the plateau regions approximately volume independent. As the employed sources are volume independent and compact on the scale of all the spatial volumes, such behavior is exceedingly unlikely. There is no indication that the values of the ground-state energies in the two-baryon correlation functions scale as a power-law with the volume, and any multi-level model of these correlation functions with the exclusion of one or more possible bound state(s) fails to reproduce the behavior shown in Fig. 7.

In summary, the “mirage plateau” issue posed by Iritani et al. Iritani et al. (2017) appears to be irrelevant to the calculations presented here. The results of the present work are consistent with the correlation functions in each of two-baryon channels relaxing into a bound state at late times, with their binding energies determined in the next section. Iritani et al. additionally question the validity of the scattering amplitudes arising from these spectral studies, but again these claims have no bearing on the current results as is shown in Appendix D (see also Ref. Beane et al. (2017), where a coherent rebuttal of Ref. Iritani et al. (2017) is presented).

### iii.3 Results and discussions

In this section, the results for the LQCD spectra will be used to: *1)* obtain the -wave^{6}*2)* constrain the ERE parametrization of the scattering amplitudes, *3)* determine bound states and their binding energies, *4)* examine the naturalness of -wave baryon-baryon interactions, and *5)* provide constraints on the leading -symmetric interactions and well as the leading -symmetric interactions in the limit of large .

#### function

Given the ten FV energy eigenvalues determined in the previous section for each two-baryon channel, each scattering amplitude can be constrained at ten kinematic points via Lüscher’s QCs, Eqs. (2) and (5). In the NR limit, the CM energy eigenvalues corresponding to a two-baryon system at rest must be identical to that of the system in motion with two units of momentum in one Cartesian direction (the direction of total spin in a spin-triplet system) Briceno et al. (2013). Therefore, two sets of energy eigenvalues obtained from and measurements on the same ensemble do not provide constraints on scattering amplitude at distinct kinematic points. Nonetheless, given that these are obtained from separate sets of measurements (they are different Fourier projections of correlation functions with the same interpolating operators), including both sets in the analysis leads to better constraints on the scattering parameters and the binding energies.

The -wave scattering amplitude of the two-baryon channels with and belonging to the irrep, e.g. , is parametrized by a single phase shift, whose value can be constrained at a given CM momentum using the QC in Eq. (2), up to contaminations from -wave interactions that are neglected. The resulting function is shown in Fig. 8 for the ten energy eigenvalues obtained in the previous section. The figure includes the corresponding functions from which is obtained, see Eqs. (2) and (3) with . The function, whose intersection with determines the location of the bound state pole in the amplitude (see Eq. (1)), is also shown in Fig. 8.

For the spin-triplet channels , and associated with the , and irreps, respectively, additional mixing into the -wave in anticipated. As discussed in Sec. II.1, in the Blatt-Biedenharn parametrization, and with the boost vectors and , the -wave phase shift can be constrained at a given CM momentum using the QC in Eq. (5), up to negligible contaminations from -wave interactions. The resulting functions are plotted in Fig. 8 for the ten energy eigenvalues obtained in the previous section in each of these channels.

#### Effective range expansion parameters

Below the start of the t-channel cut, the function for the -wave (-wave) amplitude is anticipated to be well described by an ERE, see Eq. (7). Assuming that the pion is the lightest hadron exchanged between the baryons at this value of the quark masses, the t-channel cut starts at , considerably higher than the values obtained from the FV spectra in all channels. The constrained values of as a function of can thus be fit by two and three-parameter forms in each of the two-baryon channels, and the resulting fit bands are shown in Figs. 9-12. The values at the ten kinematic points considered here are also shown. Note that the vertical and horizontal error bars are displayed for simplicity and do not reflect the strongly correlated distributions of the and results. The precise form of the uncertainties are those shown in Fig. 8. In all channels, a two-parameter ERE describes the data well. The three-parameter ERE fits provide only small improvements in the values of of the fits, with the resulting scattering lengths and effective ranges being consistent with those of the two-parameter fit but with larger uncertainties. The values of the inverse scattering lengths and effective ranges from the two and three-parameter fits, as well as the shape parameters from the three-parameter fits, are listed in Table 2. The fit parameters are correlated, with their best values described by a multi-dimensional confidence ellipsoid. The and confidence ellipses from the two-parameter ERE are shown in Fig. 13, with the values of the center of the ellipses, their semi-minor and semi-major axes, as well as the slope of the semi-major axis of each ellipse listed in Table 14 of Appendix C.

The values of the inverse scattering lengths and effective ranges of the two-parameter EREs that are tabulated in Table 2 in lattice units can be expressed in physical units:

(19) | |||

(20) | |||

(21) | |||

(22) |

The numbers in the first and second parentheses denote, respectively, the statistical uncertainty, and the systematic uncertainty propagated from the corresponding uncertainties in the energies. The uncertainty in the lattice spacing is small compared with other uncertainties. Although these calculations have been performed for heavy quark masses at the flavor-symmetric point and without QED interactions, it is still interesting to compare these parameters with those in nature. While constraints on hyperon-nucleon and hyperon-hyperon scattering are not precise enough for a useful comparison, there exist precise determinations of nucleon-nucleon scattering parameters at low energies. In particular, the experimental values of the and -wave scattering lengths and effective ranges are

(23) | |||||

(24) |

which should be compared with the scattering parameters in the and irreps above, respectively. It is observed that the ranges of interactions in both channels, as characterized by the effective range parameters, are larger in nature than they are in the present work with heavier quark masses. Furthermore, the scattering lengths are smaller in this calculation than they are in nature. A more in-depth discussion of these parameters, in particular with regard to the unnaturalness of interactions and their spin-flavor symmetries, will be presented in Secs. III.3.4 and III.3.6.

#### Bound states and binding energies

There is clear evidence for the existence of a bound state in each of the channels belonging to the , and irreps of the decomposition of the product of two octet baryons. First, as the volume is increased, the ground-state energies of two-baryon systems converge to a negatively-shifted energy far away from the two-particle threshold. Second, the analytic continuation of the amplitudes in these channels is consistent with the presence of a pole in the amplitude for , as is evident from the intersection of the ERE bands with the function in Figs. 9, 10 and 12. The location of this intersection determines the square of the binding momentum in these channels. However, given the uncertainties in the ERE fits, the infinite-volume extrapolation of negatively-shifted energies leads to a more precise determination of the binding energies. With energies determined in three volumes, a controlled extrapolation to infinite volume is possible in the present work. Fitting to the truncated form of the FV QC for negative values, Eq. (8), the infinite-volume binding momenta, , can be obtained in each channel. These results are presented in Table 3 for measurements with and , with complete agreement seen between the two determinations. The bootstrap samples of extracted values from each case can be combined to obtain a conservative estimate of the binding momenta and their uncertainties, given in the last row of Table 3. The omitted terms in the truncated form in Eq. (8) are negligible as is at most for the channels belonging to the , and irreps. The stability of the extracted binding momenta has been verified by excluding lower-order terms and by adding higher-order terms to the fits.

Table 3 also includes the values for the channels belonging to the irrep. As is seen from Fig. 11, the ground-state energy in the largest volume is close to threshold. Nonetheless, assuming that there is a bound state in this channel, a determination of based on the fit to Eq. (8) is fully consistent with the ground-state energies at the largest volume, as well as with the location of the pole in the scattering amplitude. From these results, the existence of a bound state in the irrep cannot be confirmed or excluded with statistical significance. Future calculations with higher statistics are needed in order to draw robust conclusions about the nature of the ground state in the irrep.

In physical units, the binding energies of these states are:

(25) | |||||

(26) | |||||

(27) | |||||

(28) |

where . Again, the first uncertainty is statistical and the second uncertainty encompasses both a fitting uncertainty and an uncertainty encoding variation among multiple analyses. The uncertainty in the lattice spacing is small compared with other uncertainties. These binding energies are consistent with our previous determination in Ref. Beane et al. (2013a, b), and with the binding energies obtained on the same ensembles of gauge-field configurations in Ref. Berkowitz et al. (2017) for the ground states of the two-nucleon channels in the and irreps.

#### -wave baryon-baryon interactions and naturalness

Interactions are considered unnatural if they give rise to some characteristic length scale of the system that is much larger than their range. There are at least two measures to assess naturalness in a two-particle system. For scattering states at low energies, scattering length defines a characteristic length scale, and the range of interactions can be approximated by the effective range. As an example, -wave interactions in the spin-singlet and spin-triplet two-nucleon channels in nature produce effective range to scattering length ratios, , that are and , respectively. This indicates that both channels are unnatural, particularly the spin-singlet channel. When interactions support a bound state, another characteristic length scale of the two-particle system is the inverse of the binding momentum, which defines an intrinsic size for the bound state. Considering the exchange of the pion to be the dominant contribution to the long-range part of effective interactions among two nucleons, the ratio of the binding momentum to the pion mass provides another measure of unnaturalness of interactions. In nature, and for the di-neutron and deuteron, respectively, again indication that both channels are unnatural.^{7}

The ratios of the scattering lengths to effective ranges obtained from both the two and three-parameter ERE fits are shown in Table 4 for the , , and two-baryon channels. Interestingly, these ratios are universally consistent with within uncertainties, a feature that points to a spin-flavor symmetry of interactions as will be discussed in Sec. III.3.6 (see also Ref. Beane et al. (2013b)). This value indicates that all channels are governed by -wave interactions that are only slightly less unnatural than those in the spin-triplet two-nucleon system in nature. This also implies that the -wave interactions in a spin-singlet two-nucleon state undergo a more dramatic change as a function of the quark masses and appear more finely tuned in nature, a feature that was also pointed to in Ref. Beane et al. (2013b).

The ratios of the binding momenta to the pion mass of this calculation () for the two-baryon channels in the , , and irreps are generally close to , similar to the value of for the deuteron in nature. However, a comparison between the effective range in each channel and the inverse pion mass of this calculation suggest that pion exchange may not be the dominant contribution to the long-range forces between the baryons Beane et al. (2013b). In particular, the intrinsic size of the bound state in each channel (set by the inverse binding momenta), is comparable to the corresponding effective range in each channel (with larger uncertainties in the irrep), a feature that is consistent with the results obtained for the ratio of the effective ranges to the scattering lengths, . Nonetheless, while the bound states in the , and irreps appear to have a natural intrinsic size, the scattering lengths in all channels are still large compared with the binding momenta of the bound states. Consequently, an EFT treatment of these channels at low energies with the KSW-vK power-counting scheme is justified.

#### Large- limit and leading interactions in effective field theory

The scattering parameters in two-baryon channels belonging to different irreps have similar values in the present calculations, differing by at most from an average value. This feature is illustrated in Fig. 14, in which the scattering lengths and effective ranges from the two and three-parameter ERE fits, and the shape parameters from the three-parameter ERE fits, are compared. Additionally, the correlated ratios from both ERE fits are shown for the four irreps. All of these quantities are broadly consistent between different channels, in particular between channels belonging to the and irreps. This is a manifestation of an approximate spin-flavor symmetry of nuclear and hypernuclear forces as predicted to exist in the large- limit of QCD Kaplan and Savage (1996), as introduced in Sec. II.2. Deviations from symmetry (involving the and irreps) are expected to scale as while the deviations from symmetry are expected to scale as . This is consistent with the almost identical nature of the channels belonging to the and irreps in the results presented here. Additional deviations from the spin-flavor symmetry that occur as a result of the flavor-symmetry breaking in nature are absent in the present calculations, making them an ideal testing ground for the large- relations.

Given an approximate symmetry of -wave interactions, constraints can be obtained on the coefficients, and , using Eqs. (14). Two cases are considered here: unnatural interactions and natural interactions. In the unnatural case, the leading -wave interactions are summed to all orders in perturbation theory (implementing KSW-vK power counting) introducing a UV-scale dependence in the coefficients. A convenient choice of renormalization scale is . However, any scale much above the largest inverse scattering length, but below the cut off of the theory, would lead to manifest power counting and RG-scale independence. The values of and obtained from each pair of equations in (14) are tabulated in Table 5, and are shown in Fig. 15. Within the uncertainty of each determination, these values are in agreement with each other. Note that from Eq. (14), contributions from the coefficient are suppressed by at least a factor of 3 compared with those from the coefficient, thus the rescaled coefficient is considered. A combined fit of a constant to the five determinations results in the values for and that are listed in Table 5 and shown as pink bands in Fig. 15.

Assuming the systems to be natural in the EFT analysis results in large uncertainties in the coefficient. This precludes conclusions to be drawn regarding its significance compared to the coefficient. Additionally, determinations that involve the irrep yield large uncertainties in the coefficients, signaling the inappropriate assumption of naturalness for interactions in a channel that is almost at unitarity within uncertainties. As the observations in Sec. III.3.4 point to primarily an unnatural scenario for all interactions considered, the values of the and coefficients in the unnatural scenario, defined with KSW-vK power counting, are found to be more relevant. In particular, it is observed that the value of is approximately an order of magnitude smaller than the value of . This is a signature of an accidental symmetry of nuclear and hypernuclear forces which was first predicted in Ref. Kaplan and Savage (1996), and is studied here directly with QCD for the first time.

Without constraints on the scattering parameters belonging to the and irreps, a conclusive statement regarding the spin-flavor symmetry in the interactions is not possible. However, with the observations in other irreps pointing to such a symmetry, a prediction can be made for the scattering amplitudes in the and irreps, giving , where the second uncertainty accounts for corrections to the prediction of symmetry. This result enables an extraction of all SW coefficients of the LO -symmetric interactions.