Two-Nucleon Higher Partial-Wave Scattering from Lattice QCD

We present a determination of nucleon-nucleon scattering phase shifts for l>= 0. The S, P, D and F phase shifts for both the spin-triplet and spin-singlet channels are computed with lattice Quantum ChromoDynamics. For l>0, this is the first lattice QCD calculation using the Luscher finite-volume formalism. This required the design and implementation of novel lattice methods involving displaced sources and momentum-space cubic sinks. To demonstrate the utility of our approach, the calculations were performed in the SU(3)-flavor limit where the light quark masses have been tuned to the physical strange quark mass, corresponding to m_pi = m_K ~ 800 MeV. In this work, we have assumed that only the lowest partial waves contribute to each channel, ignoring the unphysical partial wave mixing that arises within the finite-volume formalism. This assumption is only valid for sufficiently low energies; we present evidence that it holds for our study using two different channels. Two spatial volumes of V ~ (3.5 fm)^3 and V ~ (4.6 fm)^3 were used. The finite-volume spectrum is extracted from the exponential falloff of the correlation functions. Said spectrum is mapped onto the infinite volume phase shifts using the generalization of the Luscher formalism for two-nucleon systems.


I. INTRODUCTION
Understanding low-energy nuclear physics directly from the underlying theory of strong interactions, Quantum ChromoDynamics (QCD), remains a primary goal of nuclear physicists. The motivation can be broadly separated into two categories: obtaining a quantitative description of basic nuclear physics directly from QCD and probing the limits of the Standard Model and its fundamental symmetries through precision low-energy experiments in nuclear environments. In both cases, there are substantial international experimental efforts planned or underway which require a quantitative understanding of QCD.
The basic interactions of two nucleons (NN) and nuclei are well measured and have led to a variety of precise theoretical descriptions ranging from phenomenological models to effective field theories (EFT). These over-constrained NN interactions are then used to predict properties of light nuclei using a variety of methods such as multi-nucleon EFT [1][2][3], harmonic oscillator based effective theory [4,5], no-core shell model [6,7] and Green's Function Monte Carlo [8,9]. For light nuclei, the NN interactions dominate the nuclear structure, but the three-body nuclear force is necessary for accurate comparisons with the measured values [10,11]. A recent, exciting development is the use of lattice field theory to regularize the two-and three-nucleon EFT and predict properties of light nuclei, such as the recent computation of the Hoyle State [12].
All of these impressive theoretical applications rely upon experimental input of the nuclear interactions. While this is achievable precisely for the low-energy NN interactions, there are very few constraints on the three and higher nucleon forces. Determining these interactions directly from QCD is a multifaceted problem. At the core of this challenge is the non-perturbative nature of QCD which requires a numerical approach at low energy. Lattice QCD (LQCD) is the discretized version of QCD in a finite Euclidean volume. It is the only known tool to compute QCD correlation functions in the infrared for which no uncontrolled approximations are necessary. Using LQCD combined with EFT (see for example Refs. [13,14]), we will be able to determine the elusive few-body nucleon interactions directly from QCD, relevant for example for the upcoming experiments at FRIB [15], designed to study neutron-rich nuclei. We will also be able to compute the hyperon-nucleon interactions, which are extremely challenging to measure experimentally due to the rapid weak decay of the hyperons [16][17][18][19]. These calculations will be relevant for the experiments planned at the FAIR, JLab and J-PARC facilities. There are even studies of hyper-nuclei using heavy ion collisions at RHIC [20] and the LHC [21].
LQCD is also necessary to compute one, two and few-body nuclear matrix elements such as the scalar matrix elements needed for direct dark matter detection and the electroweak matrix elements which govern nuclear interactions and decays. Two recent examples of these include the N → N π parity violating process [22] and the parity conserving np → dγ rate [23]. It is known that, due to technical complications, calculations of matrix elements involving multiparticle states often receive O(1) corrections from the finite volume [24][25][26]. In order to control these corrections, the two-particle phase shifts and their derivatives must be determined. See Ref. [27] for a very nice demonstration of this technology for the case of πγ → ππ.
In this work, we present a calculation of higher partial wave scattering in the NN system. In particular, we have computed S,P,D and F partial waves in both isosinglet and isotriplet channels using the NN generalization [69] of Lüscher's formalism. Given the complex nature of this problem, this exploratory calculation was performed at the SU (3) flavor symmetric point with m π ∼ 800 MeV, enabling us to explore the implementation of our new method and demonstrate its feasibility with relatively little investment of computing time. We have also simplified the Lüscher formalism by ignoring mixing from higher partial waves contributing to a given cubic irrep. This mixing is kinematically suppressed at sufficiently low scattering energies, however, for the range of energies we explore the assumption that they do not contribute significantly may require further investigation. Some evidence that the assumption is valid can be obtained by comparing different cubic irreps coupling to the same partial wave, and is presented in Section III. This work is an extension of a previous determination of the NN S-wave interactions [53] on the same LQCD gauge configurations.

II. IMPROVED TWO-NUCLEON INTERPOLATING FIELDS
We have performed the calculations with the isotropic Wilson lattices generated by the JLab/W&M group, with a lattice spacing of b = 0.145(2) fm and two spatial extents of L/b = 24 and L/b = 32 (for more details see Ref. [57]). Since the volumes are cubic and the two-nucleon systems considered have zero total momentum, the spectra obtained are those of the irreducible representations (irreps) of the octahedral symmetry group O h .
The sink operators are defined as products of single nucleon operators in momentum space. The single nucleon operators were designed to have good overlap with the single nucleon ground states [70][71][72]. Let R be an element of O h . Let N m I 1 ms 1 (k) be a nucleon operator with spin and isospin z-components m s1 and m I1 and momentum k at time t. Due to the periodic boundary conditions, the free momenta satisfy k = 2πn/L, where n is an integer triplet. With these we can construct two nucleon operators with angular momentum Jm J and isospin Im I where and S denote the total orbital angular momentum and spin of the system and Rk is the unit vector in the Rk direction and the Y m ( Rk) are the standard spherical harmonic functions. The standard Clebsch-Gordan coefficients, C Jm J m ,Sm S = Jm j | m l , Sm s project the operators to total spin S, angular momentum J, and isospin I. The infinite volume quantum numbers are not good quantum numbers in a finite volume. Operators with different angular momentum labels will mix due to the non-spherically symmetric finite spatial boundary conditions. To project the operators above to an operator that is in a row µ of the irrep Λ of O h , we use the subduction coefficients, (1) and go directly to a set of operators belonging to definite cubic irreps, it is convenient to keep separate operators of the same cubic irrep having different [J S] labels, as these operators in some cases have very different overlap with the various excited states of a particular cubic irrep. For example, in the spin singlet, T + 1 channel, the second non-zero momentum shell of the non-interacting system contains two degenerate states, whose energies split once interactions are turned on. We find that the T + 1 operators having = 0 labels exhibit good overlap with the lower of these two states, while operators with = 2 labels overlap well onto the higher state (see Fig. 1).
Ideally, we would construct NN operators in momentum space at both the source and sink locations, as is done in two-meson calculations [28,29]. However, the computational cost of performing these calculations for two nucleons is orders of magnitude greater than for two mesons. For this reason, NN calculations typically are performed with local or volume sources for the nucleons at the source and then projected to definite momentum at the sink. The significant improvement made over previous works is the use of spatially displaced nucleon operators at the source, N † (t 0 , x 0 + r/2)N † (t 0 , x 0 − r/2). By displacing the two nucleons at the source, we find a significant increase in the overlap of the operators onto the NN A + 1 and T + 1 irreps (∼ S-wave) as compared to the local operators (see Fig. 2, left). Further, without such a displacement, for zero total momentum, the overlap of the local operators (r = 0) with the cubic irreps that contain the P, D and F waves are zero, prohibiting a determination of the spectrum in these irreps.
Displaced operators give us further freedom in designing our sources to have good overlap with the desired states by choosing from various geometries for the displacements. For example, after fixing one nucleus to a single lattice point, (0, 0, 0), one may then calculate the set of pairs where the second nucleon is displaced in all possible ways along a single axis, r = |r|(0, 0, 1) (plus all cubic rotations), or along multiple axes, such as r = |r|(0, 1, 1) (plus all cubic rotations) and r = |r|(1, 1, 1) (plus all cubic rotations). We have named these geometries "face", "edge", and "corner", respectively, and disregard more complicated geometries. Each collection of geometries may then be projected onto the desired cubic irrep as described above for the sink operators (Eq. (1)), with the momentum vectors k replaced by the set of displacements vectors, r. Note, however, that this is only a partial projection because we use a reduced set of displacement vectors compared to the full lattice volume. Example effective masses for the three types of geometries are shown in Fig. 2 (right). While the "face" sources have a reduced computational cost compared to "corner" and "edge" sources (7 inversions versus 9 ("corner") and 13 ("edge")), they have zero overlap with several channels of interest due to their simple geometrical structure. We have chosen to focus on "corner" sources for the remainder of this work, to balance good overlap with a large number of states with moderate computational cost.
In Fig. 3, we show effective mass plots for several operators in two different cubic irreps along with the resulting determination of the energy levels. We construct all sink operators that have free momenta |k|L/2π ≤ √ 6. Note that the T − 2 channel has overlap with the physical 3 P 2 , 3 F 2 and 3 F 4 channels, and the T + 2 channel has overlap with the physical 1 D 2 channel. Both channels have additional overlap with kinematically suppressed higher partial waves. The bands denote the combined statistical and systematic uncertainties obtained by performing one and two exponential

III. NN SCATTERING PHASE SHIFTS
In general, due to the reduction of rotational symmetry in a finite volume, there is not a one-to-one correspondence between the finite volume spectrum and the infinite volume scattering amplitudes. For sufficiently low energies, one can expect higher partial waves to be kinematically suppressed. Ignoring partial wave mixing and ≥ 4 waves, the In most panels we plot q 2 +1 cot δ which is used to determine the parameters in the ERE. In the 1 S0 channel, the faded points were not included in the ERE fit, but are displayed here to show consistency. In the upper right panel we also show the phase shift δ 3 P 2 as a function of the lattice momenta.
spectra of a number of cubic irreps satisfy the quantization condition [69] q cot δ Λ (q) 4π = c 00 (q 2 ) + α 4,Λ c 40 (q 2 ) q 4 + α 6,Λ where q is the on-shell relative momentum of the system, q 2 = is the scattering phase shift of the partial wave that primarily couples to the Λ irrep, α ,Λ are constants, Ref. [69], and the c m are kinematic, non-linear functions that depend solely on the momentum and the volume Employing Eqs. (3) and (4), we obtain the scattering phase shifts evaluated at the on-shell relative momenta derived from the spectrum. In Fig. 4 we give illustrative examples of the quality of the results for the spin-triplet and spinsinglet channels in various irreps that couple predominantly to a given partial wave. The bands are fits to the effective range expansion (ERE) to different orders, i.e.
where a , r , and P are the scattering length, effective range, and shape parameter for = 0 and the corresponding parameters of the ERE for > 0. Results for LO (q 0 ) fits are denoted by yellow bands, fits to NLO (q 2 ) are blue, and  NNLO (q 4 ) fits are red. The ERE parameters determined from these fits are listed in Table I. The dashed vertical line represents the t-channel cut, above which the ERE is expected to break down. While the ERE becomes formally unjustifiable past this cut, the Lüscher formalism holds for all energies below the N N π threshold, well above the energies considered. We obtain results significantly different from zero for various channels, including P and D channels. One beautiful illustration of the power of the operators and the finite volume formalism being used is that of 3 P 2 in Fig. 4 (topmiddle and top-right). In this channel the spectrum was determined with two irreps that have overlap with the same partial wave using two different volumes. The consistency of the extracted phase shifts clearly demonstrates that the generalization of the Lüscher formalism for NN-systems is working [69]. We find that over the kinematic range of our calculations, q 3 cot δ 3 P2 is consistent with a constant, even above the t-channel cut. Furthemore, we find no evidence of the t-channel cut playing an important role for any of the channels studied at these values of the quark masses.
Currently, we ignore partial wave mixing, an issue that will be addressed in subsequent publications. This formidable challenge has only been addressed in two-meson calculations [28,29,32,[46][47][48]. However, some evidence that the mixing from higher partial waves is small in at least one channel can be obtained by investigating the results for the 3 P 2 channel in Fig. 4. Again, for this channel we have two cubic irreps, T − 2 and E − , for which the lowest partial wave is 3 P 2 . Any differences between the two irreps must arise from mixing with higher partial waves. We find that the two cubic irreps give completely consistent results, indicating that this mixing must be too small to resolve within our error bars. Additionally, the first contribution from mixing in the T − 1 cubic irrep, having lowest partial wave 1 P 1 , comes from the 1 F 3 partial wave. We obtain information on the strength of the 1 F 3 phase shift independently using the A − 2 cubic irrep, and find it to be extremely small, thus it likely does not contribute to mixing in the T − 1 cubic irrep.

A. Bound states from the Effective Range Expansion
Assuming the results are within the range of convergence of the ERE, the infinite volume bound state energies may be determined by solving for poles in the derived scattering amplitude. For a single channel scattering amplitude, these satisfy which has a solution below threshold for q B = iκ B . In both S-wave channels, the ERE expansion appears to be converging well (top-left and middle-left figures in Fig. 4) with small O(q 4 ) corrections where the results exist below the t-channel cut. Using this method and the spectrum obtained in this work, we find deeply bound states in both the 1 S 0 and 3 S 1 channels, with binding energies of The first uncertainties are our fitting statistical and systematic uncertainties combined in quadrature and the second is an estimate of systematic uncertainties arising from the truncation of the ERE. In addition, in the 3 S 1 channel we find a second pole near threshold, corresponding to Corresponding to each of these poles we find finite volume states whose energies are consistent with the expected exponential volume dependence associated with bound states. However, with only two volumes we cannot definitively state whether the volume dependence is exponential or polynomial. With this precision, it is unclear whether this L unitarity bound T + 1 ( 3 S1) A  shallow bound state corresponds to a true bound state or a near-threshold scattering state. Improved analysis techniques, such as a full basis of interpolating fields in momentum space, as has been successfully used in the twomeson systems [28,29,[46][47][48], or additional statistics may be necessary to fully settle this matter. We do note, however, that all of the negative energy finite-volume shifts are larger by 4 or more standard deviations than the corresponding finite-volume energy shift for a bound state at threshold, ∆E ∼ − 3.786 M L 2 [74,75], as shown in Table II. This indicates that the interactions producing these states are more attractive than those at unitarity, and the states should therefore correspond to true bound states in the infinite volume limit. Further evidence of multiple bound states are presented in the next section.

B. Evidence for multiple negatively shifted energy states
We further elaborate on the plausibility of finding two negatively shifted energy states in both S-wave channels. The large negatively shifted energy levels, determined with the local operators, are consistent with those in Ref. [57], which also used local operators. 2 The state closer to threshold (and additionally, the negative energy state near threshold in the 1 S 0 channel) has strong overlap onto the non-local NN interpolating field, and has not been found in previous works. In Fig. 5, we plot the effective mass of the ratio correlation functions of the non-local two-nucleon interpolating fields divided by the local two-nucleon interpolating fields Both two-nucleon interpolating fields couple to the same tower of states, E Λ n = 2m N + ∆ Λ n and only differ in the relative size of their overlap factors, A Λ n (|r|). In the spin triplet, T + 1 channel, there is a clear plateau inconsistent with 0 for a relatively large time in fm (recall that each time step on these ensembles correspond to a = 0.1453(16) fm [57]), indicating that the two correlation functions each give statistically distinct plateaus. For the spin singlet, A + 1 channel, a clear plateau is only observed on the larger volume.
In the long-time limit, the effective masses for all interpolating operators in this channel must asymptote to the ground state of the system. However, if the non-local interpolating field couples strongly to the state with small negative energy shift, and weakly to the state with large negative energy shift, and conversely for the local interpolating field, then at intermediate times the dominant contribution to each interpolating field will come from the two different states respectively. This would manifest in the non-local interpolating field having an effective mass plateau at the smaller negative energy shift (the first excited state) at intermediate times. With enough statistics, one would eventually observe the effective mass "collapse" to the ground state and plateau again at the larger negative energy shift (the ground state). Note that this is also true of all excited two-nucleon elastic scattering states. Because relative momentum is not a good quantum number, the non-zero relative momentum projections serve only to enhance the overlap of the interpolating operators with excited states relative to the ground state. We find this to work quite well due to the existence of distinct plateaus at intermediate times for the different relative momentum operators, however, in the very long time limit we expect all of these correlation functions to collapse onto a single ground state plateau.
This explanation for the two states with negative energies is consistent with the intuition that the wavefunction for a shallow bound state is more extended in space, and thus will have poor overlap with an interpolating field involving only two nucleons at the same point (a smeared out delta function in our case). Conversely, the more deeply bound state would have a much more compact wave function and would thus overlap more poorly with a non-local interpolating field. We finally emphasize that these findings are consistent with the two poles found using the ERE of our phase shift, even if we disregard the finite volume negative energy results from either the local or the non-local operators. We therefore find the most plausible explanation of these results to be the existence of two distinct bound states in the T + 1 channel.

IV. SUMMARY
This work presents the implementation of new two-nucleon interpolating fields which allow, for the first time, a robust determination of > 0 scattering phase shifts in the NN sector, which we have calculated in this exploratory work using unphysically heavy quark masses. Further, this improved basis of interpolating operators are sensitive to additional states in the S-wave spectrum that were not found using only local operators and greater statistics. This has been made possible by three previously unexploited tools. First was the development of displaced two-nucleon interpolating sources. These are necessary to have appreciable overlap with partial waves beyond the S-wave as the = 0 orbital wavefunctions are zero at the origin. Second was the use of momentum space sink operators that were not restricted to the simplest cubic irreps. Finally, we applied the formalism for two-nucleon systems in a finite volume [69], ignoring higher partial wave mixing, with notable success for the 3 P 2 channel. This work represents the first crucial step towards the study of more challenging systems such as three-neutron interactions and the S → P wave parity violating pp interaction.