Emergent strongly coupled ultraviolet fixed point in four dimensions with 8 K\"ahler-Dirac fermions

The existence of a strongly coupled ultraviolet fixed point in 4-dimensional lattice models as they cross into the conformal window has long been hypothesized. The SU(3) gauge system with 8 fundamental fermions is a good candidate to study this phenomenon as it is expected to be very close to the opening of the conformal window. I study the system using staggered lattice fermions in the chiral limit. My numerical simulations employ improved lattice actions that include heavy Pauli-Villars (PV) type bosons. This modification does not affect the infrared dynamics but greatly reduces the ultraviolet fluctuations, thus allowing the study of stronger renormalized couplings than previously possible. I consider two different PV actions and find that both show an apparent continuous phase transition in the 8-flavor system. I investigate the critical behavior using finite size scaling of the renormalized gradient flow coupling. The finite size scaling curve-collapse analysis predicts a first order phase transition consistent with discontinuity exponent $\nu=1/4$ in the system without PV bosons. The scaling analysis with the PV boson actions is not consistent with a first order phase transition. The numerical data are well described by"walking scaling"corresponding to a renormalization group $\beta$ function that just touches zero, $\beta(g^2) \sim (g^2 - g^2_\star)^2$, though second order scaling cannot be excluded. Walking scaling could imply that the 8-flavor system is the opening of the conformal window, an exciting possibility that could be related to t'Hooft anomaly cancellation of the system.


INTRODUCTION
The existence of a conformal phase in 4-dimensional gauge-fermion systems is well established [1]. The upper limit of the fermion number where the conformal phase becomes infrared free can be predicted perturbatively. However, the lower limit where the QCD-like chirally broken and confining phase turns conformal requires a nonperturbative approach. The conformal phase is characterized by an infrared fixed point (IRFP) where the gauge coupling is irrelevant. In the "walking technicolor" scenario the conformal window opens up when the renormalization group (RG) β function vanishes just before chiral symmetry breaking could decouple the fermions [2][3][4][5][6]. This fixed point turns into an IRFP as the number of fermions is further increased, but at the sill the β function has a maximum that just touches zero. Alternatively, the IRFP could appear together with an ultraviolet fixed point (UVFP) as hypothesized by several authors [7][8][9][10]. At the sill of the conformal window the β function is similar to that of the walking scenario. The difference becomes evident in the conformal regime where the β function exhibits both an IRFP and a UVFP. Conformality is "lost" as the two fixed points merge and move into the complex plane. This latter scenario is particularly interesting as the emergence of a UVFP implies a second order phase transition with a new relevant operator. The corresponding phase diagram of the conformal regime was sketched in Ref. [11] and reproduced in the left panel of Fig.1. The sketch shows the phase structure and the RG flows in an extended parameter space that * Anna.Hasenfratz@colorado.edu includes a new relevant operator denoted by G. As the number of flavors increases the IRFP moves to weaker coupling. When the Gaussian FP (GFP) and the IRFP merge, the system becomes infrared free. On the other hand, when the number of flavors decreases, the IRFP approaches the UVFP. The opening of the conformal window is characterized by the merge of the UV and IR fixed points, as is sketched on the right panel of Fig.1.
The conformal IRFP does not separate phases. However, a UVFP would be associated with a continuous phase transition, indicated by the solid black curves in Fig.1. The phase transition might turn to first order before it intercepts the G = 0 action, if the two lines intercept at all (dashed black curves in Fig.1.). Phase transitions are much easier to identify in numerical simulations than RG fixed points. QCD-like, chirally broken lattice models often do not show any phase transition at zero temperature. When they do (e.g. with Wilson fermions), it is a bulk first order transition that arises from lattice artifacts. Conformal lattice systems, on the other hand, almost always exhibit a phase transition from the conformal phase at weak coupling to a confining one in the strong coupling regime. However, despite extensive numerical studies of systems with various fermion flavors, representations, and gauge colors, no continuous phase transition has been identified in four dimensional gauge-fermion systems at zero temperature. For a recent review see e.g. [12] and references therein.
In this paper, I present numerical results that indicate a continuous phase transition of an SU(3) gauge model with 8 fundamental massless fermions implemented as two sets of staggered lattice fields. This system has been investigated by several collaborations using different versions of staggered fermions [13][14][15][16][17][18][19][20][21][22][23][24]. Lattice results indicate that N f = 8 is close to the conformal sill, though Sketch of a hypothetical phase diagram of a conformal system with concomitant infrared and ultraviolet fixed points in an extended parameter space. g0 refers to the bare gauge coupling, while G denotes the coupling of the emergent relevant operator (possibly a 4-fermion Nambu-Jona-Lasinio interaction). The left panel, showing separate IR and UV fixed points in a conformal system, is reproduced from Ref. [11]. The right panel shows the "merging" of the two fixed points and corresponds to the sill of the conformal window. The relevant fermion mass is assumed to be zero, i.e. the sketches describe systems on the critical chiral surface.
there is no solid evidence that would show if it is inside or out of the conformal window. For example, finite temperature simulations at zero fermion mass do not identify a chirally broken phase before the appearance of a first order bulk transition [13,17,21]. Spectrum results in the weak coupling regime are consistent with dilaton chiral perturbation theory (dChiPT), which is designed to describe systems very close but below the conformal sill [23,[25][26][27][28]. However, dChiPT fits found that the extrapolation of the lattice data to the chiral limit is very long, and the dynamically-generated IR scale of the massless Nf=8 theory is much smaller than all the scales probed by simulations. It is also feasible that dChiPT describes systems very close but above the conformal sill [29].
Staggered fermions are equivalent to Dirac-Kähler fermions and correspond to four Dirac flavors per staggered species at the vanishing gauge coupling of the Gaussian FP [30][31][32]. The renormalized trajectory (RT) that originates from the GFP preserves this property, so N s staggered fermions describe 4N s Dirac flavors around the GFP, and even at the IRFP, if the system is conformal. Beyond the IRFP the equivalence does not have to hold.
Most of the existing studies with 2 staggered fields found a first order bulk transition, and many identified a novel strong coupling phase that exhibits spontaneous breaking of the single-site shift symmetry of staggered fermions (S4 phase) [33]. The S4 phase appears to be chirally symmetric and confining 1 , conceivably describing symmetric mass generation (SMG) [34][35][36]. Two staggered fermion species in the massless limit are equivalent to four reduced staggered fields that correspond to 16 1 The spectrum in the S4 phase has been studied extensively with N f = 12 flavors (3 staggered fermions) in Ref. [33]. In the N f = 8 case only the existence of the phase as evidenced by the S4 order parameter has been established. Simulations with finite fermion mass in the S4 phase are under way and publication describing the meson spectrum is in preparation.
Weyl spinors. This is exactly the fermion number needed to cancel all t'Hooft anomalies, a necessary condition for the fermions to become massive while remaining chirally symmetric. Thus the model could describe SMG if an appropriate four-fermion interaction is induced [35,36]. Gauge fields with momentum close to the cutoff might generate such interaction between the fermion doublers. This possibility has not been pursued since a universal continuum limit cannot be defined at first order phase transitions. The existence of a continuous phase transition would change the picture. The results I present here were obtained using uniquely improved lattice actions that contain heavy Pauli-Villars (PV) type bosonic fields [37], similar to the PV regulators in perturbation theory [38]. The mass of the PV bosons is always kept comparable to the cutoff, thus the PV fields do not influence the IR dynamics. However, they generate a local effective gauge action that shifts the bare gauge coupling towards weaker coupling where the UV fluctuations are suppressed. As the number of PV fields increases, or their mass decreases, the bulk phase transition is also shifted. Both in Ref. [37] and in the present work it is observed that the renormalized gauge coupling increases, while the discontinuity decreases at the phase transition as PV bosons are added. In the N f = 8 system the discontinuity at the phase transition eventually disappears, suggesting a continuous transition. This could imply that the first order transition observed without PV bosons at strong coupling is due to lattice fluctuations and opens the possibility that a non-perturbative, strongly coupled UVFP emerges as the cutoff effects are reduced by improved actions.
I investigate two lattice actions with different PV boson contents, in addition to the one without PV fields. I use finite size scaling to predict the critical coupling and the critical exponent at the phase transition. In the analysis I use the finite volume gradient flow renormalized coupling g 2 GF [39] at various flow times fixed by the lattice size L as √ 8t = cL, where c is an arbitrary parameter that I vary between 0.3 and 1.0, though the conclu-sion is based on the results at c = 0.45 [40]. Inspired by the arguments in Refs. [2,7,9], I also consider a scaling form corresponding to an RG β function that touches the axis at g 2 , i.e. β(g 2 ) ∼ (g 2 − g 2 ) 2 . This is reminiscent to the Berezinsky, Kosterlitz, Thouless (BKT) scaling of the 2 dimensional XY model, though with exponent ν = 1.0, instead of ν = 0.5 [41][42][43]. The finite size scaling analysis for the action without PV bosons is consistent with a first order transition with critical exponent ν close to 1/4. The PV improved actions show good curve collapse and consistent exponents for for "walking scaling" and ν ≈ 1.0, though second order scaling cannot be fully excluded either. However, numerical results of the PV improved actions are not consistent with a first order transition.
The rest of this paper is organized as follows. Sect. II A describes the lattice actions and the effect of the additional PV bosons. In Sect. II B I outline finite size scaling with the gradient flow coupling and discuss the relevant scaling forms. Sect. III A compares the phase structure of the PV improved actions and the original action without PV fields. Finally, I present the result of the finite size scaling analysis in Sect. III B. Sect. IV summarizes the results and discusses open questions.
In the simulations I set the fermion mass to am f = 0 and use symmetric L 4 volumes with L/a = 8, 10, 12, 16, and 20. The fermion boundary condition is antiperiodic in all four directions to control the fermion zero modes. The gauge fields obey periodic boundary conditions. Simulations with am f = 0 are well behaved in conformal systems or in the small volume deconfined regime of QCD-like systems. The S4 phase is confining but chirally symmetric, and the pions are not particularly light [33]. While numerical simulations are increasingly difficult at strong gauge coupling, the Hybrid Monte Carlo (HMC) updates are well controlled even in the chiral limit of the S4 phase.
For each staggered fermion field, I include N P V PV bosons. The PV bosons have the same action as the fermions but with bosonic statistics [37]. They are degenerate in mass and I choose their lattice mass am P V to be much larger than any IR scale of the system. When integrated out, the PV fields generate an effective gauge action. This is not an ultra-local action, but still local, with localization governed by the mass of the PV "pions". With heavy PV mass the bound state meson masses are approximately 2am P V , i.e. the range of the effective gauge action is a few lattice spacings if 2am P V 1.0. The leading term of the induced gauge action is a plaquette built from smeared links with coefficient where N s = 2 is the number of staggered fermion fields. Since β ind is negative, the bare gauge coupling β has to increase to keep the effective gauge coupling β eff = β + β ind and the lattice spacing constant. In general, more PV bosons and/or smaller am P V values generate a larger induced gauge action. The constant physics regime is pushed to weaker bare gauge coupling, which reduces the ultraviolet fluctuations as evidenced by larger plaquette expectation values. The effect of the PV bosons was investigated in detail in Ref. [37] for the 12-flavor SU(3) system. I want to emphasize that the effect of the PV bosons in the simulations is the generation of a local effective gauge action. The corresponding system has the same infrared properties at the Gaussian, as well as other possible fixed points as the original action, as long as am P V 1/L and am P V am f .

B. Finite size scaling with gradient flow
The nonperturbative UVFP in the 8-flavor SU(3) gauge system, if it exists, has at least two relevant parameters, the fermion mass and the new operator emerging at the UVFP. The scaling forms I consider below assumes the simplest scenario, i.e. that there are exactly two relevant parameters. If the fermion mass is set to zero, the RG flow is restricted to the critical surface of the massless theory. In that case, the UVFP has only one relevant operator, and finite size scaling methods can be used to study its properties. This is the scenario depicted in the sketches of Fig. 1.
Finite size scaling studies usually consider susceptibilities or bound state masses obtained from 2-point functions. Here I use the gradient flow coupling g 2 GF . The gradient flow (GF) [52,53]. is a smoothing transformation that has been used extensively in lattice studies GF is not an RG transformation as it lacks the essential coarse-graining step. However, coarse-graining can be implemented at the level of expectation values. If one relates the lattice GF time t/a 2 to the RG scale change as b ∝ t/a 2 , at large t/a 2 1 the GF can be interpreted as a well-defined real-space RG transformation [54]. Observables at finite flow time correspond to RG flowed observables at energy scale µ ∝ 1/ √ 8t, and exhibit hyperscaling in the vicinity of an RG fixed point.
The GF coupling was proposed in Ref. [39] for scale setting and in Ref. [40] to calculate the finite volume step scaling function in QCD-like or conformal system. It is defined as where E(t) is the energy density at bare gauge coupling β and flow time t on lattice volume (L/a) 4 2 . In lattice calculations various local operators like the plaquette or the clover operator can be used to approximate E. A frequent choice to estimate E is where Re tr U and Re tr V (t) denote the expectation value of the plaquette before GF flow and at flow time t.
The GF defines a non-linear RG transformation, therefore E(t) β,L/a does not have an anomalous dimension, i.e. the η exponent of the RG does not enter. Thus, the scaling dimension of the gradient flow coupling vanishes, its RG evolution depends only on the RG flow in the action parameter space. 3 . This property greatly simplifies the RG scaling relations, as only the scaling exponent related to the correlation length enters. One can think of g 2 GF as a scalar quantity that tracks the renormalized trajectory [11,59]. It has several advantages, most important is that it is a simple expectation value that is easy to measure with high precision. Finite size scaling relations for g 2 GF are derived in the usual way. At a second order phase transition one expects the correlation length to scale as where β denotes the critical coupling. The combination is a dimensionless scaling variable. Renormalization group scaling in finite volume therefore implies that g 2 GF (β, L; c), c = √ 8t/L, depends on (β/β − 1) L 1/ν , i.e. follows a unique but a priori unknown scaling function First order phase transitions are expected to show similar scaling behavior but with discontinuity exponent ν = 1/d = 0.25 in 4 dimensions.
where Nc refers to color, g 2 GF matches the M S coupling at tree level perturbation theory [39]. In the non-perturbative phase this normalization is arbitrary, but I keep it for consistency. 3 The use of t 2 (3 − Re tr V (t) ) to approximate g 2 GF is reminiscent of the operators used in Monte Carlo Renormalization Group (MCRG) studies [55][56][57][58]. The main difference is that GF is continuous, so one is not restricted to scale change of 2. Also, the lattice volume does not decrease during GF, so larger flow time/ RG scale change is possible.
If the RG β-function behaves as in the vicinity of the critical point, the phase transition is continuous but "walking" or BKT-type. The correlation length scales as and the corresponding scaling function depends on The 2-dimensional XY model has ν = 1/2, while the quadratic "walking" β function in Eq. 2.8 corresponds to ν = 1. Both Eqs. 2.7 and 2.10 are leading order relations valid only in the vicinity of the critical point, |β−β | 1. One also has to assure that the flow does not drive the RG blocked correlation length below the lattice spacing, t/a 2 ξ. At the same time t/a 2 has to be large enough for the RG flow to reach the renormalized trajectory. This latter condition could limit the use of the smallest volumes, since t = (cL) 2 /8 and c ≤ 0.5 is commonly preferred 4 . The RT is approached according to the largest non-leading exponent, but corrections to scaling could be important. This, however, is beyond what my present data set can describe.
The scaling functions f BKT depend on the parameter c, but should not depend on ν, ζ, and β . In addition, different lattice actions should show scaling with the same ν exponent, though β and ζ are dependent on the action. Thus testing different c values and different actions provides a consistency check and verifies the scaling form.

III. NUMERICAL RESULTS
A. The phase structure I investigated several different PV boson actions, two of them in detail. The first one contains 8 PV fields for each staggered flavor with mass am P V = 0.75 (8PV-m0.75 action), and the second one contains 4 PV fields per staggered flavor with mass am P V = 0.5 (4PV-m0.5 action). I have also looked at the action without PV bosons (0PV action). The 0PV action was studied previously in Ref. [17] at small but finite mass, and in the chiral limit but at weaker couplings in Ref. [20]. No previous work considered the 8-flavor systems at the phase 4 At larger c values the GF "wraps around" the finite volume.
That does not invalidate the use of g 2 GF , though the statistical errors increase with increasing c. A full investigation on the dependence on c, similar to Ref. [60], is a worthwhile project for the future.   2. The plaquette, the gradient flow coupling g 2 GF at c = 0.45, the S4 order parameter, and the topological susceptibility as the function of the bare gauge coupling β for the 0PV system. All simulations were done in the chiral limit on L 4 volumes. The dashed lines are guides to the eye. transition directly in the chiral limit. In contrast, all simulations I present here are performed with am f = 0.
I start by comparing the plaquette expectation value, the finite volume GF coupling, the S4 order parameter for the plaquette, and the topological susceptibility of the three actions. The top panel of Fig. 2 shows the plaquette as the function of the bare gauge coupling β on volumes L/a = 8, 10, 12, 16, and 20 for the 0PV action. The discontinuity at the phase transition is too small to be resolved on the plot. Nevertheless, I have observed tunneling and 2-state signals in several ensembles, consistent with a first order phase transition. The GF coupling g 2 GF at c ≡ √ 8t/L = 0.45 is shown on the second panel. It exhibits a rapid change around β ≈ 4.3−4.55 on these volumes, consistent with a phase transition around β 1st ≈ 4.6, the value predicted with small but finite fermion masses in Ref. [17]. The strong coupling phase exhibits spontaneous breaking of the staggered single-  site shift symmetry (S4 phase). The S4 order parameter that measures the staggered nature of the plaquette (Eq. 3. in Ref. [33]) is shown on the third panel. It is fairly noisy but useful to verify that the new phase shows the S4 symmetry breaking pattern. The S4 phase is confining but chirally symmetric, and the existence of a local order parameter ensures that it is separated from both the conformal and the chirally broken, confining phases by a phase transition [33] The bottom panel of Fig. 2 shows the topological susceptibility χ Q = Q 2 /L 4 . In the chiral am f = 0 limit of conformal and QCD-like systems one expects χ Q = 0, because the fermion determinant has a zero mode on isolated instantons. Staggered fermion lattice artifacts lift χ Q slightly at finite lattice spacing [61]. However, in the S4 phase, the topological charge is large and no longer suppressed. This is further evidence of the very different nature of this new phase. The largest GF coupling in the weak coupling phase of the 0PV model is g 2 GF 26 with L = 16 and c = 0.45. The plaquette is fairly small at the phase transition, Re tr U ≈ 1.0 (normalized to 3.0), suggesting large UV fluctuations. This system has another first order phase transition at an even stronger coupling that separates the S4 phase from a QCD-like chirally broken, confining regime [17]. I am not concerned with that transition in this work.
The panels of Figs. 3 and 4 show the same quantities as Fig. 2 but for the 8PV-m0.75 and 4PV-m0.5 actions and volumes L ≤ 24 5 . The GF coupling now appears continuous. None of my numerical simulations showed 2-state signals or phase tunneling either. The S4 order parameter becomes non-zero at the same place where g 2 GF shows a qualitative change of β dependence, around g 2 GF 35, c = 0.45. This is also where the S4 order parameter and topological susceptibility show a sudden increase. The plaquette is significantly larger at the phase transition, Re tr U ≈ 1.75 and Re tr U ≈ 1.7, respectively. The bare gauge coupling of the phase transition has also shifted to larger bare coupling, from β ≈ 4.6 of the 0PV action to β ≈ 8.7 and β ≈ 8.1. This is the consequence of the induced effective action.
The strongest GF coupling with the 0PV action is g 2 GF ≈ 26 on the volumes considered, compared to g 2 GF ≈ 35 with both PV actions. The 0PV action has large vacuum fluctuations. It appears that those trigger a first order phase transition before the GF coupling is strong enough to undergo the continuous transition seen with the 4PV-m0.5 and 8PV-m0.75 actions. This suggests that the first order phase transition could be only a lattice artifact. The picture is similar to the frequently observed first order bulk transition with Wilson fermions that are not considered physical.
If the phase transition is continuous, it should follow a scaling form like Eq. 2.7 or Eq. 2.10. A first order phase transition should follow the scaling of Eq. 2.7 with exponent ν = 1/4. The critical coupling β , exponent ν and parameter ζ can be determined based on a curve collapse analysis of g 2 GF (x; c) vs. x = |β − β| ν L or x = Lexp (−ζ|β/β − 1| −ν ), respectively. I perform this analysis similar the Nelder-Mead method as described in [62].
Specifically, I chose a reference volume L 0 and interpolate g 2 GF (β, L 0 ; c) versus β with a (smoothed) cubic spline. In the analysis below I use L 0 = 16. I fit the 4 2nd order fit parameters β , ν (and ζ) by minimizing the deviation of g 2 GF (β, L; c) from this interpolating curve at β 0 , where β 0 is the solution of the equation if fitting Eq. 2.7 or 2) if fitting Eq. 2.10 The procedure can be iterated by refining the "ideal" interpolating curve using more volumes, once the fit parameters are known approximately. In addition, I excluded those data points that are not clearly in the S4 phase 6 . However, for completeness, I show the smaller volumes and excluded data points in 6 At a first order phase transition the critical coupling of an MCRG approach can overshoot the discontinuity, especially if the transition is strongly first order [63,64]. In order to avoid data that are in the weak coupling phase I employ a conservative cut based on the S4 order parameter.   Fig. 2 and the observed phase transition in Ref. [21]. The bottom two panels show the curve collapse for the two PV improved actions. Again, only L ≥ 16 are included in the fit. In all cases χ 2 /dof is O(1), but the predicted exponent ν varies between 0.53 and 0.63, very different from the discontinuity exponent of a first order phase transition.

Second order scaling test
To test the consistency of the curve collapse, I repeat the analysis varying the parameter c = √ 8t/L between 0.3 and 1.0. I also compare fits using volumes L ≥ 16 or L ≥ 20 for the PV actions. Fig. 6 summarizes the predicted ν and corresponding β values for the three actions. The observable drift at smaller c indicates that the RG flow has not yet reached the vicinity of the renormalized trajectory. The 0PV action stabilizes at c > 0.5 and predicts ν = 0.25(1), the expected scaling exponent at a first order transition. The exponent ν for the PV improved actions drift and show significant variation depending on the action even at c ≈ 0.9. I do not find a consistent scaling fit with volumes L ≥ 16, though it is possible that on larger volumes the curve collapse analysis would predict a consistent exponent, possibly around ν = 0.5. For both PV improved actions, a first order transition is not consistent with the numerical data. The "walking" or BKT-like scaling form given in Eq. 2.10 assumes that the RG β function just touches zero.
This fit form has three free parameters, though ν = 1.0 is the theoretically motivated value [7]. Using L ≥ 16 volumes I am able to fit all three parameters. Fits with BKT-like scaling have smaller χ 2 /dof than the second-order scaling. The predicted ν ≈ 1.0 exponent and no significant dependence on c also favors the "walking" scaling scenario. However the second order scaling form discussed in Sect. III A 1 is not excluded. It is well known that corrections to BKT scaling in the 2dimensional XY model are very large [65]. In the 8-flavor system I have concentrated on fitting the leading terms only. More sophisticated curve collapse fitting, better statistics, and especially larger volumes might allow distinguishing the two cases in the future [66].

IV. SUMMARY
I have investigated the SU(3) gauge model with two sets of staggered fermions in the chiral limit. Staggered lattice fermions are equivalent to Dirac-Kähler fermions, and in the chiral limit two staggered fields correspond to 16 Weyl spinors. This system could lead to symmetric mass generation if appropriate terms are induced by the strong interaction.
I study two improved lattice actions that contain heavy Pauli-Villars type bosonic fields. Both improved actions shift the first order phase transition of the original system to weaker bare couplings and smooth out the transition. Using finite size scaling methods I showed that the phase transition with both PV improved actions appears con-  tinuous. I investigated scaling according to standard second order phase transition, and also "walking" or BKTlike scaling. First order phase transition with ν = 1/4 is not consistent with the numerical data. While I cannot exclude second order scaling, the curve collapse fit favors "walking" scaling with ν ≈ 1.0. A "walking" scaling phase transition could imply that the system is at the sill of the conformal window as sketched on the right panel of Fig. 1. If that is indeed the case, it is likely that the complete t'Hooft anomaly cancellation of the lattice formulation plays an essential role.
The weak coupling phase with all three actions appears chirally symmetric in the volumes considered, i.e. the simulations are well controlled even in the am f = 0 limit with meson spectrum showing chiral symmetry. Simulations at finite fermion mass that contrast the δ-regime rotator spectrum with the mass-deformed hyperscaling prediction might be required to establish the infrared properties of the weak coupling phase [16,67]. In this work I concentrate on the order of the phase transition from the strong coupling side and do not directly investigate the infrared properties of the weak coupling phase.
The strong coupling regime breaks the single site shift symmetry of the staggered action. The S4 phase could describe symmetric mass generation, where confinement leads to mass generation without breaking chiral symmetry. The existence of a continuous phase transition would allow taking an infinite cutoff continuum limit. This very exciting possibility requires further analytical and numerical studies.
Lattice studies with N f = 12 fundamental flavors also show an S4 phase in the strong coupling. However, my preliminary numerical studies with various PV actions did not show continuous phase transition in that system. t'Hooft anomalies do not cancel with 12 flavors (i.e. three staggered fermions), possibly explaining the difference between 12 and 8 flavors. However, further studies are needed.
If the model with 2 staggered fermions is conformal in the weak coupling regime, 8 fundamental Dirac fermions should also be conformal. The critical properties at the continuous phase transition, however, are not necessarily universal. It would be very interesting to investigate 8 fundamental flavors with improved domain wall fermions that allow the study of the strong gauge couplings where the continuous phase transition with staggered fermions occur.

ACKNOWLEDGMENTS
The phase diagram depicted in Fig. 1 came from a discussion with S. Rychkov. I am indebted for his insight that lead to the present investigation. I am grateful to Yigal Shamir for his constructive criticism and suggestions, and to Simon Catterall for discussions about t'Hooft anomaly cancellation and SMG. I have presented preliminary results to the LSD collaboration and benefited from their probing questions when finalizing the manuscript. I thank Oliver Witzel who carefully read an earlier version of this manuscript for his many comments and suggestions.
Computations for this work were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy, the RMACC Summit supercomputer [68], which is supported by the National Science Foundation (awards No. ACI-1532235 and No. ACI-1532236), the University of Colorado Boulder, and Colorado State University, and on the University of Colorado Boulder HEP Beowulf cluster. I thank the University of Colorado for providing the facilities essential for the completion of this work. My reseach is supported by DOE grant DE-SC0010005.