text:
Lattice Strong Dynamics (LSD) Collaboration

# Light Scalar Meson and Decay Constant in SU(3) Gauge Theory with Eight
Dynamical Flavors

R. C. Brower  Department of Physics and Center for Computational Science,
Boston University, Boston, Massachusetts 02215, USA  C. Culver  Department of
Mathematical Sciences, University of Liverpool, Liverpool L69 7ZL, UK  K. K.
Cushman  Department of Physics, Sloane Laboratory, Yale University, New Haven,
Connecticut 06520, USA  G. T. Fleming<EMAIL_ADDRESS>Department of
Physics, Sloane Laboratory, Yale University, New Haven, Connecticut 06520, USA
Theoretical Physics Division, Fermilab, Batavia, IL 60510, USA  A. Gasbarro
AEC Institute for Theoretical Physics, University of Bern, 3012 Bern, CH  A.
Hasenfratz  Department of Physics, University of Colorado, Boulder, Colorado
80309, USA  J. Ingoldby  Abdus Salam International Centre for Theoretical
Physics, Strada Costiera 11, 34151, Trieste, Italy  X. Y. Jin  Computational
Science Division, Argonne National Laboratory, Argonne, Illinois 60439, USA
E. T. Neil  Department of Physics, University of Colorado, Boulder, Colorado
80309, USA  J. C. Osborn  Computational Science Division, Argonne National
Laboratory, Argonne, Illinois 60439, USA  E. Owen  Department of Physics and
Center for Computational Science, Boston University, Boston, Massachusetts
02215, USA  C. Rebbi  Department of Physics and Center for Computational
Science, Boston University, Boston, Massachusetts 02215, USA  E. Rinaldi
Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS),
RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  D. Schaich  Department of
Mathematical Sciences, University of Liverpool, Liverpool L69 7ZL, UK  P.
Vranas  Physical and Life Sciences Division, Lawrence Livermore National
Laboratory, Livermore, California 94550, USA  Nuclear Science Division,
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA  E.
Weinberg  Department of Physics and Center for Computational Science, Boston
University, Boston, Massachusetts 02215, USA  NVIDIA Corporation, Santa Clara,
California 95050, USA  O. Witzel  Center for Particle Physics Siegen (CPPS),
Theoretische Physik 1, Naturwissenschaftlich-Technische Fakultät, Universität
Siegen, 57068 Siegen, Germany

###### Abstract

The SU(3) gauge theory with $N_{f}=8$ nearly massless Dirac fermions has long
been of theoretical and phenomenological interest due to the near-conformality
arising from its proximity to the conformal window. One particularly
interesting feature is the emergence of a relatively light, stable flavor-
singlet scalar meson $\sigma$ $(J^{PC}=0^{++})$ in contrast to the $N_{f}=2$
theory QCD. In this work, we study the finite-volume dependence of the
$\sigma$ meson correlation function computed in lattice gauge theory and
determine the $\sigma$ meson mass and decay constant extrapolated to the
infinite-volume limit. We also determine the infinite volume mass and decay
constant of the flavor-nonsinglet scalar meson $a_{0}$.

††preprint: FERMILAB-PUB-23-286-T

##  I Introduction

SU(3) gauge theory with $N_{f}$ flavors of massless Dirac fermions has a
conformal window for $N_{fc}\leq N_{f}\leq 16$ [1, 2]. See [3] for a review of
the early history of constraining $N_{fc}$ and see [4] for the most recent
review. While it is not known with certainty whether the massless $N_{f}=8$
theory is inside or outside the conformal window [5, 6, 7], our collaboration
has previously published results [8, 9] indicating the massless $N_{f}=10$
theory is likely inside the conformal window and for the rest of this paper we
will assume the massless $N_{f}=8$ theory is very close to the edge of the
conformal window. If it is inside the conformal window it is most likely a
very strongly coupled CFT [10]. If it is outside the conformal window,
spontaneous chiral symmetry breaking and confinement produce massless Nambu-
Goldstone bosons and a spectrum of other hadronic states which may be
different relative to QCD due to the proximity of the conformal window.

The continuum SU(3) gauge theory with $N_{f}=8$ Dirac fermions with small
vector-like mass terms is not an IR conformal theory. The small mass terms
explicitly break chiral symmetry, confinement occurs and a massive spectrum of
hadronic states is generated. Another scenario may be possible at stronger
lattice coupling but we don’t consider that here [10].

In our previous papers [11, 12, 13], we identified two specific features of
the low-energy spectrum which were different from QCD. First, the pion decay
constant $F_{\pi}$ strongly depends on the fermion mass unlike QCD where
$F_{\pi}$ is nearly constant with a small, linear correction in the fermion
mass. Second, the flavor-singlet scalar meson $\sigma$ $(J^{PC}=0^{++})$ has a
light mass $M_{\sigma}<1.5\ M_{\pi}$ in the fermion mass region where we
compute it, well below the energy threshold for decay to two pions, whereas in
QCD it is somewhat heavier $M_{\sigma}>1.9\ M_{\pi}$ [14, 15, 16, 17, 18, 19,
20, 21, 22], remaining just below decay threshold across in an equivalent
fermion mass range. The QCD picture is somewhat consistent with our earlier
$N_{f}=4$ calculation [12]. However, we also identified several features of
the $N_{f}=8$ theory which appeared similar to QCD calculations in an
equivalent range of fermion masses: the ratios $M_{\rho}/F_{\pi}$ and
$M_{\text{nucleon}}/F_{\pi}$, the decay constants $F_{\pi}$, $F_{\rho}$ and
$F_{a_{1}}$ appear consistent with QCD KSRF relations [23, 24] and the $I=2$
$\pi\pi$ scattering length $a_{\pi\pi}$ appears to agree with QCD.

In this paper, we focus larger volume calculations at the various fermion
masses (App. A) which will allow us to extrapolate our results to the infinite
volume limit, removing one potential source of systematic error (Sec. II.5).
We also introduce an improved method for analyzing $\sigma$ meson correlation
functions with a new subtraction scheme in the rest frame combined with
simultaneous analysis in several moving frames (Secs. II.2 and II.3). We rely
heavily on the method of Bayesian model averaging [25] (Sec. III). We present
new calculations of the flavor-singlet scalar decay constant $F_{S}$ and the
flavor-nonsinglet scalar meson $a_{0}$ mass and decay constant $F_{a_{0}}$
(Sec. II.4). We also comment briefly on the chiral condensate and its
contribution to the Gell-Mann-Oakes-Renner (GMOR) relation and its
generalizations (Sec. II.6).

As described in our earlier paper [12] we have chosen the bare lattice
parameter $\beta=4.8$ such that the lattice spacing $a$ is as coarse as
possible given our current action, so that we can get as close to the chiral
limit $am_{q}\to 0$ as possible with available computing resources. We are
working on calculations at $am_{q}=0.00056$ which may provide further insights
in the near future. If the massless $N_{f}=8$ theory is conformal and
sufficiently strongly coupled [10] then it is likely a new lattice action that
allows for even coarser lattice spacings will be necessary to make further
progress.

Phenomenologically, theories that exhibit approximate conformal behavior at
strong coupling are anticipated to produce large anomalous dimensions over a
wide interval of scales, which can make them attractive as candidate composite
Higgs models [26, 27, 28, 29, 30]. In particular, the
$\text{SU}(3)\,\,N_{f}=8$ theory has been used to build composite Higgs models
in [31, 32, 33]. The construction of a low energy EFT for the lightest
composites, to which the rest of the Standard Model can be coupled is a
crucial intermediate step in the creation of these models. In a separate
letter [34], we discuss various effective models that can be fit to our data.

##  II Staggered two-point correlation function construction and modeling

###  II.1 Staggered two-point correlation function construction

The continuum $N_{f}=8$ theory is approximated on a finite lattice by an
SU(2)-doublet of staggered fermion fields $(\chi_{1}\ \chi_{2})$ that carry
only an SU(3) color index at each lattice site. Each component of the doublet
represents four non-degenerate Dirac fermion “tastes” with spin and taste
degrees of freedom spread out over $2^{4}$ sites of local hypercubes. In the
continuum limit where the bare gauge coupling $g_{0}^{2}\to 0$, these tastes
become degenerate and equivalent to four Dirac flavors. Hence the doublet of
staggered fields becomes a degenerate $N_{f}=8$ theory in this continuum
limit. Staggered fermions have a remnant of chiral symmetry that can lead to
$(N_{f}/4)^{2}$ Nambu–Goldstone bosons when taking the massless chiral limit
at finite lattice spacing, assuming the chiral symmetry is spontaneously
broken by the gauge interactions. However, to recover the full flavor symmetry
it is essential to take the $g_{0}^{2}\to 0$ continuum limit prior to the
$m_{q}\to 0$ chiral limit.

In general, a staggered meson two-point correlaton function where source and
sink operators have the same quantum numbers $Q$ is (schematically)

$C_{Q}(\vec{p},|t-t_{0}|)=\left\langle\sum_{\vec{x}}e^{i\vec{p}\cdot(\vec{x}-\vec{x}_{0})}\overline{\chi}(\vec{x}+\vec{\delta}^{\prime},t)\Gamma_{Q}(\vec{x},\vec{\delta}^{\prime})\tau\chi(\vec{x},t)\
\overline{\chi}(\vec{x}_{0}+\vec{\delta},t_{0})\Gamma_{Q}(\vec{x}_{0},\vec{\delta})\tau\chi(\vec{x},t)\right\rangle$
(1)

where $\Gamma_{Q}(\vec{x},\vec{\delta})$ are phases which refer to the spin-
taste structure of the interpolating operators with quantum numbers $Q$, and
$\tau$ is either an SU(2) generator for a non-singlet correlator or a $2\times
2$ identity matrix for a singlet correlator under the SU(2) staggered doublet
symmetry. There are various phase conventions possible, one common choice is
[35]. Not shown are gauge matrices required to make the whole thing gauge-
invariant, e.g. connecting sites $(\vec{x},t)$ and
$(\vec{x}+\vec{\delta}^{\prime},t)$. Also, translation invariance of the
ensemble average $\langle\ \cdot\ \rangle$ guarantees the correlation function
depends only on the distance $|t-t_{0}|$ and not the source position
$(\vec{x}_{0},t_{0})$.

In our earlier paper [12], the LSD Collaboration constructed correlation
functions with $\vec{p}=0$ for local and point-split operators. In this study,
we focused on constructing singlet and non-singlet correlation functions of
local operators at several different momenta $\vec{p}$ with much higher
statistics. On each gauge configuration, we generate a unique set of $N$
random source points $(\vec{x}_{0},t_{0})_{n}$ and construct a primitive
staggered meson “connected” correlator

$C(\vec{x},t)=\frac{1}{N}\sum_{n}\mathcal{T}_{n}\
\text{Tr}_{\text{color}}\left[G_{F}(\vec{x}_{0},t_{0};\vec{x},t)\
G_{F}^{\dagger}(\vec{x},t;\vec{x}_{0},t_{0})\right]$ (2)

where $G_{F}(\vec{x}_{0},t_{0};\vec{x},t)$ is a $3\times 3$ color matrix of
single staggered fermion propagator from the site $(\vec{x}_{0},t_{0})$ to the
site $(\vec{x},t)$ and $\mathcal{T}_{n}$ represents the translation of the
$n$-th source location $(\vec{x}_{0},t_{0})_{n}$ to the origin $(\vec{0},0)$.
Then, we record the value of this averaged primitive correlator for every sink
point $(\vec{x},t)$ in the lattice volume. We refer to this as a connected
correlator because valence fermion lines connect the source and sink points.
With post-processing we can project this primitive correlator into eight
different non-singlet staggered meson quantum number channels of different
momenta $\vec{p}$ using Fourier transform

$C_{Q}(\vec{p},t)=\sum_{\vec{x}}e^{i\vec{p}\cdot\vec{x}}\ C(\vec{x},t)\
\phi_{Q}(\vec{x}).$ (3)

For example, if we choose the phase $\phi_{Q}(\vec{x})=1$, we get the
correlation function for the $\pi_{5}$ meson, which is the pseudo-
Nambu–Goldstone boson.

Also in our earlier work, we explained in detail how we construct a
“disconnected” correlator, where valence fermion lines do not connect the
source and sink points

$D(\vec{p},|t-t_{0}|)=\sum_{\vec{x}}\sum_{\vec{x}_{0}}e^{i(\vec{x}-\vec{x}_{0})\cdot\vec{p}}\
\text{Tr}\left[G_{F}(\vec{x}_{0},t_{0};\vec{x}_{0},t_{0})\right]\
\text{Tr}\left[G_{F}(\vec{x},t;\vec{x},t)\right]$ (4)

using a diluted noisy estimator to compute the trace at each site on the
lattice for each gauge configuration, which is again recorded as a single
value per site in the lattice volume. With post-processing, we can compute the
disconnected correlator for any spatial momentum $\vec{p}$ using FFT and the
fast convolution algorithm

$\widetilde{O}(\vec{p},\omega)=\sum_{\vec{x},t}e^{i(\vec{p}\cdot\vec{x}+\omega
t)}\ \text{Tr}\left[G_{F}(\vec{x},t;\vec{x},t)\right]$ (5)
$D(\vec{p},t)=\sum_{\omega}e^{-i\omega
t}\left|\widetilde{O}(\vec{p},\omega)\right|^{2}$ (6)

where the result is automatically invariant under any lattice translation. In
an $N_{f}$-flavor theory, the flavor-singlet scalar correlator for the
$\sigma$ meson is then

$C_{\sigma}(\vec{p},t)=\left(\frac{N_{f}}{4}\right)^{2}D(\vec{p},t)-\frac{N_{f}}{4}C_{a_{0,1}}(\vec{p},t)$
(7)

and where $C_{a_{0,1}}(\vec{p},t)$ is the flavor-nonsinglet scalar meson
correlator constructed from Eq. (3) with the appropriate choice of phases.
Note this normalization is different from [12] where we dropped an overall
factor of $N_{f}/4$.

Regarding the naming convention of mesons, we note that the continuum SU(8)
flavor representation is broken by lattice artifacts to a subgroup SU(2)
$\times$ taste, a discrete subgroup of SU(4). Meson names will follow the PDG
convention for two-flavor mesons: $\pi$, $a_{0}$, $\rho$, …plus an additional
subscript to indicate the representation under the discrete taste group:
$\pi_{5}$, $a_{0,1}$, $\rho_{i}$, …There is only one scalar meson which is a
singlet over the whole flavor group which we name $\sigma$ and no subscript is
required. The effects of taste-breaking were discussed previously [12] and we
will not expand on it, so the taste index not play a significant role here
with one important exception. In the continuum two-flavor theory, the decay
$a_{0}\to\pi\pi$ is forbidden by isospin symmetry. However, in our staggered
$N_{f}=8$ theory, the decay $a_{0,1}\to\pi_{5}\pi_{5}$ is allowed because the
$\pi_{5}$ and the $a_{0,1}$ are not in the same SU(2) flavor subgroup, as
indicated by the different taste indices. It is analogous to the decay of
$a_{0}\to KK$ in continuum three-flavor theory.

###  II.2 Model for Staggered Meson Correlation Functions

We will consider three different types of models for staggered meson two-point
correlation functions in this paper. As we are employing Bayesian model
averaging, further discussed in Sec. III, we don’t have to choose a particular
model but rely on the computed model probabilities to distinguish the most
likely models for a given correlation function. Within each model type, the
number of free parameters in each specific instance of the model will depend
upon the number of oscillating and non-oscillating states included.

The model we will use for the staggered meson correlation function in the time
domain (Model A) is

$\displaystyle C(\vec{p},t)$ $\displaystyle=$ $\displaystyle
c_{0}\delta_{\vec{p},0}+\sum_{n}\frac{c_{n}}{2\left(1-e^{-E_{n}N_{t}}\right)\sinh(E_{n})}\left[e^{-E_{n}t}+e^{-E_{n}(N_{t}-t)}\right]$
(8)
$\displaystyle+(-1)^{t}\sum_{j}\frac{c^{\prime}_{j}}{2\left(1-e^{-E^{\prime}_{j}N_{t}}\right)\sinh(E^{\prime}_{j})}\left[e^{-E^{\prime}_{j}t}+e^{-E^{\prime}_{j}(N_{t}-t)}\right]$

where we have chosen to use a particular “relativistic” normalization for the
amplitudes. As is typical for the staggered fermions, there are a set of
states labeled by $n$ whose contributions do not oscillate in time and another
set of states labeled by $j$, with different quantum numbers, that oscillate
in time with a factor $(-1)^{t}$. The energies $E_{n}$ and $E_{j}^{\prime}$
are understood to depend implicitly on the spatial momentum $\vec{p}$. We also
allow for the possibility of a $t$-independent contribution to the correlation
function, $c_{0}$, which is generally not present for flavor-nonsinglet
correlation functions due to translation invariance of the ensemble average.
But, it is the dominant contribution to the flavor-singlet $\sigma$
correlation function and must be treated carefully in order to extract
reliable estimates of model parameters. Note that the constant only
contributes to the $\vec{p}=0$ correlator, so one method of dealing with this
constant is to work with $\vec{p}\neq 0$ correlators. Given that we are
interested in the energy of the $\sigma$ meson in the rest frame,
$\lim_{\vec{p}\to 0}E_{\sigma}(\vec{p})=M_{\sigma}$, this approach requires a
good understanding of the dispersion relation on the lattice.

To motivate the normalization of amplitudes $c_{n}$ and $c_{j}^{\prime}$ in
Eq. (8), we can perform the discrete cosine transform (DCT-I) of the time-
domain correlation function into the frequency domain analytically

$\widetilde{C}(\vec{p},k)=c_{0}\delta_{\vec{p},0}\delta_{k,0}+\frac{1}{N_{t}}\sum_{n}\frac{c_{n}}{\widehat{E}_{n}^{2}+\widehat{\omega}_{k}^{2}}+\frac{1}{N_{t}}\sum_{j}\frac{c^{\prime}_{j}}{\widehat{E}^{\prime
2}_{j}+\widehat{\omega}^{\prime 2}_{k}}$ (9)

where

$\widehat{E}_{n}=2\sinh\frac{E_{n}}{2},\qquad\widehat{\omega}_{k}=2\sin\frac{2\pi
k}{2N_{t}},\qquad\widehat{\omega}^{\prime}_{k}=2\sin\left(\frac{\pi}{2}-\frac{2\pi
k}{2N_{t}}\right).$ (10)

Comparing the expression for two different spatial momenta $\vec{p}$, the
energies $E_{n}$ and $E_{j}^{\prime}$ will be different, defining some lattice
dispersion relation. But the amplitudes $c_{n}$ and $c_{j}^{\prime}$ are
momentum independent as normalized, and therefore frame-independent as
expected in a Lorentz-invariant theory, hence a “relativistic” normalization.

In our previous work, we considered another method of dealing with the
constant $c_{0}$ which was to analyze the finite difference correlation
function for the $\vec{p}=0$ $\sigma$ meson

$\Delta_{\sigma}(t)=C_{\sigma}(t+1)-C_{\sigma}(t).$ (11)

In the model, the cancellation of $c_{0}$ is exact but in our lattice
calculation there is inherent statistical noise contributing to each time
slice, so the cancellation is not exact. In this work, we propose an improved
subtraction scheme for $\vec{p}=0$ correlation functions

$\overline{C}(t)=C(t)-\frac{1}{N_{t}}\sum_{t^{\prime}=0}^{N_{t}-1}C(t^{\prime})$
(12)

for states that have a time-independent part, like the $\sigma$ meson. Given
our frequency analysis above, we can see the subtraction is the zero-frequency
component of the correlation function $\overline{C}(t)=C(t)-\widetilde{C}(0)$.
Furthermore, we know explicitly the functional form of the residual constant
that comes from the integral of the $t$-dependent part of the correlation
function

$c_{0}-\widetilde{C}(0)=-\frac{1}{N_{t}}\sum_{n=1}^{\infty}\frac{c_{n}}{\widehat{M}_{n}^{2}}-\frac{1}{N_{t}}\sum_{j=1}^{\infty}\frac{c^{\prime}_{j}}{4+\widehat{M}^{\prime
2}_{j}}$ (13)

Because some of the fit parameters appear in the residual constant, we will
include that part in the fit and shift the constant (Model B)

$\displaystyle\overline{C}(t)$ $\displaystyle=$
$\displaystyle\overline{c}_{0}+\sum_{n}^{n_{\text{max}}}\frac{c_{n}}{2\left(1-e^{-M_{n}N_{t}}\right)\sinh(M_{n})}\left[e^{-M_{n}t}+e^{-M_{n}(N_{t}-t)}\right]-\frac{c_{n}}{N_{t}\widehat{M}_{n}^{2}}$
(14)
$\displaystyle+\sum_{j}^{j_{\text{max}}}\frac{(-1)^{t}c^{\prime}_{j}}{2\left(1-e^{-M^{\prime}_{j}N_{t}}\right)\sinh(M^{\prime}_{j})}\left[e^{-M^{\prime}_{j}t}+e^{-M^{\prime}_{j}(N_{t}-t)}\right]-\frac{c_{j}^{\prime}}{N_{t}(4+\widehat{M}_{j}^{\prime
2})}$ $\displaystyle\overline{c}_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{N_{t}}\sum_{n=n_{\text{max}}+1}^{\infty}\frac{c_{n}}{\widehat{M}_{n}^{2}}-\frac{1}{N_{t}}\sum_{j=j_{\text{max}}+1}^{\infty}\frac{c^{\prime}_{j}}{4+\widehat{M}^{\prime
2}_{j}}.$

In counting free parameters, model B will have one more free parameter than
model A and the interpretation of the value of this parameter,
$\overline{c}_{0}$, will depend strongly on the choice of $n_{\text{max}}$ and
$j_{\text{max}}$. In particular, we expect $\overline{c}_{0}\to 0$ within
statistical uncertainties as the number of states included in a particular
model instance approaches the limit of available statistics to properly
constrain them.

We will also consider a modification of model B (Model C) where we constrain
$\overline{c}_{0}=0$. It will have the same number of free parameters as model
A in an instance where they include the same number of states. In the context
of Bayesian model averaging, we expect that model B will have a higher
relative probability than model C in instances where $\overline{c}_{0}$ is
statistically non-zero. But with increasing numbers of states eventually model
C should become more probable, also indicating the limit in which the power of
the available statistics to constrain parameters has been exhausted.

###  II.3 Staggered Meson Dispersion Relation

The functional momentum dependence of energies $E_{Q}(\vec{p})$ extracted from
analysis of two-point correlation functions $C_{Q}(\vec{p},t)$ is a
complicated, non-perturbative problem because Lorentz symmetry is broken by
the lattice discretization so the theory is not invariant under boosts. Still,
Lorentz symmetry is fully recovered in the continuum limit. Naively we can
expect

$a^{2}E_{Q}^{2}(\vec{p})=a^{2}M_{Q}^{2}+a^{2}p^{2}+O(a^{4}p^{4})$ (15)

where we explicitly show the lattice spacing $a$ in this dimensionless
relation and define the spatial momentum components $p_{i}=2\pi
n_{i}/(aN_{s})$ and
$n_{i}\in\left\\{-N_{s}/2+1,\cdots,0,\cdots,N_{s}/2\right\\}$ and $N_{s}$ is
the number of lattice sites in the spatial directions.

To improve upon this estimate, one would have to understand the dynamics on
the lattice of the eigenstates corresponding to these energies. This is a
challenging problem since the eigenstates are not simple single-hadron
excitations, in general, but are more likely strongly-interacting multi-hadron
states. But, the lowest energy state with given quantum numbers $Q$ may
reasonably be expected to behave like a single-hadron state, particularly if
it’s energy is well below the nearest multi-hadron threshold. In this case, we
can approximate the dispersion relation with that of a non-interacting boson
on the lattice [36]

$\widehat{E}_{Q}^{2}=\widehat{M}_{Q}^{2}+\widehat{p}^{2}+O(\widehat{p}^{4})$
(16)
$\widehat{E}_{Q}=2\sinh\frac{aE_{Q}}{2},\quad\widehat{M}_{Q}=2\sinh\frac{aM_{Q}}{2},\quad\widehat{p}_{i}=2\sin\frac{ap_{i}}{2}$
(17)

In the second equation, we have explicitly put in the lattice spacing
dependence $a$. Both lattice dispersion relations correspond to the same
continuum relation as $a\to 0$.

In either of these models, Eqs. (15) or (17), the finite size of the lattice
along spatial directions $N_{s}$ directly controls the spacing between the
discrete momenta but is not expected to appear explicitly in the finite
lattice spacing corrections $O(a^{4}p^{4})$ or $O(\widehat{p}^{4})$. When we
fit our lattice data on two or more volumes at the same value of the bare
coupling and mass, we will parameterize our fits so that the same lattice
corrections are used on all volumes.

###  II.4 Staggered Meson Decay Constants

The normalization in Eq. (8) was chosen such that $c_{n}\to\left|\left\langle
0\left|\mathcal{O}\right|n,\vec{p}=0\right\rangle\right|^{2}$ in the continuum
limit with the usual continuum relativistic normalization. Following Eq. (7.5)
of [37] we define the pion decay constant

$\sqrt{2}\widehat{F}_{\pi_{5}}\left(\widehat{E}_{\pi_{5}}^{2}-\widehat{p}^{2}\right)=2m_{q}\frac{1}{\sqrt{N_{f}}}\left\langle
0\left|P_{5}\right|\pi_{5}(\vec{p})\right\rangle\quad\implies\quad\widehat{F}_{\pi_{5}}=\frac{1}{\sqrt{2}}\
\frac{m_{q}\sqrt{|c_{\pi_{5}}|}}{\widehat{E}_{\pi_{5}}^{2}-\widehat{p}^{2}}$
(18)

where $N_{f}$ in this equation is the number of continuum flavors of a single
staggered fermion, i.e. $N_{f}=4$. Note we put the hat on the symbol for
$\widehat{F}_{\pi_{5}}$ to indicate the form of the lattice dispersion
relation used. We could have just as easily used the other form of the lattice
dispersion relation, which would lead to a slight different definition of the
decay constant. Both definitions should converge to the continuum one in the
limit of zero lattice spacing. This definition is slightly different than ones
previously used by the LSD Collaboration for the pion decay constant [12, 13],
but the difference is not statistically significant.

For the isotriplet scalar form factor, there does not seem to be a
conventional normalization [38, 39] for the decay constant in QCD as it is an
unstable resonance. See review “Scalar Mesons below 1 GeV” in [40]. In our
$N_{f}=8$ theory over the range of fermion masses we’ve studied, the non-
singlet scalar meson appears to be stable, although close in energy to its
decay threshold. We choose to normalize it analogously with the pion decay
constant

$\widehat{F}_{a_{0,1}}=\frac{1}{\sqrt{2}}\
\frac{m_{q}\sqrt{|c_{a_{0,1}}|}}{\widehat{E}_{a_{0,1}}^{2}-\widehat{p}^{2}}$
(19)

where $c_{a_{0,1}}$ is the residue of the first pole in the frequency domain
representation of the non-singlet scalar two-point correlation function, Eq.
(9).

For the isosinglet scalar decay constant, we use the normalization defined in
Eq. (72) of [41]

$\widehat{F}_{S}\
\left(\widehat{E}_{\sigma}^{2}-\widehat{p}^{2}\right)=m_{q}\left\langle
0\left|S(0,0)\right|\sigma(\vec{p})\right\rangle$ (20)

where the scalar current is defined as
$S(\vec{x},t)=\sum_{i=1}^{N_{f}/4}\overline{\chi}_{i}(\vec{x},t)\chi_{i}(\vec{x},t)$.
The two-point correlation function of this scalar current is defined in Eq.
(7) and, in terms of this correlation function, the decay constant is defined

$\widehat{F}_{S}=\frac{m_{q}\sqrt{|c_{\sigma}|}}{\widehat{E}_{\sigma}^{2}-\widehat{p}^{2}}.$
(21)

In particular, the normalization used in Eq. (7) is essential to correctly
normalizing the decay constant.

###  II.5 Finite Volume Corrections

In QCD, finite volume corrections to the pion mass and pion decay constant
extracted from a two-point correlation function calculated on a periodic torus
of spatial size $L$ can be computed in chiral perturbation theory provided
$M_{\pi}L\gg 1$ and $F_{\pi}L\gg 1$. See Eq. (6.15) of [42], for example. In
$N_{f}=8$ over the range of fermion masses for which we have relevant lattice
calculations, chiral perturbation theory is unlikely to be a good effective
description for two reasons: the strong fermion mass dependence of $F_{\pi}$
and the stable $\sigma$ meson with $M_{\sigma}\ll 4\pi F_{\pi}$. So it is not
expected that finite volume corrections computed in chiral perturbation theory
(ChiPT) will exactly match the numerical calculations. Still, it seems likely
that whatever low energy effective theory replaces ChiPT will have much the
same structure as these arise from contributions of virtual pion degrees of
freedom that probe the finite volume by wrapping the spatial cycles of the
torus, and the pion still is the lightest hadron in the eight-flavor theory.
There may be additional contributions from $\sigma$-meson degrees of freedom,
but they are expected to be sub-leading due to the somewhat heavier mass.

We will follow the approach used in [13] and use ChiPT-inspired forms to model
our finite volume corrections

$\displaystyle M_{Q}(L)$ $\displaystyle=$ $\displaystyle
M_{Q}(\infty)\left[1+\alpha_{Q}\frac{M_{\pi}^{2}}{(4\pi
F_{\pi})^{2}}\sum_{n=1}^{\infty}\frac{4\ \kappa(n)}{\sqrt{n}\
M_{\pi}L}K_{1}(\sqrt{n}\ M_{\pi}L)\right]$ (22) $\displaystyle F_{Q}(L)$
$\displaystyle=$ $\displaystyle
F_{Q}(\infty)\left[1+\beta_{Q}\frac{M_{\pi}^{2}}{(4\pi
F_{\pi})^{2}}\sum_{n=1}^{\infty}\frac{4\ \kappa(n)}{\sqrt{n}\
M_{\pi}L}K_{1}(\sqrt{n}\ M_{\pi}L)\right].$ (23)

The function $\kappa(n)$ counts the number of lattice vectors $\vec{n}$ with
integer-valued components of length $\sqrt{n}$, see Tab. 1. In QCD, it is
common to expand the sum over modified Bessel functions $K_{1}$, assuming
$M_{\pi}L\gg 1$ and keep only the leading term, particularly if
$M_{\pi}L\gtrsim 4$ in all, leading to

$\sum_{n=1}^{\infty}\frac{4\ k(n)}{\sqrt{n}\ M_{\pi}L}K_{1}(\sqrt{n}\
M_{\pi}L)\approx\frac{12\sqrt{2\pi}}{(M_{\pi}L)^{3/2}}e^{-M_{\pi}L}$ (24)

In an earlier paper [13], we also used this approximation for the finite
volume extrapolation of $M_{\pi}$ and $F_{\pi}$. We did not observe any
significant change in the result if we included more terms in the expansion.
In this analysis, we will be conservative and not expand the modified Bessel
functions and truncate the sum only after the first eight terms (up to $n=8$)
although we expect it will not make a significant difference relative to
keeping just the leading term.

$n$ | $\vec{n}$ | $|\vec{n}|$ | $\kappa({n})$
---|---|---|---
0 | (0,0,0) | 0 | 1
1 | (1,0,0) | 1 | 6
2 | (1,1,0) | $\sqrt{2}$ | 12
3 | (1,1,1) | $\sqrt{3}$ | 8
4 | (2,0,0) | 2 | 6
5 | (2,1,0) | $\sqrt{5}$ | 24
6 | (2,1,1) | $\sqrt{6}$ | 24
7 | — | — | 0
8 | (2,2,0) | $\sqrt{8}$ | 12
9 | (2,2,1) | 3 | 24
$9^{\prime}$ | (3,0,0) | 3 | 6
Table 1:  The number of lattice vectors $\vec{n}$ with integer-valued
components of length $\sqrt{n}$. Note there are no vectors of length
$\sqrt{7}$ and, starting at length 3, there may be multiple inequivalent sets
of vectors under the cubic group.

Since the infinite volume extrapolation described in this section implicitly
assumes that the pion is a pseudo-Nambu-Goldstone boson, one should use
caution when modeling the extrapolated data provided later in this paper,
particularly if one wants to explore other finite volume corrections, e.g. due
to a light isosinglet scalar. If one assumes that the massless limit of the
theory approaches a conformal fixed point, possible finite volume corrections
were discussed in [43]. In either case, one should use the finite volume data
provided in the supplementary materials [44] when performing further analysis.

###  II.6 The GMOR Relation and Near-Conformality

As a guide to constructing low energy effective descriptions $N_{f}=8$ theory,
it would be useful to characterize the extent to which one or a few light
states dominates the low-energy dynamics. An important phenomenological tool
for characterizing the degree to which the dynamics of the Nambu–Goldstone
pions dominates low-energy phenomena in QCD was first described by Gell-Mann,
Oakes and Renner (GMOR) [45]. In their original derivation, they a priori
assumed pion-pole dominance and derived the GMOR relation as a consequence.
Our derivation will not initally assume pole-dominance but start with the
integral of the axial Ward–Takehashi identity. In our notation, this can be
written

$\sum_{t=0}^{N_{t}-1}C_{\pi_{5}}(\vec{0},t)=\frac{1}{m_{q}}\text{Tr}_{\text{color}}\left[G_{F}(\vec{0},0;\vec{0},0)\right]$
(25)

which is an exact spectral identity on each gauge configuration, not just in
the ensemble average. In the chiral limit $m_{q}\to 0$ of a theory with
spontaneous chiral symmetry breaking, the trace on the right hand side will
approach a constant following the Banks–Casher relation [46, 47] and the
integrated pion correlation function will diverge due to the massless
Nambu–Goldstone pion. Using Eqs. (9) and (18) we can identify the rate of this
divergence with parameters in our fit functions

$\sum_{t=0}^{N_{t}-1}C_{\pi_{5}}(\vec{0},t)\to\frac{c_{\pi_{5}}}{\widehat{M}_{\pi_{5}}^{2}}=2\frac{\widehat{F}_{\pi_{5}}^{2}\widehat{M}_{\pi_{5}}^{2}}{m_{q}^{2}}\quad\text{as}\quad
m_{q}\to 0.$ (26)

Using the normalization of the isosinglet scalar current in Eq. (20) leads to
a generalization of the GMOR relation for general $N_{f}$

$m_{q}\left\langle
S\right\rangle=m_{q}\frac{N_{f}}{4}\left\langle\overline{\chi}\chi\right\rangle\geq\frac{N_{f}}{2}\widehat{F}_{\pi_{5}}^{2}\widehat{M}_{\pi_{5}}^{2}.$
(27)

Now, if we assume spontaneous symmetry breaking and pion pole-dominance, the
inequality becomes an equality in the limit $m_{q}\to 0$ and $\left\langle
S\right\rangle$ approaches a well-defined low energy constant, which is the
usual GMOR relation.

Patella [48] has noted that Eq. (27) should also be true in a mass-deformed
CFT with a large mass anomalous dimension $(1<\gamma^{*}<2)$ due to a large
contribution to the pion correlation function generated by the running of the
mass. They propose examining the GMOR ratio

$R_{G}(m_{q})\equiv\frac{m_{q}\left\langle\overline{\chi}\chi\right\rangle}{2\widehat{F}_{\pi_{5}}^{2}\widehat{M}_{\pi_{5}}^{2}}=\left\\{\begin{array}[]{ll}1\
,&\text{(near-conformal)}\\\ 1<R_{G}(0)<\infty\
,&(\text{CFT},1<\gamma^{*}<2)\\\ \infty\
,&(\text{CFT},0<\gamma^{*}<1)\end{array}\right.\quad\text{as}\quad m_{q}\to 0$
(28)

for an indication of whether the theory is near-conformal or conformal in the
chiral limit. In a near-conformal scenario, it is not clear at what fermion
mass $m_{q}$ one would expect to see the transition from the approximately
hyperscaling regime where $R_{G}(m_{q})>1$ to the spontaneously broken regime
where $R_{G}(m_{q})\to 1$ as $m_{q}\to 0$. Just observing $R_{G}(m_{q})>1$ at
some finite fermion mass is not sufficient to establish IR conformality. In
particular, one must follow the correct order of limits: volume to infinity,
lattice spacing to zero, and then fermion mass to zero.

##  III Bayesian Model Averaging

###  III.1 General Setup

One of the challenges observed in our previous analysis of the light meson
spectrum in the $N_{f}=8$ theory [12] were large systematic errors due to fit
parameters varying significantly over a range of different fits while
$\chi^{2}/\text{dof}$ did not. We define $\log p(D|M)$ by the usual chi-
squared prescription

$\log p(D|M)\propto-\frac{1}{2}\sum_{t,t^{\prime}\in
T_{1}}\left(C(t)-f_{M}(t)\right)\Sigma_{tt^{\prime}}^{-1}\left(C(t^{\prime})-f_{M}(t^{\prime})\right)$
(29)

where $C(t)$ is correlation function computed from our lattice ensemble $D$,
$f_{M}(t)$ is the function for model $M$ to be fitted by minimizing the log-
likelihood, $T_{1}$ is the subset of times selected for fitting and
$\Sigma_{tt^{\prime}}$ is the covariance of the correlation function $C(t)$ on
the subset $T_{1}$. Assuming all the quantities are properly estimated from
the ensemble, the log-likelihood is expected to sample the chi-squared
distribution for degrees of freedom equal to the number of times in $T_{1}$
minus the number of free parameters in $M$.

Subsequent to our earlier analysis, Jay and Neil proposed [25] a Bayesian
model averaging analysis framework which estimates $\log p(M|D)$, the
probability that a model $M$ is a good representation of the data selection
$D$. One suggested estimator of the model probability is based on the Akaike
Information Criterion (AIC), provided nuisance model parameters are assigned
to account for data subsets not included in the fit. For example, let $M$ be a
model with $N_{M}$ free parameters and the maximal data set has $N_{T}$ times
available to be fit. If we perform the fit only on a subset of times $T_{1}$
of size $N_{1}$ then the number of data points not included
$N_{0}=N_{T}-N_{1}$ must be assigned nuisance parameters. Thus, for the AIC,
the number of relevant parameters is $N_{M}+N_{0}$ and the model probability
[49] is

$\log p(M|D)\propto\log p(D|M)-(N_{M}+N_{0})$ (30)

After the model probability has been estimated for the full set of models
$\\{M\\}$ to be considered for the analysis, we normalize this set of
probabilities: $\sum_{\\{M\\}}p(M|D)=1$. In App. B, we provide some details
how we perform this sum accurately given the potential for widely varying
values of $\log p(M|D)$.

With an reasonable estimate of the model probability, it seems straightforward
to construct expectation values and variances of model parameters over the set
of possible models considered. For example, the expected value of a model
parameter is

$\displaystyle E(a)$ $\displaystyle=$
$\displaystyle\frac{1}{\Sigma_{1}}\sum_{\\{M|a\in M\\}}a_{M}\ p(M|D)\
\Theta\left[p(M|D)-p_{\text{cut}}\right]$ $\displaystyle\Sigma_{1}$
$\displaystyle=$ $\displaystyle\sum_{\\{M|a\in M\\}}p(M|D)\
\Theta\left[p(M|D)-p_{\text{cut}}\right]$ (31)

where, to make sure the notation is clear, we compute a weighted average over
only the subset of models that contain the parameter, $\\{M|a\in M\\}$ and
further consider only models where the model weight is greater than some pre-
determined minimum $p_{\text{cut}}$, as enforced by the Heavyside function
$\Theta$.

The variance of the model-averaged expectation value has two contributions.
The first, and usually dominant, contribution is the weighted average over
models of the square of the error estimate $\sigma_{a,M}$ for the parameter
$a_{M}$ in a given model $M$

$E(\sigma^{2}_{a})=\frac{1}{\Sigma_{1}}\sum_{\\{M|a\in M\\}}\sigma^{2}_{a,M}\
p(M|D)\ \Theta\left[p(M|D)-p_{\text{cut}}\right]$ (32)

The second, usually sub-dominant, contribution is the weighted variance of the
model estimates of parameters $a_{M}$, relative to the model-averaged
expectation $E(a)$

$\displaystyle\text{Var}(a)$ $\displaystyle=$
$\displaystyle\frac{\Sigma_{1}}{\Sigma_{1}^{2}-\Sigma_{2}}\sum_{\\{M|a\in
M\\}}\left(a_{M}-E(a)\right)^{2}\ p(M|D)\
\Theta\left[p(M|D)-p_{\text{cut}}\right]$ $\displaystyle\Sigma_{2}$
$\displaystyle=$ $\displaystyle\sum_{\\{M|a\in M\\}}p(M|D)^{2}\
\Theta\left[p(M|D)-p_{\text{cut}}\right]$ (33)

The final error estimate for the model average of a parameter is to add the
two contributions in quadrature

$\sigma_{a}=\sqrt{E(\sigma^{2}_{a})+\text{Var}(a)}$ (34)

Now we can discuss the motivation behind the probability cut $p_{\text{cut}}$.
In our experience, the model-averaged $E(a)$ tend to be dominated by a few
choices whose $p(M|D)\sim O(1)$. It seems reasonable to expect that
$E(\sigma_{a}^{2})$ should be similarly dominated by choices whose $p(M|D)\sim
O(1)$ and not $p(M|D)\sim O(p_{\text{cut}})$. However, we have observed cases
of overfitting for certain models where as the data selection changes such
that $p(M|D)$ decreases, $\sigma_{a,M}^{2}$ increases at a faster rate,
leading to those very unlikely model choices to dominate the model average of
the squared error $E(\sigma_{a}^{2})$. $p_{\text{cut}}$ can be adjusted to
minimize the impact of this scenario.

To understand how this can happen, we recall that uncertainty of a two-point
meson correlation function grows exponentially in Euclidean time [50]

$\text{Var}\left[C_{Q}(\vec{p},t)\right]\sim\exp\left[2\left(E_{Q}(\vec{p})-M_{\pi_{5}}\right)t\right]$
(35)

Now, for a given model function $M$ with its fixed number of exponential
terms, there is a certain $t_{\text{min}}$ for which
$-\chi^{2}/2\sim(N_{M}-t_{\text{max}}+t_{\text{min}}-1)/2$, indicating a good
fit using the usual chi-squared criteria $\chi^{2}/\text{dof}\sim 1$. For fits
on the interval $[t,t_{\text{max}}],t<t_{\text{min}}$, there will be no good
fits according to chi-squared, whereas for fits on the interval
$[t,t_{\text{max}}],t>t_{\text{min}}$, $-\chi^{2}/2$ will approximately
increase by $(t-t_{\text{min}})/2$ indicating continued goodness-of-fit.
However, as the minimum $t$ increases in a given fit, the number of times not
included in the fit also increases: $\Delta N_{0}=t-t_{\text{min}}$. The net
effect of increasing $t>t_{\text{min}}$ is to decrease
$p(M|D)\propto\exp(-(t-t_{\text{min}})/2)$. If
$E_{Q}(\vec{p})-M_{\pi_{5}}>1/4$ we expect that the uncertainties in model
parameters will grow faster than the model probability decreases as
$t>t_{\text{min}}$. Based on these considerations, we have found
$p_{\text{cut}}=10^{-3}$ is a reasonable choice for this analysis and we adopt
it throughout. While this analysis was nearing completion, an alternate
approach to dealing with these challenges was proposed [51]. It would be
interesting to compare these two approaches in future analyses.

In our analysis of $I=2$ $\pi_{5}\pi_{5}$ scattering [13], we implemented
Bayesian model averaging. As we had hoped, the systematic uncertainties for
$\pi_{5}$-related observables were greatly reduced in that paper relative to
earlier paper [12]. Also, the problem with uncertainties increasing for
$t>t_{\text{min}}$ was not apparent because we were considering primarily
$\pi_{5}$-related observables. We expect this will not be the case for
$\sigma$ and $a_{0,1}$-related observables.

###  III.2 Shrinkage Estimator of Covariance

Suppose one wants to estimate from a multivariate sample a particular element
of the covariance matrix, then one usually uses the standard unbiased sample
estimator

$\Sigma_{ij}=\frac{1}{N-1}\sum_{n=1}^{N}(x_{i}^{(n)}-\bar{x}_{i})(x_{j}^{(n)}-\bar{x}_{j})$
(36)

which is derived from the maximum likelihood estimate (MLE) of covariance of a
multivariate Gaussian distribution. By the central limit theorem, as
$N\to\infty$ the standard estimator approaches the MLE for any distribution.
To estimate a full $\mathbb{R}^{K\times K}$ covariance matrix, there are
$K(K+1)/2$ independent matrix elements that must be estimated, requiring $N$
independent samples for each one. Furthermore, accurate estimation of the
covariance is crucial when using the chi-squared prescription in Eq. (29)
since the inverse of the covariance matrix is used and the consequence of
poorly-estimated small eigenmodes is amplified. Empirically, it has been found
that approximately $50K(K+1)/2$ samples are needed in lattice QCD calculations
for the standard estimator to be sufficiently accurate for chi-squared fitting
[52, 53] .

If you only care about this particular matrix element, or perhaps one more,
then this is the optimal estimator to use. However, if you want to
simultaneously estimate three or more elements of the covariance matrix, Stein
[54] proved that this was not the optimal estimator in the sense of minimizing
the combined mean square error, i.e.
$\sum_{ij}(\Sigma_{ij}-\Sigma_{ij}^{*})^{2}$ where $\Sigma_{ij}^{*}$ is the
true but unknown covariance. This was so counterintuitive at the time, it was
called Stein’s paradox.

For our purposes, Stein’s improved estimator will take the form of the linear
shrinkage estimator of covariance

$\sigma_{ij}(\lambda)=\lambda\Sigma_{ii}\delta_{ij}+(1-\lambda)\Sigma_{ij},\quad\lambda\in[0,1]$
(37)

and for a given sample ensemble, there exists some optimal $\lambda^{*}$ that
minimizes the MSE and $\lambda^{*}\to 0$ as $N\to\infty$. Since we don’t know
the true covariance $\Sigma^{*}$ we must estimate the optimal value. Based on
work by Ledoit and Wolf [55], Schäfer and Strimmer [56] gave a fairly
straightforward estimator for optimal value of $\lambda$

$\hat{\lambda}^{*}=\frac{\sum_{i\leq
j}\widehat{\operatorname{Var}}(\Sigma_{ij})}{\sum_{i\neq j}\Sigma_{ij}^{2}}$
(38)

In App. C, we show a one-pass algorithm to compute the sample estimate of
$\widehat{\operatorname{Var}}(\Sigma_{ij})$.

The shrinkage estimator of covariance has been suggested for use in lattice
quantum field theory (LQFT) applications for some time [57, 58]. Only recently
has the shrinkage estimator been actually employed for use in published
lattice QCD analyses [59, 60, 61, 62]. Recent work by Ledoit and Wolff [63,
64] have proposed an improved non-linear shrinkage estimator. Burda and Jarosz
[65] have also developed an improved shrinkage estimator and have developed an
open-source Python library called shrinkage to assist in calculations. In this
analysis, we have conservatively chosen to use linear shrinkage rather than
one of the newer alternatives.

##  IV Detailed Example of Model-Averaging Analysis on a Single Ensemble

We will discuss in detail our analysis of the $96^{3}\times 192$, $\beta=4.8$,
$m=0.00125$ ensemble which is the larger volume companion to the $64^{3}\times
128$ ensembles discussed in our previous work [12, 13]. It will also serve as
a detailed example of how we implemented our model-averaging analysis.

###  IV.1 Data Selection

In order to compare models fit to different data subsets, we need to first
identify the maximal data set $T$ which could be considered for any model.
Although our staggered meson two-point correlation function data is computed
from $t=0$ to $t=N_{t}-1=191$, the data are first symmetrized:
$(C(t)+C(N_{t}-t))/2\to C(t)$ and now the largest possible data set is from
$t=0$ to $t=N_{t}/2=96$. As already mentioned [50], the signal-to-noise
decreases exponentially at large times, so for most correlators, particularly
at non-zero momenta $\vec{p}\neq 0$, there is insufficient signal to
reasonably include those data points in the analysis particularly since this
will exacerbate the problem of reliable covariance estimation. We will not use
data for $t=0,1$ given the difficulties of interpreting a staggered
correlation function separated by one unit in time in terms of a transfer
matrix [37]. We compute the jackknife ratio
$C_{Q}(\vec{p},t)/C_{Q}(\vec{p},1)$ and choose a minimum value for this
quantity for each state $Q$ where there is still good signal-to-noise for all
$\vec{p}$. This defines $t_{\text{max}}$ for each $Q$ and $\vec{p}$.

Figure 1: Two-point correlation functions for $\pi_{5}$ and $a_{0,1}$ mesons.
In the right panel, data points below the horizontal line at $3\times 10^{-6}$
were not included in any fits.

Fig. 1 shows examples of our procedure. On the left, for the $\pi_{5}$ meson,
we see good signal for all momenta to the middle of the lattice and we also
see nice straight lines on the log plot indicating clear signal of a single
decaying exponential. On the right, for the $a_{0,1}$ meson, the situation is
somewhat different. There does seem to be pretty good signal to the middle of
the lattice, but the nature of the signal changes at large times, with an
apparent change of slope and an oscillating signal becoming dominant. We use a
rough model to guide our choice of where to draw a horizontal line based on
the dispersion relation $E_{Q}^{2}=M_{Q}^{2}+p^{2}$ and assuming that a single
exponential dominates the correlation function at times $t_{c}$ where it
crosses the line

$\frac{C_{Q}(\vec{p},t_{c})}{C_{Q}(\vec{p},1)}=e^{-\sqrt{M_{Q}^{2}+p^{2}}(t_{c}-1)}=\text{const}\quad\implies\quad
t_{c}(\vec{p})\propto\frac{1}{\sqrt{M_{\text{eff}}^{2}+p^{2}}}$ (39)

We compare the computed values to this model and we see good agreement along
the shown cut line. If we lower the cut line, the observed values deviate from
the prediction, particularly for $p^{2}=4$, so we conservatively set the cut
line at $3\times 10^{-6}$. Note, this model will not work well as $t_{c}\to
N_{t}/2$ since it does not include the additional contribution due to periodic
boundary conditions which becomes important in that region. A modified
expression involving hyperbolic cosines can be derived but we didn’t need it
here.

Figure 2: Two-point correlation functions for $\sigma$ meson. Left panel shows
the unsubtracted $\vec{p}=0$ correlator. Center panel shows $\vec{p}=0$
subtracted correlator. Right panel data points below horizontal line for
$\vec{p}\neq 0$ not included in fits. The data from the central panel is
included on the right by shifting upwards by a sufficiently large constant
$\overline{C}(\vec{p}=0,t)+0.045$ so that the result is positive and can be
displayed on a log plot. The shifted data cannot be used in the data selection
analysis.

The situation for the $\sigma$ meson correlator is more complicated. In Fig. 2
on the left is the unsubtracted correlator $C_{\sigma}(\vec{p}=0,t)$. It
should be clear that just subtracting some constant value around
$c_{0}=262.08\cdots$ in an uncorrelated way, following Eq. (8), would be
unsatisfactory because the signal-to-noise would fall below one in a few time
units. The center panel shows $\overline{C}_{\sigma}(\vec{p}=0,t)$ and,
following Eq. (14), the previously large positive constant has been replaced
with a three orders of magnitude smaller negative constant and greatly
enhanced signal-to-noise. However, we still need to figure out at what time
$t_{c}$ the signal-to-noise of the exponentially decaying part of the
correlator falls below an acceptable level. We cannot judge this from the
central panel since the large time behavior is dominated by the integral of
the correlation function. Instead, we compute the ratio
$C_{\sigma}(\vec{p},t)/C_{\sigma}(\vec{p},1)$ for $\vec{p}\neq 0$ and rely on
our crude model Eq. (39) to extrapolate to $\vec{p}=0$, shown in right panel.
The results in the data selection procedure are summarized in Tab. 2.

| $\vec{p}=(0,0,0)$ | $\vec{p}=(1,0,0)$ | $\vec{p}=(1,1,0)$ | $\vec{p}=(1,1,1)$ | $\vec{p}=(2,0,0)$
---|---|---|---|---|---
$\pi_{5}$ | [2,96] | [2,96] | [2,96] | [2,96] | [2,96]
$a_{0,1}$ | [2,70] | [2,63] | [2,58] | [2,53] | [2,50]
$\sigma$ | [2,52] | [2,41] | [2,34] | [2,30] | [2,27]
Table 2:  Summary of maximum allowed time ranges for fitting in model
averaging procedure for the $96^{3}\times 192$, $\beta=4.8$, $m=0.00125$
ensemble.

###  IV.2 Model Averaging

As previously discussed, this analysis will use model averaging [25]. In Fig.
3, we show how varying the fitting range $t\in[t_{\text{min}},t_{\text{max}}]$
affects the relative model probabilities $p(M|D)$. We focus on the $\vec{p}=0$
mesons since those states are most susceptible to the presence of
$t$-invariant constant contribution to the correlation function. This is true
even in the case of the $a_{0,1}$ meson where the expected constant
contribution should vanish in the infinite statistics limit. The $\pi_{5}$
meson is much less affected by any such constant as can be seen by the
preference for model A fits in the model averaging.

Figure 3:  Relative model probabilities for the $\vec{p}=0$ $\pi_{5}$ $\sigma$
and $a_{0,1}$ mesons. The different models are labeled by a letter A, B, C and
integers $(n_{\text{max}},j_{\text{max}})$, the number of non-oscillating and
oscillating states, as described in Sec. II.2. The range of time values in
each fit $[t_{\text{min}},t_{\text{max}}]$ are shown in the figures. The
uppermost curves correspond to $t_{\text{max}}$ as the maximum value in Tab. 2
and the lower curves correspond to decreasing $t_{\text{max}}$ by one.

###  IV.3 Dispersive Analysis

Once the model parameters and their errors have been computed for each
correlation function computed on a given volume, at a given fermion mass, and
at a given spatial momentum $\vec{p}$, the results from various momenta can be
used to constrain the values of the parameters in the rest frame using the
dispersion relations outlined in Eqs. (15) through (17) for the rest mass
$M_{Q}$ and Eqs. (18) through (21) for the decay constants $\widehat{F}_{Q}$.
Parameter estimation is done using least-squared fitting with possible finite
lattice spacing corrections included in even powers of $\widehat{p}^{2}$ or
$(ap)^{2}$, as appropriate. Since the number of lattice correction terms
needed is unknown a priori we use model averaging to average over the
different model choices.

Figure 4: Momentum dependence of the energy $\widehat{E}_{\pi_{5}}$ and decay
constant $\widehat{F}_{\pi_{5}}$. Fits to polynomials in $\widehat{p}^{2}$ up
to quadratic order are shown.

Figure 5: Momentum dependence of the energy $\widehat{E}_{\sigma}$ and decay
constant $\widehat{F}_{S}$. Fits to polynomials in $\widehat{p}^{2}$ up to
quadratic order are shown.

Figure 6: Momentum dependence of the energy $\widehat{E}_{a_{0,1}}$ and decay
constant $\widehat{F}_{a_{0,1}}$.

This procedure is probably of marginal benefit for the $\pi_{5}$ rest mass and
decay constant since those quantities are already very accurately determined
directly in the $\vec{p}=0$ frame and the other momentum frames do not add
significant additional information, as shown in Fig. 4. However, these fits
also show how the momentum-dependence is consistent with the expected
dispersion relations up to small lattice artifacts.

For the isosinglet scalar $\sigma$ in Fig. 5 and isotriplet scalar $a_{0,1}$
in Fig. 6, we see similar consistency with the expected dispersion relations.
Now, the information from the non-zero momentum frames provides additional
significant constraints on the rest mass and decay constant resulting in
overall smaller uncertainties than if only the $\vec{p}=0$ results alone were
used. This is particularly important for the $\sigma$ channel where the
correlation function in the $\vec{p}=0$ frame has a difficult to subtract
constant which is not present in non-zero momentum frames.

##  V Infinite volume extrapolation

We repeat the steps described in detail for one ensemble in Sec. IV for all
ensembles in this study. We would like to compare the results of our
calculations with various models but those models usually only apply to the
system in an infinite volume. We will extrapolate our data to the infinite
volume limit using the model described in Sec. II.5. At each volume and
fermion mass we compute the quantity

$\xi(m_{q},L)\equiv\frac{M_{\pi_{5}}^{2}}{(4\pi\widehat{F}_{\pi_{5}})^{2}}\sum_{n=1}^{8}\frac{4\
\kappa(n)}{\sqrt{n}\ M_{\pi_{5}}L}K_{1}(\sqrt{n}\ M_{\pi_{5}}L)$ (40)

where the $m_{q}$ dependence is implicit in the relevant infinite volume
quantities $M_{\pi_{5}}$, $\widehat{F}_{\pi_{5}}$. With this computed
quantity, the analysis becomes a simple linear fit.

If we focus just on $M_{\pi_{5}}$ and $\widehat{F}_{\pi_{5}}$, we know in
chiral perturbation theory, the quantities $\alpha_{\pi_{5}}$ and
$\beta_{\pi_{5}}$ defined in Eqs. (22) and (23) appear at a specific order in
the chiral expansion and have no implicit fermion mass dependence. We use the
same finite volume model for other masses and decay constants and we will
similarly assume the parameters $\alpha_{Q}$, $\beta_{Q}$ are mass-independent
as a model choice. This means that $\alpha_{Q}$, $\beta_{Q}$ are determined by
a simultaneous fit to the data at all fermion masses and volumes.

The choice of the expansion parameter $\xi(m_{q},L)$ being defined in terms of
infinite-volume quantities might pose a chicken-and-egg problem when
attempting to extrapolate $\pi_{5}$ data since the infinite volume values are
not known a priori. In this case, we start by using the values on the largest
volume and then iterate a few times and the result converges quickly.

An earlier version of the finite volume extrapolation for $M_{\pi_{5}}$ and
$F_{\pi_{5}}$ were published previously [13] where it was observed to be a
relatively minor correction on our volumes. Our current results are consistent
with them, so we focus here on $\sigma$ channel. The fit of
$M_{\sigma}(m_{q},L)$ is shown in Fig. 7 and the fit of
$\widehat{F}_{S}(m_{q},L)$ is shown in Fig. 8. $\alpha_{Q},\beta_{Q}$ for
various channels studied in this work are summarized in Table 3. Both from the
figures and from the uncertainties on $\alpha_{\sigma}$ and $\beta_{\sigma}$
in the table, it is clear that the uncertainties in our $\sigma$ meson
observables are still too large to reliably extract the sign and magnitude of
these finite volume corrections. We hope to return to this issue in a future
publication.

Studying the other parameters in Table 3 reveals relationships between
parameters which are generated by the strong dynamics and which are
qualitatively similar to QCD. First,
$\text{sign}(\alpha_{Q})=-\text{sign}(\beta_{Q})$ is a well-known feature in
QCD. Second, the fact that
$\text{sign}(\alpha_{\pi_{5}})=-\text{sign}(\alpha_{a_{0,1}})$ is also
observed in earlier studies [3] and was previously misinterpreted as an
indication of “parity doubling” in near-conformal gauge theories because
finite volume effects would push the masses and decay constants of parity
partners $\pi_{5}$ and $a_{0,1}$ towards degeneracy. We also note that without
a proper infinite-volume extrapolation, if the mass of the $a_{0,1}$ meson
were observed to be stable but just below decay threshold, one could wonder
whether the state might become unstable in a larger volume. In our
calculations, the $a_{0,1}$ meson remains stable even after infinite volume
extrapolation as can be seen in Tab. 4.

Figure 7: Infinite volume extrapolation of the rest masses
$M_{\sigma}(m_{q})$. $\alpha_{\sigma}=1.9\pm 16.7$ and
$\chi^{2}/\text{dof}=2.05$ with 4 dof.

Figure 8: Infinite volume extrapolation of the decay constants $\widehat{F}_{S}(m_{q})$. $\beta_{\sigma}=-1.9\pm 32.0$ and $\chi^{2}/\text{dof}=0.15$ with 4 dof. | $\alpha_{Q}$ | $\chi^{2}/\text{dof}$ | $\beta_{Q}$ | $\chi^{2}/\text{dof}$
---|---|---|---|---
$\pi_{5}$ | 6.53(29) | 3.12 | -9.3(1.3) | 0.65
$\sigma$ | 2(17) | 2.05 | -2(32) | 0.15
$a_{0,1}$ | -27.4(4.6) | 2.88 | 11(11) | 0.29
Table 3: Summary of finite volume corrections $\alpha_{Q},\beta_{Q}$. All fits have 4 dof. Multiply these parameters by $3/(2\sqrt{2\pi^{3}})\approx 0.19$ to compare with [13]. $m_{q}$ | 0.00125 | 0.00222 | 0.00500 | 0.00750 | 0.00889
---|---|---|---|---|---
$M_{\pi_{5}}$ | 0.081082(32) | 0.10870(12) | 0.165691(73) | 0.205711(33) | 0.22534(13)
$\widehat{F}_{\pi_{5}}$ | 0.021677(40) | 0.02794(12) | 0.03982(10) | 0.048314(66) | 0.05262(15)
$M_{\sigma}$ | 0.1174(44) | 0.1545(79) | 0.2408(77) | 0.2744(71) | 0.301(11)
$\widehat{F}_{S}$ | 0.0254(17) | 0.0361(37) | 0.0536(32) | 0.0711(36) | 0.0787(59)
$M_{a_{0,1}}$ | 0.1536(10) | 0.2070(53) | 0.3119(28) | 0.3773(18) | 0.4193(32)
$\widehat{F}_{a_{0,1}}$ | 0.00691(14) | 0.00944(43) | 0.1480(27) | 0.01829(24) | 0.02047(42)
$M_{\rho_{i}}$ [12] | 0.1709(65) | 0.2197(37) | 0.3024(63) | 0.36962(77) | 0.4093(21)
Table 4: Final summary of infinite volume ground state rest masses and decay
constants in lattice units. Only statistical uncertainties are shown. Data for
$M_{\rho_{i}}$ copied from [12] for convenience. See Table 9 for results with
systematic uncertainties included. Supplementary results for all fit
parameters are available [44].

##  VI Systematic Error Analysis

In our previous $I=2,\pi\pi$ scattering paper [13], we made a crude estimate
of the relative systematic errors affecting our statistical determinations of
the $\pi_{5}$ meson mass $M_{\pi_{5}}(m_{q})$ and decay constant
$F_{\pi_{5}}(m_{q})$. Our current statistical-only estimate of uncertainties
for quantities like $M_{\pi_{5}}(m_{q})$ and $\widehat{F}_{\pi_{5}}(m_{q})$ as
shown in Tab. 4 are likely underestimates due to various factors: a small
number of independent samples; various modeling choices regarding dispersion
relations and finite volume effects; data quality cuts and model probability
cuts in the model averaging procedure; plus the interplay between the amount
of independent data and choices made in the rest of the analysis through the
reliability of the shrinkage estimator of covariance.

We would like to estimate how large these effects might be in terms of a
single relative systematic error parameter $\rho$ across all the ensembles. We
will estimate $\rho$ using a number of different observables and then combine
those estimates to get an average value for $\rho$. For example, if
$\sigma_{M}$ is the statistical-only estimate of the uncertainty of a given
mass $M$, we would like to estimate a relative systematic uncertainty
$\rho_{M}$ such that the total uncertainty is

$M(m_{q})\pm\sqrt{\sigma_{M(m_{q})}^{2}+M(m_{q})^{2}\rho_{M}^{2}}$ (41)

We assume that the systematic effect is similar across all the different
ensembles labeled by different fermion masses $m_{q}$ so that the parameter
$\rho_{M}$ doesn’t depend on $m_{q}$.

To estimate $\rho_{M}$, we don’t want to assume any explicit functional
dependence for $M(m_{q})$, in particular, that we would expect to be valid for
small $m_{q}$ including as $m_{q}\to 0$. Instead, we imagine that whatever the
correct function, it is relatively smooth and slowly-varying and can be
approximated by a Taylor series expansion around the midpoint $m_{0}=0.00507$
of our range of $m_{q}$ and $|m_{q}-m_{0}|\leq\Delta_{m}=0.00382$. We can fit
the data to a polynomial

$M(m_{q})\approx\sum_{n=0}^{n_{\text{max}}}a_{n}\left(m_{q}-m_{0}\right)^{n}$
(42)

Given that we have only five different $m_{q}$ values, we will compare the
$\chi^{2}$ and AIC values for $n_{\text{max}}$=2 and $n_{\text{max}}$=3 and
use those comparisons to estimate $\rho_{M}$. We will also use the ratio test
to check for convergence of the series on $m_{0}\pm\Delta_{m}$

$\frac{|a_{n+1}\Delta^{n+1}|}{|a_{n}\Delta^{n}|}<1,\quad\forall n$ (43)

Actually the ratio test only requires the ratio $<1$ as $n\to\infty$ for
convergence, but we will assume convergence if its true term-by-term up to the
largest $n$ we can fit. For this analysis, we will use the data in Tab. 4.

###  VI.1 Fits using statistical-only data

In this section, in Tab. 5 we show fits of Eq. (42) to the statistical-only
data from Tab. 4 for $n_{\text{max}}$=2, 3. We then test for convergence by
computing the ratios in Eq. (43) and collect the results in Tab. 6.

obs. | $n_{\text{max}}$ | $\chi^{2}$ | $\log p(M|D)$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$
---|---|---|---|---|---|---|---
$M_{\pi_{5}}$ | 2 | 1052.0 | -528.0 | 0.166768(60) | 18.563(11) | -1000.9(6.4) |
| 3 | 159.5 | -83.8 | 0.167345(63) | 17.472(38) | -950.4(6.6) | 100500(3400)
$\widehat{F}_{\pi_{5}}$ | 2 | 76.0 | -41.0 | 0.040144(82) | 3.967(15) | -222.3(8.4) |
| 3 | 14.3 | -11.2 | 0.040275(84) | 3.581(51) | -207.4(8.6) | 33600(4200)
$M_{\sigma}$ | 2 | 1.6 | -3.8 | 0.2382(64) | 22.9(1.3) | -2250(640) |
| 3 | 0.8 | -4.4 | 0.2384(64) | 19.5(4.1) | -2122(657) | 29(33) $\times 10^{4}$
$\widehat{F}_{S}$ | 2 | 0.6 | -3.3 | 0.0551(28) | 6.93(62) | -200(290) |
| 3 | 0.4 | -4.2 | 0.0550(28) | 6.1(2.0) | -160(300) | 6(16) $\times 10^{4}$
$M_{a_{0,1}}$ | 2 | 14.7 | -10.4 | 0.3097(22) | 33.82(36) | -1820(210) |
| 3 | 0.0 | -4.0 | 0.3136(25) | 28.0(1.5) | -1870(210) | 46(12) $\times 10^{4}$
$\widehat{F}_{a_{0,1}}$ | 2 | 1.9 | -3.9 | 0.01483(23) | 1.726(46) | -89(23) |
| 3 | 0.0 | -4.0 | 0.01489(23) | 1.50(17) | -82(23) | 19(14) $\times 10^{3}$
Table 5:  Basic fits using statistical-only data from Tab. 4 to model function in Eq. (42). Bolded entries indicate observables where model probabilities are higher for $n_{\text{max}}=2$ than $n_{\text{max}}=3$. obs. | $n_{\text{max}}$ | $\Delta_{m}\left|\frac{a_{1}}{a_{0}}\right|$ | $\Delta_{m}\left|\frac{a_{2}}{a_{1}}\right|$ | $\Delta_{m}\left|\frac{a_{3}}{a_{2}}\right|$
---|---|---|---|---
$M_{\pi_{5}}$ | 2 | 0.42521(36) | 0.2060(14) |
| 3 | 0.39884(95) | 0.2078(15) | 0.404(15)
$\widehat{F}_{\pi_{5}}$ | 2 | 0.3775(19) | 0.2141(86) |
| 3 | 0.3397(56) | 0.2212(97) | 0.618(88)
$M_{\sigma}$ | 2 | 0.368(24) | 0.37(12) |
| 3 | 0.312(68) | 0.42(15) | 0.52(66)
$\widehat{F}_{S}$ | 2 | 0.481(54) | 0.11(16) |
| 3 | 0.43(14) | 0.10(19) | 1.6(5.9)
$M_{a_{0,1}}$ | 2 | 0.4171(64) | 0.206(26) |
| 3 | 0.342(20) | 0.254(35) | 0.94(26)
$\widehat{F}_{a_{0,1}}$ | 2 | 0.445(15) | 0.197(54) |
| 3 | 0.358(46) | 0.211(64) | 0.86(72)
Table 6:  Ratios for convergence testing of fits in Tab. 5. The fit parameter covariance matrix (not shown) was used to compute these uncertainties. Bolded entries indicate observables where the convergence test may fail due to large values or uncertainties. obs. | $n_{\text{max}}$ | $f^{\prime}(m_{q})=0$ | $f^{\prime\prime}(m_{q})=0$
---|---|---|---
$M_{\pi_{5}}$ | 2 | 0.013 |
| 3 | $0.0070\pm 0.0069i$ | 0.0070
$\widehat{F}_{\pi_{5}}$ | 2 | 0.013 |
| 3 | $0.0059\pm 0.0056i$ | 0.0059
$M_{\sigma}$ | 2 | 0.0089 |
| 3 | $0.0062\pm 0.0040i$ | 0.0062
$\widehat{F}_{S}$ | 2 | 0.021 |
| 3 | $0.0046\pm 0.0054i$ | 0.0046
$M_{a_{0,1}}$ | 2 | 0.0131 |
| 3 | $0.0052\pm 0.0043i$ | 0.0052
$\widehat{F}_{a_{0,1}}$ | 2 | 0.0135 |
| 3 | $0.0053\pm 0.0050i$ | 0.0053
Table 7: Zeroes of the derivatives of fits in Tab. 5.

If we first look at the model probabilities, we see when the fit is highly-
constrained, indicated by large $\chi^{2}$ values, then the fit with
$n_{\text{max}}=3$ is preferred relative to $n_{\text{max}}=2$. This is the
expected behavior since adding extra fit parameters in a highly-constrained
fit usually reduces the $\chi^{2}$ by a sufficient amount to increase the
model probability. However, if the fit is poorly-constrained, indicated by a
small $\chi^{2}$, adding extra parameters may not increase the model
probability. Observables where this occurs are highlighted in Table 5 and
those observables are probably too noisy to help constrain the systematic
error parameter $\rho$.

Looking at the convergence test in Table 6, again we highlight examples where
data were too noisy to pass the test with confidence. Again, we will not use
those observables to help constrain $\rho$. Note also the strong overlap in
the lists of rejected observables from both tables. Finally, we don’t expect
that the functions will have extremal points in the region where it
approximates the data. The zeroes of the derivatives are shown in Tab. 7.

###  VI.2 Estimating Relative Systematic Error

To estimate the relative systematic error parameter $\rho$, from Eq. (41) as
$\rho$ increases the error bars will increase and the corresponding $\chi^{2}$
will decrease. What value of $\chi^{2}$ should we choose to determine $\rho$?
A priori, two interesting values come to mind: (I) the mean value of the chi-
squared distribution for $k$ degrees of freedom, i.e.
$\chi_{k}^{2}(\rho^{(I)})=k$; (II) the value of $\chi^{2}_{k}$ such that one
expects 68% of the time a random sample of the chi-squared distribution should
be less than or equal to that value, i.e. $\chi_{1}^{2}(\rho)=1$,
$\chi_{2}^{2}(\rho^{(II)})\approx 2.3$.

A posteriori, we noticed that from a model averaging perspective, the
4-parameter cubic polynomial fit has the higher model probability at $\rho$=0
in cases where the statistical error is small compared to the expected
systematic error. In the $\rho\to\infty$ limit, $\chi^{2}\to 0$ and the most
likely model is the one with the smallest $n_{\text{max}}$. As $\rho$
increases, there is a point where the quadratic and cubic polynomial fits have
equal probability. We define

$\text{AIC}(\rho,n_{\text{max}})=\frac{1}{2}\chi_{5-n_{\text{max}}}^{2}+n_{\text{max}}+1$
(44)

and choose a third interesting value of $\rho$: (III)
$\text{AIC}(\rho^{(III)},2)=\text{AIC}(\rho^{(III)},3)$. Note this does not
always have a solution, particularly if the quadratic fit has a lower AIC at
$\rho=0$. A posteriori we can rationalize this choice as the point where the
quadratic and cubic descriptions of the data are equally good (or bad) from an
information-theoretic perspective.

obs. | $n_{\text{max}}$ | $\rho^{(I)}$ | $\rho^{(II)}$ | $\rho^{(III)}$
---|---|---|---|---
$M_{\pi}$ | 2 | 0.0169 | 0.0158 | 0.0162
| 3 | 0.0068 | 0.0068 | 0.0162
$F_{\pi}$ | 2 | 0.0210 | 0.0196 | 0.0197
| 3 | 0.0098 | 0.0098 | 0.0197
Table 8:  Various estimates of the systematic error parameter $\rho$ as
determined by methods described in the text.

###  VI.3 Summary of Systematic Error Analysis

From Tab. 8 we can see there are eight unique $\rho$ values from approximately
0.007 to 0.021. Rather than pick just one, we consider a few summary
statistics: the arithmetic mean $\rho_{\text{a}}=0.0157$, the median
$\rho_{\text{m}}=0.0165$, or the geometric mean $\rho_{\text{g}}=0.0148$. All
give relatively similar values close to the central grouping. We make a
conservative choice and choose the largest of the three
$\rho=\rho_{\text{m}}=0.0165$. If we compare this estimate to the previous
rough guess of 0.01 quoted in [13], it is nice to see they are not too
different and that 0.01 falls within the range of estimated values. A final
summary of our results with the relative systematic error included is given in
Table 9.

$m_{q}$ | 0.00125 | 0.00222 | 0.00500 | 0.00750 | 0.00889
---|---|---|---|---|---
$M_{\pi_{5}}$ | 0.0811(13) | 0.1087(18) | 0.1657(27) | 0.2057(34) | 0.2253(37)
$\widehat{F}_{\pi_{5}}$ | 0.02168(36) | 0.02794(47) | 0.03982(66) | 0.04831(80) | 0.05262(88)
$M_{\sigma}$ | 0.1174(48) | 0.1545(83) | 0.2408(87) | 0.2744(85) | 0.301(12)
$\widehat{F}_{S}$ | 0.0254(17) | 0.0361(37) | 0.0536(33) | 0.0711(38) | 0.0787(60)
$M_{a_{0,1}}$ | 0.1536(27) | 0.2070(63) | 0.3119(58) | 0.3773(65) | 0.4193(76)
$\widehat{F}_{a_{0,1}}$ | 0.00691(17) | 0.00944(45) | 0.1480(36) | 0.01829(39) | 0.02047(54)
$M_{\rho_{i}}$ [12] | 0.1709(71) | 0.2197(52) | 0.3024(80) | 0.3696(61) | 0.4093(71)
Table 9: Final summary of infinite volume ground state rest masses and decay
constants with relative systematic error of $\rho=0.0165$ included following
Eq. (41). Data for $M_{\rho_{i}}$ derived from [12] for convenience. Results
with only statistical errors in Table 4.

In Fig. 9 we compare the previously computed results for $M_{\sigma}$ with
combined statistical and systematic errors as described in [12] with the new
results of Tab. 9 for $M_{\sigma}$ and $M_{\pi_{5}}$. The values for
$\sqrt{8t_{0}}/a$ are taken from Tab. I of [12] and the plot style is similar
to Fig. 10 of [12]. The conclusion we draw from this comparison was that the
previous analysis method for computing $M_{\sigma}$ led to systematically
lower mass values and that the method used previously to estimate systematic
errors was sufficiently conservative as to cover the the downward shift of the
result.

Figure 9:  A comparison of $M_{\sigma}$ computed previously [12] with the
results from Tab. 9.

##  VII Comments on Chiral Expansions

While the SU(3) $N_{f}=8$ theory with massive Dirac fermions is a potentially
interesting theory on its own, being a possible candidate for composite dark
matter [66], the theory closer to the chiral limit $m_{q}\to 0$ might also be
relevant for composite Higgs phenomenology. As stated in Sec. I, our results
alone are not sufficient to establish with certainty whether the massless
$N_{f}=8$ theory is inside or outside the conformal window. But, the low
energy content of the two scenarios is quite different: in one case a very
strongly coupled conformal field theory and in the other case massless Nambu-
Goldstone bosons and possible light flavor-singlet scalar resonance with a
mass of order $F_{\pi}$. Specific models will appear quite different depending
on the scenario, and when fitted to our data, those models will, in general,
be an expansion in some small parameter which vanishes in the chiral limit. We
will discuss different specific extrapolations in detail in a companion paper
[34].

Here we note that based on one’s a priori expectation for the nature of the
low-energy theory in the chiral limit, the choice of expansion parameter can
lead to very different presentations of the data. One could naively plot
results vs. the fermion mass $m_{q}$, or some power of the fermion mass
$m_{q}^{1/(1+\gamma^{*})}$, $0\leq\gamma^{*}\leq 2$ motivated by assuming
conformal symmetry in the chiral limit, or $\chi\equiv
M_{\pi_{5}}^{2}/(4\pi\widehat{F}_{\pi_{5}})^{2}$ by assuming spontaneous
chiral symmetry breaking in the chiral limit. In theory like SU(3) $N_{f}=2$
(QCD), these choices often don’t make any appreciable difference in the
presentation of the data. But in this theory, if the chiral limit is
conformal, the expansion parameter $\chi$ doesn’t vanish as $m_{q}\to 0$.
Visually, we can see the difference in Fig. 10. Since the value of
$\gamma^{*}$ is a dynamical parameter that can only be determined through a
careful extrapolation, we plot three representative values that cover weakly
and strongly-coupled CFTs plus an intermediate value. Regardless of which
value of $\gamma^{*}$ is chosen, $\chi$ varies significantly over the range of
$m_{q}$ with a fair degree of curvature which makes it difficult to estimate
how small $m_{q}$ must be before the constant term in $\chi$ dominates the
leading $m_{q}$-dependent term. Of course, if $\chi$ vanishes in the chiral
limit, then the constant term will never dominate. This suggests it will be
difficult to distinguish with much certainty given our current results whether
or not $\chi$ vanishes in the chiral limit. Calculations at smaller fermion
masses are needed.

Figure 10:  Chiral expansion parameter
$\chi=M_{\pi_{5}}^{2}/(4\pi\widehat{F}_{\pi_{5}})^{2}$ vs. other chiral
expansion parameters $m_{q}^{1/(1+\gamma^{*})}$. If the theory is conformal
$\chi$ should be non-zero in the chiral limit. If the theory is spontaneously
broken, $\chi$ should be zero in the chiral limit.

We can now present the results for the spectrum in two different ways. In Fig.
11, on the left is a presentation appropriate when assuming the theory is
conformal in the chiral limit with a mass anomalous dimension
$\gamma^{*}\approx 1$. In units of the lattice spacing $a$, the masses of all
the hadrons are expected to extrapolate to zero since any non-zero hadron mass
would break conformal symmetry. On the right is a presentation assuming the
chiral symmetry is spontaneously broken in the chiral limit and the relevant
scale of chiral symmetry breaking is set by $4\pi\widehat{F}_{\pi_{5}}$. All
the hadron masses except the pion should be non-zero in the chiral limit.
Plotted this way, the pion is shown as a simple curve since
$M_{\pi_{5}}/4\pi\widehat{F}_{\pi_{5}}=\sqrt{\chi}$. This also makes it easy
to display the energy threshold as a dotted line for decay to two pions. In
the current data set, both the flavor-singlet and non-singlet scalar mesons
appear to be unable to decay to two pions.

Figure 11:  Two different presentations of the spectrum from Tab. 9. On the
left, in units of the lattice spacing $a$ vs. a chiral expansion parameter
assuming conformal symmetry and $\gamma^{*}\approx 1$. On the right, in units
of the chiral breaking scale $4\pi\widehat{F}_{\pi_{5}}$ vs. a chiral
expansion parameter assuming spontaneous chiral symmetry breaking. The dotted
line on the right indicates the energy threshold for decays to two pions.

Focusing solely on the data presented in this section, it is still far from
clear whether or not the chiral parameter $\chi$ vanishes in the chiral limit.
On the other hand, $\chi$ varies significantly with the fermion mass which
also suggests we are far from the hyperscaling limit where $\chi$ should be a
non-zero constant. Recent numerical studies with improved gauge action [10]
suggest that the SU(3) $N_{f}=8$ system could be at the opening of the
conformal window, or at least very close to it. There are indications of an
infrared fixed point at much stronger couplings than what is probed by our
data in this paper. This is so even if $N_{f}=8$ is below the conformal
window. Therefore, corrections to scaling in the gauge coupling could be
significant. This can explain our inability to distinguish between the
conformal and chirally broken scenarios.

##  VIII Gell-Mann Oakes Renner (GMOR) Ratio Results

As discussed in Sec. II.6, numerical studies of the GMOR ratio can shed light
on the low-energy behavior of the $N_{f}=8$ theory by measuring how much the
ground-state pion pole contributes the pseudoscalar two-point correlation
function. A value close to unity indicates pion pole-dominance. Tab. 10 shows
the computed values for the chiral condensate and the integrated pseudoscalar
correlation function. Although computed by two different techniques, the
results agree extremely well with Eq. (25). Using the largest volume data at
each fermion mass for the condensate and the statistical-only data in Tab. 4,
we compute the GMOR ratio in Eq. (28) and propagate the statistical-only
errors. We then apply the relative systematic error correction estimated in
Sec. VI. The results are shown in the rightmost column of Tab. 10. The lowest
pole does not fully dominate the pion correlation function in our fermion mass
range as the result is larger than one. If we assume a mass-deformed CFT is
the correct low energy description, then the lowest pion pole will never
dominate the pseudoscalar correlation function at any fermion mass as it is
not a pseudo-Nambu-Goldstone boson, as in Eq. (28).

In the case the $N_{f}=8$ theory is outside the conformal window, one of the
poles contributing to the pion correlation function would have a pole position
at $M_{\pi_{5}}+M_{\sigma}$ and a residue proportional to
$g^{2}_{\pi\pi\sigma}$ which we will assume is
$\mathcal{O}(\widehat{F}_{\pi_{5}}^{2})$. For the leading pion pole to
dominate, $M_{\pi_{5}}^{2}/(M_{\pi_{5}}+M_{\sigma})^{2}\ll 1$. In this work,
the ratio varies from 0.167(8) — 0.183(9), which in this scenario is
interpreted as not small enough to ensure pion pole dominance. A direct
calculation of the coupling $g_{\pi\pi\sigma}$ and/or further calculations at
lighter fermion masses should shed light on this issue.

We did not perform an infinite volume extrapolation of the condensate data in
Tab. 10 similar to the ones described in Sec. II.5. The systematic effect of
this correction might be significant on the scale of the uncertainties shown
for the GMOR ratio. But, the effect is unlikely to be significant relative to
the deviation of the ratio from unity. In the future, if a detailed model is
to be fit to these data, the modeller should consider including these
neglected corrections.

$m_{q}$ | $L$ | $\left\langle\overline{\chi}\chi\right\rangle$ | $\sum_{t}C_{\pi_{5}}(\vec{0},t)$ | $R_{G}(m_{q})$
---|---|---|---|---
0.00125 | 96 | 0.0121704(53) | 9.732(28) | 2.462(42)
| 64 | 0.0121220(70) | 9.641(17) |
0.00222 | 48 | 0.019918(14) | 9.039(25) | 2.397(45)
0.00500 | 48 | 0.040808(18) | 8.164(41) | 2.344(41)
| 32 | 0.040521(21) | 8.103(12) |
0.00750 | 48 | 0.058447(12) | 7.808(20) | 2.219(37)
| 32 | 0.058445(26) | 7.786(16) |
| 24 | 0.058063(73) | 7.752(32) |
0.00889 | 32 | 0.068086(28) | 7.661(37) | 2.152(38)
| 24 | 0.067814(46) | 7.615(17) |
Table 10: Values for the staggered chiral condensate
$\left\langle\overline{\chi}\chi\right\rangle$, computed using a noisy
estimator, and the integrated pseudoscalar correlation function, computed
using a point source. The reader can verify that the columns satisfy Eq. (25).
Only statistical errors are shown. The rightmost column shows the GMOR ratio
defined in Eq. (28) with errors computed as described in the text.

##  IX Discussion

In this investigation, we have made many methodological improvements with
respect to our earlier lattice study of the $N_{f}=8$ theory [12]. In
particular, we have employed two different methods for dealing with time-
independent contributions to the flavor-singlet scalar correlator, first by
using the subtraction scheme developed in Sec. II.1 and then by working with
moving frames and applying the dispersion relation described in Sec. II.3. We
were able to substantially reduce the systematic uncertainties of our fit
results using the Bayesian model averaging approach. Additionally, we used
improved “linear” shrinkage estimators for data covariance which we found were
more reliable given the amount of statistics. There was an open question in
our previous paper whether finite volume effects could be significant even
when $M_{\pi_{5}}L\gtrsim 5.3$. Now, we can see that the finite volume effects
are mild and don’t play a significant role in the final result. We find that
$M_{\sigma}/M_{\pi_{5}}$ ranges from 1.45 – 1.34 as $M_{\pi_{5}}/M_{\rho_{i}}$
increases from 0.47 – 0.55.

We computed a new observable, the scalar decay constant $\widehat{F}_{S}$,
which, as we show in a related paper [34], provides useful independent
constraints on various low-energy effective theories. We also computed the
flavor-nonsiglet scalar meson mass $M_{a_{0,1}}$ and decay constant
$\widehat{F}_{a_{0,1}}$. The proximity of the $a_{0,1}$ to the decay threshold
suggests that a careful elastic scattering analysis might be warranted in the
future if more accurate results are desired.

###### Acknowledgements.

R.C.B. and C.R. acknowledge United States Department of Energy (DOE) Award No.
DE-SC0015845. K.C. acknowledges support from the DOE through the Computational
Sciences Graduate Fellowship (DOE CSGF) through grant No. DE-SC0019323 and
also from the P.E.O. Scholar award. G.T.F. acknowledges support from DOE Award
No. DE-SC0019061. A.D.G. is supported by SNSF grant No. 200021_17576. A.H. and
E.T.N. acknowledge support by DOE Award No. DE-SC0010005. J.I. acknowledges
support from ERC grant No. 101039756. D.S. was supported by UK Research and
Innovation Future Leader Fellowship MR/S015418/1 and STFC grant ST/T000988/1.
P.V. acknowledges the support of the DOE under contract No. DE-AC52-07NA27344
(Lawrence Livermore National Laboratory, LLNL). We thank the LLNL
Multiprogrammatic and Institutional Computing program for Grand Challenge
supercomputing allocations. We also thank Argonne Leadership Computing
Facility (ALCF) for allocations through the INCITE program. ALCF is supported
by DOE contract No. DE-AC02-06CH11357. Computations for this work were carried
out in part on facilities of the USQCD Collaboration, which are funded by the
Office of Science of the DOE, and on Boston University computers at the
MGHPCC, in part funded by the National Science Foundation (award No.
OCI-1229059). This research utilized the NVIDIA GPU accelerated Summit
supercomputer at Oak Ridge Leadership Computing Facility at the Oak Ridge
National Laboratory, which is supported by the DOE Office of Science under
Contract No. DE-AC05-00OR22725.

##  Appendix A Ensembles

A summary of the ensembles used in this paper are shown in Tables 11 and 12.

Volume | Mass | Try | MDTU | Period (MDTU) | Block (MDTU) | $N_{\text{blk}}$
---|---|---|---|---|---|---
$24^{3}\times 48$ | 0.00889 | 1 | [250,25000] | 10 | 100 | 247
$32^{3}\times 64$ |  | 1 | [1040,7000] | 40 | 80 | 75
|  | 2 | [1040,7000] | 40 | 80 | 75
|  | 3 | [1040,7000] | 40 | 80 | 75
|  | 4 | [1040,7000] | 40 | 80 | 75
|  | C |  |  | 80 | 300
$24^{3}\times 48$ | 0.0075 | 1 | [350,10000] | 10 | 90 | 107
$32^{3}\times 64$ |  | 1 | [255,1395] | 10 | 100 |
|  |  | [1400,25160] | 5 | 100 | 249
$48^{3}\times 96$ |  | 1 | [250,9990] | 10 | 70 | 139
|  | 2 | [250,9990] | 10 | 70 | 139
|  | C |  |  | 70 | 278
$32^{3}\times 64$ | 0.005 | 1 | [251,29641] | 5 | 100 | 293
|  | 2 | [20011,22815] | 2 | 100 | 28
|  | 3 | [29001,31653] | 2 | 100 | 26
|  | 4 | [10001,13293] | 2 | 100 | 32
|  | C |  |  | 100 | 379
$48^{3}\times 96$ |  | 1 | [250,4200] | 10 | 50 | 79
|  | 2 | [250,3390] | 10 | 50 | 63
|  | C |  |  | 50 | 142
Table 11:  Ensembles, or Markov chains, used in this study with $0.005\leq m_{q}\leq 0.00889$. “Try” assigns a label to each Markov chain and the label “C” indicates the combined summary for all chains at a given mass and volume. “Period” indicates how often the correlation functions were computed.  Volume | Mass | Try | MDTU | Period (MDTU) | Block (MDTU) | $N_{\text{blk}}$
---|---|---|---|---|---|---
$48^{3}\times 96$ | 0.00222 | 1 | [250,11190] | 2 | 120 | 91
|  | 2 | [1000,9930] | 2 | 120 | 74
|  | 3 | [210,1450] | 10 | 120 | 10
|  | 4 | [210,1410] | 10 | 120 | 10
|  | 5 | [210,1360] | 10 | 120 | 9
|  | 6 | [210,1290] | 10 | 120 | 9
|  | 7 | [210,1350] | 10 | 120 | 9
|  | C |  |  | 120 | 212
$64^{3}\times 128$ | 0.00125 | r0 | [200,2060] | 10 | 120 | 15
|  | r1 | [200,1990] | 10 | 120 | 15
|  | r2 | [200,2010] | 10 | 120 | 15
|  | r3 | [200,2070] | 10 | 120 | 15
|  | s0 | [8436,17088] | 12 | 120 | 72
|  | s1 | [7644,17472] | 12 | 120 | 82
|  | s2 | [7212,17412] | 12 | 120 | 86
|  | C |  |  | 120 | 300
$96^{3}\times 128$ |  | 2 | [500,3144] | 2 | 80 | 34
|  | 3 | [500,3282] | 2 | 80 | 35
|  | C |  |  | 80 | 69
Table 12:  Ensembles, or Markov chains, used in this study with $0.00125\leq
m_{q}\leq 0.00222$. “Try” assigns a label to each Markov chain and the label
“C” indicates the combined summary for all chains at a given mass and volume.
“Period” indicates how often the correlation functions were computed.

##  Appendix B Normalization of model probabilities

When performing an aggressive model averaging analysis by considering a wide
range of models $\\{M\\}$ and a wide range of data subset selections $T_{1}$
for each model, the resulting set of $\log p(M|D)$ can vary by several orders
of magnitude, making it numerically challenging to perform an accurate
calculation of $\sum_{\\{M\\}}p(M|D)$. In particular, exponentiating each
$\log p(M|D)$ and then performing the sum seems like a bad idea. So, we work
directly with $\log p(M|D)$ to compute the log of the sum. Let $\ell_{n}$ be a
sorted list of the $\log p(M|D)$:
$\ell_{1}\leq\ell_{2}\leq\cdots\leq\ell_{N}$. We can construct the partial
sums recursively

$s_{1}=\ell_{1},\qquad
s_{n}=s_{n-1}+\log\left(1+e^{\ell_{n}-s_{n-1}}\right)\quad(n>1)$ (45)

The final sum over model probabilities is $\sum_{\\{M\\}}p(M|D)=\exp s_{N}$.
The key observation is that sorting the list ensures that two wildly different
numbers are not combined at any step with accompanying large loss of
precision.

##  Appendix C  Unbiased Sample Estimator for the Variance of the Covariance

Using Mathematica’s MomentConvert[] functionality, it is a few lines of code
to express the unbiased sample estimator for
$\widehat{\operatorname{Var}}(\Sigma_{ij})$ in terms of raw moments

    centMom11Est = MomentConvert[CentralMoment[{1, 1}], "SampleEstimator"];
    bias = MomentConvert[centMom11Est, {Moment, n}];
    MomentConvert[(centMom11Est - bias)^2, {Moment, n}]

The result is

$\displaystyle\frac{n^{3}}{n-1}\widehat{\operatorname{Var}}(\operatorname{Cov}(x,y))$
(46) $\displaystyle=$
$\displaystyle(n-1)\mu_{2,2}-2(n-1)\left(\mu_{2,1}\mu_{0,1}+\mu_{1,0}\mu_{1,2}\right)+\mu_{2,0}\mu_{0,2}+(n-2)\left(\mu_{2,0}\mu_{0,1}^{2}+\mu_{1,0}^{2}\mu_{0,2}\right)$
$\displaystyle-(n-2)\mu_{1,1}^{2}+2(3n-4)\mu_{1,1}\mu_{1,0}\mu_{0,1}-2(2n-3)\mu_{1,0}^{2}\mu_{0,1}^{2}$

where we use Mathematica’s convention for raw moments

$\mu_{i,j;\mathcal{S}}=\frac{1}{n}\sum_{(x,y)\in\mathcal{S}}x^{i}y^{j}$ (47)

Following Pébay [67], we would like to construct a one-pass, parallelizable
computation. To explain the notation, $\mathcal{S}$ is a set of $n$ samples
that can be partitioned into two subsets $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$ of $n_{1}$ and $n_{2}$ samples, respectively, so
$n_{1}+n_{2}=n$. The computation can be parallelized by performing
computations on the subsets and combining the results. In the special case
where $n_{1}=n-1$ and $n_{2}=1$, the results simplify and can be used as a
one-pass algorithm

$\mu_{i,j;\mathcal{S}}=\frac{n_{1}}{n}\mu_{i,j;\mathcal{S}_{1}}+\frac{n_{2}}{n}\mu_{i,j;\mathcal{S}_{2}}=\mu_{i,j;\mathcal{S}_{1}}+\frac{n_{2}}{n}\left(\mu_{i,j;\mathcal{S}_{2}}-\mu_{i,j;\mathcal{S}_{1}}\right)$
(48)

where the first form is symmetric and more useful when $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$ are of comparable size and the second form is better suited
when $\mathcal{S}_{2}$ is a single sample $(x,y)$

$\mu_{i,j;\mathcal{S}}=\mu_{i,j;\mathcal{S}_{1}}+\frac{1}{n}\left(x^{i}y^{j}-\mu_{i,j;\mathcal{S}_{1}}\right).$
(49)

## References

  * Caswell [1974] W. E. Caswell, Asymptotic Behavior of Nonabelian Gauge Theories to Two Loop Order, Phys. Rev. Lett. 33, 244 (1974).
  * Banks and Zaks [1982] T. Banks and A. Zaks, On the Phase Structure of Vector-Like Gauge Theories with Massless Fermions, Nucl. Phys. B 196, 189 (1982).
  * Fleming [2008] G. T. Fleming, Strong Interactions for the LHC, PoS LATTICE2008, 021 (2008), arXiv:0812.2035 [hep-lat] .
  * Drach [2020] V. Drach, Composite electroweak sectors on the lattice, PoS LATTICE2019, 242 (2020), arXiv:2005.01002 [hep-lat] .
  * Ishikawa  _et al._ [2015] K. I. Ishikawa, Y. Iwasaki, Y. Nakayama, and T. Yoshie, IR fixed points in $SU(3)$ gauge Theories, Phys. Lett. B 748, 289 (2015), arXiv:1503.02359 [hep-lat] .
  * Fodor  _et al._ [2018] Z. Fodor, K. Holland, J. Kuti, D. Nogradi, and C. H. Wong, Is SU(3) gauge theory with 13 massless flavors conformal?, PoS LATTICE2018, 198 (2018), arXiv:1811.05024 [hep-lat] .
  * Fodor  _et al._ [2019] Z. Fodor, K. Holland, J. Kuti, D. Nogradi, and C. H. Wong, Case studies of near-conformal $\beta$-functions, PoS LATTICE2019, 121 (2019), arXiv:1912.07653 [hep-lat] .
  * Appelquist  _et al._ [2012] T. Appelquist _et al._ (Lattice Strong Dynamics (LSD)), Approaching Conformality with Ten Flavors, (2012), arXiv:1204.6000 [hep-ph] .
  * Appelquist  _et al._ [2021a] T. Appelquist _et al._ (Lattice Strong Dynamics (LSD)), Near-conformal dynamics in a chirally broken system, Phys. Rev. D 103, 014504 (2021a), arXiv:2007.01810 [hep-ph] .
  * Hasenfratz [2022] A. Hasenfratz, Emergent strongly coupled ultraviolet fixed point in four dimensions with eight Kähler-Dirac fermions, Phys. Rev. D 106, 014513 (2022), arXiv:2204.04801 [hep-lat] .
  * Appelquist  _et al._ [2014] T. Appelquist _et al._ (LSD), Lattice simulations with eight flavors of domain wall fermions in SU(3) gauge theory, Phys. Rev. D 90, 114502 (2014), arXiv:1405.4752 [hep-lat] .
  * Appelquist  _et al._ [2019] T. Appelquist _et al._ (Lattice Strong Dynamics), Nonperturbative investigations of SU(3) gauge theory with eight dynamical flavors, Phys. Rev. D 99, 014509 (2019), arXiv:1807.08411 [hep-lat] .
  * Appelquist  _et al._ [2022a] T. Appelquist _et al._ (Lattice Strong Dynamics (LSD)), Goldstone boson scattering with a light composite scalar, Phys. Rev. D 105, 034505 (2022a), arXiv:2106.13534 [hep-ph] .
  * Kunihiro  _et al._ [2004] T. Kunihiro, S. Muroya, A. Nakamura, C. Nonaka, M. Sekiguchi, and H. Wada (SCALAR), Scalar mesons in lattice QCD, Phys. Rev. D 70, 034504 (2004), arXiv:hep-ph/0310312 .
  * Wakayama  _et al._ [2015] M. Wakayama, T. Kunihiro, S. Muroya, A. Nakamura, C. Nonaka, M. Sekiguchi, and H. Wada, Lattice QCD study of four-quark components of the isosinglet scalar mesons: Significance of disconnected diagrams, Phys. Rev. D 91, 094508 (2015), arXiv:1412.3909 [hep-lat] .
  * Howarth and Giedt [2017] D. Howarth and J. Giedt, The sigma meson from lattice QCD with two-pion interpolating operators, Int. J. Mod. Phys. C 28, 1750124 (2017), arXiv:1508.05658 [hep-lat] .
  * Briceno  _et al._ [2017] R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Isoscalar $\pi\pi$ scattering and the $\sigma$ meson resonance from QCD, Phys. Rev. Lett. 118, 022002 (2017), arXiv:1607.05900 [hep-ph] .
  * Briceno  _et al._ [2018a] R. A. Briceno, J. J. Dudek, and R. D. Young, Scattering processes and resonances from lattice QCD, Rev. Mod. Phys. 90, 025001 (2018a), arXiv:1706.06223 [hep-lat] .
  * Briceno  _et al._ [2018b] R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Isoscalar $\pi\pi,K\overline{K},\eta\eta$ scattering and the $\sigma,f_{0},f_{2}$ mesons from QCD, Phys. Rev. D 97, 054513 (2018b), arXiv:1708.06667 [hep-lat] .
  * Guo  _et al._ [2018] D. Guo, A. Alexandru, R. Molina, M. Mai, and M. Döring, Extraction of isoscalar $\pi\pi$ phase-shifts from lattice QCD, Phys. Rev. D 98, 014507 (2018), arXiv:1803.02897 [hep-lat] .
  * Rodas  _et al._ [2023a] A. Rodas, J. J. Dudek, and R. G. Edwards, The quark mass dependence of $\pi\pi$ scattering in isospin 0, 1 and 2 from lattice QCD, (2023a), arXiv:2303.10701 [hep-lat] .
  * Rodas  _et al._ [2023b] A. Rodas, J. J. Dudek, and R. G. Edwards, Constraining the quark mass dependence of the lightest resonance in QCD, (2023b), arXiv:2304.03762 [hep-lat] .
  * Kawarabayashi and Suzuki [1966] K. Kawarabayashi and M. Suzuki, Partially conserved axial vector current and the decays of vector mesons, Phys. Rev. Lett. 16, 255 (1966).
  * Riazuddin and Fayyazuddin [1966] Riazuddin and Fayyazuddin, Algebra of current components and decay widths of rho and K* mesons, Phys. Rev. 147, 1071 (1966).
  * Jay and Neil [2021] W. I. Jay and E. T. Neil, Bayesian model averaging for analysis of lattice field theory results, Phys. Rev. D 103, 114502 (2021), arXiv:2008.01069 [stat.ME] .
  * Panico and Wulzer [2016] G. Panico and A. Wulzer,  _The Composite Nambu-Goldstone Higgs_, Vol. 913 (Springer, 2016) arXiv:1506.01961 [hep-ph] .
  * Ferretti [2016] G. Ferretti, Gauge theories of Partial Compositeness: Scenarios for Run-II of the LHC, JHEP 06, 107, arXiv:1604.06467 [hep-ph] .
  * Witzel [2019] O. Witzel, Review on Composite Higgs Models, PoS LATTICE2018, 006 (2019), arXiv:1901.08216 [hep-lat] .
  * Cacciapaglia  _et al._ [2019] G. Cacciapaglia, G. Ferretti, T. Flacke, and H. Serôdio, Light scalars in composite Higgs models, Front. in Phys. 7, 22 (2019), arXiv:1902.06890 [hep-ph] .
  * Cacciapaglia  _et al._ [2020] G. Cacciapaglia, C. Pica, and F. Sannino, Fundamental Composite Dynamics: A Review, Phys. Rept. 877, 1 (2020), arXiv:2002.04914 [hep-ph] .
  * Vecchi [2017] L. Vecchi, A dangerous irrelevant UV-completion of the composite Higgs, JHEP 02, 094, arXiv:1506.00623 [hep-ph] .
  * Appelquist  _et al._ [2021b] T. Appelquist, J. Ingoldby, and M. Piai, Nearly Conformal Composite Higgs Model, Phys. Rev. Lett. 126, 191804 (2021b), arXiv:2012.09698 [hep-ph] .
  * Appelquist  _et al._ [2022b] T. Appelquist, J. Ingoldby, and M. Piai, Composite two-Higgs doublet model from dilaton effective field theory, Nucl. Phys. B 983, 115930 (2022b), arXiv:2205.03320 [hep-ph] .
  * Appelquist  _et al._ [2023] T. Appelquist _et al._ (Lattice Strong Dynamics (LSD)), Hidden Conformal Symmetry from the Lattice, (2023), arXiv:2305.03665 [hep-lat] .
  * Daniel and Sheard [1988] D. Daniel and S. N. Sheard, PERTURBATIVE CORRECTIONS TO STAGGERED FERMION LATTICE OPERATORS, Nucl. Phys. B 302, 471 (1988).
  * Engels  _et al._ [1982] J. Engels, F. Karsch, and H. Satz, Finite Size Effects in Euclidean Lattice Thermodynamics for Noninteracting Bose and Fermi Systems, Nucl. Phys. B 205, 239 (1982).
  * Kilcup and Sharpe [1987] G. Kilcup and S. R. Sharpe, A Tool Kit for Staggered Fermions, Nucl. Phys. B 283, 493 (1987).
  * Shi  _et al._ [2000] F. Shi, T. G. Steele, V. Elias, K. B. Sprague, Y. Xue, and A. H. Fariborz, Holder inequalities and isospin splitting of the quark scalar mesons, Nucl. Phys. A 671, 416 (2000), arXiv:hep-ph/9909475 .
  * Maltman [2000] K. Maltman, Scalar meson decay constants and the nature of the $a_{0}(980)$, Nucl. Phys. A 680, 171 (2000), arXiv:hep-ph/0005155 .
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2020, 083C01 (2020).
  * Aoki  _et al._ [2017] Y. Aoki _et al._ (LatKMI), Light flavor-singlet scalars and walking signals in $N_{f}=8$ QCD on the lattice, Phys. Rev. D 96, 014508 (2017), arXiv:1610.07011 [hep-lat] .
  * Golterman [2009] M. Golterman, Applications of chiral perturbation theory to lattice QCD, in  _Les Houches Summer School: Session 93: Modern perspectives in lattice QCD: Quantum field theory and high performance computing_ (2009) arXiv:0912.4042 [hep-lat] .
  * Appelquist  _et al._ [2011] T. Appelquist, G. T. Fleming, M. F. Lin, E. T. Neil, and D. A. Schaich, Lattice Simulations and Infrared Conformality, Phys. Rev. D 84, 054501 (2011), arXiv:1106.2148 [hep-lat] .
  * Fleming  _et al._ [2023] G. T. Fleming _et al._ (Lattice Strong Dynamics (LSD)), Dataset for ”light scalar meson and decay constant in SU(3) gauge theory with eight dynamical flavors” (2023), zenodo https://dx.doi.org/10.5281/zenodo.8007955.
  * Gell-Mann  _et al._ [1968] M. Gell-Mann, R. J. Oakes, and B. Renner, Behavior of current divergences under SU(3) x SU(3), Phys. Rev. 175, 2195 (1968).
  * Banks and Casher [1980] T. Banks and A. Casher, Chiral Symmetry Breaking in Confining Theories, Nucl. Phys. B 169, 103 (1980).
  * Leutwyler and Smilga [1992] H. Leutwyler and A. V. Smilga, Spectrum of Dirac operator and role of winding number in QCD, Phys. Rev. D 46, 5607 (1992).
  * Patella [2011] A. Patella, GMOR-like relation in IR-conformal gauge theories, Phys. Rev. D 84, 125033 (2011), arXiv:1106.3494 [hep-th] .
  * Neil and Sitison [2023] E. T. Neil and J. W. Sitison, Model averaging approaches to data subset selection, (2023), arXiv:2305.19417 [stat.ME] .
  * Lepage [1989] G. P. Lepage, The Analysis of Algorithms for Lattice Field Theory, in  _Theoretical Advanced Study Institute in Elementary Particle Physics_ (1989).
  * Neil and Sitison [2022] E. T. Neil and J. W. Sitison, Improved information criteria for Bayesian model averaging in lattice field theory, (2022), arXiv:2208.14983 [stat.ME] .
  * Michael [1994] C. Michael, Fitting correlated data, Phys. Rev. D 49, 2616 (1994), arX

license: Creative Commons - Attribution - https://creativecommons.org/licenses/by/4.0/

authors:
Lattice Strong Dynamics Collaboration: R. C. Brower (1), E. Owen (1),
  C. Rebbi (1), C. Culver (2), D. Schaich (2), K. K. Cushman (3), G. T. Fleming
  (3 and 4), A. Gasbarro (5), A. Hasenfratz (6), E. T. Neil (6), J. Ingoldby
  (7), X. Y. Jin (8), J. C. Osborn (8), E. Rinaldi (9), P. Vranas (10 and 11),
  E. Weinberg (1 and 12), O. Witzel (13) ((1) Boston U., (2) Liverpool U, Dept.
  Math., (3) Yale U., (4) Fermilab, (5) Bern U., (6) Colorado U., (7) ICTP,
  Trieste, (8) Argonne, (9) Wako, RIKEN, (10) LLNL, Livermore, (11) LBNL, NSD,
  (12) NVIDIA, Santa Clara, (13) Siegen U.)

submitter: George T. Fleming


text:
iv:hep-lat/9310026 .
  * Simone [2017] J. N. Simone, Improved data covariance estimation techniques applied to lattice QCD (2017), the 35th International Symposium on Lattice Field Theory, Granada, Spain.
  * Stein [1956] C. Stein, Inadmissibility of the Usual Estimator for the Mean of a Multivariate Normal Distribution, Proc. Third Berkeley Symp. on Math. Statist. and Prob. Berkeley Symp. on Math. Statist. and Prob., 1, 197 (1956).
  * Ledoit and Wolf [2004] O. Ledoit and M. Wolf, A well-conditioned estimator for large-dimensional covariance matrices, Journal of Multivariate Analysis 88, 365 (2004).
  * Schäfer and Strimmer [2005] J. Schäfer and K. Strimmer, A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics, Statistical applications in genetics and molecular biology 4, 10.2202/1544-6115.1175 (2005).
  * Lichtl [2007] A. C. Lichtl, Model optimization in the presence of correlations, in  _Fourth International Workshop on Numerical Analysis and Lattice QCD)_ , QCDNA (Yale University, New Haven, CT, USA., 2007).
  * Cohen  _et al._ [2010] S. D. Cohen, G. T. Fleming, and H.-W. Lin, Exponential time series and lattice quantum field theory, in  _Exponential Data Fitting and Its Applications_, edited by V. Pereyra and G. Scherer (Bentham Science Publishers, 2010) pp. 71–93.
  * Rinaldi  _et al._ [2019] E. Rinaldi, S. Syritsyn, M. L. Wagman, M. I. Buchoff, C. Schroeder, and J. Wasem, Lattice QCD determination of neutron-antineutron matrix elements with physical quark masses, Phys. Rev. D 99, 074510 (2019), arXiv:1901.07519 [hep-lat] .
  * Beane  _et al._ [2021] S. R. Beane _et al._ (NPLQCD, QCDSF), Charged multihadron systems in lattice QCD+QED, Phys. Rev. D 103, 054504 (2021), arXiv:2003.12130 [hep-lat] .
  * Shanahan  _et al._ [2020] P. Shanahan, M. Wagman, and Y. Zhao, Collins-Soper kernel for TMD evolution from lattice QCD, Phys. Rev. D 102, 014511 (2020), arXiv:2003.06063 [hep-lat] .
  * Bazavov  _et al._ [2022] A. Bazavov _et al._ (Fermilab Lattice, MILC), D-meson semileptonic decays to pseudoscalars from four-flavor lattice QCD, (2022), arXiv:2212.12648 [hep-lat] .
  * Ledoit  _et al._ [2018] O. Ledoit, M. Wolf, _et al._ , Optimal estimation of a large-dimensional covariance matrix under Stein’s loss, Bernoulli 24, 3791 (2018).
  * Ledoit and Wolf [2020] O. Ledoit and M. Wolf, Analytical nonlinear shrinkage of large-dimensional covariance matrices, Ann. Statist. 48, 3043 (2020).
  * Burda and Jarosz [2022] Z. Burda and A. Jarosz, Cleaning large-dimensional covariance matrices for correlated samples, Phys. Rev. E 105, 034136 (2022), arXiv:2107.01352 [math-ph] .
  * Appelquist  _et al._ [2013] T. Appelquist _et al._ (Lattice Strong Dynamics (LSD)), Lattice Calculation of Composite Dark Matter Form Factors, Phys. Rev. D 88, 014502 (2013), arXiv:1301.1693 [hep-ph] .
  * Pébay [2008] P. Pébay, Formulas for Robust, One-Pass Parallel Computation of Covariances and Arbitrary-Order Statistical Moments, Sandia Report (2008).