Charmed baryon spectrum from lattice QCD near the physical point

We calculate the low-lying spectrum of charmed baryons in lattice QCD on the $32^3\times64$, $N_f=2+1$ PACS-CS gauge configurations at the almost physical pion mass of $\sim 156$ MeV/c$^2$. By employing a set of interpolating operators with different Dirac structures and quark-field smearings for the variational analysis, we extract the ground and first few excited states of the spin-$1/2$ and spin-$3/2$, singly-, doubly-, and triply-charmed baryons. Additionally, we study the $\Xi_c$-$\Xi_c^\prime$ mixing and the operator dependence of the excited states in a variational approach. We identify several states that lie close to the experimentally observed excited states of the $\Sigma_c$, $\Xi_c$ and $\Omega_c$ baryons, including some of the $\Xi_c$ states recently reported by LHCb. Our results for the doubly- and triply-charmed baryons are suggestive for future experiments.


I. INTRODUCTION
Recent experimental results from the LHCb Collaboration on the Ω c , Ξ c and the doubly charmed Ξ cc state have put further emphasis on the relevance of the hadron spectroscopy. There now exist 31 observed charmed baryons, 25 of which are classified with at least three stars by the Particle Data Group (PDG) [1]. Charmed baryons provide a unique laboratory to study the strong interaction and confinement dynamics due to the composition of the light and charm quarks. Studying the excited states of the charmed baryons has the potential to reveal their internal dynamics and the nature of the excitation mechanisms.
Experimentally, the singly-charmed baryon sector is most accessible. Within this sector, the Λ c channel is most established. In addition to the ground state, there are four excitations with total spin up to 5/2, although in need of a confirmation of the assigned quantum numbers. Out of the three Σ c states that are listed by the PDG, two are the lowest J P = 1/2 + and 3/2 + states and Σ c (2800) is their only observed excitation. This state has been detected in the Λ c π channel by the Belle [2] and the BABAR [3] Collaborations. Its quantum numbers are not measured. In contrast, the Ξ c sector is quite rich since it can have flavor symmetric and antisymmetric wave functions. There are up to seven Ξ c excitations observed by the Belle [4][5][6][7][8][9], the BABAR [10,11] and very recently by the LHCb [12] Collaborations in the energy range of 2920 to 3120 MeV/c 2 . The PDG considers the existence of three of them to be very likely or certain while the confidence for the other two is smaller. LHCb states are not included in the review yet. These excited states appear in the invariant mass distributions of several singly-charmed baryon B c + K or π channels depending on the strangeness number of the baryon and in the ΛD channel where the charm quark is confined in the meson system. This unique behavior makes the Ξ c system a good laboratory to study the internal excitation dynamics of the charmed baryons and the diquark correlations. The quantum numbers of these states remain undetermined. The LHCb Collaboration has also reported the precise measurements of the masses and the decay widths of five new Ω 0 c states [13], which are observed in the Ξ c K channel in the energy range from 3000 to 3120 MeV/c 2 . Their spin-parity quantum numbers remain undetermined. There are several works in the literature investigating the nature of these states and assigning conflicting spin-parity quantum numbers. It is a triumph of the experiments to identify many states in such narrow energy windows.
The lowest-lying states of the singly-charmed baryons are already established by experimental studies and the lattice QCD results agree well with those observations. The Ξ cc is the only observed doubly-charmed baryon for the time being. It was first observed by the SELEX Collaboration [14,15] but its results were not confirmed by other experiments until the LHCb Collaboration has reported the same particle with a different mass [16]. Lattice QCD predictions for the mass of the Ξ cc lie above the SELEX reported value but agree very well with the LHCb value.
The lowest-lying charmed baryon states have been studied by various lattice groups as well. Early investigations utilized the quenched approximation [58][59][60][61][62], while recent studies employ up to 2 + 1 + 1-flavor dynamical gauge configurations with several lattice spacings, volumes and light-quark masses to estimate the baryon masses at the physical point [63][64][65][66][67][68][69][70][71][72][73][74][75]. We summarize the recent studies of several lattice groups in Table I. There is a remarkable agreement between the results of the different groups utilizing different types of quark actions and approaches to the physical point. Most of those studies are motivated by the observation of the Ξ cc baryon by LHCb and thus their focus has been on the lowest-lying positive parity baryons. Extracting the excited states, however, is a challenge compared to calculating the ground states. The majority of the attention has been on the light-quark sector, especially on the Roper resonance and the Λ(1405), while there are just a few groups that have studied the excited states of the charmed baryons.
The RQCD Collaboration reported results for the singly-and doubly-charmed baryons, including excited states [69]. They employ several 2+1-flavor gauge ensembles with a fixed lattice spacing but two different volumes and varying light-quark masses with the lightest one corresponding to a pion mass of m π ∼ 260 MeV/c 2 . All the sea and valance quarks (including the charm quark) are treated via a non-perturbatively improved stout-smeared Clover action. The bare charm-quark mass is tuned to reproduce the 1S spin-averaged charmonium mass. In addition to spectrum calculations, they also investigate the light-flavor dependence of the singly and doubly charmed states. To this end, the operator set they use consists of interpolating fields based on SU (4) symmetry and heavy quark effective theory (HQET) pictures. In order to access the excited states, they perform a variational analysis over a set of interpolating fields with three different quark-field smearings. Their chiral extrapolations follow a different approach compared to the other groups since they start from an SU (3) symmetric point for the light and strange quarks and vary their masses while keeping the singlet quark mass fixed in their descent to the physical point. This leads to fits based on Gell-Mann -Okubo relations for the charmed baryons. The lowest-lying extracted states are in good agreement with the other lattice determinations and with experimental values where available.
The Hadron Spectrum Collaboration (HSC) extracts the charmed baryon spectrum including positive and negative parity baryons with total spin up to J = 7/2. They use N f = 2 + 1 anisotropic lattices generated with a treelevel tadpole-improved Clover fermion action with a pion mass of m π = 391 MeV/c 2 . The anisotropic Clover action is used for the charm quark as well with its mass parameter tuned non-perturbatively so as to reproduce the dispersion relation for the η c meson. By using a large set of continuum interpolating operators, including nonlocal covariant derivative operators, subduced to the irreducible representations of the cubic group, they form the basis for the variational correlation matrix analysis and extract the spectrum of the singly-, doubly-and triplycharmed baryons [71][72][73][74][75]. Although the systematics are left unchecked and the pion mass is unphysical, their pioneering results provide valuable insight into the charmed baryon spectrum.
In this work, we follow a conventional approach by using local operators only. Notable improvements of this study compared to the previous works that extract the excited baryon spectrum are the fully relativistic treatment of the charm quark, thus suppressing the O(am Q ) discretization errors, and working on gauge configurations with almost physical light quarks, hence eliminating the chiral extrapolation systematics. We also perform variational analyses over sets of operators with different Dirac structures and quark smearings and their combinations. Preliminary results of this work have been presented in Ref [76]. This paper is structured as follows: we outline the approach to extract the baryon energies and the formulation of the variational analysis in Section II. Details of our lattice setup, the heavy quark action that we employ, and the choice of baryon operators are given in Section III. A detailed discussion on the variational analyses and the states we extract are presented in Section IV. Section V holds the summary of our findings.

II. EXTRACTING EXCITED STATES
For a given interpolator, χ i , the two-point correlation function contains the contributions from all the states that couple to the corresponding quantum number, where (E B ) B stands for the (energy of the) baryon state. The desired parity state can be isolated by applying the parity operator, P ± C ij (t) = 1 2 (1 ± γ 4 )C ij (t). Using a set of operators that couple to the same quantum numbers, one can utilize a variational approach to extract the tower of states. One can form an N × N correlation function matrix, where each element, C ij (t), is an individual correlation function given in Equation (1). Then, by solving the generalized eigenvalue problem [77,78], one extracts the left and right eigenvectors, ψ α and φ α , and uses them to diagonalize the correlation-function matrix, to access the energies of the states, E α . One can alternatively utilize the individual eigenvalues, λ α (t, t 0 ) ∼ e −Eα(t−t0) (1 + O(e −∆Eαt )), of the left and right eigenvalue equations given in Equation (3) to extract the energies of the states. Both approaches give complementary results with some caveats [79]. We prefer the method outlined above. Note that a suitable combination of the time-slice t 0 and the time slice of the eigenvectors, t , is chosen with respect to the quality and stability of the signal. Additionally, t may or may not be chosen equal to t. Once the correlation function matrix is diagonalized, one can follow the standard techniques and perform an effective mass analysis for each state, α,

A. Quark Actions
We employ the 32 3 × 64, 2 + 1-flavor gauge configurations that are generated by the PACS-CS Collaboration [80]. These configurations are generated with the Iwasaki gauge action (β = 1.9) and with the nonperturbatively O(a)-improved Wilson (Clover) action (c sw = 1.715) for the sea quarks. We perform our simulations on the κ sea ud = 0.13781 subset, which have almost physical light quarks corresponding to m π = 156(9) MeV/c 2 as measured by PACS-CS. The hopping parameter of the strange quark is fixed to κ sea s = 0.13640, while the lattice spacing is determined to be a = 0.0907(13) fm (a −1 = 2.176 GeV).
We use the Clover action for the valence u/d and s quarks. The hopping parameter of the valence light quarks is set equal to those of sea quarks, κ val u/d = κ sea u/d . Due to an overestimation of the mass of the Ω − particle with κ val s = κ sea s , however, we re-tune the hopping parameter of the valence strange quark to κ val s = 0.13665, in order to match the physical Ω − mass on these configurations. Details of this tuning are discussed in Ref. [81].
We employ a relativistic heavy-quark action for the charm quark, where the Ψs are the heavy quark spinors and the fermion matrix is given as with the free parameters r s , ν, c B and c E to be tuned in order to remove the discretization errors appropriately. We adopt the perturbative estimates r s = 1.1881607, c B = 1.9849139 and c E = 1.7819512 [82] and the nonperturbatively tuned ν = 1.1450511 value [66]. We re-tune the charm-quark hopping parameter to κ Q = 0.10954007 non-perturbatively so as to reproduce the relativistic dispersion relation for the 1S spin-averaged charmonium state. With these parameters, masses of the η c and the J/ψ are m ηc = 2.984(2) GeV/c 2 , m J/ψ = 3.099(4) GeV/c 2 . The hyperfine splitting is estimated as ∆E (V −P S) = 116(4) MeV/c 2 , in agreement with its experimental value. Further details of our charm quark tuning can be found in Ref. [81].

B. Baryon operators
The baryon operators that we employ are tabulated in Table II in a shorthand notation while the explicit forms of the operators can be found in Table III. Note that we do not distinguish between u and d quarks since they are degenerate in our lattice setup.
For the spin-1/2 baryon, we form three individual operators by using the Dirac structures, (see Table III). An explicit example for the N -like operator is The χ 4 -type operator with the Dirac structure [Γ 1 , Γ 2 ] = [γ 5 γ 4 , 1] corresponds to the time component of an operator with [Γ 1 , Γ 2 ] = [γ 5 γ µ , γ 5 ], which couples to both spin-1/2 and spin-3/2 particles. It has been shown that projecting out the spin-1/2 component of such an operator results in two terms: a linear combination of the χ 1 and the χ 2 , and a term containing the χ 4 operator [83]. Furthermore, the χ 4 -type operator is distinct from the χ 1 and the χ 2 from a chiral transformation perspective [84], making it a viable choice for the basis set of the spin-1/2 operators. Note that we limit ourselves to only one Dirac structure for the spin-3/2 baryons, which Among the operators discussed in this section, the ones coupling to the Ξ c and Ξ c states deserve special attention. The Ξ c (Ξ c ), which belongs to an SU (3) anti-triplet (sextet) is anti-symmetric (symmetric) with respect to the exchange of s and u/d quarks, which should hold for the respective operators. For Ξ c , this can be achieved by both N -like and Λ-like operators, which will both be used in this work. Note that our N -like Ξ c operator was referred to as "HQET" in Ref. [69]. For Ξ c , we employ a different operator combination with the correct symmetry properties as shown in Table III. While Ξ c and Ξ c states decouple in the SU (3) limit, they can in principle mix in our setup due to the breaking of the SU (3) symmetry. This mixing can be studied by computing cross-correlators of Ξ c and Ξ c operators. The results of such an analysis will be discussed in Section IV.

C. Simulation details
Quark fields of the interpolating operators are Gaussian smeared in a gauge-invariant manner at the source, (x, y, z, t) = (16a, 16a, 16a, 16a), for all the baryons with three different sets of smearing parameters, corresponding to an rms radius of ∼ 0.2, 0.4 and 0.7 fms for the quark wave-functions. Sink operators are smeared in the same manner. However, we find that the signal deteriorates rapidly with increasing sink operator smearing. For this reason, we analyze the spin-1/2 baryons with smeared-source -point-sink correlation functions with a fixed source smearing for all the quark fields. Correlation functions depend mildly to the smearing of the singly represented quarks and the plateau regions become independent of the smearing after a certain number of iterations. Therefore, we apply the smearing to the quark fields depending on their flavor and quantity. We treat the u-, d-and the s-quarks on an equal footing and consider them as light quarks in comparison to the charm quark. When the interpolating operator is formed by TABLE III. Interpolating operators with generic Dirac structures for spin-1/2 and spin-3/2 baryons. C = γ2γ4 is the charge conjugation operator. [Γ1, Γ2] choices and the quark contents are given in the text and in Table II. Spin Baryon Operator two light quarks and a charm quark, we fix the smearing of the charm quark to 0.7 fms, which is the widest of the smearings that we have, to decouple its effects and perform the variational analyses over the smearings of the remaining light quarks. Smearing parameters of the individual light quarks are set to be equal. This is true for all the baryon fields with the exceptions of Ω ( * ) cc , in which case the smearing of the strange quark is fixed to 0.7 fms and the smearings of the charm quarks are varied, and Ω ccc , for which the treatment is the same as light quarks. For the spin-3/2 baryons, we use smeared-source -smeared-sink correlators to form an operator basis from an operator with fixed Dirac structure. A discussion on the operator basis is given in Section IV A 1. Parity is selected by applying the parity projection operator, P ± , to the individual correlation functions.
We bin our data with a bin size of 15 measurements to account for the autocorrelations on this ensemble and estimate the statistical errors via a single elimination jackknife analysis. We performed our computations using a modified version of the Chroma software system [85] on CPU clusters along with the QUDA library [86,87] for the valence u−/d− and s-quark propagator inversions on GPUs. The charm quark inversions are done on CPUs.

A. Variational analysis
To obtain the individual states from a set of operators, one solves the generalized eigenvalue problem on each time slice, t, against a reference time-slice, t 0 , as discussed in Section II. To ensure the consistency of this step, it is necessary to check that the solutions are stable with respect to t 0 , since it can be chosen freely. Another concern is associating the eigenvalues with the states. Eigenvalues are sorted in increasing order on each time slice. However due to the faster deterioration of the higher states' signal, their eigenvalues fluctuate heavily as time evolves and can sometimes be smaller than the eigenvalue associated with the lower state. This situation might misguide the analysis if not addressed properly. In order to make sure that the eigenvalues are associated with the correct states, we fix the time-slice of the eigenvectors, t , that is used to diagonalize the correlation function matrix, to a specific value. This procedure, however, introduces an extra parameter dependence to the analysis. We check this dependence for each channel for a range of t values. The dependencies on t 0 and t can be tracked by investigating the respective eigenvectors, whose components should be stable when changing both fictitious time parameters. We illustrate such a consistency check in Figure 1. We perform this check for each channel and select a (t , t 0 ) combination that optimizes the signal quality. A common choice is t ≥ 2a.
1. Operator dependence a. Operator basis: Having three operators with differing Dirac structures, it is possible to analyze both the full 3 × 3 correlator, but also various combinations of 2 × 2 correlators. While the full information for all of them is contained in the 3 × 3 case, the 2 × 2 correlators can provide valuable and comprehensible information about which state couples to which operator. For this purpose we here investigate the correlators with different operator sets. We find that the variational analyses over two different sets of spin-1/2 operators, namely over {χ 1 , χ 4 } and {χ 1 , χ 2 }, give two distinct second eigenvalues for the positive parity states. The {χ 4 , χ 2 } set produces similar results to that of {χ 1 , χ 2 }. For negative parity, only the {χ 1 , χ 4 } combination yields mostly wellseparated second eigenvalues, whereas the second eigenvalues of the {χ 1 , χ 2 } and {χ 4 , χ 2 } bases lie closer to the first eigenvalues. When we extend the operator basis to the {χ 1 , χ 2 , χ 4 } set and solve the corresponding 3 × 3 variational system, the 2 × 2 results are reproduced. These findings are illustrated in Figure 2 for the positive and negative parity Ξ c , Ω c , and Ξ cc baryons where we show the fit results from a plateau approach. These representative baryons are chosen such that they correspond to the different operator characteristics, i.e. Λlike, singly-charmed N -like, and doubly-charmed N -like, respectively.
b. N -like operators: Although we use the same Nlike operators for the singly-charmed and the doublycharmed spin-1/2 baryons, it is reasonable to expect a different behavior when we solve the variational system, since they belong to different layers of the mixed-flavor c. Λ-like operators: Λ c and Ξ c belong to the totally flavor-antisymmetric SU (4) anti-quadruplet and hence are studied via the flavor-octet Λ-like operators. The behaviors of these operators depicted in Figures 2c and 3 show similarities to the N -like Ξ cc case. It can naively be expected that the first term of the Λ-like operator (see Table III) would have the dominant contribution, which would mean that it is in essence the same as the N -like operator. Indeed, by rearranging the latter two terms of the Λ-like operator via Fierz transformations, one can show that the coefficient of the q T a 1 (x)Cγ 5 q b 2 (x) q c 3 (x) term of the operator is five times the other resulting terms. The same argument holds for the other Dirac structures as well. This dominance is realized in our comparisons of the Ξ c ( 1 2 + ) results illustrated in Figure 4, where we have an almost identical signal for the ground states calculated via the Λ-like and the N -like operator.
Additionally, the flavor decomposition of the Λ c studied in Ref. [88] by three of the present authors shows that the negative parity Λ c baryon consists of a mixture of flavor-singlet and flavor-octet wave-functions. The flavor-octet interpolating operator that we employ for the Λ c baryon may therefore be inadequate to resolve the lowest-lying negative parity state by itself. A similar conclusion was reached in Ref. [89]. The first excited negative parity state on the other hand, is dominated by a flavor-octet wave-function and it is possible that this state is contaminating our lowest Λ c ( 1 2 − ) signal, which could be a plausible explanation of the apparent overestimation of its mass (see Table IV and Figure 8).
We analyze the Ξ c channel with two different types of operators. One being the Λ-like, the other the Nlike operator as given in Table III. We find that both give consistent results for the positive parity case while there is a difference for negative parity. As shown in Figure 2f, the N -like operator couples to a lower-lying state for the {χ 1 , χ 4 } basis. Similar differences between these operators for the negative parity sector have been reported by the RQCD Collaboration [69].
d. Ξ c − Ξ c mixing: We perform a correlation matrix analysis consisting of the Ξ c , and N -like and Λ-like Ξ c operators in order to investigate the possible mixing between these baryons. We construct the correlationfunction matrices for this analysis in two steps. First, we solve a variational system over the {χ 1 , χ 4 } basis for each element of the correlation matrix and take the lowest lying state. We find that this approach helps to isolate the ground states better. We then solve another 2×2 correlation matrix with both Ξ c and Ξ c ground state operators to investigate the mixing effects.
For positive parity Ξ c and Ξ c , we analyze the cross correlators between the flavor-octet SU (4) Ξ c -Ξ c , and the N -like Ξ c -Ξ c individually. We find that the Ξ c and Ξ c signals separate nicely, and the N -like and Λ-like Ξ c operators produce consistent signals with negligible mixing (see Figure 4). Magnitudes of the eigenvectors also con-  Table III. Note that we only have two smearings for that case. State numbering, α = 1 − 3, follows the notation of the 3 × 3 solutions even for the 2 × 2 solutions to emphasize the coupling of certain operators to certain states. All energies are extracted via a plateau method, see main text for a discussion.  firm that the Ξ c and Ξ c states have distinct signals. In case of negative parity, there appears to be non-negligible mixing between the two states dependent on the variational parameters. Specifically, the Λ-like Ξ c has a negligible Ξ c component, while the N -like Ξ c state has up to a 10% Ξ c mixing although the effect seems to depend on the variational parameters. The reason why the negative parity Λ-like operator gives signals close to the Ξ c is understood to be related to the overestimation of the mass obtained for that operator rather than a mixing effect. The Ξ c appears to have at most a mixing of 5% with the N -like Ξ c . In all, we see that for negative parity the mixing is not completely negligible, but nevertheless quite small.

Smearing dependence
a. Spin-1/2 baryons: We observe that, evidently, the ground state signals remain stable with respect to the smearing radius. The excited-state signals on the other hand show a clear dependence to the smearing radius of the source quark fields. This is readily visible for every case given in Figure 2. For both positive and negative parity, states that are clearly separated from the ground state tend to decrease as the smearing radius increases with no apparent plateau behavior. Note that all the energies are extracted via a plateau approach, which are dependent on the choice of the fit windows. Extracting the energies from two-exponential fits are more reliable for the ∼ 0.2 and 0.4 smearings, where those fit results coincide with that of ∼ 0.7 extracted via a plateau approach or a two-exponential fit. This indicates that the signals of the widest smearing are the most reliable to estimate the energy levels.
When we enlarge the operator basis by combining two operators with two different smearings and perform a 4×4 analysis, we end up with quite noisy solutions due to the current limited statistics, which renders a conclusive analysis impossible. We, however, observe an apparent degeneracy in three out of four solutions as shown in the b. Spin-3/2 baryons: We find that solving a 3 × 3 variational system with smeared-smeared operators only, provides no additional information compared to a 2 × 2 system with the smearings at hand. One solution turns out to be indistinguishable from the other so we focus on the solutions from the two narrower smearings, which give less noisier signals.

B. Charmed baryon spectrum
The energy levels from the diagonalized correlation functions are extracted by fitting the data to the form given in Equation (4). Additional exponential terms are employed to stabilize the fits against excited-state contributions. In most of the cases, where the signal forms a plateau in the effective-mass plots, masses of the lowest states extracted from the one-exponential fits agree with the multi-exponential fit results within their error bars. Yet, a two-exponential form stabilizes the fits and improves the accuracy of the results. This is especially true when analyzing the widest smearing case. The extracted energies are compiled in Table IV. Since we are at the isospin-symmetric point, m u = m d , our results should be understood as the isospin averaged masses of the respective states.
As we have discussed in Section IV A 2, a variational analysis over a set of different smearings for a fixed operator returns solution eigenvectors that couple to the widest smearing. Therefore, we always use an operator basis with quark smearings fixed to the widest one. For the spin-1/2 cases, we perform 3 × 3 variational analyses with a fixed smearing over the operator sets {χ 1 , χ 2 , χ 4 } and extract signals of three states for each channel. The third energy level with largest energy is however usually lost to noise already at relatively early time slices or decays to the ground states due to inaccuracies in the diagonalization procedure of Equations (3) and (4). For instance, in case of the positive parity spin-1/2 Ξ c baryons, we find that the state dominantly coupling to the χ 2 operator decays to the ground state signal before showing a plateau that may be a candidate signal for an excited state (blue rectangles in the top left plot of Figure 5). Signals of possible third states for the spin-1/2, positive parity Σ c , Ξ c and Ω c channels emerge in early time slices of effective mass analyses but are quickly lost to noise. It is usually possible to identify a fit region of 2-3 points for the narrowest smearing but we find the energy extracted via this approach to be unreliable, since the fit window is very small and the smearing dependency of the state cannot be established. Positive parity spin-1/2 Ξ cc and Ω cc signals mimic the behavior of Ξ c , where there appear signals one could potentially identify as distinct states. However we find that those states are rather unstable under the change of variational parameters. In addition, extracted energies are highly dependent on the extraction method -plateau approach or a two-exponential fit. Therefore, even though we show their signals in the plots, we do not extract or report any corresponding energy values.
In general, we find that the negative parity sector appears to be richer in comparison to the positive parity case. Indeed, we could identify three distinct states for most of the negative parity spin-1/2 channels. Isolating the low-lying states via a plateau approach is a challenge here since multiple energy levels appear in a narrow energy range. Two-exponential fits are very helpful in such cases to disentangle and extract the states more accurately. Effective mass plots illustrating the above discussions are given in Figure 5.
a. Mass differences: Hyperfine splittings, the mass differences between the spin-3/2 and spin-1/2 states, of the Σ c , Ξ c , and Ω c channels are reproduced in good agreement with the experimental values. Mass differences between the positive and negative parity states also agree well with the available experimental results. The first excited states of the positive parity baryons lie quite high, 400 MeV to 1 GeV, above the ground states. A common pattern is that, more than one negative parity state for the singly-and doubly-charmed spin-1/2 baryons appear in between the positive parity ground and first-excited state. The first two negative parity states of the Σ c , Ξ c , Ω c , Ξ cc , and Ω cc channels lie close to each other. The splittings between those states are smaller for the Ω c and Ω cc baryons compared to those of Σ c , Ξ c , and Ξ cc . The situation is different for the Λ c and the Ξ c baryons where the negative parity states are roughly 300 MeV apart.
b. Scattering states: It is essential to examine the relevant thresholds for the negative parity states in order to check if they could correspond to scattering states. It is possible for the negative parity ground states to couple to the S-or D-wave scattering states of a positive parity baryon and a negative parity meson. The relevant thresholds which respect to isospin, spin, parity, strangeness and charm quantum numbers are, We plot the above two-particle thresholds together with the extracted negative parity energies in Figure 6. The two-particle scattering energies are calculated via E = M 2 1 + p 2 1 + M 2 2 + p 2 2 , where M i is the mass of the particle and p i = 2πn/L the lattice momentum. We use the π mass quoted in the PACS-CS paper [80] and the experimental K mass, since we use a strange quark mass re-tuned to its physical value via the K mass input [89], along with the positive parity baryon masses from Table IV of this work in calculating the threshold energies. The Λ + D threshold has to be estimated differently since we do not calculate the Λ baryon or the D meson in this work. In estimating the threshold, we take the experimental Λ mass and multiply it by a correction factor, Λ our c /Λ exp c , due to our overestimation of the Λ c mass. The uncertainty of this value is assumed to be same as that of Λ our c . The D meson mass is taken to be its experimental value with its uncertainty neglected.
The momenta p 1 and p 2 are set to zero. An inspection of Figure 6 shows that some of the Ξ c baryon signals may contain scattering states because of their vicinity to various thresholds. Indeed, lie close to at least one related threshold.
We also find some states that lie above the thresholds to be close to their respective boosted (n > 0) thresholds. MeV/c 2 , 40 MeV/c 2 lower than that of BABAR. It is noted in the PDG listings that the state that has been observed by BABAR might be a different Σ c excitation.
Given that these states have been seen in the Λ c π invariant mass spectra, a straightforward assignment for the quantum numbers would be J P = 1/2 − . From a quark model perspective (see paragraph f.), there are three possible low-lying negative parity spin-1/2 Σ c excitations. Two λ-modes with diquark spin j = 0 and j = 1, and a ρ-mode with diquark spin j = 1. In the heavy quark limit, the S-wave Σ c (2800) → Λ c π transitions of the j = 1 λ-and ρ-modes would be forbidden due to the violation of the spin-parity conservation of the light-quark degrees of freedom. A heavy quark effective theory calculation estimates a very large decay width, of the order of 885 MeV, for the j = 0 λmode [56], which rules out the 1/2 − quantum number for Σ c (2800). On the other hand, a D-wave transition is possible and points to the J P = 3/2 − , 5/2 − possibil-      Note that the three extracted negative parity Σ c states are well above their respective two-particle thresholds so that the two-particle contribution to the signals should be suppressed.
d. Excited Ξ c and Ξ c states: The experimental spectrum of the Ξ c and Ξ c channels consists first of the respective J P = 1/2 + ground states, and the first Ξ c ( 1 2 − ) excited state, which are all experimentally well established and which we reproduce well in our work. The energy levels above the lowest three are less well established, both experimentally and theoretically. Above 2.9 GeV/c 2 , the PDG reports the five states Ξ c (2930), Ξ c (2970), Ξ c (3055), Ξ c (3080) and Ξ c (3123), for none of which the spin and parity quantum numbers have been measured. Very recently, the spectrum of these states has received an update by a new measurement of the  LHCb Collaboration [12] in the Λ + c K − channel. According to this measurement, the Ξ c (2930) (observed earlier by the Belle [4] and the BABAR [10] Collaborations in the same channel) should be considered to be a previously unresolved combination of two independent states Ξ c (2923) and Ξ c (2939). The third observed state in Ref. [12], Ξ c (2965), corresponds either to the already seen Ξ c (2970), or is another entirely new resonance.
Let us discuss potential interpretations of our findings with regard to this rather rich experimental spectrum. We find two negative parity spin-1/2 Ξ c states in the vicinity of the lowest three (or four) states above 2.9 GeV/c 2 , Ξ c (2923), Ξ c (2939), Ξ c (2965) and potentially Ξ c (2970), which suggests that such quantum numbers can be assigned to at least two of these states. While our numerical results are not precise enough to draw any firm conclusions, our obtained spectrum is most naturally interpreted as either Ξ c (2923) or Ξ c (2939) and similarly Ξ c (2965) or Ξ c (2970) being a Ξ c ( 1 2 − ) state.
The already known Ξ c (2970) state has been observed in the Λ c Kπ channel -also proceeding approximately half of the time via the intermediate Σ c (2455)K channel -and in the Ξ c π, and Ξ c (2645)π channels by the Belle [5][6][7] and BABAR [11] Collaborations. These decay channels imply several possible quantum numbers, J P = (1/2 ± , 3/2 ± , 5/2 ± ), for this state, which is not in contradiction with the above potential assignment. The Ξ c (3055) was observed by the Belle and the BABAR Collaborations in the Σ c K channel [8,11] and in the ΛD channel only by the Belle Collaboration [9]. Finally, the Ξ c (3080) was reported by the Belle Collaboration [9] in the Σ c K, Σ * c K, and ΛD channels and by the BABAR Collaboration [11] in the Λ c Kπ channel via the Σ c (2455)K channel. Similar to the Ξ c (2970) case, these decay channels suggest several quantum numbers, such as J P = (1/2 ± , 3/2 ± , 5/2 ± ). Our second Ξ c ( 1 2 + ) state appears to be the most probable candidate for this resonance.
e. Excited Ω c states: The five new excited Ω 0 c states reported by the LHCb Collaboration [13] were seen in the Ξ c K channel. One would hence naively expect these states to have negative parity. A first dedicated lattice QCD calculation has confirmed this expectation by assigning negative parity to these states [57], with total spin ranging from J = 1/2 to 5/2. The two Ω c ( 1 We should reiterate that since we only employ local three-quark operators, we are limited in our ability to resolve all molecular, radial or orbital excitation modes of the higher lying states. Our results should hence be considered as indicative in identifying potential compact three-quark states among the experimentally observed energy levels in the Ξ c and the Ω c channels. Conversely, the levels that we are not able to reproduce, could be candidates for molecular or orbitally excited states. It is however at present too early to assign definite quantum numbers without a through scattering state analysis since some of our negative parity states lie close to the thresholds.
The values in Table IV are illustrated in Figure 7 together with the relevant experimental results. The latest Ξ c results from the LHCb Collaboration are shown as well. The similarities between the Λ c and Ξ c , and Σ c , Ξ c and Ω c are evident as expected from their flavor structures.
f. Interpretation from a quark model perspective The quark model (QM) has has been useful in giving a pictorial and intuitive interpretation of the mass spectrum obtained by lattice QCD computations. The QM derives the energy and structure of a system by considering constituent valence quarks and their interactions. For the excited states, in particular, it can clarify what the essential degrees of freedom in a specific excitation are.
For heavy-quark baryons, the heavy-quark spin symmetry plays an important role. As the coupling of a heavy quark to the magnetic component of gluons is suppressed by a 1/m Q factor, the heavy quark spin is approximately conserved. For singly charmed baryons, this symmetry is manifested by the appearance of heavy-quark spin doublets, in which spin (j −1/2, j +1/2) pair states approach each other with increasing quark masses. Here, j represents the total spin minus the heavy quark spin of the considered baryon.
We will here briefly compare the present lattice QCD results with the QM predictions and study how the essential excitation modes arise in the spectrum. Quite remarkably, multiple features of the QM predictions are confirmed in the obtained lattice QCD spectrum of the charmed baryons.
1. Our lattice QCD results for the positive parity "ground" states agree completely with the QM assignments, i.e., the spin, parity, isospin and flavor representation, and the mass orderings are consistent. The QM predictions for the splitting between the spin 1/2 and 3/2 states are also in quantitative agreement with the obtained lattice results.
2. Among the positive parity ground states, Ξ c is most interesting, because it contains three different valence quarks, c, s, and u/d. In the QM, the total spin of s and u/d can take either S = 0 (Ξ c ), or 1 (Ξ c ). The existence of two low-lying positive parity states is indeed realized in lattice QCD as well as in experiment. In the QM, the distinction of Ξ c and Ξ c is guaranteed by the flavor SU (3) symmetry, while the SU (3) breaking with m s = m u/d will mix the two Ξ c 's. The QM predicts, however, that the mixing is suppressed for the ground state due to the heavy quark spin symmetry, which is confirmed in our lattice QCD results.
3. Low-lying negative parity singly charmed baryons are described in the QM as orbital P -wave excitations. They are categorized in two classes, λ-mode and ρ-mode [17,28]. The λ-mode is characterized by the P -wave excitation between the charm quark and the center of mass of the light quarks, while the ρ-mode is given by the excitation between the light quarks. The QM predicts that the λ modes are lighter than the ρ-modes for singly heavy baryons.
The QM spectrum depends on the flavor structure: For the flavor anti-triplet Λ c and Ξ c , we find a set of (1/2 − , 3/2 − ) states in the λ-mode, and (1/2 − ), (1/2 − , 3/2 − ) and (3/2 − , 5/2 − ) states in the ρ-mode. Thus, among the three 1/2 − states, the QM predicts that one λ-mode state is lighter than the other two. This structure is indeed seen in the Λ c and Ξ c spectrum given in Table IV and Figure 7. The next 1/2 − state is about 300 MeV higher, which can be regarded as the mass splitting between the λ-and ρ-mode states.
On the other hand, the flavor 6 baryons, Σ c , Ξ c and Ω c , have two λ-mode 1/2 − states, one of them being accompanied by a 3/2 − state. In terms of the heavy-quark spin symmetry, we have a (1/2 − , 3/2 − ) spin doublet and an isolated singlet 1/2 − . The lower two λ-mode states come close in energy, but can be distinguished by the total angular momentum of the light-quark system. Thus we expect two 1/2 − and one 3/2 − states as the lowest negative parity excitations for Σ c , Ξ c and Ω c . One sees that, indeed, these three states turn out to be almost degenerate in the lattice QCD spectrum of these channels in Table IV and Figure 7. Other states are much higher in energy, which again confirms the predicted QM assignments.
In all, the low-lying spectra of both the positive and negative parity charmed baryons confirm the effectiveness of the QM in assigning the quantum numbers and symmetry properties of heavy baryons. g. Comparison to other lattice results: We compare our results to other lattice determinations and experimental values in Figure 8. Our positive parity ground states are in good agreement with the experimental results and the calculations of the other lattice groups with the exception of the Λ c , which is overestimated in our work. Taken altogether, this is a good indication that we are close to the physical point. The first excited positive parity states also mostly agree with the predictions of the HSC [71,72] and the RQCD Collaboration [69]. For negative parity, there are notable differences between our and RQCD's results, especially for the doubly-charmed baryons. For the excited states of the Ξ cc and Ω cc , there are discrepancies between our extracted spectrum and that of RQCD, while our results are similar to those obtained by the HSC [72]. Although we do not show the corresponding HSC spectrum in Figure 8, the pattern they extract in their preliminary studies for the negative parity spin-1/2 singly-charmed baryons [75] is similar to our results as well. Such a qualitative agreement for the low-lying spectrum is quite encouraging since, in contrast to the HSC, which utilizes both local and non-local operators, we only use local operators.

V. SUMMARY AND CONCLUSIONS
We have calculated the ground and the first few excited states of the charmed baryons on 2+1-flavor gauge configurations with a pion mass of ∼ 156 MeV/c 2 . The charm quark is treated relativistically by employing a relativistic heavy-quark action to remove O(am Q ) discretization errors. The states are extracted via a variational approach over a set of interpolating fields with different Dirac structures and quark-field smearings. By performing separate variational analyses with multiple subsets of the operator basis, we have studied the Dirac-structure and smearing dependence of the excited states. Our results indicate that the excited-state signals are highly susceptible to the width of the quark smearing. Additionally, solutions of a variational analysis over a set of smeared operators with fixed Dirac structure couple dominantly to the operator that is smeared the widest within our employed smearing parameter range. These results highlight the importance of forming the variational basis from different Dirac structures since relying on smeared operators only might miss some parts of the spectrum.
In comparing the operator dependence of the extracted positive and negative parity states, we have extended the SU (4) operator basis of the Ξ c baryons to include not only Λ-like, but also N -like operators. Both operators give consistent results for the positive parity case while there appears a difference for the negative parity states. We have also investigated the Ξ c -Ξ c mixing by studying the cross-correlators of this system.
Our masses of the low-lying states agree well with the available experimental results and previous lattice determinations. Consequently, the hyperfine splittings and the mass differences between the positive and negative parity states are reproduced, which is a good check of the relativistic action we employ for the charm quark. Excited states in the positive parity channel lie 400 MeV to 1 GeV above the ground states depending on the quantum numbers. One or more negative parity states appear in between. This pattern is consonant with the QM expectations. Although we identify several states that are close to observed excited Σ c , Ξ c and Ω c baryons, mostly in the negative parity channels, some of the signals are in close proximity to the related two-particle thresholds. Without a thorough scattering state analysis with multiple volumes and two-particle operators, the contamination from the thresholds remain unidentified.
From a qualitative point of view, the spectrum we extract is similar to what has been reported by the Hadron Spectrum Collaboration (HSC). This is quite encourag-