Elastic nucleon-pion scattering at $m_{\pi} = 200~{\rm MeV}$ from lattice QCD

Elastic nucleon-pion scattering amplitudes are computed using lattice QCD on a single ensemble of gauge field configurations with $N_{\rm f} = 2+1$ dynamical quark flavors and $m_{\pi} = 200~{\rm MeV}$. The $s$-wave scattering lengths with both total isospins $I=1/2$ and $I=3/2$ are inferred from the finite-volume spectrum below the inelastic threshold together with the $I=3/2$ $p$-wave containing the $\Delta(1232)$ resonance. The amplitudes are well-described by the effective range expansion with parameters constrained by fits to the finite-volume energy levels enabling a determination of the $I=3/2$ scattering length with statistical errors below $5\%$, while the $I = 1/2$ is somewhat less precise. Systematic errors due to excited states and the influence of higher partial waves are controlled, providing a pathway for future computations down to the physical light quark masses with multiple lattice spacings and physical volumes.


Introduction
Nucleon-pion (N π) scattering is a fundamental nuclear physics process. Because the pion is the lightest hadron, pion exchange between nucleons governs the long-range nuclear force and contributes to the binding of protons and neutrons into atomic nuclei. Nucleon-pion scattering also gives rise to the narrow ∆(1232) resonance which influences many nuclear processes, including lepton-nucleon and lepton-nucleus scattering relevant to a range of electron-nucleus and neutrino-nucleus scattering experiments.
While N π scattering is well understood experimentally and phenomenologically, such as through the Roy-Steiner equations [1], the ability to determine the amplitudes directly from quantum chromodynamics (QCD) is hampered by its non-perturbative nature at low energies. After QCD was established as the underlying theory of the strong nuclear force, chiral perturbation theory (χPT) [2] and chiral-EFT [3,4] were developed to systematically describe the low-energy dynamics of pions and nucleons in an effective field theory (EFT) framework. For a recent review, see Ref. [5]. While EFT methods are generally effective in treating low-energy hadron scattering processes, a number of challenges can only be addressed with first-principles QCD calculations, for which lattice QCD is an essential non-perturbative tool.
For example, many of the low-energy constants (LECs) of nuclear EFTs are difficult to determine from experimental information alone. Lattice QCD can assist in the determination of LECs by carrying out computations at a variety of quark masses and by computing processes which are difficult to measure experimentally, such as hyperon-nucleon and three-nucleon interactions, as well as short-distance matrix elements of electroweak and beyond-the-Standard Model operators. See recent reviews for further discussion [6][7][8][9][10]. The beneficial interplay between EFTs and lattice computations is already developing for meson-meson scattering [11][12][13][14][15][16], but few lattice studies of meson-baryon scattering amplitudes currently exist.
Another important issue concerns the convergence of EFTs, which are asymptotic expansions in small momenta and light quark masses with convergence not guaranteed at the physical quark masses. Lattice QCD has already provided numerical evidence that SU (2) baryon χPT is not converging at or slightly above the physical pion mass for the nucleon mass and axial coupling [17][18][19]. Including explicit ∆ degrees of freedom may improve convergence of SU (2) baryon χPT, but introduces a plethora of additional unknown LECs. Lattice QCD calculations of N π scattering at various pion masses can help verify the convergence pattern and whether it is improved with explicit ∆s [20], as well as constrain the additional LECs. N π scattering is additionally important because of the current tension between lattice QCD determinations of the nucleon-pion sigma term [21][22][23][24] σ πN and phenomenological determinations [1,25] (see Ref. [26] for a possible resolution). σ πN plays an important role in the analysis of direct dark matter detection experiments [27]. Controlled lattice QCD calculations of N π scattering may help understand this tension.
As a final example, a future prospect for lattice QCD is the determination of inputs to models of neutrino-nucleus scattering cross sections to aid next-generation experiments, like DUNE [28] and Hyper-K [29], designed to measure unknown properties associated with neutrino oscillations. The importance of lattice QCD input was recently highlighted by current lattice QCD results for elastic nucleon form factors [30]. The frontier for these lattice QCD applications is the ∆-resonance and pion-production contributions to inelastic νN structure. To carry out this program, it is essential to first demonstrate control of N π scattering, a necessary component of nucleon inelastic resonant structure.
Lattice QCD calculations of two-pion systems are well established (for a recent review, see Ref. [31]), and there are now numerous three-meson results [32][33][34][35][36][37][38]. In contrast, there are few nucleon-pion scattering studies in lattice QCD. Ref. [39] included nucleon-pion operators to extract the spectrum in the I = 1/2 channel at a single lattice spacing using a 16 3 × 32 lattice volume with m π ≈ 266 MeV and obtained an estimate of the scattering length with a significance of roughly four standard deviations. Refs. [40] and [41] each employ a single ensemble with m π 250 MeV to evaluate scattering phase shifts relevant to the ∆ resonance, but neither presented statistically significant results for the scattering lengths. There is also older work which employs the quenched approximation [42] and preliminary unpublished results for the I = 3/2 amplitudes [43][44][45]. The determination of finite-volume nucleon-pion energies in Ref. [46] is performed close to the physical quark masses, but scattering amplitudes are not computed. Lattice computations of mesonbaryon scattering lengths in other systems have also been performed [47,48].
where the errors are statistical only. The Breit-Wigner parameters for the ∆(1232)-3 resonance are also determined from the I = 3/2, J P = 3/2 + partial wave m ∆ m π = 6.257 (35), g ∆,BW = 14.41 (53), where the corresponding scattering phase shift is shown in Fig. 7. Since only a single ensemble of gauge field configurations is employed, the estimation of systematic errors due to the finite lattice size, lattice spacing, and unphysically large light quark mass is left for future work. However, systematic errors due to the determination of finitevolume energies, the reduced symmetries of the periodic simulation volume, and the parametrization of the amplitudes are addressed. The methods presented here therefore provide a step toward the lattice determination of the nucleon-pion scattering lengths at the physical point with controlled statistical and systematic errors. The remainder of this work is organized as follows. Sec. 2 discusses the effects of the finite spatial volume, including the corresponding reduction in symmetry and the relation between finite-volume energies and infinite-volume scattering amplitudes. Sec. 3 presents the computational framework, including the lattice regularization and simulation, the measurement of correlation functions, and the determination of the spectrum from them. Results for the amplitudes are presented and discussed in Sec. 4, while Sec. 5 concludes.

Finite-volume formalism
The Euclidean metric with which lattice QCD simulations are necessarily performed complicates the determination of scattering amplitudes. It was shown long ago by Maiani and Testa [59] that the direct application of an asymptotic formalism to Euclidean correlation functions does not yield on-shell scattering amplitudes away from threshold. Instead, lattice QCD computations exploit the finite spatial volume to relate scattering amplitudes to the shift of multi-hadron energies from their non-interacting values [60]. See Ref. [61] for a more complete investigation of the Maiani-Testa theorem, and Refs. [62,63] for an alternative approach to computing scattering amplitudes from Euclidean correlation functions based on Ref. [64].
This section summarizes the relationship between finite-volume spectra and elastic nucleon-pion scattering amplitudes. Due to the reduced symmetry of the periodic spatial volume, this relationship is not one-to-one and generally involves a parametrization of the lowest partial wave amplitudes with parameters constrained by a fit to the entire finitevolume spectrum. Symmetry breaking due to the finite lattice spacing is also present, but ignored. At fixed physical volume and quark masses, the continuum limit of the finite volume spectrum exists and is assumed for this discussion.
For a particular total momentum P , the relationship between the finite-volume center-of-mass energies E cm determined in lattice QCD and elastic nucleon-pion scattering amplitudes specified in the well-known K-matrix is given by the determinantal equation whereK is proportional to the K-matrix and B P (E cm ) is the so-called box matrix, using the notation of Ref. [57]. This relationship holds below the nucleon-pion-pion threshold, up to corrections which vanish exponentially for asymptotically large M L, where L is the 4 side length of the cubic box of volume L 3 and M the smallest relevant energy scale. The determinant is taken over all scattering channels specified by total angular momentum J, the projection of J along the z-axis m J , and the orbital angular momentum . For elastic nucleon-pion scattering the total spin S = 1/2 is fixed, and therefore not indicated explicitly. The K-matrix is diagonal in J and m J , and, for elastic nucleon-pion scattering, additionally diagonal in . TheK-matrix in Eq. (3) explicitly includes threshold-barrier factors which are integral powers of q cm = q 2 cm , with so thatK −1 is smooth near the nucleon-pion threshold. Each diagonal element ofK is associated with a particular partial wave specified by J P , where P is the parity, or equivalently (2J, ), so that where δ J (E cm ) is the scattering phase shift. The box matrix B P (E cm ) encodes the reduced symmetries of the periodic spatial volume, and is in general dense in all indices. The finite-volume energies used to constrain K from Eq. (3) possess the quantum numbers associated with symmetries of the box, namely a particular irreducible representation of the finite-volume little group for the total spatial momentum P = 2π L d, with d ∈ Z 3 . The matrices in Eq. (3) are therefore block-diagonalized in the basis of finite-volume irreducible representations (irreps), with each energy analyzed using a single (infinite-dimensional) block. Since the subduction from infinite-volume partial waves to finite-volume irreps is not in general one-to-one, an additional occurrence index n is required to specify the matrix elements in each block. A particular block is denoted by the finite-volume irrep Λ(d 2 ) and a row of this irrep λ. Since the spectrum is independent of the row λ, this index is henceforth omitted. For a particular block, the block-diagonalized box-matrix is denoted B J n,J n . After transforming to the block diagonal matrix, theK matrix has the form given by Eq. (35) in Ref. [57].
In practical applications, the matrices in Eq. (3) are truncated to some maximum orbital angular momentum max . Threshold-barrier arguments ensure that, at fixed E cm , higher partial waves are suppressed by powers of q cm , but systematic errors due to finite max must be assessed. The expressions for all elements of B Λ (d 2 ) relevant for this work are given in Ref. [57], although some are present already in Ref. [65]. The occurrence pattern of lowest-lying partial waves in the finite-volume irreps is given in Table 1.
Employing this formalism for nucleon-pion scattering presents additional difficulties compared to simpler scattering processes. First, due to the non-zero nucleon spin, two partial waves contribute for each non-zero , one with J = + 1/2 and the other with J = − 1/2. Secondly, off-diagonal elements of the box matrix induce mixings of different partial waves in the quantization condition. For max = 2, energies E cm in irreps with d 2 = 0 determine the quantity q 2 +1 cm cot δ J (E cm ) for sand p-waves, while these partial waves cannot be unambiguously isolated for levels in irreps with non-zero total momentum. This complication necessitates global fits of all energies to determine the desired partial waves, which are discussed in Sec. 4. 5 , (5, 2) Table 1: A list of the lowest contributing partial waves for each irrep of the finite-volume little group Λ in momentum class d employed in this work. All partial waves with ≤ max for max = 2 are shown and each partial wave is denoted by (2J, ). The superscript nocc denotes the number of multiple occurrences (subductions) of the partial wave in the irrep. The pattern of partial wave mixing is evidently more complicated for irreps with non-zero total momentum.

Spectrum computation details
This section details the numerical determination of finite-volume nucleon-pion energies used to constrain the ≤ 2 partial waves of the I = 1/2 and I = 3/2 elastic nucleon-pion scattering amplitudes. Properties of the single ensemble of gauge field configurations are given in Sec. 3.1, and computation of the nucleon-pion correlation functions from them is discussed in Sec. 3.2. The subsequent determination of the finitevolume spectra from the correlation functions is detailed in Sec. 3.3.

Ensemble details
This computation uses the D200 ensemble of QCD gauge configurations generated by the Coordinated Lattice Simulations (CLS) consortium [66], whose properties are summarized in Table 2. It was generated using the tree-level improved Lüscher-Weisz gauge action [67] and a non-perturbatively O(a)-improved Wilson fermion action [68]. Open temporal boundary conditions [69] are employed to reduce the autocorrelation of the global topological charge. However, all interpolating fields must be sufficiently far from the boundaries to reduce spurious contributions to the fall-off of temporal correlation functions. An analysis of the zero-momentum single-pion and ρ-meson correlators in Ref. [70] suggests that a minimum distance of m π t bnd 2 is sufficient to keep temporal boundary effects below the statistical errors in the determination of energies. The time ranges for the correlators employed here are such that m π t bnd 2.3.
A complete description of the algorithm used to generate the D200 ensemble is presented in Ref. [66], but some details relevant for the present work are given below. All CLS ensembles use twisted-mass reweighting [71] for the degenerate light quark doublet and the Rational Hybrid Monte Carlo (RHMC) approximation for the strange quark [72]. Both representations of the fermion determinants require reweighting factors to change the simulated action to the desired distribution. All primary observables are therefore a[fm] (L/a) 3 × T /a N meas am π am K 0.0633(4)(6) 64 3 × 128 2000 0.06617(33) 0.15644 (16) af π af K am N 0.04233(16) 0.04928(21) 0.3148(23) Table 2: Parameters of the D200 ensemble produced by the CLS consortium [66]. The lattice spacing a is from Ref. [74] (with both statistical and systematic errors) and follows the strategy of Ref. [73]. The number of gauge configurations employed here is specified by Nmeas. The pion mass mπ and nucleon mass m N determinations are discussed in Sec. 3.3. The kaon mass, denoted m K , and the pion and kaon decay constants, denoted fπ and f K , are taken from Ref. [73].

re-weighted according to
where ... W denotes the ensemble average with respect to the simulated action. W is the product of two factors W = W 0 W 1 , where W 0 and W 1 are the reweighting factors for the light and strange quark actions. They are estimated stochastically on each gauge configuration as in Ref. [66]. The lattice scale is determined at a fixed value of the gauge coupling according to the massless scheme described in Ref. [73] and updated in Ref. [74]. Specifically, the kaon decay constant f K is enforced to take its physical value at the physical point where the pion and kaon masses take their physical values. This point is identified along a trajectory in which the bare light-and strange-quark masses are varied, keeping the sum of the (renormalized) quark masses fixed. The heavier-than-physical pion mass m π = 200 MeV for the D200 therefore results in m K = 480 MeV, which is less than the physical value. In practice, the bare quark mass tuning satisfies the trajectory condition only approximately. In order to correct any mistuning a posteriori, Ref. [73] applies slight shifts to the quark masses to ensure the trajectory condition is respected in the scale determination. No such shift is applied here.
In this study, correlation function measurements are separated by four molecular dynamics units (MDU's). To check for autocorrelations, the original measurements are binned by averaging N bin consecutive gauge configurations. The dependence of the relative errors on N bin for the single-nucleon and single-pion correlators is shown in Figure  1. Although evidence of autocorrelation remains for t/a 8 − 10 between N bin = 20 and 40, these early timeslices are not used in the analysis, suggesting that N bin = 20 is sufficient to account for any autocorrelations in our energy determinations.

Correlation function construction
The determination of finite-volume nucleon-pion energies requires a diverse set of temporal correlation functions measured on the D200 gauge field ensemble. In addition to diagonal correlation functions between single-pion and single-nucleon interpolating operators, correlation matrices between all operators in each irrep are required. For the I = 3/2 irreps in Table 1 where the resonant (2J, ) = (3, 1) partial wave contributes, single-baryon operators are included in addition to nucleon-pion operators resulting in  Table 3: Parameters of the stochastic LapH implementation used to compute temporal correlators in this work. N D is the number of Dirac matrix inversions required per configuration and (ρ, nρ) the stout smearing parameters for the spatial links in the gauge-covariant Laplace operator. Nev denotes the dimension of the LapH subspace. N fix R and N rel R are the number of stochastic sources for fixed and relative quark lines, respectively. Notation used to specify the dilution scheme for each line type is explained in the text, and the number of source times on each configuration is Nt 0 .
additional valence quark-line topologies. These topologies include those with lines that start and end on the same timeslice.
Our operator construction is described in Ref. [75] and our method of evaluating the temporal correlators is detailed in Ref. [50]. Well-designed multi-hadron interpolators are comprised of individual constituent hadrons each having definite momenta. Evaluating the temporal correlations of such operators requires all-to-all quark propagators, where all elements of the Dirac matrix inverse are computed. The stochastic-LapH approach [50] enables the efficient treatment of this inverse, provided at least one of the quark fields is LapH smeared [49]. This smearing procedure is effected by a projection onto the space spanned by the N ev lowest eigenmodes of the gauge-covariant three-dimensional Laplace operator in terms of link variables which are stout smeared [76] with parameters (ρ, n ρ ). Although the N ev required to maintain a constant smearing radius grows with the spatial volume, the growth of the number of Dirac matrix inversions N D can be significantly reduced with the introduction of stochastic estimators in the LapH subspace. Such estimators are specified by the number of dilution projectors in the time ('T'), spin ('S'), and Laplacian eigenvector ('L') indices, for which 'F' denotes full dilution and 'In' some number of uniformly interlaced projectors. Different dilution schemes are used for fixed-time quark lines, denoted 'fix', which propagate between different timeslices, and relative-time lines ('rel') which start and end at the same time. In this work, the 8 relative-time quark lines were only used at the sink time, while the fixed-time lines were used for quark propagation starting and ending at the source time. Both the dilution schemes and the number of stochastic sources used for each type of line are specified in Table 3. Source times t 0 = 35, 64 were used for correlations going forwards in time, and t 0 = 64, 92 were used for correlations going backwards in time.
A beneficial property of the stochastic estimators is the factorization of the inverse of the Dirac matrix, which enables correlation construction to proceed in three steps: (1) Dirac matrix inversion, (2) hadron sink/source construction, and (3) correlation function formation. After determining the stochastically-estimated propagators in step (1), the hadron source/sink tensors are computed in step (2). These tensors are subsequently reused to construct many different correlation functions in step (3), which consists of optimized [32] tensor contractions. Averages over the N t0 = 4 different source times (two for forward propagation and two for backward propagation), all possible permutations of the available noise sources in a given Wick contraction, all total momenta of equal magnitude, and all equivalent irrep rows are performed to increase statistics.

Determination of finite-volume energies
Once the correlation functions computed as described in Sec. 3.2 are available, the determination of finite-volume energies can commence. From the (binned) correlator and reweighting factor measurements, the reweighted correlation functions are computed as secondary observables according to Eq. (6). Their statistical errors and covariances are used in fits to determine energies and estimated by the bootstrap procedure with N B = 800 samples.
In order to ensure that t bnd is sufficiently large, a maximum time separation t max = 25a is enforced globally in the analysis. Energies are determined from correlated-χ 2 fits to both single-and two-exponential fit forms, which are additionally compared to a "geometric series" form which consists of four free parameters. We also explored a "multi-exponential" variant of the geometric series, with the replacement Be −M t → n B n e −Mnt . The application of our approach to determining the nucleon and pion masses is shown in Fig. 2. As usual, a fit range is desired so that statistical errors on the energies are larger than systematic ones. This optimal range is selected according to several criteria. First, a good fit quality q 0.2 − 0.3 is enforced to ensure that the fit describes the data within the usual 68% confidence interval quoted for statistical errors. Second, the absence of any statistically significant change in the energy upon variation of t min around the chosen fit range further suggests that the asymptotic large-time behavior is applicable. Finally, consistency across different fit forms supports the hypothesis that the energy determination is statistics limited. For the pion, consistency between single-and twoexponential fits, as well as the mild variation with t min , suggests that statistical errors are dominant. As is evident in am eff N , single-exponential fits for am N are unsuitable, but consistency between the double-exponential and geometric fit forms is reassuring.
For the nucleon-pion channels, excited state energies are determined in addition to ground states. This requires a large basis of interpolating operators in each irrep and  two-point correlations between them. The resulting N op ×N op hermitian matrix, denoted C ij (t), is rotated [77] using eigenvectors v n (t 0 , t d ) of the generalized eigenvalue problem (GEVP) In our single-pivot approach, the correlation matrix is rotated by these vectors for all t, and the diagonal elements of the rotated matrix, denoted D n (t), are correlators with optimal overlap onto the lowest N op states. Although diagonalizing separately on each time-slice ensures that the optimized correlators are increasingly dominated by the desired state for asymptotically large times [78,79], in practice, the single-pivot method produces nearly identical results if (t 0 , t d ) are chosen appropriately. Systematic errors related to this are controlled by ensuring that extracted energies are insensitive to (t 0 , t d ) and N op and by ensuring that the rotated correlation matrix remains diagonal within statistical errors for all time separations t > t d . The advantages of diagonalizing on a single set of times include a better signal-to-noise ratio for large times and no need for eigenvector pinning in which eigenvectors are re-ordered for diagonalizations at different times and bootstrap samples. After forming the optimized correlators, the following ratio is taken with d 2 π and d 2 N chosen so that corresponds to the closest non-interacting energy. The ratio R n (t) is then fit to the single-exponential ansatz R n (t) = A n e −∆Ent to determine the energy shift a∆E n , from which the lab-frame energy is reconstructed aE lab n = a∆E n + aE non.int.
n . Although the ratio fits enable somewhat smaller t min when ∆E n is small, they offer little advantage for states which are significantly shifted from non-interacting levels. Nonetheless, ratio fits are employed for all levels in the nucleon-pion irreps, and are typically consistent with single-and double-exponential fits directly to D n (t).
A sample illustration of the procedure for nucleon-pion energies is shown in Fig. 3 for the second level in the I = 1/2 G(3) irrep. Due to partial wave mixing, the single-nucleon state is also present in this irrep. The GEVP is therefore required to properly isolate the desired higher-lying nucleon-pion energies. Analogous plots for all levels are given in Appendix A and Appendix B for the GEVP-and t min -stability plots, respectively.
The spectra resulting from this analysis are shown in Figs. 4a and 4b for the I = 1/2 and I = 3/2 channels, respectively.

Scattering parameter results and discussion
This section details the determination of the scattering parameters from the finitevolume energies. The parameterizations of the K matrix elements are presented, and  best-fit values for the parameters are summarized. Lastly, a comparison with chiral perturbation theory is made.

Scattering parameter determinations
The energies shown in Figs. 4a and 4b are next used to determine scattering amplitudes via the relations in Sec. 2. Although these relations are only applicable to energies below the nucleon-pion-pion threshold, the slow growth of three-body phase space near threshold suppresses corrections to Eq. (3) and the coupling of nucleon-pion-pion states to our operator basis is naively suppressed by the spatial volume, so energies somewhat above the inelastic threshold are expected to be appropriate for inclusion in our global fits. Nevertheless, we restrict our attention to center-of-mass energies below or within one standard deviation of the threshold E cm = 2m π + m N .
The goal of this analysis is a parametrization of the J P = 1/2 − partial wave for both isospins, and the 3/2 + wave with I = 3/2. As discussed in Sec. 2, energies from 12 irreps with zero total momentum directly provide points for these partial waves up to corrections from ≥ 3 contributions. However, mixing among the desired waves, as well as with others, generally occurs for energies in irreps with non-zero total momentum. The zero-momentum points are therefore a useful guide when plotted together with the partial wave fits. Each partial wave is parametrized using the effective range expansion. For the I = 3/2, J P = 3/2 + wave, the next-to-leading order is included where √ s = E cm = m 2 π + q 2 cm + m 2 N + q 2 cm , and the effective range fit parameters are reorganized to form the conventional Breit-Wigner properties of the ∆(1232) resonance, denoted g 2 ∆,BW and m ∆ . For the other waves, the leading term in the effective range expansion is sufficient where the overall √ s factors are adopted from standard continuum analysis [80], and the single fit parameter A I J P is trivially related to the scattering length Two different fit strategies are employed to determine the parameters from the finitevolume energies. The first, called the "spectrum method" [81], obtains best-fit values of the model parameters {p n } by minimizing where the q 2 cm,i are the center-of-mass momenta squared computed from lattice QCD, with covariance matrix C, and q 2,QC cm,i ({p n }) are the center-of-mass momenta squared evaluated from the model fit form for a given choice of parameters {p n }. The fact that the model depends on m N /m π , and is therefore not independent of the data to be fit, complicates the evaluation of the covariance matrix C. As discussed in Ref. [57], a simple way to avoid this complication so that C is just the covariance matrix of the data is to introduce model parameters for m π L and the ratio m N /m π , and include appropriate additional terms in the χ 2 of Eq. (14). Given the relatively small errors on m N /m π and m π L, these additional terms have little effect on the fit parameters and the resultant χ 2 , and are subsequently ignored. Note that the evaluation of the (q 2,QC cm,i /m 2 π )({p n }) for a particular choice of the parameters requires the determination of roots of Eq. (3), a procedure which can be delicate for closely-spaced energies.
The second method, called the "determinant residual" method [57], employs the determinants of Eq. (3) themselves as the residuals to be minimized. These determinants depend on the fit parameters through the K-matrix, which are adjusted to minimize the residuals and best satisfy Eq. (3). This approach avoids the subtleties associated with root-finding, but has other difficulties. For the spectrum method, the covariance between 13  Table 4: Results for the fits in the I = 3/2 channel. Npw is the number of partial waves included in the fit. Two different fit forms are included, the one denoted Npw = 2 includes only the desired partial waves, namely J P = 1/2 − and 3/2 + , while the one with Npw = 5 includes all s-, p-, and d-waves, employing the two additional energy levels in the G 1g (0) and Hu(0) irreps. For the Npw = 2 fit, results from the determinant-residual method, denoted 'DR', are shown in addition to the spectrum method, denoted 'SP'. (12) Fig. 4a. Npw is the number of partial waves included in the fit. Due to the small number of levels, all fits include only the desired J P = 1/2 − partial wave. Nonetheless, the effect of the omitted p-waves is estimated by removing the G 1 (4) level, which evidently has little influence on the result. 'SP' refers to the spectrum method, and 'DR' refers to the determinant-residual method.
the residuals is, to a good approximation, simply the covariance between the q 2 cm,i /m 2 π , which can be estimated once and does not depend on the fit parameters. Conversely, for the determinant residual method, the covariance must be re-estimated whenever the parameters are changed. Since the statistical errors on the determinant are typically larger than those on q 2 cm,i /m 2 π , this approach is less sensitive to higher partial waves, and results in a smaller χ 2 compared to the spectrum method.
For the I = 3/2 fits, the J P = 1/2 + , 3/2 − , and 5/2 − partial waves are added to the spectrum method fits along with the ground states in the G 1g (0) and H u (0) irreps. The I = 3/2 spectrum in the G 2g (0) irrep was not computed, and irreps in the I = 1/2 channel which do not contain the s-wave were also omitted. This choice was made for computational simplicity, although these irreps may be beneficial to further constrain higher partial waves in future work. The determinant residual method was found to be less able to constrain higher partial waves and was only used in fits that included just the J P = 1/2 − , 3/2 + waves. Nonetheless, the consistency between these different fitting methods, as well as those including higher partial waves, suggest that uncertainties on amplitude parameters are statistics dominated.
For the I = 1/2 channel, max = 0 is employed. Although the small number of levels precludes a sophisticated estimate of the effect of higher partial waves, the influence of the omitted p-waves can be explored by examining the influence of the highest-lying level on the fit. Table 5 indicates that the effective range is insensitive to the omission of the lowest-lying nucleon-pion level in the G 1 (4) irrep. These I = 1/2 fits are also insensitive to an additional term in the effective range expansion, and exhibit no statistically significant difference between the spectrum and determinant-residual methods.
Results from fits using both the spectrum and determinant-residual methods including 14 various partial waves are given in Tables 4 and 5 for I = 3/2 and I = 1/2, respectively. In addition to the desired partial waves, fits using the spectrum method are mildly sensitive to the J P = 1/2 + , 3/2 − , and 5/2 − waves with I = 3/2. Although not included in the table, the determination of the effective range for both isospins is robust to the addition of the next term in the effective range expansion. Results for the partial waves from the fit including only the desired partial waves are shown with the points from the total-zero momentum irreps in Figs. 5 and 6 for the I = 3/2 and I = 1/2 partial waves, respectively. The phase shift δ 3/2 + has the characteristic profile of the ∆(1232) resonance and is shown in Fig. 7. Since the scattering length is the only desired parameter from the I = 1/2 spectrum, only the lowest nucleon-pion levels from each irrep are included in the fit, as denoted by the solid symbols in Fig. 4a. Full exploration of the elastic I = 1/2 spectrum likely requires additional operators beyond the scope of this work, due to the strongly-interacting J P = 1/2 + wave containing the N (1440) Roper resonance. The spectrum method enables an additional visualization of the quality of fits to the finite-volume spectra. The residual is constructed using model values of q 2,QC cm /m 2 π which depend on the parameters and can be compared with the input data from the spectrum. Such comparisons are shown in Fig. 8 for both the I = 1/2 and I = 3/2 spectra. Although not shown explicitly on the plot, the ground states in G 1 (1), G(2), G(3), and G 1 (4) with I = 3/2 are sensitive to the J P = 3/2 + partial wave. The max = 0 approximation significantly increases the χ 2 for these levels. Conversely, these levels therefore place significant constraints on the near-threshold behavior of the 3/2 + wave, in contrast to the higher-lying levels in the H g (0), G 2 (1), F 1 (3), F 2 (3), and G 2 (4) irreps. The ground states in the G 1g (0) and H u (0) irreps are not shown on the plot, and only included in the N pw = 5 fit in Table 4.
The final results for the scattering lengths in this work are taken from the determinant residual method fit in Table 4 with N pw = 2 for I = 3/2 and the spectrum method fit which are already given in Eq. (1). In Fig. 9, the results from this work for the Breit-Wigner parameters of the ∆(1232) resonance in the I = 3/2, J P = 3/2 + partial wave are compared to the published numbers in Refs. [40] and [82] where, as is customary, the definition of the g ∆N π coupling from leading-order effective field theory is used, as defined in Eq. (39) of Ref. [82]. When considering Fig. 9, keep in mind that the quark mass trajectory used here and in Ref. [40], which fixes the sum of the quark masses, differs from that used in Ref. [82], which fixes the strange quark mass to its physical value. A comparison of the scattering lengths determined here to past lattice QCD results is also shown in Fig. 9.

Comparison with phenomenology and chiral perturbation theory
Although the results here are only at a single pion mass, it is interesting to compare our scattering lengths to those extracted phenomenologically [83,84] from pionic atoms [85][86][87], as well as to those obtained from chiral perturbation theory. The proximity of m π = 200 MeV to its physical value naively suggests that SU (2) baryon χPT may be applicable. While χPT may be poorly converging for g A and m N , the convergence pattern is an observable dependent issue which has not been explored for N π scattering.
Pionic hydrogen (πH) is sensitive to one combination of the isoscalar (a + ) and isovector (a − ) scattering lengths and pionic deuterium (πD) is sensitive to a different combination, allowing for a percent level determination [83,84]. Phenomenological values extracted from a full πN partial wave analysis have been reported in Ref. [88]. Ref. [27] determined the values of these scattering lengths in the isospin limit. While it is customary for the lattice community to use m π = 135 MeV as the pion mass in the isospin limit, it is common in the phenomenological estimates to use the charged pion mass to determine the QCD quantities in the absence of QED corrections. In order to be consistent with the phenomenological estimates, we similarly use the charged pion mass when quoting results at the physical pion mass. 16 These scattering lengths are known to fourth order in the baryon chiral expansion [89][90][91] and expressed in Appendix F of Ref. [88] and Ref. [92] in a form convenient for extrapolating LQCD results. In terms of the s-wave a ± 0 scattering lengths, the isospin 1/2 and 3/2 πN scattering lengths are given by At leading order (LO), the scattering lengths are free of LECs and given by The values of these input parameters on D200 and at the physical (charged) pion mass are A comparison of our results with the LO χPT predictions and phenomenological values in the isospin limit from Ref. [27] is presented in Table 6. Not only do our results disagree with LO χPT, but we also find the magnitude of m π a 3/2 0 exceeds that of m π a 1/2 0 , in conflict with both LO χPT and phenomenology. Note that the LO χPT predictions at     Table 6: A comparison of our N π scattering length results at mπ = 200 MeV with phenomenological values in the isospin limit and predictions from leading order chiral perturbation theory. For the χPT predictions, the first error is from uncertainties on the input parameters, π and µ, and the second error is a χPT truncation uncertainty given by π mπa I 0 [LO]. the physical point are in reasonable agreement with the phenomenological values, lying within one sigma of the estimated χPT truncation uncertainty. With only one pion mass available in this work, the reasons for the discrepancy of our results with LO χPT cannot be ascertained. Interestingly, our scattering length results can be described at next-to-leading order (NLO) using a single LEC. At NLO, one finds m π a 3/2 where Λ χ = 4πF π and we have defined the dimensionless LEC in terms of the c i LECs in the baryon chiral Lagrangian [95]. The scattering lengths in this work can be described by these NLO formulae if C is in the range 0.6-0.7. The NLO phenomenological determination finds a value of C ≈ 0.3, which is not significantly different from that needed to describe our results. However, the phenomenological extraction of the LECs in Ref. [88] is clouded by issues related to the ∆ degrees of freedom [20] and is not stable until at least next-to-next-to-next-to-leading order (N 3 LO) [88]. When results at additional pion masses, particularly lighter ones, become available, a more thorough understanding of the pion mass dependence of the scattering lengths can be achieved and a more quantitative comparison with the results from the phenomenological analysis and χPT can be performed.

Conclusion
This work presents a computation of the lowest partial waves for the elastic nucleonpion scattering amplitude on a single ensemble of gauge configurations with m π = 200 MeV. The s-wave scattering lengths are determined for both isospins I = 1/2 and I = 3/2 and compared to determinations from LO χPT and a Roy-Steiner analysis [88]. To our knowledge, this is the first (unquenched) lattice QCD determination of both nucleon-pion scattering lengths for m π < 250 MeV. The Breit-Wigner resonance parameters of the ∆(1232) in the J P = 3/2 + partial wave with I = 3/2 are determined as well.
A comparison of two different methods, the spectrum method and the determinant residual method, of extracting K-matrix information from finite-volume spectra is also performed. Although the determinant residual method avoids awkward root-finding, it was found to be less sensitive to higher partial wave contributions. Nonetheless, the consistency between these two different fitting procedures is reassuring.
These results suggest that the methods used here will prove useful for future work at the physical values of the quark masses and for other lattice spacings. Larger volumes needed at smaller quark masses will require an increase of N ev , the dimension of the LapH subspace discussed in Sec. 3.2, but not the number of Dirac matrix inversions in the stochastic-LapH algorithm for all-to-all quark propagators. Nevertheless, the increasingly severe signal-to-noise problem will likely require more configurations and source times to achieve a similar statistical precision.
This work is part of a larger effort to compute baryon scattering amplitudes on lattice QCD gauge field ensembles at quark masses in the chiral regime m π 300 MeV where effective theories may be applicable. As discussed in Sec. 3.2, the stochastic LapH approach to quark propagation enables considerable re-use of the hadron tensors in multiple multi-hadron correlation functions on the D200 ensemble employed here. Analyses are currently underway to compute the analogous amplitudes for the N Λ − N Σ and N N systems. Hopefully this exploratory computation has sufficient statistical precision to impact chiral effective theories for these baryon-baryon channels as well. Appendix A. Systematic errors from correlation matrix rotation

Acknowledgments
As discussed in Sec. 3.3, the optimized diagonal correlation functions D n (t) are obtained in this work from the GEVP using a single-pivot approach which uses one choice of (t 0 , t d ). The systematic error associated with this approach is estimated for each energy level by fixing the fit range [t min , t max ] and varying the GEVP metric and diagonalization times (t 0 , t d ) defined in Eq. (8), as well as the dimension of the input correlation matrix N op . Taking both GEVP stability and statistical precision into account, the parameters (t 0 , t d ) = (8a, 16a) are found to work well for all energies presented here. As shown in Fig. A.10, the spectrum is rather insensitive to variations in (t 0 , t d ) and N op .

Appendix B. Systematic errors from varying fit forms and time ranges
As discussed in Sec. 3.3, multiple fit ranges and fit forms are compared for every energy level to ensure systematic errors associated with excited state contamination are smaller than the statistical errors. Ultimately, single-exponential fits to the correlator ratios in Eq. (9) are chosen due to their mild sensitivity to t min and good statistical precision. The fit range is chosen to be consistent with the double-exponential t min plateau, defined as the range of t min for which the fitted energy exhibits no statistically significant variation. Most levels are additionally consistent with the single-exponential fit plateau, although as shown in Fig. 2 for m N , these fits may fail to describe correlators with significant excited-state contamination. Plots analogous to the t min -plot in