Baryon Spectroscopy in Lattice QCD

A variational analysis is performed within the framework of lattice QCD to extract the masses of the spin-3/2 positive parity (cid:2) + baryons, including radial excitations. 2+1 ﬂavour dynamical gauge-ﬁeld conﬁgurations provided by the PACS-CS collaboration via the ILDG are considered. To improve our interpolator basis, we perform an iterative process of source and sink smearing and solve a generalised eigenvalue problem using the resulting fermion operators. We obtain a clear signal for the ground and ﬁrst excited states at a light quark mass corresponding to m π = 413MeV. Furthermore, we show that one can use the eigenvectors obtained in this method to investigate the nature of these states, allowing us to classify our results as 1 s and 2 s states for the ground and ﬁrst excited states respectively. Finally, we brieﬂy highlight the method of Hamiltonian effective ﬁeld theory which can be used to make comparison with quark model expectations.


Introduction
The strong nuclear force underpins much of the world around us, giving rise to the hadronic matter in our universe through the interaction of quarks and gluons.While the equations governing this fundamental force of nature are captured by the theory of Quantum Chromodynamics (QCD), a complete understanding is elusive due to the non-perturbative nature of QCD in the low energy regime.This is precisely the realm of hadrons such as the nucleon, which is so fundamental to our understanding of the universe.
Fortunately, lattice QCD is a technique designed specifically for tackling the nonperturbative part of QCD and has seen great advances in its application alongside developments in supercomputing technology.By discretising space-time and employing a number of systematically improvable approximations, the once impossible task of computing the mass of low-lying baryons is now tractable at percent-level precisions.
The nucleon spectrum has historically presented an interesting problem for the community, focused on the N * (1440) Roper resonance [1,2].The Roper has been of particular interest in the past due to its stark juxtaposition with simple theory expectations; while it is the second state in the nucleon spectrum and has positive parity, a simple quark model without fine tuning predicts the second state should have negative parity.The closely related ∆ baryons also exhibit this anomalous energy level ordering and present a similar issue for our understanding of the strong force.
One way of understanding the Roper and other low-lying baryon resonances such as the ∆ has developed in recent years through Hamiltonian Effective Field Theory (HEFT) [2][3][4][5][6][7][8].Traditionally, this approach fits experimental data and brings this information to the finite volume of the lattice in the form of predictions of the finite volume spectra.Thus one main aim of this paper is to provide results for the ∆ baryon spectrum which can be confronted with future HEFT analyses.
Another aim for this work is to present the lattice techniques of spin and parity projection coupled with a variational analysis to isolate states in the spectrum.These techniques have been developed separately in different contexts [9,10] and we aim to show how together they allow one to isolate and classify states in, for example, the J P = 3/2 + , ∆ + spectrum.

Energies from Correlation Functions
In lattice QCD we calculate correlation functions by first calculating the relevant quark propagators.At the baryon level, these correlation functions are defined as The operator χj ν (0) acts on the vacuum to create states at the space-time point 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).T indicates time ordering of the operators.
If we are interested in the excited states of hadrons, we choose a set of operators with some non-zero overlap with the states of interest.For baryons, we have the complete set of states where B labels different baryons with momenta p and spins s.This allows us to obtain We then use the translational operator to rewrite χ i µ (x) as This allows us to simplify the expression for the correlation function to Now we set the momentum to p = 0, so that E B = M B .We also introduce the parity projection operators so one can then compute the parity projected correlator and sum over the spatial Lorentz By replacing the matrix elements in Eq. ( 5) with appropriate functions (i.e.Rarita-Schwinger spinors for the case of spin-3/2 baryons, or regular Dirac spinors for spin-1/2), one can show [11] that 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 ± .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 [9,11].
In order to extract 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 to the lowest lying state.With this preliminary discussion out of the way, we now discuss a number of techniques used to extract the various states in the ∆ spectrum.

Spin Projection
We use a standard interpolating field for the ∆ + given by This operator has overlap with both spin-1/2 and spin-3/2 states so we need to perform spin projection to guarantee we extract the masses of the desired states.The spin-1/2 and spin-3/2 projection operators are given by [10,12] Although Eq. ( 11) looks somewhat involved and cumbersome, we can do a few things to simplify these operators.First, in our lattice calculations the mass is on shell, so we have We'll consider the particles being at rest, p = 0, in which case we get p 0 = m and as usual p 2 = m 2 .Further, we also have With these results in mind, we simplify our projection operator to One can then show by using the properties of the γ-matrices and the metric, that the elements of the projectors obey (17) where m, n are spatial Lorentz indices.
We can immediately get the corresponding results for the spin-1/2 projection operator by making use of Eq. ( 12) and ( 17) and we obtain (18) With the spin projection operators in hand, a spin-s projected correlation function is then given by G µσ g σλ P s λν .
This spin projection is performed prior to the parity projection and trace in Eq. ( 8).Thus we can obtain results for the masses of both the spin-1/2 and spin-3/2 ∆ + states without needing to generate the correlation functions more than once, saving on overall compute time.

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.Smearing allows us to more accurately represent a fermion bound within a hadron by smearing out a point source/sink so it achieves some finite width.This is demonstrated in Figure (1) where higher levels of smearing lead to Gaussian distributions of larger RMS radius.This broadening of the sources and sinks leads to an improved overlap of the resulting interpolating field with the states in the spectrum.Additionally, a combination of smeared sources and sinks can also give an indication of radial excitations in the resulting spectrum.This will be discussed further in Section 8.

Variational Analysis
We now resume our main discussion of how one extracts the masses of a baryon spectrum.Recall that in Section 2 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. ( 9).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 documented variational analysis method [9,13] in order to extract the ground state mass, as well as the 1st excited state mass.We also obtain results for the 2nd and 3rd excited state masses, and will report those in a forthcoming publication.
To extract N states, we use N interpolating fields (one for each state).We generate our interpolating fields by using the same base field for the ∆ + baryon, and applying different levels of smearing.These smearing levels are chosen for their coupling to the states of interest, and the use of smearing on both sources and sinks allows us to construct a matrix of correlation functions.This correlation matrix is written as Here the λ α i and λ α j are essentially the same as seen before in Section 2, though we now use the α index to distinguish energy states.In other words, they are couplings of the interpolators χ i µ and χ j ν at the source and sink to the various energy eigenstates α = 0, . . ., N − 1.Finally, m α is the mass of the state α.
From here, we now aim to construct linear combinations of our operators 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 superposition of fields, using the basis of smeared interpolating fields.
At this point, we will attempt to construct an eigenvalue problem to solve for both u α and v α .Noting that since G ij (t) = G ji (t) in the ensemble average, we introduce an improved unbiased estimator of the correlation matrix 1/2[G ij (t) + G ji (t)].This provides us with a correlation matrix which is symmetric, so we can simultaneously compute u α and v α as discussed below.
Multiplying Eq. ( 22) 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, by premultiplying Eq. ( 22) by v α i we get from which we can arrive at our second eigenvalue equation (this time for v α ): Both Eq. ( 29) and (31) need to be solved simultaneously for each given pair of variational parameters, 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 We can clearly see that these eigenvectors have isolated a single state in the baryon spectrum, exactly as we set out to do.We then construct the effective mass function As a note, δt should not be confused with ∆t.The latter is a variational parameter which we set to allow us to get a few time slices away from the source time.The former 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 α must also be 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. ( 33) can be computed for various discrete values of t and then plotted as a function of time.As usual, one looks for time intervals over which the effective mass plot plateaus, indicating that all contamination from unwanted states has decayed away in the exponential series.We finally perform a covariance matrix analysis [14] to determine the most suitable time intervals to fit when obtaining our final masses.
We note that the principles highlighted within this section are realised only when one includes a complete set of interpolating fields effective at isolating all the states within the spectrum.This needs to include multi-particle scattering states.Our formalism is focused on the single-particle states and in particular their radial excitations.As such, Euclidean time evolution is important in suppressing contaminations from nearby scattering states.In fitting our effective mass to a plateau, we ensure single-state dominance by monitoring the χ 2 per degree of freedom (χ 2 /dof).

Lattice Details
We have simulated ∆ + baryons on a 32 3 × 64 lattice with gauge-field configurations provided by the PACS-CS Collaboration [15] through the ILDG [16].These configurations use an order a improved Wilson fermion action and an Iwasaki gauge action.The lattice is subject to periodic boundary conditions.The correlation functions were constructed using the COLA software library [17].
Our calculations are performed on gauge-field ensembles with β = 1.90 and a quark mass given by the hopping parameter κ = 0.13754, providing m π = 413 MeV.At this quark mass, and with a Sommer parameter of r 0 = 0.4921(64) fm, the lattice spacing is a = 0.0961(13) fm.
To improve the signal-to-noise ratio of our results, we perform shifts on the source insertion point and then average the resulting correlation functions.We found that at our chosen quark mass, an average over 16 source locations was sufficient.
The smearing levels used in this analysis are 16, 35, 100 and 200 sweeps.We also investigated several sets of variational parameters before settling on t 0 = 18 and ∆t = 2 relative to a source at t s = 16.These choices have been motivated in similar studies in the past [9].

Numerical Results
We present results for the ∆ 3/2 + spectrum in Figure (2).We have thus far been able to extract masses for the first 2 states in the spectrum and future work will aim to also extract the 2nd excited state.
Consider Figure (2), where our results are labelled 'CSSM' and shown as red squares alongside other contemporary results [15,18,19] having a similar quark mass.While there is some indication of the ground state in our low-lying mass, there is a significant Fig. 2 Our obtained masses for the ∆ after spin-3/2 and positive parity projection (CSSM).Resonances measured experimentally are the ∆(1232), ∆(1600) and ∆(1920), indicated by crosses along the grey dashed line at the physical point.The other studies listed here at similar quark masses to that presented here are from the PACS-CS Collaboration [15], Khan et al. 2021 [18] and the Hadron Spectrum Collaboration [19].Note that the PACS-CS result is offset from our own for clarity.
jump to get to the second state in our spectrum.This jump is more than one might expect for a state which tends to the ∆(1600) in the physical limit.In particular, we emphasise that we don't appear to see this state using a simple 3-quark interpolating field, as given by Eq. (10).
To finish this section, we have shown that a variational analysis allowed us to successfully extract ∆ masses for the ground and first excited states.We report these masses as m 0 ∆ = 1.439 ± 0.007 GeV, χ 2 /dof = 0.972 (34) where the χ 2 /dof is computed via a covariance matrix analysis [14].

Eigenvector Analysis
One can also leverage the eigenvectors u α obtained in our variational analysis.This has been shown useful in identifying radial excitations in the nucleon spectrum [9,20] and we aim to apply similar techniques here.Recall that the eigenvectors were used to construct a superposition of smeared interpolating fields in Eq. ( 24) and (26).We could more explicitly think of these superpositions as being written as where each of the χ Nsm is a fermion field which has been smeared N sm times.As usual, the superscript α is just the energy level.Now, we already identified that smearing increases the extent of our source/sink operators.Solving the variational analysis problem is equivalent then to finding a superposition of Gaussians which has largest overlap with the wave function of the desired ∆ state.
This simple reinterpretation of the superpositions gives us a powerful tool for classifying the states in the spectrum.Consider the case of superposing, for example, a narrow Gaussian with positive sign, with a wider Gaussian of negative signature.The resulting distribution will have a crossing through zero at the edge of the positive Gaussian, where the superposition becomes dominated by the negatively weighted Gaussian.In other words, if there is a relative sign between consecutive Gaussians in our superposition, this will indicate a zero crossing or a node in the radial part of the wave function.
Consider Figure (3).This illustrates the superposition of Gaussians via the u α i components obtained through the variational method.More specifically, let us focus on the subplot for State 0. All components of the eignvector are positive or approximately zero, so we would be inclined to think there is little evidence for a zero crossing.We would thus identify this state as a simple 1s state, since it has no nodes.Moving onto the State 1 plot, we see that both the 16 and 35 smearing components are centred about zero.The 100 component is clearly positive and large, and the 200 component is negative and large.We would interpret this exactly as the toy case above, where there is a narrow peak of positive signature superposed with a broader peak of negative signature.Once the positive peak becomes less prominent, the negative Gaussian will take over and there will be 1 zero crossing.We would identify this as a 2s state as the superposition has a single node.Similar analysis of the State 2 and 3 subplots shows 2 and 3 nodes respectively, so we would identify these as 3s and 4s states respectively.
A more involved analysis can be performed whereby one computes the wave function of one of the light quarks within the baryon directly using lattice QCD [20].Such a computation would be a useful point of comparison, but leveraging the variational method eigenvectors as discussed above offers an alternative for classifying states in the spectrum.

Hamiltonian Effective Field Theory
Looking ahead to future studies, we briefly allude to Hamiltonian Effective Field Theory (HEFT) and its use in studying low-lying baryon resonances [2,8].Essentially, HEFT bridges the gap between the finite volume of the lattice, and what is observed in the infinite volume of the real world (at particle collider experiments for example).By considering a Hamiltonian composed of single-particle bare states and meson-baryon multi-particle states, one is able to construct a scattering model which can be fit to experimental observables, such as phase shifts and inelasticities in, say, πN scattering.
Such a model can then be brought to the finite volume of lattice QCD and extrapolated to unphysical quark masses to confront lattice QCD data, such as that obtained in Figure (2).Having a broad quark mass range to perform fits to affords a good opportunity for such comparisons.HEFT has been shown to provide new understanding for the nature of the Roper resonance [2,6], the Λ(1405) [4,7] and for the low-lying ∆ spectrum [3,8].We aim to extend this work using the updated lattice QCD masses reported in this paper.

Conclusion
In this preliminary investigation into the ∆ 3/2 + spectrum, we have obtained ground state and first excited state finite-volume masses at a quark mass corresponding to m π = 413 MeV.Our ground state results are comparable with the earlier PACS-CS results, and our 1st excitation is similar to other recent analyses, sitting higher than the expected ∆(1600) state.
Building on previous studies using smeared interpolating fields, we have been able to identify both states in our spectrum by their nodal structure as being a 1s state and 2s state.We hope to take this work and make connection to experimental data through HEFT and thus glean some deeper understanding of resonance structure for the low energy ∆ excitations.

Fig. 1
Fig. 1 Gaussian distributions resulting from smearing a point source with 16, 35, 100 or 200 sweeps of smearing.Note the progressively wider Gaussians obtained from repeating the smearing.

Fig. 3
Fig. 3 Components of the normalised eigenvectors u α i .The figure is composed of four subplots with each labelled along the x-axis by state α.The components of the eigenvectors are labelled by the number of smearing sweeps, indicated in the legend.