Spectrum of QCD with one flavor: A window for supersymmetric dynamics

We compute the spectrum of the low-lying mesonic states with vector, scalar and pseudoscalar quantum numbers in QCD with one flavour. With three colours the fundamental and the two-index anti-symmetric representations of the gauge group coincide. The latter is an orientifold theory that maps into the bosonic sector of $\mathcal{N} = 1$ super Yang-Mills theory in the large number of colours limit. We employ Wilson fermions along with tree-level improvement in the gluonic and fermionic parts of the action. In this setup the Dirac operator can develop real negative eigenvalues. We therefore perform a detailed study in order to identify configurations where the fermion determinant is negative and eventually reweight them. We finally compare results with effective field theory predictions valid in the large $N_C$ limit and find reasonably consistent values despite $N_C$ being only three. Additionally,the spin-one sector provides a novel window for supersymmetric dynamics.


I. INTRODUCTION
Understanding the dynamics of strongly coupled gauge theories, such as QCD, has motivated the construction of several expansions complementary to the standard, perturbative, weak coupling expansion. One of the most prominent examples is the large N C limit (where N C is the number of colours), introduced by 't Hooft in Ref. [1]. In this case one keeps quarks in the fundamental representation of the gauge group SU (N C ) and organises an expansion in 1/N C using a diagrammatic approach. Several properties of QCD can then be understood in a simple way, suggesting that N C = 3 is "large". However, since quark loops are suppressed in this expansion, the properties of the η ′ -meson are not well reproduced in the 't Hooft large N C limit. Baryons also become increasingly heavy as N C grows.
Partly motivated by that, Corrigan and Ramond (CR) introduced a different large N C expansion in Ref. [2], in which quarks transform according to the two-index antisymmetric representation of the gauge group. While 't Hooft and CR expansions coincide for N C = 3, they are very different in the large N C limit. Notably, in the CR expansion, quark loops are not suppressed as N C → ∞. A simple scaling of the dimensionality of the representations of the quark fields suggests that the CR large N C limit may share non-trivial dynamical properties with supersymmetric theories. This relation has been made precise by Armoni, Shifman and Veneziano in Refs. [3,4], where a connection between the mesonic sectors of the two-index (anti-)symmetric theories and of N = 1 super Yang-Mills theory (sYM) is established. The subtle issues of the confinement properties and (in)equivalences at large N C were investigated in Ref. [5]. Further developing the correspondence, in Ref. [6] supersymmetry inspired effective Lagrangians have been constructed for gauge theories featuring one Dirac fermion transforming either in the symmetric or in the anti-symmetric twoindex representation of the gauge group SU (N C ) (orientifold theories). At leading order in the 1/N C expansion such effective theories coincide with that of supersymmetric gluodynamics restricted to its mesonic sector. These correspondences imply that non-perturbative quantities computed in orientifold theories can be related, up to 1/N C effects, to the analogous ones in sYM. By considering 1/N C supersymmetry breaking effects, including the explicit ones due to a finite quark mass, a number of predictions are made in Ref. [6] concerning the spectrum of the low-lying mesonic states. 1 In this work we confront such predictions with non-perturbative results produced by means of lattice simulations. For simplicity, in this first study we only consider N C = 3, which corresponds to one-flavour QCD. This has the advantage that available simulation packages for lattice QCD can be used without having to develop new code for handling representations of the fermionic fields different from the fundamental one. Future studies will be devoted to the extension to N C > 3. Intriguingly, by flipping the point of view (cf. Refs. [5,8]), we can use QCD results to learn about the spectrum and dynamics of supersymmetric theories, in particular N = 1 sYM. Analytic and numerical studies can now be employed to investigate several dynamical properties, including the theta-angle [6].
One-flavour QCD has been the object of several previous lattice studies. The qualitative behaviour of the theory has been discussed in Ref. [9]. In Ref. [10] the quark condensate has been computed by comparing the density of low-lying eigenvalues of the overlap Dirac operator to predictions from Random Matrix Theory [11,12]. The result is consistent with the prediction for the gluino condensate in sYM obtained in Ref. [13]. Using Wilson fermions, Ref. [14] presents a computation of the lowlying hadronic spectrum of one-flavour QCD. We improve here on that computation by considering a finer lattice spacing, larger volumes and a tree-level improved fermionic action. In Ref. [15] the one-flavour SU (2) vector gauge theory with the fermion in the fundamental representation is studied as a possible composite model for Dark Matter (DM). The Dirac operator is discretised using Wilson's regularisation. The fundamental representation of SU (2) is pseudo-real making the global symmetries and dynamics different from three colours QCD. In particular, the dark-matter model of Ref. [15] features a mass-gap with vector mesons being the lightest triplet of the enhanced SU (2) global symmetry. A similar DM model based on SU (2) gauge theory with scalar quarks was proposed in Ref. [16].
Finally, in Ref. [17] the single flavour SU (2) theory is considered with the fermion in the adjoint representation. The goal in this case is to gain insights on the emergence of the conformal window. Again the Wilson Dirac operator is used in the numerical simulations. As is highlighted by this brief review, one-flavour QCD is implemented on the lattice by adopting either overlap (or more generally Ginsparg-Wilson) or Wilson fermions. That is because in those cases the single-flavour lattice Dirac operator can 1 A string theory dual of orientifold theories has also been used in Ref. [7] to make predictions in the massless limit.
be rigorously defined. Wilson fermions are computationally cheaper but in such regularisation the spectrum of the Dirac operator may contain real negative eigenvalues for positive (but small) quark masses. That might cause a sign problem as the fermion determinant may become negative on some configurations. Following Refs. [18][19][20][21] we discuss in detail how we monitor such cases. 2 Earlier numerical investigations of orientifold theories [23,24] used the quenched approximation where the sign problem is absent.
Directly simulating supersymmetric gauge theories on the lattice has been an active research field for many years. Since the literature is vast we refer the reader to the recent review in Ref. [25] and references within for a discussion of the status and open problems.
A preliminary account of the results we present in this paper appeared in Refs. [26,27]. The latter in particular contains some algorithmic exploratory studies for N C = 4, 5 and 6.
The remainder of this paper is organised as follows. In Section II we describe our computational setup and provide algorithmic details. In Section III we investigate the consequences of the sign problem in our simulations. In Section IV we report on the correlation function fits required to extract the spectrum at non-zero quark masses, before extrapolating the meson spectrum to vanishing quark masses in Sec. V. Finally, in Section VI we confront the effective field theory predictions with our results and provide an outlook.

II. SIMULATION SETUP
For the gauge part of the action, we employ the Symanzik improved gauge action [28] with a fixed value for the gauge coupling of β = 4.5. As fermion action we use one flavour of tree-level improved Wilson fermions [29] and set the parameter of the clover term to 1. The Wilson-Dirac operator D in clover improved form is defined as follows where a is the lattice spacing, m 0 is the bare quark mass and ∇ ( * ) ν denotes the covariant forward (backward) derivative. The hopping parameter κ is related to the bare mass m 0 by 1/κ = 2(am 0 + 4). In order to map out the relevant parameter space we generated 19 gauge field ensembles covering different hopping parameters κ between 0.1350 and 0.1410 and volumes ranging from 12 3 × 64 to 32 3 × 64. An overview of the simulation parameters can be found in Table I. We measure the topological charge Q by integrating the Wilson flow [30] using a third-order Runge-Kutta scheme with a step-size of ϵ = 0.01 and 1600 integration steps. The topological charge at the largest flow time (t/a 2 = 16) is shown for all ensembles in Fig. 19 in Appendix A. The topological charge behaves as expected: its distribution is narrower for lighter quark masses and broader for larger volumes [11]. The Wilson flow further allows us to estimate the lattice spacing (via the reference flow scale t 0 ) by studying the Yang-Mills gauge action density as a function of flow-time [30]. Since our goal is to determine dimensionless quantities, we only quote the lattice spacing in order to enable qualitative comparison with other lattice calculations. As there is no reference scale for a single flavour (N f = 1), we use the average of t 0 from N f = 0 [30] and N f = 2 [31] as an estimate for the lattice spacing with N f = 1. In practice, we use a value of √ 8t 0 = 0.45 fm. This allows us to obtain an indicative value for the lattice spacing of a ≈ 0.06 fm.
All configurations are generated using the openQCD software package [32]. Since we only simulate a single fermion in the sea, it is necessary to use the rational hybrid Monte Carlo (RHMC) algorithm [33]. In the rational approximation we adopt a Zolotarev functional of degree 10. In the absence of prior knowledge about the optimal Zolotarev approximation -in particular for just one In comparison with Ref. [21] this is a rather loose approximation, which is relevant for the tunnelling between regions of configuration space with positive and negative determinants of the Dirac operator. In addition, we include frequency splitting, i.e. we factorise the Zolotarev rational into two terms, where the first factor contains the poles 1 to 5 and the second term the contribution from poles 6 to 10. Throughout the entire generation, we adopt three levels of integration schemes. The outermost employs a second-order Omelyan integrator [34] with λ = 1/6, which is used for the contributions from poles 6 to 10. For the inner two levels we use fourth-order Omelyan integrators, where the remaining fermion force is calculated in the second, and the gauge forces in the innermost level. We tune the number of fermion integration steps (ML steps) in the different levels to achieve a high acceptance (between 84% and 99.9%, c.f. Table I).
The pseudofermion actions and forces are obtained using a simple multi-shift conjugate gradient solver. For ensembles with a lighter quark mass, i.e. with larger values of κ, we take advantage of the deflated SAP [35,36] preconditioned solver given in the openQCD framework. The trajectory lengths of our ensembles are typically between 2 and 3 molecular dynamic (MD) units. In our analysis, we use every 32nd (or 40th) trajectory, which implies that configurations are at least 64 MD units apart from each other. For each ensemble the resulting number of configurations N config on which we perform all measurements is listed in Table I. To increase the amount of statistics and to utilise smaller computing resources more efficiently, we branch our simulation stream into multiple replicas after thermalisation is reached.
Since the Zolotarev approximation in the RHMC is not exact, we correct our observables by using a reweighting scheme. To achieve this, on each configuration we compute four estimators for the reweighting factors w i using code from the openQCD package. The correctly reweighted gauge average of an observable O is then given by where we define w ′ = w/ ⟨w⟩. Figure 1 shows these normalised reweighting factors w ′ as a function of the trajectory length (excluding any thermalisation times) for two representative ensembles (L/a = 32, κ = 0.1390 and L/a = 24, κ = 0.1405). In Fig. 2 we show the variation of the reweighting factors for all ensembles and observe that the fluctuations increase with volume, but are insensitive to the quark mass.
As the phase space of this theory in the regularisation we have chosen is a priori unknown, we computed the trace of the Polyakov loop. We find that the Polyakov loop vanishes within errors on each ensemble, which indicates that we are simulating in the confined phase.

III. EIGENVALUE ANALYSIS
The use of Wilson fermions for lattice QCD with an odd number of quark flavours or with non-massdegenerate (light) quarks can introduce a sign problem. This occurs because the configuration space is divided into two sectors, one associated to a positive sign of the fermion determinant and one to a negative sign. These sectors are separated by a zero of the fermionic measure. Note that the latter translates into a pole of the fermionic force in the molecular dynamics algorithm. With exact integration and an exact expression for the square root function, the negative sector cannot be reached from the positive one. In practice the algorithmic choices for the rational approximation yield a finite (rather than infinite) barrier between the two sectors.
In the thermodynamic and continuum limit the trajectory is expected to be constrained to the positive sector. However, at finite volume, the presence of the negative sector has to be accounted for by sign reweighting which requires knowledge of the sign of the fermion determinant det(D). A direct computation is numerically (prohibitively) expensive. Instead we follow a strategy in which the sign of det(D) is inferred from computing a few of the lowest eigenvalues of the Dirac operator. This can be achieved at a cost linear in the lattice volume and using the approach we will now sketch: Due to γ 5 -Hermiticity of the Wilson-Dirac operator, i.e.
the matrix Q = γ 5 D is Hermitian and its spectrum is real. Furthermore, it holds that det(D) = det(γ 5 ) det(D) = det(Q) and that a zero eigenvalue of D is also a zero eigenvalue of Q. Recalling that the eigenvalues of D come in complex conjugate pairs, for det(D) to be negative there must be an odd number of negative real eigenvalues of D.
Since the fermion determinant det(D) is assumed to be positive for large quark masses, we can infer that the determinant at the unitary mass m * 0 , used in the actual simulation, is negative if and only if there is an odd number of eigenvalues that cross zero as the mass is decreased from large quark masses to m * 0 . The idea is to locate (on each gauge configuration) the largest value m t of the quark mass such that Q(m t ), and therefore D(m t ), has a zero eigenvalue. If m * 0 is larger than this value m t then we need to determine the number of zero crossings of the lowest eigenvalue(s) λ(m 0 ) of Q(m 0 ) by varying the bare mass m 0 from above m t down to m * 0 . To that end we combine the PRIMME package with openQCD as mentioned in Ref. [21].
In practice we proceed in two steps: First we perform a preselection to identify potential candidate configurations with a negative fermion determinant and for this subset of configurations we perform a tracking analysis to identify the configurations that indeed display a negative fermion determinant.
We start the preselection by measuring the lowest O(10) eigenpairs (λ i , ψ i )(m * 0 ) and their chiralities χ i (m * 0 ), defined by where the last equality follows from the Feynman-Hellman theorem [20,21]. The chirality hence corresponds to the slope of the eigenvalue function. This allows to categorise the eigenvalues of Q into those which approach zero as m 0 is increased and those which move away from it. In Figure 3 we plot the results of the eigenvalue-chirality analysis for the four lowest lying eigenvalues of the two L/a = 16 ensembles with κ = 0.1405 (left) and κ = 0.1410 (right). If a datapoint falls into the north-east or south-west quadrant, the eigenvalue moves further away from zero when the quark mass is increased, implying that there is no zero crossing for values larger than m * 0 . This is the case for all configurations with κ = 0.1405. Conversely, if a datapoint falls into the north-west or south-east quadrant this implies that the eigenvalue approaches zero as the quark mass is increased and a zero crossing is possible. Configurations with eigenvalues which display this feature can potentially have a negative determinant and therefore require further monitoring. As can be seen in Fig. 3, on the κ = 0.1410 ensembles we find a small number of these cases for which the second step, the tracking analysis, is performed. We find that datapoints close to the horizontal axis but in the "safe" quadrants, tend to occur only for comparably large values of |λ|. Under the assumption that the chirality changes slowly in the range of masses explored, even if the sign of χ were to change, the corresponding eigenvalues are not expected to be at risk of changing sign. This assumption is justified a posteriori in the tracking analysis.
On the configurations that displayed datapoints in the north-west or south-east quadrants we now measure the lowest 20 eigenpairs for several partially quenched masses around m * 0 . The eigenvalue functions λ i (m 0 ) and the eigenbasis {ψ i } are assumed to vary slowly and continuously with m 0 . Assuming that the different partially quenched masses are sufficiently close to each other it is possible to track how a particular eigenvalue behaves as a function of the quark mass as follows. For each set of neighbouring masses m 0 and m 0 + ∆m 0 we construct the matrix M ij = ⟨ψ i (m 0 )| ψ j (m 0 + ∆m 0 )⟩ of scalar products between the ith eigenvector ψ i (m 0 ) at m 0 and the jth eigenvector ψ j (m 0 + ∆m 0 ) at m 0 + ∆m 0 . We determine the largest entry M ij and interpret this to mean that the eigenvalue i at m 0 evolves to be the eigenvalue j at m 0 + ∆m 0 . We then remove row i and column j from the matrix and iterate the procedure until each eigenpair at m 0 has been assigned a corresponding eigenpair at m 0 + ∆m 0 . Figure 4 displays a configuration of the L/a = 16 and κ = 0.1410 ensemble where a negative determinant was detected. We observe that the line connecting the red downward facing triangles does cross zero as the mass m 0 is increased from m * 0 (highlighted as the vertical dashed line). Since there is only a single eigenvalue crossing zero in the region m 0 > m * 0 , we conclude that the fermion determinant is negative on this particular configuration.
We see from the representative example shown in Figure 4 that the assumption discussed above is indeed valid and the derivatives of the eigenvalues change very little in the range of masses explored. We also see that such derivatives are either of O(1) or small. This is expected and in agreement with the discussion in Ref. [37], where approximate relations are derived between the eigenpairs corresponding to small eigenvalues of Q and those of D.
The chiralities are expected to be significantly different from zero (and in that case close to ±1) only for eigenvectors corresponding to almost real eigenvalues of D.
We performed the above analysis for the two smallest values of the quark mass corresponding to κ = 0.1405 and 0.1410 for which we each have a L/a = 16 and a L/a = 24 ensemble. As discussed above (cf. left panel in Fig. 3) we did not observe any cases of a negative determinant for κ = 0.1405 on either of the two available volumes. Since negative eigenvalues are expected to have a higher likelihood to occur at small quark masses, we did not perform this analysis for any of the remaining larger masses. At κ = 0.1410 we found 6 configurations with a negative determinant for each of the two volumes. Furthermore, we observed that the negative sector is visited at most for the Monte Carlo time corresponding to two consecutive measurements. This might be related to our choice of parameters for the rational approximation of √ D † D yielding a relatively low barrier between the two sectors. We conclude that in our computational setup the sign problem for N f = 1 QCD is mild and the relative frequency of a negative determinant of the Dirac matrix is at the sub-percent level.

IV. CORRELATOR ANALYSIS
In order to obtain the spectrum of one-flavour QCD, we create mesonic correlation functions for states with a variety of quantum numbers. We are particularly interested in states with scalar (S), pseudoscalar (P) and vector (I) quantum numbers. We employ the Laplacian Heaviside (LapH) method [38,39] which allows us to efficiently compute quark-line disconnected contributions that appear in the computation of mesonic quantities with a single flavour.

A. Construction of correlation functions
Following Ref. [38] and, where possible, using the same notation we compute the N v lowest eigenpairs (λ i , v i ) of the three-dimensional gauge-covariant Laplacian using a stout smeared gauge field. On each time slice t we arrange these eigenvectors into a matrix V s as from which we then define the Hermitian smearing matrix as a function of the number of eigenpairs that were computed as Using a low number of eigenpairs corresponds to a broad smearing profile, whereas using a large number of eigenpairs corresponds to "less" smearing and taking the limit of all eigenpairs recovers the identity. Quark lines Q are  computed as for y i β (t). This is done for each eigenvector v i (i = 1, . . . , N v ), each spin component (α = 1, . . . , 4) and each time slice (t 0 = 0, . . . , T − 1), amounting to N t × N v × 4 inversions per configuration.
In our simulation, we keep the number of eigenvalues N v = 20 fixed for all ensembles. However, from these inversions we can construct operators which use fewer than 20 eigenvalues by truncating the elements of the square matrix V † s (t)(γ 4 D) −1 V s (t 0 ). Using this we compute meson correlation functions for N v ∈ {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 15, 17, 20}, which describe the same spectrum but have different smearing functions.
In all three channels (P, S, I), we use the appropriate interpolation operator (P, S, I) in the finite volume irreducible representation. For the S-channel we additionally construct a purely gluonic operator G [40] which induces the same quantum numbers as the S operator 3 . We consider all mutual combinations of G and S in the 'scalar-glue' system.

B. Reweighting and vacuum expectation value subtraction
The vacuum subtracted correlation function C X Y can be derived from the un-subtracted correlation function C raw X Y (t) and the vacuum expectation values (vevs) v X and v Y as where ⟨·⟩ denotes the gauge average. Whilst the vev is exactly zero for the P operator and numerically zero for the I operator, it is sizable for the S and G operators. We find that the statistical signal for correlation functions including G or S deteriorates when reweighting (cf. Sec. II) is combined with the naive vacuum expectation value subtraction defined in Eq. (4.5). This is due to the fact that the product of the vevs is many orders of magnitude larger than the exponentially decaying part of the correlator and as a consequence even the little noise introduced by the reweighting factors destroys the signal for the latter almost completely.
Since the vacuum expectation value is timeindependent, an alternative way to perform the vev subtraction is to take the temporal derivative of the unsubtracted correlation function. We find that this results in a significantly better signal when combined with reweighting and are therefore utilising this. We observe that only for the earliest time slices the uncertainty of the reweighted data is limited by the accuracy of the reweighting factors.

C. Correlation function fits
For a given channel (P, I or S), the correlation function C of operators O n X with X ∈ {S, P, I, G} using n eigenvalues can be approximated by the first N states We emphasise that the induced masses m X i depend on the channel X, rather than the specific operator X , in particular all combinations of S, G induce the same spectrum m S i . We extract the three lowest-lying states of the spectrum by performing simultaneous correlated fits to the symmetrised correlation functions C n X Y (t) for several choices of n (between 2 and 4). We illustrate two such fits for the example of the vector channel in Fig. 6. We defer the discussion on the slow approach to the ground state for the bottom panel to Sec. V B. In order to assess systematic uncertainties associated with the choice of smearing radii, we vary which n enter into a particular fit. In particular, for the vector and pseudoscalar channels we perform three different fits, simultaneously fitting N v = (20,12,6), (17,10,3) or (20,15,10,5) and labelled 'fit1', 'fit2' and 'fit3', respectively. 4 For the scalar-glue basis we simultaneously fit N v = (20,3) or N v = (17, 5) ('fit1' and 'fit2') but jointly fitting C SS , C SG and C GG . In all cases, we fit three states (N = 2 in (4.6)), but only the lowest two potentially enter any subsequent analysis. We list the numerical results for the lowest two states ('gr' and 'ex', respectively) in Table II in Appendix B. In all further steps of the analysis we consider all choices of 'fit1', 'fit2' and 'fit3' to propagate any systematic uncertainties. Finally, we also compute the connected correlation function for the pseudoscalar meson, which corresponds to a non-existent state in a N f = 1 theory and in the following is therefore referred to as "fake pion". As we will discuss in the following section, m fake π → 0 can be used as a proxy for the massless limit (see also Ref. [15]). These correlation functions are generated from standard point sources and follow the same functional form as Eq. (4.6) with the replacement Z n X → π (qγ 5 q) † 0 . For these states we perform fits with N = 0 and N = 1. We note that for both κ = 0.1410 ensembles we expect large finite size effects as m fake π L < 3 and therefore discard them from the subsequent analysis.

V. ANALYSIS OF THE SPECTRUM
The goal of this section it so extrapolate the results for the meson spectrum (m P , m I and m S ) to the chiral and infinite volume limit to provide results for ratios of these masses.

A. Defining the chiral limit
We start by determining what the best proxy for the quark mass is. Figure 7 shows the lowest lying state for the pseudoscalar channel. The left panel displays this as a function of the bare quark mass, the right panel as a function of the fake pion mass. By comparing the two panels, it is evident that the fake pion mass is the more suitable choice to define the massless limit as the bare quark mass suffers from large finite volume effects. Those are due to discretisation effects in the computation of the critical parameter κ c entering the definition of the bare subtracted 5 quark mass (see Ref. [41] for a discussion in the case of QCD). In addition, in Ref. [15] it has been numerically shown, for the one-flavour SU (2) gauge theory, that the definition of the massless point from the vanishing of the fake pion mass is consistent with the rigorous definition from the continuum relation between the topological susceptibility and the quark mass 6 . In the following we therefore choose the fake pion mass to define the massless limit. 5 In other words we are saying that data should be compared at fixed bare subtracted quark mass and that differs from the bare mass by a constant related to κc, which has, at finite lattice spacing, a rather strong dependence on the volume [41]. 6 In Ref. [11] the relation ⟨ν 2 ⟩ = ΣV m is in fact established first for one-flavour QCD and then for the case of several flavours. In the equation ⟨ν 2 ⟩ is the topological susceptibility, Σ the fermion condensate and V the four-dimensional volume. We see from the plot in Appendix A that our data are in good qualitative agreement with that relation.

B. Assignment of states
To understand the behaviour of the spectrum we induced by means of our chosen interpolating fields, we investigate how the hadron masses vary as a function of quark mass and volume. We are predominantly interested in mesonic states dominated by qq contributions 7 . These are expected to display a strong quark mass dependence but at most a mild dependence on the volume, whereas any glueball state should only depend weakly on quark mass and volume. Contrary to these, states that depend mildly on the quark mass but strongly on the volume do not correspond to physical states and might be interpreted to be torelon states [42,43].
In Section V A we noted that the pseudoscalar mass is largely volume independent, but depends smoothly and strongly on the quark mass set by m fake π . We therefore identify this with the desired qq-state. In the case of the scalar and vector channels, the situation is more complicated. When comparing results of simulations at the same κ but on different volumes, there are cases that display significant volume dependence on smaller volumes.
For example, the top panel of Figure 8 shows the spectrum as a function of the inverse spatial volume but at fixed κ = 0.1390. We observe that the three largest volumes yield very consistent ground state masses. Contrary, for the two smallest volumes, we see that a lighter state is present in the spectrum, which displays a strong volume dependence. We note that the first excited state on these two volumes is numerically close to the ground state mass extracted on the larger volumes. This picture is further substantiated by investigating the behaviour of the amplitude for the matrix element as we will illustrate with the example of (Z 20 I ) i : In the bottom panel of Figure  8 we show these values for the three states we are fitting. For the three largest volumes, which are displaying a consistent ground state mass, we find that the ground state matrix element (left three magenta circles) is of similar size or larger than the other matrix elements. In contrast to this, for the smallest two volumes the situation is reversed and we find the matrix element of the lowest lying state (right two magenta circles) to be significantly smaller than that of the first and second excited states. We further note that for these two smallest volumes, the matrix element corresponding to the first excited state (rightmost two red diamonds) shows a qualitatively similar behaviour to that of the ground state for the larger volumes. In other words, for the smallest two volumes, the correlation function couples more strongly to the first excited state than the ground state. This is also the reason for the slow approach to the plateau for example in the case of the L/a = 16 and κ = 0.1390 ensemble (c.f. bottom panel of Fig. 6). The strong volume dependence and qualitatively different behaviour with respect to the matrix element indicate that the lowest lying state for the small volumes is not the qq-state we are interested in. Instead, as indicated by the values of the mass and the amplitudes we identify the first excited state with the qq state. In summary, for the vector channel at fixed κ = 0.1390, the qq state corresponds to the lowest lying state for L/a = 32, 24, 20 and to the first excited state for L/a = 16, 12. Corresponding analyses for the other quark masses yield a similar picture. Figure 9 addresses the scalar channel. The top panel shows the mass dependence at fixed volume L/a = 16. The lowest lying state is mass independent in the range of masses we simulate, but the first excited state displays a strong mass dependence. The bottom panel shows the volume dependence at fixed κ = 0.1390. Again, for small volumes, we find a state whose energy increases as the volume increases (lowest state at L/a = 12, 16), as well as a volume insensitive state (lowest state at L/a = 32, 24, 20 and first excited state at L/a = 16, 12). Furthermore, the latter coincides with the state that displayed the strong mass dependence in the top panel. In analogy with the discussion of the vector meson, we conclude that those correspond to a (mass dependent, volume independent) scalar meson state and a (mass independent, volume dependent) torelon state.
By means of similar investigations of the volume and quark mass dependence, we categorise the two lowest lying states on each ensemble and in each channel into the lowest quark mass dependent state (qq) and the remaining state, which in principle can be a torelon, an excited qq or a glueball state. Figure 10 shows the state that has been identified as the relevant qq state for the vector (top) and scalar (bottom) channels. For the large volumes, good agreement is found for all quark masses, whereas for light quark masses and small volumes finite size ef-fects are sizable. We therefore exclude the L/a = 12 and L/a = 16 from our subsequent analysis.
Summarising the discussion in this Section, the qq states we are interested in are easily identified at large volumes and small quark masses as the lowest lying states in the respective channels. Such determinations have the largest impact in the chiral and infinite volume extrapolations we discuss next. However, especially for small volumes, the identification required a more detailed study of the volume and mass dependence of both the energy levels and the overlap factors describing the correlation functions. Those are important lessons we will take into account for future studies at large values of N C .

C. Extrapolation to zero quark mass
We are interested in the spectrum at vanishing quark mass. Since we have not performed a scale setting analysis we focus on ratios of masses in the chiral limit. As discussed above, we will use the fake pion mass to define the zero quark mass limit. A completely model independent fit function would have to include even and odd powers of the fake pion mass. To give a rigorous definition of the fake pion correlator one would have to consider a partially quenched theory constructed by introducing a quark field with the same mass parameter as the original one and quenching it away by a corresponding ghost field [44]. Such a theory would be invariant (at zero quark mass) under transformations in an extended (graded) chiral symmetry group. Depending on whether the symmetry is realisedà la Wigner-Weyl orà la Nambu-Goldstone, one would obtain different relations between the fake pion mass and the quark mass. In the second case (where the symmetry is broken spontaneously by the vacuum and explicitly by the quark mass) the quark mass would turn out to be proportional to the fake pion mass squared. In this case a fit in terms of only even powers of the fake pion mass would be more appropriate.
Any such Gell-Mann-Oakes-Renner-like [45] relation is valid at low energies or very close to the massless limit and in the same limit the fake pion and the pseudoscalar masses should differ significantly, as the first is expected to vanish while the second not. Since in our data we only see small differences between such masses we cannot claim with confidence to be in the regime where such relations apply. We hence favour the more general fit function including even and odd powers of the fake pion mass for our final results. where M is either a mass (m P , m S , m I ) or ratios thereof.
In line with what we discussed above, we separately consider choices where i takes even and odd values or only even values and in both cases either leaving f 0 as a free parameter or setting it to zero. In addition to varying the fit function, we consider cuts to the data, in particular removing the smallest volumes and/or the lightest and/or heaviest masses. An example fit for the case of the pseudoscalar mass (top) and the scalar mass (bottom) is shown in Figure 11. In both of these cases we take the results obtained by 'fit1', keep f 0 as a free parameter and choose n pow = 2. Due to concerns about the finite volume effects, we exclude the smallest volumes (L/a = 12, 16) and the lightest quark mass (κ = 0.1410).
We repeat all extrapolations for the various choices of the correlation function fits, whether or not f 0 is kept as a free parameter and for different choices of n pow . For the lowest order polynomial we restrict the mass range that enters the fit. The datapoints in Fig. 12 show the results for these variations for the pseudoscalar (top) and the scalar (bottom). Only fits with an acceptable p-value of p > 0.05 are shown. The green band in these plots is derived by taking the 68th percentile of the distribution of the underlying bootstrap samples of all the fits which produced an acceptable p-value. We interpret this number to be a good approximation of systematic effects due to correlator fit choices, variations of the chiral fit ansatz and the data included in such a fit.
Ultimately we are interested in the ratio of masses in the chiral limit. We can obtain this in two ways as we will now illustrate on the example of the ratio of the pseudoscalar to the scalar mass: We can either build the ratio m P /m S at finite m fake π and then extrapolate this to the massless limit (method 1), or we can separately extrapolate the pseudoscalar and the scalar masses and then build their ratio (method 2). One example fit of the former is shown in Fig. 13. We observe that part of the mass dependence cancels in the ratios, resulting in a less steep curve than that observed in the individual fits (cf. Fig 11). The coloured stars in the left panel of Fig. 14 show different variations of the fit ansatz, analogous to for the fits of method 1 (method 2) that produced an acceptable p-value.
In general, we notice that the ratio of separate chiral extrapolations leads to larger variations than the extrapolation of the ratio of masses. This is unsurprising as, ensemble by ensemble, the underlying datapoints are statistically correlated, and therefore statistical fluctuations are reduced for the individual ratios of datapoints. Furthermore the extrapolation of the individual datapoints is more difficult to control since the slope with the fake pion mass is steeper. Our preferred number is therefore the direct extrapolation (green band in Fig. 14) whilst the orange band provides a sanity check.
We are now in a position to compare the results of the fits including even and odd powers of m fake π to those only using even powers. These two choices are compared in the two panels of Fig. 14 for the example of the ratio of pseudoscalar to scalar masses. The green bands of the two panels are in ∼ 2σ agreement, lending confidence in the results. However the errorbands of the direct and indirect methods do not overlap for the fit of the even powers only. This is even more pronounced for the case of the ratio m P /m I (c.f. the bottom panel of Fig 17). This numerical evidence further supports our preference for the more conservative fit ansatz including even and odd powers of m fake π and we therefore quote results from this choice as our final numbers.
In addition to m P , m S we have data for the vector mass m I . An example fit for the extrapolation of the vector mass is shown in Fig. 15 (cf. Fig. 11) whilst different fit variations are shown in Fig. 16 (cf. Fig. 12). Finally, we also construct the ratios m P /m I (see Fig. 17) and m I /m S (see Fig. 18) in the chiral limit via the two methods described above.

VI. DISCUSSION AND OUTLOOK
We have presented a detailed study of the spectrum of one-flavour QCD using Wilson fermions with tree-level O(a) improvement.
Results are obtained at one single lattice spacing (approximatively 0.06 fm) for different volumes (up to 32 3 × 64) and several quark masses. After extrapolating to the massless limit we obtain for the pseudoscalar to scalar meson mass ratio and for the pseudoscalar to vector ratio. In Ref.
[6] a prediction using an effective field theory approach and a 1/N C expansion was derived. In the massless limit this reads  corrections and the parameter β are small. Obviously this finding needs to be corroborated by extending our studies to larger values of N C .
Those are further motivated, firstly, by the observation about the slope in the mass for the extrapolation of m P /m S having the opposite sign compared to the prediction in Ref. [6]. Corrections to that start at O(1/N C ) and can therefore be quite large. Secondly, results at larger values of N C will allow assessing the range of validity of the two different theoretical predictions in Refs. [6] and [7]. The latter predicts a value of 1/3 for the ratio m P /m S at N C = 3 in the massless limit.
We have provided an improved estimate, concerning cutoff effects and assessment of systematic errors, compared to previous results that appeared as Proceedings in Ref. [46] (based on Ref. [14]), where a value of 0.410 (41) was found for the pseudoscalar to scalar mass ratio.
Besides having tested the predictions made in Refs. [6,7] for the spin-zero one flavour QCD mesonic state, we further provided information on the vector spectrum that can be interpreted as the leading order prediction for the N = 1 super Yang-Mills vector states.
In order to assess the size of higher order effects we are extending the computation considering N C = 4, 5 and 6. A preliminary account appeared in Ref. [27]. Table III shows values obtained for t 0 on each ensemble. We use two different action densities to compare systematic effects of setting the scale, i.e. the Wilson plaquette action (t Wilson 0 ) and the Yang-Mills action (t YM 0 ). As mentioned before, we quote these values only as reference for other lattice simulations. The quoted uncertainties originate from the numerical integration scheme and are statistical only. As an example we show the dependence of t YM 0 on the fake pion mass in Fig. 20. We note that t YM 0 displays large finite size effects for the smallest ensembles.    [30]. The quoted uncertainties are statistical only and discrete by construction, as the Wilson flow is numerically integrated and saved every step of length 0.05. We quote two estimates using the Wilson plaquette action (t Wilson 0 ) and the Yang-Mills action (t YM 0 ). The difference can be seen as an estimate for systematic uncertainties of the overall procedure.