Global Wilson-Fisher fixed points

The Wilson-Fisher fixed point with $O(N)$ universality in three dimensions is studied using the renormalisation group. It is shown how a combination of analytical and numerical techniques determine global fixed point solutions to leading order in the derivative expansion for real or purely imaginary fields with moderate numerical effort. Universal and non-universal quantitites such as scaling exponents and mass ratios are computed, for all $N$, together with local fixed point coordinates, radii of convergence, and parameters which control the asymptotic behaviour of the effective action. We also explain when and why finite-$N$ results do not converge pointwise towards the exact infinite-$N$ limit. In the regime of purely imaginary fields, a new link between singularities of fixed point effective actions and singularities of their counterparts by Polchinski are established. Implications for other theories are indicated.


I. INTRODUCTION
Wilson-Fisher fixed points provide a paradigm for continuous phase transitions and scale invariance in quantum and statistical field theory [1]. An  In specific limits such as large-N , or close to four dimensions [2], Wilson-Fisher fixed points are under exact perturbative control and allow for a complete solution of the theory [3]. At finite N , however, Wilson-Fisher fixed points are often strongly coupled and exact analytical solutions are not at hand [4,5]. Instead, one has to resort to non-perturbative tools including numerical simulations on the lattice [6], methods from conformal field theory [7] or, as will be done here, methods from the renormalisation group [8][9][10][11].
Ideas to understand interacting fixed points analytically with the help of recursive relations have recently been developed in [28,37,38]. The virtue of this approach is that all local couplings, polynomial or otherwise, can be uniquely determined in terms of a few free parameters owing to the iterative diagramatic structure of Wilsonian renormalisation group flows [49,50]. For O(N ) symmetric scalar field theories, exact recursive relations were obtained and exploited to specify the complete fixed point action using local expansions about small, large, or purely imaginary fields [28]. At infinite N , it has also been shown that local expansions are sufficient to uniquely determine the Wilson-Fisher fixed point globally, leading to closed expressions for the fixed point action. In this limit, free parameters become uniquely specified because the radii of convergence of small and large field expansions overlap. In this paper, we are interested in the global Wilson-Fisher fixed point at finite N . The main new ingredient are the fluctuations of the radial mode which modify the set-up: local recursive relations are now characterised by a larger number of free parameters, and the singularity structure in the complex field plane is altered. As a result, we observe that the radii of convergence for small and large field expansions no longer overlap and complete analytical solutions remain out of reach. Still, global fixed points can be determined accurately, and with moderate effort, by combining analytical insights with numerical integration [10,11]. Our method identifies the few free parameters numerically, and offers largely analytical results for the fixed point action for all fields within the respective radii of convergence, together with numerical solutions in those regions of field space which cannot be reached via small or large field expansions.
Our strategy is put to work for O(N ) symmetric scalar field theories covering the universality classes N = −2, −1, 0, 1, 2, 3, . . . , 10, 20, 30, . . . , 10 2 , 10 3 and 10 4 , and to leading order in the derivative expansion. We compute global fixed point potentials, leading and subleading scaling exponents, mass ratios, local fixed point coordinates, radii of convergence, as well as the parame-ters which control local expansions of the effective action. In the regime of purely imaginary fields, a novel relation between singularities of the Wilson-Polchinski RG [12] and singularities of their Legendre-transformed counterparts [13] is established. Special emphasis is put on the accuracy of numerical procedures and the validity of local expansions, both of which are monitored carefully in all steps involved.
The outline of the paper is as follows. Our setup and computational strategy is introduced in Sec. II, results are detailed in Sec. III, and Sec. IV concludes with a discussion and an outlook.

II. RENORMALISATION GROUP
In this section, we shall recall the RG equations, introduce our approximations, and explain the computational strategy.

A. Wilsonian renormalisation group
This work employs the method of functional renormalisation based on a Wilsonian version of the path integral where parts of the fluctuations have been integrated out, [12-14, 16, 51, 52]. We consider Euclidean scalar field theories with partition function where S denotes the classical action and J an external current. The Wilsonian cutoff term is given by with R k (q 2 → 0) > 0 for q 2 /k 2 → 0 and R k (q 2 ) → 0 for k 2 /q 2 → 0 to ensure that R k acts as an IR momentum cutoff [16][17][18]. As such, the partition function (1) interpolates between a microscopic (classical) theory (k → ∞) and the full physical theory (k → 0). The main object of our investigation will be the"flowing" effective action Γ k . It relates to (1) via a Legendre , where φ = ϕ J denotes the expectation value of the quantum field. The renormalisation group scale-dependence of Γ k is given by an exact functional identity [13] (see also [14,52]), It expresses the change of scale for the effective action Γ k with an operator trace over the full propagator multiplied with the scale derivative of the cutoff itself. At weak coupling, iterative solutions reproduce the perturbative loop expansion [49,50]. The exact RG flow (3) also relates to the well-known Wilson-Polchinski flow [12] by means of a Legendre transformation, and reduces to the Callan-Symanzik equation in the limit where R k (q 2 ) becomes a momentum-independent mass term [20]. The right-hand side (RHS) of the flow (3) is local in field and momentum space implying that the change of Γ k at momentum scale k is governed by fluctuations with momenta of the order of k. Optimised choices for the regulator term [5,16,17] allow for analytic flows and an improved convergence of systematic approximations [8]. Approximations of the exact flow (3) are the central input in this study.

B. Local potential approximation
The local potential approximation (LPA) expresses the O(N )-symmetric 3d effective action for real scalar fields φ by a standard kinetic term and a general effective potential, where φ denotes a vector of N scalar fields. In the infinite N limit, this approximation becomes exact. After inserting (4) into (3), the scale dependence of the dimensionless effective potential u(ρ) = V k (φ)/k 3 is determined by the partial differential equation [15], Here, t = ln k is the logarithmic scale parameter, and k the Wilsonian RG scale. The dimensionless field variable ρ = 1 2 φ 2 /k automatically takes into account the invariance of the action under φ → −φ. The four terms on the RHS originate from the canonical scaling of the potential, the field, the fluctuations of the N − 1 transversal (Goldstone) modes, and the fluctuations from the longitudinal (radial) mode, respectively. The functions I[x] parametrize threshold effects owing to the decoupling of heavy modes. In general, they take values of order unity for small argument, and vanish in the limit where field-dependent massesm 2 grow large (I[x] → 0 for x → ∞, where x =m 2 /k 2 stands for the square of the field-dependent mass in units of the RG scale). The functions I[x] also depend on the choice for the Wilsonian momentum cutoff [25]. It takes the simple analytical form for an optimised choice of the regulator, following [16,17,25]; see [28,53] for exact large-N solutions. The numerical factor A = 2/(d L d ) arises from the angular integration over loop momenta, and L d = (4π) d/2 Γ(d/2) denotes the d-dimensional loop factor [28]. At a fixed point solution, universal scaling exponents are independent of the numerical constant A. This constant can be absorbed into the potential and the fields by the rescaling which is equivalent to setting A = 1 in (6). The benefit of this normalisation is that couplings are now measured in units of the appropriate loop factors suggested by naive dimensional analysis [54]. All our numerical data below is obtained for this choice. Since the zero point energy is not determined by the RG flow (5), it is sufficient to investigate the flow of its first derivative, which is given by We are interested in the Wilson-Fisher fixed point solutions ∂ t u = 0 of (8) for all N and all fields. The fixed point potential obeys an ordinary second order non-linear differential equation where we recall to have set A = 1. In principle, (9) uniquely determines the Wilson-Fisher fixed point. In practice, integrating (9) numerically starting from small field values is difficult because the equation is stiff, e.g. [55]. Therefore we resort to a combination of analytical and numerical techniques to identify the unique fixed point solution and its universal properties.

C. Small field expansions
We briefly introduce the polynomial expansions of the scaling potential around its minimum in a phase with spontaneous symmetry breaking. It is well-known that scaling exponents can reliably be computed within a small-field expansion [25]. Here we make the polynomial ansatz [56][57][58], which we approximate at a maximum order n = M for the scaling potential in the vicinity of the potential minimum ρ 0 with u (ρ = ρ 0 ) = 0. We refer to this ansatz as expansion A. This type of ansatz has previously been studied in many works, e.g. [15,25,55,59,60]. After insertion of (10) into the flow equation (5) one obtains a set of coupled ordinary differential equations for the couplings λ n and ρ 0 . The fixed point conditions become ∂ t ρ 0 = 0 together with ∂ t λ n ≡ β n ({λ i }) = 0 for all n. For the expansion A, these equations can be solved algebraically, leading to the exact recursive relation [28] where we use the short-hand notation When solved iteratively, the result (11) provides us with expressions for all couplings {λ n , n ≥ 3} in (10) as explicit functions of the vacuum expectation value and the quartic self-coupling, the two couplings which remain undetermined by the recursive solution. Results for the two couplings (12) which remain undetermined by (11) are given below, for all N . The series (10) is closely linked to an expansion in the "symmetric phase" about vanishing field which we denote as expansion B. Expressions similar to (11) are found for the fixed point couplings in (13). The relevant expressions read [28] λ s,n+1 = 1 When solved iteratively, the expressions (14) provide us with explicit results for all couplings {λ s,n , n ≥ 2} as functions of the mass term at vanishing field, Results for the critical mass parameter (15) for all N will be reported below. Notice that (14) only depends on one unknown parameter (15), whereas (11) depends on two unknowns (12). In expansion point in field space parameter minimum vanishing field large real large imaginary info infinite N : none m 2 γ ζ Ref. [28] finite N : ρ 0 and λ 2 m 2 γ ζ andζ this work either case, the remaining parameters must be determined by some other means. However, it has been observed previously that the expansion (10) converges more rapidly than expansion (13), see [25,55,[59][60][61]. For this reason, in the remaining part of this paper, we will focus on the expansion (10) with (11). Local polynomial expansions such as (10) and (13) necessarily have a finite radius of convergence R, typically of the order of the vacuum expectation value [25,28] In order to extend the fixed point solution for all fields, additional expansions are required.

D. Large field expansions
To obtain the fixed point solution for all fields, we additionally consider expansions for asymptotically large real (ρ → ∞) and imaginary (ρ → −∞) field amplitudes. In the limit ρ → ∞, the flow (8) admits an expansion of the form [28] u (ρ) = γρ 2 which we approximate at some finite maximum order n = M . We refer to this ansatz as expansion C. Interestingly, for all N and using the beta functions for the couplings γ n , we can solve the fixed point condition recursively to provide us with unique functions for the couplings [28] γ n = γ n (γ) (18) in terms of the leading order coefficient γ, Of this one-parameter family of solutions, only a single one is related to the Wilson-Fisher fixed point. Below, the values (19) are determined numerically for each universality class. We are also interested in the behaviour of the Wilson-Fisher fixed point solution for large negative field squared ρ → −∞, corresponding to purely imaginary fields. In this limit, the fixed point potential has an expansion of the form [28] which we approximate at a finite maximal order m = M . We refer to this ansatz as expansion D.
In [28], it has been established by exploiting the beta-functions for all couplings ζ n,m recursively that the putative fixed point couplings are expressed as functions of two independent couplings ζ m,n = ζ m,n (ζ,ζ) .
Of this two-parameter family of fixed point candidates, only a few, including the Gaussian and the Wilson-Fisher fixed point, extend over all field space. The independent parameters in (21) are given by the coefficients see (20). Values for these will be determined numerically for all N .

E. Computational strategy
A full, global fixed point solution for all values of the field amplitude ρ and predictions for the critical exponents can be obtained with the following two-step strategy. In Step 1, we determine all fixed point couplings in a local approximation, using expansion A owing to its fast convergence. For a given order in the polynomial approximation M around the minimum (10), we obtain the set of couplings (11). Next, we must determine the remaining free parameters (12) numerically. To order M in the polynomial expansion, we impose additional boundary conditions on some of the higher order couplings which are not contained in the approximation up to order M . Specifically, we use see [10,25], which provides us with two additional equations for the two unknowns (12). The constraint (23) is equivalent to omitting all couplings λ n with n ≥ M + 1 in the effective action from the outset. This strategy has the form of a bootstrap, where (23) serves as additional input for closure. Note also that the full non-perturbative fixed point will not obey (23) exactly. Hence, it is important to establish stability of our fixed point search under variations in the boundary condition and extensions of the approximation order M → M + 1.
The numerical accuracy of the fixed point solution, for all fields, is monitored by the quantity which provides a measure for the deviation from scaling. We have N acc → ∞ for the exact solution, for all fields. Finally, the universal critical exponent ν and the sub-leading scaling exponents ω n are extracted as the eigenvalues of the stability matrix at criticality, This technique has been successfully applied previously [25] including for anti-symmetric corrections to scaling [61], and extended for a high-accuracy determination of Ising exponents in [10].
Here we adopt it to compute the scaling exponents for all N ranging between N = −2 up to N = 10 4 . The accuracy in a scaling exponent X (or a critical coupling) with increasing order in the approximation is monitored via On the right-hand side, X n denotes the nth order approximation for the quantity X, and X M denotes the highest order value thereof. Hence, (26) counts the number of decimal places which agree between the approximated and the highest order result. By construction we expect that X M approaches the exact solution for M → ∞ with N X → ∞.
In order to connect the local fixed point at small fields with the asymptotic expansions at large fields, we have to determine the values for those parameters which have remained undetermined within the large real (19) or large imaginary field solution (21), see Tab. 1. Provided the radii of convergence of the small field polynomial expansion u small and the large field expansion u large overlap, it is possible to fix the coefficients (27) in the respective overlapping regimes by requiring u small − u large = 0 to within the desired accuracy. In general, we find that the asymptotic expansions for large fields converge, albeit slowly. While the radii of convergence do overlap in the limit of infinite N [28], this cannot be guaranteed at finite N . Furthermore, with decreasing N , a very high order in the expansion would be required to potentially achieve overlapping radii of convergence with the small field regime and good estimates for the unknown parameters (27). Therefore, as Step 2, we integrate the differential equation (9) at the fixed point numerically to bridge the gap between the small and the large field expansions. This will also provide numerical values for the parameters (27). Integrating the stiff differential equation (9) from the vacuum expectation value towards large fields is numerically unstable and can terminate e.g. in a moving singularity. The instability can be dealt with by integrating into the opposite direction (from large to small fields), which requires a "shoot and relax"-type strategy. To that end, we choose field values ρ p1 and ρ p2 with 0 < ρ p1 < ρ p2 and well within the range of validity of the small (ρ p1 ) and large field expansion (ρ p2 ), respectively. Using initially a trial value for the parameter γ, we choose the starting point (u large , u large ) at ρ = ρ p2 to integrate the differential equation (9) from ρ p2 down to ρ p1 , thereby bridging the two asymptotic expansions. In this direction the integration is stable. The boundary condition for the parameter γ at ρ = ρ p2 is then adapted iteratively until the result from the numerical integration u num at ρ p1 agrees with the result u small from the polynomial small field approximation to within the desired accuracy, say at the matching point. To ensure stability in the result, we varied all matching points over a wide range and confirmed that our final values for γ are independent of the procedure. In the opposite regime, for large negative ρ, we determine the coefficients of the asymptotic approximation ζ 1,0 and ζ 4,0 in a two-step procedure. We first chose a point ρ m < 0 for which the polynomially approximated potential u small is reliable with, typically, N acc (ρ m ) > 10. The starting values (u small , u small ) at ρ m are used for the numerical integration of (9) towards large negative values ρ < ρ m . In this direction, the numerical integration is stable. At sufficiently large negative field values, the leading order behaviour of the solution is solely governed by the coefficient ζ 1,0 .
The sensitivity to the coefficient ζ 4,0 is suppressed as (−ρ) −3/2 . Hence the numerical solution is used to first determine ζ 1,0 from a fit to the asymptotic expansion at −ρ ≈ 10 10 . We have confirmed that the result is not sensitive against variations of the matching point.
Finally, the value of the subleading coefficient ζ 4,0 is determined by choosing a second matching point at a significantly smaller purely imaginary field value, together with the numerical result for u (ρ) and the prior determination of ζ 1,0 . Again, we confirm that the result is independent of the details of the procedure, as long as the second matching point is still located within the range of validity of the asymptotic expansion. Alternatively, we note that ζ 4,0 can also be extracted from the asymptotic expansion of the function u + 2ρ u . Here, the dependence on the coefficient ζ 1,0 drops out allowing for a direct access of ζ 4,0 at large negative ρ. We have checked that both procedures give the same result.

III. RESULTS
This section summarises our results for the scaling solutions, the scaling exponents and various other universal and non-universal quantities and their accuracy, the large-N scaling of the theory's parameters, and links with the Wilson-Polchinski RG.

A. Potential
At a fixed point of the renormalisation group, quantities are said to be non-universal provided that their values depend on technical parameters such as the RG scheme. Non-universal quantities cannot be measured in any experiment. This applies to the global shape of scaling potentials and the values of couplings at fixed points, which may vary depending on the choices for the Wilsonian cutoff and the RG scale. Quantities are said to be universal provided that their values solely depend on the characteristics of the universality class such as the dimensionality of space-time, the short-or long-range nature of interactions, and the number and type of degrees of freedom. Universal quantities such as scaling exponents can be measured in experiments. For the 3d O(N ) symmetric scalar theories considered here, the different universality classes are then specified by N . Finally, quantities are said to be superuniversal provided they take identical values irrespective of the universality class.
Plots of the non-universal global potential and its derivatives are given in Figs. 1,2 for a representative selection of values for N . The first plot in Fig. 1 (left panel) shows the results for the complete scaling potential u rescaled as in (7) with A = ρ 0 (N ), for N = 10 n , n = 0, 1, 2, 3 and 4 (the intensity of the line colour increases in this order). Notice that the x-and y-axes are rescaled as x → x/(2 + |x|) and y → y/(2 + y) to allow a display of the global fixed point potential for all real and imaginary fields. The potentials show a single global minimum at ρ 0 for all N and all fields. The plots at finite N are complemented by the solution in the infinite-N limit (red dashed line, cf. [28]). With increasing but finite N , the solutions seem to approach the infinite-N result smoothly (note that the red dashed line is nearly on top of the large-N solution for u/ρ 0 for all ρ and is therefore hardly visible) -see Sec. III C for a detailed discussion of the N -dependence. Our results can further be used to derive the global equation of state and universal amplitude ratios for all N (for an overview, see [4]), following the lines of [62] for the Ising universality class at LPA and beyond (e.g. including corrections through anomalous dimensions).   (31). For display purposes, the x-and y-axes are stretched as x → x/(2 + |x|) and y → y/(2 + y).

B. Mass scales
Next we discuss derivatives of the scaling potential. The RG flow is driven by the two relevant mass scales of the model, the field-dependent Goldstone massm 2 ≡ m 2 k 2 and the field-dependent mass of the radial modeM 2 ≡ M 2 k 2 which, in units of the RG scale, are related to the potential as Either mass term contributes to the flow except at infinite N (N = 1) where the radial (Goldstone) contribution is absent. As can be seen from the flow equation (8) together with (29), the limits signal a putative singularity [63,64]. Dynamically, however, the limit (30) is achieved smoothly without encountering any singularity, as follows from Figs. 1,2 and from the analytical result (20), see [28]. Also, the mass scale ordering m 2 < M 2 for large positive ρ > 0 becomes inverted −1 < M 2 ≤ m 2 once ρ < 0. With increasing N , the radial mass fully dominates the solution for imaginary fields.
Notice that (31) is consistent with the radius of convergence R D ≈ −1.5 determined from the expansion (20) at infinite N [28]. For fields below (31), ρ < ρ c , the infinite-N solution no longer fullfills M 2 ≥ −1. Hence, (31) denotes the smallest possible value for ρ down to which finite-N results may converge pointwise towards the infinite N expressions. Furthermore, from the data, we also observe that the scaling solutions do approach the infinite N result for all fields with ρ > ρ c . Taking the limit N → ∞ for the solutions of the finite-N RG equations (8), we observe that M 2 → −1 in the entire regime ρ < ρ c . This is the convex hull of all functions M 2 for any finite N , defined via a differential equation for u which is solved by for all ρ < ρ c . The coefficient ζ is given by ζ 1,0 in (20), in the limit N → ∞. We observe the absence of subleading terms in (32), indicating that all subleading coefficients ζ m,n with m ≥ 2 in the general result (20) become suppressed parametrically for sufficiently large N .
We conclude that the large-N limit of the flow (8) at finite N does not commute for all fields with the integration of the flow at infinite N . Rather, this limit is qualitatively different from the solution of the large-N equations in the regime for purely complex fields. The culprit for these differences are the flucutations of the longitudinal (radial) mode, which at any finite N interfere with those of the transversal (Goldstone) modes and modify the analytic structure of the fixed point solution.

D. Mass differences
Next, we investigate in more detail the region indicated by the arrow in Fig. 2 (right panel), where the field-dependent mass difference Y ≡ M 2 − m 2 is displayed as a function of X ≡ ρ/ρ 0 . We use the vacuum expectation value ρ 0 to set the "scale" in field space, with X measuring the magnitude of fields in units of the vev, per universality class. Although the functions M 2 and m 2 take quite different shapes for different N , we observe from the data that the functions Y (X) seem to coincide close to the region where X ≈ −2, Fig. 2. An exact coincidence for all N would be an indication for superuniversal behaviour. However, a closer inspection shows that this is not quite the case: for each pair of universality classes (N 1 , N 2 ) with N 1 = N 2 we have computed the field ratio X such that Y (N 1 ) = Y (N 2 ), leading to a collection of coincidence points (X, Y ). For large N 1 and N 2 , the data points show a mild dependence on N with values in the range X ≈ [−1.7, −2.2] and Y ≈ [−0.27, −0.28]. We conclude that within the approximations used here the mass difference does not show superuniversal behaviour other than the trivial one where (X, Y ) = (1, 0) for any N , reflecting equality of masses at vanishing field. Note also that we have excluded the exact infinite N result from the analysis, because the pointwise convergence of finite-N results towards their large-N limit no longer coincide with the exact infinite N result for the range of X values found above (X < ρ c /ρ 0 ), see Sec. III C. It would be interesting to check corrections to this picture beyond the local potential approximation.

E. Field ratios
At the Wilson-Fisher fixed point, the physical masses (29) vanish identically in the infra-red limit k → 0. Furthermore, away from the fixed point both masses always coincide for vanishing field, m 2 (ρ = 0) = M 2 (ρ = 0) < 0. At k = 0, and with increasing field values, both masses increase, with M 2 > m 2 for all fields ρ > 0. The square of the Goldstone mass changes sign at the potential minimum, ρ = ρ 0 ≡ 1 2 φ 2 0 /k 2 and the square of the radial mass changes sign at the inflection point ρ = ρ M ≡ 1 2 φ 2 M /k 2 . This makes their ratio a potentially universal quantity worth being explored, also in the light of our results in the previous section. The ratio also gives access to global aspects of the fixed point solution. At the Wilson-Fisher fixed point in the infinite N limit, we find the ratio of field values at the inflection point and at the potential minimum analytically. It reads where x = −0.26763 · · · is determined as the unique solution of the transcendental equation 0 = 8 + 25x + 15x 2 + 15 √ x(1 + x) 2 arctan √ x, leading to Our numerical results for finite N are displayed in Fig. 5. It is noteworthy that the asymptotic limit (35) is an excellent approximation for (33) for most N down to N ≈ 4, suggesting that (33) is largely super-universal, i.e. independent even of the universality class. Quantitatively, the N -dependence of (33) is much weaker than the N -dependence observed for the universal scaling exponents themselves, and much less so than the fixed point coordinates.

F. Barrier height
In Fig. 5 we also display the "barrier height" -refering to the difference between the potential at vanishing field and at the global minimum of the fixed point solution -as a function of the universality class and in units of the vacuum expectation value The dimensionless parameter (36) supplies us with global information of the fixed point potential. It is of order unity for all N , and characterises the universality class. Note that the physical barrier height is given by k 3 (u(ρ 0 ) − u(ρ)), and vanishes in the infrared limit k → 0 as it must at a second order phase transition. In the infinite-N limit, we have  Table 2. Numerical results for the Wilson-Fisher fixed point for all N parametrised in terms of the (dimensionless) vacuum expectation value ρ 0 , the first few polynomial couplings at the potential minimum, and the mass term squared at vanishing field. For N = −2, couplings have been scaled with appropriate powers of (N + 2). The quoted digits are significant according to the criterion detailed in Sec. III G. We have cut the number of quoted digits at a maximum of 12.
where x = −0.388347 · · · uniquely solves the transcendental equation H(x) = −1, with Quantitatively, we therefore have As can be seen from Fig. 5, the N -dependence of (36) is much more pronounced than the one for the scaling exponents and for the field ratio, in particular for small N .

G. Polynomial couplings
The coordinates of the vacuum expectation value ρ 0 , the dimensionless mass parameter at the potential minimum m 2 = u (0), and the polynomial fixed point couplings λ i are summarised in  Tab. 2 (see also Fig. 3). We restrict ourselves to presenting only the leading and the first few subleading couplings and recall that perturbative loop factors have been scaled into the couplings. The order M up to which we approximate the expansion around the minimum, is also given in the table. The quoted digits are those which did not change for at least four preceeding orders of the approximation (that is, between M − 4 and M ). In the last row of Tab. 2 we also provide the N → ∞ limits, empirically determined for the couplings λ i based on (40) in the case of the scaling exponents, and based on [28] for m 2 . For N = −2 the potential's minimum resides exactly at the origin, ρ 0 = 0. As N increases ρ 0 moves towards larger values, asymptotically reaching ρ 0 = N for 1/N → 0. Hence, ρ 0 (N ) sets a "scale" for the overall normalisation of couplings with N . For this reason, it is natural to use ρ 0 (N ) as a normalisation factor for the results. For the purpose of Tab. 2 we therefore have rescaled all couplings with (appropriate powers of) N + 2. In this way all couplings achieve a finite limit even for large N . 1 By definition, the scaled couplings vanish in the limit N → −2. Quartic couplings are of order unity. By naive dimensional analysis [54], the fixed point is boarderline strongly coupled for most N . We also observe that couplings reach maxima for moderate values of N , and that they approach their finite large-N limits from above.

H. Exponents
Exact results for the scaling exponents in 3d are known in the infinite-N limit [65,66] where the universal eigenvalues θ approach those of the spherical model. One finds where n = 0, 1, 2, · · · takes integer non-negative values. The sole negative eigenmode is related to the exponent ν by ν = −1/θ 0 , and the sub-leading corrections-to-scaling exponents are ω n = θ n . In the vicinity of N = −2, the exact result reads [66,67] Numerical results for exponents at finite N have been given in [18,25] with an accuracy up to six significant digits to leading order in the derivative expansion. 2 In [10], the accuracy has been extended up to 14 significant digits for the Ising universality class. The present study provides the exponents for any N within −2 < N < ∞, and enhances the accuracy up to about 30 digits in the result. Our numerical results all smoothly approach their known asymptotic value as N is increased. This is illustrated in Figs. 3, 4. The data underlying these plots is summarised in Tab. 3. It is worth pointing out that the Wilsonian regulator function R k in (3) is key to help enhancing convergence and stability of results [16,17,25,66]. Moreover, universal exponents show a characteristic "cusp-like" structure as a function of the RG scheme, with best results located exactly at the cusp [68]. Incidentally, cusp-like structures have also been observed within the conformal boostrap technique [7]. For our models, the cusp is achieved for the optimised cutoff R k ∼ (k 2 −q 2 )θ(k 2 −q 2 ) [17]. Hence, results are "as good as it gets" to leading order in a derivative expansion, in exact agreement with results from Polchinski's equation [8,12]. For completeness, we compare our findings to leading order in a derivative expansion with the best estimates to date based on Monte Carlo (MC) simulations and the conformal bootstrap where impressive levels of accuracy have been achieved as of late [69][70][71][72]. Specficially, Ising exponents ν = 0.629971(4) and η = 0.0362978 (20) have been found recently with the conformal bootstrap technique [71]. MC simulations have found ν = 0.58759700(40) for entangled polymers (N = 0) [69]; ν = 0.63002(10) and η = 0.03627(10) for Ising universality (N = 1) [70], and ν = 0.6717(1) and η = 0.0381(2) for the XY model (N = 2) [72]. LPA results for ν agree with these to within 0.75% (N = 0), 3% (N = 1), and 5% (N = 2), respectively. Subleading scaling exponents such as ω [66] or antisymmetric corrections to scaling ω 5 [61] are more sensitive to higher derivative interactions [5] and deviate more strongly in LPA [66] (often in the 25%-40% range) provided they are correlated with ν (if not, much smaller deviations are observed [66]). Sub-percent accuracy of RG results is achieved through suitably "optimised" higher orders in the derivative and/or vertex expansion; see [5,[73][74][75] for results up to fourth order in the derivative expansion.

I. Large-field asymptotics
In Tab. 4, we present our results for the expansion coefficients (27), as obtained from the numerical integration of the flow equation, using the techniques discussed in Sec. II E. For the determination of the large field coefficient γ, we use an asymptotic expansion (17) up to order  Table 4. Numerical results for the asymptotic expansion coefficients (27) for all N (ρ 0 is given in Tab. 2). We quote only those digits which remain stable under variations of the matching points and under a variation of the corresponding approximation order by ∆M = 2.
ρ −30 . For the coefficients ζ 1,0 and ζ 4,0 we use an asymptotic expansion up to order ( corresponding to about a hundred terms in the series (20). Note also that we have rescaled the coefficients (27) with ρ 0 (N ) to factor-out redundant N -dependences. Again, this is equivalent to a rescaling with (7) of (8) for the choice A = ρ 0 (N ). The plot in Fig. 6 illustrates that in this normalisation the coefficient γ(N ) smoothly approaches its infinte-N asymptotic value given by in accord with the findings in [28]. Similarly, for large negative ρ, corresponding to purely imaginary fields and assuming that ρ 0 is the only global "scale" of the fixed point solution, naive dimensional analysis suggests that the coefficients ζ 1,0 (N ) scale proportionally to √ ρ 0 (N ), and that ζ 4,0 (N ) should scale proportionally to ρ 2 0 (N ), modulo logarithmic corrections. Our results for these are shown in Tab. 4 and in Fig. 6. The expected behaviour is confirmed for ζ 1,0 (N ). Unlike the coefficient γ(N ), however, it has not yet fully settled to its asymptotic value. The reason for this is that the limit N → ∞ is continuous only for real fields with ρ > ρ c , but not so for purely imaginary fields with ρ < ρ c , see Sec. III B. In fact, the asymptotic solution for large negative ρ at N = ∞ is different from the expansion at any finite N due to the absence of the radial mode fluctuations. We also find that ζ 4,0 (N ) roughly scales as ∝ ρ 0 (N ), rather than ∝ ρ 2 0 (N ). For large-N , it is therefore smaller by a power of N than naively expected. This behaviour indicates that subleading terms are more strongly suppressed than the leading one, for large N , in accord with (32).

J. Radii of convergence
In order to further quanitfy the above we define an empirical radius of convergence for the small-field expansion around ρ = ρ 0 based on the data, see Fig. 7. To that end, we determine a maximal (minimal) field value ρ r (ρ m ) such that the accuracy of the fixed point solution N acc (ρ), defined in (24), remains larger than a pre-set desired accuracy N da for all fields ρ 0 < ρ < ρ r (ρ m < ρ < ρ 0 ). Our results are illustrated in Fig. 7. The blue circles (red stars) represent for each value of N the results for ρ r (upper lines) and ρ m (lower lines), for a polynomial approximation at order M = 30 (M = 20). The colour intensity increases with increasing "desired accuracy" N da = 0, 2, 5, 6. The data points above and below the 0-axis correspond to the radius of convergence above and below the potential minimum, respectively, indicating approximately symmetric convergence properties in these two domains, ρ r − ρ 0 ≈ ρ 0 − ρ m . The common radius of convergence is read off as In order to guide the eye we have connected those data points for M = 30 where the difference between the full result and the polynomial approximation in terms of N acc is of order one; this line is our estimate for the radius of convergence. Broadly speaking, we find that for N > 1. It is intuitively clear that the dark points lie closer to the origin. In fact, the requirement that N acc = 6 is a much more stringent choice for the definition of the radius of convergence. This empirical definition for the radius of convergence for intermediate values of N is compatible with (16). The results for M = 20, marked by red stars, behave in a very similary way. In fact, the results agree very well over a wide range of N indicating that an excellent accuracy is already achieved at M = 20 irrespective of the universality class. Exemplarily, we evaluate the radius of convergence with higher accuracy for the Ising universality class N = 1, and for the infinite-N limit. In the Ising case, we use the Mercer-Roberts test [28] for the numerical expansion coefficients to find the estimate for N = 1, which is in good agreement with our findings in Fig. 7 and the estimate of the radius based on the present study, displayed in Fig. 8. At infinite N , the radius of convergence is related to the unique pair of complex conjugate solutions z p and z * p of the transcendental equation [28] H H −1 (z p ) = 0 , where the function H is given in (38) and H −1 denotes its inverse, H −1 (H(z)) = z. From (46), one then determines the radius of convergence in units of the vacuum expectation value as R A = |z p | × ρ 0 . Quantitatively, which is in good agreement with the numerical estimates achieved for large N , see Fig. 7.
K. Accuracy of global fixed point In Fig. 8 we illustrate exemplarily for N = 1 and N = 100 the accuracy of the fixed point solution to leading order in the derivative expansion, covering the entire field space ρ ∈ [−∞, ∞] of real and purely imaginary fields. We also highlight the regime where we used a small field (i) At infinite N , only Goldstone mode fluctuations contribute and the radial mode has decoupled. We observe that the local expansion A has an overlapping radius of convergence with those of expansions C and D, see Fig. 9. Consequently, the parameters m 2 , γ and ξ are uniquely determined through the couplings of expansion A, see Tab. 1. This provides us with a global and analytical solution u * (ρ) for all fields [28].
(ii) At finite N , the fluctuations of the radial mode contribute. We observe that large and small field expansions no longer have overlapping radii of convergence, illustrated in Fig. 8 through the intermediate (green) regions E and F . The gap in field space is narrower for real (E) than for imaginary fields (F ), and mildly depends on N and the competition between Goldstone and radial mode fluctuations. The fixed point equation, albeit stiff, admits a stable numerical integration in the intermediate regimes E and F . The accuracy in region E is rather homogeneous. In region F (where the integration starts at large negative ρ) the accuracy slowly decreases inwards, with decreasing −ρ. In either case the overall accuracy is solely limited by the accuracy of the adopted Figure 9. The global Wilson-Fisher fixed point u * (ρ) in the infinite N limit, showing that the small field expansion (A), the large real field expansion (C) and the large imaginary field expansion (D) have overlapping radii of convergence. Axes are rescaled as in Fig. 1 for better display. numerical algorithm. 3 (iii) The accuracy of all local expansions decreases with increasing distance from the expansion point, as is evident from the decrease of N acc in the bottom-half of the plots, see Fig. 8. This is a necessary fingerprint of local expansions. We have checked that local expansions agree with the numerical integration in the regions where they overlap. This way, an overall accuracy of 10 − 20 relevant digits is achieved.

L. Accuracy of universal exponents
In Fig. 10, the accuracy and convergence to leading order in the derivative expansion of the four leading scaling exponents are shown as a surface plot for all N and all approximation orders M . All plots show a monotonic increase in the number of significant digits N Xn as the order of the polynomial approximation is increased, M → M + 1. This means that the scaling potential can indeed be improved systematically, and for all N . Interestingly, the accuracy in the result is systematically higher for larger values of N . This result establishes that the numerical strategy used here is applicable, and that the boundary condition (23) is self-consistent. We have also found indications that an improved boundary condition, informed by semi-analytical input for the higher order couplings, improves the convergence even further. Given the quality in the results, we have not attempted to implement more refined strategies here.  Figure 10. Accuracy of scaling exponents, characterised by the number of stable relevant digits N X (26) for the exponents X = ν, ω, ω 1 , ω 2 (from top left to bottom right, respectively) for all universality classes studied here, and as a function of the order of the polynomial approximation M . We note that the accuracy at fixed order in the approximation increases rapidly with increasing N (see text).

M. Wilson-Polchinski flow
We finally discuss a link between optimised and Wilson-Polchinski flows related to singularities of their fixed point solutions in the regime of purely imaginary fields. In [8], it has been conjectured that the optimised RG flows [16,17] studied here are equivalent to the Wilson-Polchinski RG flow [12,51] in the same approximation. In fact, both of these RG flows are related by a Legendre transform. By now, their equivalence has been established both numerically [8,10], by evaluating universal scaling exponents in either version, and analytically [10,76], with the help of the Legendre transform.
Here, we briefly discuss how the singularity for large and purely imaginary fields translates to the corresponding fixed point solution of the Wilson-Polchinski flow. For asympotically large imaginary fields the first derivative of the scaling potential takes the form (20), up to next-to-next-to-leading order corrections, see (20), (22). We note that since the variable ρ ∝ φ a φ a takes all real values on the fixed point solution, the fields φ a exhaust all values on the imaginary field axis. In particular, u * remains finite throughout. We now turn to the Wilson-Polchinski flow for the potential v(ρ). In d = 3 dimensions and for N = 0 it can be written as [8] where ρ = 1 2 φ a φ a /k is the square of the field φ a in units of the RG scale k. At the Wilson-Fisher fixed point v * (ρ), and in the regime of purely imaginary fields and negative ρ, it has been observed that [77] v where the parameter ρ p (N ) > 0 depends on the universality class N . In contrast to (48), the result (50) shows that the integrated Wilson-Polchinski flow exhibits a simple pole at finite fields on the imaginary field axis, with ϕ = |φ a |. Performing a Legendre transform following [10] provides us with relations between the effective potential v(ρ) of the Wilson-Polchinski flow and the potential u(ρ) from the optimised RG. This leaves us with an explicit link between the parameters ρ p in (50) and the coefficient ζ in (48), Hence, although the singularities in the complexified field plane in either setting are qualitatively different, as evidenced in (48), (50), the behaviour of the Wilson-Fisher fixed point for asymptotically large and purely imaginary fields in one setting dictates the location of a pole in the other, both given by the same parameter ζ. In Tab. 5 we confirm the result (52) numerically, comparing our results for ζ from (48) with those derived numerically from (50) in [77]. Besides offering a consistency check for the analytical results of [28], the quantitative agreement also confirms equivalence of renormalisation group flows even for purely imaginary fields.

IV. DISCUSSION
The combined use of analytical and numerical methods is often vital for the understanding of non-perturbative phenomena in quantum and statistical field theory. Here, in combination with the analytical studies put forward in [28], we have performed a detailed numerical analysis of O(N ) symmetric scalar field theories at their Wilson-Fisher fixed point in three dimensions. We have exploited polynomial expansions up to very high order, showing that these deliver reliable and rapidly converging results for scaling exponents when combined with expansions about asymptotically large, or purely imaginary fields, and numerical integration. In this manner we have achieved global fixed point solutions for the quantum effective action with moderate numerical effort, also obtaining high numerical accuracy for scaling exponents to leading order in the derivative expansion, critical couplings, field ratios, and a number of other universal, super-universal, and non-universal features of scaling solutions. Our findings strengthen the view that scaling exponents can efficiently be deduced from data encoded in the small field region of the fixed point effective action [25].
Furthermore, our study showed that large-and small-field expansions often do not yield overlapping radii of convergence owing largely to the fluctuations of the radial mode. This affects the physically most interesting universality classes where the number of fields, parametrised by N , is finite and small. Interestingly, this phenomenon is absent at infinite N where the radial mode decouples and the Goldstone modes dominate instead. The competition between longitudinal and transversal fluctuations appears to be most pronounced for moderate values of N where the Wilson-Fisher fixed point is more strongly coupled. Quantitatively, however, we found that the non-perturbative gap in field space is narrower for real than for purely imaginary fields, and that the intermediate field region is very efficiently covered by straightforward and stable numerical integration. The procedure appears to be sufficient for all technical purposes.
We have also provided a detailed analysis covering the entire field space including the region of purely imaginary fields. In the latter region, qualitative differences between longitudinal and transversal mode fluctuations are particularly pronounced. The strict convexity bound (30) is relaxed for the radial modes at infinite N without causing a singularity, but not so at any finite N where the bound remains firmly in place. Consequently, the N -dependence of the effective potential is discontinuous at infinite N . No such discontinuity is visible in the region where fields are real. As a by-product, we also observed that while fixed point solutions of (5) are free of singularities for any real or purely imaginary field, those of the Polchinski flow (49) unavoidably develop singularities at finite and purely imaginary field. Interestingly, the behaviour of the Wilson-Fisher fixed point solution for large and purely imaginary fields in the former setting dictates the location of the singularity in the latter, for any universality class. Overall, we have thus established a comprehensive global picture covering all the salient features of non-perturbative scaling solutions.
Although we have limited ourselves to the leading order in the derivative expansion, our line of reasoning and methodology can be used beyond the leading order as well as for other systematic expansions. At second and fourth order in the derivative expansion, some focus has been given to optimised choices for the Wilsonian momentum cutoff [5,16,[73][74][75]. It is conceivable that predictions for universal exponents can be further improved by exploiting the structure of equations in the complexified field plane [28], possibly supported by conformal mappings [27] and the close links with the Polchinski flow. Some of the qualitative and quantitative insights of our study will also prove useful for fixed point searches in other theories, e.g. 3d Wess-Zumino models [40,41], or 4d quantum gravity [38]. While the gravitational renormalisation group is more involved than the equations explored here, the interplay between singularities in the complexified field plane and radii of convergence is operative in gravity as well [78].