The elastic $I=3/2$ $p$-wave nucleon-pion scattering amplitude and the $\Delta(1232)$ resonance from $N_{\mathrm{f}}=2+1$ lattice QCD

We present the first direct determination of meson-baryon resonance parameters from a scattering amplitude calculated using lattice QCD. In particular, we calculate the elastic $I=3/2$, $p$-wave nucleon-pion amplitude on a single ensemble of $N_{\mathrm{f}}=2+1$ Wilson-clover fermions with $m_{\pi}=280\mathrm{MeV}$ and $m_{K}=460\mathrm{MeV}$. At these quark masses, the $\Delta(1232)$ resonance pole is found close to the $N-\pi$ threshold and a Breit-Wigner fit to the amplitude gives $g^{\mathrm{BW}}_{\Delta N\pi}=19.0(4.7)$ in agreement with phenomenological determinations.


INTRODUCTION
Accurate and precise predictions of hadron-hadron scattering amplitudes from first principles are desirable for many phenomenological applications. While lattice QCD has been successful in calculating many singlehadron properties, hadron-hadron scattering amplitudes have been more difficult. Since lattice QCD simulations are performed in Euclidean time and finite volume, realtime infinite-volume scattering amplitudes cannot be calculated directly [1]. Instead, the finite volume may be exploited to determine scattering amplitudes using the shift of interacting two-hadron energies from their noninteracting values [2].
However, these calculations have been hampered by the difficulty in evaluating correlation functions with twohadron interpolating operators. Thanks to algorithmic advances in the treatment of all-to-all quark propagators [3,4] and increasing computational resources, lattice QCD studies of scattering amplitudes have undergone substantial recent progress. As reviewed in (e.g.) Ref. [5], many calculations of resonant meson-meson amplitudes have been performed. Elastic meson-meson scattering amplitudes are therefore quickly entering an era of high precision, while first progress has been made on amplitudes with multiple coupled meson-meson scattering channels and on amplitudes coupled to external currents .
The situation with meson-baryon scattering amplitudes is considerably less advanced. There are calculations of non-resonant amplitudes and scattering lengths [36][37][38][39], as well as first steps towards resonant N − π amplitudes [40,41]. Nonetheless, to date published determinations of meson-baryon resonance parameters from amplitudes calculated using lattice QCD are lacking. Unpublished preliminary progress toward a calculation of the amplitude considered in this work was communicated privately in Fig. 17 of Ref. [42].
The ∆(1232) is the lowest-lying baryon resonance, but remains phenomenologically interesting. For instance, as discussed in Ref. [43] nucleon-∆(1232) transition form factors are an important phenomenological input for neutrino-nucleus scattering experiments such as NOνA and DUNE. The scattering amplitude calculated in this work is a required first step in the calculation of such form factors using lattice QCD.
There are many difficulties associated with mesonbaryon scattering amplitudes. First, the exponential degradation in the signal-to-noise ratio is typically worse in correlation functions containing baryon interpolators than in the pure-meson sector. Second, the additional valence quark results in increased computational and storage costs, as well as a proliferation of the necessary Wick contractions. Furthermore, resonant meson-baryon amplitudes require 'annihilation diagrams' which are present in correlation functions between single-baryon and meson-baryon interpolators. Finally, nonzero baryon spin complicates the construction of irreducible mesonbaryon operators [44] and the extraction of scattering amplitudes [45].
Despite these difficulties, we present here the first lattice QCD calculation of the resonant I = 3/2, N − π amplitude in the elastic region on a single ensemble of gauge field configurations with N f = 2 + 1 dynamical flavors of Wilson-clover fermions generated as part of the Coordinated Lattice Simulations (CLS) consortium [46]. Although this ensemble is at unphysically heavy (degenerate) light quark mass corresponding to m π = 280MeV, we observe an analogue of the ∆(1232) resonance close to N − π threshold. Using a variety of moving frames [47,48], we employ finite-volume energies to determine the amplitude at six points in the elastic region. The energy dependence of the resultant amplitude is well-described by a Breit-Wigner resonance shape, yielding m ∆ = 1344(20)MeV and g BW ∆N π = 19.0(4.7) already at our heavier-thanphysical light quark masses.
This letter is organized as follows. We first detail the lattice QCD ensemble employed here, our method for extracting the spectrum, and the subsequent determination of the amplitude. This is followed by the presentation and analysis of the results. Finally, we close with conclusions and an outlook.

LATTICE QCD METHODS
Ensemble details: The gauge field ensemble employed here is from the CLS consortium [46,49], which has generated a large set of N f = 2+1 flavor ensembles at several lattice spacings and quark masses. This single ensemble is detailed in Tab. I and does not have the quark masses set to their physical values but belongs to a quark mass trajectory where 2m +m s is kept fixed as m = m u = m d is lowered to its physical value. Therefore we have both m π > m phys π and m K < m phys K . This ensemble also employs the open temporal boundary conditions of Ref. [51]. In order to ensure a Hermitian matrix of correlation functions, our interpolating operators are always separated from the temporal boundaries by at least t bnd , where m π t bnd = 3.5. Using the zeromomentum single-pion correlation function, which is the most precisely determined correlation function, we have demonstrated that this separation is sufficient to reduce temporal boundary effects below the statistical precision. Temporal boundaries are therefore neglected in all subsequent analysis.
Correlation functions: In order to efficiently treat the all-to-all quark propagators required in two-hadron correlation functions, we employ the stochastic LapH method [4]. While brute-force calculation of the entire quark propagator is intractable, this method projects it onto a low-dimensional subspace spanned by N ev eigenmodes of the stout link-smeared [52] gauge-covariant 3-D Laplace operator. This projection is a form of quark smearing, a common technique used to reduce unwanted excited state contamination in temporal correlation functions.
The stochastic LapH method [4] then introduces stochastic estimators for the smeared-smeared quark propagator Q(x, y) in this subspace spanned by time ('T'), spin ('S'), and Laplacian eigenvector ('L') indices. The variance of these estimators may be improved via dilution [53]. In each index we shall either consider full dilution, denoted 'F', or n uniformly interlaced dilution projectors, denoted 'In'. Furthermore, it is beneficial to employ different dilution schemes for 'fixed' quark lines, where x 0 = y 0 , and for 'relative' quark lines, where x 0 = y 0 . Fixed and relative dilution schemes are denoted by the subscripts 'F' and 'R', respectively. Information on the stochastic LapH implementation is given in Tab. II. This scheme together with the N t0 = 2 source times for our fixed quark lines results in N D = 1152 Dirac matrix inversions per configuration. Using this algorithm, all required Wick contractions are evaluated as described in Ref. [40] while only a single permutation of the stochastic quark line estimates is employed. To increase statistics we average over all irrep rows and a subset of equivalent total momenta P.
Energy calculation: Shifts of the finite-volume N − π energies from their noninteracting values are calculated directly by fitting the ratios [35] to the ansatz R n (t) = Ae −∆Ent . In Eq. 1 C(t) is a correlator matrix in a particular irreducible representation (irrep) and v n (t 0 , t d ) an eigenvector from the generalized eigenvalue problem (GEVP) C π and C N are single-pion and single-nucleon correlation functions (respectively) with momenta equal to those in the noninteracting N − π level closest to E n . While the GEVP enables the extraction of excited states, it is also practically advantageous to enhance ground-state overlap in correlators with significant mixing between the operators and eigenstates. We include one single-site (smeared) ∆ interpolating operator and several nucleon-pion interpolators in the GEVP for each irrep resulting in correlation matrices of dimension N op < ∼ 5. We employ the ground state in each irrep in our analysis, as well as a single precisely determined excited state. With our current level of statistics, other levels in the elastic region have insufficient statistical precision to constrain the amplitude.
Effects due to variation of (t 0 , t d ) and N op are not visible with our current statistical precision. Furthermore, as seen in Fig. 1, we choose fit ranges [t min , t max ] with t min large enough so that the systematic error due to unwanted excited state contamination is smaller than the statistical error. It should be noted that the excited state contamination in R(t) may be non-monotonically decreasing, leading to 'bumps' in the t min -plots shown in Fig. 1. Nonetheless, the overall magnitude of the excited state contamination is considerably smaller than in single-exponential fits to just the numerator of Eq. 1. Furthermore, the chosen t min values lie in the plateau region for individual effective masses of both numerator and denominator, so we are confident that R n (t) behaves asymptotically for t ≥ t min .
Although multihadron correlation functions containing baryons are more computationally intensive than those with just mesons, the overall measurement cost is still dominated by the Dirac matrix inversions, which we perform efficiently using the DFL SAP GCR solver in openQCD [54]. However, the baryon functions defined in Eq. 23 of Ref. [4] are the dominant storage cost.
Amplitude calculation: A variant of the methods of Refs. [2,47] detailed in Refs. [45,48] is applied to relate finite-volume N − π energies to the infinite-volume elastic scattering amplitude. For each total momentum P and irrep Λ, these relations are given as determinant conditions of the form det(K −1 − B (P,Λ) ) = 0, (2) which hold up to exponentially suppressed residual finite volume effects. The determinant is taken over indices corresponding to the total angular momentum J, total orbital angular momentum , total spin S, and an occurrence index n labelling multiple occurrences of the partial wave in the irrep. The (infinite dimensional) matrix B depends on the irrep and encodes the reduced symmetries of the finite volume. It is diagonal in S but (in general) dense in all other indices. Expressions for all required elements of B up to J = 13/2 and = 6 are given in Ref. [45].K is diagonal in J, equal to the identity in n, and is related to the usual K-matrix via S; S where q cm is the center-ofmass momentum.
For this first calculation we only include irreps in which the J η = 3/2 + , p-wave is the lowest contributing partial wave [48], namely (P 2 , Λ) = {(0, H g ), (1, G 2 ), (3, F 1 ), (3, F 2 ), (4, G 2 )}. In addition to ignoring the exponential finite volume effects in Eq. 2, contributions from higher > 1 are expected to be negligible based on threshold angular-momentum suppression. We assess the effect of this truncation to = 1 by performing a fit also including all = 2 contributions, namely the J η = 3/2 − and 5/2 − partial waves. For this fit with the = 2 waves, we additionally include the ground state energy in the (0, H u ) irrep where the 3/2 + wave does not contribute, but both the 3/2 − and 5/2 − are present.

RESULTS
Results for the I = 3/2, p-wave elastic N −π scattering amplitude are presented in Fig. 2, where the (rescaled) real part of the inverse amplitude (q cm /m π ) 3 cot δ 3 2 1 is shown as a function of the center-of-mass energy E cm . This quantity is smooth near the elastic N − π threshold and, unlike the scattering phase shift, can describe both near-threshold resonances and bound states. However it is a highly nonlinear function of E cm , so that conventional horizontal and vertical error bars would be significantly correlated. Instead of these, in Fig. 2 for each energy we display one point for each of the central 68% of bootstrap samples. This therefore gives a visual representation of the 1-σ confidence interval for each point in this two-dimensional plot. Finally, in Fig. 2 E cm = (E cm − m N )/m π is shown on the horizontal axis so that the elastic region is given by 1 <Ẽ cm < 2.
We describe the energy dependence of this amplitude with a Breit-Wigner shape with fit parameters m ∆ /m π and g BW ∆N π . The fit is performed using the method of [45] in which the residuals in the correlated-χ 2 are taken to be where A =K −1 − B (P,Λ) from Eq. 2. We take µ = 1, although fit parameters do not vary outside their statistical errors when going from µ = 1 to µ = 10.
The results of this fit (which neglects > 1 partial waves) are m ∆ m π = 4.738 (47), g BW ∆N π = 19.0(4.7), where the errors are statistical only. While our small number of data points makes fits to other parametrizations difficult, we can attempt to describe this partial wave in a nonresonant manner by truncating the effective range expansion at leading order, yielding the oneparameter fit form This fit gives (m π a 3 2 1 ) −3 = −0.099 (14) with χ 2 /d.o.f. = 2.50, indicating a poorer description of the data compared to the Breit-Wigner form of Eq. 3.
We can also assess the impact of the 3/2 − and 5/2 − dwaves which are present in the irreps with nonzero total momenta. In addition to the six energies included in the previous fits, we add the ground state in the total zero momentum H u channel, where these two waves are the lowest contributing partial waves. Although we only have seven energy levels, we nonetheless perform a fourparameter fit including the leading term in the effective range expansion for each of these additional waves The values for m ∆ and g BW ∆N π are consistent with those obtained from truncating to = 1, confirming our insensitivity to these = 2 waves.
Since m N /m π = 3.732 (56), there is no significant difference between m ∆ and the elastic threshold at E th /m π = 1 + m N /m π . By employing the scale determination of Ref. [50] we obtain a mass in physical units of m ∆ = 1344 (20)MeV, where the error on the scale has been combined in quadrature. It is worth emphasising that the quark masses for this ensemble are tuned to satisfy TrM = 2m l +m s = (TrM ) phys as m l is lowered to its physical value, in contrast with the more standard choice where m s = m phys s for all values of m l . Comparison of g BW ∆N π can be made to experiment using the experimental values m exp ∆ ≈ 1232MeV and Γ exp ≈ 117MeV from Ref. [55] and the relation Γ BW = , where q ∆ is the center-of-mass momentum corresponding to the resonance mass. Such a comparison yields g BW,exp ∆N π ≈ 16.9, in agreement with our result. An alternative convention for the ∆N π-coupling is provided by leading-order chiral effective theory [56], which defines Γ as where E N and E π are the energies of the nucleon and pion, respectively, with momenta equal to q ∆ . Using our calculated values for the resonance parameters and m N /m π gives g LO πN ∆ = 37.1(9.2).

CONCLUSIONS AND OUTLOOK
This work presents the first lattice determination of meson-baryon resonance parameters directly from the scattering amplitude. It builds on demonstrably successful algorithms from the meson-meson sector [35], in particular the stochastic LapH method [4], which reduces computational and storage costs for multihadron correlation functions containing baryons significantly compared to the distillation approach of Ref. [3].
The I = 3/2, p-wave, elastic N − π scattering amplitude is calculated here on a single ensemble of gauge configurations, and thus the usual lattice systematic errors due to finite lattice spacing and unphysical quark masses are not addressed. Furthermore while the magnitude of the exponentially suppressed finite volume effects indicated in Eq. 2 is presumably insignificant, this has not been checked explicitly.
This first analysis also avoids the influence of = 0 partial wave mixing in Eq. 2 by judiciously choosing irreps where this wave is absent and the J = 3/2 p-wave is the leading contribution. Future work will include also irreps where the corresponding s-wave is present, which can be analyzed as described in Ref. [45]. These additional finite volume energies will better constrain the energy dependence of the amplitude and enable a more precise analysis of higher partial wave contributions. Furthermore, this calculation employs only a single permutation of stochastic quark line estimates. Additional 'noise orderings' may significantly improve the statistical precision.
Furthermore, these results will soon be complemented by measurements on other CLS ensembles. This will not only enable a check of the lattice spacing and (exponential) finite volume effects, but also elucidate the quark-mass dependence by using ensembles along the TrM = const. trajectory down to m π < ∼ 200MeV. While we have employed the ground states in each irrep and a single excited state from the (3, F 2 ) irrep here, at lighter pion masses, as m ∆ moves further above the elastic threshold, more excited states will also be included to provide additional points.
This first elastic ∆(1232) → N π calculation may be viewed as a stepping stone in several respects. First, techniques for computing correlation functions, extracting finite volume energies, and analysing determinant conditions may be extended to other resonant meson-baryon systems. Such systems of interest may present additional complications like coupled scattering channels, as are present when studying the Λ(1405) resonance. Finally, building upon this calculation of ∆(1232) → N π will ultimately enable lattice QCD calculations of ∆(1232) transition form-factors, for which the theoretical foundations can be found in Ref. [60].
We acknowledge helpful discussions with Harvey Meyer and Daniel Mohler. This work was supported by the U.S. National Science Foundation under Award No. PHY-1613449. Some of our computations are performed using the CHROMA software suite [61]. Computing facilities were provided by the Danish e-Infrastructure Cooperation (DeIC) National HPC Centre at the University of Southern Denmark. We are grateful to our CLS colleagues for sharing the gauge field configurations on which this work is based. We acknowledge PRACE for awarding us access to resource FERMI based in Italy at CINECA, Bologna and to resource Super-MUC based in Germany at LRZ, Munich. Furthermore, this work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s384. We are grateful for the support received by the computer centers. The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS share of the supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC). GCS is the alliance of the three national supercomputing centres HLRS (Universität Stuttgart), JSC (Forschungszentrum Jülich), and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF).