Infrared fixed point in the massless twelve-flavor SU(3) gauge-fermion system

We present strong numerical evidence for the existence of an infrared fixed point in the renor-malization group flow of the SU(3) gauge-fermion system with twelve massless fermions in the fundamental representation. Our numerical simulations using nHYP-smeared staggered fermions with Pauli-Villars improvement do not exhibit any first-order bulk phase transition in the investigated parameter region. We utilize an infinite volume renormalization scheme based on the gradient flow transformation to determine the renormalization group β function. We identify an infrared fixed point at g 2GF ⋆ = 6 . 69(68) in the GF scheme and calculate the leading irrelevant critical exponent γ ⋆g = 0 . 204(36). Our prediction for γ ⋆g is consistent with available literature at the 1-2 σ level.

Most lattice studies are affected by the presence of a bulk first-order phase transition [41][42][43][44][45][46].Such unphysical phase transitions are triggered by strong ultraviolet fluctuations in the fermion sector that prevent lattice simulations from reaching deep into the infrared regime.Even when strong couplings are reached, lattice cutoff effects make it difficult to take the proper continuum limit, leading to inconsistent results between different lattice formulations.It is imperative to reduce the ultraviolet fluctuations that trigger first-order bulk phase transitions.Ref. [47] suggested to include unphysical heavy Pauli-Villars fields to achieve the necessary improvement.Lattice Pauli-Villars fields are similar to their continuum analogue -they have the same action as the fermions, but possess bosonic statistics.Their mass is at the level of the cutoff.Therefore, they decouple in the infrared limit, γ g = 0.199 (32) Non-pert.continuum 1-loop univ.

3-loop GF
FIG. 1.Our continuum prediction for βGF g 2 GF as a function of g 2 GF (maroon band) juxtaposed against the 1-(dashed), 2-(dotted) and 3-loop (dashed-dotted) gradient flow β function from perturbation theory [50].The width of the maroon band indicates the error.Also given is our prediction for the leading irrelevant critical exponent at the IRFP, γ ⋆ g = 0.199 (32).
while in the ultraviolet they compensate for cutoff effects introduced by the fermions.This idea has been tested in simulations of the SU(3) gauge-fermion system with ten massless fundamental Dirac fermions (flavors) and the SU(4) gauge-fermion system four massless fundamental and four massless two-index Dirac fermions [48,49].Both studies extended the reach into the infrared regime significantly, and both present clear evidence for infrared conformality in those systems.See also Refs.[44,45] for an alternative proposal to remove unphysical bulk phase transitions.
In this work, we utilize Pauli-Villars (PV) improvement to study the infrared properties of the massless SU(3) gauge-fermion system with N f = 12 fundamental flavors.We calculate the renormalization group (RG) β function, defined as the logarithmic derivative of the renormalized running coupling g 2 (µ) with respect to an energy scale µ as If the β function exhibits a fixed point at some coupling arXiv:2402.18038v2[hep-lat] 13 Mar 2024 g 2 ⋆ , the gauge coupling becomes irrelevant and the system is infrared conformal.The value of g 2 ⋆ depends on the renormalization scheme, but the existence of the infrared fixed point (IRFP) does not.
On the lattice, it is convenient to use gradient flow (GF) transformation [51][52][53] as a continuous smearing operation to define an infinite-volume renormalization scheme [28,[54][55][56][57][58].In the GF scheme, the renormalized coupling g 2 GF (t) in infinite volume at flow time t ∝ 1/µ 2 is defined in terms of the Yang-Mills energy density E(t) as where N = 128π 2 /(3N 2 c − 3) is a constant chosen to match g 2 GF to g 2 MS at tree-level [52].The corresponding RG β function is To calculate the infinite volume gradient flow β function β GF g 2 GF from finite-volume simulations, we utilize the continuous β function method (CBFM) proposed in Refs.[28,54,57] and deployed extensively in Refs.[48,49,[59][60][61][62] to a variety of strongly-coupled gauge-fermion systems.In this paper we follow the steps described in [60,62] with additional extensions for improved error estimation.We discuss the continuous β function method and its implementation in further detail in Sec.III.
As a preview, we show our nonperturbative prediction for β GF g 2 GF as a function of g 2 GF in Fig. 1.The predicted β function converges to the universal 1-/2-loop and 3-loop gradient flow perturbative β functions at small g 2 GF [50].Around g 2 GF⋆ = 6.60 (62), the nonperturbative β function unambiguously exhibits an infrared fixed point.From the slope of β GF g 2 GF at g 2 GF⋆ , we calculate the leading irrelevant critical exponent γ ⋆ g = 0.199 (32) and find that it is consistent with the perturbative calculations of Refs.[4,6] at the 1σ level and the lattice calculation of Ref. [21] at the 2σ level.We control for systematic errors in the infinite volume extrapolation step of the CBFM using Bayesian model averaging.Additionally, our Pauli-Villars improved simulations offer tight control over systematics in the continuum extrapolation step of the CBFM.Our data is publicly available at Ref. [63].
This paper is laid out as follows.In Sec.II, we summarize details of our numerical simulations.In Sec.III, we review the continuous β function method and discuss our analysis.We explain our calculation of the leading irrelevant critical exponent in Sec.IV, and wrap up in Sec.V with conclusions.

II. NUMERICAL DETAILS
We simulate the massless twelve-flavor SU(3) gaugefermion system using an adjoint-plaquette gauge action with β F /β A = −0.25 (β F ≡ β b ) and a massless (am f = 0.0) nHYP-smeared staggered fermion action with four massive (am PV = 0.5) Pauli-Villars (PV) fields per staggered fermion [47,[64][65][66].The "pions" of these PV fields have mass am (PV) PS ≳ 1.03 and generate gauge loops in the effective gauge action with a size that decays exponentially with am (PV) PS [47].As long as the volume is much larger than 1/m the PV fields decouple in the infrared.Their only effect is a modified, but local, gauge action.One of the goals of the present work is to illustrate the validity of this expectation.
We use antiperiodic boundary conditions in all four directions for both the staggered fermion fields and PV fields.Our numerical simulations are performed using the hybrid Monte Carlo algorithm [67] implemented in a modified version of the MILC library4 and the Quantum EXpressions (QEX) library5 [68].We set the molecular dynamics trajectory length to τ = 1.0.Our configurations are separated by ten trajectories (10 molecular dynamics time units).We perform our simulations at fourteen bare gauge couplings (9.20 ≤ β b ≤ 14.6) and five symmetric volumes (24 ≤ L/a ≤ 40).In Table I, we list the total number of thermalized configurations on each ensemble.
Our gradient flow measurements are performed using either the modified MILC or QEX libraries [68].We flow our configurations using Wilson flow [52,53], integrating the gradient flow equations using the 4th-order Runge-Kutta algorithm discussed in Ref. [52] with time step dt/a 2 = 0.02 for t/a 2 ≤ 5.0 and dt/a 2 = 0.1 for t/a 2 > 5.0.At each integration step, we measure the Yang-Mills energy density E(t) using the Wilson (W) and clover (C) discretizations.In the rest of this paper, we refer to results based on Wilson flow and Wilson operator as WW, while we refer to results based on Wilson flow and clover operator as WC.Our data for the Yang-Mills energy density E(t) from the Wilson and clover operator is available at Ref. [63].
In Fig. 2 we plot the expectation value of the magnitude of the gradient-flowed Polyakov loop |P | at 8t/a 2 ≈ (L/2a) 2 against the bare gauge coupling β b for all ensembles in Table .I. Errors are estimated using the "Γmethod" technique implemented by the pyerrors package [69][70][71][72][73]; however, they are likely underestimated.Nonetheless, the Polyakov loop suggests that the system is not confining in the range of couplings and volumes that we use in the present study.

III. NONPERTURBATIVE β FUNCTION
We measure the gradient flow coupling in finite volume in terms of the Yang-Mills energy density E(t) as where δ(t, L) corrects for gauge zero modes [74].From g 2 GF t; L, g 2 0 , we calculate the gradient flow β function in finite volume as where we discretize d/dt with a 5-point stencil.The autocorrelation time for g 2 GF t; L, g 2 0 and β GF t; L, g 2 0 is typically between 20-80 molecular dynamics time units (MDTUs), with occasional jumps to 120-200 MDTUs.
To extract the continuum β GF g 2 GF as a function of g 2 GF , we follow the CBFM procedure outlined in Refs.[60,62].
1. Take the infinite volume limit by independently extrapolating both g 2 GF t; L, g 2 0 and β GF t; L, g 2 0 linearly in (a/L) 4 → 0 at fixed t/a 2 and β b .
3. Take the continuum limit by extrapolating β GF t; g 2 0 linearly in a 2 /t → 0 at fixed g 2 GF (t).
Correlated uncertainties are propagated throughout our analysis using the automatic error propagation tools provided by the gvar library [75].Fits are performed using the SwissFit library, which integrates directly with gvar [76].The steps of the CBFM are detailed in the rest of this section.

A. Infinite volume extrapolation
Because the Yang-Mills energy density is a dimension-4 operator, leading finite-volume corrections to g 2 GF t; L, g 2 0 and β GF t; L, g 2 0 are expected to be O t 2 /L 4 .Therefore, we extrapolate both g 2 GF t; L, g 2 0 and β GF t; L, g 2 0 to a/L → 0 by independently fitting them to the ansatz at fixed t/a 2 and β b .This analysis strategy was first outlined in Ref. [62] and subsequently applied in Refs.[60,61].Alternative methods are discussed in Refs.[28,54,57,59].We account for the systematic uncertainty that is associated with choosing a particular subset of volumes for the infinite volume extrapolation using Bayesian model averaging [77][78][79].We do so by first fitting over all possible subsets of volumes L/a ∈ {24, 28, 32, 36, 40} with at least three volumes in each subset η.We calculate the model weight w η for a particular subset η as where χ 2 η is the χ 2 statistic of fit η and d η is the number of data points not included in fit η from the full set of volumes {24, 28, 32, 36, 40}.Denoting the mean of k i from fit η as k (η) i , our model-averaged prediction for the mean k i of k i is (a/L) 4 3. Result of our infinite volume extrapolation of g 2 GF t; L, g 2 0 (left panels) and βGF t; L, g 2 0 (right panels) for the Wilson (W) operator at β b = 9.60 (top panels), 9.80 (middle panels) and 10.2 (bottom panels).Black (×) markers with error bars are the data included in our extrapolation.Extrapolations with errors that are predicted from Bayesian model averaging are indicated by multi-colored bands at t/a 2 = 2.5 (red), 3.5 (light green), 4.5 (cyan), and 6.0 (light purple).We do not show the infinite volume extrapolation of βGF t; L, g 2 0 at t/a 2 = 4.5 for visualization purposes.
where the weights w η have been normalized such that where i } i=1,2 from fit η.In Figs. 3 and 4, we show the result of our modelaveraged infinite volume extrapolation for the W and C discretization of E(t), respectively, over a range of flow times 2.5 ≤ t/a 2 ≤ 6.0 (different colors).The left panels of Figs.3-4 show our infinite extrapolation of g 2 GF t; L, g 2 0 , while the right panels show our infinite volume extrapolation of β GF t; L, g 2 0 .The bare gauge couplings that we chose for these plots are in the vicinity where the continuum β function predicts an IRFP.For all three bare gauge couplings shown in Figs.3-4, the model average is dominated by subsets containing L/a ∈ {28, 32, 36, 40}, as L/a = 24 often deviates from the linear trend in a 4 /L 4 , particularly as the flow time increases.This is reflected in the model average, as fits including L/a = 24 possess small model weights w η and contribute negligibly to the model average.Result of our infinite volume extrapolation of g 2 GF t; L, g 2 0 (left panels) and βGF t; L, g 2 0 (right panels) for the clover (C) operator at β b = 9.60 (top panels), 9.80 (middle panels) and 10.2 (bottom panels).Black (×) markers with error bars are the data included in our extrapolation.Extrapolations with errors that are predicted from Bayesian model averaging are indicated by multi-colored bands at t/a 2 = 2.5 (red), 3.5 (light green), 4.5 (cyan), and 6.0 (light purple).We do not show the infinite volume extrapolation of βGF t; L, g 2 0 at t/a 2 = 4.5 for visualization purposes.

B. Intermediate interpolation
The continuum limit a 2 /t → 0 of β GF t; g 2 0 is taken at fixed g 2 GF .We predict pairs t/a 2 , β GF t; g 2 0 at a set of fixed g 2 GF for the continuum extrapolation by interpolating β GF t; g 2 0 in g 2 GF t; g 2 0 at fixed t/a 2 using the ansatz At each t/a 2 , we account for the uncertainty in g 2 GF t; g 2 0 by including the mean and covariance of g 2 GF t; g 2 0 as a Gaussian prior.We also set a Gaussian prior on each coefficient p n with zero mean and a width of 0.1, which helps stabilize the fit.We choose N = 4, as it is the lowest value of N that fits the data well.In Fig. 5  are indicated by colored bands, with t/a 2 = 2.5 (red), 3.5 (light green), 4.5 (cyan), and 6.0 (light purple).The width of the band indicates the error.The data contributing to each interpolation is indicated by an open circular marker with both xand y-errors.We compare our interpolation against the continuum 1-(dashed), 2-(dotted), and 3-loop (dashed-dotted) gradient flow β function from perturbation theory [50].
at weak/strong coupling.Therefore, use N = 4 for our central analysis.We will discuss the systematic effect that is associated with the order N in our estimate of g 2 GF⋆ and γ ⋆ g in Sec.IV.

C. Continuum extrapolation
The final step is the continuum (a 2 /t → 0) limit over a set of fixed g 2 GF that predicts β GF (g 2 GF ) as a function of g 2 GF .The range of [t min /a 2 , t max /a 2 ] used in the continuum extrapolation must be chosen with care.The value of t min /a 2 must be large enough for the RG flow to reach the renormalized trajectory.Once this is the case, finitecutoff effects are O(a 2 /t).In practice, one can identify when t min /a 2 is close enough to the renormalized trajectory by the overlap between the continuum prediction for β GF (g 2 GF ) from both operators.Because finite-volume effects are expected to be O t 2 /L 4 , t max /a 2 must be chosen such t/a 2 < t max /a 2 has a reliable infinite volume  WW WC 1-loop univ.2-loop univ.
3-loop GF step-scaling [21] FIG. 7. Our continuum prediction for βGF g 2 GF as a function of g 2 GF for the W operator (gold band) and C operator (maroon band).The width of the band indicates the error.The nonperturbative results are juxtaposed against the 1-(dashed), 2-(dotted), and 3-loop (dashed-dotted) gradient flow β function from perturbation theory [50].Also shown is the step-scaling β function in the c = 0.25 scheme from Ref. [21] as a grey band.
extrapolation.Ideally, we would apply Bayesian model averaging to the continuum extrapolation to automatically account for systematic effect that is associated with making a particular choice in [t min /a 2 , t max /a 2 ].However, at this point in the analysis, we no longer have access to the full covariance matrix, which means that we no longer have access to a reliable estimate of the model weights.To estimate the error in our continuum extrapolation, we use the half-difference of the prediction from the continuum extrapolation performed at ±1σ.This approach was also taken in Ref. [60].
Fig. 6 shows examples of the continuum extrapolation performed in the range 2.0 ≤ g 2 GF ≤ 6.0 (different colors) using [t min /a 2 , t max /a 2 ] = [3.5, 6.0].For t/a 2 ≲ 3.5, β t; g 2 0 has a slight curvature in a 2 /t, indicating emerging higher-order cutoff effects for both the W and C operator.For t/a 2 ≳ 6.0, the data begins to deviate from a linear trend in a 2 /t, indicating that the infinite volume extrapolation is getting unreliable.Our choice of [t min /a 2 , t max /a 2 ] = [3.5, 6.0] avoids both of these two regimes.In Sec.IV, we discuss the sensitivity of our prediction for g 2 GF⋆ and γ ⋆ g to our choice of t min /a 2 , t max /a 2 .We show our prediction for the continuum β GF g 2 GF in Fig. 7 from the WW (gold band) and WC (maroon band) combination.The continuum predictions from both operators are consistent with one another across the entire range of investigated renormalized couplings g 2 GF .At small g 2 GF , the continuum β GF g 2 GF appears to converge to the 1-, 2-, and 3-loop perturbative gradient flow β function [50].At g 2 GF⋆ ≈ 6.60, our continuum β function predicts an infrared fixed point.The location of the fixed point is slightly below the predicted IRFP from the stepscaling calculation of Ref. [21].Note that, because the calculation in Ref. [21] was done in a different gradientflow-based renormalization scheme, the predicted g 2 GF⋆ values do not have to agree.Our final result for β GF g 2 GF from both operators is provided as an ASCII file.

IV. THE IRFP AND ITS LEADING IRRELEVANT CRITICAL EXPONENT
In the vicinity of the RG fixed point where γ ⋆ g is the universal critical exponent of the irrelevant gauge coupling.The factor of 1/2 is chosen to match the convention of Refs.[6,21].We estimate g 2 GF⋆ and γ ⋆ g via the following procedure.
1. Interpolate the central value of β g 2 GF and the central value of β g 2 GF ± 1σ in g 2 GF using a monotonic spline.
2. Estimate the central value of g 2 GF⋆ from the root of the spline interpolation of β g 2 GF in g 2 GF .
3. Estimate the central value of γ ⋆ g from the derivative of the spline at g 2 GF⋆ .4. Repeat Steps (2)-( 3) with a spline interpolation of β g 2 GF ± 1σ in g 2 GF .
5. Estimate the error in g 2 GF⋆ and γ ⋆ g from the half difference of their predictions from the interpolations in Step (4).
Steps (1)-( 5) yield g 2 GF⋆ = 6.69 (35), 6.60 (36) and γ ⋆ g = 0.206 (19), 0.199 (18) for WW and WC, respectively.In Fig. 8 we look at how the continuum limit predictions for g 2 GF⋆ (top panels) and γ ⋆ g (bottom panels) vary with our choice of t min /a 2 (x-axes) and t max./a 2 (different colors) for the W (left panels) and C (right panels) operators.The central values for both quantities and both operators are stable; they vary well within error.The stability in t min /a 2 is attributed to the linearity of the continuum extrapolation over a wide range of a 2 /t, while the stability in t max /a 2 is likely attributed to our control over the infinite volume extrapolation.
We take the result for γ ⋆ g from the WC combination with the value for [t min /a 2 , t max /a 2 ] = [3.5, 6.0] from Sec. V as our central result.We estimate additional systematic errors by varying our analysis as follows.
• Choosing a higher-order polynomial for the intermediate interpolation in Sec.III B. The highestorder interpolation in N that we can use before we lose control over our continuum extrapolation due to overfitting is N = 6.This shifts the value of γ ⋆ g by ≈ 0.001, and we take the latter difference as an estimate for the systematic error that is associated with our choice of N for the intermediate interpolation.
• Choosing a different t min /a 2 , t max /a 2 in the continuum extrapolation.We estimate the systematic error that is associated with our choice of t min /a 2 , t max /a 2 by the difference in the most extreme values of γ ⋆ g in our variations illustrated in Fig. 8.This yields a systematic error of ≈ 0.006.
• Choosing instead the prediction for γ ⋆ g = 0.206(19) from the WW combination as our central result.We take the difference in these predictions (≈ 0.007) as an estimate of the systematic error associated with making a particular choice of flow/operator combination.
To be conservative, we combine the error in our analysis of γ ⋆ g with the systematic error estimates above linearly.This yields the final of prediction γ ⋆ g = 0.199 (32).Repeating the same exercise for g 2 GF⋆ yields a systematic error of ≈ 0.12 from the interpolation order, ≈ 0.05 from the continuum extrapolation, and ≈ 0.09 from the flow/operator combination.Including the systematic error in g 2 GF⋆ linearly yields a final prediction of g 2 GF⋆ = 6.60 (62) [21] (teal error bar) and Ref. [6] (dark gold error bar).The smaller error bar on our result indicates the error without accounting for systematic effects; the larger error bar indicates our error after accounting for systematic effects.We indicate our total error with a grey band for visualization.
In Fig. 9 we compare our prediction for γ ⋆ g with those available in the literature [6,21].Our result is plotted as a maroon star with errors indicated by an error bar.The smaller error bar is our error estimate before accounting for systematic effects and the larger error bar includes systematic effects.The result for γ ⋆ g from the perturbative calculation of Ref. [6] (dark gold error bar) is within 1σ of our estimate for γ * g .The lattice calculation of γ ⋆ g from Ref. [21] (cyan error bar) is within 2σ of our result.However, it is important to note that the lattice calculation in Ref. [21] uses smaller volumes and very coarse lattices in comparison to the present work.It is also worth noting the "scheme-independent" prediction of γ ⋆ g ≈ 0.228 from Ref. [4] is also within 1σ of our predicted γ ⋆ g ; we don't show this result in Fig. 9 because no estimate of the systematic error in this result is available.

V. CONCLUSIONS
We have calculated the non-perturbative β function of the SU(3) gauge-fermion system with twelve massless fundamental fermions using a Pauli-Villars improved lattice action.We find strong evidence for an infrared fixed point at g 2 GF⋆ = 6.60(62) from our gradient-flow-based renormalization scheme.Our study utilizes a wide range of couplings and volumes.In particular, we can reach renormalized gauge coupling values well above the predicted IRFP without the interference of a bulk phase transition.We include systematic effects from the infinite volume extrapolation directly into our analysis using Bayesian model averaging.Our data exhibits cutoff effects that are consistent with the leading O(a 2 /t) form over a wide range t/a 2 .The consistency between the W and C operators further supports the leading order scaling behavior.We believe the improved scaling is due to the additional PV bosons that reduce cutoff effects.In contrast, we found significantly larger cutoff effects when we re-analyzed the data that were generated without PV fields and used in Ref. [21].Our data is publicly available at Ref. [63].
Overall, the systematics of our continuum extrapolation are well controlled.Based on the continuum prediction for β g 2 GF in g 2 GF , we estimate the leading irrelevant critical exponent γ ⋆ g = 0.199 (32).This estimate includes conservative systematic errors from various choices in our analysis.Our result for γ ⋆ g agrees with Refs.[4,6] at the 1σ level and Ref. [21] at the 2σ level.
FIG.3.Result of our infinite volume extrapolation of g 2 GF t; L, g 2 0 (left panels) and βGF t; L, g 2 0 (right panels) for the Wilson (W) operator at β b = 9.60 (top panels), 9.80 (middle panels) and 10.2 (bottom panels).Black (×) markers with error bars are the data included in our extrapolation.Extrapolations with errors that are predicted from Bayesian model averaging are indicated by multi-colored bands at t/a 2 = 2.5 (red), 3.5 (light green), 4.5 (cyan), and 6.0 (light purple).We do not show the infinite volume extrapolation of βGF t; L, g 2 0 at t/a 2 = 4.5 for visualization purposes.
FIG.4.Result of our infinite volume extrapolation of g 2 GF t; L, g 2 0 (left panels) and βGF t; L, g 2 0 (right panels) for the clover (C) operator at β b = 9.60 (top panels), 9.80 (middle panels) and 10.2 (bottom panels).Black (×) markers with error bars are the data included in our extrapolation.Extrapolations with errors that are predicted from Bayesian model averaging are indicated by multi-colored bands at t/a 2 = 2.5 (red), 3.5 (light green), 4.5 (cyan), and 6.0 (light purple).We do not show the infinite volume extrapolation of βGF t; L, g 2 0 at t/a 2 = 4.5 for visualization purposes.

FIG. 6 . 2 0
FIG. 6. Illustration of our continuum extrapolation of βGF t; g 2 0 at fixed g 2 GF = 2.0 (teal), 4.0 (dark orange), 6.0 (magenta), and 8.0 (forest green).Data contributing to our extrapolation with the W operator are shown as error bars with triangular markers and the C operator are shown as error bars with circular markers.Our extrapolations are shown as colored bands, where the error is indicated by the width of the band. 2 .