Searching for the first radial excitation of the $\Delta(1232)$ in lattice QCD

We present a lattice QCD analysis of the $ \Delta $-baryon spectrum, with the goal of finding the position of the $ 2s $ radial excitation of the $ \Delta(1232) $ ground state. Using smeared three-quark operators in a correlation matrix analysis, we report masses for the ground, first and second excited states of the $ J^P = 3/2^+ $ spectrum across a broad range of $ m_\pi^2 $. We identify the lowest lying state as being a $ 1s $ state, consistent with the well known $ \Delta(1232) $. The first excitation is identified as a $ 2s $ state, but is found to have a mass of approximately 2.15~GeV on our $ \sim3 $ fm lattice, which does not appear to be associated with the $ \Delta(1600) $ resonance in a significant manner. We also report on the spin-$ 1/2 $ and odd-parity states accessible via our methods. The large excitation energies of the radial excitations provide a potential resolution to the long-standing missing baryon resonances problem.


I. INTRODUCTION
Short-lived hadronic states, otherwise known as resonances, present various challenges in our understanding of the strong force.While the quark model performs well for describing long-lived ground state mesons and baryons [1,2], the complex dynamics of resonances and the possibility of several chains of decays makes modelling and analysis of these and their excited states non-trivial.The foremost example is the N * (1440) Roper resonance [3][4][5], with similar issues associated with the ∆(1600) resonance [6,7].A recent description, arising through Dyson-Schwinger methods, is that the first excitation in both of these channels is a radial excitation of the ground state [5,6].
The only currently tenable method for resolving the properties of QCD at hadronic energy scales is lattice QCD.Lattice QCD is well established as a non-perturbative framework within which many physical quantities can be extracted such as particle masses, decay constants, form factors, etc.The systematics to do with the discretisation of spacetime and using heavier-than-physical quark masses are now readily controlled and understood, leading to subpercent level errors [8].
Light baryon resonances are of particular interest in the context of lattice QCD studies for several reasons.Firstly, as energy increases beyond relevant thresholds, state of the art calculations need to consider contributions from a number of open many-body channels, as in Refs.[9][10][11].This necessitates a robust basis of interpolators that couple to these scattering states, and so provide a spectrum of energies which can be connected to infinite-volume resonances by methods such as that of Lüscher [12][13][14].
One such approach, known as Hamiltonian Effective Field Theory (HEFT), has been implemented in the case of the Roper resonance, leading to strong evidence that it is not a simple radial excitation of the ground state nucleon [15][16][17][18][19].This approach draws on lattice QCD finite-volume energies across a broad range of m 2 π , and so, in anticipation of applying this formalism to extract the nature of the ∆(1600) resonance, the current work aims to provide a similar foundation of masses in the ∆ regime.
In studying the Roper, it was found that a basis of smeared three-quark interpolating fields provides a method of classifying the associated eigenstates in terms of radial excitations [20], based on the number of nodes in the wave function [21,22].We apply the same approach here in the ∆ spectrum, and seek to identify the approximate position of the first radial excitation of the ∆(1232).This will provide some initial insight into the nature of the ∆(1600), and whether it is indeed the radial excitation of the ground state ∆.
In constructing a basis of three-quark interpolating fields, one finds that various combinations of total angular momentum and parity can be projected out, and so it is of interest to also explore beyond the ∆(3/2 + ) and consider the 3/2 − , 1/2 + and 1/2 − states in the ∆ spectrum.Resonances of interest, as well as their PDG classification in terms of how well studied they are, are given in Table I.
The paper is organised as follows.Section II describes the methodology used to compute the ∆ spectrum.We outline the variational method used for extracting excited states in the spectrum as well as a spin projection technique.Finally, we present techniques used for identifying radial excitations in the spectrum based on nodes generated in the superposition of interpolating fields.Our results for the J P = 3/2 + spectrum are discussed in detail in Section III.While our focus is on the 3/2 + spectrum, Section IV contains results for the 3/2 − and 1/2 ± states that we are also able to extract using our methods.These act to disclose the role of wave function nodes in the excitations, with additional comparison available from the work of other recent studies.

II. LATTICE QCD CALCULATIONS
In this work we use well established tools for extracting baryon spectra within the framework of lattice QCD.This involves constructing interpolating fields which couple to single-particle, three-quark dominated states on the lattice with quantum numbers of the ∆-baryon, and constructing 2-point correlation functions from these.These correlation functions can then be spin-and parity-projected to isolate different states sharing the same quark flavours, whose masses can then be extracted using the variational method [20,24,25].The details of each of these steps will be explored in the subsections below.

A. Baryon Spectroscopy
The masses of states in the energy spectrum of a particular baryon species can be simply derived from 2-point correlation functions.At the hadronic level, a general 2point correlation function is defined as The indices i, j refer to choices of interpolating field χ, with Greek indices µ, ν labelling the Lorentz degrees of freedom.
T indicates time ordering of the operators.Putting all this together, the operator χj ν (0) acts on the vacuum to create states at the space-time point (0, 0).These states then propagate through Euclidean time t before being annihilated at a new space-time point x = (t, x) by the operator χ i µ (x).The Fourier transform projects the baryon momentum p.
For the case of a particular baryon species, we have the complete set of states where B labels the baryon species with momenta p ′ and spins s.This allows us to obtain In order to have meaningful overlap with this complete set of baryon states, the form of the interpolating fields should be chosen so that they share the same quantum numbers and symmetry properties as the states we wish to extract.
We then use the translational operator to rewrite χ i µ (x) as This allows us to simplify the expression for the correlation function to At this point it becomes important to distinguish the states allowed for the system in question.Given our choice of three-quark operator for the interpolating field coupling to the ∆ ++ we need to consider both spin-3/2 and spin-1/2 contributions.
Thus, for our case of interest we have the following available matrix elements, with the associated decomposition into Dirac spinors u(p, s) for the spin-1/2 components, and Rarita-Schwinger spinors u µ (p, s) for the spin-3/2 components: Here, the λ, α and β factors are various couplings and M is the mass of the state with relativistic energy E. We note that as both the Rarita-Schwinger spinor and the interpolating field transform under parity as a pseudovector, Eq. ( 9) includes a γ 5 next to the Dirac spinor.Including the corresponding expressions for the χ matrix elements from Eq. ( 5), one can use the Dirac spinor spin-sum relation and the analogous Rarita-Schwinger spin-sum relation to compute the sum over s and write the contributions from each of the required products of matrix elements.At p = 0, so that E B = M B , and with µ = ν = n = 1, 2, 3, we have where we've relabelled β 1/2 ± → λ 1/2 ± for simplicity.The notation λ is used to denote the overlap of the interpolator at the source.
Hence we can see that the even-parity states exist in the (1, 1) and (2, 2) Dirac components of the correlation function, while the odd-parity states live in the (3, 3) and (4, 4) Dirac components.These are then readily accessed by applying the appropriate projection operators and taking the trace over Dirac indices.
In particular, we note that the parity information of our correlation functions is governed by the upper and lower components of the spinors describing the baryon states.This will be of importance in Section IV in discussing the identification of radial excitations for odd-parity states.
One can then define the parity projected correlator as a sum over the spatial Lorentz indices The result is a series of decaying exponentials governed by the mass of the baryon states where λ iB ± and λ jB ± are coupling strengths between the interpolating fields χ i µ and χ j ν and the parity projected baryon states B ± , up to some overall constant factors.We note that these coupling strengths can be taken to be real by considering both the original gauge-field links and their complex conjugates, weighted equally in the ensemble average [20,26].
The parity projected correlation function still contains various spin states (such as spin-1/2 and spin-3/2 in our case) as well as a tower of excited states.In order to extract, say, the ground state mass of a particular baryon state, one takes the long-time limit in which all the excited states have decayed off.Explicitly, where the λ ± i0 and λ ± j0 are couplings of baryon interpolators at the source and sink to the lowest lying state.In order to isolate spin states one can employ further projection techniques as discussed in the following section.A discussion of how one can access excited energy states is presented in Section II D.

B. Interpolating Operators and Spin Projection
To simulate ∆ baryons on the lattice, we start from the well established form for the ∆ ++ [27, 28] The superscript T denotes the transpose, a, b, c are colour indices and u(x) is the Dirac spinor field for the u-quark.
Note that this interpolating field couples to both spin-1/2 and spin-3/2 states (by coupling of the spins of three spin-1/2 quarks).Hence we require additional spin projection to distinguish the relevant contributions to the overall correlation function.
The general form of a spin-s projected correlator G s in our formalism is G µσ g σλ P s λν , where P s µν are spin-s projectors for s = 1/2 or 3/2.The form of these is given in Appendix A. Once a correlator is spin-projected, it can then be subsituted directly into Eq.( 18) with the calculation proceeding as shown.

C. Source and Sink Smearing
In order to improve the overlap of the interpolating fields with the states of interest, we apply Gaussian smearing to the spatial components of the interpolating fields.The general procedure is to take some fermion field ψ i (t, x) and iteratively apply a smearing function F (x, x ′ ).Explicitly, this takes the form where the smearing function is We take the smearing parameter to be α = 0.7 in our calculations.The use of repeated applications of the smearing function controls the width of our source by gradually smearing out an initial point source.This is visualised in Fig. 1 where, for 16, 35, 100 and 200 smearing sweeps, the amplitude of the distribution for U µ (x) = I is plotted on the x-y plane with the third spatial dimension fixed at the centre of the source.These choices of smearing sweeps are motivated by considering the condition number of the correlation matrix, as in Ref. [20].

D. Variational Method
Recall that in Section II A we showed how one can readily obtain the mass of baryon ground states by simply taking the leading order contribution to the correlation function as in Eq. (20).In order to isolate states of higher energies, a more nuanced approach is required, since these states are at sub-leading order in the exponential series.We make use of the well-known variational analysis method [20,24] in order to extract the ground state mass as well as the excited state masses.
To extract N states, we require N interpolating fields (one for each state).We generate our interpolating fields by using the same base field for the ∆ + baryon, obtained from Eq. ( 21) and applying different levels of smearing consisting of 16, 35, 100 or 200 sweeps of smearing.These smearing levels are chosen both for their coupling to the states of interest and to produce a well defined matrix of correlation functions after applying these smearing levels to both the source and sink [15,20].This correlation matrix is written as Here the λ α i and λ α j are essentially the same as seen before in Section II A, though we now use the α index to distinguish between energy states.In other words, they are couplings of the smeared interpolators χ i µ and χ j ν at the sink and source, respectively, to the various energy eigenstates α = 0, . . ., N − 1. m α is the mass of the state α.
From here, we now aim to construct linear combinations of our smeared interpolating fields to cleanly isolate the N states in the baryon spectrum.Labelling these baryon states |B α ⟩, we thus wish to construct the superpositions such that where u µ (α, p, s) is a Rarita-Schwinger spin vector.Here, the z α and z α are the couplings of the superpositions ϕ α µ and ϕ α ν to the state |B α ⟩.The u α i and v α i are simply the weights for the basis of smeared interpolating fields.
At this point, we construct an eigenvalue problem to solve for both u α and v α .Noting that since G ij (t) is real and symmetric, G ij (t) = G ji (t), in the ensemble average we introduce an improved unbiased estimator of the correlation matrix [G ij (t) + G ji (t)]/2.This provides us with a correlation matrix which is symmetric, so we can simultaneously compute u α and v α as discussed below.
Multiplying Eq. ( 25) on the right by u α j we obtain Then, since the exponential is the only time dependent part of the correlation function, we can form a recurrence relation at some time after source insertion by introducing the variational parameters t 0 and ∆t: Then, multiplying on the left by the inverse [G ij (t 0 )] −1 and suppressing the indices i and j gives which we recognise as an eigenvalue equation for the vector in interpolator space u α .Similarly, premultiplying Eq. ( 25) by v α i we get from which we arrive at our second eigenvalue equation (this time for v α ): Both Eq. ( 31) and (33) need to be solved simultaneously for each given pair of variational parameters t 0 and ∆t, and we do so using a generalised eigenvalue problem solver.Solving for these eigenvectors automatically gives us the weights for the superpositions of interpolating fields, by construction.
Finally, the eigenstate and parity projected correlation function is then taken to be The eigenvectors are essentially used to isolate particular states in the baryon spectrum, exactly as we set out to do.
We then construct the effective mass function δt is typically taken to be small and is set independently of the variational parameters.We take δt = 2 in our calculations.
It is also worth noting that the eigenvectors u α and v α are equal since G ij (t) is a real symmetric matrix.From here on, we will refer to the u α vector, for simplicity.
The effective mass defined in Eq. ( 35) can be computed for various discrete values of t and then plotted as a function of Euclidean time.As usual, one looks for time intervals over which the effective mass plot plateaus, indicating that all contamination from unwanted excited states has decayed away in the exponential series.We then perform a covariance matrix analysis of the χ 2 per degree of freedom [29] to determine the most suitable time intervals to fit when obtaining our final masses.The details of this procedure for calculating our fit uncertainties is given in Appendix B.
With a basis of 4 smeared interpolating fields, we have access to 4 energy eigenstates in principle.However, while the third excited state is available to us, it is susceptible to excited state contamination.In some cases the signalto-noise ratio degrades too rapidly for us to reliably report results, and these are omitted.Hence we only report the masses of the ground, first and second excited states.
As a final note, the principles highlighted within this section are fully realised only when one includes a complete set of interpolating fields effective at isolating all the states within the spectrum.In principle this needs to include multi-particle scattering states as highlighted in, for example, Ref. [30].However, our aim is more modest; to determine the approximate position of the 2s excitation of the ∆(1232) in the spectrum.Our formalism is suited to exciting these single-particle states and in particular iden-  tifying their radial node structure.As such, Euclidean time evolution is important in suppressing contamination from nearby states.In fitting our effective mass to a plateau, we enhance single-state dominance by monitoring the χ 2 per degree of freedom (χ 2 /dof).We enforce an upper limit of χ 2 /dof ≤ 1.2 for our fits.As we will see, the 2s excitation sits at approximately 2.2 GeV presenting a formidable challenge to future calculations aiming to include all the multiparticle interpolating fields required below this energy.While our results will miss the subtle energy shifts associated with avoided level crossings of nearby scattering states, knowledge of these subtleties is not required to make progress in understanding the gross structure of ∆-baryon resonances.

E. Simulation Details
Our calculations are based on the PACS-CS Collaboration's gauge-field ensembles [31], made available via the ILDG.These configurations make use of the Iwasaki gauge action and an O(a)-improved Wilson quark action with Clover coefficient c SW = 1.715, in full 2 + 1 flavour dynamical QCD.The simulations use β = 1.90 on a 32 3 × 64 lattice.
The masses of the up and down quarks are taken to be degenerate, and range over 5 values down to near the physical point.Meanwhile, the strange quark mass is held fixed.All masses are described by the Hopping parameter expansion with the Sommer parameter r 0 = 0.04921(64) fm used to set the physical scale, and the characteristics of each ensemble can be seen in Table II.The source position is held fixed at t 0 = 16, one-quarter of the lattice extent in the time direction.This allows for sufficient Euclidean time evolution of states to reduce excited state contamination, while reducing the impact of possible backwards propagating states reflected from the fixed boundary condition in the time direction.
For each mass ensemble, we apply a number of circular shifts to the gauge-fields, effectively allowing the source location to roam around the lattice N src times.We then simply take an average over the resulting N src sets of correlators for each configuration, allowing us to reduce the impact of the signal-to-noise ratio problem for baryons [32][33][34][35].Jackknife uncertainties are calculated in a second-order single-configuration elimination calculation.

III. RESULTS:
We obtain results for each of the J P combinations 1/2 + , 3/2 + , 1/2 − , 3/2 − , by simply applying the corresponding pairs of spin and parity projection operators to the correlation functions of the correlation matrix.
A wide range of variational parameters was considered including t 0 = 17, 18, 19, 20, 21 and ∆t = 2, 3, 4. While early values of t 0 and ∆t produce smaller uncertainties, better plateaus are obtained for larger values.However, large time considerations suffer from both exponential suppression of the excited-state information and the onset of large uncertainties.In general, we find that a choice of t 0 = 18 and ∆t = 2 produces the most stable plateau fits across all 5 quark masses.This allows us to eliminate the variational parameters as having contributed any systematic uncertainty to the calculations as we compare across quark masses.We first present and discuss the results for the effective masses in the 3/2 + channel, and follow this with identification of the radial excitations.

A. Spectrum
The results of our computations are given in Table III as well as in Fig. 2 where data from other recent studies are shown for comparison [31,36,37].Note that the studies labelled HSC and Khan et al. use different sets of interpolating fields and so are able to extract more states.We compare with the other studies to illustrate the systematic uncertainties associated with the calculations presented.All these simulation results producing large excitation energies necessarily do not include the multitude of multiparticle interpolatiing fields that can participate at these energies.However, the three-quark interpolators considered are constructed using diverse techniques and agreement among the various collaborations brings credence to the results presented.Still, a variety of volumes are combined in these comparison plots and as such the avoided level crossings of the volume-dependent energy eigenstates will differ.These more subtle shifts are to be expected and can be used to inform the level of systematic uncertainty in the results presented.As our goal is to determine the approximate position of the 2s excitation of the ∆(1232), the node information provided in this study is key.
Starting with the lowest-lying state, we see that our results are consistent with the original PACS-CS results (offset slightly from ours in Fig. 2 for simplicity) and our lightest ensemble point improves on the original PACS-CS result in terms of the associated uncertainty.This can be attributed to our having averaged over 64 source locations as compared to their 4 source locations.The other lattice results shown in Fig. 2 reflect similar values of the ground state mass.As expected, there is a nonlinear dependence in the light quark mass region of low m 2 π [38][39][40].Moving on to the first excitation results, we see resonable agreement between our results and the nearby data from other studies employing three-quark interpolators.We observe several features of note.Firstly, our excited state appears almost constant in pion mass, and secondly it lies well above the region that one would typically associate with the ∆(1600).We do note however that the analysis of Khan et al. [36] identifies lower-lying states with a large hybrid operator component (highlighted in Fig. 2) and therefore these are not associated with the three-quark excitation of our calculation.Their next states in the spectrum are dominated by three-quark operators and compare favourably with our results.In summary, it is reassuring that despite differences in methodology, lattice volumes and scale setting, our results are in approximate agreement with those of the other collaborations.
The first low-lying excitations reported by the Hadron Spectrum Collaboration (HSC) [37] next to their lowestlying states are not observed in our analysis.A forthcoming analysis [41] using Hamiltonian effective field theory to bring πN scattering data in the ∆-resonance channel to the finite-volume of the lattice indicates that the only finitevolume states in this energy region on a 2 fm lattice are the lowest-lying states.Thus, these states are not associated with the ∆(1600) but rather should be degenerate with the other lowest-lying HSC states presented.We defer further discussion on these particular HSC states to Ref. [41].
These considerations lead to a quite profound qualitative result; implementing smeared three-quark operators for describing the states in the ∆(3/2 + ) spectrum on the lattice appears insufficient for exciting states in the ∆(1600) resonance region.The question arises then; is the ground state a simple 1s state as one expects, and is the first excitation that we have calculated near ∼ 2.15 GeV the corresponding 2s excitation?This would give some tension with the long-held stance that the ∆(1600) is the radial excitation of the ∆(1232).The identification of radial excitations is explored below.

B. Radial Excitations
With the lattice masses calculated, we now identify these in terms of radial excitations.Our approach to this is based on counting the number of nodes in the wave function, as is discussed in detail in [20][21][22].We present here a brief summary of the techniques involved before discussing our results.
In solving the GEVP of the variational method, we build up a description of the wave function for each state in the spectrum as a superposition of Gaussian distributions.An important feature of this is the possible presence of nodes in the wave function.
Consider, for example, the case of a narrow Gaussian of positive signature, superposed with a broader Gaussian of negative signature.As we consider points further away from the centre of the distribution, the positive Gaussian will diminish and the negative Gaussian will begin to dominate the superposition.As a result, there will be some radius at which the superposition crosses through zero.The associated probability amplitude will thus possess a single radial node, corresponding to a 2s state.For a superposition with additional Gaussians which alternate in signature as the widths increase, one can produce more nodes in the wave function.
These techniques of superposing smeared quark sources and sinks through a correlation matrix analysis were used to study the wave functions of even-parity nucleon excitations in Refs.[21,22].There, the linear combination of smeared sources defined in the GEVP analysis was used at the source to excite the energy eigenstates.After fixing the gauge fields to Landau gauge, the quark fields in the sink operator were given individual spatial positions such that they could be pulled apart.For example, fixing the two u quarks of the proton at the centre of the lattice, the wave function of the d-quark was resolved by studying the strength of the baryon annihilation matrix element as a function of the d-quark position.Upon observing nodes in the d-quark wave function, it was noted that this information was already apparent in the correlation matrix eigenvectors for the energy eigenstates, manifested by sign alternating eigenvector components as the Gaussian widths of the sources increase.It is this method that is used here to identify nodes.
We can view our GEVP procedure as constructing the superposition of Gaussians which best describes the wave function of the states in the spectrum, that is, solving the coefficients u α i in such that ϕ α is a good approximation to the wave function for energy level α.The values for the various u α i are displayed in Fig. 3.
To reiterate then, we start by considering the u 16 components (the weights for the narrowest Gaussian) and observe whether the weight changes sign as we move on to the u 35 .If it does, we would associate this with a node.We repeat this process as we continue onto the u 100 and u 200 .
For the eigenvector components to have a physical interpretation, a normalisation prescription is required for the baryon interpolating fields of various smearing levels.We proceed by introducing a real scaling factor, s i , for each baryon interpolating field, χ i → s i χ i , and normalising the diagonal spin-parity projected correlation functions of Eq. ( 18) to one at one time slice after the source, i.e.G ii (t − t 0 = 1, p = 0) = 1 for each i = 16, 35, 100, 200.This normalisation condition for the diagonal correlators having the same smeared source and sink provides a constraint for the determination of the scaling factor s i .With the χ i suitably normalised, the coefficients u α i may then be applied to normalised Gaussians to gain insight into the nature of the interpolating field wave function.
We can visualise these superpositions by simply adding up normalised Gaussian distributions weighted by u α i for each state α.Taking a 2-dimensional slice of our lattice in the x − y plane, the square of this superposition can then be plotted as a probability amplitude as shown in Fig. 4 for the m π = 413 MeV ensemble.Indeed, one can see the presence of nodes in both Figs. 3 and 4.
Consider Fig. 3 first, and in particular the middle quark mass results for the u vector components.Small eigenvector components will introduce ripples in the wave functions but are not large enough to change the sign within the lattice volume.In this light, State 0 has all of its large components being non-negative, so there is no node visible in the probability amplitude illustrated in the left-hand plot of Fig. 4 2. J P = 3/2 + ∆ spectrum with CSSM lattice results from this investigation along with other comparable studies in the recent literature: PACS-CS [31], HSC [37], and Khan et al. [36].In reporting results from Khan et al. we have differentiated states dominated by conventional three-quark operators and those having significant hybrid components.We omit the results for State 2, these are given in Table III.with u 100 > 0 and u 200 < 0. This case is effectively the same as the toy example above where a broad Gaussian is being subtracted from a narrow Gaussian.The result is that we expect a single node in the probability amplitude and this is reflected in the right-hand plot of Fig. 4. Turning our attention to State 2 in Fig. 3, we note the u 16 values are small relative to the u 35 values and do not create a node.However sign oscillations from u 35 to u 100 to u 200 , all of which have significant magnitudes, induce 2 nodes in the interpolator wave function.Similarly, State 3 has 3 nodes.
At this point we have identified the masses of three states in the spectrum, and have further classified these as being a 1s, 2s and 3s state.Interestingly, while the 1s state aligns well with the ∆(1232), the 2s state has a mass around ∼ 2.15 GeV so cannot be associated in a significant manner with the ∆(1600).Thus it appears that this state is associated with 2-and 3-body scattering channels.The use of Hamiltonian Effective Field Theory (HEFT) as a Lüscher adjacent method will allow us to provide insight in a forthcoming paper into which channels may be useful in future lattice studies for examining the ∆(1600) resonance.
IV. ADDITIONAL RESULTS: While the spin-3/2 even-parity case is of most interest in the current study, we also report results for the spin-1/2 states as well as odd-parity states for both spins.The spectra for these cases are given in Figs. 5, 6 and 7, with numerical values for our fits given in Tables IV, V and VI.While we are not able to extract excited states for all of these cases, we can observe a distinct trend in the ground state masses: • the 3/2 + has the lowest mass.
• the 3/2 − and 1/2 − are a little higher in energy, with similar masses.
• the 1/2 + has a mass well above these states.This is consistent with previous CSSM results using quenched QCD [42], and has good alignment with the resonances observed in experiment [23].

A. Spectra
Considering the resonance positions in Table I, we note alignment between our results and those from the PDG, in that the ground state 1/2 − and 3/2 − states are nearby to one another, and likewise for the corresponding excited states.These appear to correspond approximately to the ∆(1620)(1/2 − ) and ∆(1700)(3/2 − ), and the ∆(1900)(1/2 − ) and ∆(1940)(3/2 − ), respectively.This suggests these resonances may be single-particle dominated states dressed by their decay channels.
Furthermore, there is excellent agreement across all our lowest-lying masses with those reported by HSC, as well as Khan et al.Our 3/2 − results are competitive with those of the other studies, and show fair agreement even at the first excitation level.
The 1/2 + case presents the greatest challenge to our approach.Even the lowest-lying state could not be extracted at the lightest pion mass, and no excited states were able to be extracted at any pion mass.This is due to the initial interpolating field we use being constructed for coupling primarily to spin-3/2 states.Once we project both the spin and parity, there is little signal left for the spin-1/2 case.This is particularly problematic in the 1/2 + case since the lowest-lying state is expected to lie much higher in energy than the other cases.
While the 1/2 + results are clearly not competitive with those of the other groups, there is approximate agreement with their findings.A more suitable choice of interpolator basis is required for further study of this channel.

B. Radial Excitations
Similarly to the 3/2 + results, we can also consider the u vector components for each of these states, and visualise the resulting descriptions of the wave functions.However, there is a subtle point worth mentioning here, specifically how the odd-parity properties are accessed.In particular, as was outlined in Section II A, the parity nature of our correlators is dictated by the presence of a γ 5 matrix in the overlap of interpolating fields and baryon states.We thus expect that odd-parity character should be contained in only the spinor part of the wave function.This means that even for the odd-parity baryons, we should expect a radially symmetric description of the spatial wave function, just as was obtained in Fig. 4.
The visualisations of the interpolator wave functions for the 3/2 − and 1/2 − states are shown in Figs. 8 and 9.Note that while we are able to cleanly extract all u vector components for all 4 energy states in each spin-parity pairing, we only show the wavef function plots for which we were also able to extract masses.The 1/2 + case can also be assigned a Gaussian superposition for each energy state, but since we were only able to extract a ground state energy in this case, we omit the plots of the associated Gaussian distributions.
Once again, we see that the first excited state involves a node in its structure.For the 3/2 − case, our results suggest that the ∆(1940)(3/2 − ) may not be a simple excitation of the ∆(1700)(3/2 − ), since it lies well below the excitation energy determined herein.This resonance may be associated with one or more scattering channels in a non-trivial way, as has been suggested for the ∆(1600)(3/2 + ) [43].
For the 1/2 − case, there is doubt whether the ∆(1900)(1/2 − ) or even the ∆(2150)(1/2 − ) are simple radial excitations of the ∆(1620)(1/2 − ), owing to a substantial energy gap between these and the observed 2s excitation energy on the lattice.Again, a coupled-channel analysis approach may shed some light on the structure of these resonances.

V. CONCLUSION
We have presented an analysis of the low-lying ∆-baryons using techniques within the framework of lattice QCD.Our focus is on identifying the position of the 2s excitation of the ∆(1232)(3/2 + ) in the spectrum.Our technique of superposing smeared fermion fields in constructing our baryon sources lends itself to identifying the nodes of wave functions, thereby identifying radial excitations of states obtained by the variational method.This has allowed us to conclude the following for the energy eigenstates observed in our calculations: 1.The lowest-lying state having an energy in the ∆(1232)(3/2 + ) resonance region is identified as a 1s state.
2. The first excitation has an energy of 2.14(5) GeV in the lightest mass ensemble, with only weak dependence on quark mass evident across the considered ensembles.This state is identified as a 2s state, and thus the radial excitation of the ∆(1232)(3/2 + ).There is a substantial energy gap between this state and the ∆(1600)(3/2 + ) and it is unlikely to be strongly associated with this resonance.We note that this situation is similar to that observed in Refs.[18,21] where the 2s excitation was resolved at 1.90(6) GeV, well above the position of the Roper resonance at 1.44 GeV.
3. The second excitation is identified as a 3s state with two nodes in the wave function.The energy of this state on our 3 fm lattice is large at 3.10(17) GeV.
When compared with the 1s and 2s energies of 1.26(2) and 2.14( 5) GeV respectively, one observes a mass splitting of ∼ 0.9 GeV between the radial excitations, very large compared to conventional quark models.
We conclude that the missing baryon resonances problem of the last three decades follows from incorrectly tuning quark model parameters to reproduce the N (1440) Roper resonance and the ∆(1600) resonance as 2s radial excitations.We now know the 2s radial excitations sit at 1.9 and 2.1 GeV for the N and ∆ respectively.Thus there is an important opportunity for the development of modern quark models whose spectra are compatible with this modern understanding of quantum field theory.
Aside from the main results for the 3/2 + channel, we have also reported results for the 3/2 − , 1/2 + and 1/2 − spectra.Here we resolved a node in the first excitations of the 1/2 − and 3/2 − states illustrating the internal structure of these states for the first time.Our spectra are qualitatively similar to other collaborations pursuing these high-energy excitations.
Since a simple three-quark dominated model appears to be insufficient to describe the structure of the ∆(1600), attention must turn to interactions in multi-particle mesonbaryon channels.It is expected that the dominant scattering channels of πN and π∆ will play a key role in describing this resonance.Here a coupled-channel analysis in the framework of Hamiltonian Effective Field Theory holds the promise to resolve the structure of this resonance [41].
For convenience we leave all time dependence implicit from now on and just note that the subsequent process is to be applied to each time slice.
We now introduce the jackknife estimator for the mean and uncertainty in the mean.This is a necessary step since an individual Mone Carlo sample is not an estimate of the mean.We should instead look to use averaged quantities when calculating uncertainties.
The i-th first order jackknife sub-ensemble is constructed by This is motivated by considering the average of N con − 1 samples after removing a single configuration.We can then estimate the population mean by taking an average of the N con jackknife sub-ensembles The standard deviation is then computed in terms of the jackknife sub-ensemble averages and the average of these, One can reduce this to the usual expression for the standard deviation of normally distributed data by making the assumption that each C k is a good approximate of C and thus arrive at By taking a jackknife estimate of the samples, we are using sub-ensemble means Ck which are (N − 1) 2 times more accurate than the C k when estimating the population mean C.
Up until now we have disregarded any correlation between configurations on different time slices.Since the importance sampling of the lattice action forms dependencies between the time slices, configurations on nearby time slices are in fact correlated.If we wish to perform a fit to C(t) over a time-interval, then we need a more nuanced approach to account for these cross correlations.
We hence introduce the covariance matrix which is formed by considering the deviation of the jackknife ensemble average from the average of the jackknife ensemble averages across two different time slices.Explicitly, For a fit window of T f points, this provides us with a T f ×T f matrix, with entries along the diagonal V (t i , t i ) equal to σ 2 c (t i ) using Eq.(B4).In the absence of correlations the off-diagonal factors are independent and vanish in the sum.
We can finally compute the χ 2 value for a fit T (t) to the mean, given by where the t i and t j range over all values in the interval of T f points.The inverse of the covariance matrix is computed via the singular value decomposition algorithm.
Finally, we obtain the degrees of freedom by counting the number of time slices in the fit interval T f , and subtracting both the number of parameters in the theoretical model and the number of singular values found when inverting the covariance matrix V (t i , t j ).

FIG. 1 .
FIG. 1. Plots of point sources after applying the smearing operator 16, 35, 100 and 200 times.The height of the peak is kept fixed in the visualisation to illustrate the broadening of the quark distribution.

FIG. 3 .
FIG.3.Plot of the u eigenvector components for each of the states in our 4×4 correlation matrix analysis.For each state, the quark mass increases from left to right, and the different smeared Gaussian components are denoted by the various coloured symbols.Noting the small values near zero introduce ripples but do not change the sign of the wave function within the lattice volume, one can infer the presence of nodes by counting the number of crossings through zero at a particular quark mass as the smearing number increases from 16 to 200.We find state α has α nodes when counting from zero.

FIG. 7 .
FIG. 7. J P 1/2 ∆ spectrum.CSSM results from this study are compared with the HSC [37] and Khan et al. [36].Note that the HSC and Khan et al. results include 3 states at each pion mass.

FIG. 9 . 1 / 2 −
FIG. 9. 1/2 − probability amplitudes, based on the superposition of smeared interpolating fields on the middle quark mass ensemble.These symmetric wave functions multiply the lower components of the Dirac spinors to generate odd-parity states.Nodes in the wave function are indicated by radially symmetric regions of dark blue away from the edges of the volume.

TABLE II .
Table of parameters for our lattice ensembles, characterised by the Hopping parameter κ.
. In the case of State 1, we start with u 16 ∼ u 35 ∼ 0 and then the remaining components are both far from zero, FIG.

TABLE VI .
Masses for the 1/2 + ∆ spectrum.Note that we are only able to identify plateaus in the effective mass series for the lowest-lying state at the 4 heavier quark masses.