Low-energy Scattering and Effective Interactions of Two Baryons at $m_{\pi}\sim 450$ MeV from Lattice Quantum Chromodynamics

The interactions between two octet baryons are studied at low energies using lattice QCD (LQCD) with larger-than-physical quark masses corresponding to a pion mass of $m_{\pi}\sim 450$ MeV and a kaon mass of $m_{K}\sim 596$ MeV. The two-baryon systems that are analyzed range from strangeness $S=0$ to $S=-4$ and include the spin-singlet and triplet $NN$, $\Sigma N$ ($I=3/2$), and $\Xi\Xi$ states, the spin-singlet $\Sigma\Sigma$ ($I=2$) and $\Xi\Sigma$ ($I=3/2$) states, and the spin-triplet $\Xi N$ ($I=0$) state. The $s$-wave scattering phase shifts, low-energy scattering parameters, and binding energies when applicable, are extracted using L\"uscher's formalism. While the results are consistent with most of the systems being bound at this pion mass, the interactions in the spin-triplet $\Sigma N$ and $\Xi\Xi$ channels are found to be repulsive and do not support bound states. Using results from previous studies at a larger pion mass, an extrapolation of the binding energies to the physical point is performed and is compared with experimental values and phenomenological predictions. The low-energy coefficients in pionless EFT relevant for two-baryon interactions, including those responsible for $SU(3)$ flavor-symmetry breaking, are constrained. The $SU(3)$ symmetry is observed to hold approximately at the chosen values of the quark masses, as well as the $SU(6)$ spin-flavor symmetry, predicted at large $N_c$. A remnant of an accidental $SU(16)$ symmetry found previously at a larger pion mass is further observed. The $SU(6)$-symmetric EFT constrained by these LQCD calculations is used to make predictions for two-baryon systems for which the low-energy scattering parameters could not be determined with LQCD directly in this study, and to constrain the coefficients of all leading $SU(3)$ flavor-symmetric interactions, demonstrating the predictive power of two-baryon EFTs matched to LQCD.


2
lower maximum values for neutron star masses. While experimental data on scattering cross sections in the majority of the Y N channels are scarce, there are reasonably precise constraints on the interactions in the ΛN channel from scattering and hypernuclear spectroscopy experiments [2,3], and they indicate that the interactions in this channel are attractive. Given that the Λ baryon is lighter than the other hyperons, it is likely the most abundant hyperon in the interior of neutron stars. However, models of the EoS including Λ baryons and attractive ΛN interactions [4] predict a maximum neutron star mass that is below the maximum observed mass at 2M [5][6][7][8][9]. 1 Several remedies have been suggested to solve this problem, known in the literature as the "hyperon puzzle" [11][12][13]. For example, if hyperons other than the Λ baryon (such as Σ baryons) are present in the interior of neutron stars and the interactions in the corresponding Y N and Y Y channels are sufficiently repulsive, the EoS would become more stiff [14,15]. Another suggestion is that repulsive interactions in the Y N N , Y Y N and Y Y Y channels may render the EoS stiff enough to produce a 2M neutron star [4,14,[16][17][18]. Repulsive density-dependent interactions in systems involving the Λ and other hyperons have also been suggested, along with the possibility of a phase transition to quark matter in the interior of neutron stars, see Refs. [11][12][13] for recent reviews. Given the scarcity or complete lack of experimental data on Y N and Y Y scattering and all three-body interactions involving hyperons, SU (3) flavor symmetry is used to constrain effective field theories (EFTs) and phenomenological meson-exchange models of hypernuclear interactions. In this way, quantities in channels for which experimental data exist can be related via symmetries to those in channels which lack such phenomenological constraints [19,20]. For example, the lowest-order effective interactions in several channels with strangeness S ∈ {−2, −3, −4} were constrained using experimental data on pp phase shifts and the Σ + p cross section in the same SU (3) representation in the framework of chiral EFT (χEFT) in Refs. [21][22][23]. However, only a few of the SU (3)-breaking low-energy coefficients (LECs) of the EFT could be constrained [23]. 2 To date, the knowledge of these interactions in nature remains unsatisfactory, demanding more direct theoretical approaches. 3 Building upon our previous works [24][25][26][27][28][29][30][31], we present further studies which constrain hypernuclear forces in nature by direct calculations starting from the underlying theory of the strong interactions, Quantum Chromodynamics (QCD). To this end, the numerical technique of lattice QCD (LQCD) is used to obtain information on the low-energy spectra and scattering on twobaryon systems, which can be used to constrain EFTs or phenomenological models of two-baryon interactions. In recent years, LQCD has allowed a wealth of observables in nuclear physics, from hadronic spectra and structure to nuclear matrix elements [32][33][34], to be calculated directly from interactions of quarks and gluons, albeit with uncertainties that are yet to be fully controlled. In the context of constraining hypernuclear interactions, LQCD is a powerful theoretical tool because the lowest-lying hyperons are stable when only strong interactions are included in the computation, circumventing the limitations faced by experiments on hyperons and hypernuclei. Nonetheless, LQCD studies in the multi-baryon sector require large computing resources as there is an inherent signal-to-noise degradation present in the correlation functions of baryons [25,26,[35][36][37][38], among other issues as discussed in a recent review [34]. Consequently, most studies of two-baryon systems to date [24,[26][27][28][29][30][31][39][40][41][42][43][44][45][46][47][48][49][50][51] have used larger-than-physical quark masses to expedite computations, and only recently have results at the physical values of the quark masses emerged [52][53][54][55], making it possible to directly compare with experimental data [56]. The existing studies are primarily based on two distinct approaches. In one approach, the low-lying spectra of two baryons in finite spatial volumes are determined from the time dependence of Euclidean correlation functions computed with LQCD, and are then converted to scattering amplitudes at the corresponding energies through the 3 use of Lüscher's formula [57,58] or its generalizations [59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75]. In another approach, non-local potentials are constructed based on the Bethe-Salpeter wavefunctions determined from LQCD correlation functions, and are subsequently used in the Lippmann-Schwinger equation to solve for scattering phase shifts [55,76,77]. Given that Lüscher's formalism is model-independent below inelastic thresholds, it is this approach that is used in the present study as the basis to constrain scattering amplitudes and their low-energy parametrizations in a number of two-(octet)baryon channels with strangeness S ∈ {0, −1, −2, −3, −4}.
While LQCD studies at unphysical values of the quark masses already shed light on the understanding of (hyper)nuclear and dense-matter physics, a full account of all systematic uncertainties, including precise extrapolations to the physical quark mass, is required to further impact phenomenology. Additionally, LQCD results for scattering amplitudes can be used to better constrain the low-energy interactions within given phenomenological models and applicable EFTs. In the case of exact SU (3) flavor symmetry and including only the lowest-lying octet baryons, there are six two-baryon interactions at leading order (LO) in pionless EFT [78,79] that can be constrained by the s-wave scattering lengths in two-baryon scattering [80]. LQCD has been used in Ref. [31] to constrain the corresponding LECs of these interactions by computing the s-wave scattering parameters of two baryons at an SU (3) flavor-symmetric point with m π ∼ 806 MeV. Strikingly, the first evidence of a long-predicted SU (6) spin-flavor symmetry in nuclear and hypernuclear interactions in the limit of a large number of colors (N c ) [81] was observed in that study, along with an accidental SU (16) symmetry. This extended symmetry has been suggested in Ref. [82] to support the conjecture of entanglement suppression in nuclear and hypernuclear forces at low energies, pointing to intriguing aspects of strong interactions in nature.
The objective of this paper is to extend our previous study to quark masses that are closer to their physical values, corresponding to a pion mass of ∼ 450 MeV and a kaon mass of ∼ 596 MeV, and further to study these systems in a setting with broken SU (3) flavor symmetry as is the case in nature. The present study provides new constraints that allow preliminary extrapolations to physical quark masses to be performed, and complements previous independent LQCD studies at nearby quark masses [24, 26-29, 39, 42, 45, 46, 83, 84]. In particular, predictions for the binding energies of ground states in a number of Y N and Y Y channels based on the results of the current work and those of Ref. [31] at larger quark masses are consistent with experiments and phenomenological results where they exist. Our LQCD results are used to constrain the leading SU (3) symmetry-breaking coefficients in pionless EFT. This EFT matching enables the exploration of large-N c predictions, pointing to the validity of SU (6) spin-flavor symmetry at this pion mass as well, and revealing a remnant of an accidental SU (16) symmetry that was observed at a larger pion mass in Ref. [31]. Strategies to make use of the QCD-constrained EFTs to advance the ab initio many-body studies of larger hypernuclear isotopes and dense nuclear matter are beyond the scope of this work. Nevertheless, the methods applied in Refs. [85][86][87] to connect the results of LQCD calculations to higher-mass nuclei can also be applied in the hypernuclear sector using the results presented in this work. This paper is organized as follows. Sec. II presents a summary of the computational details (Sec. II A), followed by the results for the lowest-lying energies of two-baryon systems from LQCD correlation functions, along with a description of the method used to obtain these spectra (Sec. II B), a determination of the s-wave scattering parameters in the two-baryon channels that are studied, along with the formalism used to extract the scattering amplitude (Sec. II C), and finally the binding energies of the bound states that are identified in various channels, including an extrapolation to the physical point (Sec. II D). Sec. III discusses the constraints that these results impose on the low-energy coefficients of the next-to-leading order (NLO) pionless EFT Lagrangian, including neutron stars, can be used to indirectly probe the strangeness content of dense matter, and provide complementary constraints on models of hypernuclear interactions [11]. TABLE I. Parameters of the gauge-field ensembles used in this work. L and T are the spatial and temporal dimensions of the hypercubic lattice, β is related to the strong coupling, b is the lattice spacing, m l(s) is the bare light (strange) quark mass, N cfg is the number of configurations used and N src is the total number of sources computed. For more details, see Ref. [42]. only SP correlation functions are produced for computational expediency. Table I summarizes the parameters of these ensembles. Correlation functions are constructed by forming baryon blocks at the sink [39]: where S (f ),n n is a quark propagator with flavor f ∈ {u, d, s} and with combined spin-color indices n, n ∈ {1, . . . , N s N c }, where N s = 4 is the number of spin components and N c = 3 is the number of colors. The weights w B i j k are tensors that antisymmetrize and collect the terms needed to have the quantum numbers of the baryons B ∈ {N, Λ, Σ, Ξ}. The sum over the sink position x projects the baryon blocks to well-defined three-momentum p. In particular, two-baryon correlation functions were generated with total momentum P = p 1 + p 2 , where p i is the three-momentum of the i-th baryon taking the values p i = 2π L n with n ∈ {(0, 0, 0), (0, 0, ±1)}. Therefore, P = 2π L d, with d ∈ {(0, 0, 0), (0, 0, ±2)}. 4 Additionally, two baryon correlation functions with back-toback momenta were generated at the sink, with momenta p 1 = −p 2 = 2π L n. This latter choice provides interpolating operators for the two-baryon system that primarily overlap with states that are unbound in the infinite-volume limit, providing a convenient means to identify excited states as well. The construction of the correlation functions continues by forming a fully-antisymmetrized quark-level wavefunction at the location of the source, with quantum numbers of the two-baryon system of interest. Appropriate indices from the baryon blocks at the sink are then contracted with those at the source, in a way that is dictated by the quark-level wavefunction, see Refs. [30,102] for more detail. The contraction codes used to produce the correlation functions in this study are the same as those used to perform the contractions for the larger class of interpolating operators used in our previous studies of the SU (3) flavor-symmetric spectra of nuclei and hypernuclei up to A = 5 [30], and two-baryon scattering [31,41,42]. 5 .
In this study, correlation functions for nine different two-baryon systems have been computed, ranging from strangeness S = 0 to −4. Using the notation ( 2s+1 L J , I), where s is the total spin, L is the orbital momentum, J is the total angular momentum, and I is the isospin, the systems are: Under strong interactions, these channels do not mix with other two-baryon channels or other hadronic states below three-particle inelastic thresholds. In the limit of exact SU (3) flavor symmetry, the states belong to irreducible representations (irreps) of SU (3): 27 (all the singlet states), 10 (triplet N N ), 10 (triplet ΣN and ΞΞ), and 8 a (triplet ΞN ). In the rest of this work, the isospin label will be dropped for simplicity.

B. Low-lying finite-volume spectra of two baryons
The two-point correlation functions constructed in the previous section have spectral representations in Euclidean spacetime. Explicitly, the correlation function CÔ ,Ô (τ ; d) formed using the 6 source (sink) interpolating operatorsÔ (Ô ) can be written as: where all quantities are expressed in lattice units. E (i) is the energy of the ith eigenstate |E (i) , , and V = L 3 . The boost-vector dependence of the energies, states, and overlap factors is implicit. The lowest-lying energies of the one-and two-baryon systems required for the subsequent analyses can be extracted by fitting the correlation functions to this form. To reliably discern the first few exponents given the discrete τ values and the finite statistical precision of the computations is a challenging task. In particular, a well-known problem in the study of baryons with LQCD is the exponential degradation of the signal-to-noise ratio in the correlation function as the source-sink separation time increases-an issue that worsens as the masses of the light quarks approach their physical values. First highlighted by Parisi [35] and Lepage [36], and studied in detail for light nuclei in Refs. [25,26], it was later shown that this problem is related to the behavior of the complex phase of the correlation functions [38,104]. Another problem that complicates the study of multi-baryon systems is the small excitation gaps in the finite-volume spectrum that lead to significant excited-state contributions to correlation functions. To overcome these issues, sophisticated methods have been developed to analyze the correlation functions, such as Matrix Prony [37] and the generalized pencilof-function [105] techniques, as well as signal-to-noise optimization techniques [106]. Ultimately, a large set of single-and multi-baryon interpolating operators with the desired quantum numbers must be constructed to provide a reliable variational basis to isolate the lowest-lying energy eigenvalues via solving a generalized eigenvalue problem [107], as is done in the mesonic sector [75]. Such an approach is not yet widely applied to the study of two-baryon correlation functions, given its computational-resource requirement, but progress is being made. In Ref. [50], a partial set of two-baryon scattering interpolating operators were used to study the two-nucleon and H-dibaryon channels with results that generally disagreed with previous works [27,30,49]. Investigations continue to understand and resolve the observed discrepancies [31,33,34,88,[108][109][110]. For the present study, in which only up to two types of interpolating operators were computed, no variational analysis could be performed. Instead, we have developed a robust automated fitting methodology to sample and combine fit range and model selection choices for uncertainty quantification.
Given that the correlation functions are only evaluated at a finite number of times and with finite precision, to fit Eq. (2) the spectral representation is truncated to a relatively small number of exponentials and fitted in a time range {τ min , τ max }, where τ max is set by a threshold value determined by examining the signal-to-noise ratio, and τ min is chosen to take values in the interval [2, τ max − τ plateau ]. Here, τ plateau = 5 is chosen to be the minimum length of the fitting window (numbers are expressed in units of the lattice spacing). A scan over all possible fitting windows is treated as a means to quantify the associated systematic uncertainty. With a fixed window, a correlated χ 2 -function is minimized to obtain the fit parameters Z i Z * i and E i for i ∈ {0, 1, . . . , e}, where e+1 is the number of exponentials in a given fit form. Variable projection techniques [111,112] are used to obtain the value of the overlap factors for a given energy, since they appear linearly in CÔ ,Ô (τ ; d). Furthermore, given the finite statistical sampling of correlation functions, shrinkage techniques [113] are used to better estimate the covariance matrix. The number of excited states included in the fit is decided via the Akaike information criterion [114]. The confidence intervals of the parameters are estimated via the bootstrap resampling method. For a fit to be included in the set of accepted fits (later used to extract the fit parameters and assess the resulting uncertainties), several checks must be passed, including χ 2 /N dof being smaller than 2, and different optimization algorithms leading to consistent results for the parameters (within a tolerance)-see Ref. [115] for further details on this and other checks. The accepted fits are then combined to give the final result for the mean value of the energy, with weights ω f that are chosen to be the following combination of the p-value, p f , and the uncertainty of each fit, δE f : see Ref. [116] for a Bayesian framework. Here, the indices f, f run over all the accepted fits. The statistical uncertainty is defined as that of the fit with the highest weight, while the systematic uncertainty is defined as the average difference between the weighted mean value and each of the accepted fits: It should be noted that instead of fitting to the correlation function, the effective energy function can be employed, derived from the logarithm of the ratio of correlation functions at displaced times, where τ J is a non-zero integer that is introduced to improve the extraction of E (0) (for a detailed study, see Ref. [37]). Consistent results are obtained when either correlation functions or the effective energy functions are used as input.
In order to identify the shift in the finite-volume energies of two baryons compared with noninteracting baryons, the following ratio of two-baryon and single-baryon correlation functions can be formed with an associated effective energy-shift function, Constant fits to ratios of correlation functions can be used to obtain energy shifts ∆E (0) = E (0) − m 1 −m 2 (where m 1 and m 2 are the masses of baryon B 1 and B 2 , respectively), requiring time ranges such that both the two-baryon and the single-baryon correlation functions are described by a singlestate fit. However, if both correlation functions are not in their ground states, cancellations may occur between excited states (including the finite-volume states that would correspond to elastic scattering states in the infinite volume), either in correlation function or in ratios of correlation functions, producing a "mirage plateau" [88]. Despite this issue, as demonstrated in Ref. [108], our previous results, such as those in Refs. [30,31] are argued to be free of this potential issue (similar discussions can be found in Ref. [110]). To determine ∆E (0) in this work, the two-baryon and singlebaryon correlation functions are fit to multi-exponential forms (which account for excited states) within the same fitting range, and afterwards the energy shifts are computed at the bootstrap level, in such a way that the correlations between the different correlation functions are taken into account. The use of correlated differences of multi-state fit results is convenient in particular for automated fit range sampling, since the number of excited states can be varied independently for one-and two-baryon correlation functions, unlike fits to the ratio in Eq. (7). Consistent results were obtained via fitting the ratio in Eq. (7) in the allowed time regions. The effective mass plots (EMPs) for the single-baryon correlation functions, and for each of the ensembles studied in the present work, are displayed in Fig. 22 of Appendix D. The bands shown in the figures indicate the baryon mass which results from the fitting strategy explained above, with the statistical and systematic uncertainties included, and the corresponding numerical values listed in Table II. The table also shows the baryon masses extrapolated to infinite volume, obtained by fitting the masses in the three different volumes to the following form: where M (∞) B and c B are the two fit parameters. This form incorporates leading-order (LO) volume corrections to the baryon masses in heavy-baryon chiral perturbation theory (HBχPT) [117]. As evident from the m π L values listed in Table I, the volumes used are large enough to ensure small volume dependence in the single-baryon masses. 6 This is supported by the observation that the value M While the Λ baryon is not relevant to subsequent analysis of the two-baryon systems studied in this work, the centroid of the four octet-baryon masses is used to define appropriate units for the EFT LECs, hence M Λ is reported for completeness.
The results for the two-baryon energy shifts are shown in Fig. 1. 7 For display purposes, the effective energy-shift functions, defined in Eq. (8), are shown in Figs. 23-31 of Appendix D, along with the corresponding two-baryon effective-energy functions, defined in Eq. (6). The associated numerical values are listed in Tables XVII-XXV of the same appendix. In each subfigure of Figs. 23-31, two correlation functions are displayed: the one yielding the lowest energy (labeled as n = 1 in Tables XVII-XXV) corresponds to having both baryons at rest or, if boosted, with the same value of the momentum, and the one yielding a higher energy (labeled as n = 2 in the tables) corresponds to the two baryons having different momenta, e.g., having back-to-back momenta or one baryon at rest and the other with non-zero momentum. While the first case (n = 1) couples primarily to the 9 ground state, the latter (n = 2) is found to have small overlap onto the ground state, and gives access to the first excited state directly.
As a final remark, it should be noted that the single-baryon masses and the energies extracted for the two-nucleon states within the present analysis are consistent within 1σ with the results of Ref. [42], obtained with the same set of data but using different fitting strategies. Despite this overall consistency, the uncertainties of the two-nucleon energies in the present work are generally larger compared with those reported in Ref. [42] for the channels where results are available in that work. The reason lies in a slightly more conservative systematic uncertainty analysis employed here. The comparison between the results of this work and that of Ref. [42] is discussed extensively in Appendix B.
C. Low-energy scattering phase shifts and effective-range parameters Below three-particle inelastic thresholds, Lüscher's quantization condition [57,58] provides a means to extract the infinite-volume two-baryon scattering amplitudes from the energy eigenvalues of two-baryon systems obtained from LQCD calculations, e.g., those presented in Sec. II B. This condition holds if the range of interactions is smaller than (half of) the spatial extent of the cubic volume, L, and the corrections to this condition scale as e −mπL for the two-baryon systems. Such corrections are expected to be small in the present work given the m π L values in Table I. The quantization conditions are those used in Refs. [31,42]: in the case of spin-singlet sates, only the s-wave limit of the full quantization condition is considered. For coupled 3 S 1 − 3 D 1 states, in which the Blatt-Biedenharn parametrization [118] of the scattering matrix can be used, only the α-wave approximation of the quantization condition is considered [74]. In both cases, and denoting the (s-wave or α-wave) phase shift by δ, the condition can be written as [63]: where k * is the center-of-mass (c.m.) relative momentum of each baryon, d is the total c.m. momentum in units of 2π/L, and c d lm is a kinematic function related to Lüscher's Z-function, Z d lm : with γ = E/E * being the relativistic gamma factor. Here, E and E * are the energies of the system in the laboratory and c.m. frames, respectively. The three-dimensional zeta-function is defined as where r =γ −1 (n − αd) and α = 1 2 1 + (m 2 1 − m 2 2 )/E * 2 , 8 with m 1 and m 2 being the masses of the two baryons. The factorγ −1 acting on a vector u modifies the parallel component with respect to d, while leaving the perpendicular component invariant, i.e.,γ −1 u = γ −1 u + u ⊥ . Convenient expressions have been derived to exponentially accelerate the numerical evaluation of the function in Eq. (12) [64,[119][120][121], and the following expression is used in the present analysis: 8 Here α should not be confused with α-wave mentioned above. FIG. 1. Summary of the energy shifts extracted from LQCD correlation functions for all two-baryon systems studied in this work, together with the non-interacting energy shifts defined as ∆E = m 2 1 + |p 1 | 2 + m 2 2 + |p 2 | 2 − m 1 − m 2 , where |p 1 | 2 = |p 2 | 2 = 0 corresponds to systems that are at rest (continuous line), |p 1 | 2 = |p 2 | 2 = ( 2π L ) 2 corresponds to systems which are either boosted or are unboosted but have back-toback momenta (dashed line), and |p 1 | 2 = 0 and |p 2 | 2 = ( 4π L ) 2 corresponds to boosted systems where only one baryon has non-zero momentum (dashed-dotted line). The points with no boost have been shifted slightly to the left, and the ones with boosts have been shifted to the right for clarity. Quantities are expressed in lattice units.
The values of k * cot δ at given k * 2 values are shown for all two-baryon systems in Fig. 2, and the associated numerical values are listed in Tables XVII-XXV of Appendix D. The validity of Lüscher's quantization condition must be verified in each channel, in particular in those that have exhibited anomalously large ranges, such as ΣN ( 3 S 1 ), in previous calculations. The consistency between solutions to Lüscher's condition and the finite-volume Hamiltonian eigenvalue equation using a LO EFT potential was established in Ref. [29] for the ΣN ( 3 S 1 ) channel and at values of the quark masses (m π ∼ 389 MeV) close to those of the current analysis. The conclusion of Ref. [29], therefore, justifies the use of Lüscher's quantization condition in the current work for this channel. The energy dependence of k * cot δ can be parametrized by an effective range expansion (ERE) below the t-channel cut [122][123][124] where a is the scattering length, r is the effective range, and P is the leading shape parameter. These parameters can be constrained by fitting k * cot δ values obtained from the use of Lüscher's quantization condition as a function of k * 2 . To this end, one could use a one-dimensional choice of the χ 2 function, minimizing the vertical distance between the fitted point and the function, where 10 f (a −1 , r, P, k * 2 ) corresponds to the ERE parametrization given by the right-hand side of Eq. (14), and the sum runs over all extracted pairs of {k * 2 i , (k * cot δ) i }, where the compound index i counts data points for different boosts, n values of the level, and different volumes. Each contribution is weighted by an effective variance that results from the combination of the uncertainty in both k * 2 with δx being the mid-68% confidence interval of the quantity x. The uncertainty on the {k * 2 i , (k * cot δ) i } pair can be understood by recalling that each pair is a member of a bootstrap ensemble with the distribution obtained in the previous step of the analysis. To generate the distribution of the scattering parameters, pairs of {k * 2 i , (k * cot δ) i } are randomly selected from each bootstrap ensemble and are used in Eq. (15) to obtain a new set of {a −1 , r, P } parameters. This procedure is repeated N times, where N is chosen to be equal to the number of bootstrap ensembles for {k * 2 i , (k * cot δ) i }. This produces an ensemble of N values of fit parameters {a −1 , r, P }, from which the central value and the associated uncertainty in the parameters can be determined (median and mid-68% intervals are used for this purpose).
Alternatively, one can use a two-dimensional choice of the χ 2 function. 11 Knowing that k * cot δ values must lie along the Z-function, as can be seen from Eq. (10) and Fig. 2, one could take the distance between the data point and the point where the ERE crosses the Z-function along this function (arc length) in the definition of χ 2 . Explicitly, where σ 2 i is now defined as and D Z [{x 1 , y 1 }, {x 2 , y 2 }] denotes the distance between the two points {x 1 , y 1 } and {x 2 , y 2 } along the Z-function. The quantity K * 2 is the point where the ERE (f in Eq. (16)) crosses the Z-function.
To obtain this point, and given the large number of discontinuities present in the Z-function, Householder's third order method can be used as a reliable root-finding algorithm [125]: where the starting point is set to be K * 2 0 = k * 2 i , and the number of primes over 1/F indicates the order of the derivative computed at the point K * 2 m . The stopping criterion is defined as |K * 2 m+1 − K * 2 m | < 10 −6 , which occurs for m ∼ O (10). Since the extraction of this point requires knowledge of scattering parameters, the minimization must be implemented iteratively. This second choice of χ 2 function has been used in the main analysis of this work, however, the use of the one-dimensional  For a precise extraction of the ERE parameters, a sufficient number of points below the t-channel cut must be available, for positive or negative k * 2 . In general, for the channels studied throughout this work, there are only a few points in the positive k * 2 region below the t-channel cut (starting at k * 2 t-cut ∼ 0.018 l.u.). For a non-interacting system, states above scattering threshold have c.m. energies m 2 1 + k * 2 + m 2 2 + k * 2 , with the c.m. momenta roughly scaling with the volume as k * 2 ∼ (2π|n|/L) 2 . With the minimum value of |n| 2 used in this work (|n| 2 = 1), only the states from the ensemble with L = 48 are expected to lie below the t-channel cut (4π 2 /48 2 < k * 2 t-cut ). This behavior is consistent with the data. Comparing with the results of the analysis at m π ∼ 806 MeV in Ref. [31], where lattice configurations of comparable size (in lattice units) were used, the larger value of the pion mass resulted in the position of the t-channel cut being moved further away from zero, and the majority of the lowest-lying states extracted in that study remained inside the region where the ERE parametrization could be used. Therefore, with only ground-state energies available for the analysis of the ERE in the ensembles with L ∈ {24, 32}, the precision in the extraction of scattering parameters is noticeably reduced compared with the study at m π ∼ 806 MeV in Ref. [31]. Inclusion of the shape parameter, P , does not improve the fits, and although the scattering lengths remain consistent with those obtained with a two-parameter fit, the effective ranges are larger in magnitude, and the uncertainties in the scattering parameters are increased. Moreover, the central values of the extracted shape parameters are rather large, bringing into question the assumption that the contribution of each order in the ERE should be smaller than the previous order. However, uncertainties on the shape parameters are sufficiently large that no conclusive statement can be made regarding the convergence of EREs. In one case, i.e., the ΣN ( 1 S 0 ) channel, the threeparameter ERE fit is not performed given the large uncertainties. For these reasons, while the TABLE III. The values of the inverse scattering length a −1 , effective range r, and shape parameter P determined from the two-and three-parameter ERE fits to k * cot δ versus k * 2 for various two-baryon channels (see Fig. 3). Quantities are expressed in lattice units.  scattering parameters are reported for both the two-and three-parameter fits in this section, only those of the two-parameter fits will be used in the EFT study in the next section.
Fits to k * cot δ as a function of k * 2 in various two-baryon channels are shown in Fig. 3, along with the correlation between inverse scattering length and effective range in each channel depicted in Fig. 4 using the 68% and 95% confidence regions of the parameters. The areas in the parameter space that are prohibited by the constraints imposed by Eq. (10) are also shown in Fig. 4, highlighting the fact that the two-parameter ERE must cross the Z-functions for each volume in the negative-k * 2 region. For fits including higher-order parameters, these constraints are more complicated and are not shown. For the ΣN ( 3 S 1 ) and ΞΞ ( 3 S 1 ) channels, the ground-state energy is positively shifted, i.e., ∆E 0, and only the values of k * 2 associated with the ground states are inside the range of validity of the ERE. As a result, no extraction of the ERE parameters is possible in these channels given the number of data points. Results for the scattering parameters obtained using two-and three-parameter ERE fits in the other seven channels are summarized in Table III, and are shown in Fig. 5 for better comparison in the case of two-parameter fits.
The inverse scattering lengths extracted for all systems are compatible with each other (albeit within rather large uncertainties), providing a possible signature that SU (3) flavor symmetry is approximately respected in these two-baryon systems at this pion mass and at low energies, see Sec. III. The effective range in most systems is compatible with zero. 12 Furthermore, the ratio r/a can be used as an indicator of the naturalness of the interactions; for natural interactions, r/a ∼ 1, while for unnatural interactions r/a 1. At the physical point, both N N channels are unnatural and exhibit large scattering lengths, with r/a being close to 0.1 for the spin-singlet channel, and 0.3 for the spin-triplet channel. From Table IV, the most constrained ratios are obtained for the ΣΣ ( 1 S 0 ), ΞΞ ( 1 S 0 ), and ΞN ( 3 S 1 ) channels, for which r/a ∼ 0.2 − 0.3, indicating unnatural interactions at low energies. For other channels, the larger uncertainty in this ratio precludes drawing conclusions about naturalness. Alternatively, naturalness can be assessed by considering the ratio of the binding momentum to the pion mass, as this quantity is better constrained in this study, see Table VI in the next subsection. The values for κ (∞) /m π in each of the bound two-baryon channels are between 0.2 and 0.4, indicating that the range of interactions mediated by the pion exchange is not the only characteristic scale in the system, suggesting unnaturalness. However, at larger-than-physical quark masses, pion exchange may not be the only significant contribution to the long-range component of the nuclear force, as is discussed in Ref. [41]. For these reasons, both natural and unnatural interactions are considered in the next section when adopting a powercounting scheme in constraining the LECs of the EFT. 13 TABLE V. The values of the parametersã −1 ,r,P from a two-or three-parameter polynomial fit for twobaryon channels that exhibit smooth and monotonic behavior in k * cot δ as a function of k * 2 beyond the t-channel cut. Quantities are expressed in lattice units.

Two-parameter polynomial fit
Three-parameter polynomial fit

36
(+08)(+08) (−11)(−04) 12 In Appendix B, results for N N channels are compared with the previous scattering parameters obtained in Ref. [42] using the same correlation functions, as well as with the predictions obtained from low-energy theorems in Ref. [90]. Through a thorough investigation, the various tensions are discussed and resolved. 13 For a detailed discussion of naturalness in EFTs, see Ref. [126].
FIG. 6. k * cot δ values as a function of the c.m. momenta k * 2 , along with the bands representing the twoand three-parameter polynomial fits for two-baryon systems under the assumption that there is a smooth and monotonic behavior in k * cot δ as a function of k * 2 beyond the t-channel cut. Quantities are expressed in lattice units.
Although the ERE is only valid below the t-channel cut, one may still fit the k * cot δ values beyond this threshold using a similar polynomial form as the ERE in Eq. (14). To distinguish the "model" fit parameters from those obtained from the ERE, two-and three-parameter polynomials are characterized by two {ã −1 ,r} or three {ã −1 ,r,P } parameters. Such forms are motivated by the fact that in most channels, k * cot δ values as a function of k * 2 exhibit smooth and monotonic behavior beyond the t-channel cut, as is seen in Fig. 2. The only exceptions are the spin-singlet N N and ΣN channels, for which such a polynomial fit will not be performed. The results of this fit, using the same strategy as described above for ERE fits, are shown in Table V and Fig. 6. In the next section, the EFTs and approximate symmetries of the interactions will be utilized to make predictions for the inverse scattering length in channels for which ERE fits could not be performed, i.e., ΣN ( 3 S 1 ) and ΞΞ ( 3 S 1 ) channels, and in those cases, the scattering length is found consistent with theã −1 values obtained from this model analysis. It should be emphasized that such a polynomial fit beyond the t-channel cut is only one out of many applicable parametrizations of the amplitude, and a systematic uncertainty associated with multiple model choices and modelselection criteria needs to be assigned to reliably constrain the energy dependence of the amplitude at higher energies. 14

D. Binding energies
A negative shift in the energy of two baryons in a finite volume compared with that of the noninteracting baryons may signal the presence of a bound state in the infinite-volume limit. However, to conclusively discern a bound state from a scattering state, a careful inspection of the volume dependence of the energies is required. Lüscher's quantization condition can be used to identify the volume dependence of bound-state energies. 15 Explicitly, the infinite-volume binding momentum κ (∞) can be determined by expanding Eq. (10) in the negative-k * 2 region [63]: where Z 2 is the residue of the scattering amplitude at the bound-state pole. In this study, the boost vectors are d = (0, 0, 0) and d = (0, 0, 2), and the values of γ deviate from one at the percent level. 16 Therefore, all systems considered are non-relativistic to a good approximation. Only the first few terms in the sum in Eq. (19), corresponding to |m| ∈ {0, 1, √ 2}, are considered in the volume extrapolation performed below, with corrections that scale as O(e −2κ (∞) L ).
Alternatively, one can compute κ (∞) by finding the pole location in the s-wave scattering amplitude: To obtain κ (∞) , the scattering amplitude has to first be constrained using Lüscher's quantization condition as discussed in the previous subsection, and then be expressed in terms of an ERE expansion. This approach, therefore, requires an intermediate step compared with the first method, but does not require a truncation of the sum in Eq. (19). Results for the infinite-volume binding momenta κ (∞) are shown in Table VI. The columns labeled as d = (0, 0, 0) and d = (0, 0, 2) correspond, respectively, to fitting separately the values of k * with no boost, or with boost d = (0, 0, 2), using Eq. (19). The column labeled as d = {(0, 0, 0), (0, 0, 2)} is the result of fitting both sets of k * values simultaneously, i.e., imposing the same value for κ (∞) and Z 2 in both fits. The last column shows the κ (∞) values obtained using Eq. (20), with the parameters listed in Table III as obtained with a two-parameter ERE fit to k * cot δ. The results obtained with the different extractions of κ ∞ are seen to be consistent with each other within uncertainties. The largest difference observed is in the ΞΞ ( 1 S 0 ) channel, with a difference between the volume-extrapolation and pole-location results of around 1.5σ. The agreement between the two approaches suggests that the higher-order terms neglected in the sum in Eq. (19) are not significant.
The binding energy, B, is defined in terms of the infinite-volume baryon masses and binding momenta as where M (∞) i is the infinite-volume mass of baryon i obtained from Eq. (9). This quantity is computed for all systems that exhibit a negative c.m. momentum squared in the infinite-volume limit, i.e., those listed in Table VI. The binding energies in physical units are listed for these systems in Table VII. The binding energies of the two-nucleon systems computed here are consistent within 1σ with the values published previously in Ref. [42] using the same LQCD correlation functions. The same two-baryon systems studied here were also studied at m π ∼ 806 MeV in Ref. [31], and were found to be bound albeit with larger binding energies. While the results at m π ∼ 806 MeV were inconclusive regarding the presence of bound states in the 10 irrep, the ΣN ( 3 S 1 ) and ΞΞ ( 3 S 1 ) systems are found to be unbound at this pion mass. The results obtained in the present work can be TABLE VI. The infinite-volume binding momenta κ (∞) for bound states obtained by either using the extrapolation in Eq. (19) or from the pole location of the scattering amplitude as in Eq. (20). Quantities are expressed in lattice units.
TABLE VII. Binding energies for bound states in MeV. The values are obtained using κ (∞) from the volumeextrapolation method with a combined fit to d = (0, 0, 0) and d = (0, 0, 2) data. The uncertainty from scale setting is an order of magnitude smaller than the statistical and systematic uncertainties quoted.
13.1 combined with those of Ref. [31] obtained at m π ∼ 806 MeV to perform a preliminary extrapolation of the binding energies to the physical pion mass. 17 This enables a postdiction of binding energies in nature in cases where there are experimental data, and a prediction for the presence of bound states and their binding in cases where no experimental information is available. For systems with non-zero strangeness, experimental knowledge is notably limited in comparison to the nucleon-nucleon sector, and almost all phenomenological predictions are based on SU (3) flavor-symmetry assumptions as discussed in the introduction. There is a significant body of work devoted to building phenomenological models of two-baryon interactions based on one-bosonexchange potentials, such as the Nijmegen hard-core [129][130][131], soft-core (NSC) [132][133][134] and extended-soft-core (ESC) [135][136][137][138][139][140][141][142][143] models, as well as the Jülich [144][145][146] and Ehime [147,148] models. EFTs [19,[21][22][23][149][150][151][152][153][154] and quark models [155][156][157] have also been used to construct two-baryon potentials. A short summary of the results in the literature for the relevant channels with non-zero strangeness is as follows: • The 1 S 0 and 3 S 1 ΣN channels do not exhibit bound states in any of the models listed above.
The spin-singlet state behaves in a similar way to N N ( 1 S 0 ), and the interactions are slightly attractive, while those in the spin-triplet channel are found to be repulsive.
• For the ΞN ( 3 S 1 ) system, almost all the models find that the interactions are slightly attractive, but only a few exhibit a bound state. 18 Among the most recent results are "ESC08a" [138] which gives B = 0.9 MeV ‡ , and "ESC08c1" [139] which gives B = 0.5 MeV ‡ . There is one LQCD calculation of this system near the physical values of the quark masses performed by the HAL QCD collaboration [158] using a different method than the current work, and no bound state is observed.
• The "NSC97" model [134] finds a bound state for the ΣΣ ( 1 S 0 ) channel, with binding energies ranging from 1.53 to 3.17 MeV. χEFT at NLO [23] finds a binding energy between 0 and 0.01 MeV (no bound state is found with ESC or quark models in this channel).
• The ΞΣ ( 1 S 0 ) system is found to be bound in the "NSC97" model [134], with a binding energy between 3.02 − 16.5 MeV, and by χEFT [22,23], with a binding energy between 2.23 − 6.18 MeV at LO and 0.19 − 0.58 MeV at NLO. With the quark model "fss2" [156], although the interaction in this system is found to be attractive, no bound state is predicted (similar to the "ESC08c1" model [139]).
• Using one-boson-exchange potentials, with "NSC97" [134] the ΞΞ ( 1 S 0 ) state is bound with a binding energy between 0.1 − 15.8 MeV, and with "Ehime" [148] between 0.23 − 0.71 MeV (no bound state is found with "ESC08c1" [139]). χEFT [22,23] also finds this state to be bound with a binding energy of 2.56 − 7.27 MeV at LO and 0.40 − 1.00 MeV at NLO. The quark model "fss2" [156] does not find a bound state. In the ΞΞ ( 3 S 1 ) channel, no bound state is found with one-boson-exchange potentials, except for "Ehime" that finds a deeply bound state with a binding energy of 9 − 15 MeV. "fss2" [156] finds this channel to be repulsive.
The quark-mass dependences of multi-baryon spectra have not been studied extensively in the literature. For the octet-baryon masses, it was found that LQCD calculations performed with 2+1 dynamical fermions are consistent with a linear dependence on the pion mass at unphysical values of the quark masses, compared to the HBχPT prediction of quadratic dependence at LO [159][160][161]. Nonetheless, recent precision studies near the physical values of the quark masses appear to be more consistent with chiral predictions [162]. In the two-baryon sector the situation is more complicated. On the theoretical side, χEFT was used in Ref. [163] to extrapolate LQCD results to the physical point, assuming no dependence on the light quark masses for the LECs of the EFT (at a fixed order). The same premise was taken in Ref. [29] to determine the I = 3/2 ΣN interaction at LO, which was used to address the possible appearance of Σ − hyperons in dense nuclear matter. In the absence of a conclusive form for the quark-mass extrapolation of two-baryon binding energies, two naive expressions with linear and quadratic m π dependence were used in Ref. [40] to extrapolate the binding energy of H-dibaryon to its physical value. In Refs. [164][165][166], under the assumption that the H-dibaryon is a compact 6 valence-quark state (and not a two-baryon molecule), χEFT was used to extrapolate the binding energies, resulting in an unbound state.
Two analytical forms with different m π dependence are used here to obtain the binding energies at the physical light-quark masses, using the results presented in Ref. [31] at m π ∼ 806 MeV and those listed in Table VII for m π ∼ 450 MeV: where lin and B (1) quad are parameters to be constrained by fits to data. These fits are shown in Fig. 7, along with the experimental value and predictions at the physical point. The binding energies extrapolated to the physical point, i.e., B lin (m phys π ) and B quad (m phys π ), are summarized in Table VIII. 19 These extrapolations highlight some interesting features. The values obtained at the physical point are consistent with the experimental values for the N N channels. The rest of the binding predictions are at the same level of precision as the phenomenological results. The ΞΞ ( 1 S 0 ) and ΞN ( 3 S 1 ) channels are more consistent with being bound than the other channels, using both extrapolation functions. Moreover, the ΣN ( 3 S 1 ) channel was found not to support a bound state in this study, a conclusion that is in agreement with phenomenological models. The same conclu- sion holds for ΞΞ ( 3 S 1 ), noting that only in one model, namely "Ehime", a different conclusion is reached [148]. The spread of results and some contradictory conclusions in the models motivate the need for LQCD studies of these states at near-physical values of the quark masses in the upcoming years.

III. EFFECTIVE LOW-ENERGY INTERACTIONS OF TWO BARYONS
A. Leading and next-to-leading order interactions in the EFT Even though SU (3) flavor symmetry is explicitly broken in this study by the different values of the light-and strange-quark masses, it is still useful to classify the different two-(octet)baryon channels according to the SU (3) irrep that they belong to. In the spin-flavor decomposition of the product of two octet baryons with J P = 1 2 + , the 64 existing channels can be grouped into: The Note that the structure of the LQCD interpolating operators used in this study, i.e., single-point quark-level wavefunctions at the source, does not allow accessing channels in the 8 s irrep. Moreover, the state in the 1 irrep is a coupled flavor channel, which is excluded from this study given that a large number of kinematic inputs are required to constrain the corresponding coupled-channel scattering amplitudes. 20 The Lagrangian for the low-energy interactions of two octet baryons was first constructed in Ref. [80] using the heavy-baryon chiral EFT (HBχEFT) formalism, and consists of two-baryon contact operators at LO. These interactions have also been studied in chiral perturbation theory (χPT) in Refs. [19,167], where in addition to the momentum-independent operators at LO, the pseudoscalar-meson exchanges are included in the interacting potential. At LO, all terms in both HBχEFT and χPT are SU (3) symmetric. At NLO, there are two types of contributions: the SU (3)-symmetric interactions, obtained by the addition of derivative terms to the LO Lagrangian, and the SU (3) symmetry-breaking interactions, denoted by SU (3) in the following, that arise from the inclusion of the quark-mass matrix. The NLO extension of the two-baryon potential within χPT was first presented in Refs. [20,151], and includes interactions in higher partial waves.
In this paper, two-baryon systems are analyzed at low energies, therefore only s-wave interactions are considered. The LO Lagrangian of Ref. [80] is used, and the NLO contributions are formed to follow the organization of the LO terms. In other words, the same spin-flavor operator structure is preserved in the NLO Lagrangian, up to the inclusion of derivative operators and the quarkmass matrix. The EFT considered is therefore a pionless EFT [78,79] in the hypernuclear sector. The LO coefficients are known as Savage-Wise coefficients in the literature. This organization is different from that of Petschauer and Kaiser in Ref. [20], and while the notation used here to label the NLO LECs is the same as in Ref. [20], their meaning is different. The differences between the two organizations and the relations between both sets of SU (3) coefficients are presented in Appendix C. The full pionless EFT Lagrangian, up to NLO, is written as: with where only terms that contribute to s-wave interactions are included in the NLO Lagrangian L and χ is the quark-mass matrix, which can be written in terms of the meson masses using the Gell-Mann-Oakes-Renner relation [168]: where the constant B 0 is proportional to the quark condensate. In order to constrain the values of the LECs c i ,c i , and c χ i , the LO and NLO EREs of the inverse scattering amplitudes in the s-wave can be used. It is known that if the interactions between octet baryons are unnatural, that is r/a 1, a better justified power-counting scheme in the EFT is the KSW-vK (Kaplan, Savage and Wise [169,170] and van Kolck [78]) scheme, where at LO in the scattering amplitude, the contributions from LO momentum-independent operators are summed to all orders. With natural interactions, a power-counting scheme based on naive dimensional analysis is used and the expansion of the amplitude remains perturbative in the interaction strength, including for the LO interaction. As mentioned in Section II C, given the large uncertainties in the scattering parameters (in particular in the effective range), the ratio r/a shown in Table IV is not well constrained, and does not conclusively prove unnaturalness in all channels. Since in at least two channels the interactions seem unnatural, in the following both the natural and unnatural cases will be considered in expressing relations between LECs and the scattering parameters. These relations for each two-baryon channel can be separated into those which are momentum-independent, with contributions from LO and NLO SU (3) terms in the Lagrangian, and momentum-dependent, with only contributions from NLO SU (3) terms: where c (irrep) stands for the appropriate linear combinations of the c i LECs defined in the Lagrangian in Eq. (25). These relations are given in Table IX The same relations hold forc (irrep) , replacing and c χ B 1 B 2 are linear combinations of the c χ i LECs as given in Table IX. The variables a B 1 B 2 and r B 1 B 2 are the scattering length and effective range of the channel B 1 B 2 , and M B 1 B 2 is the reduced mass of that system. The renormalization scale µ depends on the naturalness of the interactions. For the natural case µ = 0, and Eqs. (31) and (32) correspond to a tree-level expansion of the scattering amplitude. For the unnatural case, the expansion does not converge for momenta larger than a −1 , and in the KSW-vK scheme µ is introduced as a renormalization scale for the s-channel two-baryon loops appearing in the all-orders expansion of the amplitude with LO interactions. Since a pionless EFT is used, a convenient choice is µ = m π (where m π ∼ 450 MeV is the mass of the pion obtained with the quark masses used in the LQCD study).
Two sets of inputs can be used to constrain the numerical values for the LECs: 1) the scattering parameters {a −1 , r} obtained from two-parameter ERE fits in Section II C, tabulated in Table III, can be used to compute LECs of both momentum-independent and momentum-dependent operators (method I), and 2) the binding momenta from Section II D can be used to compute the corresponding scattering length, related at LO by −a −1 + κ (∞) = 0, and this single parameter can be used to constrain the LECs of momentum-independent operators (method II). This second method is motivated by the fact that κ (∞) is extracted with higher precision than the parameters from the ERE fits, therefore enabling tighter constraints on the LECs of momentum-independent operators. The results for both types of LECs are presented in Table X, and are depicted in Fig. 8.
Results are presented in units of 2π/M B for the momentum-independent operators and 4π 2 /M 2 B for the momentum-dependent operators, where M B is the centroid of the octet-baryon masses, As can be seen from the values of the LECs that are obtained, the NLO SU (3) coefficients have large uncertainties, and are mostly consistent with zero, because the effective ranges used to constrain them have rather large uncertainties. Another feature of the results is that assuming the interactions to be unnatural leads to better-constrained parameters in general, as a non-zero scale µ in the left-hand side of Eqs. (31) and (32) reduces the effect of uncertainties on the scattering lengths (this was also observed in Ref. [31] for systems at m π ∼ 806 MeV). Furthermore, as is expected, the values obtained with method II have smaller uncertainties than the ones obtained from method I, given the more precise scattering lengths, although the method is limited to LO predictions. Another anticipated feature is that in the cases where the effective range is resolved from zero TABLE X. LECs of the momentum-independent and momentum-dependent operators as they appear in Table IX for m π I 11.9  Unnatural case (µ = m ⇡ ) ], where M B is the centroid of the octet-baryon masses. The gray-circle markers denote quantities that are extracted using the ERE parameters (method I), while black-square markers are those obtained from scattering lengths that are computed from binding momenta (method II).
within uncertainties (e.g., in the ΞΞ ( 1 S 0 ) channel), the values from method II are slightly different from those obtained from method I, indicating the non-negligible effect of the NLO effective-range contributions that are neglected with this method.
It should be noted that the input for scattering parameters is not sufficient to disentangle the LO SU (3) and NLO SU (3) coefficients in general, hence the c (irrep) + c χ B 1 B 2 entry in Table X   Method c (27)  The results are shown in Table XI and Fig. 9, along with the result for the ΣΣ channel for comparison purposes, as there is no contribution from SU (3) interactions for this channel at this order. From these results, it can be seen that the values of the symmetry-breaking coefficients c χ 3 − c χ 4 and c χ 1 − c χ 2 + c χ 11 − c χ 12 are compatible with zero. Together with the observation that the scattering lengths and binding energies in all of the systems are similar within uncertainties, it appears that the SU (3) flavor symmetry remains an approximate symmetry at the quark masses used in this study. These observations in the two-baryon sector are consistent with those in the single-baryon sector as presented in Ref. [42] at the same quark masses. There, it was found that the quantity , which is a measure of SU (3) flavor-symmetry breaking, is an order of magnitude smaller than its experimental value. 22 In Appendix C, the full list of relations needed to independently constrain all 24 different LECs that appear at LO and NLO are shown, demonstrating that the proper combinations of 18 twobaryon flavor channels are sufficient to extract all these LECs. These channels will be the subject of upcoming LQCD studies toward the physical values of the quark masses.

B. Compatibility with large-N c predictions
In the limit of SU (3) flavor symmetry and large N c , two-baryon interactions are predicted to be invariant under an SU (6) spin-flavor symmetry, with corrections that generally scale as 1/N c [81]. In the two-nucleon sector, this encompasses the SU (4) spin-flavor Wigner symmetry [173][174][175], with corrections that scale as 1/N 2 c . Under SU (6) group transformations, the baryons transform as a three-index symmetric tensor Ψ µνρ , where each SU (6) index is a pair of spin and flavor indices (iα). At LO, only two independent terms contribute to the interacting Lagrangian of two-baryon 22 The violation of the Gell-Mann-Okubo mass relation [171,172] where the baryon tensor can be expressed as a function of the octet-baryon matrices, B: 23 Here, α, β, γ, ω are flavor indices, i, j, k are spin indices, and the Levi-Civita tensor is either in flavor or spin space depending on the type and number of indices. A priori, the relative size of the Kaplan-Savage coefficients, a and b, is unknown, and only experimental data or LQCD input may constrain these LECs. As is seen in Eqs. (36) below, the contribution from the b coefficient to the LO amplitude is parametrically suppressed compared with that of the coefficient a. As a result, if b in Eq. (34) is comparable or smaller than a, there remains only one type of interaction that contributes significantly to the scattering amplitude, a situation that would realize an accidental SU (16) symmetry of the nuclear and hypernuclear forces. The first evidence for SU (16) symmetry in the two-(octet)baryon sector was observed in a LQCD study at a pion mass of ∼ 806 MeV [31], and the goal of the present study is to examine these predictions at smaller values of the light-quark masses. Such a symmetry is suggested in Ref. [82] to be consistent with the conjecture of maximum entanglement suppression of the low-energy sector of QCD. As in Sec. III A, the a and b coefficients can be matched to scattering amplitudes in a momentum expansion at LO. Since at least some of the SU (3) symmetry-breaking LECs c χ i were found to be consistent with zero in this study, one can assume an approximate SU (3) symmetry in general, and relate the SU               23 Those components of the field Ψ that correspond to decuplet baryons [81] have been neglected as they are not relevant to the low-energy scattering of two octet baryons. given irreps: In order to extract a and b, states in the 27 and 10 irreps can be combined with those in the 8 a irrep, allowing for six possible extractions. 24 The results are shown in Table XII and Fig. 10. As seen in Eqs. (36), the contributions from the b coefficient are suppressed by at least a factor of 3 compared with those from the a coefficient, thus the rescaled coefficient b/3 is considered. Considering that the results presented should be valid only up to corrections that scale as 1/N c , individual values of the coefficients a and b/3 obtained from different pairs of channels exhibit remarkable agreement, indicating that the SU (6) spin-flavor symmetry is a good approximation at these values of the quark masses. A correlated weighted average 25 of the results is obtained, following the procedure 24 Note that the ERE parameters were obtained in the previous section only for two-baryon channels belonging to the {27, 10, 8a} irreps. 25 The average of a series of values {xi} with uncertainties {σi} is computed as: where, since the different values of xi (and their uncertainties) are correlated, a 100% correlation is assumed when   introduced by Schmelling [176] and used by the FLAG collaboration [177], and is shown as the pink bands in Fig. 10. Given the uncertainty in b/3, no conclusion can be drawn about the relative importance of a and b/3. We will return to the question of the presence of an accidental SU (16) symmetry shortly. Given the extracted values of a and b/3, several checks can be performed, and several predictions can be made. The simplest check is to compute all of the LO SU (3) LECs, c (irrep) , using the relations in Eqs. (36). The results are shown in the first rows of Table XIII and the upper panels of Fig. 11. Columns with hashed backgrounds are the coefficients whose values were used as an input to make predictions for other coefficients, presented in panels with solid colored backgrounds. These input coefficients (c (27) , c (10) , and c (8a) ) can be reevaluated using the average values of a and b/3, which therefore gives back consistent values but with different uncertainties (for c (27) , the average of the values given in Table XI is computed). The large uncertainties in the c (8s) , c (1) , and c (10) coefficients are due to the fact that b/3, with a larger uncertainty than a, is numerically more important in these cases, see Eqs. (36). Additionally, the Savage-Wise coefficients c i can be computed by inverting the relations in Eqs. (33), and the resulting values are presented in the last rows of Table XIII and the lower panels of Fig. 11. Due to large uncertainties in the natural case, no conclusions can be made regarding the relative size of the coefficients. In the unnatural case and at the chosen value of the renormalization scale, the c 5 coefficient has a larger value than the rest of the coefficients. The relative importance of c 5 is a remnant of an accidental approximate SU (16) symmetry of s-wave two-baryon interactions that is more pronounced in the SU (3)-symmetric study with m π ∼ 806 MeV in Ref. [31]. It will be interesting to explore whether the remnant of this symmetry remains visible in studies closer to the physical quark masses.
The values of the SU (6) coefficients a and b allow predictions to be made for the scattering lengths of the systems that could not be constrained in this study by an ERE fit, namely the ΣN ( 3 S 1 ) and ΞΞ ( 3 S 1 ) channels. Using the c (irrep) coefficients computed previously, the relations in Eq. (31) can be inverted to obtain a −1 , assuming that the values of c χ i are negligible compared computing σ 2 average . For asymmetric uncertainties in xi, the following procedure is used to symmetrize them: a value xi = c  The gray-circle symbols denote quantities that have been extracted using the scattering parameters obtained from the ERE fit (method I), while black-square symbols denote those that are obtained from scattering lengths constrained by binding momenta (method II). The hashed background in the upper panels denotes coefficients whose values were used to constrain a and b, and hence are not predictions.
with those of c (irrep) (an observation that is only confirmed for given linear combinations of these LECs but is assumed to hold in general given the hints of an approximate SU (3) symmetry in this study). This exercise leads to consistent results for the inverse scattering length for systems for which the ERE allowed a direct extraction of this parameter, while it provides predictions for the channels shown in Table XIV. For the case of natural interactions, the scattering lengths are not constrained well, although they are consistent within uncertainties with those in the unnatural case, demonstrating the renormalization-scale independence of the scattering length. For the unnatural case, both methods are consistent and give rise to inverse scattering lengths that are positive and larger than those obtained for the rest of the systems studied in this work. This is in agreement with the parameters found when fitting the results for k * cot δ in these channels beyond the t-channel cut, see Table V.

IV. CONCLUSIONS
Nuclear and hypernuclear interactions are key inputs into investigations of the properties of matter, and their knowledge continues to be limited in systems with multiple neutrons or when hyperons are present. In recent years, LQCD has reached the stage where controlled first-principles studies of nuclei are feasible, and may soon constrain nuclear and hypernuclear few-body interactions in nature. The present work demonstrates such a capability in the case of two-baryon interactions, albeit at an unphysically large value of the quark masses corresponding to a pion mass of ∼ 450 MeV. It illustrates how Euclidean two-point correlation functions of systems with the quantum numbers of two baryons computed with LQCD can be used to constrain a wealth of quantities, from scattering phase shifts to low-energy scattering parameters and binding energies, to EFTs of forces, or precisely the LECs describing the interactions of two baryons. This same approach can be expected to be followed in upcoming computations with the physical quark masses, and its output, both in form of finite-volume energy spectra and constrained EFT interactions, can serve as input into quantum many-body studies of larger isotopes, both at unphysical and physical values of quark masses, see e.g., Refs. [85][86][87] for previous studies in the nuclear sector. By supplementing the missing experimental input for scattering and spectra of two-baryon systems, such LQCD analyses can constrain phenomenological models and EFTs of hypernuclear forces. In summary, the present paper includes a computation of the lowest-lying spectra of several twooctet baryon systems with strangeness ranging from 0 to −4. These results have been computed in three different volumes, using a single lattice spacing, and with unphysical values of the light-quark masses. The finite-volume nature of the energies provides a means to constrain the elastic scattering amplitudes in these systems through the use of Lüscher's formalism. Assuming small discretization artifacts given the improved LQCD action that is employed, our results reveal interesting features about the nature of two-baryon forces with larger-than-physical values of the quark masses. In particular, the determination of scattering parameters of two-baryon systems at low energies has enabled constraints on the LO and NLO interactions of a pionless EFT, both for the SU (3) flavorsymmetric and symmetry-breaking interactions. While the two-baryon channels studied in this work only allowed two sets of leading SU (3) symmetry-breaking LECs to be constrained, and those values are seen to be consistent with zero, the present study is the first such analysis to access these interactions, extending the previous EFT matching presented in Ref. [31] at an SU (3)symmetric point with m π = m K ∼ 806 MeV. Given the limited knowledge of flavor-symmetrybreaking effects in the two-baryon sector in nature, this demonstrates the potential of LQCD to improve the situation. Finally, the observation of an approximate SU (3) symmetry in the two-baryon systems of this work led to an investigation of the large-N c predictions of Ref. [81], through matching the LQCD results for scattering amplitudes to the EFT. In particular, the swave interactions at LO are found to exhibit an SU (6) spin-flavor symmetry at this pion mass, as also observed in Ref. [31] at a larger value of the pion mass. Both of the two independent spin-flavorsymmetric interactions at LO are found to contribute to the amplitude. Nonetheless, the extracted values of the coefficients of the LO SU (3)-symmetric EFT suggest a remnant of an approximate accidental SU (16) symmetry observed in the SU (3) flavor-symmetric study at m π ∼ 806 MeV [31]. It will be interesting to examine these symmetry considerations in the hypernuclear forces at the physical values of the quark masses, particularly given the conjectured connections between the nature of forces in nuclear physics and the quantum entanglement in the underlying systems [82]. While no attempt is made in the current work to constrain forces within the EFTs at the physical point, a naive extrapolation is performed using the results of this work and those at m π ∼ 806 MeV, with simple extrapolation functions, to make predictions for the binding energies of several two-baryon channels. The results for ground-state energies of two-nucleon systems are found to be compatible with the experimental values. Furthermore, stronger evidence for the existence of bound states in the ΞΞ ( 1 S 0 ) and ΞN ( 3 S 1 ) channels is observed compared with other two-baryon systems. Such predictions are in agreement with current phenomenological models and EFT predictions, and can be improved systematically as LQCD studies of multi-baryon systems progress toward physical values of the quark masses in the upcoming years.

ACKNOWLEDGMENTS
We would like to thank Joan Soto and Isaac Vidaña for enlightening discussions on the effective field theory formalism and on the physics of neutron stars and the hyperon-puzzle, respectively. We would like to also thank them, as well as André Walker-Loud, for comments on the first version of this manuscript.  [178] and QUDA [179,180] software suites. We thank André Walker-Loud and Thomas Luu for contributions during initial stages of production prior to 2014.
After this manuscript was completed, an additional preprint [181] appeared regarding twonucleon scattering, using a set of scattering interpolating operators similar to Ref. [50].  -Consistency between ERE parameters for k * 2 < 0 and k * 2 > 0: In the two-baryon channels studied in this work, there are not sufficient data points for k * cot δ below the t-channel cut to extract precise scattering parameters, as pointed out in Sec. II B. Nonetheless, for the cases for which two sets of data at positive and negative values of k * 2 are available, the ERE fits obtained by fitting to all k * 2 versus only fitting to k * 2 < 0 values are fully consistent with each other, as is shown in Fig. 13.
-Non-singular scattering parameters: None of the scattering parameters extracted show singular behavior, as can be seen from the values in Table III. -Requirement on the residue for the scattering amplitude at the bound-state pole: In order to support a physical bound state, the slope of the ERE as a function of k * 2 must be smaller than the slope of the − √ −k * 2 at the bound-state pole. The two slopes and associated uncertainty bands are depicted in Fig. 14 for all two-baryon channels and the two-parameter EREs obtained, demonstrating that the needed inequality is satisfied. The values of binding momenta used in this analysis are taken from Table VI (the d = {(0, 0, 0), (0, 0, 2)} column).
-The absence of more than one bound state with an ERE parametrization of amplitudes: None of the systems analyzed exhibit more than one bound state, i.e., the ERE does not cross the − √ −k * 2 curve more than once. Therefore, applying the ERE parametrization of the s-wave scattering amplitude in all channels appears to be justified.
-Constrained range for ERE parameters in the presence of a bound state: If the system presents a bound state, the ratio r/a must be smaller than 1/2 for the two-parameter ERE to cross the − √ −k * 2 function once from below, which is the condition for a physical bound state. Moreover, the ERE must cross the Z-functions corresponding to different volumes to satisfy Lüscher's quantization condition, introducing more constraints on scattering parameters. With the use of the two-dimensional χ 2 in this work to fit the k * cot δ values, the confidence region of the ERE parameters does not cross these prohibited areas, as was demonstrated in Fig. 4.  Table VI. Quantities are expressed in lattice units.

Appendix B: Comparison with previous LQCD results and those obtained from low-energy theorems
A subset of the correlation functions used in this work has already been analyzed in Ref. [42], where the N N ( 1 S 0 ) and N N ( 3 S 1 ) channels were studied. In the following, we present the outcome of a careful comparison of the results obtained using both analyses, along with a comparison of the updated scattering parameters from this work and those obtained from low-energy theorems in Ref. [182].

Differences in the fitting strategy
The ground-state and first excited-state energies obtained in this work and those from Ref. [42] are shown in Fig. 15. While all numbers are in agreement within uncertainties, it is clear that, examining the individual fits from all accepted time windows. These are shown in Fig. 16 for the N N ( 1 S 0 ) L = 24 and L = 32 ground states, sorted by their weight, w f , as defined in Eq. (4). As can be seen, there are cases for which the size of the uncertainty is similar to or smaller than that presented in Ref. [42]. However, the final combined uncertainty, represented by the band in Fig. 16, is larger. This can be understood as using a more conservative procedure for quantifying the systematic uncertainty, as well as a more thorough one: not only are variations of the fitting range considered, but also variations in the fitting form, including forms with multiple exponentials, see Sec. II B. Next, the implications of using the HL estimator (instead of the mean) on the individual SP and SS correlation functions are analyzed. When correlations are fully taken into account, the covariance matrix associated with the HL estimator is computed with the Median Absolute Deviation (MAD): whereC(τ ) is the bootstrap ensemble computed with the HL estimator of the original correlation function C(τ ). However, in some cases the resulting covariance matrix is found not to be positive semi-definite, and it only becomes well behaved when a single type of correlation function is used (or a linear combination of several) in the form of an effective-(mass) energy function. To illustrate this, Fig. 17 shows the normalized inverse covariance matrix, C −1 (τ, τ )/ C −1 (τ, τ )C −1 (τ , τ ), for the N N ( 1 S 0 ) ground state with L = 24 and L = 32 for all possible choices, i.e., HL estimator versus mean and correlation function versus the effective-energy function. Therefore, in order to incorporate the HL estimator into the fitting strategy used here, only the fully uncorrelated covariance matrix can be used, and this leads to results which are compatible with the ones presented here using the mean. In Fig. 18, the effective-energy functions computed with the mean and HL are compared for the N N ( 1 S 0 ) first excited states, showing agreement within uncertainties.
To understand the ill-behaved behavior of some of the HL correlation functions, is important to recall that baryonic correlation functions exhibit distributions that are largely non-Gaussian with heavy tails, and the mean becomes Gaussian only in the limit of large statistics. However, at late times, the signal-to-noise degradation worsens, and outliers occur more frequently in the distribution. For the L = 32 and L = 48 cases, the point at which the HL estimator gives different results compared with the usual estimator (mean and standard deviation), which would indicate a deviation from Gaussian behavior, occurs at a much later time compared with the maximum time included in the fits using the automated fitter of this work. For the L = 24 case, the data is more noisy than on the other two ensembles, showing non-Gaussianity at earlier times. To illustrate the different behavior between the L = 24 and L = 32 ensembles, the second and third cumulants of C(τ ), defined as: with n ∈ {2, 3}, respectively, are shown in Fig. 19 for the two ensembles in the case of the N N ( 1 S 0 ) first excited state. Looking at the second cumulant (variance), κ 2 , it is clear that L = 24 is more noisy than L = 32, and looking at the third cumulant (skewness), κ 3 , it is clear that L = 24 deviates from zero, an indication of the non-Gaussian behavior. The use of robust estimators is, therefore, questionable in this case. This is the main reason for abandoning the use of the HL estimator in the analysis of correlation functions in the present study.

Differences in the scattering parameters
The 68% confidence region of the scattering parameters from a two-parameter ERE extracted in this work and in Ref. [42] are shown in Fig. 21. It can be seen that the values of the parameters obtained in the two analyses do not fully agree at the 1σ level, although the uncertainties are rather large. There are two significant differences between the two analysis: 1) the use of the new definition for the χ 2 function (2D-χ 2 ) in the present work, as opposed to the usual χ 2 function (1D-χ 2 ) used in Ref. [42], and 2) the use of the L-dependent ground-state k * 2 values in the fits to ERE in the present work, instead of using only the infinite-volume extrapolated value, κ (∞) , used in Ref. [42]. To see the effects of each, a comprehensive analysis has been performed, the results of which are shown in Fig. 20. Here, four different possibilities, corresponding to the types of the χ 2 function (1D or 2D) and the use of ground-state k * 2 data (L-dependent or extrapolated), are tested using the lowest-lying spectra obtained in Ref. [42] and those in the present work.
From these tests, several interesting features are observed. First, the use of the 1D-χ 2 , either with the L-dependent k * 2 or the extrapolated one, is insensitive to the conditions imposed by Lüscher's  quantization condition, and as a result, the confidence regions of the scattering parameters could lie on top of the prohibited regions. This is because the distance minimized in the 1D-χ 2 is the vertical one, and not the one along the Z-function, so the ERE is not forced to cross it. Second, when the 2D-χ 2 is used with the extrapolated k * 2 value, κ (∞)2 , the only region that is avoided is the one corresponding to L = ∞ in the figures, which is expected: with the value of the pole position given by Eq. (20), the function k * cot δ| k * =iκ (∞) equals − √ −k * 2 and the ERE crosses the − √ −k * 2 function, imposing the r/a < 1/2 constraint on the scattering parameters. Third, it is reassuring that the regions obtained using the two different energy inputs, from this work or from Ref. [42], are always overlapping.
Perhaps the most significant observation is that the choice of including the points in the negative k * 2 region in the fit, i.e., the infinite-volume extrapolated value of the momenta versus the Ldependent values, has far more impact on the differences observed than which χ 2 function is used. What the new χ 2 function does is to move the scattering parameters to the allowed region by the Z-functions. Furthermore, with the new fitting methodology, several questions raised about the validity of the ERE fits are addressed, as was presented in Appendix A. An important one is that the updated results of this work recover the position of the bound state pole obtained via the infinite-volume extrapolation of the energies, and do not yield a second pole near threshold, which would be incompatible with the use of the ERE. As a final remark, it should be noted that the data fitted to extract these parameters are highly non-Gaussian, as can be seen from the correlation between k * 2 and k cot δ in Fig. 2, and exhibit large uncertainties. This can be compared with the results of Refs. [31,41] at m π ∼ 806 MeV, where more finite-volume energy eigenvalues, with better precision, could be used in the ERE fitting. As a result, it has been verified that either the L-dependent or the infinite-volume extrapolated value of k * 2 in the ERE fitting give compatible scattering parameters.
In Ref. [90], low-energy theorems [182] were used to compute the scattering parameters from the binding energies of the N N systems obtained in Ref. [42], and it was pointed out that there were some tensions with the scattering parameters obtained from the LQCD data using Lüscher's method, i.e., those reported in Ref. [42]. Since the binding energies obtained in this work are in full agreement with those obtained in Ref. [42], the results obtained in Ref. [90] can be compared with the updated scattering parameters of this work. As is depicted in Fig. 21, the tension has reduced considerably. For the two-parameter ERE results, the scattering length is now completely consistent with the low-energy theorem predictions, both at LO and NLO. For the effective range, since the NLO predictions of the low-energy theorems enter the prohibited region for the twoparameter ERE, the comparison may only be made with the LO results. As is seen, for both the 1 S 0 and 3 S 1 channels, the effective ranges are also in agreement (with the 1 S 0 state having a better overlap).
Appendix C: On leading flavor-symmetry breaking coefficients in the EFT Table 10 of Ref. [20] lists the SU (3) flavor-symmetry-breaking LECs c i χ for all of the two-(octet) baryon channels. These coefficients are a combination of different terms in the Lagrangian shown in Table 9 of the same reference (terms [29][30][31][32][33][34][35][36][37][38][39][40]. The relations between the c i χ from Ref. [20] and the ones in Eq. (28) introduced in the present work are presented in Table XV. Instead of the ( 2s+1 L J , I) notation, the channels are labeled as ( 2I+1 2s+1 ) for brevity, as L = 0 in all cases. In Table XVI, a list of the two-baryon channels one needs to study in order to obtain independently all the LECs of this work is provided. There are 6 LO and 12 NLO symmetrybreaking coefficients that are referred to as momentum independent in this paper, as well as 6 NLO momentum-dependent coefficients, making a total of 24 parameters that need to be constrained in a more exhaustive study in the future. For the momentum-independent coefficients, the choice of the systems is not unique, as there are 37 different channels that can be used to constrain only 18 parameters (assuming SU (2) flavor symmetry and no electromagnetic interaction). For the momentum-dependent coefficients, no extra channels are needed besides those used for the momentum-independent coefficients. For simplicity, only channels that do not change the baryon content are used (e.g. ΣN → ΣN , denoted as ΣN in short).

Appendix D: Supplementary figures and tables
This appendix contains all the figures omitted from the main body of the paper for ease of presentation. These include the effective-mass plots of the single baryons in Fig. 22, and the effective energy and effective energy-shift plots for the two-baryon systems in Figs. 23-31. In Fig. 22, the thin horizontal line and the horizontal band surrounding it represent, respectively, the central value of the baryon mass at each volume, and the associated statistical and systematic uncertainties combined in quadrature, obtained with the fitting procedure described in Sec. II B. Similarly, in Figs. 23-31 the line and the band represent, respectively, the central value of the two-baryon energy shifts compared to non-interacting baryons at rest (bottom panels) for each volume, and the associated statistical and systematic uncertainties combined in quadrature.
The appendix also contains the numerical results that were omitted from the main body. These include the energy shifts, ∆E, of the two-baryon systems, the c.m. momenta, k * 2 , and the value of k * cot δ for all the systems in Tables XVII-XXV. In these tables, the values in the first and second parentheses correspond to statistical and systematic uncertainties, respectively, while those in the upper and lower parentheses are, respectively, the right and left uncertainties when the error bars are asymmetric, as is generally the case for the k * cot δ values. When there is a dash sign in the tables, it indicates that the quantity k cot δ diverges due to the singularities in the Z d 00 function. All quantities in the plots and tables are expressed in lattice units.  [20] for the two-baryon channels for which only one c i χ appears in that reference.