# Scale setting the Möbius domain wall fermion on gradient-flowed HISQ action
using the omega baryon mass and the gradient-flow scales $t_{0}$ and $w_{0}$

Nolan Miller  Department of Physics and Astronomy, University of North
Carolina, Chapel Hill, NC 27516-3255, USA  Logan Carpenter  Department of
Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA  Evan
Berkowitz  Department of Physics, University of Maryland, College Park, MD
20742, USA  Institut für Kernphysik and Institute for Advanced Simulation,
Forschungszentrum Jülich, 54245 Jülich Germany  Chia Cheng Chang (張家丞)
Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS),
RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Nuclear Science Division,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA  Department of
Physics, University of California, Berkeley, CA 94720, USA  Ben Hörz  Nuclear
Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
USA  Dean Howarth  Physics Division, Lawrence Livermore National Laboratory,
Livermore, CA 94550, USA  Nuclear Science Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  Henry Monge-Camacho  Escuela de Físca,
Universidad de Costa Rica, 11501 San José, Costa Rica  Department of Physics
and Astronomy, University of North Carolina, Chapel Hill, NC 27516-3255, USA
Enrico Rinaldi  Arithmer Inc., R&D Headquarters, Minato, Tokyo 106-6040, Japan
Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS),
RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  David A. Brantley  Physics
Division, Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
Christopher Körber  Department of Physics, University of California, Berkeley,
CA 94720, USA  Nuclear Science Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  Chris Bouchard  School of Physics and
Astronomy, University of Glasgow, Glasgow G12 8QQ, UK  M.A. Clark  NVIDIA
Corporation, 2701 San Tomas Expressway, Santa Clara, CA 95050, USA  Arjun
Singh Gambhir  Design Physics Division, Lawrence Livermore National
Laboratory, Livermore, CA 94550, USA  Nuclear Science Division, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, USA  Christopher J. Monahan
Department of Physics, The College of William & Mary, Williamsburg, VA 23187,
USA  Theory Center, Thomas Jefferson National Accelerator Facility, Newport
News, VA 23606, USA  Amy Nicholson  Department of Physics and Astronomy,
University of North Carolina, Chapel Hill, NC 27516-3255, USA  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Pavlos Vranas  Physics Division, Lawrence Livermore National Laboratory,
Livermore, CA 94550, USA  Nuclear Science Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  André Walker-Loud  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Physics Division, Lawrence Livermore National Laboratory, Livermore, CA 94550,
USA  Department of Physics, University of California, Berkeley, CA 94720, USA

( - 1:33)

###### Abstract

We report on a subpercent scale determination using the omega baryon mass and
gradient-flow methods. The calculations are performed on 22 ensembles of
$N_{f}=2+1+1$ highly improved, rooted staggered sea-quark configurations
generated by the MILC and CalLat Collaborations. The valence quark action used
is Möbius domain wall fermions solved on these configurations after a
gradient-flow smearing is applied with a flowtime of $t_{\rm gf}=1$ in lattice
units. The ensembles span four lattice spacings in the range $0.06\lesssim
a\lesssim 0.15$ fm, six pion masses in the range $130\lesssim m_{\pi}\lesssim
400$ MeV and multiple lattice volumes. On each ensemble, the gradient-flow
scales $t_{0}/a^{2}$ and $w_{0}/a$ and the omega baryon mass $am_{\Omega}$ are
computed. The dimensionless product of these quantities is then extrapolated
to the continuum and infinite volume limits and interpolated to the physical
light, strange and charm quark mass point in the isospin limit, resulting in
the determination of $\sqrt{t_{0}}={0.1422(14)}$ fm and $w_{0}={0.1709(11)}$
fm with all sources of statistical and systematic uncertainty accounted for.
The dominant uncertainty in both results is the stochastic uncertainty, though
for $\sqrt{t_{0}}$ there are comparable continuum extrapolation uncertainties.
For $w_{0}$, there is a clear path for a few-per-mille uncertainty just
through improved stochastic precision, as recently obtained by the Budapest-
Marseille-Wuppertal Collaboration.

††preprint: LLNL-JRNL-816949, RIKEN-iTHEMS-Report-20, JLAB-THY-20-3290

##  I Introduction

Lattice QCD (LQCD) has become a prominent theoretical tool for calculations of
hadronic quantities, and many calculations have reached a level of precision
to be able to supplement and/or complement experimental determinations Aoki
_et al._ (2020). Precision calculations of Standard Model processes, for
example, are crucial input for experimental tests of fundamental symmetries in
searches for new physics.

Lattice calculations receive only dimensionless bare parameters as input, so
the output is inherently dimensionless. In some cases, dimensionless
quantities or ratios of quantities may be directly computed without the need
to determine any dimensionful scale. Calculations of $g_{A}$ and
$F_{K}/F_{\pi}$ are examples for which a precise scale setting is not
necessary to make a precise, final prediction. However, there are many
quantities for which a precise scale setting is desirable, such as the hadron
spectrum, the nucleon axial radius, the hadronic contribution to the muon
$g-2$ Aoyama  _et al._ (2020) and many others.

In these cases, a quantity which is dimensionful (after multiplying or
dividing by an appropriate power of the lattice spacing) is calculated and
compared to experiment, following extrapolations to the physical point in
lattice spacing, volume, and pion mass. Because the precision of any
calculations of further dimensionful quantities is limited by the statistical
and systematic uncertainties of this scale setting, quantities which have low
stochastic noise and mild light quark mass dependence, such as the omega
baryon mass $m_{\Omega}$, are preferred. The lattice spacing on each ensemble
may then be determined by comparing the quantity calculated on a given
ensemble to the continuum value.

However, the most precise quantities one may calculate are not necessarily
accessible experimentally. For example, the Sommer scale $r_{0}$ Sommer (1994)
has been one of the most commonly used scales. This scale requires a
determination of the heavy-quark potential which is susceptible to fitting
systematic uncertainties. More recently, the gradient flow scales $t_{0}$
Lüscher (2010) and $w_{0}$ Borsanyi  _et al._ (2012) have been used for a more
precise determination of the lattice spacing Deuzeman and Wenger (2012); Bruno
and Sommer (2014); Dowdall  _et al._ (2013); Bornyakov  _et al._ (2015);
Bazavov  _et al._ (2014a); Blum  _et al._ (2016); Bazavov  _et al._ (2016);
Bruno  _et al._ (2017); Borsanyi  _et al._ (2021). In this case, a well-
controlled extrapolation of these quantities to the physical point is also
necessary.

In this paper we present a precision scale setting for our mixed lattice
action Berkowitz  _et al._ (2017a) which uses $N_{f}=2+1+1$ highly improved,
rooted staggered sea-quark (HISQ) configurations generated by the MILC Bazavov
_et al._ (2013) and CalLat Collaborations and Möbius domain wall fermions for
the valence sector. We compute the dimensionless products
$\sqrt{t_{0}}m_{\Omega}$ and $m_{{\Omega}}w_{0}$ on each ensemble and
extrapolate them to the physical point resulting in the determinations

$\displaystyle\sqrt{t_{0}}m_{\Omega}$
$\displaystyle={1.2051(82)^{s}(15)^{\chi}(46)^{a}(00)^{V}(21)^{\rm
phys}(61)^{M}}$ $\displaystyle={1.205(12)}\,,$
$\displaystyle\frac{\sqrt{t_{0}}}{\rm fm}$
$\displaystyle={0.1422(09)^{s}(02)^{\chi}(05)^{a}(00)^{V}(02)^{\rm
phys}(07)^{M}}$ $\displaystyle={0.1422(14)}\,,$ (1) $\displaystyle
w_{0}m_{\Omega}$
$\displaystyle={1.4483(82)^{s}(15)^{\chi}(45)^{a}(00)^{V}(26)^{\rm
phys}(18)^{M}}$ $\displaystyle={1.4483(97)}$ $\displaystyle\frac{w_{0}}{\rm
fm}$
$\displaystyle={0.1709(10)^{s}(02)^{\chi}(05)^{a}(00)^{V}(03)^{\text{phys}}(02)^{M}}$
$\displaystyle={0.1709(11)}\,,$ (2)

with the statistical ($s$), chiral ($\chi$), continuum-limit ($a$), infinite
volume ($V$), physical-point (phys), and model selection uncertainties ($M$).

We then perform an interpolation of the values of $t_{0}/a^{2}$ and $w_{0}/a$
to the physical quark-mass limit and extrapolation to infinite volume which
allows us to provide a precise, quark mass independent scale setting for each
lattice spacing, with our final results in Table 5. In the following sections
we provide details of our lattice setup, our methods for extrapolation, and
our results with uncertainty breakdown. We conclude with a discussion in the
final section.

##  II Details of the lattice calculation

###  II.1 MDWF on gradient-flowed HISQ

Table 1:  Input parameters for our lattice action. The abbreviated ensemble name Bhattacharya  _et al._ (2015) indicates the approximate lattice spacing in fm and pion mass in MeV. The S, L, XL which come after an ensemble name denote a relatively small, large and extra-large volume with respect to $m_{\pi}L=4$.  Ensemble | $\beta$ | $N_{\rm cfg}$ | Volume | $am_{l}$ | $am_{s}$ | $am_{c}$ | $L_{5}/a$ | $aM_{5}$ | $b_{5},c_{5}$ | $am_{l}^{\rm val}$ | $am^{\rm res}_{l}\times 10^{4}$ | $am_{s}^{\rm val}$ | $am^{\rm res}_{s}\times 10^{4}$ | ${\sigma}$ | $N$ | $N_{\rm src}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
a15m400111Additional ensembles generated by CalLat using the MILC code. The m350 and m400 ensembles were made on the Vulcan supercomputer at LLNL while the a12m310XL, a12m180L, a15m135XL, a09m135, and a06m310L ensembles were made on the Sierra and Lassen supercomputers at LLNL and the Summit supercomputer at OLCF using QUDA Clark  _et al._ (2010); Babich  _et al._ (2011). These configurations are available to any interested party upon request, and will be available for easy anonymous downloading—hopefully soon. | 5.80 | 1000 | $16^{3}\times 48$ | 0.0217 | 0.065 | 0.838 | 12 | 1.3 | 1.50, 0.50 | 0.0278 | 9.365(87) | 0.0902 | 6.937(63) | 3.0 | 30 | 8
a15m35011footnotemark: 1 | 5.80 | 1000 | $16^{3}\times 48$ | 0.0166 | 0.065 | 0.838 | 12 | 1.3 | 1.50, 0.50 | 0.0206 | 9.416(90) | 0.0902 | 6.688(62) | 3.0 | 30 | 16
a15m310 | 5.80 | 1000 | $16^{3}\times 48$ | 0.013 | 0.065 | 0.838 | 12 | 1.3 | 1.50, 0.50 | 0.0158 | 9.563(67) | 0.0902 | 6.640(44) | 4.2 | 45 | 24
a15m310L11footnotemark: 1 | 5.80 | 1000 | $24^{3}\times 48$ | 0.013 | 0.065 | 0.838 | 12 | 1.3 | 1.50, 0.50 | 0.0158 | 9.581(50) | 0.0902 | 6.581(37) | 4.2 | 45 | 4
a15m220 | 5.80 | 1000 | $24^{3}\times 48$ | 0.0064 | 0.064 | 0.828 | 16 | 1.3 | 1.75, 0.75 | 0.00712 | 5.736(38) | 0.0902 | 3.890(25) | 4.5 | 60 | 16
a15m135XL11footnotemark: 1 | 5.80 | 1000 | $48^{3}\times 64$ | 0.002426 | 0.06730 | 0.8447 | 24 | 1.3 | 2.25, 1.25 | 0.00237 | 2.706(08) | 0.0945 | 1.860(09) | 3.0 | 30 | 32
a12m40011footnotemark: 1 | 6.00 | 1000 | $24^{3}\times 64$ | 0.0170 | 0.0509 | 0.635 | 8 | 1.2 | 1.25, 0.25 | 0.0219 | 7.337(50) | 0.0693 | 5.129(35) | 3.0 | 30 | 8
a12m35011footnotemark: 1 | 6.00 | 1000 | $24^{3}\times 64$ | 0.0130 | 0.0509 | 0.635 | 8 | 1.2 | 1.25, 0.25 | 0.0166 | 7.579(52) | 0.0693 | 5.062(34) | 3.0 | 30 | 8
a12m310 | 6.00 | 1053 | $24^{3}\times 64$ | 0.0102 | 0.0509 | 0.635 | 8 | 1.2 | 1.25, 0.25 | 0.0126 | 7.702(52) | 0.0693 | 4.950(35) | 3.0 | 30 | 8
a12m310XL11footnotemark: 1 | 6.00 | 1000 | $48^{3}\times 64$ | 0.0102 | 0.0509 | 0.635 | 8 | 1.2 | 1.25, 0.25 | 0.0126 | 7.728(22) | 0.0693 | 4.927(21) | 3.0 | 30 | 8
a12m220S | 6.00 | 1000 | $24^{4}\times 64$ | 0.00507 | 0.0507 | 0.628 | 12 | 1.2 | 1.50, 0.50 | 0.00600 | 3.990(42) | 0.0693 | 2.390(24) | 6.0 | 90 | 4
a12m220 | 6.00 | 1000 | $32^{3}\times 64$ | 0.00507 | 0.0507 | 0.628 | 12 | 1.2 | 1.50, 0.50 | 0.00600 | 4.050(20) | 0.0693 | 2.364(15) | 6.0 | 90 | 4
a12m220ms | 6.00 | 1000 | $32^{3}\times 64$ | 0.00507 | 0.0304 | 0.628 | 12 | 1.2 | 1.50, 0.50 | 0.00600 | 3.819(26) | 0.0415 | 2.705(20) | 6.0 | 90 | 8
a12m220L | 6.00 | 1000 | $40^{3}\times 64$ | 0.00507 | 0.0507 | 0.628 | 12 | 1.2 | 1.50, 0.50 | 0.00600 | 4.040(26) | 0.0693 | 2.361(19) | 6.0 | 90 | 4
a12m180L11footnotemark: 1 | 6.00 | 1000 | $48^{3}\times 64$ | 0.00339 | 0.0507 | 0.628 | 14 | 1.2 | 1.75, 0.75 | 0.00380 | 3.038(13) | 0.0693 | 1.888(11) | 3.0 | 30 | 16
a12m130 | 6.00 | 1000 | $48^{3}\times 64$ | 0.00184 | 0.0507 | 0.628 | 20 | 1.2 | 2.00, 1.00 | 0.00195 | 1.642(09) | 0.0693 | 0.945(08) | 3.0 | 30 | 32
a09m40011footnotemark: 1 | 6.30 | 1201 | $32^{3}\times 64$ | 0.0124 | 0.037 | 0.44 | 6 | 1.1 | 1.25, 0.25 | 0.0160 | 2.532(23) | 0.0491 | 1.957(17) | 3.5 | 45 | 8
a09m35011footnotemark: 1 | 6.30 | 1201 | $32^{3}\times 64$ | 0.00945 | 0.037 | 0.44 | 6 | 1.1 | 1.25, 0.25 | 0.0121 | 2.560(24) | 0.0491 | 1.899(16) | 3.5 | 45 | 8
a09m310 | 6.30 | 780 | $32^{3}\times 96$ | 0.0074 | 0.037 | 0.44 | 6 | 1.1 | 1.25, 0.25 | 0.00951 | 2.694(26) | 0.0491 | 1.912(15) | 6.7 | 167 | 8
a09m220 | 6.30 | 1001 | $48^{3}\times 96$ | 0.00363 | 0.0363 | 0.43 | 8 | 1.1 | 1.25, 0.25 | 0.00449 | 1.659(13) | 0.0491 | 0.834(07) | 8.0 | 150 | 6
a09m13511footnotemark: 1 | 6.30 | 1010 | $64^{3}\times 96$ | 0.001326 | 0.03636 | 0.4313 | 12 | 1.1 | 1.50, 0.50 | 0.00152 | 0.938(06) | 0.04735 | 0.418(04) | 3.5 | 45 | 16
a06m310L11footnotemark: 1 | 6.72 | 1000 | $72^{3}\times 96$ | 0.0048 | 0.024 | 0.286 | 6 | 1.0 | 1.25, 0.25 | 0.00617 | 0.225(03) | 0.0309 | 0.165(02) | 3.5 | 45 | 8

The lattice action we use is the mixed-action Renner  _et al._ (2005); Bar
_et al._ (2003) with Möbius Brower  _et al._ (2017) domain wall fermions
Kaplan (1992); Shamir (1993); Furman and Shamir (1995) solved on $N_{f}=2+1+1$
highly improved staggered quarks Follana  _et al._ (2007) after they are
gradient-flow smeared Narayanan and Neuberger (2006); Lüscher and Weisz
(2011); Lüscher (2013) (corresponding to an infinitesimal stout-smearing
procedure Morningstar and Peardon (2004)) to a flow-time of $t_{\rm
gf}/a^{2}=1$ Berkowitz  _et al._ (2017a). The choice to hold the flow-time
fixed in lattice units is important to ensure that as the continuum limit is
taken, effects arising from finite flow-time also extrapolate to zero.

This action has been used to compute the nucleon axial coupling, $g_{A}$, with
a 1% total uncertainty Bouchard  _et al._ (2017); Berkowitz  _et al._ (2017b);
Chang  _et al._ (2018); Berkowitz  _et al._ (2018a), the
$\pi^{-}\rightarrow\pi^{+}$ matrix elements relevant to neutrinoless double
beta-decay Nicholson  _et al._ (2018) and most recently, $F_{K}/F_{\pi}$
Miller  _et al._ (2020). Our calculation of $F_{K}/F_{\pi}$ was obtained with
a total uncertainty of 0.4% which provides an important benchmark for our
action, as the result is consistent with other determinations in the
literature Follana  _et al._ (2008); Blossier  _et al._ (2009); Durr  _et al._
(2010); Dowdall  _et al._ (2013); Blum  _et al._ (2016); Carrasco  _et al._
(2015); Bazavov  _et al._ (2014b); Dürr  _et al._ (2017); Bornyakov  _et al._
(2017); Bazavov  _et al._ (2018) (and the FLAG average Aoki  _et al._ (2020)),
and also contributes to the test of the universality of lattice QCD results in
the continuum limit.

Our plan to compute the axial and other elastic form factors of the nucleon
with this mixed-action, as well as other quantities, leads to a desire to have
a scale setting with sufficiently small uncertainty that it does not increase
the final uncertainty of such quantities. It has been previously observed that
both $w_{0}$ Borsanyi  _et al._ (2012); Bazavov  _et al._ (2016) and the omega
baryon mass Allton  _et al._ (2008); Walker-Loud  _et al._ (2009); Lin  _et
al._ (2009); Durr  _et al._ (2008); Alexandrou  _et al._ (2009, 2014);
Brantley  _et al._ (2016); Borsanyi  _et al._ (2021) have mild quark mass
dependence and that they can be determined with high statistical precision
with relatively low computational cost. The input parameters of our action on
all ensembles are provided in Table 1.

###  II.2 Correlation function construction and analysis

Figure 1:  Stability plots of the ground state omega baryon mass on the three
physical pion mass ensembles. The left plots show the effective mass data (in
lattice units) and reconstructed effective mass from the chosen fit for both
the SS and PS correlation functions. The dark gray and colored band are
displayed for the region of time used in the analysis and an extrapolation
beyond $t_{\rm max}$ is shown after a short break in the fit band. The
horizontal gray band is the prior used for the ground state mass. The right
plots show the corresponding value of $E_{0}$ as a function of $t_{\rm min}$
and the number of states $n_{s}$ used in the analysis, as well as the
corresponding $Q$ value and relative weight as a function of $n_{s}$ for a
given $t_{\rm min}$, where the weight is set by the Bayes Factor. See Appendix
A for more detail on the selection of the final fit. The chosen fit is denoted
with a filled black symbol and the horizontal band is the value of $E_{0}$
from the chosen fit. The $y$-range of the upper panel of the stability plots
is equal to the prior of the ground state energy (the horizontal gray band in
the left plot).

For the scale setting computation, we have to determine four or five
quantities on each ensemble, the pion, kaon and omega masses, the gradient-
flow scale $w_{0}$ and the pion decay constant $F_{\pi}$. For $m_{\pi}$,
$m_{K}$ and $F_{\pi}$, we take the values from our $F_{K}/F_{\pi}$ computation
for the 18 ensembles in common. For the four new ensembles in this work
(a15m310L, a12m310XL, a12m220ms, a12m180L), we follow the same analysis
strategy described in Ref. Miller  _et al._ (2020).

The a12m220ms ensemble is identical to a12m220 except that the strange quark
mass is roughly 60% of the physical value rather than being near the physical
value. The a15m310L ensemble has identical input parameters as the a15m310
ensemble but $L=24$ (3.6 fm) instead of $L=16$ (2.4 fm), while the a12m310XL
ensemble is identical to the a12m310 ensemble but with $L=48$ (5.8 fm) instead
of $L=24$ (2.9 fm). The a12m180L and a12m310XL ensembles have a lattice volume
that is the same size as a12m130 but pion masses of roughly $m_{\pi}\simeq
180$ and 310 MeV. These new ensembles provide important lever arms for the
various extrapolations. The a12m220ms provides a unique lever arm for varying
the strange quark mass significantly from its physical value, the a15m310L and
a12m310XL provide other pion masses where we can perform a volume study and
the a12m180L ensemble provides an additional light pion mass ensemble to help
with the physical pion mass extrapolation. The first of these is important for
this scale setting while the latter three will be more important for future
work.

The omega baryon correlation functions are constructed similarly to the pion
and kaon. A source for the propagator is constructed with the gauge invariant
Gaussian smearing routine in Chroma Edwards and Joo (2005)
(GAUGE_INV_GAUSSIAN). Then, correlation functions are constructed using both a
point sink as well as the same gauge invariant Gaussian smearing routine with
the same parameters as the source. The values of the “smearing width”
(${\sigma}$) and the number of iterations ($N$) used to approximate the
exponential smearing profile are provided in Table 1. The correlation
functions constructed with the point sink are referred to as PS and those with
the smeared sink as SS.

Local spin wave functions are constructed following Refs. Basak  _et al._
(2005a, b). Both positive- and negative-parity omega-baryon correlation
functions are constructed with the upper and lower spin components of the
quark propagators in the Dirac basis. The negative-parity correlation
functions are time-reversed with an appropriate sign flip of the correlation
function, effectively doubling the statistics with no extra inversions. The
four different spin projections of the omega are averaged as well to produce
the final spin and parity averaged two-point correlation functions.

The reader will notice that the values of ${\sigma}$ and $N$ do not follow an
obvious pattern. This is because in our first computations of $g_{A}$ Bouchard
_et al._ (2017); Berkowitz  _et al._ (2017b), we applied an “aggressive”
smearing with a larger value of ${\sigma}$ and correspondingly larger number
of iterations, which led to a large suppression of excited states, but also
showed evidence of “over smearing” such that the non-positive-definite PS
correlation functions displayed symptoms of having a relatively large negative
overlap factor for excited states (there were wiggles in the PS effective
masses). In a subsequent paper studying the two-nucleon system on the a12m350
ensemble Berkowitz  _et al._ (2017c), where we utilized matrix Prony Beane
_et al._ (2009) to form linear combinations of PS and SS nucleons to construct
a “calm” nucleon which is ground-state dominated earlier in time, we observed
that using a milder smearing with smaller width and fewer iterations provided
a much more stable extraction of the ground state and did not show signs of
large negative overlap factors. Hence, many but not all of the ensembles have
been rerun with our improved choices of ${\sigma}$ and $N$. We have observed
the choice ${\sigma}=3.0$ and $N=30$ works well for the a15 and a12 ensembles
and that ${\sigma}=3.5$ with $N=45$ works well for the a09 and a06 ensembles.

In order to determine the omega baryon mass on each ensemble, we perform a
stability analysis of the extracted ground state mass as a function of $t_{\rm
min}$ used in the fit as well as the number of states used in the analysis.
The correlation functions are analyzed in a Bayesian framework with
constraints Lepage  _et al._ (2002). We choose normally distributed priors for
the ground-state energy and all overlap factors, and log-normal distributions
for excited-state energy priors. The ground-state energy and overlap factors
are motivated by the plateau values of the effective masses with the priors
taken to be roughly 10 times larger than the stochastic uncertainty of the
respective effective mass data in the plateau region. The excited-state energy
splittings are set to the value of two pion masses with a width allowing for
fluctuations down to one pion mass within one standard deviation.

In Fig. 1 we show sample extractions of the ground state mass on our three
physical pion mass ensembles. In the left plot, we show the effective mass
data from the two correlation functions. The weights are normalized on a given
time slice by the largest Bayes factor at that $t_{\rm min}$ value. We have
not implemented a more thorough algorithm to weight fits against each other
that utilize different amounts of data, as described for example in Ref. Jay
and Neil (2020). Rather, we have chosen a fit for a given ensemble (the filled
black symbol in the right panels highlighted by the horizontal colored band)
that has a good fit quality, the maximum or near maximum relative weight, and
consistency with the late-time data. We tried to optimize this choice over all
ensembles simultaneously, with $t_{\rm min}$ held nearly fixed for a given
lattice spacing, rather than hand-picking the optimal fit on each ensemble
separately, in order to minimize the possible bias introduced by analysis
choices. Good fits are obtained on all ensembles with $n_{s}=2$, simplifying
the model function and reducing the chance of overfitting the correlation
functions, which is most relevant on ensembles with the more aggressive
choices of smearing parameters. In Appendix C, we show the corresponding
stability plots for all remaining ensembles. In Table 2 we show the resulting
values of $am_{\Omega}$ on all ensembles used in this work.

###  II.3 Calculation of $t_{0}$ and $w_{0}$

In order to efficiently compute the value of $t_{0}/a^{2}$ and $w_{0}/a$ on
each ensemble, we have implemented the Symanzik flow in the QUDA library
Collaboration ; Clark  _et al._ (2010); Babich  _et al._ (2011). We used the
tree-level improved action and the symmetric, cloverleaf definition of the
field-strength tensor, following the MILC implementation mil . We used a step
size of ${\epsilon}=0.01$ in the Runge-Kutta algorithm proposed by Lüscher
Lüscher (2010), which leads to negligibly small integration errors. The scales
$t_{0}$ and $w_{0}$ are defined by the equations

$\displaystyle t^{2}\langle E(t)\rangle\Big{|}_{t=t_{0,\rm orig}}=0.3\,,$
$\displaystyle W_{\rm orig}\equiv t\frac{d}{dt}(t^{2}\langle
E(t)\rangle)\Big{|}_{t=w_{0,\rm orig}^{2}}=0.3\,,$ (3)

where $\langle E(t)\rangle$ is the gluonic action density at flow time $t$. In
Fig. 2, we show the determination of the $w_{0,\rm orig}/a$ on the two
physical pion mass ensembles that we have generated. The uncertainties are
determined by observing a saturation of the uncertainty as the bin-size is
increased when binning the results from configurations close in Monte Carlo
time. These uncertainties were cross-checked with an autocorrelation study
using the ${\Gamma}$-method Wolff (2004) implemented in the unew Python
package De Palma  _et al._ (2019).

Figure 2:  Determination of $w_{0,\rm orig}$ on the two new physical pion mass
ensembles. Table 2:  The omega baryon mass ($m_{\Omega}$) and gradient flow
scales ($t_{0,\rm orig}$, $t_{0,\rm imp}$, $w_{0,\rm orig}$ and $w_{0,\rm
imp}$) determined on each ensemble are listed as well as their dimensionless
products. Additionally, in the bottom panel, we list the parameters used to
control the physical point extrapolation:
$l_{\Omega}^{2}=m_{\pi}^{2}/m_{\Omega}^{2}$,
$s_{\Omega}^{2}=(2m_{K}^{2}-m_{\pi}^{2})/m_{\Omega}^{2}$,
$l_{F}^{2}=m_{\pi}^{2}/(4\pi F_{\pi})^{2}$,
$s_{F}^{2}=(2m_{K}^{2}-m_{\pi}^{2})/(4\pi F_{\pi})^{2}$, $m_{\pi}L$,
${\epsilon}_{a}=a/(2w_{0,\rm orig})$ and ${\alpha}_{S}$. The values of
${\alpha}_{S}$ are taken from Table III of Ref. Bazavov  _et al._ (2016) which
were determined with a heavy-quark potential method Davies  _et al._ (2003).
An HDF5 file is provided with this publication which includes the resulting
bootstrap samples of all these quantities which can be used to construct the
correlated uncertainties.

Ensemble | $am_{\Omega}$ | $t_{0,\rm orig}/a^{2}$ | $t_{0,\rm imp}/a^{2}$ | $\sqrt{t_{0,\rm orig}}m_{\Omega}$ | $\sqrt{t_{0,\rm imp}}m_{\Omega}$ | $w_{0,\rm orig}/a$ | $w_{0,\rm imp}/a$ | $w_{0,\rm orig}m_{\Omega}$ | $w_{0,\rm imp}m_{\Omega}$
---|---|---|---|---|---|---|---|---|---
a15m400 | 1.2437(43) | 1.1905(15) | 0.9563(11) | 1.3570(47) | 1.2162(42) | 1.1053(12) | 1.0868(14) | 1.3747(49) | 1.3516(50)
a15m350 | 1.2331(31) | 1.2032(14) | 0.9660(11) | 1.3526(35) | 1.2120(32) | 1.1154(11) | 1.0988(13) | 1.3754(37) | 1.3550(38)
a15m310L | 1.2287(31) | 1.2111(07) | 0.9720(05) | 1.3521(34) | 1.2113(30) | 1.1219(06) | 1.1066(07) | 1.3784(35) | 1.3596(35)
a15m310 | 1.2312(36) | 1.2121(09) | 0.9727(07) | 1.3555(40) | 1.2143(36) | 1.1230(07) | 1.1079(09) | 1.3826(41) | 1.3640(41)
a15m220 | 1.2068(26) | 1.2298(08) | 0.9861(07) | 1.3383(30) | 1.1984(27) | 1.1378(07) | 1.1255(08) | 1.3731(31) | 1.3582(31)
a15m135XL | 1.2081(19) | 1.2350(04) | 0.9897(03) | 1.3425(21) | 1.2018(19) | 1.1432(03) | 1.1319(04) | 1.3811(22) | 1.3674(22)
a12m400 | 1.0279(25) | 1.6942(11) | 1.4172(10) | 1.3380(33) | 1.2237(30) | 1.3616(07) | 1.3636(08) | 1.3997(35) | 1.4017(35)
a12m350 | 1.0139(26) | 1.7091(14) | 1.4298(12) | 1.3255(35) | 1.2124(32) | 1.3728(10) | 1.3755(11) | 1.3918(38) | 1.3946(38)
a12m310XL | 1.0072(41) | 1.7221(07) | 1.4407(06) | 1.3217(54) | 1.2089(50) | 1.3831(05) | 1.3865(06) | 1.3930(58) | 1.3964(58)
a12m310 | 1.0112(32) | 1.7213(19) | 1.4398(17) | 1.3267(42) | 1.2134(39) | 1.3830(13) | 1.3863(12) | 1.3985(46) | 1.4019(46)
a12m220ms | 0.8896(92) | 1.7891(15) | 1.4977(13) | 1.190(12)0 | 1.089(11)0 | 1.4339(12) | 1.4406(12) | 1.276(13)0 | 1.282(13)0
a12m220S | 0.9970(26) | 1.7466(20) | 1.4614(17) | 1.3177(35) | 1.2053(32) | 1.4021(15) | 1.4069(16) | 1.3980(39) | 1.4027(39)
a12m220L | 0.9944(30) | 1.7489(09) | 1.4633(08) | 1.3150(40) | 1.2028(37) | 1.4041(06) | 1.4090(07) | 1.3962(43) | 1.4010(43)
a12m220 | 0.9924(60) | 1.7498(14) | 1.4641(12) | 1.3127(80) | 1.2007(73) | 1.4047(10) | 1.4096(11) | 1.3940(85) | 1.3988(86)
a12m180L | 0.9924(26) | 1.7553(05) | 1.4686(05) | 1.3148(35) | 1.2026(32) | 1.4093(05) | 1.4145(05) | 1.3985(38) | 1.4037(38)
a12m130 | 0.9801(26) | 1.7628(07) | 1.4749(06) | 1.3013(34) | 1.1903(31) | 1.4155(05) | 1.4211(06) | 1.3873(37) | 1.3928(37)
a09m400 | 0.7716(23) | 2.9158(42) | 2.6040(33) | 1.3176(41) | 1.2451(38) | 1.8602(26) | 1.8686(27) | 1.4353(48) | 1.4418(48)
a09m350 | 0.7561(35) | 2.9455(37) | 2.6301(34) | 1.2977(61) | 1.2262(58) | 1.8810(25) | 1.8900(26) | 1.4222(69) | 1.4291(70)
a09m310 | 0.7543(36) | 2.9698(32) | 2.6521(29) | 1.2998(63) | 1.2283(59) | 1.8970(22) | 1.9066(22) | 1.4308(71) | 1.4381(71)
a09m220 | 0.7377(30) | 3.0172(16) | 2.6952(15) | 1.2814(53) | 1.2111(50) | 1.9282(12) | 1.9388(12) | 1.4224(59) | 1.4302(60)
a09m135 | 0.7244(25) | 3.0390(12) | 2.7147(11) | 1.2629(44) | 1.1936(42) | 1.9450(10) | 1.9563(11) | 1.4091(50) | 1.4172(50)
a06m310L | 0.5069(21) | 6.4079(45) | 6.0606(44) | 1.2830(54) | 1.2478(52) | 2.8958(20) | 2.9053(20) | 1.4678(62) | 1.4726(62)

Ensemble | $l_{F}^{2}$ | $s_{F}^{2}$ | $l_{\Omega}^{2}$ | $s_{\Omega}^{2}$ | $m_{\pi}L$ | ${\epsilon}_{a}^{2}$ | ${\alpha}_{S}$
---|---|---|---|---|---|---|---
a15m400 | 0.09216(33) | 0.2747(10) | 0.05928(42) | 0.1767(12) | 4.85 | 0.20462(44) | 0.58801
a15m350 | 0.07505(28) | 0.2915(09) | 0.04609(25) | 0.1790(09) | 4.24 | 0.20096(39) | 0.58801
a15m310L | 0.06018(23) | 0.2984(11) | 0.03630(18) | 0.1800(09) | 5.62 | 0.19864(20) | 0.58801
a15m310 | 0.06223(17) | 0.3035(09) | 0.03675(23) | 0.1792(11) | 3.78 | 0.19825(26) | 0.58801
a15m220 | 0.03269(11) | 0.3253(09) | 0.01877(09) | 0.1868(08) | 3.97 | 0.19311(23) | 0.58801
a15m135XL | 0.01319(05) | 0.3609(11) | 0.00726(03) | 0.1986(06) | 4.94 | 0.19129(10) | 0.58801
a12m400 | 0.08889(30) | 0.2648(10) | 0.05610(28) | 0.1671(08) | 5.84 | 0.13484(15) | 0.53796
a12m350 | 0.07307(37) | 0.2810(13) | 0.04454(24) | 0.1713(09) | 5.14 | 0.13266(19) | 0.53796
a12m310XL | 0.05904(23) | 0.2909(12) | 0.03506(29) | 0.1727(14) | 9.05 | 0.13069(10) | 0.53796
a12m310 | 0.05984(25) | 0.2933(12) | 0.03482(23) | 0.1707(11) | 4.53 | 0.13071(24) | 0.53796
a12m220ms | 0.03400(16) | 0.2000(09) | 0.02229(46) | 0.1311(27) | 4.25 | 0.12159(20) | 0.53796
a12m220S | 0.03384(19) | 0.3210(19) | 0.01849(13) | 0.1754(09) | 3.25 | 0.12716(27) | 0.53796
a12m220L | 0.03289(15) | 0.3195(14) | 0.01816(12) | 0.1765(11) | 5.36 | 0.12681(11) | 0.53796
a12m220 | 0.03314(15) | 0.3202(15) | 0.01831(23) | 0.1769(22) | 4.30 | 0.12670(18) | 0.53796
a12m180L | 0.02277(09) | 0.3319(13) | 0.01220(07) | 0.1779(10) | 5.26 | 0.12588(08) | 0.53796
a12m130 | 0.01287(08) | 0.3429(14) | 0.00687(05) | 0.1832(10) | 3.90 | 0.12477(09) | 0.53796
a09m400 | 0.08883(32) | 0.2638(09) | 0.05512(34) | 0.1637(10) | 5.80 | 0.07225(20) | 0.43356
a09m350 | 0.07256(32) | 0.2827(11) | 0.04358(42) | 0.1698(16) | 5.05 | 0.07066(19) | 0.43356
a09m310 | 0.06051(22) | 0.2946(10) | 0.03481(34) | 0.1695(16) | 4.50 | 0.06947(16) | 0.43356
a09m220 | 0.03307(14) | 0.3278(13) | 0.01761(15) | 0.1746(14) | 4.70 | 0.06724(08) | 0.43356
a09m135 | 0.01346(08) | 0.3500(17) | 0.00674(05) | 0.1752(12) | 3.81 | 0.06608(07) | 0.43356
a06m310L | 0.06141(35) | 0.2993(17) | 0.03481(29) | 0.1696(14) | 6.81 | 0.02981(04) | 0.29985

Reference Fodor  _et al._ (2014) determined the tree-level in lattice
perturbation theory improvement coefficients for the determination of these
gradient flow scales through $\mathrm{O}(a^{8}/t^{4})$ for various choices of
the gauge action, the gradient flow action and the definition of the field-
strength tensor. As defined in Ref. Fodor  _et al._ (2014), we have
implemented the Symanzik-Symanzik-Clover (SSC) scheme with the relevant
improvement coefficients (see Table 1 of Ref. Fodor  _et al._ (2014))

$\displaystyle C_{2}=-\frac{19}{72}\,,$ $\displaystyle
C_{4}=\frac{145}{1536}\,,$ $\displaystyle C_{6}=-\frac{12871}{276480}\,,$
$\displaystyle C_{8}=\frac{52967}{1769472}\,.$ (4)

One can then determine improved scales, $t_{0,\rm imp}$ and $w_{0,\rm imp}$ in
which one has perturbatively removed the leading discretization effects in
these flow observables,

$\displaystyle\frac{t^{2}\langle
E(t)\rangle}{1+\sum_{n}C_{2n}\frac{a^{2n}}{t^{n}}}\bigg{|}_{t=t_{0,\rm
imp}}=0.3\,,$ $\displaystyle t\frac{d}{dt}\left(\frac{t^{2}\langle
E(t)\rangle}{1+\sum_{n}C_{2n}\frac{a^{2n}}{t^{n}}}\right)\bigg{|}_{t=w_{0,\rm
imp}^{2}}=0.3\,,$ (5)

In the present work, we explore using both the original and improved versions
of $t_{0}$ and $w_{0}$ when performing our scale setting analysis. For the
improved versions, we have implemented the fourth order improvement [up to and
including the $C_{8}(a^{2}/t)^{4}$ correction]. This is the same
implementation performed by MILC in Ref. Bazavov  _et al._ (2016).

In Table 2, we list the values of $t_{0}/a^{2}$ and $w_{0}/a$ for the original
and improved definitions as well as the dimensionless products of
$\sqrt{t_{0}}m_{\Omega}$ and $w_{0}m_{\Omega}$.

##  III Extrapolation Functions

This work utilizes 22 different ensembles, each with $\mathrm{O}(1000)$
configurations (Table 1), to control the systematic uncertainties in the LQCD
calculation of the scale. This allows us to address:

  1. 1.

The physical light and strange quark mass limit;

  2. 2.

The physical charm quark mass limit;

  3. 3.

The continuum limit;

  4. 4.

The infinite volume limit.

####  III.0.1 Physical light and strange quark mass limit

The ensembles have a range of light quark masses which correspond roughly to
$130\lesssim m_{\pi}\lesssim 400$ MeV. We have three lattice spacings at
$m_{\pi}\simeq m_{\pi}^{\rm phys}$ such that the light quark mass
extrapolation is really an interpolation. On all but one of the 22 ensembles,
the strange quark mass is close to its physical value, allowing us to perform
a simple interpolation to the physical strange quark mass point. One ensemble
has a strange quark mass of roughly 2/3 its physical value (a12m220ms),
allowing us to explore systematics in this strange quark mass interpolation.

To parametrize the light and strange quark mass dependence, we utilize two
sets of small parameters

$\displaystyle{\Lambda}=F:$ $\displaystyle
l_{F}^{2}=\frac{m_{\pi}^{2}}{{\Lambda}_{\chi}^{2}},$ $\displaystyle
s_{F}^{2}=\frac{2m_{K}^{2}-m_{\pi}^{2}}{{\Lambda}_{\chi}^{2}},$ (6)
$\displaystyle{\Lambda}={\Omega}:$ $\displaystyle
l_{\Omega}^{2}=\frac{m_{\pi}^{2}}{m_{\Omega}^{2}},$ $\displaystyle
s_{\Omega}^{2}=\frac{2m_{K}^{2}-m_{\pi}^{2}}{m_{\Omega}^{2}},$ (7)

where we have defined

${\Lambda}_{\chi}\equiv 4\pi F_{\pi}\,.$ (8)

Using the Gell-Mann–Oakes–Renner relation Gell-Mann  _et al._ (1968), the
numerators in these parameters correspond roughly to the light and strange
quark mass, $m_{\pi}^{2}=2B\hat{m}$ and $2m_{K}^{2}-m_{\pi}^{2}=2Bm_{s}$,
where $\hat{m}=\frac{1}{2}(m_{u}+m_{d})$. The first set of parameters, Eq.
(6), is inspired by $\chi$PT and commonly used as a set of small expansion
parameters in extrapolating LQCD results. The second set of small parameters,
Eq. (7), is inspired by Ref. Lin  _et al._ (2009). In Fig. 3, we plot the
values of these parameters in comparison with the physical point. Since we are
working in the isospin limit in this work, we define the physical point as

$\displaystyle m_{\pi}^{\rm phys}$
$\displaystyle=\bar{M}_{\pi}=134.8(3)\textrm{ MeV},$ $\displaystyle m_{K}^{\rm
phys}$ $\displaystyle=\bar{M}_{K}=494.2(3)\textrm{ MeV},$ $\displaystyle
F_{\pi}^{\rm phys}$ $\displaystyle=F_{\pi^{+}}^{\rm phys}=92.07(57)\textrm{
MeV},$ $\displaystyle m_{\Omega}^{\rm phys}$
$\displaystyle=1672.43(32)\textrm{ MeV}\,,$ (9)

with the first three values from the FLAG report Aoki  _et al._ (2020) and the
omega baryon mass from the PDG Tanabashi  _et al._ (2018). The values of
$l_{F,{\Omega}}$ and $s_{F,{\Omega}}$ are given in Table 2 for all ensembles.

Figure 3:  The parameter space of $s_{F}^{2}$ versus $l_{F}^{2}$ (left) and
$s_{\Omega}^{2}$ versus $l_{\Omega}^{2}$ used in this calculation. The
vertical and horizontal gray lines represent the physical point defined by Eq.
(III.0.1).

####  III.0.2 Physical charm quark mass limit

The FNAL and MILC Collaborations have provided a determination of the input
value of the charm quark mass that reproduces the “physical” charm quark mass
for each of the four lattice spacings used in this work. The mass of the
$D_{s}$ meson was used to tune the input charm quark mass until the physical
$D_{s}$ mass was reproduced (with the already tuned values of the input
strange quark mass), defining the “physical” charm quark masses Bazavov  _et
al._ (2014b),

$\displaystyle a_{\rm 15}m_{c}^{\rm phys}=0.8447(15),$ $\displaystyle a_{\rm
12}m_{c}^{\rm phys}=0.6328(8),$ $\displaystyle a_{\rm 09}m_{c}^{\rm
phys}=0.4313(6),$ $\displaystyle a_{\rm 06}m_{c}^{\rm phys}=0.2579(4).$ (10)

Comparing to Table 1, the simulated charm quark mass is mistuned by less than
2% of the physical charm quark mass for all ensembles used in this work except
the a06m310L ensemble, whose simulated charm quark is almost 10% heavier than
its physical value.

In order to test the sensitivity of our results to this mistuning of the charm
quark mass, we perform a reweighting Ferrenberg and Swendsen (1988) study of
the a06m310L correlation functions and extracted pion, kaon and omega baryon
masses. While the relative shift of the charm quark mass is small, this shift
is approximately equal to the value of the physical strange quark mass,

$\displaystyle{\delta}a_{06}m_{c}$ $\displaystyle=a_{06}m_{c}^{\rm
phys}-a_{06}m_{c}^{\rm HMC}$ $\displaystyle=-0.0281(4)\,.$ (11)

As the reweighting factor is provided by a ratio of the charm quark fermion
determinant, it is an extensive quantity, and the relatively large volume we
have used to generate the a06m310L ensemble causes some challenges in
accurately determining the reweighting factor.

The summary of our study is that our scale setting is not sensitive to this
mistuning of the charm quark mass, in line with prior expectation. For
example, we find

$\displaystyle a_{\rm
06}m_{\Omega}=\left\\{\begin{array}[]{cl}0.5069(21)\,,&\textrm{unweighted}\\\
0.5065(29)\,,&\textrm{reweighted}\\\ \end{array}\right.\,,$ (14)
$\displaystyle a_{\rm 06}(m_{\Omega}^{\rm reweighted}-m_{\Omega}^{\rm
unweighted})=-0.0004(34)\,,$ (15)

where the splitting is determined under bootstrap. We provide extensive
details in Appendix A.

####  III.0.3 Continuum limit

In order to control the continuum extrapolation, we utilize four lattice
spacings ranging from $0.057\lesssim a\lesssim 0.15$ fm. For most of the pion
masses, we have three values of $a$ with four values at $m_{\pi}\sim 310$ MeV
and one value at $m_{\pi}\sim 180$ MeV. The parameter space is depicted in
Fig. 4. The small dimensionless parameter we utilize to extrapolate to the
continuum limit is

${\epsilon}_{a}^{2}=\frac{a^{2}}{(2w_{0,\rm orig})^{2}}\,.$ (16)

As noted in Ref. Miller  _et al._ (2020), this choice is convenient as the
values of ${\epsilon}_{a}^{2}$ span a similar range as $l_{F}^{2}$. This
allows us to test the ansatz of our assumed power counting that treats
corrections of $\mathrm{O}(l_{F}^{2})=\mathrm{O}({\epsilon}_{a}^{2})$ which we
found to be the case for $F_{K}/F_{\pi}$ Miller  _et al._ (2020).

An equally valid way to define the small parameter characterizing the
discretization corrections is to utilize the gradient flow scale that is also
used to define the observable $y$ being extrapolated. The following
normalizations are comparable to our standard choice, Eq. (16)

${\epsilon}_{a}^{2}=\left\\{\begin{array}[]{rl}\frac{a^{2}}{(2w_{0,\rm
orig})^{2}}\,,&y=w_{0,\rm orig}m_{\Omega}\\\
\vspace{1pt}\frac{a^{2}}{(2w_{0,\rm imp})^{2}}\,,&y=w_{0,\rm imp}m_{\Omega}\\\
\vspace{1pt}\frac{a^{2}}{4t_{0,\rm orig}}\,,&y=\sqrt{t_{0,\rm
orig}}m_{\Omega}\\\ \vspace{1pt}\frac{a^{2}}{4t_{0,\rm
imp}}\,,&y=\sqrt{t_{0,\rm imp}}m_{\Omega}\\\ \end{array}\right.\,.$ (17)

While these choices do not exhaust the possibilities, they are used to explore
possible systematic uncertainties arising from this choice. Unless otherwise
noted, the fixed choice in Eq. (16) is used in subsequent results and plots.

Figure 4:  Parameter space of pion mass and lattice spacing utilized in this
work expressed in terms of $l_{F}^{2}$ and ${\epsilon}_{a}^{2}$.

####  III.0.4 Infinite volume limit

The leading sensitivity of $m_{\Omega}$, $t_{0}$ and $w_{0}$ to the size of
the volume is exponentially suppressed for sufficiently large $m_{\pi}L$
Luscher (1986). We have ensembles with multiple volumes at a15m310, a12m310
and a12m220 to test the predicted finite volume corrections against the
observed ones. We derive the predicted volume dependence of $w_{0}m_{\Omega}$
to the first two nontrivial orders in Sec. III.2.

###  III.1 Light and strange quark mass dependence

The light and strange quark mass dependence of the omega baryon has been
derived in $\mathrm{SU}(3)$ heavy baryon $\chi$PT (HB$\chi$PT) Jenkins and
Manohar (1991a, b) to next-to-next-to-leading order (N2LO) which is
$\mathrm{O}(m_{\pi,K,\eta}^{4})$ Lebed (1994); Lebed and Luty (1994); Tiburzi
and Walker-Loud (2005). It has been shown that $\mathrm{SU}(3)$ HB$\chi$PT
does not produce a converging expansion at the physical quark masses Walker-
Loud  _et al._ (2009); Torok  _et al._ (2010); Ishikawa  _et al._ (2009);
Jenkins  _et al._ (2010); Walker-Loud (2012), and so using these formulas to
obtain a precise, let alone subpercent, determination, at the physical pion
mass is not possible when incorporating systematic uncertainties associated
with the truncation of $\mathrm{SU}(3)$ HB$\chi$PT.

However, many LQCD calculations, including this one, keep the strange quark
mass fixed near its physical value. Therefore, a simple interpolation in the
strange quark mass is possible. Further, as the omega is an isosinglet, it
will have a simpler, and likely more rapidly converging chiral expansion of
the light-quark mass dependence than baryons with one or more light valence
quarks. This has motivated the construction of an $\mathrm{SU}(2)$ HB$\chi$PT
for hyperons which considers only the pion as a light degree of freedom Beane
_et al._ (2007); Tiburzi and Walker-Loud (2008); Jiang and Tiburzi (2009);
Jiang  _et al._ (2011); Jiang and Tiburzi (2010). In particular, the chiral
expansion for the omega baryon mass was determined to
$\mathrm{O}(m_{\pi}^{6})$ Tiburzi and Walker-Loud (2008)

$\displaystyle m_{\Omega}$
$\displaystyle=m_{0}+{\alpha}_{2}\frac{m_{\pi}^{2}}{{\Lambda}_{\chi}}+\frac{m_{\pi}^{4}}{{\Lambda}_{\chi}^{3}}\left[{\alpha}_{4}{\lambda}_{\pi}+{\beta}_{4}\right]$
$\displaystyle\phantom{=}+\frac{m_{\pi}^{6}}{{\Lambda}_{\chi}^{5}}\left[{\alpha}_{6}{\lambda}_{\pi}^{2}+{\beta}_{6}{\lambda}_{\pi}+{\gamma}_{6}\right]\,,$
(18)

where

${\lambda}_{\pi}={\lambda}_{\pi}(\mu)=\ln\left(\frac{m_{\pi}^{2}}{\mu^{2}}\right)\,,$
(19)

${\alpha}_{n}$, ${\beta}_{n}$ and ${\gamma}_{6}$ are linear combinations of
$\mu$-dependent dimensionless low energy constants (LECs) of the theory, and
$m_{0}$ is the mass of the omega baryon in the $\mathrm{SU}(2)$ chiral limit
at the physical strange quark mass. The renormalization group Weinberg (1979)
restricts the coefficient of the $\ln^{2}$ term to be linearly dependent on
${\alpha}_{2}$ and ${\alpha}_{4}$, with the relation provided in Ref. Tiburzi
and Walker-Loud (2008). In standard HB$\chi$PT power counting, in which the
expansion includes odd powers of the pion mass, this order would be called
next-to-next-to-next-to-next-to-leading order (N4LO), where leading order (LO)
is the $\mathrm{O}(m_{\pi}^{2}/{\Lambda}_{\chi})$ contribution, next-to-
leading order (NLO) would be an $\mathrm{O}(m_{\pi}^{3}/{\Lambda}_{\chi}^{2})$
contribution, which vanishes for $m_{\Omega}$, etc.

The light quark mass dependence for $t_{0}$ and $w_{0}$ has also been
determined in $\chi$PT through $\mathrm{O}(m_{\pi}^{4})$ Bar and Golterman
(2014) which is N2LO in the meson chiral power counting. For example

$w_{0}=w_{0,\rm{ch}}\bigg{\\{}1+k_{1}\frac{m_{\pi}^{2}}{{\Lambda}_{\chi}^{2}}+k_{3}\frac{m_{\pi}^{4}}{{\Lambda}_{\chi}^{4}}+k_{2}\frac{m_{\pi}^{4}}{{\Lambda}_{\chi}^{4}}{\lambda}_{\pi}\bigg{\\}}\,,$
(20)

where the LO term, $w_{0,\rm{ch}}$, is the value in the chiral limit and the
$k_{i}$ are linear combinations of dimensionless LECs. The expression for
$t_{0}$ is identical in form and will have different numerical values of the
LECs.

From these expressions, we can see both $m_{\Omega}$ and $t_{0}$ and $w_{0}$
depend only upon even powers of the pion mass through the order we are
working: $m_{\Omega}$ receives a chiral correction that scales as
$\mathrm{O}(m_{\pi}^{7})$ from a double-sunset two-loop diagram Tiburzi and
Walker-Loud (2008) and the next correction to $t_{0}$ and $w_{0}$ will appear
at $\mathrm{O}(m_{\pi}^{6})$. We can multiply these expressions together, Eqs.
(III.1) and (20), in order to form an expression describing the light-quark
mass dependence of $w_{0}m_{\Omega}$. As the characterization of the order of
the expansion with respect to the order of $m_{\pi}^{2}$ is not the same for
$w_{0}$ and $m_{\Omega}$, we define the contributions to $w_{0}m_{\Omega}$ as

$\displaystyle w_{0}m_{\Omega}$
$\displaystyle=c_{0}+{\delta}_{ls,{\Lambda}}^{\rm
NLO}+{\delta}_{ls,{\Lambda}}^{\textrm{{N${}^{2}$LO}}}+{\delta}_{ls,{\Lambda}}^{\textrm{{N${}^{3}$LO}}}\,,$
(21)

with a similar expression for $\sqrt{t_{0}}m_{\Omega}$. We add polynomial
terms in $s_{{\Lambda},{\Omega}}^{2}$ such that

$\displaystyle{\delta}_{ls,{\Lambda}}^{\rm NLO}$
$\displaystyle=l_{\Lambda}^{2}c_{l}+s_{\Lambda}^{2}c_{s}\,,$
$\displaystyle{\delta}_{ls,{\Lambda}}^{\textrm{{N${}^{2}$LO}}}$
$\displaystyle=l_{\Lambda}^{4}(c_{ll}+c_{ll}^{ln}{\lambda}_{\pi})+l_{\Lambda}^{2}s_{\Lambda}^{2}c_{ls}+s_{\Lambda}^{4}c_{ss}\,,$
$\displaystyle{\delta}_{ls,{\Lambda}}^{\textrm{{N${}^{3}$LO}}}$
$\displaystyle=l_{\Lambda}^{6}(c_{lll}+c_{lll}^{ln}{\lambda}_{\pi}+c_{lll}^{ln^{2}}{\lambda}_{\pi}^{2})+l_{\Lambda}^{4}s_{\Lambda}^{2}{\lambda}_{\pi}c_{lls}^{ln}$
$\displaystyle\phantom{=}+l_{\Lambda}^{4}s_{\Lambda}^{2}c_{lls}+l_{\Lambda}^{2}s_{\Lambda}^{4}c_{lss}+s_{\Lambda}^{6}c_{sss}\,.$
(22)

We will consider both ${\Lambda}=F$ and ${\Lambda}={\Omega}$ for the two
choices of small parameters. For convenience, we set $\mu={\Lambda}_{\chi}$
and $\mu=m_{\Omega}$ respectively for these choices. For a detailed discussion
how one can track the consequence of such a quark mass dependent choice for
the dim-reg scale, see Ref. Miller  _et al._ (2020).

###  III.2 Finite volume corrections

The finite-volume (FV) corrections for $m_{\Omega}$ are determined at one loop
through the modification to the tadpole integral Gasser and Leutwyler (1987);
Colangelo and Haefeli (2004)

$\frac{(4\pi)^{2}}{m^{2}}{\mathcal{I}}^{\rm
FV}=\ln\left(\frac{m^{2}}{\mu^{2}}\right)+4k_{1}(m_{\pi}L)$ (23)

where $k_{1}(x)$ is given by

$\displaystyle k_{1}(x)=\sum_{|\mathbf{n}|\neq
0}c_{n}\frac{K_{1}(x|\mathbf{n}|)}{x|\mathbf{n}|}\,.$ (24)

$K_{1}(x)$ is a modified Bessel function of the second kind and $c_{n}$ are
multiplicity factors for the number of ways the integers $(n_{x},n_{y},n_{z})$
can form a vector of length $|\mathbf{n}|$; see Table 3 for the first few.

Table 3:  Multiplicity factors for the finite volume corrections of the first ten vector lengths, $|\mathbf{n}|$.  $|\mathbf{n}|$ | 1 | $\sqrt{2}$ | $\sqrt{3}$ | $\sqrt{4}$ | $\sqrt{5}$ | $\sqrt{6}$ | $\sqrt{7}$ | $\sqrt{8}$ | $\sqrt{9}$ | $\sqrt{10}$
---|---|---|---|---|---|---|---|---|---|---
$c_{n}$ | 6 | 12 | 8 | 6 | 24 | 24 | 0 | 12 | 30 | 24

At N3LO, the finite volume corrections for $m_{\Omega}$ are also trivially
determined, as the only two-loop integral that contributes is a double-tadpole
with un-nested momentum integrals; see Fig. 2 of Ref. Tiburzi and Walker-Loud
(2008). The N3LO correction to $w_{0}$ is not known. However, the isoscalar
nature of $w_{0}$ means that at the two-loop order, just like the correction
to $m_{\Omega}$, it will only receive contributions from trivial two-loop
integrals with factorizable momentum integrals. Therefore, the N3LO FV
correction can also be determined from the square of the tadpole integral

$\frac{(4\pi)^{4}}{m^{4}}\left[\left({\mathcal{I}}^{\rm
FV}\right)^{2}-\left({\mathcal{I}}^{\infty}\right)^{2}\right]\\\
=8{\lambda}_{\pi}k_{1}(m_{\pi}L)+16k_{1}^{2}(m_{\pi}L)\,,$ (25)

resulting in

$\displaystyle{\delta}_{L,F}^{\textrm{{N${}^{2}$LO}}}(l_{F},m_{\pi}L)$
$\displaystyle=c_{ll}^{ln}l_{F}^{4}\,4k_{1}(m_{\pi}L)$
$\displaystyle{\delta}_{L,F}^{\textrm{{N${}^{3}$LO}}}(l_{F},m_{\pi}L)$
$\displaystyle=c_{lll}^{ln^{2}}l_{F}^{6}\ln(l_{F}^{2})\,8k_{1}(m_{\pi}L)$
$\displaystyle\phantom{=}+c_{lll}^{ln}l_{F}^{6}\,16k_{1}(m_{\pi}L)\,.$ (26)

The FV correction through N3LO arising from loop corrections to $w_{0}$ and
$m_{O}$ is given by

${\delta}_{L,F}={\delta}_{L,F}^{\textrm{{N${}^{2}$LO}}}+{\delta}_{L,F}^{\textrm{{N${}^{3}$LO}}}\,,$
(27)

with a similar expression for ${\delta}_{L,{\Omega}}$.

We are neglecting a few FV corrections at
${\delta}_{L,F}^{\textrm{{N${}^{3}$LO}}}$. The NLO $l_{F,{\Omega}}^{2}$
correction to the omega mass is from a quark mass operator, which has been
converted to an $m_{\pi}^{2}$ correction,
$2B\hat{m}_{l}=m_{\pi}^{2}+\mathrm{O}(m_{\pi}^{4}/{\Lambda}_{\chi}^{2})$. This
choice for organizing the perturbative expansion induces corrections in what
we have called N2LO and N3LO. At N2LO, the corrections arise from single
tadpole diagrams and so the FV corrections are accounted for through Eq.
(III.2). At the next order, the corrections to the pion self-energy involve
more complicated two-loop diagrams Bijnens and Rössler (2015) and so the FV
corrections arising from these are not captured in our parametrization.
Similar corrections arise from expressing $4\pi F_{0}=4\pi
F_{\pi}+\mathrm{O}(m_{\pi}^{2}/{\Lambda}_{\chi}^{2})$ when using $l_{F}^{2}$
to track the light-quark mass corrections [$F_{0}$ is $F_{\pi}$ in the
$\mathrm{SU}(2)$ chiral limit].

While we have neglected these contributions, the FV corrections to
$m_{\Omega}$ are suppressed by an extra power in the chiral power counting
compared with many observables, beginning with an
$m_{\pi}^{4}/{\Lambda}_{\chi}^{4}$ prefactor, Eq. (III.2). In Fig. 5, we show
the predicted FV correction along with the results at three volumes on the
a12m220 ensembles. As can be observed, the predicted FV corrections are very
small and consistent with the numerical results.

Figure 5:  Predicted finite volume corrections (band) compared with the
results on the a12m220 ensembles. The N3LO result is from Eq. (27) using
$l_{F}^{2}$ as the small parameter.

###  III.3 Discretization corrections

A standard method of incorporating discretization effects into the
extrapolation formula used for hadronic observables is to follow the strategy
of Sharpe and Singleton Sharpe and Singleton (1998):

  1. 1.

For a given lattice action, one first constructs the Symanzik effective theory
(SET) by expanding the discretized action about the continuum limit. This
results in a local effective action in terms of quark and gluon fields
Symanzik (1983a, b);

  2. 2.

With this continuum effective theory, one builds a chiral effective theory by
using spurion analysis to construct not only operators with explicit quark
mass dependence, but also operators with explicit lattice spacing dependence.

Such an approach captures the leading discretization effects in a local
hadronic effective theory. At the level of the SET, radiative corrections
generate logarithmic dependence upon the lattice spacing. The leading
corrections can be resummed such that, for an $\mathrm{O}(a)$ improved action,
the leading discretization effects then scale as
$a^{2}{\alpha}_{S}^{n+\hat{\gamma}_{1}}$ Balog  _et al._ (2009, 2010) where
$n=0$ for an otherwise unimproved action, $n=1$ for a tree-level improved
action and $n=2$ for a one-loop improved action. The coefficient
$\hat{\gamma}_{1}$ is an anomalous dimension which has been determined for
Yang-Mills and Wilson actions Husung  _et al._ (2020).222See also the
presentation by N. Husung at the MIT Virtual Lattice Field Theory Colloquium
Series, http://ctp.lns.mit.edu/latticecolloq/.

For mixed-action setups Renner  _et al._ (2005); Bar  _et al._ (2003) such as
the one used in this work, a low-energy mixed-action effective field theory
(MAEFT) Bar  _et al._ (2003, 2004, 2005); Golterman  _et al._ (2005); Tiburzi
(2005); Chen  _et al._ (2006, 2007); Orginos and Walker-Loud (2008); Chen  _et
al._ (2009a, b) can be constructed to capture the manifestation of infrared
radiative corrections from the discretization333While this might seem
counterintuitive, it is analogous to the infrared sensitivity of hadronic
quantities to the Higgs vacuum expectation value (vev): hadronic quantities
have infrared (logarithmic) sensitivity to the pion mass from radiative pion
loops, and the squared pion mass is proportional to the light quark mass which
is proportional to the Higgs vev hig (2005).. Corrections come predominantly
from a modification of the pseudoscalar meson spectrum as well as from
“hairpin” interactions Bernard and Golterman (1992) that are proportional to
the lattice spacing in rooted-staggered Bernard and Golterman (1994) and
mixed-action theories Golterman  _et al._ (2005); in partially quenched
theories, these hairpins are proportional to the difference in the valence and
sea quark masses Sharpe (1997); Sharpe and Shoresh (2000, 2001).

In our analysis of $F_{K}/F_{\pi}$, we observed that the use of continuum
chiral perturbation theory with corrections polynomial in ${\epsilon}_{a}^{2}$
was highly favored over the use of the MAEFT expression, as measured by the
Bayes-Factor, though the results from both were consistent within a fraction
of 1 standard deviation Miller  _et al._ (2020). Similar findings have been
observed by other groups for various quantities; see for example Refs.
Borsanyi  _et al._ (2013); Colquhoun  _et al._ (2016); Dowdall  _et al._
(2019). Therefore, in this work, we restrict our analysis to a continuum-like
expression enhanced by polynomial discretizaton terms.

The dynamical HISQ ensembles have a perturbatively improved action such that
the leading discretization effects (before resumming the radiative corrections
Balog  _et al._ (2009, 2010); Husung  _et al._ (2020)) scale as
$\mathrm{O}({\alpha}_{S}a^{2})$ Follana  _et al._ (2007). The MDWF action, in
the limit of infinite extent in the fifth dimension, has no chiral symmetry
breaking other than that from the quark mass. Consequently, the leading
discretization corrections begin at $\mathrm{O}(a^{2})$ Sharpe (2007); Aoki
_et al._ (2011). For finite $L_{5}$, the $\mathrm{O}(a)$ corrections are
proportional to $am_{\rm res}$ which is sufficiently small that these terms
are numerically negligible. Therefore, we parametrize our discretizaton
corrections with the following terms where we count ${\epsilon}_{a}^{2}\sim
l_{\Lambda}^{2}\sim s_{\Lambda}^{2}$:

$\displaystyle{\delta}_{a,{\Lambda}}$
$\displaystyle={\delta}_{a,{\Lambda}}^{\rm
NLO}+{\delta}_{a,{\Lambda}}^{\textrm{{N${}^{2}$LO}}}+{\delta}_{a,{\Lambda}}^{\textrm{{N${}^{3}$LO}}}\,,$
$\displaystyle{\delta}_{a,{\Lambda}}^{\rm NLO}$
$\displaystyle=d_{a}{\epsilon}_{a}^{2}+d^{\prime}_{a}{\alpha}_{S}{\epsilon}_{a}^{2}\,,$
$\displaystyle{\delta}_{a,{\Lambda}}^{\textrm{{N${}^{2}$LO}}}$
$\displaystyle=d_{aa}{\epsilon}_{a}^{4}+{\epsilon}_{a}^{2}\left(d_{al}l_{\Lambda}^{2}+d_{as}s_{\Lambda}^{2}\right)\,,$
$\displaystyle{\delta}_{a,{\Lambda}}^{\textrm{{N${}^{3}$LO}}}$
$\displaystyle=d_{aaa}{\epsilon}_{a}^{6}+{\epsilon}_{a}^{4}(d_{aal}l_{\Lambda}^{2}+d_{aas}s_{\Lambda}^{2})$
$\displaystyle\phantom{=}\quad+{\epsilon}_{a}^{2}(d_{all}l_{\Lambda}^{4}+d_{als}l_{\Lambda}^{2}s_{\Lambda}^{2}+d_{ass}s_{\Lambda}^{4})\,.$
(28)

##  IV Extrapolation Details and Uncertainty Analysis

We perform our extrapolation analysis under a Bayesian model-averaging
framework as described in detail in Refs. Chang  _et al._ (2018); Miller  _et
al._ (2020); Kass and Raftery (1995), which is more extensively discussed for
lattice QFT analysis in Ref. Jay and Neil (2020). We consider a variety of
extrapolation functions by working to different orders in the power counting,
using the $l_{F},s_{F}$ or $l_{\Omega},s_{\Omega}$ small parameters, by
including or excluding the chiral logarithms associated with pion loops, and
by including or excluding discretization corrections scaling as
${\alpha}_{S}a^{2}$. The resulting Bayes factors are then used to weight the
fits with respect to each other and perform a model averaging. In this
section, we discuss the selection of the priors for the various LECs and then
present an uncertainty analysis of the results.

###  IV.1 Prior widths of LECs

In our $F_{K}/F_{\pi}$ analysis Miller  _et al._ (2020), we observed that
using ${\epsilon}_{\pi}^{2}=l_{F}^{2}$, ${\epsilon}_{K}^{2}$ and
${\epsilon}_{a}^{2}$ as the small parameters in the expansion,444The small
parameter
${\epsilon}_{K}^{2}=\frac{m_{K}^{2}}{{\Lambda}_{\chi}^{2}}=\frac{1}{2}(s_{F}^{2}+l_{F}^{2})$.
the LECs were naturally of $\mathrm{O}(1)$. We therefore have a prior
expectation that this may hold for $\sqrt{t_{0}}m_{\Omega}$ and
$w_{0}m_{\Omega}$ as well.

Let us use $w_{0,\rm orig}m_{\Omega}$ to guide the discussion. The $a\sim
0.12$ fm ensembles were simulated with a fixed strange quark mass. Therefore,
the entire change in $w_{0,\rm orig}m_{\Omega}$ between the a12m130 and
a12m180L ensembles can be attributed to the change in $l_{F}$. This allows us
to “eyeball” the $c_{l}$ prior to be $c_{l}\simeq 1$ if we assume the dominant
contribution comes from the NLO $c_{l}l_{F}^{2}$ term.555This is analogous to
using the effective mass and effective overlap factors to choose conservative
priors for the ground state parameters in the correlation function analysis.
Motivated by $\mathrm{SU}(3)$ flavor symmetry considerations, we can roughly
expect $c_{s}\sim c_{l}$. In order to be conservative, we set the prior for
these LECs as

$\tilde{c}_{l}=\tilde{c}_{s}=N(\mu=1,\sigma=1)\,,$ (29)

where $N(\mu,\sigma)$ denotes a normal distribution with mean $\mu$ and width
$\sigma$. A similar observation is made for $w_{0,\rm imp}m_{\Omega}$ and the
original and improved values of $\sqrt{t_{0}}m_{\Omega}$. We observe (with a
full analysis) that the log-Bayes-Factor (logGBF) prefers even tighter priors,
with logGBF continuing to increase as the width is taken down to 0.1 on these
NLO LECs.666For fixed data, $\exp\\{\rm logGBF\\}$ provides a relative weight
of the likelihood of one model versus another.

The observation that $m_{\Omega}$ increases with increasing values of $l_{F}$
and $s_{F}$ (normalized by any and all gradient flow scales) allows us to
conservatively estimate the LO prior,

$\tilde{c}_{0}=N(1,1)\,.$ (30)

We then conservatively estimate the priors for all of the higher order $l_{F}$
and $s_{F}$ LECs to be

$\tilde{c}_{i}=N(0,1)\,.$ (31)

We observe, with a full analysis, that this choice is near the optimal value
as measured by the logGBF weighting.

For the discretization corrections, see Fig. 7 in Sec. IV.2, as we change the
gradient flow scale from $w_{0,\rm orig}$ to the improved version to using the
original and improved versions of $\sqrt{t_{0}}$, the approach to the
continuum limit can change sign. We also observe, the convexity of the
approach to the continuum limit (the ${\epsilon}_{a}^{4}$ contributions) can
change sign. Therefore, we perform a prior-optimization study for the
discretization LECs in which we change the prior width of the NLO and N2LO
LECs in concert, such that the priors are given by

$\tilde{d}_{i}=N(0,{\sigma}_{a}^{\rm opt})\,,$ (32)

with the optimized values

${\sigma}_{a}^{\rm opt}=\left\\{\begin{array}[]{cl}1.2,&w_{0,\rm
orig}m_{\Omega}\\\ 1.4,&w_{0,\rm imp}m_{\Omega}\\\ 1.8,&\sqrt{t_{0,\rm
orig}}m_{\Omega}\\\ 1.4,&\sqrt{t_{0,\rm imp}}m_{\Omega}\\\
\end{array}\right.\,,$ (33)

maximizing the respective logGBF factors. We set the N3LO priors for
discretization LECs to

$\tilde{d}_{i}=N(0,1)\,.$ (34)

Although we find using tighter priors increases the logGBF, the final results
are unchanged.

When we add the ${\alpha}_{S}{\epsilon}_{a}^{2}$ term in the analysis, this
introduces a fourth class of discretization corrections. As we only have four
lattice spacings in this work, we perform an independent prior-width
optimization for this LEC. For all four choices of gradient-flow scales, we
find the choice,

$\tilde{d}_{a}^{\prime}=N(0,0.5)\,,$ (35)

to be near-optimal, with three of the analyses preferring an even tighter
prior width. An empirical Bayes study Lepage  _et al._ (2002) in which the
widths of all the chiral and all the discretization priors are varied together
at a given order leads to similar choices of all the priors.

In Table 4 we list the values of all the priors used in the final analysis.
The full analysis demonstrates that these choices also result in no tension
between the priors and the final posterior values of the LECs, further
indication that our choices are reasonable.

When we use $l_{\Omega}^{2}$ and $s_{\Omega}^{2}$ as the small parameters
instead of $l_{F}^{2}$ and $s_{F}^{2}$, we note that since
$(m_{\Omega}/{\Lambda}_{\chi})^{2}\sim 2$, we can use the same prescription,
except to double the mean and width of all the NLO priors (which scale
linearly in $l_{\Omega}^{2}$ and $s_{\Omega}^{2}$), set the widths to be 4
times larger for the N2LO priors and 8 times larger for the N3LO LECs. The
mixed contributions which scale with some power of ${\epsilon}_{a}^{2}$ and
$l_{\Omega}^{2}$ and $s_{\Omega}^{2}$ are scaled accordingly; see also Table
4.

Table 4:  The values of the priors used in the analysis.  LEC | ${\Lambda}=F$ | ${\Lambda}={\Omega}$
---|---|---
$c_{0}$ | $N(1,1)$ | $N(1,1)$
$c_{l},c_{s}$ | $N(1,1)$ | $N(2,2)$
$c_{ll},c_{ll}^{ln},c_{ls},c_{ss}$ | $N(0,1)$ | $N(0,4)$
$c_{lll},c_{lll}^{ln},c_{lll}^{ln^{2}},c_{lls},c_{lls}^{ln},c_{lss},c_{sss}$ | $N(0,1)$ | $N(0,8)$
$d_{a},d_{aa}$ | $N(0,{\sigma}_{a}^{\rm opt})$ | $N(0,{\sigma}_{a}^{\rm opt})$
$d^{\prime}_{a}$ | $N(0,0.5)$ | $N(0,0.5)$
$d_{al},d_{as}$ | $N(0,{\sigma}_{a}^{\rm opt})$ | $N(0,2{\sigma}_{a}^{\rm opt})$
$d_{aaa}$ | $N(0,1)$ | $N(0,1)$
$d_{aal},d_{aas}$ | $N(0,1)$ | $N(0,2)$
$d_{all},d_{als},d_{ass}$ | $N(0,1)$ | $N(0,4)$

###  IV.2 Extrapolation analysis

For each of the four quantities, $w_{0,\rm orig}m_{\Omega}$, $w_{0,\rm
imp}m_{\Omega}$, $\sqrt{t_{0,\rm orig}}m_{\Omega}$ and $\sqrt{t_{0,\rm
imp}}m_{\Omega}$, we consider several reasonable choices of extrapolation
functions to perform the continuum, infinite volume and physical quark-mass
limits. The final result for each extrapolation is then determined through a
model average in which the relative weight of each model is given by the
exponential of the corresponding logGBF value. The various choices we consider
in the extrapolations consists of

$\displaystyle\begin{array}[]{rl}\text{Include the $\ln(m_{\pi})$ terms or
counterterm only}:&\times 2\\\ \text{Expand to {N${}^{2}$LO}, or
{N${}^{3}$LO}}:&\times 2\\\ \text{Include/exclude finite volume
corrections}:&\times 2\\\ \text{Include/exclude the $\alpha_{S}a^{2}$
term}:&\times 2\\\ \text{Use the ${\Lambda}=F$ or ${\Lambda}={\Omega}$
expansion}:&\times 2\\\ \hline\cr\text{total choices}:&32\end{array}$

We find that there is very little dependence upon the particular model chosen.
In Fig. 6, we show the stability of the final result of $\sqrt{t_{0,\rm
imp}}m_{\Omega}$ and $w_{0,\rm imp}m_{\Omega}$ as various options from the
above list are turned on and off. In addition, we show the impact of including
or excluding the a12m220ms ensemble, whose strange quark mass is $m_{s}\sim
0.6\times m_{s}^{\rm phys}$, as well as the impact of including the a06m310L
ensemble. We observe a small variation of the result when either of these
ensembles is dropped, but the results are still consistent with our final
result (top of the figure). Using the fixed definition of
${\epsilon}_{a}^{2}$, Eq. (16), we find

$\displaystyle\sqrt{t_{0,\rm orig}}m_{\Omega}$
$\displaystyle={1.2107(78)^{s}(14)^{\chi}(40)^{a}(00)^{V}(22)^{\rm
phys}(07)^{M}}\rightarrow\frac{\sqrt{t_{0}}}{\rm
fm}={0.1429(09)^{s}(02)^{\chi}(05)^{a}(00)^{V}(03)^{\rm phys}(01)^{M}}\,,$
$\displaystyle\sqrt{t_{0,\rm imp}}m_{\Omega}$
$\displaystyle={1.2052(73)^{s}(13)^{\chi}(40)^{a}(00)^{V}(22)^{\rm
phys}(08)^{M}}\rightarrow\frac{\sqrt{t_{0}}}{\rm
fm}={0.1422(09)^{s}(02)^{\chi}(05)^{a}(00)^{V}(03)^{\rm phys}(01)^{M}}\,,$
$\displaystyle w_{0,\rm orig}m_{\Omega}$
$\displaystyle={1.4481(81)^{s}(14)^{\chi}(46)^{a}(00)^{V}(26)^{\rm
phys}(11)^{M}}\rightarrow\frac{w_{0}}{\rm
fm}={}0.1709(10)^{s}(02)^{\chi}(05)^{a}(00)^{V}(03)^{\rm phys}(01)^{M}\,,$
$\displaystyle w_{0,\rm imp}m_{\Omega}$
$\displaystyle={1.4467(83)^{s}(15)^{\chi}(44)^{a}(00)^{V}(24)^{\rm
phys}(10)^{M}}\rightarrow\frac{w_{0}}{\rm
fm}={0.1707(10)^{s}(02)^{\chi}(05)^{a}(00)^{V}(03)^{\rm phys}(01)^{M}}\,,$
(36)

with the statistical ($s$), chiral interpolation ($\chi$), continuum-limit
($a$), infinite-volume ($V$), physical-point (phys), and model selection
uncertainties ($M$). The conversion to physical units is performed with Eq.
(III.0.1).

Figure 6:  Model breakdown and comparison using the scale-improved, gradient-
flow derived quantities (i.e., $\sqrt{t_{0,\text{imp}}}m_{\Omega}$,
$w_{0,\text{imp}}m_{\Omega}$). The vertical band is our model average result
to guide the eye. $\chi$pt-full: $\chi$PT model average, including
$\ln(m_{\pi}^{2}/\mu^{2})$ corrections. $\chi$pt-ct: $\chi$PT model average
with counterterms only, excluding $\ln(m_{\pi}^{2}/\mu^{2})$ corrections.
N3LO/N2LO: model average restricted to specified order. $\Lambda_{\chi}$:
model average with specified chiral cutoff. incl./excl. $\alpha_{S}$: model
average with/without $\alpha_{S}$ corrections. fixed/variable
$\epsilon_{a}^{2}$: model average with Eqs. (16) or (17), respectively. excl.
a06m310L: model average excluding $a=0.06$ fm ensemble (a06m310L). excl.
a12m220ms: model average excluding small strange quark mass ensemble
(a12m220ms). Below solid line: results from other collaborations for various
numbers of dynamical fermions: BMWc [2020] Borsanyi  _et al._ (2021), MILC
[2015] Bazavov  _et al._ (2016), HPQCD [2013] Dowdall  _et al._ (2013), CLS
[2017] Bruno  _et al._ (2017), QCDSF-UKQCD [2015] Bornyakov  _et al._ (2015),
RBC [2014] Blum  _et al._ (2016), HotQCD [2014] Bazavov  _et al._ (2014a),
BMWc [2012] Deuzeman and Wenger (2012) and ALPHA [2013] Bruno and Sommer
(2014).

Figure 7:  Example extrapolations of $\sqrt{t_{0,\rm orig}}m_{\Omega}$ and
$w_{0,\rm orig}m_{\Omega}$ versus $l_{F}^{2}$ (top) and ${\epsilon}_{a}^{2}$
(bottom) using the N3LO $\chi$PT analysis. For the continuum extrapolation
plots (bottom) we also show the results using the “improved” determinations of
the scales, Eq. (II.3). The numerical results have been shifted from the
original data points as described in the text.

As discussed in Sec. III.0.3, we explore potential systematics in the use of
one definition of the small parameter used to characterize the discretization
corrections, Eq. (16), versus another equally valid choice, Eq. (17). For each
choice of $w_{0}m_{\Omega}$ and $\sqrt{t_{0}}m_{\Omega}$ that is extrapolated
to the physical point, we repeat the above model averaging procedure, but we
also include the variation of using these two definitions of
${\epsilon}_{a}^{2}$ for which we find

$\displaystyle\sqrt{t_{0,\rm orig}}m_{\Omega}$
$\displaystyle={1.2050(87)^{s}(16)^{\chi}(46)^{a}(00)^{V}(21)^{\rm
phys}(61)^{M}}\rightarrow\frac{\sqrt{t_{0}}}{\rm
fm}={0.1422(10)^{s}(02)^{\chi}(05)^{a}(00)^{V}(02)^{\rm phys}(07)^{M}}\,,$
$\displaystyle\sqrt{t_{0,\rm imp}}m_{\Omega}$
$\displaystyle={1.2052(77)^{s}(14)^{\chi}(46)^{a}(00)^{V}(21)^{\rm
phys}(10)^{M}}\rightarrow\frac{\sqrt{t_{0}}}{\rm
fm}={0.1422(09)^{s}(02)^{\chi}(05)^{a}(00)^{V}(02)^{\rm phys}(01)^{M}}\,,$
$\displaystyle w_{0,\rm orig}m_{\Omega}$
$\displaystyle={1.4481(81)^{s}(14)^{\chi}(46)^{a}(00)^{V}(26)^{\rm
phys}(11)^{M}}\rightarrow\frac{w_{0}}{\rm
fm}={}0.1709(10)^{s}(02)^{\chi}(05)^{a}(00)^{V}(03)^{\rm phys}(01)^{M}\,,$
$\displaystyle w_{0,\rm imp}m_{\Omega}$
$\displaystyle={1.4485(82)^{s}(15)^{\chi}(44)^{a}(00)^{V}(25)^{\rm
phys}(18)^{M}}\rightarrow\frac{w_{0}}{\rm
fm}={0.1709(10)^{s}(02)^{\chi}(05)^{a}(00)^{V}(03)^{\rm phys}(02)^{M}}\,.$
(37)

For all choices of the gradient-flow scale besides $\sqrt{t_{0,\rm orig}}$,
the average over the choice of how to define ${\epsilon}_{\alpha}^{2}$ has
minimal impact on the final result, as can be seen comparing Eq. (IV.2) and
(LABEL:eq:final_t0w0_avg_eps2a). In the case of $\sqrt{t_{0,\rm orig}}$, the
two choices for the cutoff-effect expansion parameter lead to a slight
difference in the continuum-extrapolated value, which is reflected in the
model-averaging uncertainty. 777Rather than performing a model average over
the two definitions of $\epsilon^{2}_{a}$ as defined in Eqs. (16) and (17),
one might instead consider a model average over the choices for fixed
$\epsilon^{2}_{a}$, i.e.
$\epsilon^{2}_{a}=\left\\{\left(\frac{a}{2w_{0,\text{orig}}}\right)^{2}\,,\left(\frac{a}{2w_{0,\text{impr}}}\right)^{2}\,,\frac{a^{2}}{4t_{0,\text{orig}}}\,,\frac{a^{2}}{4t_{0,\text{impr}}}\right\\}\,.$
Performing the model average in this manner instead yields results for
$\sqrt{t_{0}},w_{0}$ consistent within a fraction of a sigma of Eqs. (I) and
(I).  In all cases, the dominant uncertainty is statistical, suggesting a
straightforward path to reducing the uncertainty to a few per-mille.

To arrive at our final determination of $\sqrt{t_{0}}$, Eq. (I) and $w_{0}$,
Eq. (I), we perform an average of the results in Eq.
(LABEL:eq:final_t0w0_avg_eps2a). As the data between the two choices differ
slightly, we can not perform this final averaging step under the Bayes model-
averaging procedure; instead we treat each result with equal weight. We would
add half the difference between the central values as an additional
discretization uncertainty, but as is evident from Eq.
(LABEL:eq:final_t0w0_avg_eps2a), the central values are essentially the same.

In Fig. 6, we also compare our result with other values in the literature. All
the results, except the most recent one from BMWc Borsanyi  _et al._ (2021),
have been determined in the isospin symmetric limit. Our results are in good
agreement with the more recent and precise results, though one notes, there is
some tension in the values of $\sqrt{t_{0}}$ and $w_{0}$ reported.

In Fig. 7, we show the resulting extrapolation of $\sqrt{t_{0,\rm
orig}}m_{\Omega}$ and $w_{0,\rm orig}m_{\Omega}$ projected into the
$l_{F}^{2}$ plane using the N3LO analysis including the $\ln(m_{\pi})$ type
corrections. The finite lattice spacing bands are plotted with a value of
${\epsilon}_{a}^{2}$ taken from the near-physical pion mass ensembles from
Table 2, a15m135XL, a12m130 and a09m135. For the $a\sim 0.06$ fm band, we use
the value,

$\frac{w_{0,{\rm orig}}}{a_{06}}=3.0119(19)\,,$ (38)

from Table IV of Ref. Bazavov  _et al._ (2016) with
$m^{\prime}_{l}/m^{\prime}_{s}=1/27$ and use this to construct
${\epsilon}_{a}^{2}$. The data points are plotted after being shifted to the
extrapolated values of all the parameters using the posterior values of the
LECs from the N3LO fit.

The lower panel of Fig. 7 is similarly constructed by shifting all the data
points to the infinite volume limit, $l_{F}^{\rm phys}$, $s_{F}^{\rm phys}$
and the value of ${\epsilon}_{a}^{2}$ from the particular ensemble with the
corresponding band in this same limit and only varying ${\epsilon}_{a}^{2}$.
We plot the continuum extrapolation of both the original and improved values
to demonstrate the impact of the improvement at finite lattice spacing, noting
the agreement in the continuum limit.

For $w_{0}$, there is very little difference between the original and improved
values with very similar continuum extrapolations. In contrast, there is a
striking difference between the original and improved values using
$\sqrt{t_{0}}$, though they agree in the continuum limit. We also observe that
the use of $\sqrt{t_{0,\rm orig}}$ is susceptible to larger model-
extrapolation uncertainties arising from different choices of parametrizing
the continuum extrapolation; see Eq. (LABEL:eq:final_t0w0_avg_eps2a).
Additional results at $a\lesssim 0.06$ fm will be required to control the
continuum extrapolation using $\sqrt{t_{0}}$ in order to obtain a few-per-
mille level of precision.

###  IV.3 Interpolation of $t_{0}$ and $w_{0}$

With our determination of $t_{0}$ and $w_{0}$, Eqs. (I) and (I), we can
determine the lattice spacing for each bare coupling. We could use the near-
physical pion mass ensemble values of the gradient-flow scales, or
alternatively, we could interpolate the results to the physical quark mass
point using the predicted quark mass dependence Bar and Golterman (2014). The
interpolation can be performed for each lattice spacing separately, or in a
combined analysis of all lattice spacings simultaneously. The latter is
preferable in order for us to determine the lattice spacing $a_{06}$ as we
only have results at a single pion mass at this lattice spacing. To perform
the global analysis, we use an N2LO extrapolation function (which has the same
form for $t_{0}$),

$\displaystyle\frac{w_{0}}{a}$ $\displaystyle=\frac{w_{0,{\rm
ch}}}{a}\bigg{\\{}1+k_{l}l_{F}^{2}+k_{s}s_{F}^{2}+k_{a}{\epsilon}_{a,{\rm
ch}}^{2}$
$\displaystyle\qquad+k_{ll}l_{F}^{4}+k_{lln}l_{F}^{4}\ln(l_{F}^{2})+k_{ls}l_{F}^{2}s_{F}^{2}+k_{ss}s_{F}^{4}$
$\displaystyle\qquad+k_{aa}{\epsilon}_{a,{\rm
ch}}^{4}+k_{al}l_{F}^{2}{\epsilon}_{a,{\rm
ch}}^{2}+k_{as}s_{F}^{2}{\epsilon}_{a,{\rm ch}}^{2}\bigg{\\}}\,,$ (39)
$\displaystyle{\epsilon}_{a,{\rm ch}}$
$\displaystyle\equiv\frac{1}{(2w_{0,{\rm ch}}/a)}\,.$ (40)

This global analysis treats the value of LO parameter $\frac{w_{0,{\rm
ch}}}{a}$ for each lattice spacing as a separate unknown parameter, and then
assumes that the remaining dimensionless LECs are shared between all lattice
spacings. We use this LO parameter to also construct ${\epsilon}_{a,{\rm ch}}$
which controls the discretization corrections rather than using
${\epsilon}_{a}$, as ${\epsilon}_{a}$ is half the inverse of the left-hand
side of Eq. (IV.3).

It is tempting to think of this as a combined chiral and continuum
extrapolation analysis of $w_{0}$, but it is not as one normally thinks of
them. Because we do not know the lattice spacings already, there remains an
ambiguity in the interpretation of $w_{0,\rm ch}/a$ and the LECs $k_{a}$,
$k_{aa}$, $k_{la}$ and $k_{sa}$: we are not able to interpret $w_{0,\rm ch}/a$
as the chiral limit value of $w_{0}$ divided by the lattice spacing. We
perform this analysis for all four gradient-flow scales with independent LECs
for each scale as well as a similar parametrization of the ${\epsilon}_{a,{\rm
ch}}$ parameter.

When we perform the interpolation for each lattice spacing separately, we
utilize this same expression except that we set all parameters proportional to
any power of ${\epsilon}_{a,{\rm ch}}$ to zero. When the individual
interpolations are used, the resulting values of $a_{15}$, $a_{12}$ and
$a_{09}$ are compatible with those from the global analysis well within 1
standard deviation. In Fig. 8, we show an example interpolation of $w_{0,\rm
orig}/a$ using the global analysis of all ensembles to the $l_{F}^{\rm phys}$
and $s_{F}^{\rm phys}$ point for each lattice spacing. The open circles show
the raw values of $w_{0}/a$ while the open squares show the values shifted to
$\frac{w_{0}}{a}(l_{F}^{\rm ens},s_{F}^{\rm phys},{\epsilon}_{a}^{\rm ens})$
using the resultant parameters determined in the global analysis.

Figure 8:  Interpolation of $w_{0}/a$ values to the physical $l_{F}$ and
$s_{F}$ point. The open circles are the raw data from Table 2 and the open
squares have been shifted to the infinite volume and $s_{F}^{\rm phys}$ point.
The vertical gray line is at $l_{F}^{2,\rm phys}$.

A quark-mass independent determination of the lattice spacing can be made by
using the determination of $\sqrt{t_{0}}$, Eq. (I) or $w_{0}$, Eq. (I) at the
physical point, combined with the physical-quark mass interpolated values of
$t_{0}/a^{2}$ or $w_{0}/a$ from either the original or improved values of
these gradient flow scales. Each choice represents a different scheme for
setting the lattice spacing. The continuum extrapolated value of some
observable quantity, using any of these schemes, should agree in the continuum
limit, while at finite lattice spacing, the results can be substantially
different, as is evident in Fig. 7.

In Table 5, we provide the determination of the lattice spacing for each bare
coupling, expressed in terms of the approximate lattice spacing. It is
interesting to note that the determination of the lattice spacing with
$\sqrt{t_{0,\rm orig}}$ is substantially different than with the other
gradient-flow scales, while the scale determined with from the three remaining
auxiliary scales are very similar.

Table 5:  The value of the lattice spacing in fm using the interpolation of the four different determinations of the gradient-flow scales, combined with the determination of $\sqrt{t_{0}}$ and $w_{0}$ at the physical point, Eq. (I) and Eq. (I) respectively. Any one scheme for determining the lattice spacing can be picked as a quark-mass independent scale setting that can be used to convert results from dimensionless lattice units to MeV. The different schemes will result in different approaches to the continuum; however all choices should agree in the continuum limit.  Scheme |  $a_{15}$/fm |  $a_{12}$/fm |  $a_{09}$/fm |  $a_{06}$/fm
---|---|---|---|---
$t_{0,\rm orig}/a^{2}$ | 0.1284(10) | 0.10788(83) | 0.08196(64) | 0.05564(44)
$t_{0,\rm imp}/a^{2}$ | 0.1428(10) | 0.11735(87) | 0.08632(65) | 0.05693(44)
$w_{0,\rm orig}/a$ | 0.1492(10) | 0.12126(87) | 0.08789(71) | 0.05717(51)
$w_{0,\rm imp}/a$ | 0.1505(10) | 0.12066(88) | 0.08730(70) | 0.05691(51)

##  V Summary and Discussion

We have performed a precise scale setting with our MDWF on gradient-flowed
HISQ action Berkowitz  _et al._ (2017a) achieving a total uncertainty of $\sim
0.6\%$–$0.8\%$ for each lattice spacing, Table 5. The scale setting was
performed by extrapolating the quantities
$\sqrt{t_{0}}m_{\Omega}(l_{F},s_{F},{\epsilon}_{a},m_{\pi}L)$ and
$w_{0}m_{\Omega}(l_{F},s_{F},{\epsilon}_{a},m_{\pi}L)$ to the continuum
(${\epsilon}_{a}\rightarrow 0$), infinite volume ($m_{\pi}L\rightarrow\infty$)
and physical quark mass limits ($l_{F}\rightarrow l_{F}^{\rm phys}$ and
$s_{F}\rightarrow s_{F}^{\rm phys}$), and using the experimental determination
of $m_{\Omega}$ to determine the scales $\sqrt{t_{0}}$ and $w_{0}$ in fm. The
values of $\sqrt{t_{0,\rm orig}}/a$, $\sqrt{t_{0,\rm imp}}/a$, $w_{0,\rm
orig}/a$ and $w_{0,\rm imp}/a$ were interpolated to the infinite volume and
physical quark mass limits for each lattice spacing, allowing for the quark-
mass independent determination of $a$ for each bare coupling $\beta$,
expressed in terms of the approximate lattice spacing; see Table 5.

Of note, the approach to the continuum limit of $\sqrt{t_{0,\rm
orig}}m_{\Omega}$ and $\sqrt{t_{0,\rm imp}}m_{\Omega}$ are quite different,
Fig. 7, with the use of $\sqrt{t_{0,\rm imp}}$ leading to an almost flat
continuum extrapolation. The two different extrapolations agree quite nicely
in the continuum limit, as they must if all systematic uncertainties are under
control. In contrast, the use of the original and improved values of $w_{0}$
leads to very similar continuum extrapolations of $w_{0,\rm orig}m_{\Omega}$
and $w_{0,\rm imp}m_{\Omega}$, which also agree very nicely in the continuum
limit.

We also observe that the use of $l_{\Omega}$ and $s_{\Omega}$ as small
parameters to control the quark-mass interpolation are relatively heavily
penalized as compared to the use of $l_{F}$ and $s_{F}$; see Table 7 for an
example. We observe the same qualitative weighting with all choices of the
gradient-flow scale. Perhaps this is an indication that this parametrization
is suboptimal.

Our final uncertainty using $w_{0}$ is dominated by the stochastic
uncertainty, Eq. (I), providing a clear path to reducing the uncertainty by
almost a factor of 3 before an improved understanding of the various
systematic uncertainties becomes relevant. At such a level of precision, a
systematic study of the effect of isospin breaking on the scale setting, as
has been performed by BMWc Borsanyi  _et al._ (2021), is likely required to
retain full control of the uncertainty. For $\sqrt{t_{0}}$, we observe the
model-selection uncertainty is comparable to the stochastic uncertainty, Eq.
(I), which arises from the different ways to parameterize the continuum
extrapolation; see Eq. (LABEL:eq:final_t0w0_avg_eps2a). Therefore, additional
results at $a\lesssim 0.06$ fm will be required to obtain a fer-per-mille
precision with $\sqrt{t_{0}}$.

The pursuit of our physics program of determining the nucleon elastic
structure functions and improving the precision of our $g_{A}$ result Chang
_et al._ (2018); Berkowitz  _et al._ (2018a) will naturally lead to an
improved scale setting precision. The current precision is already expected to
be subdominant for most of the results we will obtain, but a further improved
precision is welcome.

###### Acknowledgements.

We thank M. Goltermann and O. Bär for useful discussion and correspondence
about the chiral corrections to $w_{0}$. We thank Andrea Shindler for helpful
discussions regarding reweighting. We thank Peter Lepage for helpful
correspondence regarding lsqfit and the log Gaussian Bayes factor. We thank K.
Orginos for use of wm_chroma that was used to compute some of the correlation
functions used in this work. We thank the MILC Collaboration for providing
some of the HISQ configurations used in this work, and A. Bazavov, C. Detar
and D. Toussaint for guidance on using their code to generate the new HISQ
ensembles also used in this work. We thank R. Sommer for helpful
correspondence and encouragement to determine $t_{0}$ as well as $w_{0}$.
Computing time for this work was provided through the Innovative and Novel
Computational Impact on Theory and Experiment (INCITE) program and the LLNL
Multiprogrammatic and Institutional Computing program for Grand Challenge
allocations on the LLNL supercomputers. This research utilized the NVIDIA GPU-
accelerated Titan and Summit supercomputers at Oak Ridge Leadership Computing
Facility at the Oak Ridge National Laboratory, which is supported by the
Office of Science of the U.S. Department of Energy under Contract No. DE-
AC05-00OR22725 as well as the Surface, RZHasGPU, Pascal, Lassen, and Sierra
supercomputers at Lawrence Livermore National Laboratory. The computations
were performed utilizing LALIBE lal  which utilizes the Chroma software suite
Edwards and Joo (2005) with QUDA solvers Clark  _et al._ (2010); Babich  _et
al._ (2011) and HDF5 The HDF Group (1997-NNNN) for I/O Kurth  _et al._ (2015).
They were efficiently managed with METAQ Berkowitz (2017); Berkowitz  _et al._
(2018b) and status of tasks logged with EspressoDB Chang  _et al._ (2020). The
HMC was performed with the MILC Code mil , and for the ensembles new in this
work, running on GPUs using QUDA. The final extrapolation analysis utilized
gvar v11.2 Lepage (2020a) and lsqfit v11.5.1 Lepage (2020b). The analysis and
data for this work can be found at this git repo: https://github.com/callat-
qcd/project_scale_setting_mdwf_hisq. This work was supported by the NVIDIA
Corporation (MAC), the Alexander von Humboldt Foundation through a Feodor
Lynen Research Fellowship (CK), the DFG and the NSFC Sino-German CRC110 (EB),
the RIKEN Special Postdoctoral Researcher Program (ER), the U.S. Department of
Energy, Office of Science, Office of Nuclear Physics under Awards No. DE-
AC02-05CH11231 (CCC, CK, BH, AWL), No. DE-AC52-07NA27344 (DAB, DH, ASG, PV),
No. DE-FG02-93ER-40762 (EB), No. DE-AC05-06OR23177 (CM); the Nuclear Physics
Double Beta Decay Topical Collaboration (DAB, HMC, AN, AWL); the U.K. Science
and Technology Facilities Council Grants No. ST/S005781/1 and No. ST/T000945/1
(CB); and the DOE Early Career Award Program (CCC, AWL).

##  Appendix A CHARM QUARK MASS REWEIGHTING

The use of reweighting Ferrenberg and Swendsen (1988) to estimate a
correlation function with a slightly different sea-quark mass than the one
simulated with is very common in LQCD; see for example Refs. Aoki  _et al._
(2010); Luscher and Schaefer (2013); Blum  _et al._ (2016). A nice discussion
of mass reweighting, including single flavor reweighting is found in Refs.
Finkenrath  _et al._ (2013); Leder  _et al._ (2014).

In our case, we are interested in reweighting the computation from the hybrid
Monte Carlo (HMC) simulated charm quark mass, $m_{c}^{\rm HMC}$, to the
physical quark mass, $m_{c}^{\rm phys}$ which requires an estimate of the
ratio of the fermion determinant with the physical mass to the determinant
with the HMC mass. If the mass shift is $a\;{\delta}m_{c}=m_{c}^{\rm
phys}-m_{c}^{\rm HMC}$ then up to $\mathcal{O}(\delta m^{2})$ this ratio can
be written (including the quarter-root arising from rooted-staggered fermions)

$w^{4}_{i}=\det\left[\mathbbm{1}+\delta m\;(D[U_{i}]+m_{c}^{\rm
phys})^{-1}\right]$ (41)

for each configuration $U_{i}$ and observables may be computed using the
weight $w$,

$\langle O\rangle_{\rm phys}=\frac{\langle Ow\rangle_{\rm HMC}}{\langle
w\rangle_{\rm HMC}}.$ (42)

We can use two methods to stochastically estimate $w$ for each configuration.
First, by rewriting the determinant as the exponential of a trace-log, one
finds

$w_{i}^{\rm tr}=\exp\frac{\delta m}{4}\text{tr}\left[(D[U_{i}]+m_{c}^{\rm
phys})^{-1}\right]+\mathcal{O}(\delta m^{2})$ (43)

and we can use vectors of complex Gaussian noise $\eta$ to estimate the trace,

$\displaystyle\text{tr}\left[(D[U_{i}]+m_{c}^{\rm phys})^{-1}\right]$
$\displaystyle\approx\frac{1}{N_{\eta}}\sum_{i}\eta_{i}^{\dagger}(D[U_{i}]+m_{c}^{\rm
phys})^{-1}\eta_{i}\,,$ (44) $\displaystyle\eta$
$\displaystyle\sim\frac{1}{\pi^{V}}\exp-\eta^{\dagger}\eta\,,$ (45)

where $V$ is the size of each $\eta$ vector.

Alternatively, we may estimate the determinant in the reweighting factor (41)
using the identity,

$\det
A=\pi^{-V}\int\mathcal{D}\eta\;\mathcal{D}\eta^{\dagger}\;e^{-\eta^{\dagger}A^{-1}\eta},$
(46)

which is often used to implement pseudofermions. Up to $\mathcal{O}(\delta
m^{2})$, this becomes

$w^{4}_{i}=\pi^{-V}\int\mathcal{D}\eta\;\mathcal{D}\eta^{\dagger}\;e^{-\eta^{\dagger}\eta}e^{\delta
m\;\eta^{\dagger}(D[U_{i}]+m_{c}^{\rm phys})^{-1}\eta}\,,$ (47)

which tells us to draw $\eta$ according to the same Gaussian (45) and estimate

$w^{\rm ps}_{i}=\left(\frac{1}{N_{\eta}}\sum_{j}e^{\delta
m\;\eta_{j}^{\dagger}(D[U_{i}]+m_{c}^{\rm phys})^{-1}\eta_{j}}\right)^{1/4}.$
(48)

Both the trace method (43) and the pseudofermion method (48) are only valid to
$\mathcal{O}({\delta}m^{2})$; when they agree we assume those corrections are
under control. In order to stabilize the numerical estimate of the reweighting
factors, it is also common to split the reweighting factor into a product of
reweighting factors where each is computed with a fraction of the full mass
shift Hasenbusch (2001); Hasenfratz  _et al._ (2008); Liu  _et al._ (2013).
For example, with a simulated mass of $m_{1}$ and target mass of
$m_{1}+{\Delta}m$, one could use two steps of ${\Delta}m/2$ and estimate the
reweighting factor with the trace method,

$\displaystyle w_{i}^{\rm tr}$ $\displaystyle=w^{\rm tr}_{i,1}w^{\rm
tr}_{i,2}\,,$ (49) $\displaystyle w^{\rm tr}_{i,1}$
$\displaystyle=\exp\frac{\Delta
m/2}{4}\frac{1}{N_{\eta}}\sum_{j}\eta^{\dagger}_{j}\left(D[U_{i}]+m_{1}+\frac{\Delta
m}{2}\right)^{-1}\eta_{j}\,,$ $\displaystyle w^{\rm tr}_{i,2}$
$\displaystyle=\exp\frac{\Delta
m/2}{4}\frac{1}{N_{\theta}}\sum_{j}\theta^{\dagger}_{j}\left(D[U_{i}]+m_{1}+\Delta
m\right)^{-1}\theta_{j}\,,$

using independently sampled complex Gaussian noise $\eta$ and $\theta$. Of
course, one may split the shift ${\Delta}m$ into finer steps if needed, for
increased computational cost.

The reweighting factor accounts for a change in the action and is exponential
in the spacetime volume. This can lead to numerical under- or overflow. As a
cure, we factor out the average reweighting factor. Recognizing the trace of
the inverse Dirac operator on a configuration $U_{i}$ as the scalar quark
density times the lattice volume,

$\left(V\bar{q}q\right)_{i}=\text{tr}\left[(D[U_{i}]+m_{q})^{-1}\right]\,,$
(50)

we can rescale $w$, shifting by the ensemble average $V\langle\bar{c}c\rangle$
computed via (44). For example, rescaling the trace method (43) gives

$w_{i}^{\rm
tr}=\exp\frac{{\delta}m}{4}\left(\text{tr}\left[(D[U_{i}]+m_{c}^{\rm
phys})^{-1}\right]-V\langle\bar{c}c\rangle\right);$ (51)

such a rescaling cancels exactly in the reweighting procedure (42). A similar
rescaling cures the pseudofermion method (48). If we split the mass shift as
in (49), each factor of the weight may be independently so stabilized.

On the a06m310L ensemble the lattice volume and shift in the mass are
$V=72^{3}\times 96$ and

$\displaystyle{\delta}a_{06}m_{c}$ $\displaystyle=a_{06}m_{c}^{\rm
phys}-a_{06}m_{c}^{\rm HMC}$ $\displaystyle=-0.0281(4)\,.$ (52)

While the shift in the mass is only about 10% of the physical charm quark
mass, it is of the order of the physical strange quark mass. In order to
stabilize the numerical estimate of the reweighting factors, we split this
mass shift into ten equal steps, and for each step, we used $N_{\eta}=128$
independent Gaussian random noise vectors. For each step in the reweighting,
we used the same Naik value of ${\epsilon}_{N}=-0.0533$ as was used in the
original HMC. This ensures that the Dirac operator only differs from one mass
to the next by the quark mass itself. As the Naik term is used to improve the
approach to the continuum limit, this is a valid choice to make as it results
in a slightly different approach to the continuum than if one simulated at the
physical charm quark mass with the optimized Naik value for that mass.

In Fig. 9, we show the reweighting factors for each of the mass steps, with
the bottom panel having a mass $a_{06}m_{c}=0.28319$ closest to the HMC mass
and the second from top panel having $a_{06}m_{c}^{\rm phys}=0.2579$, scaled
for numerical stability. The top panel is the resulting product reweighting
factor normalized by the average reweighting factor. One observes that there
are a few large reweighting factors of $\mathrm{O}(100)$. We have verified
that the trace estimation method (43) and the pseudofermion method (48)
produce comparable normalized reweighting factors. The large reweighting
factors are likely due to the parent HMC distribution of configurations having
a suboptimal overlap with the physical, target distribution.

Figure 9:  Distribution of reweighting factors, $r$, for the a06m310L
ensemble. The reweighting was performed with ten equal steps of the mass
difference from $a_{06}m_{c}^{\rm HMC}=0.286$ to $a_{06}m_{c}^{\rm
phys}=0.2579$. The reweighting factors from each step are shown in the bottom
ten subpanels with the product reweighting factor shown in the top panel. The
a06m310L ensemble was generated with two different streams of equal length.
For more information on the ensembles, see Ref. Miller  _et al._ (2020).

###  A.1 Reweighted spectrum

The next task is to understand how this reweighting impacts the extracted
spectrum. To aid in this discussion, we reiterate our strategy for fitting the
correlation functions. Because the noise of the omega baryon correlation
function grows in Euclidean time, it is the most challenging to fit and so we
focus our discussion on the omega. Our strategy is to find a good quality fit
to the correlation function for which the extracted ground state energy is
stable against the number of states and the time-range used in the analysis.
For a given $t_{\rm min}$, we opt to chose the simplest model which satisfies
this criteria, which amounts to picking the minimum number of excited states
possible.

The SS correlation functions are positive definite, therefore implying that
the excited state contamination of the effective mass must come from above.
When examining the SS omega-baryon effective mass on the a06m310L ensemble,
one observes that around $t=25$, the effective mass stops decreasing, and even
increases a little. Because this is not allowed for a positive definite
correlation function, we can conclude this must be due to a correlated
stochastic fluctuation; see Fig. 10. In the rewighted effective mass, one
observes more dramatic behavior of the effective mass beginning around the
same time. To be conservative, we set $t_{\rm max}=30$ in our analysis as this
allows the analysis to be sensitive to these stochastic fluctuations, which
fluctuate in the opposite direction between the unweighted and reweighted
configurations.

As a comparison, we also show the reweighting factors and reweighted omega
baryon effective mass on the a12m130 ensemble (Fig. 11) where the charm quark
was 2% different from its physical value. In this case, the reweighting
factors are much easier to estimate, and we do not observe any large values.

The lower panels in Fig. 10 show the extracted ground state mass as a function
of $t_{\rm min}$ and $n_{\rm state}$. We observe that an $n_{\rm state}=2$ fit
beginning at $t_{\rm min}=19$ and $15$ for the unweighted and reweighted
correlation functions satisfies our optimization criteria. These lead to the
estimate value of $m_{\Omega}$ given in Eq. (14). In Table 6, we also show the
reweighted value of $m_{\pi}$ and $m_{K}$ on this a06m310L ensemble. While the
pion and kaon have a statistically significant shift from the reweighting,
when we use the reweighted values of $m_{\pi}$, $m_{K}$ and $m_{\Omega}$ from
this ensemble, the final extrapolated value of $w_{0}m_{\Omega}$ is within 1
standard deviation of the completely unweighted analysis. As the a06m310L
ensemble has the largest potential change from reweighting, we conclude that
at the level of precision we currently have, our results are not sensitive to
the slight mistuning of the charm quark mass from its physical value on each
of the configurations.

Figure 10:  Same as Fig. 1 for the a06m310L ensemble.

Figure 11:  The reweighting factors for the a12m130 ensemble (left) and the reweighted omega baryon effective mass (right).  Table 6:  The pion, kaon and omega masses on the a06m310L ensemble with and without reweighting as well as the correlated difference.  State | Unweighted | Reweighted | Difference
---|---|---|---
$am_{\pi}$ | $0.09456(06)$ | $0.09419(14)$ | $-0.00037(14)$
$am_{K}$ | $0.16204(07)$ | $0.16173(16)$ | $-0.00032(15)$
$am_{\Omega}$ | $0.5069(21)$ | $0.5065(29)$ | $-0.0004(34)$

##  Appendix B MODELS INCLUDED IN AVERAGE

We use $w_{0,\rm orig}m_{\Omega}$ as an example to demonstrate the model
averaging and naming conventions for the various models. The other gradient-
flow scale studies have identical in form, extrapolation functions. Table 7
provides the various models and their relative weights for $w_{0,\rm
orig}m_{\Omega}$. Of note, in all four gradient-flow scale extrapolations, the
extrapolations which make use of $l_{\Omega}$ and $s_{\Omega}$ versus $l_{F}$
and $s_{F}$, all have relatively smaller weights in the averaging, such that
in all cases, these extrapolations could be dropped from the model averaging
without effecting the final result. A reprodution of this model averaging as
well as for the other three gradient-flow scales can be obtained by
downloading and running the associated analysis software at
https://github.com/callat-qcd/project_scale_setting_mdwf_hisq.

Table 7:  Models included in final model average of $w_{0,\rm orig}m_{\Omega}$. Moddel | chi2/dof | $Q$ | logGBF | weight | $w_{0}m_{\Omega}$ |  $w_{0}$[fm]
---|---|---|---|---|---|---
xpt_n3lo_alphaS_F | 1.065 | 0.378 | 64.919 | 0.077 | 1.448(10) | 0.1709(12)
xpt_n3lo_F | 1.070 | 0.372 | 64.916 | 0.077 | 1.4467(91) | 0.1707(11)
taylor_n3lo_alphaS_F | 1.089 | 0.350 | 64.863 | 0.073 | 1.449(10) | 0.1709(12)
xpt_n3lo_FV_alphaS_F | 1.068 | 0.374 | 64.861 | 0.073 | 1.448(10) | 0.1709(12)
xpt_n3lo_FV_F | 1.074 | 0.367 | 64.858 | 0.073 | 1.4467(91) | 0.1707(11)
taylor_n3lo_FV_alphaS_F | 1.088 | 0.350 | 64.857 | 0.073 | 1.449(10) | 0.1709(12)
taylor_n3lo_F | 1.094 | 0.343 | 64.857 | 0.073 | 1.4470(90) | 0.1707(11)
taylor_n3lo_FV_F | 1.094 | 0.343 | 64.851 | 0.072 | 1.4470(90) | 0.1707(11)
xpt_n2lo_alphaS_F | 1.124 | 0.310 | 64.610 | 0.057 | 1.449(10) | 0.1710(12)
xpt_n2lo_F | 1.130 | 0.304 | 64.600 | 0.056 | 1.4475(91) | 0.1708(11)
xpt_n2lo_FV_alphaS_F | 1.128 | 0.306 | 64.551 | 0.054 | 1.449(10) | 0.1710(12)
xpt_n2lo_FV_F | 1.134 | 0.299 | 64.540 | 0.053 | 1.4475(91) | 0.1708(11)
taylor_n2lo_alphaS_F | 1.161 | 0.272 | 64.410 | 0.047 | 1.450(10) | 0.1711(12)
taylor_n2lo_FV_alphaS_F | 1.161 | 0.272 | 64.403 | 0.046 | 1.450(10) | 0.1711(12)
taylor_n2lo_F | 1.168 | 0.265 | 64.392 | 0.046 | 1.4480(90) | 0.1708(11)
taylor_n2lo_FV_F | 1.167 | 0.266 | 64.386 | 0.045 | 1.4480(90) | 0.1708(11)
xpt_n2lo_FV_alphaS_O | 0.971 | 0.498 | 59.864 | 0.000 | 1.441(14) | 0.1700(16)
xpt_n3lo_FV_alphaS_O | 0.937 | 0.545 | 59.813 | 0.000 | 1.440(13) | 0.1699(16)
taylor_n3lo_alphaS_O | 1.089 | 0.349 | 59.792 | 0.000 | 1.444(13) | 0.1704(15)
xpt_n2lo_alphaS_O | 0.985 | 0.480 | 59.789 | 0.000 | 1.441(13) | 0.1701(16)
xpt_n3lo_FV_O | 0.937 | 0.545 | 59.515 | 0.000 | 1.440(12) | 0.1699(15)
taylor_n3lo_FV_alphaS_O | 1.089 | 0.349 | 59.394 | 0.000 | 1.444(13) | 0.1704(15)
xpt_n3lo_alphaS_O | 0.949 | 0.529 | 59.380 | 0.000 | 1.441(13) | 0.1700(16)
xpt_n2lo_FV_O | 0.971 | 0.498 | 59.234 | 0.000 | 1.441(12) | 0.1700(14)
xpt_n2lo_O | 0.985 | 0.480 | 59.172 | 0.000 | 1.441(13) | 0.1700(15)
taylor_n3lo_O | 1.090 | 0.348 | 58.874 | 0.000 | 1.444(12) | 0.1703(14)
taylor_n3lo_FV_O | 1.090 | 0.349 | 58.798 | 0.000 | 1.444(12) | 0.1703(14)
taylor_n2lo_alphaS_O | 1.131 | 0.303 | 58.742 | 0.000 | 1.445(13) | 0.1705(15)
taylor_n2lo_FV_O | 1.131 | 0.302 | 58.607 | 0.000 | 1.444(12) | 0.1704(14)
xpt_n3lo_O | 0.949 | 0.528 | 58.329 | 0.000 | 1.440(12) | 0.1699(14)
taylor_n2lo_FV_alphaS_O | 1.131 | 0.303 | 58.292 | 0.000 | 1.445(13) | 0.1705(15)
taylor_n2lo_O | 1.132 | 0.302 | 58.070 | 0.000 | 1.444(12) | 0.1704(14)
Model average |  |  |  |  | 1.4481(98)(11) | 0.1709(12)(01)

A few example extrapolation formula are given below to demonstrate the naming
convention:

$\displaystyle w_{0}m_{\Omega}$
$\displaystyle=\underbrace{c_{0}}_{\text{LO}}+\underbrace{{\delta}^{\text{NLO}}(l_{F},s_{F})}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral
NLO}}+\underbrace{{\delta}^{\text{NLO}}_{a,F}}_{\text{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}
disc NLO}}$ [xpt_nlo_fv_F] (53a) $\displaystyle w_{0}m_{\Omega}$
$\displaystyle=\underbrace{c_{0}}_{\text{LO}}+\underbrace{{\delta}^{\text{NLO}}(l_{F},s_{F})}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral
NLO}}+\underbrace{{\delta}^{\text{{N${}^{2}$LO}}}(l_{F},s_{F})}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral
{N${}^{2}$LO}}}+\underbrace{{\delta}^{\text{{N${}^{2}$LO}}}_{\ln}}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral log {N${}^{2}$LO}}}$ [xpt_n2lo_fv_F]
$\displaystyle\phantom{=}+\underbrace{{\delta}^{\text{NLO}}_{a,F}}_{\text{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}
disc
NLO}}+\underbrace{{\delta}^{\text{{N${}^{2}$LO}}}_{a,F}}_{\text{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}
disc
{N${}^{2}$LO}}}+\underbrace{{\delta}^{\text{{N${}^{2}$LO}}}_{L,F}}_{\text{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}
FV {N${}^{2}$LO}}}$ (53b) $\displaystyle w_{0}m_{\Omega}$
$\displaystyle=\underbrace{c_{0}}_{\text{LO}}+\underbrace{{\delta}^{\text{NLO}}(l_{F},s_{F})}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral
NLO}}+\underbrace{{\delta}^{\text{{N${}^{2}$LO}}}(l_{F},s_{F})}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral
{N${}^{2}$LO}}}+\underbrace{{\delta}^{\text{{N${}^{3}$LO}}}(l_{F},s_{F})}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral
{N${}^{3}$LO}}}+\underbrace{{\delta}^{\text{{N${}^{2}$LO}}}_{\ln}}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral log
{N${}^{2}$LO}}}+\underbrace{{\delta}^{\text{{N${}^{3}$LO}}}_{\ln}}_{\text{\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}
chiral logs {N${}^{3}$LO}}}$ [xpt_n3lo_fv_F]
$\displaystyle\phantom{=}+\underbrace{{\delta}^{\text{NLO}}_{a,F}}_{\text{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}
disc
NLO}}+\underbrace{{\delta}^{\text{{N${}^{2}$LO}}}_{a,F}}_{\text{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}
disc
{N${}^{2}$LO}}}+\underbrace{{\delta}^{\text{{N${}^{3}$LO}}}_{a,F}}_{\text{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}
disc
{N${}^{3}$LO}}}+\underbrace{{\delta}^{\text{{N${}^{2}$LO}}}_{L,F}}_{\text{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}
FV
{N${}^{2}$LO}}}+\underbrace{{\delta}^{\text{{N${}^{3}$LO}}}_{L,F}}_{\text{\color[rgb]{1,.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,.5,0}
FV {N${}^{3}$LO}}}$ (53c)

where the  chiral,  discretization and  finite volume corrections are defined
in Eqs. (III.1), (III.3) and (III.2) respectively.

## References

  * Aoki  _et al._ (2020) S. Aoki _et al._ (Flavour Lattice Averaging Group), “FLAG Review 2019: Flavour Lattice Averaging Group (FLAG),” Eur. Phys. J. C 80, 113 (2020), arXiv:1902.08191 [hep-lat] .
  * Aoyama  _et al._ (2020) T. Aoyama _et al._ , “The anomalous magnetic moment of the muon in the Standard Model,” Phys. Rept. 887, 1–166 (2020), arXiv:2006.04822 [hep-ph] .
  * Sommer (1994) R. Sommer, “A New way to set the energy scale in lattice gauge theories and its applications to the static force and alpha-s in SU(2) Yang-Mills theory,” Nucl. Phys. B 411, 839–854 (1994), arXiv:hep-lat/9310022 .
  * Lüscher (2010) Martin Lüscher, “Properties and uses of the Wilson flow in lattice QCD,” JHEP 1008, 071 (2010), arXiv:1006.4518 [hep-lat] .
  * Borsanyi  _et al._ (2012) Szabolcs Borsanyi _et al._ , “High-precision scale setting in lattice QCD,” JHEP 09, 010 (2012), arXiv:1203.4469 [hep-lat] .
  * Deuzeman and Wenger (2012) Albert Deuzeman and Urs Wenger, “Gradient flow and scale setting for twisted mass fermions,” PoS LATTICE2012, 162 (2012).
  * Bruno and Sommer (2014) Mattia Bruno and Rainer Sommer (ALPHA), “On the $N_{f}$-dependence of gluonic observables,” PoS LATTICE2013, 321 (2014), arXiv:1311.5585 [hep-lat] .
  * Dowdall  _et al._ (2013) R. J. Dowdall, C. T. H. Davies, G. P. Lepage, and C. McNeile, “Vus from pi and K decay constants in full lattice QCD with physical u, d, s and c quarks,” Phys. Rev. D88, 074504 (2013), arXiv:1303.1670 [hep-lat] .
  * Bornyakov  _et al._ (2015) V.G. Bornyakov _et al._ , “Wilson flow and scale setting from lattice QCD,” (2015), arXiv:1508.05916 [hep-lat] .
  * Bazavov  _et al._ (2014a) A. Bazavov _et al._ (HotQCD), “Equation of state in ( 2+1 )-flavor QCD,” Phys. Rev. D 90, 094503 (2014a), arXiv:1407.6387 [hep-lat] .
  * Blum  _et al._ (2016) T. Blum _et al._ (RBC, UKQCD), “Domain wall QCD with physical quark masses,” Phys. Rev. D 93, 074505 (2016), arXiv:1411.7017 [hep-lat] .
  * Bazavov  _et al._ (2016) A. Bazavov _et al._ (MILC), “Gradient flow and scale setting on MILC HISQ ensembles,” Phys. Rev. D93, 094510 (2016), arXiv:1503.02769 [hep-lat] .
  * Bruno  _et al._ (2017) Mattia Bruno, Tomasz Korzec, and Stefan Schaefer, “Setting the scale for the CLS $2+1$ flavor ensembles,” Phys. Rev. D 95, 074504 (2017), arXiv:1608.08900 [hep-lat] .
  * Borsanyi  _et al._ (2021) Sz. Borsanyi _et al._ , “Leading hadronic contribution to the muon 2 magnetic moment from lattice QCD,” Nature (2021), 10.1038/s41586-021-03418-1, arXiv:2002.12347 [hep-lat] .
  * Berkowitz  _et al._ (2017a) E. Berkowitz, C. Bouchard, C.C Chang, M.A. Clark, B Joo, T. Kurth, C. Monahan, A. Nicholson, K. Orginos, E. Rinaldi, P. Vranas, and A. Walker-Loud, “Möbius Domain-Wall fermions on gradient-flowed dynamical HISQ ensembles,” Phys. Rev. D96, 054513 (2017a), arXiv:1701.07559 [hep-lat] .
  * Bazavov  _et al._ (2013) A. Bazavov _et al._ (MILC), “Lattice QCD ensembles with four flavors of highly improved staggered quarks,” Phys. Rev. D87, 054505 (2013), arXiv:1212.4768 [hep-lat] .
  * Bhattacharya  _et al._ (2015) Tanmoy Bhattacharya, Vincenzo Cirigliano, Saul Cohen, Rajan Gupta, Anosh Joseph, Huey-Wen Lin, and Boram Yoon (PNDME), “Iso-vector and Iso-scalar Tensor Charges of the Nucleon from Lattice QCD,” Phys. Rev. D92, 094511 (2015), arXiv:1506.06411 [hep-lat] .
  * Clark  _et al._ (2010) M.A. Clark, R. Babich, K. Barros, R.C. Brower, and C. Rebbi, “Solving Lattice QCD systems of equations using mixed precision solvers on GPUs,” Comput.Phys.Commun. 181, 1517–1528 (2010), https://github.com/lattice/quda, arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2011) R. Babich, M.A. Clark, B. Joo, G. Shi, R.C. Brower, _et al._ , “Scaling Lattice QCD beyond 100 GPUs,” (2011), arXiv:1109.2935 [hep-lat] .
  * Renner  _et al._ (2005) Dru Bryant Renner, W. Schroers, R. Edwards, George Tamminga Fleming, Ph. Hagler, John W. Negele, K. Orginos, A. V. Pochinski, and D. Richards (LHP), “Hadronic physics with domain-wall valence and improved staggered sea quarks,”  _Lattice field theory. Proceedings, 22nd International Symposium, Lattice 2004, Batavia, USA, June 21-26, 2004_ , Nucl. Phys. Proc. Suppl. 140, 255–260 (2005), [,255(2004)], arXiv:hep-lat/0409130 [hep-lat] .
  * Bar  _et al._ (2003) Oliver Bar, Gautam Rupak, and Noam Shoresh, “Simulations with different lattice Dirac operators for valence and sea quarks,” Phys. Rev. D67, 114505 (2003), arXiv:hep-lat/0210050 [hep-lat] .
  * Brower  _et al._ (2017) Richard C. Brower, Harmut Neff, and Kostas Orginos, “The möbius domain wall fermion algorithm,” Comput. Phys. Commun. 220, 1–19 (2017), arXiv:1206.5214 [hep-lat] .
  * Kaplan (1992) David B. Kaplan, “A Method for simulating chiral fermions on the lattice,” Phys. Lett. B288, 342–347 (1992), arXiv:hep-lat/9206013 [hep-lat] .
  * Shamir (1993) Yigal Shamir, “Chiral fermions from lattice boundaries,” Nucl. Phys. B406, 90–106 (1993), arXiv:hep-lat/9303005 [hep-lat] .
  * Furman and Shamir (1995) Vadim Furman and Yigal Shamir, “Axial symmetries in lattice QCD with Kaplan fermions,” Nucl. Phys. B439, 54–78 (1995), arXiv:hep-lat/9405004 [hep-lat] .
  * Follana  _et al._ (2007) E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD), “Highly improved staggered quarks on the lattice, with applications to charm physics,” Phys. Rev. D75, 054502 (2007), arXiv:hep-lat/0610092 [hep-lat] .
  * Narayanan and Neuberger (2006) R. Narayanan and H. Neuberger, “Infinite N phase transitions in continuum Wilson loop operators,” JHEP 0603, 064 (2006), arXiv:hep-th/0601210 [hep-th] .
  * Lüscher and Weisz (2011) M. Lüscher and P. Weisz, “Perturbative analysis of the gradient flow in non-abelian gauge theories,” JHEP 1102, 051 (2011), arXiv:1101.0963 [hep-th] .
  * Lüscher (2013) Martin Lüscher, “Chiral symmetry and the Yang–Mills gradient flow,” JHEP 1304, 123 (2013), arXiv:1302.5246 [hep-lat] .
  * Morningstar and Peardon (2004) Colin Morningstar and Mike J. Peardon, “Analytic smearing of SU(3) link variables in lattice QCD,” Phys.Rev. D69, 054501 (2004), arXiv:hep-lat/0311018 [hep-lat] .
  * Bouchard  _et al._ (2017) Chris Bouchard, Chia Cheng Chang, Thorsten Kurth, Kostas Orginos, and Andre Walker-Loud, “On the Feynman-Hellmann Theorem in Quantum Field Theory and the Calculation of Matrix Elements,” Phys. Rev. D96, 014504 (2017), arXiv:1612.06963 [hep-lat] .
  * Berkowitz  _et al._ (2017b) E. Berkowitz, C. Bouchard, D. B. Brantley, C.C Chang, M.A. Clark, N. Garron, B Joó, T. Kurth, C. Monahan, H. Monge-Camacho, A. Nicholson, K. Orginos, E. Rinaldi, P. Vranas, and A. Walker-Loud, “An accurate calculation of the nucleon axial charge with lattice QCD,” (2017b), arXiv:1704.01114 [hep-lat] .
  * Chang  _et al._ (2018) C. C. Chang, A. Nicholson, E. Rinaldi, E. Berkowitz, N. Garron, D. A. Brantley, H. Monge-Camacho, C. J. Monahan, C. Bouchard, M. A. Clark, B. Joó, T. Kurth, K. Orginos, P. Vranas, and A. Walker-Loud, “A per-cent-level determination of the nucleon axial coupling from quantum chromodynamics,” Nature 558, 91–94 (2018), arXiv:1805.12130 [hep-lat] .
  * Berkowitz  _et al._ (2018a) Evan Berkowitz, M. A. Clark, A. Gambhir, K. McElvain, A. Nicholson, E. Rinalid, P. Vranas, A. Walker-Loud, C. C. Chang, B. Joó, T. Kurth, and K. Orginos, “Simulating the weak death of the neutron in a femtoscale universe with near-Exascale computing,”  (2018a), 10.1109/SC.2018.00058, arXiv:1810.01609 [hep-lat] .
  * Nicholson  _et al._ (2018) A. Nicholson, E. Berkowitz, H. Monge-Camacho, D. Brantley, N. Garron, C. C. Chang, E. Rinaldi, M. A. Clark, B. Joó, T. Kurth, B. Tiburzi, P. Vranas, and A. Walker-Loud, “Heavy physics contributions to neutrinoless double beta decay from QCD,” Phys. Rev. Lett. 121, 172501 (2018), arXiv:1805.02634 [nucl-th] .
  * Miller  _et al._ (2020) Nolan Miller _et al._ , “$F_{K}/F_{\pi}$ from M\”{o}bius Domain-Wall fermions solved on gradient-flowed HISQ ensembles,” Phys. Rev. D 102, 034507 (2020), arXiv:2005.04795 [hep-lat] .
  * Follana  _et al._ (2008) E. Follana, C.T.H. Davies, G.P. Lepage, and J. Shigemitsu (HPQCD, UKQCD), “High Precision determination of the pi, K, D and D(s) decay constants from lattice QCD,” Phys. Rev. Lett. 100, 062002 (2008), arXiv:0706.1726 [hep-lat] .
  * Blossier  _et al._ (2009) B. Blossier _et al._ (ETM), “Pseudoscalar decay constants of kaon and D-mesons from $N_{f}=2$ twisted mass Lattice QCD,” JHEP 07, 043 (2009), arXiv:0904.0954 [hep-lat] .
  * Durr  _et al._ (2010) S. Durr, Z. Fodor, C. Hoelbling, S.D. Katz, S. Krieg, T. Kurth, L. Lellouch, T. Lippert, A. Ramos, and K.K. Szabo, “The ratio FK/Fpi in QCD,” Phys. Rev. D 81, 054507 (2010), arXiv:1001.4692 [hep-lat] .
  * Carrasco  _et al._ (2015) N. Carrasco _et al._ , “Leptonic decay constants $f_{K},f_{D},$ and $f_{{D}_{s}}$ with $N_{f}=2+1+1$ twisted-mass lattice QCD,” Phys. Rev. D91, 054507 (2015), arXiv:1411.7908 [hep-lat] .
  * Bazavov  _et al._ (2014b) A. Bazavov _et al._ (Fermilab Lattice, MILC), “Charmed and light pseudoscalar meson decay constants from four-flavor lattice QCD with physical light quarks,” Phys. Rev. D90, 074509 (2014b), arXiv:1407.3772 [hep-lat] .
  * Dürr  _et al._ (2017) Stephan Dürr _et al._ , “Leptonic decay-constant ratio $f_{K}/f_{\pi}$ from lattice QCD using 2+1 clover-improved fermion flavors with 2-HEX smearing,” Phys. Rev. D 95, 054513 (2017), arXiv:1601.05998 [hep-lat] .
  * Bornyakov  _et al._ (2017) V.G. Bornyakov, R. Horsley, Y. Nakamura, H. Perlt, D. Pleiter, P.E.L. Rakow, G. Schierholz, A. Schiller, H. Stüben, and J.M. Zanotti (QCDSF–UKQCD), “Flavour breaking effects in the pseudoscalar meson decay constants,” Phys. Lett. B 767, 366–373 (2017), arXiv:1612.04798 [hep-lat] .
  * Bazavov  _et al._ (2018) A. Bazavov _et al._ , “$B$\- and $D$-meson leptonic decay constants from four-flavor lattice QCD,” Phys. Rev. D98, 074512 (2018), arXiv:1712.09262 [hep-lat] .
  * Allton  _et al._ (2008) C. Allton _et al._ (RBC-UKQCD), “Physical Results from 2+1 Flavor Domain Wall QCD and SU(2) Chiral Perturbation Theory,” Phys. Rev. D78, 114509 (2008), arXiv:0804.0473 [hep-lat] .
  * Walker-Loud  _et al._ (2009) A. Walker-Loud _et al._ , “Light hadron spectroscopy using domain wall valence quarks on an Asqtad sea,” Phys. Rev. D79, 054502 (2009), arXiv:0806.4549 [hep-lat] .
  * Lin  _et al._ (2009) Huey-Wen Lin _et al._ (Hadron Spectrum), “First results from 2+1 dynamical quark flavors on an anisotropic lattice: Light-hadron spectroscopy and setting the strange-quark mass,” Phys. Rev. D 79, 034502 (2009), arXiv:0810.3588 [hep-lat] .
  * Durr  _et al._ (2008) S. Durr _et al._ , “Ab-Initio Determination of Light Hadron Masses,” Science 322, 1224–1227 (2008), arXiv:0906.3599 [hep-lat] .
  * Alexandrou  _et al._ (2009) C. Alexandrou, R. Baron, J. Carbonell, V. Drach, P. Guichon, K. Jansen, T. Korzec, and O. Pene (ETM), “Low-lying baryon spectrum with two dynamical twisted mass fermions,” Phys. Rev. D 80, 114503 (2009), arXiv:0910.2419 [hep-lat] .
  * Alexandrou  _et al._ (2014) C. Alexandrou, V. Drach, K. Jansen, C. Kallidonis, and G. Koutsou, “Baryon spectrum with $N_{f}=2+1+1$ twisted mass fermions,” Phys. Rev. D 90, 074501 (2014), arXiv:1406.4310 [hep-lat] .
  * Brantley  _et al._ (2016) David A. Brantley, Balint Joo, Ekaterina V. Mastropas, Emanuele Mereghetti, Henry Monge-Camacho, Brian C. Tiburzi, and Andre Walker-Loud, “Strong isospin violation and chiral logarithms in the baryon spectrum,” (2016), arXiv:1612.07733 [hep-lat] .
  * Edwards and Joo (2005) Robert G. Edwards and Balint Joo (SciDAC Collaboration, LHPC Collaboration, UKQCD Collaboration), “The Chroma software system for lattice QCD,” Nucl.Phys.Proc.Suppl. 140, 832 (2005), arXiv:hep-lat/0409003 [hep-lat] .
  * Basak  _et al._ (2005a) S. Basak, R. G. Edwards, G. T. Fleming, U. M. Heller, C. Morningstar, D. Richards, I. Sato, and S. Wallace, “Group-theoretical construction of extended baryon operators in lattice QCD,” Phys. Rev. D72, 094506 (2005a), arXiv:hep-lat/0506029 [hep-lat] .
  * Basak  _et al._ (2005b) Subhasish Basak, Robert Edwards, George T. Fleming, Urs M. Heller, Colin Morningstar, David Richards, Ikuro Sato, and Stephen J. Wallace (Lattice Hadron Physics (LHPC)), “Clebsch-Gordan construction of lattice interpolating fields for excited baryons,” Phys. Rev. D72, 074501 (2005b), arXiv:hep-lat/0508018 [hep-lat] .
  * Berkowitz  _et al._ (2017c) Evan Berkowitz, Amy Nicholson, Chia Cheng Chang, Enrico Rinaldi, M. A. Clark, Bálint Joó, Thorsten Kurth, Pavlos Vranas, and André Walker-Loud, “Calm Multi-Baryon Operators,” in  _35th International Symposium on Lattice Field Theory (Lattice 2017) Granada, Spain, June 18-24, 2017_ (2017) arXiv:1710.05642 [hep-lat] .
  * Beane  _et al._ (2009) Silas R. Beane, William Detmold, Thomas C. Luu, Kostas Orginos, Assumpta Parreno, Martin Savage, Aaron Torok, and Andre Walker-Loud, “High Statistics Analysis using Anisotropic Clover Lattices: (I) Single Hadron Correlation Functions,” Phys.Rev. D79, 114502 (2009), arXiv:0903.2990 [hep-lat] .
  * Lepage  _et al._ (2002) G.P. Lepage, B. Clark, C.T.H. Davies, K. Hornbostel, P.B. Mackenzie, C. Morningstar, and H. Trottier, “Constrained curve fitting,” Nucl. Phys. B Proc. Suppl. 106, 12–20 (2002), arXiv:hep-lat/0110175 .
  * Jay and Neil (2020) William I. Jay and Ethan T. Neil, “Bayesian model averaging for analysis of lattice field theory results,” (2020), arXiv:2008.01069 [stat.ME] .
  * (59) The MILC Collaboration, “QUDA: a library for performing calculations in lattice qcd on graphics processing units (gpus), leveraging nvidia’s cuda platform.” https://github.com/lattice/quda.
  * (60) “MILC Collaboration code for lattice qcd calculations,” https://github.com/milc-qcd/milc_qcd.
  * Wolff (2004) Ulli Wolff (ALPHA), “Monte Carlo errors with less errors,” Comput. Phys. Commun. 156, 143–153 (2004), [Erratum: Comput.Phys.Commun. 176, 383 (2007)], arXiv:hep-lat/0306017 .
  * De Palma  _et al._ (2019) Barbara De Palma, Marco Erba, Luca Mantovani, and Nicola Mosco, “A Python program for the implementation of the $\Gamma$ -method for Monte Carlo simulations,” Comput. Phys. Commun. 234, 294–301 (2019), arXiv:1703.02766 [hep-lat] .
  * Davies  _et al._ (2003) C. Davies _et al._ , “The Determination of alpha(s) from lattice QCD with 2+1 flavors of dynamical quarks,” Nucl. Phys. B Proc. Suppl. 119, 595–597 (2003), arXiv:hep-lat/0209122 .
  * Fodor  _et al._ (2014) Zoltan Fodor, Kieran Holland, Julius Kuti, Santanu Mondal, Daniel Nogradi, and Chik Him Wong, “The lattice gradient flow at tree-level and its improvement,” JHEP 09, 018 (2014), arXiv:1406.0827 [hep-lat] .
  * Gell-Mann  _et al._ (1968) Murray Gell-Mann, R. J. Oakes, and B. Renner, “Behavior of current divergences under SU(3) x SU(3),” Phys. Rev. 175, 2195–2199 (1968).
  * Tanabashi  _et al._ (2018) M. Tanabashi _et al._ (Particle Data Group), “Review of Particle Physics,” Phys. Rev. D 98, 030001 (2018).
  * Ferrenberg and Swendsen (1988) A.M. Ferrenberg and R.H. Swendsen, “New Monte Carlo Technique for Studying Phase Transitions,” Phys. Rev. Lett. 61, 2635–2638 (1988).
  * Luscher (1986) M. Luscher, “Volume Dependence of the Energy Spectrum in Massive Quantum Field Theories. 1. Stable Particle States,” Commun. Math. Phys. 104, 177 (1986).
  * Jenkins and Manohar (1991a) Elizabeth Ellen Jenkins and Aneesh V. Manohar, “Baryon chiral perturbation theory using a heavy fermion Lagrangian,” Phys. Lett. B255, 558–562 (1991a).
  * Jenkins and Manohar (1991b) Elizabeth Ellen Jenkins and Aneesh V. Manohar, “Chiral corrections to the baryon axial currents,” Phys. Lett. B259, 353–358 (1991b).
  * Lebed (1994) Richard F. Lebed, “Baryon decuplet mass relations in chiral perturbation theory,” Nucl. Phys. B 430, 295–318 (1994), arXiv:hep-ph/9311234 .
  * Lebed and Luty (1994) Richard F. Lebed and Markus A. Luty, “Baryon masses at second order in chiral perturbation theory,” Phys. Lett. B 329, 479–485 (1994), arXiv:hep-ph/9401232 .
  * Tiburzi and Walker-Loud (2005) Brian C. Tiburzi and Andre Walker-Loud, “Decuplet baryon masses in partially quenched chiral perturbation theory,” Nucl. Phys. A748, 513–536 (2005), arXiv:hep-lat/0407030 [hep-lat] .
  * Torok  _et al._ (2010) Aaron Torok _et al._ (NPLQCD), “Meson-Baryon Scattering Lengths from Mixed-Action Lattice QCD,” Phys. Rev. D81, 074506 (2010), arXiv:0907.1913 [hep-lat] .
  * Ishikawa  _et al._ (2009) K. I. Ishikawa _et al._ (PACS-CS), “SU(2) and SU(3) chiral perturbation theory analyses on baryon masses in 2+1 flavor lattice QCD,” Phys. Rev. D80, 054502 (2009), arXiv:0905.0962 [hep-lat] .
  * Jenkins  _et al._ (2010) Elizabeth Ellen Jenkins, Aneesh V. Manohar, John W. Negele, and Andre Walker-Loud, “A Lattice Test of 1/N(c) Baryon Mass Relations,” Phys. Rev. D81, 014502 (2010), arXiv:0907.0529 [hep-lat] .
  * Walker-Loud (2012) Andre Walker-Loud, “Evidence for non-analytic light quark mass dependence in the baryon spectrum,” Phys. Rev. D86, 074509 (2012), arXiv:1112.2658 [hep-lat] .
  * Beane  _et al._ (2007) Silas R. Beane _et al._ (NPLQCD), “Hyperon nucleon scattering from fully-dynamical lattice QCD,” Nucl. Phys. A794, 62–72 (2007), arXiv:hep-lat/0612026 .
  * Tiburzi and Walker-Loud (2008) Brian C. Tiburzi and Andre Walker-Loud, “Hyperons in Two Flavor Chiral Perturbation Theory,” Phys.Lett. B669, 246–253 (2008), arXiv:0808.0482 [nucl-th] .
  * Jiang and Tiburzi (2009) Fu-Jiun Jiang and Brian C. Tiburzi, “Hyperon Axial Charges in Two-Flavor Chiral Perturbation Theory,” Phys. Rev. D 80, 077501 (2009), arXiv:0905.0857 [nucl-th] .
  * Jiang  _et al._ (2011) Fu-Jiun Jiang, Brian C. Tiburzi, and Andre Walker-Loud, “Kaon Thresholds and Two-Flavor Chiral Expansions for Hyperons,” Phys. Lett. B 695, 329–336 (2011), arXiv:0911.4721 [nucl-th] .
  * Jiang and Tiburzi (2010) Fu-Jiun Jiang and Brian C. Tiburzi, “Hyperon Electromagnetic Properties in Two-Flavor Chiral Perturbation Theory,” Phys. Rev. D 81, 034017 (2010), arXiv:0912.2077 [nucl-th] .
  * Weinberg (1979) Steven Weinberg, “Phenomenological Lagrangians,” Physica A96, 327–340 (1979).
  * Bar and Golterman (2014) Oliver Bar and Maarten Golterman, “Chiral perturbation theory for gradient flow observables,” Phys. Rev. D 89, 034505 (2014), [Erratum: Phys.Rev.D 89, 099905 (2014)], arXiv:1312.4999 [hep-lat] .
  * Gasser and Leutwyler (1987) J. Gasser and H. Leutwyler, “Light Quarks at Low Temperatures,” Phys. Lett. B184, 83–88 (1987).
  * Colangelo and Haefeli (2004) Gilberto Colangelo and Christoph Haefeli, “An Asymptotic formula for the pion decay constant in a large volume,” Phys. Lett. B 590, 258–264 (2004), arXiv:hep-lat/0403025 .
  * Bijnens and Rössler (2015) Johan Bijnens and Thomas Rössler, “Finite Volume at Two-loops in Chiral Perturbation Theory,” JHEP 01, 034 (2015), arXiv:1411.6384 [hep-lat] .
  * Sharpe and Singleton (1998) Stephen R. Sharpe and Robert L. Singleton, Jr, “Spontaneous flavor and parity breaking with Wilson fermions,” Phys. Rev. D58, 074501 (1998), arXiv:hep-lat/9804028 [hep-lat] .
  * Symanzik (1983a) K. Symanzik, “Continuum Limit and Improved Action in Lattice Theories. 1. Principles and phi**4 Theory,” Nucl. Phys. B226, 187–204 (1983a).
  * Symanzik (1983b) K. Symanzik, “Continuum Limit and Improved Action in Lattice Theories. 2. O(N) Nonlinear Sigma Model in Perturbation Theory,” Nucl. Phys. B226, 205–227 (1983b).
  * Balog  _et al._ (2009) Janos Balog, Ferenc Niedermayer, and Peter Weisz, “Logarithmic corrections to O(a**2) lattice artifacts,” Phys. Lett. B 676, 188–192 (2009), arXiv:0901.4033 [hep-lat] .
  * Balog  _et al._ (2010) Janos Balog, Ferenc Niedermayer, and Peter Weisz, “The Puzzle of apparent linear lattice artifacts in the 2d non-linear sigma-model and Symanzik’s solution,” Nucl. Phys. B 824, 563–615 (2010), arXiv:0905.1730 [hep-lat] .
  * Husung  _et al._ (2020) Nikolai Husung, Peter Marquard, and Rainer Sommer, “Asymptotic behavior of cutoff effects in Yang–Mills theory and in Wilson’s lattice QCD,” Eur. Phys. J. C 80, 200 (2020), arXiv:1912.08498 [hep-lat] .
  * Bar  _et al._ (2004) Oliver Bar, Gautam Rupak, and Noam Shoresh, “Chiral perturbation theory at $\mathcal{O}(a^{2})$ for lattice QCD,” Phys. Rev. D70, 034508 (2004), arXiv:hep-lat/0306021 [hep-lat] .
  * Bar  _et al._ (2005) Oliver Bar, Claude Bernard, Gautam Rupak, and Noam Shoresh, “Chiral perturbation theory for staggered sea quarks and Ginsparg-Wilson valence quarks,” Phys. Rev. D72, 054502 (2005), arXiv:hep-lat/0503009 [hep-lat] .
  * Golterman  _et al._ (2005) Maarten Golterman, Taku Izubuchi, and Yigal Shamir, “The Role of the double pole in lattice QCD with mixed actions,” Phys. Rev. D 71, 114508 (2005), arXiv:hep-lat/0504013 .
  * Tiburzi (2005) Brian C. Tiburzi, “Baryons with Ginsparg-Wilson quarks in a staggered sea,” Phys. Rev. D72, 094501 (2005), [Erratum: Phys. Rev.D79,039904(2009)], arXiv:hep-lat/0508019 [hep-lat] .
  * Chen  _et al._ (2006) Jiunn-Wei Chen, Donal O’Connell, Ruth S. Van de Water, and Andre Walker-Loud, “Ginsparg-Wilson pions scattering on a staggered sea,” Phys. Rev. D73, 074510 (2006), arXiv:hep-lat/0510024 [hep-lat] .
  * Chen  _et al._ (2007) Jiunn-Wei Chen, Donal O’Connell, and Andre Walker-Loud, “Two Meson Systems with Ginsparg-Wilson Valence Quarks,” Phys. Rev. D75, 054501 (2007), arXiv:hep-lat/0611003 [hep-lat] .
  * Orginos and Walker-Loud (2008) Kostas Orginos and Andre Walker-Loud, “Mixed meson masses with domain-wall valence and staggered sea fermions,” Phys. Rev. D77, 094505 (2008), arXiv:0705.0572 [hep-lat] .
  * Chen  _et al._ (2009a) Jiunn-Wei Chen, Donal O’Connell, and Andre Walker-Loud, “Universality of mixed action extrapolation formulae,” JHEP 04, 090 (2009a), arXiv:0706.0035 [hep-lat] .
  * Chen  _et al._ (2009b) Jiunn-Wei Chen, Maarten Golterman, Donal O’Connell, and Andre Walker-Loud, “Mixed Action Effective Field Theory: An Addendum,” Phys. Rev. D79, 117502 (2009b), arXiv:0905.2566 [hep-lat] .
  * hig (2005) (2005), Example given by M. Savage in a discussion on the inclusion of discretization effects in an effective field theory.
  * Bernard and Golterman (1992) Claude W. Bernard and Maarten F.L. Golterman, “Chiral perturbation theory for the quenched approximation of QCD,” Phys. Rev. D 46, 853–857 (1992), arXiv:hep-lat/9204007 .
  * Bernard and Golterman (1994) Claude W. Bernard and Maarten F. L. Golterman, “Partially quenched gauge theories and an application to staggered fermions,” Phys. Rev. D49, 486–494 (1994), arXiv:hep-lat/9306005 [hep-lat] .
  * Sharpe (1997) Stephen R. Sharpe, “Enhanced chiral logarithms in partially quenched QCD,” Phys. Rev. D 56, 7052–7058 (1997), [Erratum: Phys.Rev.D 62, 099901 (2000)], arXiv:hep-lat/9707018 .
  * Sharpe and Shoresh (2000) Stephen R. Sharpe and Noam Shoresh, “Physical results from unphysical simulations,” Phys. Rev. D62, 094503 (2000), arXiv:hep-lat/0006017 [hep-lat] .
  * Sharpe and Shoresh (2001) Stephen R. Sharpe and Noam Shoresh, “Partially quenched chiral perturbation theory without Phi0,” Phys. Rev. D64, 114510 (2001), arXiv:hep-lat/0108003 [hep-lat] .
  * Borsanyi  _et al._ (2013) Szabolcs Borsanyi, Stephan Durr, Zoltan Fodor, Stefan Krieg, Andreas Schafer, Enno E Scholz, and Kalman K Szabo, “SU(2) chiral perturbation theory low-energy constants from 2+1 flavor staggered lattice simulations,” Phys. Rev. D88, 014513 (2013), arXiv:1205.0788 [hep-lat] .
  * Colquhoun  _et al._ (2016) B. Colquhoun, R.J. Dowdall, J. Koponen, C.T.H. Davies, and G.P. Lepage, “$B\rightarrow\pi\ell\nu$ at zero recoil from lattice QCD with physical u/d quarks,” Phys. Rev. D 93, 034502 (2016), arXiv:1510.07446 [hep-lat] .
  * Dowdall  _et al._ (2019) R.J. Dowdall, C.T.H. Davies, R.R. Horgan, G.P. Lepage, C.J. Monahan, J. Shigemitsu, and M. Wingate, “Neutral B-meson mixing from full lattice QCD at the physical point,” Phys. Rev. D 100, 094508 (2019), arXiv:1907.01025 [hep-lat] .
  * Sharpe (2007) Stephen R. Sharpe, “Future of Chiral Extrapolations with Domain Wall Fermions,” in  _Workshop on Domain Wall Fermions at Ten Years_ (2007) arXiv:0706.0218 [hep-lat] .
  * Aoki  _et al._ (2011) Y. Aoki _et al._ (RBC, UKQCD), “Continuum Limit Physics from 2+1 Flavor Domain Wall QCD,” Phys. Rev. D83, 074508 (2011), arXiv:1011.0892 [hep-lat] .
  * Kass and Raftery (1995) Robert E. Kass and Adrian E. Raftery, “Bayes factors,” Journal of the American Statistical Association 90, 773–795 (1995).
  * (115) “LALIBE: CalLat Collaboration code for lattice qcd calculations,” https://github.com/callat-qcd/lalibe.
  * The HDF Group (1997-NNNN) The HDF Group, “Hierarchical Data Format, version 5,” (1997-NNNN), http://www.hdfgroup.org/HDF5/.
  * Kurth  _et al._ (2015) Thorsten Kurth, Andrew Pochinsky, Abhinav Sarje, Sergey Syritsyn, and Andre Walker-Loud, “High-Performance I/O: HDF5 for Lattice QCD,”  _Proceedings, 32nd International Symposium on Lattice Field Theory (Lattice 2014): Brookhaven, NY, USA, June 23-28, 2014_ , PoS LATTICE2014, 045 (2015), arXiv:1501.06992 [hep-lat] .