Analytic Coupling Structure of Large $N_f$ (Super) QED and QCD

We study the analytic properties of the 't Hooft coupling expansion of the beta-function at the leading nontrivial large-$N_f$ order for QED, QCD, Super QED and Super QCD. For each theory, the 't Hooft coupling expansion is convergent. We discover that an analysis of the expansion coefficients to roughly 30 orders is required to establish the radius of convergence accurately, and to characterize the (logarithmic) nature of the first singularity. We study summations of the beta-function expansion at order $1/N_f$, and identify the physical origin of the singularities in terms of iterated bubble diagrams. We find a common analytic structure across these theories, with important technical differences between supersymmetric and non-supersymmetric theories. We also discuss the expected structure at higher orders in the $1/N_f$ expansion, which will be in the future accessible with the methods presented in this work. Understanding the structure of the large-$N_f$ expansion is an essential step towards determining the ultraviolet fate of asymptotically non-free gauge theories.


I. INTRODUCTION
The discovery of four-dimensional asymptotically safe quantum field theories [1] has spurred recent phenomenological and theoretical interest. The original result made use of the Veneziano limit, in which one considers a large number of both colors and flavors. These theories feature perturbative safety and contain not only gauge and fermion degrees of freedom but also scalars. It is therefore theoretically and phenomenologically important to investigate the ultraviolet fate of nonasymptotically free gauge theories featuring a small number of colors but still a large number of flavors. In particular, one wishes to either exclude or demonstrate that a large number of flavors can lead to an asymptotically safe scenario in gauge-fermion theories.
This quest has revitalized the study of quantum field theories at a large number of flavors N f . The timeliness of our investigation is further corroborated by the fact that the large-N f nonasymptotically free regime of gauge-fermion theories is being, for the first time, investigated via firstprinciple lattice simulations, where we expect the first results to appear soon [2]. An intriguing property of this limit is that, at each order in the 1=N f expansion, only a finite number of underlying topologies contributes, where each gauge line is dressed with matter loops. Correspondingly, at fixed order in 1=N f , the number of diagrams grows polynomially with the loop order, suggesting that a closed-form resummed result with a finite radius of convergence may be achievable.
The first notable study is large-N f quantum electrodynamics (QED) [3], while large-N f quantum chromodynamics (QCD) is considered later in Ref. [4]. A historical summary of the techniques and earlier results can be found in Refs. [5,6]. The generalization to a wide class of semisimple gauge-Yukawa theories appeared only recently in Refs. [7][8][9][10][11][12]. For gauge theories with different fermion matter representations, the new phase diagram as a function of the number of flavors and colors was put forward in Refs. [13,14], and it was termed Conformal Window 2.0, extending and generalizing the original phase diagram of Refs. [15,16] to contain, besides an infrared conformal window, also an ultraviolet (safe) one.
In the limit of a large number of matter fields, it is natural to introduce the 't Hooft coupling with the gauge coupling g and the Dynkin index S 2 ðRÞ, normalized to 1=2 for the fundamental representation. The generic beta function has a formal expansion as an inverse series in N f : where each β ðkÞ ðKÞ has itself a perturbative expansion in the 't Hooft coupling K. Similar expansions hold for anomalous dimensions and other critical quantities. Asymptotic freedom is lost for theories at a finite number of colors and a large number of flavors, and therefore such theories can only be fundamental if they develop an interacting fixed point in the ultraviolet. This cannot occur in perturbation theory without Yukawa interactions [1,17], but it may occur nonperturbatively above a critical number of flavors [9]. To see how this might work, let us schematically consider the leading nontrivial order-1=N f beta function, which up to a normalization reads At this order, the function β ð1Þ ðKÞ=K 2 must develop a singular behavior for the beta function to develop a zero as N f → ∞. This indeed happens for QED and QCD, as summarized in Refs. [5,9,14].
In this paper, we investigate this phenomenon further and make a systematic study of the analytic structure of the 't Hooft coupling expansion of the leading large-N f beta function for QED, QCD, super QED (SQED), and super QCD (SQCD). We discover that for each theory, the 't Hooft coupling expansion is convergent, but a large number of expansion coefficients are needed in order to determine accurately the radius of convergence and to extract the logarithmic nature of the first singularity of the theory. Additionally, by a detailed investigation of the summation properties of the beta function at a leading order of 1=N f , we identify the physical origin of the singularities from the iterated self-energy diagrams. We find a universal analytic structure across the theories investigated here, while being able to resolve important physical differences between supersymmetric and nonsupersymmetric theories.
The paper is organized as follows: In Sec. II, we investigate large-N f QED and introduce the relevant mathematical tests and tools that we use for the various theories. These include the asymptotic analysis of the expansion coefficients and Padé approximants. We then identify the physical origin of the poles. We extend this analysis to QCD, SQED, and SQCD in Sec. III. There we also elucidate and highlight the crucial differences among the various theories. We present our conclusions in Sec. IV. In the Appendix A, we briefly review Darboux's theorem, relevant for the large-order behavior of the expansion coefficients, and in Appendix B, we describe on a technical level how we extracted the numerical coefficients of the beta function.

II. LARGE-N F QED
QED is structurally the simplest gauge theory, but it still has a rich perturbative and nonperturbative structure, which we probe here in the large-N f limit. The QED beta function has been computed in Ref. [3] at the leading nontrivial order in the 1=N f expansion: Here, the integrand function for QED is This resummed beta function is shown in Fig. 1 compared with the supersymmetric version of the model. The first few terms in the 1=N f term of Eq. (4) read In the coefficients c ð1Þ n , we recognize characteristic π powers and zeta values, familiar from algebraic properties of Feynman perturbation theory and harmonic polylogarithms [18,19].
Before discussing the analytic structure of the integral representation in Eq. (4), consider the following pragmatic question: supposing, as is often the case, one were given only a finite number of terms of the expansion in Eq. (6), what could we learn about the physical nature of the expansion? There is a well-developed formalism to address such a question [20,21]. The first observation is that the expansion is convergent. This can be confirmed by a variety of simple ratio tests: for example, the radius of convergence c ð1Þ Ã can be deduced from the limit c ð1Þ Ã ¼ lim n→∞ jc ð1Þ n j −1=n , or from the limit 1=c ð1Þ Ã ¼ lim n→∞ jc ð1Þ n =c ð1Þ n−1 j. However, more information about the physics of the expansion can be obtained by applying Darboux's theorem [20][21][22], which relates the rate of growth of the perturbative expansion coefficients to the behavior of the expansion about the leading and subleading singularities. Note that this is a stronger statement than simply saying that the location of the nearest singularity determines the radius of convergence. The expansion coefficients also encode further information about the nature of the singularity. The general argument is summarized in Appendix A.

A. Asymptotic analysis of expansion coefficients
Our goal in this section is to deduce physical information from a finite number of expansion coefficients c ð1Þ n in Eq. (6). We studied these coefficients up to order M ¼ 60, and from 60 terms we obtain a great deal of asymptotic information. Using Richardson extrapolation [23] with these 60 coefficients, we learn that as n → ∞, where R 0 ¼ 0.063044292, R 1 ¼ −0.013027009, and R 2 ¼ 0.0033170626. These numbers can be fit to R 0 ¼ 28 45π 2 , R 1 ¼ − 9 70π 2 , and R 2 ¼ 11 336π 2 , identifications that can be confirmed to higher precision using higher-order Richardson extrapolations. We explain the origin of these coefficients below in Eq. (12).
Thus, using Darboux's theorem (see Appendix A), from these 60 perturbative expansion coefficients we learn that (i) the radius of convergence is 15=2; (ii) the leading singularity of β there are no higher-order corrections associated with this singularity; and (iv) there are higher-order corrections associated with further singularities at K ¼ 21 2 and K ¼ 27 2 . Interestingly, we need approximately M ¼ 30 terms of the expansion to be able to deduce precise information about the leading singularity. With fewer than M ¼ 30 terms, even identifying the radius of convergence to be 15=2 is noisy; see Thus, the leading behavior of the 1=N f correction to the beta function as K approaches the radius of convergence is This implies that in order to obtain a zero of the beta function in the large-N f limit, we must approach a nonperturbative fixed point at [1] This physical information has been deduced from a finite number of terms in the perturbative expansion of β ð1Þ ðKÞ. However, since we have an all-orders integral representation [3] in Eq. (4), we can probe the analytic structure more precisely by studying the properties of the integrand function F QED ðxÞ defined in Eq. (5). The singularities of the integrand are simple poles at x n ¼ 15 2 þ 3n, for n ≥ 0, generated by Γð 5 2 − x 3 Þ. These are the only singularities, as can be seen from the decomposition where the incomplete gamma function Γð 5 2 − x 3 ; 1Þ is regular. The potential poles at x ¼ 3 − 3n, with n ≥ 0, coming from the denominator in Eq. (5) are in fact canceled by the sinð πx 3 Þ factor in the numerator. Alternatively, one can rewrite the integrand using the gamma-function reflection formula as  7). From n ≈ 30 onwards, the coefficients agree with the expectation from the large-order behavior.
from which we see that the only singularities come from the secð πx 3 Þ factor, with the poles at x ¼ 3 2 and x ¼ 9 2 excluded. Therefore, the positions and residues of the (simple) poles of the integrand are These coincide precisely with the numerical values extracted from the asymptotic analysis in Eq. (7). Furthermore, the noisiness of the expansion coefficients at low order can be traced to the oscillatory nature of the sinð πx 3 Þ tanð πx 3 Þ factor in Eq. (11). Recall that the poles in Eq. (12) are simple poles of the integrand of the beta function in Eq. (4). After integration over x, these poles translate into logarithmic branch points of the beta function, which were found above [see, e.g., Eq. (8)] by a numerical Darboux analysis of the coefficients of the perturbative expansion of the beta function to finite order.

B. Padé approximations
Padé approximation is a commonly used method for studying perturbative expansions in physical systems [23,24]. Given the integrand F QED ðxÞ in Eq. (5) expressed in terms of gamma functions, there is a unique analytic continuation beyond its radius of convergence. However, if we only had a finite number of terms of the expansion, not its full analytic form, we could still probe beyond the radius of convergence using Padé approximation.
Padé approximants construct analytic continuations of truncated Taylor series (i.e., polynomial) approximations to functions, expressing the given polynomial as a ratio of two polynomials of lower order, with coefficients determined purely algorithmically. Padé approximants thus convert a polynomial to a rational function, which can also be expressed as a partial fraction expansion, whose residues and poles are determined by the coefficients of the original truncated Taylor series. This means that Padé approximants tend to be quite good at representing functions with poles, while they are less good at representing functions with branch cuts [23,24].
The conversion of a truncated Taylor series to a Padé approximant, is algorithmic, leading to a ratio of two polynomials P R ðxÞ and Q S ðxÞ, of orders R and S, respectively, where R þ S ¼ M. It is, in fact, a built-in function in symbolic mathematics languages such as Maple or Mathematica.
We took up to 60 terms of the expansion about x ¼ 0 of the integrand F QED ðxÞ and converted it to a diagonal Padé approximant P ½M=2;M=2 ðxÞ, for various values of M. In Fig. 3, we display the function F QED ðxÞ together with the diagonal Padé approximants starting from M ¼ 20, 30, 40, 50, 60 coefficients. With M ¼ 20 coefficients, we do not even "see" the first pole. With M ¼ 30, coefficients we accurately probe the first pole, but not the second pole. For the second pole, we need approximately M ¼ 40 coefficients, while with M ¼ 50 coefficients we accurately resolve the third pole. These numbers are consistent with the number of coefficients required in the ratio test and asymptotic analysis of the beta-function coefficients in Sec. II A, to resolve the logarithmic singularities of β ð1Þ ðKÞ.
In fact, a full Padé analysis constructs the "Padé table" of all Padé approximants P ½R;S ðxÞ, with R þ S ¼ M. It turns out that certain off-diagonal approximants are even better at representing the integrand function F QED ðxÞ. This can be understood from the analytic representation of the integrand in Eq. (11). Given that trigonometric and gamma functions have well-known product formula representations, we see that the Γð1 þ x 3 Þ=Γð 1 2 þ x 3 Þ factor in Eq. (11) is naturally represented as a near-diagonal Padé approximant, but because of the sin 2 ð πx 3 Þ= cosð πx 3 Þ factor, there is effectively one extra trigonometric factor in the numerator. Thus, a Padé representation whose numerator is a higher-order polynomial than the denominator polynomial will represent the analytic structure of the actual function F QED ðxÞ more accurately. We have confirmed that this is the case, starting from the 60 expansion coefficients, but we note that the simple diagonal Padé representations shown in Fig. 3 are already remarkably precise.

C. Physical origin of the poles
We have seen that the finite radius of convergence, K Ã ¼ 15 2 , of the expansion of the 1=N f beta function β ð1Þ QED ðKÞ can be traced directly to the leading pole of the Γð 5 2 − x 3 Þ factor in the integrand function F QED ðxÞ. This gamma factor arises because it enters the leading 1=N f computation via iteration of the basic building block of the one-bubble self-energy diagram, whose amplitude Π 0 is given by In the resummation, this amplitude typically enters the full beta function as its inverse, 1=Π 0 , and its argument is rescaled with the value of the 1=ϵpole 1 [6]. For the QED computation, the value of the 1=ϵ-pole in (14) is 2 3 . Consequently, we expect the resummed 1=N f beta function to contain the factor Indeed, this factor appears in the integrand function F QED ðxÞ and governs the pole structure underlying the asymptotics of the perturbative expansion coefficients, and the structure of the Padé approximations to the integrand function F QED ðxÞ. Knowing this, one can devise improved expansions in which this Π −1 0 ð 2 3 xÞ factor is factored out, with only the remaining factors needing to be analyzed. Not surprisingly, this leads to noticeable improvement of the resulting Padé approximations, and a much faster approach to the asymptotic behavior of the expansion coefficients.
The general idea of using Padé approximants to study the behavior of beta functions is significant for analyzing higher orders in the 1=N f expansion, for which no closed-form resummation formula is currently known. In particular, these methods may allow us to access the leading-order pole structure at higher orders in the 1=N f expansion, if enough coefficients can be extracted from the relevant diagrams. This information has quantitative implications for the stability, size, and structure of the asymptotically safe conformal window. One can imagine the following possible scenarios at higher orders in the 1=N f expansion: (i) A new singular structure may emerge closer to the origin, de facto disconnecting the putative fixed point in Eq. (9) from the Gaussian fixed point at the origin. The detailed structure of the new singularity would determine whether or not the theory remains UV safe to this order. Alternatively, if the radius of convergence of the series keeps shrinking as the order in 1=N f increases, the UV fixed point could eventually disappear. (ii) The current singular structure, and its location, could be further reinforced by higher-order corrections. This possibility is partially supported by the fact that the fermion self-energy amplitude is responsible for the singular structure of the theory. The order of the pole in the integrand might become stronger because n bubble chains appear in diagrams at the order Oð1=N n f Þ. In this case, the ultimate UV fate of the theory will depend on the character, sign, and strength of the reinforced singular structure. (iii) No further singularities emerge, or a new singular structure appears further away from the leadingorder one. This would be an indication that the putative fixed point in Eq. (9) is indeed physical. For example, the leading isolated pole, like the one that we will see appearing in QCD in the next section, is a candidate for this scenario, since it is not due to the fermion self-energy amplitude. Of course, even the ultimate confirmation of a nonperturbative zero in a generic beta function away from the origin is typically insufficient to establish the existence of a physical conformal field theory. Other critical quantities such as the variation of the a function or anomalous dimensions can potentially violate physical bounds [13,25].

III. COMPARING DIFFERENT THEORIES AND THEIR PHYSICS
The beta function at order 1=N f for an SUðN c Þ gauge theory was first calculated in Ref. [4] and written in a closed integral form in Ref. [5]. The result is The integrand function is now 1 This is the case for diagrams containing one resummed gauge chain.

ANALYTIC COUPLING STRUCTURE OF LARGE N f (SUPER) …
PHYS. REV. D 100, 015013 (2019) 015013-5 where dðGÞ, dðRÞ are the dimensions of the group G and the representation R, and similarly for the quadratic Casimirs C 2 . This result is very similar to that for QED and agrees with it in the Abelian limit C 2 ðGÞ → 0 and dðGÞ=dðRÞ → 1. The gamma factor Γð 5 2 − x 3 Þ produces the same pole pattern as in QED. However, a new isolated simple pole appears at x Ã ¼ 3, leading to a smaller radius of convergence, K Ã ¼ 3. This effect is purely due to the non-Abelian nature of the theory, as can be seen also from its residue: By an argument similar to that in Sec. II C, we can identify this pole with the gluon bubble loop rather than the fermion bubble loop. Since this diagram does not appear iterated in a chain, it does not result in an entire series of poles. The simple pole of F QCD ðxÞ at x Ã ¼ 3 leads to a logarithmic behavior of β ð1Þ QCD ðKÞ: This implies that in order to obtain a zero of the beta function in the large-N f limit, we must approach a nonperturbative fixed point at [1,13] Since the leading singularity for QCD is closer to the origin than for QED, fewer perturbative orders are required to resolve it using an asymptotic or Padé analysis. For example, for N c ¼ 3 and with fermions in the fundamental representation, the leading residue can be extracted with Oð10 −3 Þ accuracy from the asymptotic expansion of the coefficients already at ∼14th order. Retaining up to 15th order in the expansion of the integrand, the Padé approximant P ½7;7 ðxÞ gives a good reconstruction of the integrand within the radius of convergence. A similar analysis can be carried out for the other poles in the integrand: the results for the theories considered here are summarized in Fig. 4. Note that our result for QED is compatible with Ref. [26], where an expansion up to the 24th order was not sufficient to find a stable zero in the beta function.
Our results across the various theories indicate that the main factor determining the number of coefficients needed to resolve a given pole is the distance of the latter to the origin. 2 As discussed in the end of Sec. II C, the behavior near the leading pole is associated with the amplitude factor Π −1 0 , so we expect a similar relation between the number of coefficients and the location of the poles at the next orders in the 1=N f expansion. Furthermore, since no closed-form resummed perturbative expressions are known at higher orders in 1=N f , this motivates the importance of a similar Padé analysis of the perturbative series at 1=N 2 f .

B. Large-N f supersymmetric results
In this section, we review and discuss the results obtained in Refs. [27,28] for large-N f N ¼ 1 supersymmetric QED and QCD. The results are obtained in dimensional reduction (DRED) in d ¼ 4 − 2ϵ. For SQED, one finds while for SQCD, Notice again that the integrand expression for QCD in Eq. (17) agrees with the one for QED in Eq. (5) at the Abelian limit. For each of the SUSY beta functions in Eqs. (21) and (22), the integrand function FðxÞ has its first singularity as a simple pole at x ¼ 3 with a negative residue 3 ; see Table I. For example, this opposite sign explains the different behavior of the beta function near the first singularity, as shown in Fig. 1 for QED and SQED. Due to this negative sign for SQED and SQCD, the associated logarithmic singularity of the beta function cannot provide a cancellation between the first two orders in the large-N f expansion [Eq. (3)], and therefore no nonperturbative fixed point arises.
It is interesting to note that this conclusion holds also in the Novikov-Shifman-Vainshtein-Zakharov (NSVZ) scheme, which can be related to DRED by an order-by-order coupling redefinition [29]; see also Ref. [30] for details on such a relation. The well-known NSVZ beta function (see Ref. [31] for a recent discussion) is It admits a zero where the anomalous dimension takes the value We dropped the NSVZ label, as this quantity is scheme independent at the alleged fixed point. In the limit N f ≫ N c , the theory has lost asymptotic freedom, and therefore such a fixed point has to develop in the UV. However, due to the violation of the a theorem [25,32,33], it is disconnected from the IR Gaussian fixed point. The absence of an UV fixed point smoothly connected to the origin agrees with the large-N f result in the DRED scheme, where, in fact, no UV fixed point is seen to complete the Gaussian.

IV. CONCLUSIONS
Our analysis of the convergence properties of the leading 1=N f behavior of QED, QCD and their supersymmetric cousins has revealed several interesting features. We observe the emergence of a common analytic structure stemming from the leading 1=N f corrections, with the important difference that the coefficient of the logarithmic branch singularity is positive for QED and QCD, but it switches sign for their supersymmetric counterparts. The sign plays a crucial role when considering the UV fate of these theories. For example, for the supersymmetric case it implies that the theories are not fundamental, in agreement with other nonperturbative analyses.
We have demonstrated, by direct comparison with the full result, that the analysis of the large-order behavior of the 't Hooft coupling expansion is able to identify the location and nature of the leading logarithmic singularities, including the overall sign and magnitude of their coefficients. This suggests that a large-order analysis can be used in the near future to tackle the next-to-leading order in the 1=N f expansion, in the absence of a closed-form result. These corrections are crucial to test the singular structure of the leading 1=N f result, with important consequences for the UV fate of these theories.
This argument can be run in reverse, so that an analysis of the large-order behavior of the coefficients b n enables us first to determine the radius of convergence z 0 , and then also the nature p of the singularity-for example, a pole, or a branch cut, and what type of branch cut. The overall factor determines ϕðz 0 Þ, and further subleading information determines higher orders of the expansion of ϕðzÞ about z 0 . If the singularity is logarithmic, where ϕðzÞ and ψðzÞ are analytic near z 0 , then the Taylor expansion coefficients of fðzÞ near the origin have largeorder growth Once again, the large-order behavior of the convergent expansion coefficients determines the nature of the singularity and the fluctuations about it.

APPENDIX B: OBTAINING THE NUMERICAL COEFFICIENTS OF THE BETA FUNCTIONS
The Oð1=N f Þ beta functions discussed in this work are known in their resummed form, so we can reexpand them to 60 orders and perform a Padé analysis. The motivation for this is to obtain an estimate of how many perturbative terms are required in order to identify both the location and nature of the leading singularity, with a view towards a direct perturbative computation of the Oð1=N 2 f Þ beta functions, for which no resummed version is currently known. It is not a priori clear whether one might need 10 terms, or several hundred. Our work has shown that at Oð1=N f Þ roughly 30-40 perturbative terms are required, due to the lower-order oscillatory behavior, which we have associated with the appearance of the amplitude factor in Eq. (14). Since these factors appear also in an Oð1=N 2 f Þ computation, we expect that at least the same number of perturbative coefficients would be necessary in such a computation at a given higher order in the 1=N f expansion. We now comment briefly on the steps required to make such a computation feasible at higher orders in the 1=N f expansion. To this end, we first describe how the Oð1=N f Þ perturbative expansion is obtained in the diagrammatic approach.
The diagrams that contribute to the beta function at the order 1=N f in QED and QCD are displayed, e.g., in Ref. [9]. For the contraction of the diagrams, one can use the Mathematica package FeynCalc [34], which performs the trace over Lorentz and Dirac indices in d dimensions. Complicated diagrams can be traced with the symbolic manipulation system FORM [35,36], as well as the Mathematica package FormTracer [37]. The contracted diagrams can be evaluated with standard multiloop techniques along the lines of Refs. [38][39][40][41]. The diagrams contain fully dressed gauge propagators, and thus one has to apply reduction formulas that only hit nondressed propagators.
This procedure can be extended to higher orders in the 1=N f expansion to determine the perturbative coefficients up to arbitrary order in the coupling but subleading in 1=N f . At higher order in 1=N f , the loop order and the number of dressed propagators increases. Thus, at higher order in 1=N f , we expect that more efficient reduction formulas will be required.
Another complication at higher order in 1=N f is that at any order, the integrated diagrams contain gamma and hypergeometric functions, which need to be expanded in ϵ ¼ d − 4. For the nth-loop coefficient, we need to expand these functions up to ϵ n−1 . These expansions are slow, and at higher order in 1=N f , there will be more such factors to expand. Analytic expansions are only known for specific cases or only to low order. Thus, in our Oð1=N f Þ computation, we used numerical expansion methods-in particular, the package NumExp [42]. This numerical precision will need to be balanced with the precision required for the subsequent Padé analysis.