Elastic Nucleon-Pion scattering amplitudes in the ∆ channel at physical pion mass from Lattice QCD

.


I. INTRODUCTION
The precise treatment of nucleon resonances is still a formidable task in lattice QCD.While Lüscher's method [1][2][3] is the established theoretical basis for connecting observed lattice QCD spectra to scattering amplitudes, thus allowing the investigation of properties of bound states and resonances from first principles, in practice the study of meson-baryon, two-hadron states remains challenging.This is especially so when using simulations with physical values of the light quark mass, which carry increased statistical errors.Nevertheless, the ab initio computation of low-energy elastic pion-nucleon (πN ) scattering from lattice QCD is essential for the study of nucleon interactions, and any such treatment necessarily starts with the lowest-lying meson-baryon resonance, namely the I (J P ) = 3  2 ( 3 2 ) + ∆ P -wave resonance.The ∆ resonance governs nucleonpion, nucleon-photon, and nucleon-neutrino scattering as the dominant channel.Within the lattice QCD formalism, excited state contributions from pion-nucleon scattering states dominate the spectrum in nucleon form factor calculations in a finite volume for gauge ensembles generated with close to physical pion mass [4,5].The interaction of nucleon and pion has been studied by various approaches in lattice QCD in the past, using gauge ensembles generated with heavier-than-physical pion masses.Nucleon-pion scattering amplitudes and the ∆ in particular have been the subject of Refs.[6][7][8][9][10][11][12][13] using the Lüscher method.Refs.[14,15] used an alternative method based on Refs.[16,17].There are also studies using ensembles generated with quark masses at which the ∆ is a stable state [18][19][20][21][22][23].With this work, we extend the lattice calculation of the ∆ from heavier-than-physical to physical pion mass, and explore the application of the Lüscher method to the πN − ∆ channel at the physical point.We estimate the ∆ resonance pole in the P -wave channel as well as the Swave isospin-3/2 scattering length, which experimentally enters the evaluation of the pion-nucleon σ-term using the Roy-Steiner-equations.For this first physical point calculation, we use a single ensemble with two degenerate light quarks, strange, and charm quarks (N f = 2+1+1).The paper is organized as follows: In Sec.II, we describe our application of the Lüscher method, with further details on the implementation of the lattice spectroscopy in Sec.III.Secs.IV and V present our analysis of the lattice data for the πN − ∆ spectrum, and subsequent fits to the finite-volume quantization condition.In addition, we use the threshold expansion of the Lüscher quantization condition to determine the S-wave scattering length.A discussion of the results and concluding remarks are given in Sec.VI.In appendices A, B, and C we include extended figures and tables of our results, as well as some additional details of our analysis.

II. LATTICE QCD FORMALISM FOR PION-NUCLEON SCATTERING
We consider elastic scattering of two particles of nonequal mass and with spin in a hypercubic box of spatial extent L. The constraints on infinite-volume scattering amplitudes from the finite volume lattice spectrum for this setup are given in Lüscher's original work [1][2][3] and a series of extensions, to moving frames [24,25], to non-degenerate particles [26] and to particles with spin [27][28][29][30].While the formalism for multiple coupled two-particle decay channels is also known [31,32], for the purposes of this work we perform a single decay channel analysis, since the expected branching fraction ∆ → πN is almost 100% making this channel a prototype application for lattice QCD.The finite-volume method for three particles has been further developed for spin-less particles [33][34][35][36].For the N ππ system a formalism is not yet available, and we thus only consider the πN interaction here.We will define an upper limit to the spectrum entering our elastic scattering amplitude fits, such that the impact of the three-particle scattering is negligible.
For the pion-nucleon system, even and odd partial waves mix.Moreover, the spin of the nucleon couples to the orbital angular momentum, such that in the finite-volume analysis the most relevant partial wave amplitudes are J = 1/2, containing S-and P -wave, and J = 3/2, with P -and D-wave.The J = 3/2 P -wave amplitude corresponds to the ∆ resonance and is expected to dominate, with a sub-leading contribution from J = 1/2 S-wave, while amplitudes from the corresponding higher ℓ values are suppressed by angular momentum barrier.
To the extent described above, the Lüscher quantization conditions have been given in detail in Refs.[30,37] with the master equation given by det where M ⃗ P ,Λ denotes the reduced Lüscher finite volume matrix for a reference frame with total momentum ⃗ P of the πN system, and for irreducible representation (irrep) Λ of the little group LG( ⃗ P ) ⊆ O D h .Here we use the reduced Lüscher matrices and quantization conditions from Refs.[30,37].The determinant is taken in the linear space of total angular momentum J, J ′ , orbital angular momentum ℓ, ℓ ′ and multiplicity (or occurrence) of the irrep Λ n, n ′ .For ensembles with physical value of the pion mass, the thresholds for elastic nucleon-pion scattering are given by E 2−thr = m N + m π ≈ 1080 MeV and the three-particle threshold E 3−thr = m N + 2m π ≈ 1220 MeV.The correspondingly narrow window in the center of mass energy for elastic πN scattering does not cover the expected resonance region and puts prohibitive cuts on the usable finite-volume lattice spectrum.However, based on experimental observations, the πN scattering amplitude in the I = 3/2 channel is vastly dominated by the elastic two-particle interaction up to E ′ 3−thr = m ∆ + m π as a proxy three-particle threshold, where the ∆ and the pion can propagate on-shell [38].We thus consider lattice energy levels up to √ s ≲ E ′ 3−thr ≈ 1360 MeV.The partial wave amplitudes are parameterized via the associated phase shifts.In particular, for the considered leading phase shifts δ Jℓ , we employ the analytic Breit-Wigner form for the resonant ∆ channel and the leading order effective range expansion with isospin I = 3/2 Swave scattering length a 0 .
The center of mass momentum (q cmf ) of the πN system is given by and the resonance decay width has invariant mass dependence parameterized by

III. CORRELATION MATRIX CONSTRUCTION A. Interpolating operators
To constrain pion-nucleon scattering amplitudes using the Lüscher method, we determine the low-lying finitevolume energy spectrum of the lattice Hamiltonian with isospin I = 3/2.The operator basis for the correlation matrices of two-point functions is constructed from single-and two-hadron timeslice interpolating operators.We consider the case of maximal isospin, i.e.I 3 = +3/2, meaning we use the ∆ ++ , and the proton and charged pion (N + and π + ).The single-hadron, quark model ∆-type interpolating operator reads The two-hadron interpolators are generated by products of nucleon and pion interpolators given by with where u and d are up-and down-quark spinor fields, C = γ 4 γ 2 is the charge-conjugation matrix, and α denotes the spinor index.The total momentum of such interpolating fields is ⃗ P = ⃗ p N + ⃗ p π .
The above operators are projected to irreducible representations of the lattice rotation symmetry groups (O D h for ⃗ P = 0 or little group LG( ⃗ P )) for the rest and moving frames.The values of the total momenta ⃗ P and lattice irreps used are given in Table I, where we also indicate the subduced angular momenta.
The projection to irrep Λ, row r for occurrence n, follows from the group theory master formula, namely for the single hadron operators O ∆ and the πN operators O πN we use where the group element G ∈ LG( ⃗ P ) means either a proper rotation R or a rotation-inversion IR, such that R ⃗ P = − ⃗ P and concatenation with spatial inversion I leaves invariant ⃗ P .Analogously, U (J) (G) denotes the SU (2) spin-J representation matrix of the proper rotation or rotation-inversion group element.The rotation matrix U (1) acting on the momentum vector ⃗ p in Cartesian basis, is denoted as R for simplicity.Moreover, ∆ α,k and N α are Dirac four-spinors and the rotation matrix for the four-component spinors is denoted as U (1/2⊕1/2) .The space inversion operation (⃗ x, t) → (−⃗ x, t) is represented by π + (⃗ x, t) → −π + (−⃗ x, t) for the pseudoscalar pion field, and by ∆ α,k (⃗ x, t) → (γ 4 ) αβ ∆ β,k (−⃗ x, t) for the four-component nucleon and ∆ spinors.
The constructed set of irrep-projected operators O Λ,r,n ∆; α,k ( ⃗ P ), O Λ,r,n N π; α ( ⃗ P , ⃗ p) is still linearly dependent and we extract a basis set by the Gram-Schmidt procedure.The orthogonalization is done with respect to the tensor components α, k and momentum vector ⃗ p.

B. Equivalent moving frames
The explicit subduction coefficients for the operators, which are obtained based on the convention in Ref. [39], pertain to the reference moving frames with total momentum ⃗ P /(2π/L) = (0, 0, 1), (0, 1, 1), (1, 1, 1).We include lattice correlation functions also from all other equivalent moving frames of the same orbit as ⃗ P under discrete rotations.to their reference directions.The first set of two columns correspond to reference moving frame ⃗ P /(2π/L) = (0, 0, 1), the second to (0, 1, 1), and the third set to (1, 1, 1).The superscript i denotes the inverse group element.
Our method to compute correlation functions of singleand two-particle operators efficiently is based on the factorization of the quark flow diagrams [7].We show the diagrams that contribute to the N π−N π two-point function in Fig. 1.The colors used for the different propagators shown in Fig. 1 denote the method used to evaluate them.We denote spinor indices by lower case Greek indices, whereas color indices are lower case Latin indices.
In particular, point-to-all propagators (S) denoted with red lines read, where x i , β, and b are the coordinates, Dirac index, and color index of the source, respectively.Sequential propagators (T ) are denoted with blue lines and are given by, constructed from an inversion of the Dirac operator using a point-to-all propagator as a source vector (also referred to as the sequential source) with support on timeslice t seq and with a vertex given by Dirac structure Γ seq and three-momentum insertion ⃗ p seq .The stochastic sources (η) and propagators (ϕ) are denoted with green single and double lines respectively, where η r α,a (x) are sources with Z 2 × iZ 2 independent and identically distributed (iid) noise of zero mean and unit variance, The expectation value E[.] is taken over the stochastic noise index r.In practice, we use time-slice noise sources, i.e. sources with support on a single timeslice, η (r;t0) α,a (x) = η r α,a (x)δ tx,t0 .The gray, one-end-trick propagators, are constructed using spin-diluted Z 2 × iZ 2 stochastic time-slice sources, appropriately multiplied by a momentum phase, η (r;t0;µ;⃗ p) α,a (x) = η r α,a (x)δ tx,t0 δ α,µ e i⃗ p⃗ x and solution vectors ϕ indicated as ϕ , Analogous to Eq. ( 14), the source components are iid with Z 2 ×iZ 2 noise and have zero mean and unit variance.
The diagrams in Fig. 1 are further split into products of building blocks.These cuts are illustrated by the dashed lines.The quark connected diagrams are split at the source point by virtue of spin-color dilution, and by the stochastic decomposition of unity with the stochastic sources and propagators in Eqs. ( 13) and ( 15), schematically η η † ≈ 1.The factors are two-and three-fold propagator products, which are partially reduced in spin-color space and momentum projected at the pion vertex labeled f 2 in Fig. 1 and the nucleon vertex f 1 .Each type of the four diagrams depicted represents several combinations of contractions, and thus by the factorization multiple contractions benefit from sharing a small number of common factors.
The Wick contractions for ∆ − ∆ and ∆ − N π two-point functions are not found to significantly benefit from this factorization and are therefore calculated from point-toall and point-to-all plus sequential propagators, respectively.
Both the computation of building blocks and their subsequent recombination to full correlation functions is carried out on GPUs using the PLEGMA software package [40].

D. Quark field smearing
We apply Gaussian smearing [41] to the source and sink for all quark propagators.The smearing parameters are N Gauss = 140 steps with weight α Gauss = 1.The gauge links entering the Gaussian smearing kernel are APE-smeared [42] with N APE = 50 steps and weight α APE = 0.4.The parameters α Gauss and N Gauss are tuned in order to approximately give a smearing radius for the nucleon of 0.5 fm.∆ ground state energy in the center-of-mass frame, as a function of the smearing parameters confirmed that these values were appropriate.

E. Gauge ensemble and statistics
We use a gauge ensemble generated by the Extended Twisted Mass Collaboration [43] with two degenerate light quarks with twisted mass parameter tuned to reproduce the pion mass and the strange and charm quarks (N f = 2 + 1 + 1) with masses tuned close to their physical values by matching the kaon meson mass and the ratio of renormalized quark masses m c /m s | MS,µ=2 GeV , respectively.The fermion action is given by twisted mass fermions at maximal twist with the addition of a Sheihkoleslami-Wohlert "clover" term.The gluon action is the Iwasaki gauge action.The parameters of the ensemble, denoted as "cB211.072.64", are collected in Table III.We use N conf = 400 well-separated gauge configurations and on each configuration generate correlation functions We build real symmetric correlation matrices where O i (t) are the correlation functions after projecting to the lattice symmetry group, Eq. ( 10).
For extracting the energy levels from the correlation matrices, we use four methods, which we detail in what follows.

Generalized Eigenvalue Problem (GEVP)
In what we will refer to as the GEVP method, we solve the so-called generalized eigenvalue problem (GEVP) [44,45] for the matrix of correlation functions where v n j (t) is the j th component of the n th eigenvector on time-slice t and λ n (t, t 0 ), the n th eigenvalue of this GEVP, referred to as the principal correlator.In order to obtain energy eigenvalues, we fit the principal correlator to a single-exponential form An example of the analysis using GEVP is shown in Fig. 2, where we plot the effective masses m i eff = log Cii(t)  Cii(t+1) of the diagonal elements of the correlation matrix for the case of the center-of-mass frame H g .The five interpolating operators, that include four momentum combinations of πN and one ∆, are indicated in the legend.Operator: The effective mass from each diagonal correlator in the center-of-mass frame Hg correlation matrix, using a 5×5 basis of interpolating operators listed in the legend.The gray bands show the energy levels obtained by the GEVP analysis.

Prony Generalized Eigenvalue Method (PGEVM)
We apply the Prony Generalized Eigenvalue Method [46,47] directly on the principal correlators obtained via the GEVP method.In PGEVM, we solve the second-level GEVP for the correlation matrix C (2) defined by λ n as While with the GEVP method the stability of the ground state from λ n is tested by a conservative choice of t min , with PGEVM a second, consecutive ground state projection is applied, which is expected to lead to an earlier onset of ground state dominance and thus smaller statistical errors.

Athens Model Independent Analysis Scheme (AMIAS)
In the AMIAS method [48] we perform multi-state fits directly to the correlation matrix.Namely, the spectral decomposition of the correlation functions where the amplitudes, A (jk) n and the energies, E n , are fit parameters and n max is used to truncate the spectral expansion.The probability distribution function (PDF) for the complete set of parameters is , where the normalization factor Z = nmax n=1 dA n dE n e −χ 2 /2 .The estimates for the values of the fit parameters and their uncertainties are then obtained as the expectation values and the standard devi-ations of the corresponding PDF, These integrals are computed using standard Monte Carlo methods.In AMIAS, one investigates the behavior of the distributions of the fit parameters of the lower states of interest as the truncation parameter n max is increased.FIG.3: Results on the energy spectrum for the G1 irrep using the AMIAS method.We show the distributions for the energies when tmin/a=11 (blue curves), 12 (red curves), and 13 (green curves).The correlation matrix is the same as in the 5 × 5 problem used in Fig. 2.
As demonstrated in previous applications of this method [48][49][50], at large values of n max , the additional parameters added, to which χ 2 is insensitive, become irrelevant in the integrals of Eq. ( 21) and thus the distributions of the energies of interest converge, without loss of accuracy.In this way larger fit intervals [t min , t max ] (i.e.smaller t min ) can be probed.An example application of AMIAS is shown in Fig. 3, where the distributions of the energy levels for the case of the G 1 irrep are plotted.As can be seen, using n max = 5 all five energy levels are clearly distinguishable.A similar analysis is carried out for the amplitudes A (ij) n and for all irreps considered, to obtain the mean values and errors of these parameters via Eq.(21).All results in this work quoted as using the AMIAS method used up to n max =4.

Ratio Method
Following Ref. [51] we take the principal correlator from the GEVP analysis, and fit the energy shift with respect to a given non-interacting energy level.In practice this is done by taking the ratio of the principal correlator λ n (t, t 0 ) from Eq. ( 16) to the product of single hadron correlation functions with the appropriate momenta, given by We will refer to this approach as the Ratio method.The leading energy dependence of the ratio in Eq. ( 22) follows as where E ⃗ pN,Λ N ,0 and E ⃗ pπ,Λπ,0 are the ground state energy of the nucleon and pion two-point function in the relevant irreps, respectively.Taking the ratio reduces substantially the excited state contribution in a given principal correlator, and especially for small energy shifts ∆E ⃗ P ,Λ,n (⃗ p N , ⃗ p π ) allows to determine the latter with significantly increased statistical precision.For large t the logarithm of C R (t) will converge to a linear function, with the slope corresponding to the energy shift ∆E An example of such an analysis is shown in Fig. 4, for the case of the H g irrep, which corresponds to the center-ofmass frame with total momentum zero (see Table I).The ratio is applied to the π N case for three different relative momenta ⃗ p N , ⃗ p π between the nucleon and pion.From the energy shift obtained by the slope, we reconstruct the interacting two-hadron energy level by shifting back the energies using the continuum dispersion relation Eq. ( 25) has the added advantage, that high-precision estimates for the nucleon and pion mass can be employed.
Using different single-hadron momenta for pion and nucleon with same total momentum to determine the same energy shift and interacting energy level is part of our systematic error analysis.

Correlation matrix basis selection
For our analysis, we progressively add interpolating fields to the correlation matrix used, selecting the most appropriate basis by checking the stability of the spectrum.In particular, the steps we follow are:  1.We include all the relevant single-nucleon and pion momentum combinations and occurrences and perform a GEVP analysis.
2. Based on the eigenvectors obtained from the full GEVP analysis, we restrict the GEVP to using the interpolating operators that dominate the first few energy levels.
3. Starting from the smaller GEVP of the previous step, we gradually extend the basis.The interpolating operators to be added are chosen by observing the effective mass of their diagonal correlators, and whether they yield higher energy states than already seen in the smaller GEVP.As we add these interpolating operators, we check that the statistical errors of the lower-lying energy spectrum do not deteriorate.
In Fig. 5, we illustrate the basis selection process using an example taken from the GEVP method and the center-of-mass frame irrep H g .What is plotted are the five components of each eigenvectors obtained via the GEVP, which is solved on each time-slice.These correspond to the overlaps of each interpolating field used with the lowest-lying energy state.For a certain choice of eigenvector v n obtained via the GEVP, and interpolating operator number α, the overlap is defined as The eigenvectors are normalized to unity,  We note that this same basis was used for Fig. 2 discussed earlier, where we see that the pion-nucleon interpolating field decreases the ground state energy in the particular channel we consider here.We use this basis to extract energy levels for the analysis that follows.

B. Lattice spectrum and stability test
The four methods detailed above are complementary in the way the excited state contamination is treated.Thus, by comparing the results obtained among them, we can check the robustness of our observed energy levels.The comparison is carried out for all irreps considered, observing the stability of the results as we increase the initial fit-range (t min ).An example is shown in Fig. 6, where the four methods are compared for the specific case of the H g irrep.As can be seen in this plot, the PGEVM and ratio methods provide stable results at smaller t min compared to the standard GEVP.Furthermore, the statistical errors carried by the ratio method confirm our expectation that the correlations between numerator and denominator in Eq. ( 22) help in reducing the statistical fluctuations in the energy shift from the non-interacting energy.This same analysis is carried out for all irreps, with the corresponding plots given in appendix A. By observing the stability of the fitted masses as t min is increased, as well as the χ 2 /d.o.f of the fit, we indicate our selected values for each level and for each method with the bands and in appendix B collect the results for GEVP, PGEVM, and AMIAS in Table VIII and for the ratio method in Table IX, where we also include results for two larger t min values that we use in a model averaging for our final result in Sec.V 1.
In Fig. 7, we collect all the energy levels extracted from all irreps considered, using all four methods.The πN threshold (E 2−thr ) and ππN (E 3−thr ) are shown for comparison at 1080 MeV and 1220 MeV respectively.Furthermore, for each irrep we indicate the permitted noninteracting energy levels that correspond to the lattice volume of the ensemble used.For GEVP and PGEVM, the band, indicating the uncertainty of the energy levels, contains both statistical and a systematic error from the fit range variation.For AMIAS the systematic error from varying the fit range is negligible.For the ratio method we include only statistical errors here.A dedicated discussion of how we obtain the systematic errors of the ratio method is given in Sec.V 1.
As can be seen, the GEVP, PGEVM, and AMIAS methods yield comparable statistical errors, while the ratio method yields smaller statistical errors, as expected given the previous discussion of Fig. 6.In general, we can identify several energy levels that are incompatible with noninteracting energy levels, however the statistical errors carried by our results, combined with the proximity of the ππN threshold and the first non-interacting energy, make clearly identifying the ∆ resonance rather challenging, even with the large statistics of O(10 4 ) used in this work.Our analysis, when compared to the analogous mesonic system of the ρ resonance [52], highlights the increased requirements for extracting resonance parameters of systems that include baryons and when using physical point ensembles.
All data plotted in Fig. 7 are included in Tables VIII and IX of appendix B.
1.  7: The πN interacting two-hadron energy levels obtained by our analysis.For each irrep and total angular momentum J, indicated on the x-axis, we include results using the four methods employed, namely, from left to right, the GEVP method, the PGEVM method, AMIAS, and the ratio method, with the band height indicating our estimated uncertainties as explained in the text.On the left y-axis we indicate the energy levels in physical units, while on the right y-axis in lattice units.The gray dashed and dotted lines spanning the entire x-axis correspond to the πN and ππN thresholds, namely 1080 and 1220 MeV, respectively.The green, thicker dashed lines correspond to the non-interacting πN energy levels permitted for each irrep and for the volume used in this work.

V. SCATTERING PARAMETERS
The scattering amplitude near the resonance is well described by a Breit-Wigner type resonance.We thus parameterize the ∆-resonance phase shift, which via the Lüscher quantization condition then predicts the finite volume energy spectrum.The lattice spectrum we have determined is then fitted to the prediction, and from the minimization of χ 2 we extract the optimal parameters of the resonance given our spectrum results.In practice, we determine the roots of the quantization condition of Eq. ( 1) for the given set of resonance parameters and construct a correlated, non-linear χ 2 , given by where w ij is the covariance matrix between the lattice data for aE i,lat cms and aE j,lat cms estimated using jackknife resampling.To determine the errors of the fit parameters via the jackknife procedure, we perform the minimization of χ 2 in each jackknife sample.In Table V, we collect the results obtained for the scattering parameters when using the energy levels determined via the GEVP, PGEVM, or AMIAS methods.For the ratio method, which we will use to quote our final values, we carry out a more thorough analysis of the errors that we describe later in this manuscript.For the results in Table V, we either restrict to using only P -wave dominated irreps, with partial wave J = 3/2 and ℓ=1, which involves including five levels in the χ 2 minimization, or we use S-and P -wave dominated irreps, with two partial waves (J, ℓ)=(3/2, 1) or (1/2, 0), thus including in total 14 energy levels.For the latter case, we estimate the scattering length from the combined S-and P -wave fits via Eqs.( 2) and (3).In this work we restrict to providing the scattering length only for this channel, leaving the isospin 1/2 case for a future publication.As can be seen from Table V, results when using the three methods are overall compatible and within statistical errors are consistent with the experimental determinations of the ∆ mass and width.However, the statistical errors for the resonance width are large and do not permit a significant comparison with experiment.We, therefore, opt to using the ratio method, presented below, to quote our final results for the resonance parameters.We note TABLE V: Results for the scattering parameters, namely the resonance mass, MR, resonance width, ΓR, and scattering length, Mπa0, using the Lüscher quantization condition and energy levels determined via the GEVP, PGEVM and AMIAS methods.First three rows when using P -wave only and last three rows when using S − P wave dominant irreps.

Method
Breit-Wigner parameters Mπa0 that for the results in Table V, an interpolation method was used to accelerate the minimization of χ 2 in each jackknife bin, described in detail in appendix C.

Results using the ratio method
The most accurately determined energy levels are obtained using the ratio method and are employed to obtain our final values of the resonance parameters.Given these smaller errors, a thorough analysis of the sources of systematic errors is merited, and we, therefore, consider the following in our fits when using the ratio method: 1. We consider two different ranges of center-of-mass energies, namely (a) including only energy levels below N ππ threshold (E 3−thr ), which leads to including 12 energy levels, and (b) energy levels up to ∆ π (E ′ 3−thr ), which leads to including 14 levels.The latter is defined by the onset of the rise of the inelasticity in the J = 3/2, ℓ = 1 pion nucleon scattering channel [38].
2. We explore the partial wave dependence, i.e. we consider fits with energy levels coming only from P -wave dominant irreps, which leads to including 5 levels, and energy levels having also an Swave contribution, which leads to the combinations mentioned in the previous item.We attempted to include higher partial waves but this led to prohibitively large statistical errors on the parameters.
3. The S-wave contribution is entirely parameterized by the scattering length, which can be obtained directly from the ratio of correlators used in the ratio method, as presented in more detail in Sec.V 2 below.We either perform fits using this direct determination of the scattering length or leaving the scattering length free as an additional fit parameter.
4. We use three different fit ranges for the energy levels, i.e. three values of t min , that indicatively span between t min ≃ 0.3 fm and 1 fm.The smallest t min is determined from the onset of the plateau in the ratio method and is that corresponding to the band in Fig. 6 and Figs.11-18.The largest t min is determined from the onset of the plateau in the single exponential fit to the principal correlator from the GEVP.An intermediate t min is also used between these two values.
5. We vary the non-interacting energy level, i.e. that of the denominator in Eqs. ( 22), (23) in the ratio method.We found in our fits, that the energy levels most sensitive to this variation are the ones for the H g irrep.We, therefore, perform a separate analysis for the H g energy levels using the ratio method with pion-nucleon states with three values of back-to-back momenta, namely , as in Fig. 4.
Considering these variations, our analysis yields 45 results, over which we quantify our systematic uncertainty.
The results are tabulated in FIG. 8: The P -wave phase-shift as a function of the invariant mass Ecms = √ s.The error band is determined using jackknife resampling.The points with horizontal errorbars show each fitted energy level included its jackknife errorbar.
The results of the different fits are averaged according to Ref. [53] to derive a combined statistical and systematic error.As an example, in Fig. 8, we show our result for the P -wave phase-shift for one of the 45 fits that has a significant contribution to the model average.The cumulative distribution function (CDF) for the resonance parameters obtained from the 45 different fits is shown in Fig. 9. Our final result for the resonance mass and width obtained through model averaging are M R = 1269 (39) Stat.(45) Total MeV Γ R = 144 (169) Stat.(181) Total MeV, respectively, where the first error is the statistical error and the second is the total, combining statistical and systematic errors.For the individual contributions of the different sources of systematic error, we refer to Table VI.The systematic error for each source of uncertainty (a) is estimated via where O i and p i are averaged over all other systematics other than a and N a is the number of variations for the given systematic error.For the I = 3/2 S-wave scattering length, we obtain from the combined S-and P -wave fit of the lattice spectrum M π a 0 = −0.16(11).
For further illustration, in Fig. 9, we show separately the systematics that originate from varying the fit ranges, from using different non-interacting levels in the ratio method and considering the two center-of-mass energy ranges.The fact that these curves collapse onto the curve that corresponds to the total statistical plus systematic error indicates that our dominant source of error is statistical.

Direct extraction of the scattering length
The scattering amplitude in the S-wave around the π N threshold can be well described by the leading order effective range expansion using a single parameter, the scattering length.As an alternative to the analysis of the previous section, by which the scattering length is extracted from the phase shift, we can extract this quantity from the energy shift ∆E obtained directly from the ratio method.In particular, correlation functions computed with increased statistics, indicated in the first row of Table IV, allow for a high-statistics calculation of the levels in the G 1u irrep.Using the first level of this irrep and the effective range expansion, we determine the scattering length from [1] ∆E FIG.9: The cumulative probability distribution function for the model averaging carried out to obtain the resonance mass MR (top) and the resonance width ΓR (bottom).We distinguish between systematics arising from the energy-level fit range choice (blue triangles), considering energy levels also slightly above N ππ threshold (brown squares) and including all five sources of systematic errors as explained in the text (red small circles).For the resonance mass all three curves fall on top of each other, while for the width the yellow squares are on top of the red circles.Each dark magenta point corresponds to the central value and statistical error of one particular fit with y-axis corresponding to the numerically determined CDF.The vertical band shows the range of MR (top) and ΓR (bottom) between 16% and 84% which we take as our total error.
This extraction of ∆E from the ratio method benefits from the smaller statistical errors associated with the energy shifts.The data obtained for the ratio and the resulting fit are shown in Fig. 10.
A comparison of our result to other work in the literature is provided in Table VII.In particular, we compare our results to a recent lattice calculation using N f = 2 + 1 clover improved fermions at the heavier-thanphysical pion mass of 200 MeV.We also compare to three other determinations, namely from pionic atoms [54][55][56][57]  using the values of the scattering lengths updated in Ref. [58], from unitarized chiral perturbation theory [59], and from a global fit to low-energy pion-nucleon crosssection data [60].Within our quoted statistical uncertainty, our result is consistent with the latter three results.
TABLE VII: Scattering length in the isospin 3/2 pion-nucleon channel.We compare our result, obtained via Eq.( 29), to a lattice calculation using 200 MeV pion mass [6], a calculation using pionic atom data [58], a calculation using unitarized ChPT [59], and a phenomenological determination via global fits to pion-nucleon cross-section data [60].The error in our determination is only statistical.

VI. CONCLUSIONS
Using an ensemble of twisted mass fermions simulated with two degenerate light quarks, and strange and charm quarks with masses tuned to their physical values, we determine the Breit-Wigner resonance parameters of the lowest lying resonance in the π N I(J ) channel.To our knowledge, this is the first such lattice study using physical point simulations.We use a large number of measurements in order to tackle the expected increase of the statistical uncertainties in the meson-baryon correlation functions.To determine the energy levels, we form correlation matrices of one and two-particle correlation functions after an appropriate group-theoretic projection to the relevant lattice irreducible representations.We employ four methods to extract the energy levels from the correlation functions, including the stan-dard GEVP method and variants thereof.While the four methods yield consistent results, the ratio method, where we determine the energy gap between interacting and non-interacting energy levels via appropriate ratios of two-hadron to single-hadron correlation functions, yields considerably smaller statistical errors.
Restricting our analysis to energy levels obtained via the ratio method, we solve the Lüscher quantization condition to extract the resonance mass and width, varying our fits to probe systematic uncertainties.For the scattering length, we use a direct approach, forming a ratio between correlation functions with increased statistics to obtain the energy shift between non-interacting and interacting states explicitly.For the resonance parameters we find the values, for the scattering length in the 3/2 isospin channel, the resonance mass, and the resonance width, respectively.Our result for the scattering length compares well with phenomenological determinations [61] and determinations from ChPT [5,[58][59][60].A recent lattice study using simulations with 200 MeV pions [6] yields a value between 1.4 and 2.7 times larger than our value when taking into account our statistical uncertainties, a factor which is compatible with the ratio of pion masses used.Our result for the resonance mass is compatible with the expected 1230 -1234 MeV values quoted by the PDG [62], while our value for the resonance width has large uncertainties, requiring more lattice input to be determined with significance.As a first application of the Lüscher method to the πN − ∆ channel at the physical point, this calculation was restricted to a single ensemble, which does not allow for a complete assessment of lattice systematic errors, such as cut-off effects.
The current study, and in particular the level of precision obtained given the statistics used, paves the way for future calculations with multiple physical point ensembles allowing for controlled continuum and infinite volume extrapolations.
In Figs.We list the energy levels obtained using the GEVP method, the PGEVM method, and AMIAS in Table VIII.
The fit range given is the optimal one, along with the associated reduced χ 2 .For GEVP and PGEVM, the quoted uncertainty of the energy levels contains both statistical and a systematic error from the fit range variation.For AMIAS the systematic error from varying the fit range is negligible.For the ratio method, the error is purely statistical.
The ratio method energy levels that were included in   IX.The quoted uncertainty is statistical only and the closest non-interacting levels were used to obtain the energy shift.See Sec.V 1 and appendix D for details on how the systematic error is computed for the ratio method, which we quote in our final results.
TABLE X: Results for the scattering parameters, namely the resonance mass, MR, resonance width, ΓR, and scattering length, Mπa0, using the Lüscher quantization condition and energy levels determined from the ratio method.We delineate two of the five fit variations listed in Sec.V 1 via the table subheadings, namely the variation arising from 1.) varying the number of energy levels included and 2.) whether considering only P -or P -and S-wave.For 3.), we give the value of the scattering length a0 when this is left as a free parameter and omit it when it is fixed via the direct evaluation.The first column labels the remaining two variations, namely 4.) that obtained by varying the fit-range in three ways (first index) and 5.) by varying the back-to-back momenta of the noninteracting pion-nucleon operator used in the denominator of the ratio method (second index).

α |v n α | 2
= 1 for all n.From Fig. 5, we see that the first two eigenvectors of the GEVP are largely dominated by the amplitudes from two interpolators, O Hg ∆ and O Hg N π (⃗ p, −⃗ p) with ⃗ p 2 = 1.

FIG. 10 :
FIG.10:The ratio method as applied to the first level of the G1u irrep.With red band we indicate our fit in the range t a ∈[3,17].

Fig. 7
Fig. 7 are listed in TableIX.The quoted uncertainty is statistical only and the closest non-interacting levels were used to obtain the energy shift.See Sec.V 1 and appendix D for details on how the systematic error is computed for the ratio method, which we quote in our final results.

TABLE IV :
Statistics used in this work.We indicate the number of configurations used (N conf ), and per configuration, the number of source points (Nsrc) for point-to-all propagators, the number of stochastic timeslice sources (N stoch ), and the number of stochastic one-end-trick sources (Noet).
IV. SPECTRAL ANALYSISA.Correlator matrices, fits and excited state identification FIG.6: Energy levels for the case of the Hg irrep.The set of points, from left to right, correspond to results from the GEVP method, the PGEVM method, AMIAS, and the ratio method.Dashed lines indicate non-interacting levels for comparison.The red bands correspond to values listed in Table VIII for GEVP, PGEVM, and AMIAS, and to the values with the smallest tmin listed in TableIXfor the ratio method.
Table X of appendix D.

TABLE VI
11,12,13,14,15,16,17, and 18 we plot the extracted energy levels in all irreps for all four methods as a function of the lower fit-range.The figures follow the conventions of Fig.6in the main text.N π energy level fits using the GEVP, PGEVM, AMIAS, and Ratio methods for the irrep G1u.The colored bands correspond to our final selection that is used in Fig.7.FIG.13:Same as in Fig. 11 but for irrep G1.FIG.14: Same as in Fig. 11 but for irrep G2.FIG.15: Same as in Fig. 11 but irrep 2G.FIG.16: Same as in Fig. 11 but for irrep 3G.
h , H g FIG.12: Same as in Fig.11but for irrep Hg.