Universal location of Yang-Lee edge singularity in classic O(N) universality classes

Employing the functional renormalization group approach at next-to-leading order of the derivative expansion, we refine our earlier findings for the location of the Yang-Lee edge singularity in classic O(N) universality classes. For the universality classes of interest to QCD, in three dimensions, we found $|z_c|/R_\chi^{1/\gamma} = 1.612(9),\ 1.597(3)$ for $N=2$, $4$ correspondingly. We also established $|z_c| = 2.04(8),\ 1.69(3)$ for $N=2$, $4$ albeit with greater systematic error.

Due to the divergence of the correlation length near a second-order phase transition, the dynamics of the system becomes independent of the microscopic details and only reflects the grand properties -the dimensionality and global symmetries. This allows one to collect systems of varied microscopic origin into a limited number of universality classes. By studying one member of such a class, the emergent universal behavior allows one to establish the properties of many different systems regardless of their microscopic complexity.
Thus, it is not surprising that for the most ubiquitous classic universality classes of O(N ) systems many universal properties (such as critical exponents, critical universal amplitudes, critical equations of state) are known with extreme precision, see e.g. Refs. [1][2][3][4][5] and references therein. One notable exception is the universal location of the Yang-Lee edge singularity, which was only recently determined in Refs. [6,7]. In this paper, we continue to refine these results.
Lee and Yang demonstrated an intimate connection between the analytical structure of the equation of state and the phase structure [8,9]. Specifically, in the symmetric phase, An interesting property of the YLE critical point is that it is characterized by a ϕ 3 theory and consequently its upper critical dimension is 6 [14]. Therefore, the conventional ε expansion near four dimensions applied to study the Wilson-Fisher critical point of the underlying universality class has only limited predictive power for locating the YLE singularity.
We come back to this in more detail in Sect. III B. with the realistic critical exponents and the location of the singularity obtained |z c | from Ref. [7] complemented by the value of R χ = 1.72 from Ref. [2]. For the mean-field equation of state and for the large N limit in d = 3, the spinodals are located on the real h axis due to integer values of 2βδ.
Only positive values for real values of h are shown.
The numerical calculations in this paper are performed using the Functional Renormalization Group (FRG) approach, see Ref. [15] for a review. We extend the results of our previous work (see Refs. [6] and [7]) significantly. First, we improve the truncation scheme going to the (truncated) first order derivative expansion and including the dependence of the wave function renormalizations on the field for N > 1. In our original study [6], the calculations were performed in the so-called LPA' approximation which assumes a field-independent wave function renormalization, while in Ref. [7] we only investigated the Ising universality class, N = 1. Second, we accounted for the residual dependence on the regulator by performing a minimal sensitivity analysis [16] which was motivated by minimizing the sensitivity to nonphysical parameters in conventional perturbation theory with different renormalization schemes [17].
The paper is organized as follows. We start by defining a required set of universal quantities, functions and non-universal metric factors in Sect. II. We then review analytical results for the location of the YLE singularity in Sect. III: the large N limit and for the number of spatial dimensions close to 4. For the number of components N ̸ = 1, we discuss the behavior of the singularity near two dimensions. In Section IV, we turn to FRG calculations where we extract the location of the singularity for various N in three spatial dimensions.
We end with conclusions in Sect. V.

II. SCALING EQUATION, CRITICAL AMPLITUDES AND EXPONENTS
Consider a system near a critical point with two relevant parameters t and h introduced in a such a way as to detune the system from criticality which occurs at t = h = 0. We will refer to t as the temperature. Its defining property is that non-zero values of t do not explicitly break any symmetries of the system. However, a non-zero t may lead to a spontaneous symmetry breaking either for positive or negative t. Conventionally we assign positive values of t to when the spontaneous symmetry breaking is not possible -in other words, t > 0 defines the symmetric phase of the system. In contrast to t, non-zero values of h, to which we will refer to as the external magnetic field, break the symmetry explicitly. We quantify the system's response to t and h by measuring the order parameter, to which we also will refer to as magnetization M .
The renormalization group analysis (see e.g. Ref. [18]) demonstrates that the equations of state describing the dependence of the magnetization on the parameters t and h has a homogeneous form and can be written as where β and δ are universal critical exponents, and f (x) is a universal scaling function.
In general, the parameters t and h are related to the physical parameters of the system through two non-universal proportionality coefficients, also called metric factors. The metric factors are usually chosen in a such a way as to satisfy two normalization conditions for the function f : f (−1) = 0.
The above form of the scaling equation of state was suggested by Widom [19]. One of its advantages is that it can be straightforwardly derived using the ε expansion. Its disadvantage is that it leads to an implicit dependence of M on t and h. An alternative form solves this issue. Here we have introduced the so-called gap critical exponent, ∆ = βδ. The function f G is a function of one variable; it encodes most of the critical statics. It has to satisfy the normalization conditions lim to be consistent with the Widom scaling. As we alluded to before, the set of the normalization conditions requires the redefinition of the non-universal parameters t and h. Generically near a critical point we have The normalization conditions require us to define t and h in a such a way as to absorb the prefactors B and B c : M = (−t) β , h = 0 and t < 0 .
In Sect. IV, we formulate the FRG approach to locating the YLE singularity. As it is our primary objective, our truncation method is optimized to perform simulations in the symmetric phase. Calculations in the broken phase are possible in a different truncation scheme; however, we want to extract all required quantities within one scheme to avoid introducing systematic errors by mixing different truncations in the simulations. We thus strive to avoid the broken phase. This motivates us to introduce another universal quantity where the universal ratio R χ is defined by the limit and γ is the critical exponent connected to δ and β through the scaling relation: The asymptotic behaviour of the function f G (z) at large argument has a simple physical origin: the magnetization has to be a linear function of h in the symmetric phase t > 0.
From Eq. (4) follows that the scaling function f G (z) therefore has to go like z −(∆−β) = z −γ leading to the identity Eq. (13). Note that working in the symmetric phase allows us to directly extract the critical exponent γ through the scaling for the magnetic susceptibility Using this expression it is straightforward to show that (15) and that the introduced ζ is independent of the amplitude B: thus explicitly demonstrating that in order to extract the location of the YLE singularity in ζ we do not need to perform simulations in the broken phase.
We stress that ζ and z are related through a universal number R χ and universal critical exponent γ. On one hand, R χ can be computed in the FRG approach 1 but would require probing the broken phase and thus switching the FRG truncation scheme used in this paper.
The associated systematic error is difficult to assess. On the other hand, for applications to lattice QCD, |z c | is often considered. We thus will provide a separate set of results for |z c | using R χ obtained in Ref. [20].
Finally, for the purpose of the next sections, we also introduce the anomalous dimension critical exponent η. It describes the power law dependence of the static correlation function on the distance at the critical point. In d spatial dimensions: The anomalous dimension satisfies the following scaling relation [18] 2 − η = d δ − 1 δ + 1 .
1 Precision calculations were performed in a state of the art study in Ref. [20].

III. ANALYTICAL RESULTS FOR LOCATION OF YANG-LEE EDGE SINGU-LARITY
We remind the reader that the Yang-Lee edge singularities are branch points of the function f G (z) in the symmetric phase t > 0. They can be determined by finding zeros of the inverse magnetic field susceptibility. Most generally, the Lee-Yang theorem [8,9] implies that they have to be located on the imaginary h axis. Thus the argument of the singularity z c (and its complex conjugate) is fully determined by the critical exponents of the underlying O(N ) universality class: Note that the argument of ζ c coincides with that of z c .
As far as the absolute value of the location of the YLE singularity |z c | is concerned, O(N ) universality classes with N > 1 do not enjoy many (neither exact nor approximate) analytical results. Notable exceptions are N → ∞ limit and the theory near four dimensions. We detail both analytical results below.
The large N scaling equation of state can be readily computed; for a review, see Ref. [21], where the Widom scaling relation and critical exponents were derived: and The function f defines the value of γ = 2/(d − 2) at asymptotically large values of x and R χ = 1 following Eq. (12). To find the position of the YLE singularity, it is convenient to consider the inverse magnetic field susceptibility χ −1 = ∂h ∂M t . Its zeroes define the values of x at the YLE, x c . In terms of the function f (x), we have This leads to Now we can proceed with finding z c . For this we express z c in terms of x and f : leading to At large N , R χ = 1, thus |ζ c | = |z c |. To compare to the result of the next section, we perform the expansion near four dimensions d = 4 − ε: B. The ε-expansion Although N → ∞ limit provides an analytic result for any d in the range 2 < d < 4, it is not well suited to describe phenomenologically relevant universality classes. Specifically for finite temperature QCD we are interested in the Ising universality class N = 1 and the Heisenberg model with N = 2 (due to lattice discretization artifacts, see e.g. Ref. [22]) and N = 4. There is another analytic limit in which one can perform the calculation -near the upper critical dimension d = 4 − ε. As we document below, as far as the position of the YLE singularity is concerned, the utility of this approach is somewhat restricted and it cannot be systematically improved to yield a reliable result in three dimensions. However, it provides some useful information on the location of the YLE singularity near four dimensions. It is also not limited to N → ∞. Moreover, it serves an estimate of the value of N at which one can apply a large N approximation for the purpose of locating the YLE singularity.
In the conventional ε expansion, see e.g. Ref. [18], near the upper critical dimension d = 4 − ε, the critical exponents are The same method yields the universal amplitude ratio, see Ref. [23], Note that, in the large N limit, this is consistent with the previous subsection, R χ = 1.
To the linear order in ε [24], the scaling function f has the following form We now turn to the location of the YLE singularity. At the leading order ε-expansion, we obtain x (0) c = −3 for the solution of Eq. (22). For our purpose, the exact expression forx(ε) is of no importance as will be demonstrated below. Moreover the first correctionx(ε) to the leading order x (0) c is already non-perturbative (See Appendix A for details). To the first order in ε, Herex(ε) is the leading correction to ε = 0 value of f c . It seems that the presence of the non-perturbative contribution would prevent us from finding the ε order correction to the location of the YLE singularity. This is, however, not the case. Indeed, after expressing z in we obtain where the terms proportional tox(ε) cancel. At higher orders of the ε expansion, this cancellation does not happen. This prevents us from extracting corrections beyond the linear order. For the special case of N = 1, Eq. (35) was previously derived in Ref. [7].
For the absolute value of z, we get Note that, the slope of the function |z c (ε)| is negative for N < 1 + 27 ln 3 ln 2 − 1 ≈ 16.8. It demonstrates that, at least as long the slope of the d dependence is concerned, to reproduce the result of N → ∞ limit one has to consider rather large values of N ≫ 17. Note that the ε-expansion at any given ε leads to a monotonic dependence of the location on N . As we demonstrate in Sect. IV I, for the physical point d = 3 or ε = 1, this dependence is non monotonic.
Using Eq. (31) for the universal ratio R χ leads to In N → ∞ limit, this results reproduces the leading order expansion near four dimension of Sect. III A in d = 4 − ε, see Eq. (27).

C. Behaviour near/at d = 2
For the Ising universality class, N = 1, d = 2 and d → 1 + are analytically treatable. We refer the reader to Refs. [25,26] for d = 2 and to Ref. [7] for d → 1 + . In the latter reference, the two-dimensional results were also presented in the same normalization scheme as in this paper.
The case of d = 2, N ̸ = 1 deserves a special attention.
• N > 2: The perturbative analysis of the non-linear sigma model concludes that the theory near its lower critical dimension d = 2 has a stable ultra-violet fixed point for N > 2 with N = 2 corresponding to Berezinskii-Kosterlitz-Thouless phase transition [28][29][30]. Complemented by the scaling relations this analysis also reveals the full set of critical exponents near two dimensions d = 2 +ε [31,32] η =ε Using Eq. (18), we can establish that to the leading order in ε where N = −2 cos 2π v and 1 < v < 2, smoothly traverses through N = 1 (which 2 See also FRG studies at fractional N and d = 2 in Refs. [34][35][36]. In the presence of ∆ S k [ϕ], the partition function reads (42) and thus becomes scale-dependent. Usually, the following choice is considered In order the match the symmetry of the system, the mass-like regulator function R k (x, y) is chosen to be invariant under rotations (including the internal space rotations) and translations, i.e. R k (x, y) = R k (x − y). Furthermore, in order to suppress modes with p ≲ k while leaving modes with p ≳ k intact, the following must hold for the Fourier transform of the regulator: R k (p) adds a mass of order k 2 to the low-energy modes, thereby suppressing their contributions to the path integral.
The effective action, Γ k [ϕ] is obtained via a modified Legendre transform The functional Γ k [ϕ] satisfies the Wetterich equation [39,44,45],  [65] has been shown to have a finite radius of convergence and a systematic error estimate is possible [63,66]. We use the next-to-leading order of this expansion, i.e. first order in momentum-squared, in this work. This is elaborated in the next sections.

B. First order derivative expansion
In order to solve the flow equation (47) numerically we will consider a constant field configuration around small momentum. The latter is only required to extract equations for the wave function renormalization. Specifically, we use the following anzatz for the scale-dependent effective action where ρ = 1 2 ϕ a ϕ a . The above expression contains all possible terms up to ∂ µ ∂ µ in the derivative expansion. This approximation is appropriate for describing long-wavelength excitations in the critical region.
Consider small deviations around the homogeneous field background, chosen to be non zero for the i = 1 field component then the wave-function renormalization for transverse and radial modes are C. Flow equations for the next-to-leading order in the derivative expansion The flow equations for the potential can be obtained by substituting the truncation Eq. (48) to Eq. (47) evaluated for a constant field configuration where For the integral measure, here and below we use Introducing the tilde differential operator (see Ref. [15]) allows us to write the equations in a succinct diagrammatic manner, and In the diagrams, the internal lines represent the scale-dependent Green functions for the transverse (solid line) and radial (dashed line) modes. The vertices describing interaction relevant for the above flow equations are Using these vertices and applying the tilde derivative on the resulting expressions one can easily derive the flow equations for Γ (2),∥ k (p) and Γ (2),⊥ k (p): and Finally by taking the derivative of ∂ t Γ (2),∥ k (p) and ∂ t Γ (2),⊥ k (p) with respect to p 2 and evaluating at zero p 2 we arrive to the flow equations for the wave-function renormalizations: where we introduced short-hand notations for as well as to simplify the expression we denoted G ′ = ∂G ∂q 2 and G ′′ = ∂ 2 G (∂q 2 ) 2 . In order to obtain Eqs. (66) and (67) we applied the identity and used the expansion In this paper, we consider the so-called strict derivative expansion. The logic behind this approximation is straightforward, see Ref. [63]. Below we rephrase it for the truncation of interest. At the ∂ 2 -order, the momentum-dependent contributions of order q 4 to Γ (4) (p, q, q) are neglected; this justifies neglecting similar terms originating from the square of three-field vertex, (Γ (3) (p, q)) 2 . Dropping the corresponding terms we end up with where γ 0 i = γ i (q = 0). The final form of equations used for the flows of the expansion functions U k , Z ⊥ k , and Z || k in (52) and (74)-(75), while written above in terms of the fields ϕ, are re-expressed in terms of ρ when probing the Wilson-Fisher point. In this form, the regularity of the flows at ρ = 0 becomes apparent. With that said, the expressions in terms of ϕ are also necessary to probe the Yang-Lee edge singularity, as will be discussed.

D. Regulator and wave function renormalization
There are many different choices for the function R k (p). In this work, we consider Litim Here a is a parameter to be varied and optimized under the principle of minimal sensitivity [16], we come back to it in Sect. IV F. We note that we included the radial wave function renormalization at a given field background ϕ 0 , Z ∥ k = Z ∥ k (ϕ = ϕ 0 ), in the regulator. In this, we deviate from the conventional way when Z ⊥ (k) is introduced in the regulator. Our choice is shaped by the problem we are solving: at the YLE, it is expected that the transverse degrees of freedom decouple (the YLE is at a finite imaginary value of the magnetic field, see Fig. 1) while the radial mode is massless and thus dominant. It is convenient to explicitly separate Z ∥ k from the field-dependent Z ∥ k (ϕ). Moreover, we also normalize Z ⊥ k (ϕ) by the same factor, that is To simplify the equation, it is also convenient to introduce the "renormalized" field This enables us to rewrite the Green functions in the following form and At the fixed points, the anomalous dimension is related to the wave-function normalization The anomalous dimension for the transverse component can be defined analogously, but it is not required for our needs.
Note that and by analogy to Eq. (77), it is convenient to introduce

E. Truncation: notation and methodology
Equation (47) is a differential (in k) and functional-differential (in the field space) equation.
Its solution is not known. Without introducing a truncation scheme, this equation cannot be treated numerically.
We thus perform a Taylor series expansion of the functions z to perform an expansion in terms of the "renormalized" field ρ r = ϕ 2 r /2, see Eq. (78). To simplify the notation, we omit the subscript and imply ρ r → ρ unless indicated explicitly.
We have In order to find the location of the YLE singularity, we also will perform the expansion in ϕ, that is We also omitted the subscript k, but emphasize that all parameters are running here. In this paper, we require the truncation orders to satisfy the requirement (N U , N Z ∥ , N Y ) = 2(N U , N Z ∥ , N Y ). This guarantees that the truncations are consistent and one can perform switching between the variables without losing information about the corresponding function.
The FRG flow equations for the coefficients can be readily derived starting from Eqs. (52), (74) and (75). For renormalized quantities, we get a set of equationṡ where we introduced the polarization index i = (∥, ⊥). The coefficients z n ⊥ are related to y (n) through z n ⊥ = z n ∥ − ny (n−1) − ρ 0 y (n) . In order to find the fixed point solutions, we also need to determine the flow equations for the dimensionless renormalized coefficients. They can be obtained by computing the expansion coefficient of dimensionless quantities, e.g. U/k d and z ∥ , as functions of dimensionless We haveu These equations can be rewritten in terms of the expansion coefficients of the series of ϕ without many modifications. It amounts to replacing z i n , u n and ρ 0 with z i n , u n and ϕ 0 and η with η/2 in the right hand sides of Eqs. The conventional expansion scheme, which is used in the majority of applications, follows the physical point defined by the minimum of the effective potential. Since the YLE singularity is located at the (imaginary) magnetic field h c where the radial mass at vanishes, finding it requires numerically expensive fine-tuning with this scheme. The advantage of our expansion scheme is that we can directly follow the flow of the edge singularity. Since the magnetic field enters as a linear term in the effective action, it is not renormalized in the flow equation (47). We can therefore simply read-off h c as the magnetic field which turns our expansion point ϕ 0 into the physical point in the IR. The Lee-Yang theorem guarantees that we only have to resolve the effective action for purely real or purely imaginary fields in the symmetric phase, so full information in the complex plane is not required.
A downside of our expansion scheme is that numerical computations in the broken phase are more challenging. The reason is that in this case at any finite k our expansion point ϕ 0 lies on the real axis below the physical point. The FRG flow flattens the potential in this region since it has to be convex in the deep IR. This convexity-restoring flow is driven by near-singular propagators and therefore numerically challenging to resolve, see, e.g., Refs. [68,69]. However, as we have argued in Sect. II, we do not need to compute in the broken phase. The critical exponent η fully defines δ through the scaling relation Eq. (18). Moreover, calculating the stability matrix at the fixed point solution allows one to find the critical exponent ν and thus the gap critical exponent through the relation ∆ = ν 2 (d + 2 − η). This motivates our strategy in defining the parameter a as the extremum of the function ∆(a) where a enters to the regulator through Eq. (76). This fixes the regulator parameter a = a ∆ that we use for calculating the metric factors (B c and C + ), critical exponents (δ and γ), and the location of the YLE singularity. By choosing the extremum as a function of a, it is guaranteed that among the family of regulators defined in Eq. (76), we use the one where the regulator dependence, and hence the systematic error, is minimal at least for ∆. An alternative approach would be to use minimal sensitivity analysis for all universal quantities, this, however, is not feasible for the location of the YLE singularity as it is defined through the (non-universal) metric factors and the critical exponents.
We show the dependence of the gap critical exponent and the anomalous dimension on the regulator parameter a in Fig. 3 for a few values of N . These calculations were performed using the strict derivative expansion and the truncation scheme (N U , N Z ∥ , N Y ) = (6, 3, 2).
As can be seen in the figure, for N of phenomenological interest, the location of the extrema of ∆(a) and η(a) are fairly close at this truncation.
We obtained reasonable values of the anomalous dimension critical exponent displayed in Fig. 4. The non-monotonic dependence of η on the number of components N is expected from the ε-expansion.
The minimal sensitivity analysis for different quantities does not necessarily result in the same regulator parameter a. In order to estimate the corresponding systematic uncertainty we also performed the analysis for the value of the magnetic field at the YLE singularity, see Table I. There is one subtle point related to how one approaches this fixed point. We remind the reader that at the Wilson-Fisher fixed point, the FRG equations are k-independent for properly scaled variables as we introduced in Sec. IV E. The k independence implies that the fixed point can be reached at a finite value of k. To the contrary, when we start from the general equations for the multi-component field theory, one cannot expect k-independence for the YLE fixed point, as the complete separation of the transverse degrees of freedom is only possible at asymptotically small k. Thus, when finding the algebraic equation for the YLE fixed point, one additionally has to take the limit k → 0 for the terms that involve transverse degrees of freedom. In effect, this amounts to taking the limit of the dimensionless renormalized Goldstone mass to infinity. By computing this limit, we were able to show explicitly that the fixed point equations of our theory reduced to those of the single component theory, that is, the N -dependent terms originating from Eqs. (66) and (52)   As was demonstrated in Ref. [6], this approach is parametrically slow. arbitrary N , coincides with that performed in Ref. [7].

I. Location of Yang-Lee edge singularity
To locate the YLE singularity, we compute in the symmetric phase t > 0 and at vanishing renormalized mass of the radial excitations. Practically we use small but nonzero values m R = 0+. As we are interested in extracting the universal location, we have to consider rather small positive t to minimize non-universal contributions. Similarly to Ref. [7] we perform the switch of the parametrization of the solution from ρ to ϕ at some small but non-zero negative value of ρ s 0 . We checked that our results are in sensitive to the variation of the value of ρ s 0 . Solving the flow numerically yields the determination of h YLE = U ′ (ϕ 0 ) in the infrared.
Performing the mapping to ζ c (see Ref. [7] for details) we obtained the results presented in Table II and illustrated in Fig. 5. The error was computed as follows. First, to estimate the error due to the truncation of the field dependence, ∆ tr , we compare the results in the (5, 2, 1) and (5, 3, 2) truncation schemes to the ones obtained for (6,3,2). We use the result of the highest truncation (6,3,2) for central points and maximal of the absolute values of the two differences |ζ c | (6,3,2) − |ζ c | (5,3,2) and |ζ c | (6,3,2) − |ζ c | (5,2,1) for the error estimate due to truncation. Second, to evaluate the uncertainty associated with the regulator dependence ∆ reg , we perform calculations of |ζ c | at two values of a: a ∆ determined by the minimal sensitivity analysis applied to the critical exponent ∆ at the Wilson-Fisher fixed point (see Sect. IV F) and a h determined by the dependence of the magnetic field at the YLE singularity h YLE . The difference between the corresponding values of |ζ c | determines ∆ reg .
The numerical values for the regulators determined by both schemes are listed in Table I. Both our errors are measures for the convergence of our truncation within the next-to-leading order of the derivative expansion. A meaningful estimate for the systematic error of the derivative expansion itself requires us to go to next-to-next-to-leading order [63,66] 3 . While this is required for a precision calculation, this is beyond the scope of the present work. The next step is to perform the transformation from |ζ c | to |z c |. This step requires the determination of R χ . Our current set up does not allow us to compute this quantity as it necessitates solving the FRG equation in the broken phase which cannot be done using an expansion point defined by equation ∂ t m R = 0 at finite k. Fortunately, high precision calculations of R χ were performed recently in Ref. [20] for N = 2, 3, 4 and 5. We will use   these results together with the value of R χ computed for N = 1 in Ref. [70]. Reference [70] does not provide systematic uncertainty on the value of R χ ; we estimated it by comparing to earlier calculations of R χ in the LPA' FRG of Ref. [46]. We list the results in Table III, where to find R 1/γ χ the value of γ was taken from Ref. [63]. With this we can perform the transformation to |z c |. The result is presented in Table IV. Within the systematic uncertainty the values are consistent to our previous calculations in LPA', see Ref. [6].
Tables II and IV constitute the main results of this paper.

V. CONCLUSIONS
Using the Functional Renormalization Group, we extended our previous results, see Refs.
[6] and [7], for the universal location of the Yang-Lee edge singularity in the most two-dimensional Ising model is taken from Ref. [26], see also Ref. [7] for reparametrization to |z c |.
We expect that, at d = 2, going to fractional values of N would fill in the gap from N = 1 to N = 2.
The line corresponding to N = 1 reaches maximum at some value of d in the range 1 < d < 2 and then drops to 1 at the lower critical dimension d = 1.
important three-dimensional classic O(N ) symmetric universality classes to the (truncated) next-to-leading order in the derivative expansion. Furthermore, we used the prescription of the principle of minimal sensitivity to minimize the regulator dependence of our results.
Our method is best suited for investigating the symmetric phase, and thus our main results are that for |ζ c | = |zc| R 1/γ χ . We used the high precision FRG calculations of R χ and γ from Refs. [20,63] for N =2, 3, 4, and 5 in order to find the location |z c |. For N = 1, R χ was obtained from Refs. [46,70]. See Tables II and IV for  Complemented by recent ideas and methods of Refs. [71][72][73][74][75][76][77][78] our findings might help to establish the existence and potentially the location of the QCD critical point. VI.

ACKNOWLEDGEMENT
We are grateful to B. Friman, T. Schaefer, and, especially, S. Mukherjee for stimulating discussions leading to this work. We thank N. Wschebor for an illuminating discussion on strict derivative expansion. We acknowledge the computing resources provided on Henry2, We see that despite the loop counting, both contributions are of order ϵ (that is the same as the one loop term)! Higher order loop terms also contaminate ϵ 1 order and similar terms involving nested logarithms. This necessitates all loop order resummation and thus brings us to the main conclusion of this appendix that the corrections to the location of YLE singularity are not perturbative in ϵ. Note that this does not prevent us from extracting the leading order correction to z c . This correction is of order of ϵ; higher order contributions suffer from the non-perturbative contribution.