Nonperturbative investigations of SU(3) gauge theory with eight dynamical flavors

We present our lattice studies of SU(3) gauge theory with $N_f$ = 8 degenerate fermions in the fundamental representation. Using nHYP-smeared staggered fermions we study finite-temperature transitions on lattice volumes as large as $L^3 \times N_t = 48^3 \times 24$, and the zero-temperature composite spectrum on lattice volumes up to $64^3 \times 128$. The spectrum indirectly indicates spontaneous chiral symmetry breaking, but finite-temperature transitions with fixed $N_t \leq 24$ enter a strongly coupled lattice phase as the fermion mass decreases, which prevents a direct confirmation of spontaneous chiral symmetry breaking in the chiral limit. In addition to the connected spectrum we focus on the lightest flavor-singlet scalar particle. We find it to be degenerate with the pseudo-Goldstone states down to the lightest masses reached so far by non-perturbative lattice calculations. Using the same lattice approach, we study the behavior of the composite spectrum when the number of light fermions is changed from eight to four. A heavy flavor-singlet scalar in the 4-flavor theory affirms the contrast between QCD-like dynamics and the low-energy behavior of the 8-flavor theory.


I. INTRODUCTION
The discovery of a Higgs particle at the Large Hadron Collider (LHC) [1,2] was a major step towards the longstanding goal of determining the mechanism of electroweak symmetry breaking. The properties of this particle, which are so far consistent with the Standard Model [3], could also result from new strong dynamics at or above the TeV scale. Lattice gauge theory is an indispensable tool to study the relevant strongly coupled systems, which will generally differ qualitatively from QCD in order to be phenomenologically viable.
The flavor-singlet scalar meson (labelled interchangeably by σ or 0 ++ in the following) is of particular interest as a potential composite Higgs candidate. The strength of current LHC experimental bounds on new particles favors a "little hierarchy" [20] between the Higgs boson and the other states arising from this new sector, which can arise dynamically when the 0 ++ state is found to be light compared to the rest of the spectrum. To assess the viability of models of dynamical electroweak symmetry breaking built on gauge theories with a particularly light 0 ++ meson, detailed knowledge of the appropriate lowenergy effective theory in the presence of this light state is required. Lattice calculations provide crucial input into discriminating candidate low-energy effective theories.
In this paper we continue to explore these issues in the context of SU(3) gauge theory with N f = 8 light fermions in the fundamental representation. Previous lattice studies of this system have identified several features quite distinct from QCD, which make it a particularly interesting representative of the broader class of near-conformal gauge theories. These features include slow running of the gauge coupling (a small β function) [21,22], a reduced electroweak S parameter [23], a slowly evolving mass anomalous dimension γ m [24], and changes to the composite spectrum including a light flavor-singlet 0 ++ scalar [8, 11, 15-17, 23, 25, 26] and a heavier flavor-singlet 0 −+ pseudoscalar [27] (these states are referred to as the σ and η mesons in QCD.) Several lattice groups continue to investigate the 8-flavor theory in order to learn more about its low-energy dynamics and relate it to phenomenological model building. ods, including the computation of the running coupling and its discrete β function [21,22,28,29], exploration of the phase diagram through calculations at finite temperature [30][31][32][33][34][35][36], analysis of hadron masses and decay constants [8, 11, 15-17, 23, 25-27, 37-40], study of the eigenmodes of the Dirac operator [24,[37][38][39]41], and more [42][43][44][45][46][47][48][49][50]. While most of these studies obtain results consistent with spontaneous chiral symmetry breaking in the massless limit for N f = 8 [8, 15-17, 21-26, 28-37, 40], this has not yet been established definitively and some recent works favor the existence of a conformal infrared fixed point (IRFP) [46,47,49]. For example, although all lattice studies of the 8-flavor discrete β function obtain monotonic results, with no non-trivial IR fixed point where β(g 2 ) = 0, it remains possible that an IRFP could exist at some stronger coupling that these works were not able to access.
This possibility can be tested through complementary studies of phase transitions at finite temperature T = 1/(aN t ), where 'a' is the lattice spacing and N t is the temporal extent of the lattice. In a chirally broken system such as QCD, the bare (pseudo-)critical couplings g cr (N t ) of these transitions must move to the asymptotically free UV fixed point g cr → 0 as the UV cutoff a −1 → ∞ and N t → ∞ holding T (N t ) = T cr fixed. In an IR-conformal system, in contrast, the finite-temperature transitions must accumulate at a finite bare coupling as N t → ∞, so that T cr is independent of N t , and remain separated from the weak-coupling conformal phase by a bulk transition.
Unlike running coupling studies, finite-temperature lattice calculations use non-zero bare fermion mass am to give mass aM to the pseudo-Nambu-Goldstone bosons (PNGBs) which appear in the chirally broken phase. If the Compton wavelength of the PNGBs ∼ 1/(aM ) is not small relative to the spatial extent of the lattice, significant finite-volume effects will occur. Results must be extrapolated to the am → 0 chiral limit to ensure that the chiral symmetry breaking is truly spontaneous. Although previous works observed QCD-like scaling of g cr for N f = 8 with sufficiently large masses am 0.01 [30,31,33,34], this did not persist at smaller am ≤ 0.005, where the finite-temperature transitions merged with a bulk transition into a lattice phase.
In Section II we revisit this finite-temperature analysis, employing larger N t than those previous works. Although these larger lattices allow us to consider smaller masses down to am = 0.0025, we find that the finitetemperature transitions still run into a strongly coupled lattice phase before reaching the chiral limit. That is, we are not able to directly confirm spontaneous chiral symmetry breaking. The details of the lattice phase depend on our lattice action, which we also review in Section II. We use improved nHYP-smeared staggered fermions, which conveniently represent N f = 8 continuum flavors as two (unrooted) lattice fields. Staggered fermions (with or without various forms of improvement) are also used by almost all of the other studies summarized above, with the exceptions of Refs. [45,46] (Wilson fermions) and Refs. [23,49] (domain wall fermions).
Beginning in Section III we turn to the main topic of this paper, large-scale studies of the 8-flavor composite spectrum at zero temperature. Section III describes the zero-temperature ensembles that we have generated for this work, which reach the lightest masses considered to date. We use the Wilson flow to set the scale of these ensembles, and observe the Wilson flow scale to be more sensitive to the fermion mass than would be expected for lattice QCD. In Section IV we review the details of our staggered spectrum analyses, separately considering the flavor-singlet scalar σ that involves contributions from fermion-line-disconnected diagrams.
Our results for the 8-flavor spectrum are discussed in Section V. Their most significant feature is the presence of a remarkably light flavor-singlet scalar particle (σ), which remains degenerate with the PNGBs (π) down to the lightest fermion masses reached so far by lattice calculations. This is a significant contrast with QCD, which we strengthen by carrying out a similar 4-flavor spectrum calculation. Using the same lattice action and analysis procedure, our QCD-like 4-flavor results do not produce a light scalar, as expected. The comparison is reported in Fig. 1 where the masses of composite states are normalized by the π decay constant (F π ).
However, other aspects of our 8-flavor spectrum results are qualitatively similar to QCD. At the most basic level, the ratio of the vector mass to the pseudoscalar mass steadily increases as we approach the chiral limit, providing indirect indication that the theory exhibits spontaneous chiral symmetry breaking. We also find the ratio of the vector mass to the pseudoscalar decay constant to be comparable to its QCD value, M ρ /F π ≈ 8, and rather constant as we decrease the masses. In the context of models of dynamical electroweak symmetry breaking, this suggests (via the Kawarabayashi-Suzuki-Riazuddin-Fayyazuddin (KSRF) relations [51,52]) a multi-TeV scale vector resonance with a large decay width Γ/M 0.2, comparable to that of the QCD ρ meson. This is broader than the typical width assumed in past LHC searches for such states [20]; dedicated searches for broad resonances, although challenging, are well-motivated by the lattice results.
We summarize our conclusions and prospects for further progress in Section VI. In particular, we focus on the issue of the appropriate low-energy effective field theory (EFT) to describe the 8-flavor spectrum we observe. A consequence of the light flavor-singlet scalar is that we cannot expect to carry out chiral extrapolations by fitting our data to chiral perturbation theory (χPT), which assumes that the PNGBs are much lighter than all other particles. Finally, in the appendices we provide additional information about auto-correlations, topological charge evolution, more technical details about fitting correlation functions for the flavor-singlet scalar, and studies of finite-volume and discretization effects.

II. LATTICE ACTION AND FINITE-TEMPERATURE PHASE DIAGRAM
Our numerical calculations use improved nHYPsmeared staggered fermions [53,54] with smearing parameters α = (0.5, 0.5, 0.4), and a gauge action that includes both fundamental and adjoint plaquette terms with couplings β F and β A , respectively, related by β A /β F = −0.25 [38]. This lattice action was used in several previous studies of the 8-flavor system, including explorations of the phase diagram [33,34,38], the composite spectrum [26], the discrete β function [21] and the scale-dependent mass anomalous dimension γ m (µ) [24]. Using the same lattice action for all of these complementary investigations makes it easier to compare their results and thereby gain more comprehensive insight into the dynamics of N f = 8.
The first work using this action observed a strongly coupled " S 4 " lattice phase in which the single-site shift symmetry (S 4 ) of the staggered action is spontaneously broken [38]. In the massless limit, a first-order bulk (zero-temperature) transition around β F ≈ 4.6 separates the S 4 phase from the weak-coupling phase where the continuum limit is defined. At even stronger couplings there is a second bulk transition into a chirally broken lattice phase. A similar phase structure has been seen by other many-flavor lattice investigations using different improved staggered actions [55,56]. 1 However, the characteristics of these strong-coupling phases are not universal and depend on the details of the lattice action. Although in this section we scan the lattice phase diagram, including the transition into the S 4 phase, our zero-temperature calculations reported in the rest of the paper will consider a coupling β F = 4.8 safely on the weak-coupling side of this bulk transition.
The presence of the S 4 phase prevents lattice investigations from reaching arbitrarily strong couplings. For example, Ref. [21] was only able to determine the continuum-extrapolated discrete β function for renormalized couplings up to g 2 c 14 (in finite-volume Wilson flow renormalization schemes introduced by Ref. [59]). As summarized in Section I, although this β function is monotonic throughout the accessible range of couplings, this does not guarantee that the 8-flavor theory exhibits spontaneous chiral symmetry breaking. It remains possible that, at stronger couplings, the β function might reach an extremum and then return to β(g 2 ) = 0 at some large g 2 15. (Indeed, this happens in four-loop perturbation theory in the MS scheme, which predicts g 2 ≈ 19.5 [60,61], but perturbation theory seems unlikely to be reliable at such strong couplings.) In the remainder of this section we present a complementary search for spontaneous chiral symmetry breaking in the 8-flavor system, by studying its finitetemperature phase diagram. Initial results from this work appeared in Ref. [36]. As described in Section I, in order to establish spontaneous chiral symmetry breaking, the finite-temperature transitions we observe at non-zero bare fermion mass am must persist in the chiral limit. Previous finite-temperature studies with am ≥ 0.005 and N t ≤ 20 instead found that these transitions merged with the bulk transition into the S 4 phase [33,34]. Here we move to larger 40 3 × 20 and 48 3 × 24 lattice volumes, which allows the exploration of smaller masses, 0.0025 ≤ am ≤ 0.01.
For each combination of lattice volume and mass am, we generate ensembles of gauge configurations with β F ranging from the strong-coupling S 4 phase to the deconfined phase at weak coupling. To generate these ensembles we use either QHMC/FUEL [62] or a modified version of the MILC code, 2 in both cases employing the hybrid Monte Carlo (HMC) algorithm with a second-order Omelyan integrator [63] accelerated by additional heavy pseudofermion fields [64] and multiple time scales [65]. (The only exception is the N f = 4 zerotemperature ensemble with am f = 0.003, which used a approximate force-gradient integrator with θ = 0.109 [66,67].) We monitor several observables to identify both bulk and finite-temperature transitions, including ψψ , the Polyakov loop, S 4 order parameters introduced in Ref. [38], and the massless Dirac eigenvalue spectrum ρ(λ). Of these observables, the most useful are the Polyakov loop and ρ(λ), for which representative results are shown in Figs. 2 and 3, respectively. To improve the Polyakov loop signal we compute it after applying the Wilson flow, a continuous transformation that smooths lattice gauge fields to systematically remove short-distance lattice cutoff effects [68]. For sufficiently large flow time t the Wilson-flowed Polyakov loop P L W shows a clear contrast between confined systems with vanishing or small P L W 1 and deconfined systems with large P L W ∼ 1, even for large N t ≥ 20 which tend to produce noisy results for the unimproved Polyakov loop. This observable is a modern variant of the RG-blocked Polyakov loop investigated in previous studies, which showed that the improvement in the signal does not affect the location of the transition [33,34].
In Fig. 2 we consider the real part of the Wilson-flowed Polyakov loop computed at flow time t = 4.5 for 40 3 ×20 lattices with am = 0.01, 0.005 and 0.0025 vs. the bare lattice coupling β F . This flow time was chosen to produce √ 8t/N t = 0.3, which we keep fixed by considering t = 6.48 for N t = 24. As the fermion mass decreases the finite-temperature transitions in Fig. 2 steadily sharpen and move to stronger coupling. At am = 0.0025 the finite-temperature transition merges with the bulk transition into the S 4 phase, implying that even larger volumes are required to establish whether chiral symmetry breaking occurs spontaneously in the massless limit. Our N t = 24 results behave similarly, as we discuss further below.
First, we review how we determine the bulk transition into the S 4 phase. Ref. [38] identified two order parameters of the single-site staggered shift symmetry, which take the form of differences between local observables on neighboring lattice sites, ∆P µ = ReTr n,µ − ReTr n+ µ,µ nµ even (1) Here ReTr n,µ is the average real trace of the six plaquettes that include the gauge link U µ (n) connecting sites n and n + µ. χ(n) is the staggered fermion field, while α µ (n) = (−1) ν<µ nν is the usual staggered phase factor. Finally, the expectation value · · · nµ even is taken only over sites whose µ component is even.
In addition, Ref. [38] found that the eigenvalue spectrum of the massless Dirac operator exhibits an unusual 'soft edge' [69][70][71] in the S 4 phase, ρ(λ) ∝ √ λ − λ 0 with λ 0 > 0. This allows the spectral density to distinguish between all three phases of interest, as illustrated in Fig. 3 for 40 3 ×20 ensembles with am = 0.005. At weak couplings, including β F = 4.8 in the figure, the system is deconfined and chirally symmetric, with ρ(0) = 0 and a gap below the smallest eigenvalue. The gap becomes larger for weaker couplings omitted from the plot. At intermediate couplings such as β F = 4.7 we observe the expected chiral symmetry breaking, with ρ(0) = 0 and a small slope dρ dλ . Finally, once we enter the S 4 phase at stronger couplings β F 4.6 the spectral density shows the onset of a soft edge with ρ(λ) ∼ √ λ. This is most clear at β F = 4.5, where a gap has reopened and ρ(0) = 0 indicates the restoration of chiral symmetry.
The four curves shown in Fig. 3 therefore indicate that the N t = 20 chiral transition with am = 0.005 is between 4.7 < β (c) F < 4.8, consistent with the Wilsonflowed Polyakov loop shown in Fig. 2. In addition, we can see that the transition into the S 4 phase is between 4.6 < β (c) F < 4.7. While the β F = 4.6 spectral density does not show a clear gap and may possess a non-zero ρ(0) in the infinite-volume limit, its dependence on λ is consistent with the square-root behavior of the soft edge, in contrast to the curves with β F ≥ 4.7. Since the S 4 order parameters discussed above indicate that our β F = 4.6 ensemble is in the S 4 phase, it seems most likely that the lack of a clear gap is related to the nonzero sea mass am = 0.005. Based on the observables described above, we identify the N t = 20 and 24 finite-temperature transitions shown in Fig. 4. For N t = 20 and am = 0.0025 we find the system at β F = 4.65 to be in the weak-coupling phase while that at β F = 4.6 is in the S 4 phase. Although a finer scan of intermediate values 4.6 < β F < 4.65 might still reveal a narrow chirally broken phase, we conclude that the N t = 20 finite-temperature transitions have effectively merged with the bulk transition by am ≈ 0.0025. For the larger N t = 24 we see that the finite-temperature transitions at a fixed mass move to weaker couplings. While this is the expected behavior for a chirally broken system, it is not sufficient to conclusively establish that the 8-flavor theory exhibits spontaneous chiral symmetry breaking in the massless limit. Even though the Finite-temperature transitions from lattices with temporal extents Nt = 20 and 24, with lines connecting points to guide the eye. The region above these lines is confined and chirally broken, while the region below is deconfined and chirally symmetric. The left edge of the plot indicates the bulk transition into the S 4 lattice phase. The finite-temperature transitions merge with this bulk transition at am > 0, preventing a direct confirmation of spontaneous chiral symmetry breaking.
FIG. 5. Considering fixed βF = 4.7 and 4.8, the Nt dependence of the finite-temperature transitions in Fig. 4 suggests that lattice volumes with temporal extent roughly Nt 34 are needed in order to establish spontaneous chiral symmetry breaking in the massless limit at these values of βF a. N t = 24 transitions manage to reach slightly smaller am before running into the S 4 phase, it is clear that these transitions will also merge with the bulk transition at a non-zero mass.
Although our new 40 3 × 20 and 48 3 × 24 finitetemperature investigations do not suffice to establish spontaneous chiral symmetry breaking, it is still significant that the finite-temperature transitions move to smaller masses as the temporal extent of the lattice increases. If we assume that this trend continues for larger N t , then our results allow us to estimate the lattice volumes that would be needed in order for the transitions to reach the chiral limit without running into the S 4 phase.
From Fig. 4 we can read off the masses at which the N t = 20 and 24 transitions cross a fixed value of β F = 4.7 only slightly weaker than the bulk transition into the S 4 phase. These two points, plotted vs. N t in Fig. 5, suggest that the finite-temperature transitions for N t 34 should reach the am = 0 chiral limit rather than running into the S 4 phase. This figure also considers the somewhat weaker coupling β F = 4.8 used in our main zero-temperature investigations. For this β F the transitions are also projected to reach the chiral limit around roughly the same N t 34.
In principle it might be possible to construct an alternate lattice action that would allow us to reach stronger couplings before encountering a lattice phase. Then we would be able to obtain comparable results from smaller lattices. However, an attempt to do this by adding a second nHYP smearing step was not successful [21]. We do not currently plan to generate the larger-volume 64 3 ×32 and 72 3 × 36 lattice ensembles that would be needed to pursue a more definitive demonstration of spontaneous chiral symmetry breaking using our current action. While such lattice volumes are within the reach of existing algorithmic and computing technology (cf. the zero-temperature 64 3 × 128 investigations presented below), they would consume significant resources that we prefer to invest in more promising directions discussed in Section VI.

III. ZERO-TEMPERATURE LATTICE ENSEMBLES AND SCALE SETTING
We now focus on our main N f = 8 investigations, which determine the composite spectrum of the theory at zero temperature. On the basis of the phase diagram discussed above, we carry out our computations at a relatively strong coupling β F = 4.8 that is still safely on the weak-coupling side of the S 4 phase. This relatively strong coupling, in combination with large lattice volumes up to 64 3 ×128, allows us to consider the lightest masses yet to be reached by lattice studies of the 8-flavor theory. Table I summarizes the ensembles of gauge configurations that we have generated using the same software and algorithm as described in the previous section. These include four 4-flavor ensembles, two of which (with lattice volume 24 3 ×48) are matched in alternate ways to the 8flavor 24 3 ×48 ensemble with the largest am f = 0.00889: one (with β F = 6.6) matches the π and ρ masses in lattice units, while the other (with β F = 6.4) matches the π mass and √ 8t 0 in lattice units.
In addition to recording the lattice volume, fermion mass and available statistics, Table I also reports values for the reference scale √ 8t 0 introduced in Ref. [72]. This reference scale is defined through the Wilson flow (discussed in the previous section), by requiring that FIG. 6. Four determinations of t 2 E(t) (using clover-and plaquette-based definitions of the energy density, with and without perturbative improvement) illustrate discretization artifacts for the N f = 8 32 3 ×64 ensemble with am f = 0.005. These artifacts are most significant for small t 1, and are ameliorated by the corresponding tree-level finite-volume perturbative corrections, which reduce the larger plaquette-based results while increasing the smaller clover-based results. For future EFT analyses we want √ 8t0 < 1/Mπ, corresponding to t smaller than the vertical red line at t = 1/(8M 2 π ).
is evaluated after flow time t using the standard cloverleaf construction of F µν . The factor of δ(L, t) is the treelevel finite-volume perturbative correction introduced by Ref. [73], which reduces discretization artifacts in the energy density. The specific form of this correction depends on the gauge action, on the (Wilson) action used in the gradient flow transformation, and on the (clover) operator used to define E(t). At tree level our fundamentaladjoint plaquette gauge action is equivalent to the Wilson gauge action, so in the terminology of Ref. [73] we use the WWC scheme. In Fig. 6 we confirm that this perturbative correction does indeed reduce discretization artifacts, by comparing our clover-based results against the Wilson plaquette operator is the trace of the plaquette normalized to 3. We consider the N f = 8 32 3 ×64 ensemble with m = 0.005; the other ensembles show similar behavior. The two lattice definitions of the energy density differ by discretization artifacts, which are most significant for small t 1 where E plaq can be more than twice as large as the clover-based result. Appropriately, the perturbative correction reduces the plaquettebased results while increasing the clover-based results, reducing the overall discretization artifacts. To simplify comparisons with other groups' results we do not include the t-shift improvement of Ref. [74] in our analyses.
In addition to quantifying discretization artifacts, Fig. 6 also indicates that these artifacts are more se-   Perturbatively improved clover-based t 2 E(t) (Eq. (3)) for all 8-flavor ensembles in Table I,  vere for E plaq than for the clover discretization of the energy density. For sufficiently large t the discretization artifacts are removed by the Wilson flow and t 2 E(t) becomes approximately linear. The perturbatively improved clover-based curve is roughly linear already for t 1, where significant non-linearities are still visible in the plaquette-based results. This motivates our choice of the perturbatively improved clover discretization in Eq. (3). Figure 7 shows the corresponding results for all 8-flavor ensembles in Table I. We note that the con-dition t 2 E(t) t=t0 = 0.3 produces √ 8t 0 < 1/M π for all these ensembles, with the separation growing as the fermion mass decreases. Investigation of the behavior of the scale √ 8t 0 under different EFT descriptions of this theory, as has been carried out for chiral perturbation theory [75], may shed light on this relation and allow discrimination between different possible EFTs.
Finally, in Fig. 8 we plot our results in Table I for the 4-and 8-flavor Wilson flow scale √ 8t 0 vs. the fermion mass am f . The figure shows that the N f = 8 scale increases rapidly as am f decreases, with significantly more sensitivity to the fermion mass than is seen for N f = 4 or for lattice QCD (cf. , for example, Fig. 4 of Ref. [76]).

IV. SPECTRUM ANALYSIS DETAILS
In this section we discuss our analysis of hadronic two-point correlation functions to extract estimates for hadron masses and decay constants of the 8-flavor theory. In the continuum theory the global symmetries after spontaneous symmetry breaking are SU(8) V ×U(1) B . On the lattice, the staggered discretization further reduces the flavor symmetry to U(2) × SW 4 ⊂ SU(8) V × U(1) B [77]. The continuum theory has 63 massless Nambu-Goldstone bosons (NGBs) in the chiral limit whereas the staggered theory has 4 massless NGBs, corresponding to the unbroken U(2) subgroup, plus 59 light PNGBs which become massless in the continuum limit. We use a typical QCD naming scheme for the flavor- nonsinglet mesons, including the NGB pseudoscalar (π), the vector (ρ) and the axial-vector (a 1 ) mesons, to indicate the J P C quantum numbers of the mesons. In the continuum limit, these mesons will transform as the adjoint representation of SU(8) V but at finite lattice spacing they will transform under representations of the U(2) × SW 4 staggered flavor group. For three states we compute both the mass (M π5 , M ρs , M a1,5 ) and the decay constants (F π , F ρ , F a1 ). We use the decay constant normalization corresponding to the QCD value F π = 92.2(1) MeV (cf. section 5.1.1 of Ref. [78]). We also compute the masses of other flavor-nonsinglet mesons, several baryons including the lightest nucleon (M N ), and the flavor-singlet scalar meson (M σ ). A complete list of computed masses is given in Table II. The flavor-singlet σ state is the most challenging state to study because it mixes with the vacuum and receives contributions from fermion-line-disconnected diagrams.
The use of staggered fermions requires some discussion of the staggered flavor (or taste) quantum number. A complete analysis of the staggered meson correlators constructed from fermion bilinear interpolating operators re-quires considering operators distributed, at a minimum, over the 2 4 unit hypercube. In this paper we only consider the operators listed in Table II, all of which are at most one-link separated. Below we will discuss the spin and taste structures of our interpolating operators in more detail.
We split the discussion of our spectrum analyses into three parts. In the next subsection we consider the computation and analysis of the connected spectrum, including decay constants. We then separately discuss the additional disconnected observables needed for the computation of the flavor-singlet scalar σ correlator, and the corresponding analysis. Finally, we explain our procedure for numerical analysis of the correlators including determination of systematic errors associated with the choice of fit range [t min , t max ], which is common to all of the spectroscopy. Although we have also analyzed gluonic operators that couple to the scalar channel, these produce much noisier data from which we could not extract meaningful results, so we will not discuss them further.

A. Connected spectrum
To extract the mass of the π (both PNGB and tastesplit), ρ, a 1 and N , we use Coulomb gauge-fixed wall sources [79]. To extract the decay constants F π , F ρ and F a1 , we additionally use the conserved current interpolating operators A 4 , V i and A i , respectively, where i indexes spatial directions. Special to the PNGB π 5 , we additionally use the point pseudoscalar operator, P , and a random Gaussian wall source thereof. In all cases we use point operators at the sink. We compute all twopoint functions using N t /8 sources distributed evenly in Euclidean time, averaging over those data to obtain a single estimate in each channel from each analyzed configuration.
As summarized in Table I, we bin all estimates of observables to reduce the effects of auto-correlations, setting the bin size to be of the same order as the autocorrelation time of the chiral condensate ψψ , which we determine with UWerr [80]. In Appendix A we provide more information about auto-correlations.
In general, a meson two-point function where both the source and sink interpolating operators have the same quantum numbers is parameterized as List of the states we analyze on our 8-flavor lattice ensembles. Each row indicates a type of hadronic correlator and the columns indicate the spin and flavor (taste) quantum numbers of the states in the direct (non-oscillating) and oscillating channels, cf. Eq. (4). The source labels are "W " for wall sources, "C" for conserved currents, and "P " for the point pseudoscalar operator; the 0 ++ analysis uses a combination of regular and stochastic sources as described in the text. In the last column we also identify the states we used to extract the decay constants.
where Γ( δ) refers to the spin-taste structure of the interpolating operator in channel X, and τ a is an appropriate generator of the flavor symmetry (since we are computing flavor-nonsinglet states). In the staggered formulation a meson operator with quantum numbers Γ will also couple to a "parity partner" state with quantum numbers Γγ 4 γ 5 ξ 4 ξ 5 [81]. A single interpolating operator will also couple to multiple states with the same quantum numbers. We parameterize the direct-state energies and amplitudes as E X,i and A X,i , respectively, with oscillatingstate E X,i and A X,i . In both cases 0 ≤ i ≤ N exc for N exc excited states in addition to the i = 0 ground states; our central analysis sets N exc = 1 for all states. With this in mind, E X,0 refers to the ground-state energy of the parity partner state in a given channel. This contribution is non-zero except for the Goldstone channel, where the would-be parity partner has the quantum numbers of the staggered vector charge. The sign of the term e −E X,i (Nt−t) is positive for all correlators except A 4 (0)P (t) , where it is negative due to odd timereversal symmetry. Since meson correlators are symmetric or anti-symmetric about the middle timeslice N t /2, the data are 'folded' by averaging Staggered baryons have a more complicated structure, and similarly more complicated interpretation of parity partners. Ref. [82] discusses the parameterization of baryon two-point functions. For baryon correlators we do not include the backwards-propagating e E X,i (Nt−t) states and fit to C(t) only for t ≤ 3N t /8. The signal-to-noise in our correlators is poorest for t > 3N t /8, so we do not expect this restriction to significantly reduce the precision of our results.
For the pseudoscalar, vector, and axial-vector channels, we carry out joint fits to two correlators, one containing the appropriate conserved-current source and the other using a wall source with the same quantum numbers. This joint fit improves the precision of the ground-state energy and thus our determinations of the mass and decay constant in each channel. The particular operators used for the determination of the decay constants are noted in Table II.
We define the continuum decay constants via for the pseudoscalar, vector, and axial-vector channels, respectively. Here i , with i = 1, 2, 3, are polarization vectors. Because we use conserved staggered currents no renormalization factors appear (i. e. the Z-factors are equal to 1 exactly.) We can also define F π from the Goldstone pseudoscalar interpolating operator [77] where m f is the fermion mass. Since staggered fermions preserve an exact chiral symmetry this fermion mass receives no additive renormalization.
Since each staggered species corresponds to multiple continuum Dirac fermions, there is a non-trivial conversion factor between the continuum matrix elements and the lattice matrix elements [77]: where 1/ √ 4 comes from the four continuum flavors per staggered fermion, s refers to one of the meson states π, ρ, a 1 , and O is the corresponding interpolating operator.

B. Disconnected spectrum
Unlike the connected spectrum, estimating the mass of the flavor-singlet 0 ++ (also denoted σ here) requires computing fermion-line-disconnected diagrams. We only use one interpolating operator for the 0 ++ meson, the zero-momentum local operator which has a spin-taste-singlet, as well as flavor-singlet, structure. Because the operator is a flavor singlet, there are two ways to contract the staggered fermion fields, which differ by a minus sign because of anticommutation relations, corresponding to two different estimates of the staggered Dirac propagator G F : The two types of contractions split the computation of the 0 ++ correlator into two pieces: a disconnected (double-trace) piece D(t) and a connected (single-trace) piece −C(t). The factor of N f /4 is a relative loop counting factor between the two pieces, where the factor of four corresponds to the four continuum Dirac flavors encoded by each staggered lattice field.
Exactly computing the disconnected piece requires evaluating all diagonal elements of the inverse of the Dirac matrix. This is computationally impractical. Instead of explicitly computing the diagonal of the inverse, we use an improved stochastic trace estimator with dilution, which we now review step by step.
First, we consider a stochastic estimate for any element of the inverse of a complex matrix S(x, y), where x and y index rows and columns, respectively. We can construct a set of noise vectors, η i (x), satisfying where I(x, y) is the identity matrix. Some common choices for η i (x) are U(1), Z 2 , Z 4 , or an appropriately scaled normal distribution. In this paper we always use either U(1) or Z 4 noise. We use N i = 2 or 3 noise sources for L = 64 and N i = 6 for all L < 64. For each η i (x) we define the vector A stochastic estimate of the inverse matrix S −1 (x, y) is given by For finite N i the error of this approximation decreases ∝ 1/ √ N i as expected.
In our context S(x, y) is the staggered Dirac matrix. The inverse of the Dirac matrix is the Green's function G F ( x, t; y, t ), where we suppress color indices. The trace that enters the disconnected (double-trace) piece of the O 0 ++ correlator can be estimated by When forming the double-trace correlator, we need to avoid taking the product of two estimates from the same source. More explicitly, the double-trace correlator must be computed as where we impose j = i to avoid quadratic noise contributions. This is modified, and to some extent relaxed, in the context of dilution. Away from the infinite-source limit, stochastic sources can suffer from noisy couplings to other states. In particular, stochastic estimates of single-trace correlation functions can be sensitive to states with lighter masses. This is solved by dilution, where instead of computing the propagator from a full-volume noise source, we first partition the source into separate pieces and then compute propagators for each individual piece [83]. The ultimate effect of this partitioning is that some off-diagonal components of the sum in Eq. (12) are exactly zero even with finite N i . In practice, appropriate use of dilution improves convergence: instead of the O(1/ √ N i ) convergence of stochastic methods, dilution can give an effective O(1/N i ) convergence for an equal amount of computational work.
We take advantage of dilution in time, color, as well as even/odd in space to individually compute the singletrace pieces of the disconnected correlator. We can then construct the disconnected correlator, D(t), by computing the all-to-all correlator of these single-trace operators. Dilution in time allows us to relax the constraint j = i in Eq. (16) for t = 0: the contributions from t = 0 come from different noise sources because the diluted η sources are non-zero on different timeslices. The issue at t = 0 can be avoided by imposing j = i uniquely for that separation. The construction of the disconnected correlator can be efficiently done in Fourier space via the convolution theorem. Because the 0 ++ channel has the same quantum numbers as the vacuum, the disconnected piece suffers from a vacuum contamination that must be removed, as we discuss below.
Instead of using Eq. (15) for O 0 ++ (t) we employ the improved stochastic estimator which follows from the Ward identity for the staggered chiral symmetry [77]. With this definition, the stochastic estimate for O 0 ++ (t) can be written as This is an improved estimator because it involves only positive-definite contributions (corresponding to the positive Goldstone propagator in Eq. (17)).
The connected piece −C(t) is much simpler to compute. While point sources could be used for this computation, in this work we reuse the stochastic propagators discussed above, simply contracting them in a different way as described by Ref. [83]. Whereas the full correlator S(t) couples to the 0 ++ meson as its lightest singleparticle state, the positive-definite correlator −C(t) by itself instead couples to the a 0 meson.
Neglecting excited states and scattering states, we parameterize S(t) and −C(t) by The staggered oscillating partner of the 0 ++ meson is a pseudoscalar which is a singlet under the exact U(2) flavor subgroup but a non-singlet under the SW 4 staggered flavor group, so we label it η 45 (cf. Tab. II).
On the other hand, D(t) is not a well-defined correlator in the sense that it does not uniquely couple to one set of quantum numbers. Numerically, we can still pa-rameterize it as a linear combination of S(t) and −C(t), For M 0 ++ < M a0 , we can extract the mass of the 0 ++ from D(t) alone for asymptotically large t [8]. This correlator appears to have smaller excited-state contamination than S(t), leading to longer plateaus in effectivemass plots. On the other hand, consistent with the construction of D(t) as a difference between correlators, the excited states are not positive definite, and cause the correlator to curve down instead of up at small t, as shown in Fig. 9.
In practice, instead of fitting S(t) or D(t) directly to the functional forms shown above, we fit to the finite difference S(t + 1) − S(t). This has the advantage that the vacuum contribution a vac is removed entirely before we attempt any nonlinear fits. Even with the finite difference, it is difficult to extract the ground state from the individual correlators S(t) and −C(t): both channels suffer from strong excited-state contamination, especially because, in the infinite stochastic source limit, we are effectively using the point, or local, source and sink operators. Instead, we carry out a joint fit to S(t) and D(t) simultaneously. The opposite sign of the excitedstate contamination in D(t) helps to counterbalance the significant positive contamination in S(t), leading to a more robust joint estimate of M 0 ++ .
In Appendix B we give further numerical details of our extraction of M 0 ++ , including comparisons between our main joint-fit-based analysis and individual fits to the D(t) and S(t) finite-difference correlators.

C. Fit selection and fit-range systematic error
To determine the masses and decay constants, we carry out fully correlated fits over a fit range [t min , t max ] for all of our correlators, subject to the models described above. Those restrictions allow for more reliable determination of the correlator covariance matrix with finite statistics, and also suppress contributions from higher excited states at very small t. We make use of the lsqfit Python package [84] for our fits. All results presented here include one excited state in each of the oscillating and non-oscillating channels, i.e., we include four states in total to describe any correlator with an oscillating parity-partner contribution. We include priors on the energies of all states for numerical stability only: the prior values are aE = 1(5) for ground states, and log(∆E) = −1(4) for the splitting between each excited state and its respective ground state.
For many lattice studies, it is common practice to fix [t min , t max ] based on some empirical analysis of the cor- relators. We select our central fits based on a version of the empirical 'MILC criterion' [85] Here p is the unconstrained p-value of the fit, N dof is the number of degrees of freedom in the fit, and σ n /µ n is the best-fit uncertainty on the nth fit parameter normalized by its mean value. We include only the ground-state parameters in the sum in the denominator. The choice of a specific fit window in principle introduces a systematic effect into the analysis if the resulting masses or decay constants have some residual dependence on [t min , t max ]. To account for this possibility, we assign an additional fit-range systematic error. Over all fits considered for a particular correlator, we compare the central result (mass or decay constant) as obtained for all fits passing a minimal fit quality cut, in our case p ≥ 0.1, and passing the other cuts described in the next paragraph. We further discard all fits which agree with the central result at 1σ, as well as fits which are consistent with zero at the 3σ level, to avoid unnecessarily inflating our errors due to outlier fits with large uncertainties. The fit-range systematic is then taken to be half of the maximum difference between central values of fits passing all cuts.
To carry out the analysis described above, we consider all possible combinations of [t min , t max ] subject to the following constraints. First, we require t max −t min < 2 √ N s , where N s is the number of independent Monte Carlo samples; this cut is placed to avoid including fits where the data covariance matrix is ill-determined by our statistics. We also require that N dof > 0, which imposes a minimum separation t max − t min . To exclude data with poor signalto-noise which could inflate the overall error, we require t max ≤ 3N t /8, and to avoid a regime where the ground state is poorly determined we require t max ≥ N t /4. Finally, we fix t min ≥ 4, as the data at smaller t suffer from substantial excited-state contamination A representative plot showing the determination of this fit-range systematic error is shown in Fig. 10. For each fixed t min , the point with solid error bars represents the best fit, while the dashed error bars, if any, represent the fit-range systematic determined by varying t max at that particular value of t min . The colorized point shows the best fit as determined from the MILC criterion, and the dashed line on that point represents the total fit-range systematic from consideration of all [t min , t max ].

V. SPECTRUM RESULTS AND DISCUSSION
Our results for the spectrum of light bound states in the 8-flavor and 4-flavor theories are collected in Tables III and IV, and plotted in Figs. 1, and 11, and 12. The errors shown for the results include both statistical and fit-range systematic errors, combined in quadrature. No systematic uncertainty for finite volume or lattice discretization is included. We believe that the former effect is negligible for our results, but discretization effects are difficult to estimate since we work at a single value of the bare coupling β F . Appendix C contains a more detailed study of these systematic effects.
The most striking feature apparent in the N f = 8 spectrum is the relative lightness of the flavor-singlet scalar state, σ. On most ensembles it is degenerate with the pion within our uncertainty; in all cases, in the range of fermion masses we study both σ and π are much lighter than the next heaviest hadron, which is the ρ meson. This result stands in contrast with QCD, where phenomenological estimates place the f 0 flavor-singlet scalar between 400-550 MeV [19], above the threshold for two-pion decay. Lattice QCD calculations with heavierthan-physical pion masses show that the f 0 state remains heavy relative to the π states, eventually becoming heavier than M ρ [86][87][88].
Our study of the 4-flavor theory provides a direct test for systematic effects in our analysis of the σ scalar that might give an artificially small value for its mass. The spectrum results for N f = 4 shown in Fig. 12 indicate that the σ mass in this case is qualitatively heavier, closer to the ρ meson than to the pions. We surmise that the appearance of a light σ scalar meson is not an artifact of our analysis procedure, and the additional fermion species in the N f = 8 theory seem to be important for σ to be light.
Decay constants for the eight-flavor theory are plotted  in Fig. 13. The decay constants F ρ and F a1 suffer from significant fit-range systematic errors in our analysis. We can observe qualitatively that they tend to be larger than F π and roughly degenerate with each other.
Here we do not attempt to extrapolate our results to the chiral limit am f → 0, as this requires the choice of a specific chiral EFT to model the data, and the presence of the light σ state indicates that a careful study of possible EFT descriptions is needed. We defer such a study to future work, and here present some ratios of quantities at finite m f that may provide useful insights.
We first present the ratio of meson mass M ρ /M π in Fig. 14. This ratio should diverge in the m f → 0 limit for a theory with spontaneous chiral symmetry breaking, as the pion mass vanishes while the ρ mass remains finite. On the other hand, in either an infrared-conformal theory whose masses exhibit hyperscaling [89] or in a theory with very heavy fermion masses which exceed the binding energy, this ratio should remain roughly constant. Our results in the fermion mass range studied show growth of this ratio with decreasing m f , but no clear indication of divergence. This effect has been argued [47] to be due to finite lattice volume, which we discuss in Appendix C.
The KSRF relations [51,52] are a pair of equations relating the vector and PNGB masses and decay constants, based on current algebra and vector-meson dominance.
In our conventions, they read These relations, and their implications for the strong decay process ρ → ππ, were discussed previously in [15] in more detail. Figure 15 shows our numerical results for the KSRF relations. The top panel shows F ρ /F π , which is again found to be relatively independent of m f and consistent with the KSRF prediction of √ 2. The bottom panel assumes the validity of the second KSRF relation and estimates the coupling g ρππ , which is again relatively mass-independent and qualitatively similar to the value inferred from QCD experiment. This implies a broad decay width Γ/M ∼ 0.2 for a vector resonance in any model of electroweak symmetry breaking built on this theory, given the assumptions discussed in [15].

VI. CONCLUSIONS
The presence of a light, unflavored stable scalar meson σ in the spectrum of the SU(3) N f = 8 theory which is approximately degenerate with the PNGBs π when M π /M ρ < 0.5 is perhaps the most dramatic difference between this theory and QCD where M π /M ρ ≈ 0. Decay constants of the 8-flavor theory, defined in equations (5)- (7). Both the vector and axial-vector decay constants are fairly imprecise, suffering from significant fitrange systematic effects in our analysis. above decay threshold [19]. Our result remains consistent with studies of the same theory reported by the LatKMI collaboration [16] at heavier fermion masses. There is growing evidence from QCD that for heavier PNGB masses M π /M ρ 0.46 the f 0 may become stable with a mass just below the decay threshold [88,90,91], but this is still very different from the SU(3) gauge theory with larger number of flavors or other near-conformal theories with light scalars [7-10, 13, 14].
However, there are other significant differences between N f = 8 and QCD which have been noted in the past [23] and remain present in this work. Of particular interest is the relatively steep slope of the pion decay constant F π vs. am f in the light-fermion-mass regime (cf. Fig. 13). The leading-order prediction from chiral perturbation theory is that F π should be independent of the fermion mass, so that significant am f dependence indicates the presence of large higher-order contributions; past studies of SU (3) theories with N f > 2 have noted the rapid growth of NLO contributions with N f [92][93][94]. The poor convergence of χPT signaled by the large size of these higher-order contributions may be related to the existence of a light scalar which has been omitted from the EFT. It would be interesting to explore this possibility, and the influence of the light scalar's presence on other physical quantities of interest such as the chiral condensate.
On the other hand, we also find some qualitative sim-ilarities between the N f = 8 theory and QCD. In particular, the ratios of ρ and π masses and decay constants appearing in the KSRF relations (Fig. 15) show only mild mass dependence and indicate a vector meson with properties quite similar to that of the QCD ρ, including a large strong decay width. This is an important qualitative feature which should be accounted for in dedicated LHC searches for heavy vector resonances associated with composite Higgs models built on theories such as this one.
Another set of differences between N f = 8 and QCD has been long predicted under the name parity doubling. The expectation is that the ρ and a 1 meson masses will become approximately degenerate as they are parity partners. Similarly the a 0 meson is expected to become substantially lighter as a parity partner of the π [95,96]. In the range of masses we study, we find no evidence for parity doubling between ρ and a 1 . This is in contrast to previously reported results for the N f = 6 theory [96], where degeneracy between the ρ and a 1 masses appeared to set in at light fermion masses. We have not reported any results for the mass of the a 0 state, due to concerns about large systematic effects in our attempts to determine it; such a study is deferred to future work.
Ultimately, a deeper understanding of these results requires investigation of candidate low-energy EFTs appropriate for describing the observed spectrum with a light σ state. Several proposals for such EFTs have been made [97][98][99][100][101][102][103][104][105][106][107][108]; distinguishing between these proposals requires investigation of quantities beyond the spectrum of light states, both on the EFT and lattice sides. Although not strictly an EFT, comparison to holographic predictions for the spectrum [109] would be interesting to expore as well. We have initiated new studies of PNGBs 2 → 2 elastic scattering [17] and form factors, with the goal of exploring which of these EFTs provide the best description of the 8-flavor theory at low energies.
Appendix A: Auto-correlations and topological charge evolution The disconnected correlator D(t) in Eq. (16) is effectively the time correlation of the ψψ and we study its autocorrelation function along the measurement time series at a fixed time separation. In Fig. 16 we show the D(t = 9) history and histogram for the first stream of the am f = 0.00125 N f = 8 ensemble, together with its autocorrelation function. The autocorrelation function is compatible with zero within two standard deviations (dashed grey line) when the distance between measurements is ≈ 20, and the integrated autocorrelation time is ≈ 10. We notice that the autocorrelation function is significantly smaller when we consider the finite time difference of correlators D(dt) = D(t + 1) − D(t) that we fit to extract the 0 ++ meson mass, as can be seen in Fig. 17, with an integrated autocorrelation time ≈ 3. We find that similar observations hold for all ensembles. On all the configurations where we measure the spectrum, we also measure the topological charge, using the clover definition of the discretized F µνF µν gauge tensor. For each ensemble we use gradient flow to smooth the configurations and obtain the topological charge Q at flow time t w satisfying √ 8t w = L/2, where L is the number of points in the spatial directions. For the lightest mass ensembles with N f = 8 and N f = 4, which also have the smallest lattice spacing (cf. Fig. 8), we see a well balanced topological charge, with frequent fluctuations and no sign of topological freezing. The autocorrelation function for the topological charge is also compatible with zero (within two standard deviations) after a distance of a few configurations (grey dashed lines in the plots). The topological charge history, histogram and autocorrelation functions for the individual Monte Carlo streams at am f = 0.00125 for N f = 8 are shown in Fig. 18 to Fig. 20 Here we give additional details on our numerical procedures and results for determination of the 0 ++ σ meson mass. Our procedure for fitting to determine the mass and estimating the fit-range systematic error is described in detail in Secs. IV B and IV C above. Briefly, our procedure for the σ meson is the same as that for any other state, except for two key differences: first, we fit finite time-differences of the correlation functions in order to remove the vacuum contribution; second, we carry out a joint fit of the D(t) and S(t) correlators in order to obtain more stable results for M σ . As a test of the joint fit procedure, we can compare our central results for M σ to determinations from the individual correlators S(t) and D(t). These results are tabulated in Table V and plotted in Fig. 22 and Fig. 23.
There is a clear systematic trend for M σ,S > M σ,D , which is expected due to the fact that contamination from states other than the σ enters with positive sign in S(t) but negative sign in D(t); see Eqs. (19) and (21). The results for M σ from the joint fit are broadly compatible with the individual S(t) and D(t) fits, but offer improved precision compared to the more conservative possibility of taking the difference between M σ,S and M σ,D as a systematic error estimate. Figures 24-28 show detailed numerical results for the 0 ++ joint fit on each 8-flavor ensemble. The top panel shows the best fit versus the "effective correlator" C eff (t) ≡ C(t) A0 e +E0t , where A 0 and E 0 are the best-fit ground state amplitude and energy respectively. The solid line shows the range of t values used in the best fit, while the dashed lines show the extrapolation of the fit. These plots demonstrate that our best-fit models describe the correlator data well, even when extrapolated beyond the fit range. The lower panel in each figure shows the t min dependence of the joint fit, used to determine the fit-range systematic error. Finally, we note that the comparison between the 4-flavor and 8-flavor theories provides one of the more important cross-checks on our 0 ++ mass extraction. The qualitatively different results at N f = 4, with the σ found to be much heavier than the π, are obtained with the same procedure, lattice action, etc. This rules out the possibility that due to some unaccounted-for systematic effect, we are simply finding M π itself from our measurements in the σ channel. Comparison in the 4-flavor theory of the 0 ++ mass spectrum as determined through joint fits to S(t) and D(t) as in our main analysis, and determined with individual fits to S(t) and D(t). A horizontal offset is added to the S(t) and D(t) results to make the plot more easily readable.           In this appendix, we give further detail on possible systematic effects due to the finite lattice spacing and lattice volume in our calculations. We begin with the lattice spacing. For each of the N f = 8 and N f = 4 theories, our lattice calculations are carried out at a single value of the bare gauge coupling. If we adopt a mass-independent lattice scale setting prescription, i.e., a = a(β), then we are unable to carry out a continuum extrapolation. Adopting a mass-dependent scheme a = a(β, m f ) instead, for example using the significant m f dependence observed in t 0 (see Fig. 8), is another possibility; in this case the continuum and chiral (m f → 0) extrapolations are intertwined. We defer such an analysis to future work in which the mass dependence is more closely investigated, since this requires the choice of an appropriate chiral effective theory. A more direct test for discretization effects involves "taste breaking", splitting between the masses of staggered hadrons due to the non-zero lattice spacing. No significant taste breaking is seen in the ρ mesons. The pions, which are particularly sensitive to chiral symmetry breaking, show more significant taste-breaking effects, with masses on the order of 20%-30% heavier than the Goldstone pion. These results are comparable to the 15%-20% splittings seen by the MILC collaboration in QCD for a lattice spacing of a ≈ 0.12 fm [110] using asqtad-improved staggered fermions. We note the partial restoration of taste symmetry at finite a predicted by Lee and Sharpe [111] is present in our results, with the π i5 and π 45 tastes having degenerate masses. We now turn to finite-volume corrections. In order to explicitly test for such effects, we analyze a number of additional ensembles at N f = 8 which are matched to specific ensembles in our main analysis but have smaller or larger lattice volumes. These ensembles are specified in Table VII. Results of our spectrum analysis on these ensembles are given in Table VIII and shown in Figs. 35, 36, and 37. The disconnected diagrams needed to reconstruct the σ correlator were measured only on a subset of these ensembles, so in some cases we do not compute its mass. No significant finite-volume dependence is seen for any of the states in our spectrum, including the σ light scalar, although we only have an explicit test of the latter at the heavy fermion mass am f = 0.00889. In general, direct finite-volume tests are only available on the ensembles with heavier fermion mass, but we note that the volume of the lighter-mass ensembles scales up such that the Taste-split spectrum for the 4-flavor theory in both the pseudoscalar and vector channels. The results are qualitatively similar to what is observed for the 8-flavor theory (Fig. 33), although the pion taste splitting is somewhat smaller here. ensembles used in the main analysis.