Towards the QED beta function and renormalons at $1/N_f^2$ and $1/N_f^3$

We determine the $1/N_f^2$ and $1/N_f^3$ contributions to the QED beta function stemming from the closed set of nested diagrams. At order $1/N_f^2$ we discover a new logarithmic branch-cut closer to the origin when compared to the $1/N_f$ results. The same singularity location appears at $1/N_f^3$, and these correspond to a UV renormalon singularity in the finite part of the photon two-point function.


I. INTRODUCTION
The discovery of asymptotically safe quantum field theories in four dimensions [1,2] triggered renewed interest in studying the ultraviolet fate of quantum field theories once asymptotic freedom is lost. The original proof of asymptotic safety made use of the Veneziano-Witten large number of flavors and colors limit for a class of gauge-Yukawa theories that displayed perturbatively trustable ultraviolet fixed points. Without scalars, it is impossible to analytically disentangle the ultraviolet fate of asymptotically non-free gauge-fermion theories. Nevertheless one can make progress by analyzing the large N f dynamics of these theories at finite number of colors [3][4][5][6][7][8][9] including again certain type of Yukawa interactions [10][11][12][13]. These studies make use of the large N f resummation techniques to derive the all orders in the 't Hooft coupling beta functions of these theories at order 1/N f . This large N f beta function has several interesting properties including the emergence of singularities undermining the consistency of the expansion, whose physical interpretation remains still to be clarified [14]. In the meantime first principle lattice simulations have begun to explore the large N f dynamics in a systematic manner [15]. It is therefore highly desirable to gain insight into the sub-leading 1/N 2 f corrections. This task, however, turns out to be challenging. The present work constitutes a step forward in this direction by determining these sub-leading corrections for a closed class of diagrams in QED.
To achieve our goal we will make use of the technologies developed in our recent work [16] according to which it is shown that it is possible to reconstruct the 1/N f beta function and its properties using a finite number of coefficients of the perturbative series. We determined the stability of the series and showed that about thirty terms were needed to properly reconstruct the 1/N f beta * Email: dondi@cp3.sdu.dk; ORCID: 0000-0002-6971-2028 † Email: gerald.dunne@uconn.edu; ORCID: 0000-0003-1338-339X ‡ Email: reichert@cp3.sdu.dk; ORCID: 0000-0003-0736-5726 § Email: sannino@cp3.sdu.dk; ORCID: 0000-0003-2361-5326 function up to the leading singularity. The technology includes Padé methods, combined with the study of the large order growth of the perturbative series.
In this paper, we employ this technology to deduce the first complete set of 1/N 2 f and 1/N 3 f corrections for the QED nested diagrams. We discover the emergence of a novel singularity at order 1/N 2 f within this subset of diagrams. The latter appears closer to the origin when compared to the original 1/N f singularity of the full beta function. We refrain from speculating about the physical content of this singularity given that the remaining diagrams are still to be computed. The nature of the singularity is captured by a novel logarithmic branch-cut at a value of the 't Hooft coupling where the beta function remains finite. We also find an intriguing correspondence between the UV renormalons in the finite part of the photon two-point function, appearing at multiples of 3 in the Borel plane, and the leading singularities of the divergent part of the photon two-point function, at order 1/N 2 f and 1/N 3 f . The paper is organized as follows. In Sec. II we fix the notation and introduce the basic building blocks for our computation. This is followed by Sec. III, in which we present details of the computation and uncover the beta function contribution and its leading singularity of the nested diagrams. In Sec. IV we study the appearance of renormalons in the finite part of the photon two-point function. We offer our conclusions in Sec. V. The details of the various computations can be found in the appendices.

II. LARGE N f QED SETUP
In QED with a large number of flavors N f , it is natural to introduce the 't Hooft coupling which we keep fixed when sending N f to infinity. This allows organizing the beta function as a series in 1/N f where each β (k) (K) constitutes itself a perturbative expansion in the 't Hooft coupling K. With this counting, a fermion bubble is of order one and each photon line in a diagram is dressed with n fermion bubbles as depicted in the following diagram We indicate an n-loop bubble chain with a grey blob with an index n referring to the number of fermion bubbles in the photon chain. At the zeroth order in the expansion only the single fermion bubble contributes to the beta function, which reads The first order is given by the diagrams We indicate the sum of these diagrams by a grey square labelled by the same n. The diagrams in (5) were computed for the first time in [3]. The analogous contribution for QCD was computed in [4,6], see [9] for a review. Since the diagrams only contain a single bubble chain, the resulting beta function can be resummed and expressed by a closed integral representation The integrand function is given by This beta function has logarithmic branch cuts at K n = 15 2 + 3n for n ≥ 0. The leading behavior with which it approaches the radius of convergence is Thus, the behavior near the first branch cut is negative, which allows the 1/N f contribution, (6), to cancel the leading contribution, (4). This leads to a zero in the beta function at this order in 1/N f , which has triggered speculations about the existence of a UV fixed point. Similarly in QCD, a negative logarithmic branch cut at K = 3 allows for a zero in the beta function. The potential fixed point in QED has a diverging fermion anomalous mass dimension [17] and is thus considered unphysical. At the analogous fixed point in QCD, the fermion anomalous mass dimension is instead vanishing [17], but the fixed point suffers from glueball operators with diverging anomalous dimensions [18]. In Ref. [18], the authors argued that this fact can be interpreted as an operator decoupling and thus the interacting fixed point might still be physical. The potential existence of the fixed point has triggered already many phenomenological studies [10,19,20]. The viability of the fixed points has been studied on the lattice [15] and with critical point methods [14,21].
Here we go beyond the state of the art by computing part of the full beta function at orders 1/N 2 f and 1/N 3 f . We are interested in whether new singularities can appear at this order that could shrink the overall radius of convergence. The complete knowledge of the full beta function at these orders would be ideal to test the physical nature of the potential fixed points. Given the complexity of the task, we focus here on QED and determine the contributions coming from diagrams of the nested type (to be defined in the next section) to the 1/N 2 f and 1/N 3 f order. We will see that these contributions alone show rather interesting features.

III. NESTED DIAGRAMS
In [16], we have laid the foundations to access crucial information regarding the singular structure of the beta function of the theory by knowing finitely many coefficients of the perturbative series in the coupling K. We showed that to extract the precise radius of convergence and uncover the first singularity in K, roughly the first thirty coefficients of the perturbative expansions are needed. Clearly, the task to extract so many coefficients becomes progressively more demanding when going beyond the leading result. Therefore, although, the procedure can, in principle, be applied to the full 1/N 2 f and 1/N 3 f beta function, already to O(1/N 2 f ) it requires to determine four-loop diagrams and non-planar three-loop diagrams, for which the master integrals with dressed propagators are not known, see Fig. 9 in the appendix for the full set of diagrams.
Fortunately there is one closed set of diagrams, the nested ones, which is gauge and RG-scale independent, and which can be tackled. The associated Feynman diagrams are obtained iterating the 1/N f topologies and are given by We represent the sum of these contributions with a grey hexagon and three indices labelling the number of fermion bubbles on each photon propagator. The full amplitude from these diagrams is given in (A10) of App. A in terms of discrete convolutions of the 1/N f amplitude.    nested (K) reveals that the radius of convergence is K = 3. A Richardson extrapolation is used to accelerate the convergence of the series. Right: The prefactor of the leading large-order behavior is determined to be − 1 2n(n−1)3 n , which leads to a much faster convergence than − 1 2n 2 3 n . See App. B.
There is an additional contribution from a counterterm that makes the result gauge and RG scale independent.
We computed these contributions to the beta function separately up to The full amplitudes for these diagrams are now given in (A11) and (A12) of App. A. Again, together with the counter term contributions, this subset of diagrams is gauge and RG scale independent. We determined these contributions to the beta function up to K 32 and report the coefficients in App. D. Determining the coefficients from the diagrams in (9) and (10) required a significant computational effort. We use the Mathematica package FeynCalc [22] to contract the diagrams and standard multi-loop techniques to evaluate them, see, for example, [23][24][25][26]. The most computation power is needed for the numerical extraction of the divergent and finite part at each loop order. The Mathematica package NumExp [27], which numerically expands hypergeometric functions, turned out to be very useful in this context.
We are now ready to analyze the large-order behavior of the expansion coefficients of β (2) nested (K), in order to extract physical information concerning possible physical singularities. We apply ratio tests, and Darboux's theorem [28][29][30], as well as Padé methods to access in-formation about the leading singular structure of the associated beta function. To keep the presentation light we report the details of the methods in App. B and App. C.
A. Leading singularity of β To extract the leading singularity we must first and foremost demonstrate that the number of terms at our disposal is sufficient to see convergence. This is performed by running the ratio test b n+1 /b n , with b n the coefficients of the β (2) nested series. For sufficiently large n the ratio approaches the inverse radius of convergence. We report the results in the left panel of Fig. 1 and give the detailed analysis in App. B. To accelerate the convergence we further employed Richardson extrapolation. From the left plot of Fig. 1, one learns that about thirty coefficients are sufficient to approach the convergence of the series.
One learns that the radius of convergence is K = 3. From the right plot of Fig.1 we further learn that in total the coefficients have leading decay b n ∼ −1/(2 · 3 n n 2 ). This allows us to subtract from each b n its leading large n contribution, allowing us to determine its sub-leading large n behavior which is discovered to go as −1/(2 · 3 n n 3 ). This trend repeats after each subtraction and therefore we can determine the large-order behavior of the nested QED beta function .
It is interesting that with just 30 expansion coefficients we can clearly distinguish the correct sub-leading large n behavior, as shown in the right plot of Fig. 1. Note that the behavior in (11)  singularity once we resum the series This implies that in the vicinity of the leading singularity at K = 3, where the sub-leading terms are analytic at K = 3. This is a remarkable result for a number of reasons: i) The nested QED beta function at O(1/N 2 f ) has a logarithmic branch cut at K = 3 while remaining finite there. ii) Once the large order behavior is subtracted the remaining contribution at K = 3 is regular. iii) The singularity occurs at a value of K which is smaller than the leading singularity of QED occurring for K = 15/2 at the first order in 1/N f . iv) The singularity occurs at the same value as the leading one for QCD.
It is worth mentioning that because of the simple structure of the function multiplying the logarithmic singularity in (12), it is possible to confirm this behavior by analyzing the series obtained by a second derivative with respect to K of the nested beta function. This is so because the logarithmic singularity turns into a simple pole that it is more easily accessed by the test. In fact, in this case, the onset of the converge occurs already for ∼ K 28 as detailed in Fig. 10 of App. B.
We can now use Padé approximants to deduce the full form of the nested beta function up to the singularity, which is depicted in Fig. 2. More information about the Padé analysis is provided in App. C. [21,21] P [22,22] CP [21,21] CP [22,22] (14). The Padé approximants seem to hint towards a pole or a branch-cut at K = 15 2 .
B. Sub-leading singularity of β We now turn to the sub-leading singular behavior of β (2) nested . To that end we subtract the leading logarithmic branch cut behavior from the nested beta functioñ With the coefficients at hand no further singular behavior is revealed by the ratio test. From the naive expectation that the next singularity arises at K = 15 2 , which is the point where the leading order full beta function is singular, we estimate that roughly 55 coefficients would be needed, which goes beyond the scope and resources of this investigation.
Since the ratio test is not sufficiently precise to display any sub-leading singular behavior, we use Padé methods instead. The three highest Padé approximants are displayed in Fig. 3. They converge well up to K ≈ 7 and thus we know that there is no singular behavior for K < 7. For K > 7, some approximants show singularities and we suspect to find the next singular behavior at K = 15 2 , since this was the location of the logarithmic branch cut at 1/N f .
Having a hint for the location of the next singularity, we now employ the conformal Padé method which is expected to be more accurate in this case. The conformal Padé takes as input the location of the singularity and maps the series to the unit disc. After the conformal transformation, one re-expands the function and applies the standard Padé method and, in the end, one inverts back the conformal transformation. The Padé-Conformal method is described in detail in App. C. In [31], the improvement by a suitable conformal transformation was studied on the example of the Painlevé I equation. We tested the conformal Padé method for the leading singularity at 1/N 2 f and it indeed gives improved results, for details see App. C, and Fig. 3 and Fig. 11. [14,14] P [15,15] CP [14,14] CP [15,15]  For the sub-leading singularity, the conformal Padé results are displayed in Fig. 3. We therefore feel confident to have captured both the leading branch-cut singularity at K = 3, and the sub-leading one occurring at K = 15/2. However, for the latter the available data (i.e. the available expansion coefficients) is not sufficient to determine the precise nature of this sub-leading singularity.

C. Leading singular behavior of β (3) nested
Here we use directly the Padé method to infer the location of the first singularity given the fewer computable coefficients than needed for the ratio test. The highestorder Padé approximants are displayed in Fig. 4 showing convergence up to K ≈ 2.7 and hence we can exclude any singular behavior in this range. It is reasonable to expect the first singularity to occur for K = 3, since all Padé approximants have a pole shortly after K = 3. We therefore apply the conformal Padé method (see App. C) to get a more accurate representation of the associated beta function up to the singularity, and plot it in Fig. 4.

IV. FINITE PARTS AND RENORMALONS
So far, we discussed the divergent part of the photon two-point correlation function because it is directly related to the beta function of the theory. The finite part, however, also has an interesting story to tell, the renormalon story [32]. Indeed, the class of diagrams contributing to the 1/N f correlator are the ones that were originally considered as a renormalon source. Renormalons emerge as singularities of the Borel transform of the finite part of the correlation function and produce a factorially growing series not associated with diagram proliferation. Resummation methods for the finite parts in QED at order 1/N f have been explored in [33] and shown not to  Figure 6. The 1/N 2 f coefficients of the 1/ part of the nested diagrams and corresponding counter term are each factorially divergent. As these are not separately RG-independent we choose µ 2 = −p 2 /4π, p 2 being the external momentum. Their sum produces a RG-independent convergent series according to the large order behavior given in (11). is confirmed at the 1/N f order.
Remarkably, even with our limited number of expansion coefficients, we observe the same factorial growth and singularity structure at sub-leading orders 1/N 2 f and 1/N 3 f of the divergent part of the nested diagrams. See Fig. 6 for the rate of growth of the expansion coefficients in the 1/N 2 f case. The coefficients for the nested diagrams and the corresponding counterterm both grow factorially fast, but with the same rate and opposite signs, in such a way that the factorial growth cancels, leaving coefficients of a rapidly convergent expansion, as shown in (11).
Furthermore, in Fig. 7 and Fig. 8 we plot the real parts of the Borel transform of the finite parts at orders 1/N 2 f and 1/N 3 f , respectively, as shown in (A22). These plots should be compared with the 1/N f result in Fig. 5. At orders 1/N 2 f and 1/N 3 f we do not have the luxury of closed-form expressions, but with the limited number of expansion coefficients (see App. D) our Borel transforms, after conformal mapping and Padé approximation (see App. C), clearly reveal singularities at t = +3 and t = −6, with strong indications of a further singularity at t = −9. With more coefficients, one would be able to resolve even more Borel singularities. Note that without the conformal map, the Padé approximation to the Borel transform cannot see any physical singularities beyond the leading ones, because Padé tries to represent the leading branch cut with an array of poles and zeros, which have no physical content beyond a crude representation of the cut, and these unphysical poles therefore obscure further physical singularities. On the other hand, the Padé approximation to the conformally-mapped expansion, as described in App. C, does not place unphysical singularities on the cut [31,34], so higher physical singu- larities can be seen. These results suggest a relation between the leading order renormalon factorial growth and the singularities in the divergent part of the nested diagrams. This QED Borel structure suggests singularities on the positive real axis associated with UV renormalons, and singularities on the negative real axis associated with IR renormalons [32].

V. CONCLUSIONS AND OUTLOOK
In this paper we have determined the contribution to the QED beta function stemming from the gauge and RG-scale independent class of nested diagrams to order 1/N 2 f and 1/N 3 f , resolving their leading singularity structure. We have shown that: 1) The nested beta function at O(1/N 2 f ) has a new logarithmic branch cut at K = 3, coinciding with the QCD branch cut at O(1/N f ). The nested beta function is finite at the branch cut.
2) The next singularity of the nested beta function appears at K = 15 2 . However, we do not have enough perturbative data to fully characterize its nature.
3) The first singularity of the nested beta function at O(1/N 3 f ) appears at K = 3, but its nature remains to be determined. An important message from this analysis is that it is indeed feasible, as proposed in [16], to use a finite number of perturbative expansion terms to probe certain nonperturbative properties at higher orders of the large N f expansion; we do not require the closed-form expressions which are available at leading order. This suggests a new strategy for studying physical properties of large N f expansions.
Acknowledgements We would like to thank Anders Thomsen, Tommi Alanne and Simone Blasi for discussions on the finite parts resummations. This work is partially supported by the Danish National Research Foundation grant DNRF:90, and is based upon work supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics under Award Number DE-SC0010339.
Appendix A: Renormalization procedure and computation of nested diagrams In this appendix, we detail the applied renormalization procedure. We apply dimensional regularisation in d = 4 − dimensions. The 1PI photon two-point function is parameterised by We expand the renormalization of the coupling Z K = K/K 0 , where K 0 is the bare 't Hooft coupling, as well as Π in orders of N f Each Z n can be written as a series in 1/ The simple pole in given by Z n determines the beta function at each order in 1/N f . The latter is given by We use a minimal subtraction scheme and thus the renormalization condition is We expand this equation in orders of 1/N f . From this, also using the fact that Π 0 is linear in the bare coupling, we obtain Here, Π 0 is precisely the single fermion bubble and thus Z 0 = 1 − 2K 3 . Π 1 is given by the diagrams displayed in (5). In Z 2 , the first term contains the factor Π 2 , which is precisely the diagrams displayed in Fig. 9. The second term can be viewed as the 1/N f diagrams (5) with a 1/N f counter term insertion. In Z 3 , the first term is again given by 1/N 3 f diagrams, while the second and third term can be viewed as lower-order diagrams with counter term insertions.

Nested diagrams
We now display the structure of the nested diagrams. For this it is useful to write the 1/N f contribution to the photon two-point function, i.e., the diagram given in (5), in an expansion in loop orders Here, Π (n) 1 corresponds to the contribution with n inserted fermion bubbles, µ is the RG scale, and G 0 is the one-bubble amplitude given by This notation allows us to write down the nested amplitudes in a convenient way. For the 1/N 2 f nested diagrams, displayed in (9), we assign n and m fermion bubbles to the outer photon propagators and fermion bubbles to the inner photon propagator. Then the amplitude is given by where we used m,n=0 f (m + n) = k=0 (k + 1)f (k). In straight analogy, we write down the nested amplitudes for 1/N 3 f . For the two diagrams in the first line of (10), the amplitude reads where we used ,m,n=0 f ( + m + n) = k=0 1 2 (k + 1)(k + 2)f (k). The two diagrams in the second line of (10) result in the amplitude Π 3,nested,diag2 = K 6 0 ∞ l,m,n,p,q=0 We close this appendix with a short discussion of the diagrams at O(1/N 2 f ) which are not computed in this paper. The full set of diagrams contributing to the beta function at O(1/N 2 f ) is displayed in Fig. 9. For all diagrams in the first line of Fig. 9, the corresponding master integral is known: they all are topologically still two-loop or even one-loop diagrams. All diagrams in the second line are topologically three-loop diagrams, except the last one, which is a topological four-loop diagram. The most challenging diagrams are the last two diagrams in the second line of Fig. 9: the first is a non-planar topological three-loop diagram with two bubble chains, while the second one is a topological four-loop diagram with three bubble chains.

Finite parts
We now detail the computation of the finite part of the regularised two-point function Z 0 Π 1 and its corresponding Borel transform. We write the amplitude of the two-point function schematically as Where the function H is regular in n for constant and in for constant n . Note that this is a different representation of Π 1 than in (A8). In the following, we use the expansion of H in n as well as in , which we denote by We plug this into (A13) and also expand K 0 in . We obtain where we introduced This is computed as We denote the finite part, i.e., the limit → 0, of the amplitude (A15) as F (1) . This can be computed as The first part is a convergent series, while the second one is asymptotic. For this reason, only the second part can contribute to singularities in the corresponding Borel transform. We define the Borel transform here by We then obtain This is the same formula as in [33] adapted to our notation. The function H is given by with H(0, 0) = 1 and M = − 4πµ 2 p 2 . Note that the singularities at t = 0 and t = 3 in the first and third term cancel. Only the term with the hypergeometric function is contributing to the singularity at t = 3. We display the Borel transform of the finite part in Fig. 5.
The finite parts of the 1/N 2 f and 1/N 3 f contributions are defined as No closed-form resummation is possible for these contributions. We perform the series expansion in K numerically and provide the coefficients in App. D. We computed the nested part of F (2) up to K 32 and the nested part of F (3) up to K 28 . The Borel transforms are obtained by dividing out the leading factorial growth. This leads to the following definition of the respective Borel transforms These Borel transforms are analytically continued using the Conformal-Padé method described in App. C and plotted in Fig. 7 and Fig. 8. where φ(z) and ψ(z) are analytic near z 0 , then the Taylor expansion coefficients of f (z) at the origin have large-order growth These results can be used in reverse to find the singularity location z 0 , the exponent p (or to detect logarithmic behavior), and properties of the coefficient function φ(z), from the large-order growth of the expansion coefficients at the origin. We implemented this strategy on the perturbative expansion of β (2) nested (K), with expansion coefficients b n . A simple ratio test suggests that b n+1 /b n → 1/3, see Fig. 1. This can be refined using Richardson extrapolation to accelerate the convergence of the ratio test. Richardson extrapolation is based on the ansatz [35] a n = a + A n where a is the anticipated convergent value. First-order Richardson extrapolation is obtained by setting all parameters beyond 1/n to zero, i.e., B = C = . . . = 0. The evaluation at n and n + 1 yields R (1) a n ≡ a = (n + 1)a n+1 − na n .
In Fig. 1 we display the effects of the second-order Richardson extrapolation on the enhancement of the convergence of the ratio test series, clearly indicating convergence to 1/3, indicating the existence of a singularity at K * = 3, and hence a radius of convergence equal to 3. Given z 0 , we can now fit the growth of the coefficients b n to the branch cut forms in (B2) and (B4). This can be done by studying the sub-leading behavior of the ratio test. This reveals that the ratio behaves as where the precise sub-leading coefficient, − 2 3 , can be extracted using Richardson acceleration once again. This indicates logarithmic behavior, as in (B4). Now we can probe this further to deduce information about the analytic function φ(z) multiplying the logarithmic branch cut. The result (B8) implies that φ(3) = 0 and φ (3) = 1/6. Analysis of further sub-leading corrections indicate that all higher derivatives of φ(z) vanish at z = 3. This leads to the result (2) nested exact P [15,15] CP [15,15]  for the logarithmic branch cut in (12). To confirm this result, we plot in Fig. 1 a precise test of the deduced largeorder behavior of the b n coefficients in (11). The agreement is excellent. Note that with 44 coefficients we can clearly distinguish between b n ∼ − 1 2 1 3 n 1 n(n−1) , and the cruder estimate b n ∼ − 1 2 1 3 n 1 n 2 . An interesting further consistency test of the logarithmic form of the singularity is to differentiate (twice) the nested beta function β (2) nested (K), and then apply the Darboux analysis. The resulting function has a simple pole, which is easy to detect with a ratio test, see Fig. 10 for the convergence of the ratio test in this case.

Appendix C: Padé versus Padé-Conformal
Padé approximants provide well-known analytic continuations of truncated series expansions and are widely used in physical applications [35]. It has further been observed empirically that combining Padé approximants with conformal maps often yields further improved precision [31,36,37]. This improved precision is explained and quantified in [34]. The Padé-Conformal analytic continuation procedure for a truncated series in the presence of a branch cut is: (i) first, make a conformal transformation from the cut complex plane to the unit disk; (ii) second, re-expand to the same order inside the conformal disk; (iii) third, make a Padé approximation to the resulting series inside the disk; (iv) finally, map back to the original cut plane with the inverse conformal transformation. This procedure is algorithmically straightforward and is provably exponentially more precise than just Padé if there is a cut [31,34].
In the presence of a single cut (interestingly, it does not matter what the precise nature of the cut is, just where it is), the explicit conformal map from the K plane cut along the positive real axis with a branch point at K * , together with its inverse are: The branch cut itself is mapped to the unit circle in the complex z plane. Given K * , which we have determined to be 3, it is now a completely algorithmic procedure to implement this Padé-Conformal extrapolation. The result is much more precise than just making a Padé approximation, especially in the vicinity of the branch point and branch cut. The results are shown in Fig. 3 and Fig. 11 for β (2) nested (K), and in Fig. 4 for β nested (K). In the presence of two cuts along the real axis, as occurs for the Borel analysis in Sec. IV, we use the conformal map from the Borel t plane cut along the positive real axis t ∈ [b, ∞) and along the negative real axis t ∈ (−∞, −a], to the unit disk in the z plane: For the expansions of the finite parts in App. A 2, the leading singularities are at t = +3 and t = −6, so we choose b = 3 and a = 6. We map the Borel transform to the unit conformal disk in the z plane, re-expand, and map back again to the Borel t plane. The resulting plots for the 1/N 2 f and 1/N 3 f Borel transforms are shown in Fig. 7 and Fig. 8, respectively. Note that the conformal mapping is crucial for revealing the existence of sub-leading Borel singularities [31,34].
The coefficients of all finite parts grow factorially, as expected. From the above expressions for the finite parts we can already see that the factorial rate of growth of the coefficients is comparable at O(1/N f ), O(1/N 2 f ), and O(1/N 3 f ). Moreover, since all the coefficients are positive, we deduce that the leading Borel singularity should be on the positive real axis. Indeed, our further analysis shows that this leading singularity is at t = 3, for each order of the large N f expansion. See Fig. 5, Fig. 7 and Fig. 8.