The Positive-Parity Baryon Spectrum and the Role of Hybrid Baryons

We calculate the low-lying spectra for the positive-parity $\Delta$ and $N$ at two pion masses of 358 and 278 MeV using an isotropic clover action with two degenerate light-quaark and one strange-quark flavors through the application of the generalized variational method within the distillation framework. The spectrum exhibits the general features observed in previous calculations using an anisotropic clover lattice, with a counting of states at least as rich as the quark model. Furthermore, we identify states that are hybrid in nature, where gluonic degrees of freedom play a structural role, indicatinng that such states appear a feature of the excited baryon spectrum, irrespective of the lattice action, or the precise details of the smearing of the lattice interpolating operators used to identify such states.


I. INTRODUCTION
Lattice QCD (LQCD) provides a powerful numerical approach to solve QCD from the first principles, and has been successfully applied to address a range of key quantities in high-energy and nuclear physics, from the calculation of the ground-state spectrum, to nuclear charges and key measures of hadron structure. The calculation of the excited-state spectrum of QCD presented a particular challenge, in that the formulation of lattice QCD in Euclidean space precludes the direct calculation of scattering amplitudes.
The starting point for a study of the excited-state spectrum from lattice QCD is the extraction of the discrete, low-lying energy levels on the Euclidean lattice where lattice QCD is formulated. The most straightforward, albeit naïve, approach is to identify those energy levels with the single-or multi-particle states in the spectrum. Whilst the resulting energy spectrum should be independent of the basis of interpolating operators employed, it has the potential to be highly sensitive to that basis. Thus the first attempts to compute the spectrum employed baryon interpolating operators with a straightforward three-quark structure, mirroring the valence structure [1][2][3][4][5][6][7][8][9][10].
A focus of several studies has been the nature of the Roper resonance, the lowest-lying positive-parity excitation of the N , and whether the spectrum computed on the lattice exhibits the ordering of states revealed in nature whereby the Roper is of lower mass than the lowest-lying negative-parity excitation [11][12][13][14][15]. Here there are indications that the inverse ordering observed at the relatively large quark masses where the calculations have been performed, is reversed in the approach to the physical quark masses [13,14]. The need for a faithful representation of chiral symmetry to obtain a low-lying Roper mass has been argued in ref [15], though a calculation employing the overlap and clover fermion actions on the same set of gauge configurations suggests that some other mechanism may be at play [16]. For a more complete picture, it is necessary to include multi-hadron interpolators that would couple to states such as N π, N σ, ∆π, N ρ, and N ππ that are present in the spectrum. The calculation of the two-point functions employing such operators is computationally challenging due to the many additional Wick contractions that enter into their construction. Despite this, considerable progress has been made at including such operators in the analysis [17][18][19], though the evidence of the emergence of a low-lying Roper as a N π scattering state in these works remains scant.
The excited states are of course resonances, and a revolution over the last few years has arisen from the realization that the energy shifts at finite volume could be related to infinite-volume scattering amplitudes [20] which has transformed our ability to study the resonance spectrum, and the interaction of hadrons, from lattice QCD. The application of these methods to meson spectroscopy has been extensive 1 , and the first steps have been taken in applying them to the more computationally demanding baryon spectrum [18,19,22]. Given the computational demands that such calculations entail, a finite-volume Hamiltonian approach has been developed to relate the energies levels computable in currently practical computations to physical scattering parameters [23], and applied to the baryon spectrum in ref. [24].
The focus of our work here is on the higher states in the spectrum. Considerable insight into the excited baryon spectrum of QCD has been obtained, most notably, the extracted spectrum is found to be at least as rich as the quark model [25][26][27], and exhibits a counting of energy levels commensurate with SU(6) ⊗ O(3) spin-flavor symmetry [28]. Indeed, the ordering of the states was found to be in agreement with the predictions of the constituent quark model which can be considered as an intermediate phenomenological model consistent with the experimental results. However, such a model interprets the Roper resonance predominantly as a qqq state augmented by a meson cloud [29], and fails to find the low-mass Roper seen in experiments.
Moreover, the positive-parity excited baryon spectrum reveals suggestions of "hybrid" states, that is those in which the gluonic degrees of freedom play an essential, structural role, beyond those of the quark model, with a common mechanism with comparable states in the meson section [30]. For the case of mesons, there can exist "exotic" states that have quantum numbers J P C not available within a regular qq valence structure. Such states may have a dominant qqqq component, so-called tetraquarks [31,32], or be predominantly "hybrid" qqg states [33][34][35][36][37], with a manifestation of the gluonic valence component. Thus it is straightforward to separate an "exotic" meson state from a regular one. But for baryons, the regular qqq states can have all the J P values. So, hybrid or so-called pentaquark states will always have to"share" quantum numbers with regular states, thus making them very difficult to identify.
There has been a number of models proposed to calculate the spectrum of hybrid baryons. There is the bag model [38] where the quark and gluon fields are confined within a cavity with the fields satisfying appropriate boundary conditions at the wall of the cavity. In the flux-tube model [39,40], the quarks sit at the ends of a string-like structure. A meson contains a single flux tube between the quark and the anti-quark; in a hybrid meson, this string is excited by a transverse oscillation. For the case of a baryon, there are three tubes which either meet at a junction or form a triangle. There are also QCD sum-rule methods [41] and quark potential models [42] that make predictions about hybrid baryons.
For the case of hybrid baryons, and indeed for the non-exotic hybrids catalogued from lattice calculations in the meson sector [43], their identification has proceeded through observing that the dominant interpolating operators are"hybrid" in nature, in the sense that the operators vanish for trivial unit gauge configurations. Such an identification by its nature introduces a degree of model dependence, and therefore it is important to study the robustness of such an identification. The aim of this paper is precisely to test such robustness by performing a calculation of the low-lying positive-parity baryon spectrum at pion masses lower than those of ref. [30], using different gauge and fermion actions at a finer spatial lattice spacing.
The remainder of the paper is organized as follows. In section II, we describe the baryon interpolating operators used in our calculation, and briefly outline the distillation methodology used to construct the correlation functions and our implementation of the variational method. Section III contains details of our calculation, beginning with the parameters of the ensembles used in the calculation, a description of our fitting procedure, the robustness of the spin identification on our lattices and the stability of the fits under the variation of the parameters of distillation process. Section IV contains our results for the lowlying positive-parity ∆ and N spectrum, together with the quark-gluon assignment of the states, including those we identify as hybrid baryons. Section V contains a comparison between the results presented here, and earlier calculations using the anisotropic lattices [28,30]. In section VI, we summarize our work and outline future avenues for research.

II. COMPUTATIONAL STRATEGY
Since the focus of our calculation is the low-lying positive-parity spectrum, we employ a basis of interpolating operators that have been found to have the dominant overlaps with those states. The construction of the interpolating operators, and the identification of the operators that couple primarily to the low-lying spectrum, has been described in detail in refs. [28] and [30], so we only summarize the salient elements here. The interpolating operators follow a continuum construction, and are expressed as a product of terms describing the flavor structure, Dirac spin and orbital angular momentum implemented through derivatives: where B, S and D denote the flavor, Dirac spin and orbital angular momentum, L, components respectively, and Σ F , Σ S , Σ D are the corresponding permutation symmetries. The resulting operators are projected through suitable Clebsch-Gordon coefficients to total spin J; the label n distinguishes different combinations that have the same spin structure, while the d is the order of the gaugecovariant derivative. For this work, our basis comprises the non-relativistic operators constructed from the upper components, in a suitable γ representation of the Dirac spinors, with up to two covariant derivatives, allowing the operators with up to two units of orbital angular momenta. We also include additional operators containing the commutator of two covariant derivatives acting on the same quark field, corresponding to the chromomagnetic components of the gluonic field-strength tensor and denoted by D [2] L=1,M in the following; it is these operators, which vanish for a unit gauge configuration, that are referred to as "hybrid" operators, and for which a dominant overlap with a given state, we treat as the signature of the hybrid nature of that state, as we discuss below. Finally, as the calculations are done on a discretized lattice, the operators are subduced from the continuous Hilbert space onto the different lattice irreps. H g , G 1g and G 2g , where the subscript g denotes positive parity, such that the correlators constructed using these operators, receive contributions in the forward direction only from the positive-parity states. As a consequence, for total angular momentum J = 5 2 and higher, the continuum operators are subduced onto multiple irreps, as detailed in Table I.
The operators created directly from the fields of the lattice Lagrangian couple to states at all scales, thus making the extraction of the lightest states in the spectrum difficult. In order to solve this problem, a linear operator is applied on the quark fields on appropriate time-slices and operators are built from those "smeared" fields. In this work, the smearing method used is known as Distillation TABLE I: The numbers of ∆ and N interpolating operators used in the calculation, together with their subductions onto the irreps. of the cubic group. In the final two columns, the number in brackets denotes the number of hybrid operators within each irrep.. [44]. The distillation operator is defined as: x (t) is the k th eigenvector of the second-order three-dimensional differential operator, 2 evaluated on the background of the spatial gauge fields of time-slice t, once the eigenvectors have been sorted by the ascending order of the eigenvalues; N D is the dimension of the distillation space.
The reasons for adopting distillation in our calculation are two-fold. First, the computationally demanding parallel transporters of the theory, the perambulators, depend only on the gauge field, and not on the interpolating operators. So, we can calculate the perambulators on an ensemble of gauge field once, and then reuse them for an arbitrary basis of operators at both the source and the sink, and indeed for a range of calculations. Second, the method provides a more complete sampling of each gauge configuration through imposing a momentum projection at both the source and the sink correlation time-slices. The expectation is that N D should scale as the physical volume, and the cost of computing the corresponding correlation functions scales as N 4 D for the case of baryons. Thus there is a computational imperative to use as small a distillation space as possible while still providing a faithful description of the physics, and we will investigate the sensitivity of our results to N D below. An approach that aims to overcome this scaling with the spatial volume is to stochastically sample the eigenvector space, a method known as LaPH [45], which we do not pursue here.
A variety of approaches have been developed to obtain the energies and operator overlaps for the excited baryon spectrum [6,46]. Here we use the variational method as implemented in ref. [47]. Our starting point is the generalized eigenvalue equation (GEV) for the two-point where without loss of generality we take the source interpolating operator to be at time slice t = 0, and where i, j label the operators in a given representation of the cubic group. The GEV equation is expressed as where u α are the generalized eigenvectors which satisfy the orthonormality condition u † α C(t 0 ) u β = δ αβ , and the corresponding principle correlators are λ α (t, t 0 ) behaving as Here m α is the energy of the state labeled by α and δm represents the contributions from other states. Our subsequent results are derived from two-state fits to the principle correlators of this form. Furthermore, C ij (t) can be decomposed into the form, where the overlap factor Z α i = 0|O i |α can be written as, The matrix U is formed using the generalized eigenvectors u α as its columns. The overlaps can thereby be obtained from the solution of the generalized eigenvector matrix.

III. COMPUTATIONAL DETAILS
Earlier calculations using this basis of operators were performed on anisotropic clover lattices, with a spatial lattice spacing of around a 0.123 fm, and an anisotropy, ξ ≡ a s /a t 3.5 [48,49] with two mass-degenerate light-quark flavors and a strange-quark. Here we use an isotropic clover action at a smaller lattice spacing a 0.094 fm, determined using the w 0 scale [50], and likewise with 2 + 1 flavors. We use ensembles at two values of the light-quark masses, corresponding to pion masses of m π = 358 and 278 MeV respectively. Details of the parameters of our ensemble are listed in Table II. All the gauge-links entering in the operator constructions are stout-smeared [51]. In order to achieve the best possible sampling of the lattice, we evaluate the two-point correlators from each time-slice on the lattice, and on each configuration, average the correlators over the different time sources to account for the correlations along the temporal direction of the lattices. We compute the perambulators for N D = 64 eigenvectors.
FIG. 1: Fits to the principal correlators for the low-lying positive-parity spectrum of the H g irrep. of the ∆ on the ensemble a094m358, using t 0 = 5. The plots show λ α (t, t 0 ).e mα(t−t0) data on the y-axes and the lattice time-slices on the x-axes; the curves are two-exponential fits as described in the text. In each panel, the mass corresponding to the leading exponential state is labelled by m and given in lattice units.   [50], and N cfg is the number of configurations.

A. Fitting Procedure
The fitting procedure we employ to extract the mass spectra and the overlap factors is discussed in detail in ref. [34], and we only summarize the procedure here. The GEV of eqn. 2 is solved over a range of t 0 , and for each t 0 , the resulting principle correlators are fit to the twoexponential form, We restrict the fitting range such that, for each principle correlator, we only include time-slices for which the noiseto-signal ratio is less than 0.05; in practice, this restricts the largest value of t included in the fits to be around 8. Furthermore, we only include the correlators with source-sink separation greater than two lattice units to avoid possible contact terms. For each t 0 , our fit to each principle correlator is based on an acceptable χ 2 /dof. We also require that the coefficient A in eqn. 7 is less than around 0.1, and that, for a fit to an N -dimensional matrix of correlators, the matrix is largely saturated by the lowest-lying N states reflected in that, for each α, m α is larger than the lowest-lying masses, m α obtained for all the principle correlators. Finally, we compute the overlap factors of eqn. 6 using the eigenvectors at a reference time-slice, t z . In Figure 1, we show fits to the leading the principal correlators for the H g irrep. of the ∆. For each panel, the curve is the reconstruction from the fitted parameters with the purple region indicating the data points included in the fit. The approach of the plateaux close to unity at large times is indicative of the small value of A in the fits, and the small contribution of the second excited state to each principal correlator. As anticipated, the uncertainties on the leading mass increase as the mass increases.

B. Spin Identification
The breaking of rotational symmetry induced by the discretization onto the lattice renders the determination of the spin corresponding to the different energy levels within the irreps. less than straightforward. In the case of the glueball spectrum in pure Yang-Mills theory [52], the identification of the spins was accomplished by the identification of degeneracies across different lattice irreps in the approach to the continuum limit. This requires the generation of ensembles at several lattice spacings, a formidable task once quark degrees of freedom are included, and further requires statistical precision far beyond that attainable with reasonable computational cost to delineate overlapping energies within the spectrum. We need a spin identification method which uses data obtained from only a single lattice spacing, albeit one sufficiently fine that it preserves the rotational symmetry to a sufficient degree at the hadronic scale.
FIG. 2: Histogram plot of the operator overlaps Z for the ∆ on the a094m358 ensemble, normalized such that, for a given operator, the largest overlap across all the states is unity. The overlaps are obtained from a variational analysis across all operators within a given lattice irrep, irrespective of the continuum spins from which they are derived. The histogram plot of each state is accompanied by its mass in lattice units. The asterisks denote the hybrid-type operators, and the energy levels identified with them.
Here we use the method introduced in ref. [35], and applied for the baryon spectrum in refs. [28,30] whereby the operator overlap factors are used to identify the spin of a state. It relies on the observation that each operator used in the calculation carries an essence of the continuum spin of the operator from which it is subduced, and therefore we would expect an operator subduced from, say angular momentum J, to have large overlaps only with states of the same continuum angular momentum J. Positive-parity states corresponding to the continuum angular momentum J = 5 2 and J = 7 2 will appear in the spectrum of the H g and G 2g irreps, and of the H g , G 1g and G 2g irreps respectively, and we would expect overlaps to be dominated by the operators subduced from the same continuum operators across those irreps. This is indeed what we observe, as can be seen in Figure 2 for the ∆ spectrum, where the overlaps are obtained from a variational analysis using all the operators within a given lattice irrep. Further, we find the resulting energies are degenerate, with, for states of spin 5 2 and 7 2 , the energies obtained in the H g irrep. within 1% of the values obtained in the G 1g and, for the case of spin 7 2 , G 2g irreps. Rather than applying the variational method to a basis comprising all the operators within a lattice irrep., we instead apply the method to a more restricted basis of operators comprising those operators within an irrep. derived from a given continuum J. In Figure 3, the ∆ spectrum obtained by analysing all the operators within a given lattice irrep. is compared with that where we apply the varational method in each lattice irrep. to only those operators derived from a given continuum spin. The comparison reveals that there are no significant differences between these two spectra, prompting us to analyse the operators of each angular momentum separately as this requires calculating the two-point correlators using a smaller basis of operators at one time; since the computational cost of computing the full correlation matrix goes as the square of the operator basis, this reduces the computational cost significantly.

C. Stability under variation of Distillation Space
The previous studies of the low-lying baryon spectrum using this implementation of distillation employed N D = 56 distillation eigenvectors, on a 16 3 spatial lattice with a 0.123 fm. With the expected scaling of the number of eigenvectors with physical volume discussed earlier, that would suggest that as many as 230 eigenvectors might be needed to capture the same physics on the ensembles employed here, with in excess of three times the physical spatial volume. In this paper, we have generated perambulators and the baryon elementals that encode the operators for N D = 64 eigenvectors, and begin our discussion by examining the sensitivity of the extracted spectra to the variation of N D . A study of the various charges of the N , both for a state at rest [53], and at non-zero spatial momentum [54], performed on the same lattices as those employed here, suggests that the ground-state properties can indeed be resolved with this number of eigenvectors. However, ascertaining the sensitivity of our results for the spectrum and overlaps of the excited states is an important prerequisite for our subsequent discussion.
In Figure 4, we show the lowest energy levels in the positive-parity H g irreducible representation of the ∆ as we reduce the number of eigenvectors down to N D = 24. While the ground state is indeed reliably extracted with only a minimal number of eigenvectors, it is only when we reach N D = 56 eigenvectors that the lowest five states are obtained with acceptable uncertainties, with consistency between the N D = 56 and N D = 64 determinations. Since a major aim of this work is establishing evidence for the properties of the extracted states, as evidence through the operator overlaps, we likewise need to ensure that these features are robust under the rank of the distillation space. In Figure 5, we show the operator overlaps corresponding to the energies in Figure 4; as for the case of the energies, the overlaps, and in particular the dominant operators corresponding to each state, show stability between N D = 56 and N D = 64, but with some qualitative differences, notably in the ordering of the states, for N D = 48. We will therefore use N D = 64 in the remainder of this paper.

IV. RESULTS
The low-lying positive-parity spectra of the ∆ and N for both the a094m278 and a094m358 ensembles using this fitting procedure are shown in Figures 6 and 7, respectively. For the spin 5 2 and spin 7 2 energy levels, the splittings between the values obtained in the H g and G 2g , and in the H g , G 2g and G 1g irreps, respectively, are remarkably small, reflecting the partial O(a 2 ) breaking of rotational symmetry, and the smaller spatial lattice spacing than that used in comparable studies using an anisotropic lattice. As expected, the quality of the spectrum is somewhat worse at the lighter value of the quark mass, and the spin identification procedure less convincing for the highest states in the N spectrum. We emphasise that the qualitative properties of the spectrum, and in particular the counting of states, is consistent with that obtained on the anisotropic lattices at a coarser value of the spatial lattice spacing, but a considerably finer temporal lattice spacing. Notably, our calculation does not exhibit the low-lying Roper resonance, in accord with calculations using the anisotropic action, and indeed most calculations using a Wilson-type action.

Hybrid States
As we noted in the introduction, in contrast to the case of the meson spectrum, "exotic" baryons cannot be distinguished through their quantum numbers. Therefore, the identification of baryons as "hybrid" in nature FIG. 6: The low-lying positive-parity ∆ spectrum in lattice units on the a094m278 (left) and a094m358 (right) ensembles, using the fitting procedure described in the text. For the states identified as spin 5 2 and 7 2 , the boxes contain the energy levels obtained after the subduction onto the different lattice irreps. Energy levels identified as those of hybrid baryons are denoted by the green asterisks.
inevitably involves a degree of interpretation. Here we identify the hybrid states as those whose overlap, defined through eqn. 6, is dominated by the hybrid-type operators, that is those that would vanish for the case of a trivial gauge configuration [30]. For the case of the ∆, this identification is very apparent, as can be seen in Figure 2 for the a094m358 ensemble, where we find one hybrid state in the J = 1 2 channel and one in J = 3 2 channel. For the N spectrum on the a094m358 ensemble, we likewise find clear evidence for hybrid-baryon states through the nature of their overlaps, where we identify two states in the J = 1 2 channel, two states in the J = 3 2 channel and one state in the J = 5 2 channel. On the a094m278 ensemble, the identification and multiplicities of the hybrid baryons follow those of the a094m358 ensemble except for the J = 1 2 channel in the N spectrum, where there is no obvious candidate for a hybrid baryon using the criterion of the operator overlaps. In spite of this, the multiplicity in both the ∆ and N spectrum confirm the findings in the earlier studies using the anisotropic lattice [30], with a multiplicity of states at least as rich as the quark model, and the presence of additional states that appear to be hybrid in nature.

V. DISCUSSION
We now compare our results with those of previous works, and in particular the previous calculation of the low-lying positive-parity baryon spectrum obtained on the heavier of the two anisotropic clover lattices employed in ref. [30]. To facilitate this comparison, we consider the excitation energy with respect to the ground-state Nucleon mass, in units of the Ω mass, a quantity that is somewhat insensitive to the light-quark masses. In Figures 8 and 9, we show the comparison among these lattices for the J = 1 2 and J = 3 2 channels for the ∆ and N , respectively. Also shown are the lowest-lying noninteracting two-particle energy levels.
A notable feature of most states for both the ∆ and N is that the splitting with respect to the ground-state Nucleon mass shows only a weak dependence on the quark masses, while the energies of the non-interacting twoparticle states exhibit a far stronger dependence. This suggests that we are observing predominantly "singlehadron" states rather than multi-hadron states, and leads further support to our assertion that the three-quark operators used in this study couple only weakly to the multi-hadron states. These observations are more prominent for the hybrid baryons whose masses, with respect to the ground-state Nucleon mass, remain more or less the same irrespective of the quark masses. However, there is one qualification to this observation, namely that the first excited-state energy seen in the N 1 2 channel exhibits a stronger dependence on the light-quark masses, and is indeed consistent with that of N (0) π(0) π(0) and N (1) π(−1) multi-hadron states.
A focus of this paper is whether the identification of hy-FIG. 7: The low-lying positive-parity N spectrum in lattice units on the a094m278 (left) and a094m358 (right) ensembles, using the fitting procedure described in the text. For the states identified as spin 5 2 and 7 2 , the boxes contain the energy levels obtained after the subduction onto the different lattice irreps. Energy levels identified as those of hybrid baryons are denoted by the green asterisks. brid baryons is indeed robust. In Figures 10 and 11, we show that the dominant operators for each of the states of all three ensembles for J = 3 2 in the case of the ∆, and for J = 1 2 in the case of the N . For the ∆, the hybrid baryon in each of the ensembles has almost identical overlap distribution across the operators, with the hybrid operator having the predominant overlap. The ground states also have a comparable distribution across the three ensembles, though we note that the work here includes an additional operator whose orbital structure is of the form D [2] L=0,S that can be interpreted as an operator of additional width with respect to the S-type orbital D L=0,S , and therefore of the same orbital structure. For the case of the N , the identification and ordering of the hybrid baryons is comparable for the heavier a094m358 ensemble and the anisotropic ensemble, but for out lighter a094m278 ensemble, that identification is less obvious in spite of having significant overlap from the hybrid operators. As in the case of the ∆, there is consistency in the overlaps for the ground state and the first two excited states across all three ensembles.

VI. CONCLUSIONS
In this work, we have computed the positive-parity ∆ and N spectra using an isotropic clover action. Our results support the observations in earlier works at heavier pion masses, and using the anisotropic clover action at a coarser spatial lattice spacing, but finer temporal lattice spacing. In particular, we find that rotational symmetry is largely observed at the hadronic scale, enabling us to reliably identify the spins of the states through their predominant overlap of operators derived from continuum operators of definite spin. However, the most significant outcome of this work is that we find that the spectra exhibit a counting of states in line with that of the quark model, but with additional states that we can identify as "hybrid" in nature, with the gluonic degrees of freedom playing a structural role. The means used to identify such hybrids through the predominant overlap of a class of "hybrid" operators, pioneered in ref. [30], must inevitably raise the issue of the operator dependence of such an identification. Here we use a different action, with a different lattice spacing and essentially different interpolating operators implemented through the variation of the number of distillation eigenvectors. Thus the identification of hybrid-type states in the spectrum is indeed robust.
This work has important limitations in its use of "single-hadron" operators which do not fully capture the low-lying energy levels in the finite-volume spectrum. The next step in the investigation of the nature of "hybrid" baryons would be to include the multi-hadron operators, and subsequently to compute the infinite-volume momentum-dependent phase shifts. Such a study could also reveal the decay modes of such states, and indeed the first study of the decays modes of the exotic 1 −+ hybrid FIG. 8: The left and right panels show the excitation energies for the ∆ with respect to the ground-state Nucleon on each of our lattices vs. m 2 π , together with the corresponding result from ref. [30], for the J = 1 2 and J = 3 2 channels respectively. We use the Ω mass to set the scale. The higher excited states are displaced for clarity. The N (0) π(0) π(0) and N (1) π(−1) energy thresholds are identified by horizontal dashes.
has recently been performed [55]. This is more computationally challenging for baryons than for mesons through the increased cost of Wick contractions, the scaling of the number of distillation eigenvectors with increasing volume, and the numerous final states to which they can decay. Nonetheless, the advent of the exascale era of computation makes such computations increasingly realizable. Ultimately, a largely model-independent determination of the quark and gluon content of such resonances will be achieved by the probing of their structure through external currents, and the theoretical framework for such studies is an area of rapid development [56][57][58] and application [59,60].

VII. ACKNOWLEDGMENTS
We thank Jozef Dudek, Robert Edwards, Archana Radhakrishnan and Christopher Johnson for useful discussions, and for the use of the reconfit fitting package. This work is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177. Computations for this work were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy. This work was performed in part using computing facilities at The College of William and Mary which were provided by contributions from the National Science Foundation (MRI grant PHY-1626177), and the Commonwealth of Virginia Equipment Trust Fund. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. Specifically, it used the Bridges system, which is supported by NSF award number ACI-1445606, at the Pittsburgh Supercomputing Center (PSC) [61,62]. In addition, this work used resources at NERSC, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract #DE-AC02-05CH11231, as well as resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. #DE-AC05-00OR22725. The software codes Chroma [63], QUDA [64,65] and QPhiX [66] were used in our work. The authors acknowledge support from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research and Office of Nuclear Physics, Scientific Discovery through Advanced Computing (SciDAC) program, and of the U.S. Department of Energy Exascale Computing Project. TK was support in part by the Cen-FIG. 9: The left and right panels show the excitation energies for the N with respect to the ground-state Nucleon on each of our lattices vs. m 2 π , together with the corresponding result from ref. [30], for the J = 1 2 and J = 3 2 channels respectively. We use the Ω mass to set the scale. The higher excited states are displaced for clarity. The N (0) π(0) π(0) and N (1) π(−1) energy thresholds are identified by horizontal dashes.