Asymptotic safety in the Litim-Sannino model at four loops

We consider a four-dimensional $SU(N_c)$ gauge theory coupled to $N_f$ species of color fermions and $N_f^2$ colorless scalars. The quantum field theory possesses a weakly interacting ultraviolet fixed point that we determine from beta functions computed up to four-loop order in the gauge coupling, and up to three-loop order in the Yukawa and quartic scalar couplings. The fixed point has one relevant direction giving rise to asymptotic safety. We compute fixed point values of dimensionless couplings together with the corresponding scaling exponents up to the first three nontrivial orders in Veneziano parameter $\epsilon$, both for infinite and finite number of colors $N_c$. We also consider anomalous dimensions for fields, scalar mass squared, and a class of dimension-three operators. Contrary to previous studies, we take into account possible mixing of the latter and compute eigenvalues of the corresponding matrix. Further, we investigate the size of the conformal window in the Veneziano limit and its dependence on $N_c$.


I. INTRODUCTION
The study of asymptotic behavior of the dimensionless couplings in quantum field theory (QFT) provides important information both for the Standard Model (SM) and Beyond the Standard Model (BSM) scenarios.One of these behaviors is known as asymptotic freedom [1,2], which is a defining feature of quantum chromodynamics.This behavior entails a decrease in the value of a coupling with the energy scale.Thus, in the deep ultraviolet (UV), this coupling tends to approach the Gaussian noninteractive fixed point (FP).
Asymptotic safety (AS) is an extension of the concept of asymptotic freedom, as outlined in the work of S. Weinberg [3].In AS, the coupling in the deep UV also reaches a fixed point, but unlike in asymptotic freedom, the fixed point value is not zero, which means that the theory remains interactive.Such theories are referred to as asymptotically safe.
The concept of asymptotic safety was initially introduced by S. Weinberg in the late 1970s as a means of achieving nonperturbative renormalizability for the four-dimensional theory of gravity [3].However, in recent years, AS has been widely applied in the context of gauge theories to address issues with U (1) gauge couplings (Landau pole) by stabilizing them at an interactive fixed point at a certain scale.Indications of AS has been found in simple [4,5], semisimple [6], and supersymmetric gauge theories coupled to matter [7][8][9].From the phenomenological point of view, properties of such UV interactive fixed points for matter fields can be transmitted down to the low-energy regime and lead us to some phenomenological predictions, see, e.g., recent reviews [10,11].
An example of a UV-complete particle theory with a weakly interacting fixed point is a model of N f fermions coupled to SU (N c ) gauge fields and elementary scalars through gauge and Yukawa interactions [5,12].In the large-N Veneziano limit, the fixed point can be systematically studied in perturbation theory using a small control parameter ϵ, allowing for the extraction of specific details of theory.Previous studies [5,[13][14][15] have identified critical couplings and universal exponents up to second order in ϵ, including finite N corrections.Phenomenological applications of the model and its extensions can be found in Refs.[15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33].It is also worth mentioning that the model was the first nonsupersymmetric theory investigated in the large global charge limit [34,35].The holographic description of the model was considered in Ref. [36].
In this paper, we confirm the study [37] of the UV critical theory and provide the fixed point couplings and conformal data up to third order in ϵ (both for infinite and finite-N c scenarios) by considering four-loop gauge, three-loop Yukawa, and quartic β functions.We also determine three-loop anomalous dimensions for dimension-three operators and discuss peculiarities in their computations as compared to Ref. [37].
The paper is organized as follows.Section II provides the main information about the considered model and operators.In Sec.III we give the details of computation of the renormalization-group (RG) functions.In Sec.IV we demonstrate our results for fixed points, anomalous dimensions and scaling exponents for finite N c .The bounds on the conformal window are given in Sec.V.We conclude in Sec.VI.The full expressions for beta functions and various Representations of matter fields under gauge SU (Nc) and flavor UL(N f ) and UR(N f ) groups.
anomalous dimensions can be found in Appendixes A 1, A 2, and A 3. App.B contains additional information related to finite-N c results.

II. MODEL DESCRIPTION
We consider four-dimensional theory with SU (N c ) gauge fields coupled to N f massless Dirac fermions and a scalar singlet field H.The last is uncharged under the gauge group and carries two flavor indices, such that it can be written as a N f × N f complex matrix.The corresponding Lagrangian is where F a µν is the field strength of the gauge bosons G a µ with a = 1, . . ., N 2 c − 1.The trace in Eq. ( 1) runs over both color and flavor indices and ψ = ψ L + ψ R are fermions with P L/R = 1  2 (1 ± γ 5 ).In what follows we use a linear R ξ -gauge with L gf = − 1 2ξ (∂ µ G a µ ) 2 together with the corresponding ghost Lagrangian L gh .The theory (1) is invariant under global G = U L (N f ) × U R (N f ) "flavor" symmetry corresponding to independent unitary rotations of left-and right-handed chiral fermions.The matrix scalar field H is colorless but transform under G (see Table I).In this paper we also consider a class of operators that breaks the flavor symmetry down to diagonal U (N f ).One can introduce independent couplings for these operators resulting in the following additional contribution to the Lagrangian (1) The choice of the G-breaking terms is dictated by the fact that the operator O 1 = ψψ, also considered in Ref. [37], mixes under renormalization with the two operators coupled to h 2,3 .The model has four dimensionless couplings: gauge coupling g, the Yukawa y, and two quartic scalar couplings u and v.One usually introduces a set of rescaled couplings [38] The latter has been done since we consider the Veneziano limit [39] with N f , N c → ∞.One also introduces a parameter that becomes continuous and may take any value between (− 11 2 , ∞).The benefit of the Veneziano limit is that it allows systematic expansions in a small parameter.In our work we suppose that and treat it as a small control parameter for perturbativity.In order to study dimension-three operators in the Veneziano limit, we rescale the corresponding couplings and the operators as This allows one to absorb all corrections with positive powers of N c and N f appearing in the beta functions for ⃗ κ into the rescaled couplings given in Eqs. ( 3) and (6).It is worth noticing that the parameters m ϕ , h 1 and h 2 are rescaled in the same way as the dimensionless couplings y, u, and v, respectively.

III. CALCULATION METHODS
The beta functions and anomalous dimensions in the MS scheme can be computed by standard methods, so we omit full description of the techniques refereeing to appropriate literature.In this section, we only discuss peculiarities of the current calculation.Let us mention here that we do not rescale the couplings and operators according to Eqs. ( 3) and (6) in our explicit computation.The transition to the Veneziano-limit normalization is carried out at the final stage.Nevertheless, we present all our results in terms of α g,y,u,v and for the rescaled operators (7).
In order to compute the required RG functions, we rewrite the Lagrangian (1) in terms of real scalars ϕ a utilizing a decomposition (a = 1, . . ., with and satisfying the following identities Such a decomposition gives rise to the Feynman rules (see Fig. 1) for the vertices involving ϕ a .In Fig. 1 we also indicate the flavor "flow" for the double-trace and single-trace scalar couplings that can be associated with the left (L) and right (R) chiral fermions.In the absence of U L (N f ) × U R (N f ) breaking, the "left-" and "right-handed" flows are "conserved" separately.
We implemented the rules in the DIANA [41] package based on QGRAF [42] and utilize FORM [43] to deal with index contractions and to compute Feynman integrals via the MATAD [44] code.In order to derive RG equations for u, v and y, we generate Green functions corresponding to radiative corrections for the tree-level Yukawa and quartic vertices up to three loops and extract local divergent terms by applying suitable projectors.Our explicit calculations heavily rely on well-known infrared rearrangement (IRR) trick [45,46], which allows one to deal only with fully massive vacuum integrals.The fermion-fermion-scalar interaction involves the γ 5 -matrix, which requires special treatment in dimensional regularization (see, e.g., Ref. [47]).In this paper, we restrict ourselves to the seminaive approach [48,49] and by explicit computation we prove that potential ambiguities do not appear in the final result for the RG functions, thus, providing an independent cross-check of the results obtained in Ref. [37].As for the four-loop gauge-coupling beta function the γ 5 -ambiguity can be fixed by means of Weyl consistency conditions [50][51][52][53].Moreover, we do not carry out explicit computations here, but use the RGBeta code [54] extended to 432-order in Refs.[55,56].
Let us now switch to the discussion of a family of dimension-three operators that includes O 1 = ψψ.In dimensionally regularized theory2 (d = 4 − 2ε) within the MS-scheme, we have a relation between bare and renormalized quantities where all couplings in ⃗ κ R , and all operators in ⃗ O R have in our case the mass dimension one and three, respectively.The corresponding renormalization matrices Z κ and Z O are defined as The renormalization matrices involve poles in ε and depend on the renormalization scale µ together with the running dimensionless couplings from (1) denoted collectively by α(µ).They satisfy In what follows, we will routinely use Z and Z to denote renormalization constants with and without explicit dependence on the renormalization scale µ.The diagonal matrices involving powers of µ ε account for mass dimensions of bare couplings and operators.We include them in the definition of Zs for convenience, since we can write the beta functios of ⃗ κ R in a compact form and relate the matrix anomalous dimension γ κ to the anomalous dimensions of the dimension-three operators Both γ O (α) and γ κ (α) should be finite in the limit ε → 0, which serves as a welcome check of the computation.As a consequence, in d = 4 we have This relation can be used in two ways.Given the beta functions of (a closed set of) the dimension-1 couplings, one can extract the matrix anomalous dimension for the corresponding operators.In case we know γ O , it is possible to reconstruct the MS beta functions for the operator couplings, e.g., via FIG. 2. Feynman rules for dimension-three operators.All the operators break the flavor symmetry G by "flipping" the "chirality" of the flavor flow, which in the case of scalar operators we indicate by a box with a cross inside.
In this paper we explicitly carry out the renormalization of dimension-three operators and compute γ O up to three loops.In addition, we also cross-check our results at lower loops by adding (2) to the Litim-Sannino model (1) implemented in public computer codes ARGES [57] and RGBeta3 [54] and computing ⃗ β κ .To find three-loop corrections to Z O , we consider insertions of the operators O 1−3 into one-particle irreducible (1PI) Green's functions.The corresponding Feynman rules are given in Fig. 2. The renormalized 1PI Green functions are finite and are related to the bare ones via Here Õj is a product of external fields at different space-time points (either bare or renormalized) corresponding to the local operator O j with a property that ⟨O i • Õj ⟩ tree ∝ δ ij .The factors Z j,ext account for the external fields renormalization entering Õj , e.g., for ψ 0 = Z ψ • ψ, and H 0 = √ Z H • H.We determine the matrix elements of (Z O ) ij in the MS scheme order-by-order from the requirement that there are no ε poles in the rhs (21).
It is important to stress that if we ignore O 2,3 (or the h 2,3 couplings in the Lagrangian (2)), we can compute (Z O ) 11 and (γ O ) 11 (the latter is denoted by γ m ψ in Ref. [37]) by considering ⟨O 1 (x)ψ(y) ψ(z)⟩ 1PI up to the two-loop order without hitting any difficulties.At three loops there are diagrams (see,e.g., Fig. 3a), which require an insertion of the O 2 operator as a counterterm (Fig. 3b).Because of this, we are forced to include O 2 in the game.The latter mixes with O 1 already at the one-loop order4 .Moreover, starting from one loop, the O 3 operator is needed to account for all the divergences appearing in four-point functions ⟨O 2 ϕ a ϕ b ϕ c ⟩ 1PI .
We also compute the divergences of the two-point functions ⟨O i (x)ϕ a (y)⟩ 1PI for i = 1, 2, 3 and extract the mixing of O i with O 4 = 1/2(∂ 2 Tr(H) + h.c.), thus, modifying (13) as [ This gives rise to a 4 × 4 mixing renormalization matrix Z entering where i = 1, . . ., 4 and ZO is schematically depicted in Fig. 4. If we take into account the equations of motion (EOM) for the bare operators where Z y,u,v are renormalization constants for the corresponding couplings, and express O 4 as a linear combination of ⃗ O, we obtain another 3 × 3 matrix schematically presented in Fig. 5.The corresponding anomalous dimension is given by O .This can be interpreted as the fact that only two of the considered dimension-three eigenoperators are independent, while the remaining one is a descendent of Tr(H) + h.c.
Before discussing results, let us note how our approach differs from that of Ref. [37].The authors of Ref. [37] utilized the same IRR technique and the seminaive treatment of γ 5 in their three-loop calculations.However, instead of DIANA and MATAD, they used a private framework to deal with Feynman diagrams generated by QGRAF.In our study we do not rely on the dummy-field method (see, e.g., Ref. [58] and references therein) to derive the beta functions of dimensional couplings, but extract them from the Green's functions with insertions of the corresponding operators.Moreover, the authors of Ref. [37] considered only the O 1 = ψψ operator and did not account for possible mixing with O 2,3 .Because of this, the anomalous dimension5 denoted by γ m ψ is not an eigenvalue of the corresponding matrix and, if computed at a fixed point, does not represent a correction to canonical scaling.
In what follows, we present all our results in terms of rescaled couplings and for the rescaled operators O ′ (7).All expressions are available in computer-readable form as Supplemental Material.

IV. RESULTS
In this section, we summarize our results for β-functions for (in)finite N c and anomalous dimensions.Moreover, we determine fixed points and universal scaling dimensions up to the third nontrivial order in the Veneziano parameter also for both cases.
One-particle-irreducible Green functions used to extract matrix elements of ZO.The factor Green's functions used to extract anomalous dimensions γO of the dimensions-3 operators with the account of equations of motion, expressing O4 as a linear combination of O1−3.

Fixed points
As it is well known in perturbation theory the β-functions can be given as where denotes the n-th loop contribution, and x = {g, y, u, v}.In Ref. [37] they were found for the first time and provided both in the Veneziano limit and with finite-N c corrections.We recomputed them independently and provide the lengthy expressions in Appendix A 1 and as Supplemental Material.
Using the 433-order β functions and solving β i (α * j ) = 0 systematically, we determine interacting fixed points up to complete third order in the small parameter ϵ.We expressed fixed points as a series expansion α * x = c (1)  x f (1)  x (N c )ϵ + c (2)  x f (2)  x (N c )ϵ 2 + c (3)  x f (3)  x (N c )ϵ 3 + O(ϵ 4 ), (29) where x = {g, y, u, v}.We have 6 possible solutions with fixed points, however, following Refs.[5,6,14], we choose a fully interacting UV fixed point (α * g,y,u,v ̸ = 0) that exhibits asymptotically safe behavior.The coefficients c (i) x of series expansion can be found by solving the beta functions in this limit order by order for the stationary point.We recomputed the exact expressions for c (i) x and found agreement with Ref. [37] c (1)   For convenience, we provide the numerical results in the Veneziano limit [37] α * g = 0.456ϵ + 0.781ϵ 2 + 6.610ϵ 3 +24.137ϵ 4 , (39) where we also include subleading terms O(ϵ 4 ) terms that will be modified when the 544-result will be available [37].
In addition, Table II shows various β (l) x evaluated at the fixed point (α * g , α * y , α * u , α * v ) given above.One sees the typical size of the coefficients that cancel at each order of ϵ-expansion up to O(ϵ 5 ) in β g , and up to O(ϵ 4 ) in other beta functions.Subleading terms that do not add up to zero due to missing high-order corrections are also indicated.
The functions f (i) x first introduced in Ref. [14] capture the dependence on N c and can be computed in the same manner as the c (i) x .It should be noted that lim Nc→∞ f (i) x ≡ 1, which means that the FPs found for in(finite) cases are connected continuously.Therefore, for the infinite-N c case we have only c x , which are exact numbers.However, the full form of f x is complicated 6 .Because of this, following the ideas of [57], we fitted all the finite-N c corrections in the range N c ∈ [3, 100] as ratios of two second-order polynomial up to ϵ 3 .The results can be found in Appendix B. Then, we carried out another fit when ratios of two fourth-order polynomials were considered and obtained the following expressions: which turn out to be more accurate. .The maximal value of the mean squared error (MSE) for the f x with our fourth (second) order polynomial fits in each order of ϵ is as follows: 10 −12 (4•10 −11 ), 10 −9 (8•10 −5 ), 8•10 −5 (4•10 −3 ).While the fit with the second-order polynomials seems to provide good approximation we restrict ourselves to a more precise four-order result.
At the end, we carry our numerical comparison of the two types of fits and exact results.The green line in Fig. 6 corresponds to the second order fit (see Appendix B), and purple line is our fourth order fit (43), (44), (45), (46).One can see that the green line "misses" exact points in Fig. 6 in the presented range of N c , while the purple line goes through the dots.

Critical exponents
Let us study the universal critical exponents.The latter can be obtained as the eigenvalues of the stability matrix which again we expand as a power series and find perfect agreement with [37] (note that for θ 1 expansion starts at ϵ 2 so c (1)  43),( 44),( 45),( 46), the green line is 2nd order fit (B1).The black dots correspond to exact values.It should be noted that we omit the point Nc = 3, since it lies higher in the FP's scale, therefore the difference between the fits is not clearly visible.

c
(1) Numerical evaluation of these coefficients gives in the Veneziano limit [37] and for not too large ϵ (see below) there is only one relevant direction corresponding to θ 1 < 0 giving rise to an asymptotically-safe scenario.In Eq. ( 58) we again highlight terms that are not determined precisely in the 433 approximation and will be modified by high-order terms.
Finite-N c corrections are incorporated in the factors, which we approximate as7 : The maximal value of MSE for the given (second-order result from the appendix) approximations are 10 −19 (10 −10 ), 10 −10 (10 −3 ), and 10 −7 (0.15) for first, second and third nontrivial orders in ϵ.

Other anomalous dimensions
In this subsection we provide results for the scalar and fermion anomalous dimensions.In addition, we consider the operator Tr(H † H) coupled to m 2 in (1) together with the anomalous dimensions corresponding to dimension-three eigenoperators discussed earlier.The full expressions for anomalous dimensions beyond the Veneziano limit can be found in Appendix A 3 and are also available in the Supplemental Material.
At the fixed point the anomalous dimensions are given as a series expansion in ϵ.In the Veneziano limit we have gauge-independent coefficients c (1) for the scalar fields, and c which is valid in the Landau gauge ξ = ξ * = 0 (see also Ref. [37]).The gauge-independent result for m 2 can be cast into c Substituting the coefficients into the power series, we have in the Veneziano limit and (for ξ = 0) [37] γ H = 0.2105ϵ + 0.4625ϵ 2 + 2.4711ϵ 3 , (67) The factors that account for the finite-N c corrections are again approximated as where the maximal MSE is no worse than 3 • 10 −7 in comparison to 10 −1 for the approximations given in Eq. (B2).Finally, the eigenvalues for the dimension-three operator 4 × 4 anomalous dimension matrix are given by with Evaluating the coefficients in the Veneziano limit, we obtain One can see that all but γ 2 is positive for ϵ > 0 and does not pose a problem to unitarity.Moreover, γ 2 = −γ H corresponds to the linear combination of operators 9 O 1−4 that vanish due to the equations of motion, while γ 1 = γ H corresponds to a linear combination of O 1−3 that becomes a descendant of Tr(H) + h.c.when EOMs are imposed.
As for the finite-N c factors, we provide the following approximate expressions 9 This fact can be deduced by considering the renormalization of local operator The maximal value of MSE in these fits are 10 −9 , while the corresponding value for the ratio of two second-order polynomials (B4) reaches 4 • 10 −2 .

V. CONFORMAL WINDOW
In this section we investigate the size of the UV conformal window for the asymptotically safe theory with action equation ( 1) using perturbation theory.

A. Constraints on the UV conformal window
We can find the UV conformal window directly from the expressions for fixed points and scaling exponents given in the previous sections.In the case when we retain only first three non-vanishing powers of ϵ, we will call the bound on ϵ strict.The reason for this is that the higher-order ("subleading") coefficients in the power expansion ( 29), (48) are not (yet) accurately determined due to the absence of higher loop terms in β-functions.This scheme is dictated first by the following constraints • no fixed point merger [14] (the collision of the UV fixed point with an IR fixed point corresponds to θ = 0).
The second strategy employs the approximation, where we retain subleading terms in ϵ, so we refer to its bounds as subleading.There we also take into account all the above mentioned constraints from couplings, vacuum stability and critical exponents.

B. Investigation of the UV conformal window
The UV conformal window can be investigated using the above mentioned constraints.To do so, we first can equate the perturbative expressions for FP couplings to unity, (α * u + α * v ) and scaling exponents to zero and choose the smallest positive solution for ϵ.Second, we can find the Padé approximants for these constraints and make the same manipulations.We represent our results as α * = ϵP ij and θ = ϵ 2 P ij , where P ij are Padé approximants and i + j = 2(3) for the strict (subleading) case.However, we cannot confidently trust the obtained results, because these approximants contain non-physical poles.Nevertheless, they give tighter constraints on the conformal window and we provide all approximants that can be constructed from available series together with the corresponding bounds.
• From perturbative expansion of couplings ( 39)-( 42), we note that the tightest bound on ϵ strict arises from the gauge coupling.First, equating the gauge coupling to unity, we can find ϵ strict and ϵ subl where the first line correspond to the strict case, the second to subleading, and find the bounds • In the same manner we can investigate the bounds arising from the vacuum stability condition [60]: To provide more stronger constrains for ϵ, we use the Padé approximants: The UV conformal window in this case is constrained as • After the calculation of scaling exponents (58), we notice the behavior of θ 2,3,4 with same-sign corrections at every order.However, the sign of leading ϵ 2 coefficient for the relevant scaling exponent θ 1 differs from other loop terms, which is the indication of possible FP merger.Thus, we can extract the constraints from the relevant scaling exponent, solving θ 1 = 0: And using the Padé approximation: we get It should be noted, that at 433 order we have only one fixed point merger.However, if the subleading tendency continues in the 544 order, it will lead to the additional θ 3 merger ( 57), which we do not study in this paper.Moreover, we illustrate our results in the Fig. 7.Here we show up the bounds for strict and subleading approaches (full lines) (88),( 91),(94) and their Padé approximants (dashed lines) that provide tightest constraints(90),( 93),(96) at 433 order.In addition, we add 211 and 322 orders for comparison.From these figures and from (88),( 90),( 91),( 93),(94), and (96), we can deduce that the ϵ subl bound is systematically tighter than the ϵ strict bound, which arise from the last summands (39), (55).The same situation was obtained at 322 order, see Ref. [13].
At the end, we demonstrate the size of the UV conformal window in Fig. 8.The left panel includes ϵ strict bounds arising from vacuum stability (91),(93).For comparison we have also indicated the previous bound at 322 order [13].The right panel illustrates the boundaries for finite values of N c .From this picture we can see that all constraints (88),(91),(94) share roughly the same rate of convergence.In addition, it is clear, that finite-N c corrections contract the conformal window (full lines) in comparison to infinite results (dashed lines).
The dashed lines represent the asymptotic value; the full lines represent the upper boundary for the ϵ as functions of Nc.

C. (De)stabilizing fluctuations
Following the ideas of Refs.[13,14], let us consider into which direction the higher loop corrections shift the beta functions.To do this we substitute the fixed points to order O(ϵ) 3 back to the beta functions and take series expansion in ϵ up to the first non-vanishing order.Keeping only highest available terms in RG functions, we obtain (c.f., also Table II) Negative shifts to the beta functions (∆β < 0) are supposed to stabilize the UV fixed point [13,14], while ∆β > 0, conversely, shift the zero towards larger values, which could even destabilize it.Thus, the obtained leading shifts for finite N c gives us a qualitative picture of the trend from higher-order loop contributions.It should be noted that these shifts have changed their signs compared to the previous results obtained in Refs.[13,14].
In the same manner, we insert the fixed points up to ϵ 3 to the full available beta functions and find the subleading shifts for finite N c (see Table II  (99) We see that the results for β g,y,u change signs as compared to (98).However, they retained their behavior as in the case of Refs.[13,14].Therefore, we can expect, that in the 544 order we will obtain similar results.

VI. CONCLUSION
The availability of interacting UV fixed points in particle physics presents numerous prospects for constructing models, see, e.g., Ref. [20].However, comprehending the size of the corresponding conformal window is equally crucial for any real-world applications.In this paper we investigated the Litim-Sannino model with action (1) at the 433 order.Extending the findings of [5,13,14] and confirming the results of [37], we have performed a full search for interacting fixed points and computed scaling exponents and anomalous dimensions up to third order in the small parameter ϵ (4).It was also noted that their full expressions for finite N c are complicated, therefore, we performed fourth(second) order fits and provide the corresponding approximate expressions together.By comparing means square errors, we conclude that the ratio of fourth-order polynomials provide better approximation in certain cases; see, e.g., Fig. 6.
Moreover, we studied the size of the conformal window imposing conditions on the fixed point values of the couplings and scaling exponents.We used different approximation orders and estimate the effect of subleading corrections.We compared various restrictions coming from the perturbativity of the strong coupling (88), vacuum instability (91) and possible fixed point merger (94).Despite their qualitatively different origins, constraints are quantitatively similar, with vacuum stability offering the tightest bound (91).In addition, following Ref.[37] we tried to resum the ϵ series by means of Padé approximation and find that it gives even stronger constraints (90),( 93),(96).However, the presence of non-physical poles in the approximants used to derive the bounds undermines our confidence in their accuracy.Perhaps, this should be explored using other types of approximations and when high-order contributions in the 544-scheme will be available.We summarize our results in Fig. 8, where we illustrated the asymptotic safety regime together with the size of the UV conformal window.At the end, the asymptotically safe quantum field theories, which lie within the allowed conformal window were also demonstrated in Eq.(97) and Fig. 8.
Furthermore, we notice that the authors of Ref. [37] mentioned that the anomalous dimension γ m ψ of the fermion mass, which they computed in their paper, is negative to the leading order in ϵ, while all the next-to-leading order terms are positive.Because of this, γ m ψ can become negative and potentially pose a threat to unitarity [37].We argue that this negative leading-order result is nothing else but the leading-order contribution to (−γ H ). Positive higherorder terms computed in Ref. [37] cannot be trusted when evaluating scaling dimensions, since at the two-loop level the mixing comes into play.We account for this mixing in the present study and compute the anomalous-dimension matrix eigenvalues, one of which should be (−γ H ) at any loop.But this is not the end of the story.We also argue that if EOMs are taken into account the anomalous dimension matrix of dimension-three operators is modified γ O → γO such that γO has the same eigenvalues as γ O but with the flipped sign of the "dangerous" eigenvalue (−γ H ) → γ H .This can be anticipated and represents correct scaling of the operator that enters EOMs alongside with O 4 ∝ ∂ 2 Tr(H).Our analysis shows that all eigenvalues of the reduced anomalous dimension matrix are positive so dimension-three operators that break G flavor symmetry do not spoil unitarity.
We believe that our findings can be used in a more elaborate analysis of the asymptotic safety in the Veneziano limit and beyond the latter.It is interesting how the results will be modified when the 544 order beta functions will be available, e.g., in connection with possible additional FP merger due to potentially negative contribution of O(ϵ 4 ) to θ 3 .

FIG. 3 .
FIG.3.A three-loop diagram (a) with a O1 = ψψ insertion that requires a two-loop counterterm due to O2 (b).The grey blob gives rise to a one-loop contribution to (ZO)21.
and is different from γ O(16).At a fixed point, γ * O (γ * O ) has a eigenvalue γ * H (−γ * H ). The other two eigenvalues of γ * O coincide with those of γ *

FIG. 6 .
FIG.6.The dependence of couplings from the number of colors Nc for fixed ϵ = 0.09.Here the purple line is our fourth-order polynomial fit (43),(44),(45),(46), the green line is 2nd order fit (B1).The black dots correspond to exact values.It should be noted that we omit the point Nc = 3, since it lies higher in the FP's scale, therefore the difference between the fits is not clearly visible.

FIG. 8 .
FIG.8.Left panel: conformal window with asymptotic safety (yellow band), also showing regimes with asymptotic freedom (green) and effective theories (grey).Dots indicate the first few integer solutions (97) at 433 order.Moreover, this panel contains the new upper bound on ϵ at 433 order (91), and the tighter Padé approximant bound (dashed line) (93).Previous upper bounds at 322[13] order also indicated.The right panel compares different schemes (88),(91),(94) as given by ϵstrict.The dashed lines represent the asymptotic value; the full lines represent the upper boundary for the ϵ as functions of Nc.

TABLE II .
The model beta-functions at the fixed point as series in ϵ.Cancellations between contributions from different loop orders are presented up to O(ϵ 5 ) for βg and up to O(ϵ 4 ) for βy,u,v.Higher-order terms marked by the red color are not reliably calculated and are not summed up to zero, but (the sum) can be used to estimate the size of typical high-order contribution.