Operator product expansion coefficients from the nonperturbative functional renormalization group

Using the nonperturbative functional renormalization group (FRG) within the Blaizot-M\'endez-Galain-Wschebor approximation, we compute the operator product expansion (OPE) coefficient $c_{112}$ associated with the operators $\mathcal{O}_1\sim\varphi$ and $\mathcal{O}_2\sim\varphi^2$ in the three-dimensional $\mathrm{O}(N)$ universality class and in the Ising universality class ($N=1$) in dimensions $2 \leq d \leq 4$. When available, exact results and estimates from the conformal bootstrap and Monte-Carlo simulations compare extremely well to our results, while FRG is able to provide values across the whole range of $d$ and $N$ considered.


I. INTRODUCTION
The nonperturbative functional renormalization group (FRG) provides us with a versatile technique to study strongly correlated systems. It has been used in many models of quantum and statistical field theory ranging from statistical physics and condensed matter to highenergy physics and quantum gravity [1][2][3]. Besides the interest in models where perturbative approaches or numerical methods are difficult for various reasons, there is an ongoing effort to characterize and quantify the efficiency and the accuracy of the FRG approach by considering well-known models of statistical physics. It is now proven that the FRG yields very accurate values of the critical exponents associated with the Wilson-Fisher fixed point of O(N ) models [4,5], comparable with the best estimates from field-theoretical perturbative RG [6,7], Monte Carlo simulations [8][9][10][11][12][13] or conformal bootstrap [14][15][16][17]. The FRG also allows the computation of universal quantities defined away from the critical point, such as universal scaling functions [18][19][20][21] or universal amplitude ratios [22], again in remarkable agreement with Monte Carlo simulations when available.
On the other hand, the operator product expansion (OPE) has received little attention in the framework of the FRG until recently [23][24][25][26][27][28][29]. Wilson and Kadanoff suggested independently that in a quantum field theory the product of two operators in the short distance limit is equivalent to an infinite sum of operators multiplied by possibly singular functions when inserted in any correlation function [30][31][32][33]. The validity of the OPE has been proven to all orders in pertubation theory [34] and can be established in full generality in the case of conformal field theories [35]. Indeed, the OPE has been fundamental in the study of conformal field theories in two and higher dimensions [36,37]. In this context, the conformal bootstrap program [37][38][39][40] has lead to a large number of precise results. The OPE has been instrumental as well in studies regarding quantum chromodynamics [41] and condensed matter, where it has been used to derive the thermodynamic properties of quantum gases [42,43].
Despite the fact that both the FRG formalism and the OPE offer non-perturbative approaches to quantum field theory, it is not yet clear to what extent these two aspects can be usefully combined to extract information regarding the non-perturbative regime of a field theory.
From the perspective of perturbation theory, the FRG provides a useful framework that allows one to prove the existence of the OPE perturbatively [23][24][25][26][27]. Moreover, by following the proposal of Cardy relating the OPE coefficients to the second order terms in the expansion of the beta functions around a fixed point [44], the standard perturbative renormalization group has been used to derive certain OPE coefficients within the expansion [45]; we refer to [46] for an FRG perspective on these issues based on a geometric approach to theory space.
In principle, one may reconstruct from the FRG the full operator product and express the latter as an OPE [28,29]. However, this may be rather cumbersome in practice. For a conformally invariant fixed point theory [47], a further possibility explored in [29] consists in extracting the OPE coefficient from three-point functions. It has been shown that within this approach it is possible to calculate the OPE coefficients in the epsilon expansion.
The main quantities of interest in the FRG are the effective action, defined as the Legendre transform of the free energy, and the one-particle irreducible (1PI) vertices. Taking the Wilson-Fisher fixed point of the O(N ) model as an example, we show how the OPE coefficient c 112 associated with the operators O 1 ∼ ϕ and O 2 ∼ ϕ 2 can be deduced from a small number of low-order 1PI vertices. One difficulty in the computation of OPE coefficients is that the latter are determined by the full momentum dependence of the vertices in the critical regime. For this reason, one has to go beyond the derivative expansion in order to accurately determine the OPE coefficients. The latter can be computed in the so-called Blaizot-Méndez-Galain-Wschebor (BMW) approximation that enables the determination of the momentum dependence of the correlation functions [48][49][50]. This approximation scheme has been used in the past to obtain the spectral function of the "Higgs" amplitude mode in the (2 + 1)-dimensional O(N ) model [51] providing an estimate of the Higgs mass that has been confirmed by subsequent numerical simulations of lattice models [52][53][54].
The outline of the paper is as follows. In Section II we recall the relation between the OPE coefficients and the two-and three-point functions in momentum space, focusing on the coefficient c 112 in the d-dimensional O(N ) model. We then show how to relate c 112 to the 1PI vertices. Finally we briefly describe the nonperturbative FRG formalism and the BMW approximation. In Section III the results obtained from a numerical solution of the flow equations are discussed for the three-dimensional O(N ) model and the Ising university class (N = 1) in dimensions 2 ≤ d ≤ 4, and compared with exact values in some particular cases and estimates from conformal bootstrap and Monte Carlo as well as and large-N expansions.

A. Correlation functions in momentum space
We consider a critical, conformally invariant, theory. For fields O a (x) (be them composite or not) with scaling dimensions ∆ a , the two-and three-point correlation functions are given by and

13
(2) where x 12 = |x 1 − x 2 |, etc. Equation (1) assumes a proper normalization of the fields and the coefficient c abc in (2) can be identified with the OPE coefficient [55]. Since in practice we shall work in momentum space, it is convenient to consider the Fourier transformed correlation functions. For the two-point one, where with Γ(x) the gamma function and d the dimension. The is given by a complicated expression but it is sufficient to consider the limit to extract the coefficient c abc [29]. Equation (5) entails that the OPE coefficient c abc can be deduced from the three-point function (5) once the fields have been properly normalized in order to satisfy (3).

B. O(N ) model and Wilson-Fisher fixed point
In the following we consider the O(N ) model in d dimensions defined by the action and regularized by a UV momentum cutoff Λ. ϕ = (ϕ 1 , · · · , ϕ N ) is a N -component field. The model can be tuned to its critical point by varying r 0 . The correlation functions are then scale and conformal invariant in the momentum range |p| is the Ginzburg scale. In the following we shall only be interested in the critical point and the scaling limit |p| p G ; we refer to [56] for an overview of the various regimes of the O (N ) model.
We focus on the operators  (3) in the scaling limit |p| p G . Even though O 2 is not, strictly speaking, a scaling operator, it can be expressed as a linear combination of scaling operators, among which that associated with ∆ 2 . In the scaling limit, corrections due to higher scaling dimension operators are suppressed and we neglect them; a more detailed explanation regarding this point is provided at the end of Section II D. The scaling dimension is related to the anomalous dimension η while where ν is the correlation-length exponent.
To deal with the composite field O 2 , in addition to the linear source J we introduce in the partition function a source h coupled to ϕ 2 , The correlation functions of interest, besides the propa- for the connected correlation function), are the scalar susceptibility and the three-point function Here and in the following, there is no implicit summation over the index i. The computation of G(p) and χ s (p) allows us to determine the normalization constants N 1 and N 2 since at criticality for p → 0. The knowledge of χ(p 1 , p 2 , −p 1 − p 2 ) then yields the OPE coefficient c 112 using (5).

C. Effective action
The effective action is defined as the Legendre transform of the free energy [57]. The source J and the order parameter field φ are related by All correlation functions for h = 0 can be obtained from the one-particle irreducible (1PI) vertices where, assuming the absence of spontaneously broken symmetry, we have set φ = 0. In particular, the propagator is related to the inverse of the two-point vertex computed with a vanishing source h = 0. The other two correlation functions of interest are given by [51] where we have used the fact that Γ (1,1) i vanishes when evaluated for φ = 0. To alleviate the notations we do not write the last argument of the three-point vertices, e.g., Γ We are now in a position to relate the OPE coefficient c 112 to the 1PI vertices at criticality. From Eqs. (13) we obtain the normalization constants Considering (5) in the limit p 2 = 0 and p 1 → 0, we finally deduce Equations (19) to (21) are the basic ingredients to determine the OPE coefficient c 112 in the effective action formalism. In Appendices A and B, we show that they yield the known results for the free theory (u 0 = 0), and in the large-N limit to the leading order and for d < 4, in agreement with the literature [58,59].

D. FRG formalism and BMW approximation
The nonperturbative FRG allows one to compute the effective action beyond standard perturbation theory [1][2][3]. Fluctuations are regularized by the infrared regulator term where the momentum scale k varies from the UV cutoff Λ down to zero. A possible choice for the cutoff function with the function r(y) taken to be for instance r W (y) and r E (y) define respectively the so-called Wetterich and exponential regulators. In either case, α is a constant of order one and Z k a field renormalization factor which varies as k −η at criticality [60]. Thus the regulator suppresses fluctuations with momenta |q| k but leaves unaffected those with |q| k. The partition function is now k dependent. The scale-dependent effective action (28) is defined as a slightly modified Legendre transform which includes the subtraction of ∆S k [φ]. Assuming that for k = Λ the fluctuations are completely frozen by the regulator term, On the other hand, the effective action of the O(N ) model (6) is given by Γ k=0 since R k=0 vanishes. The FRG approach aims at determining Γ k=0 from Γ Λ using Wetterich's equation [61][62][63] The infinite hierarchy of flow equations satisfied by the k-dependent 1PI vertices Γ (n,m) k can be obtained from (30) by taking functional derivatives wrt φ and h. The presence of the source h in addition to the field φ allows one to follow the flow of composite fields, an approach which proved to be useful in tackling a wide range of issues [51,[64][65][66][67][68][69][70][71][72].
In the BMW approximation [48][49][50], one considers the flow equations of the 1PI vertices in a uniform field φ even if one is eventually interested in the vanishing field configuration. These equations are shown diagrammatically in Fig. 1  . Noting then that a vertex with a vanishing momentum can be related to a lower-order vertex, e.g., we obtain a closed set of equations satisfied by Γ The reader may wonder if the simple choice of bare operators made in (7) is suitable to address our objectives. Thus, let us expound on some fundamental properties of composite operators within the FRG formalism. In the FRG framework a composite operator can be defined by differentiating the effective action Γ k wrt an external source. For instance, the running composite operator corresponding to O 2 is defined as φ 2 (k) ≡ − δΓ k δh [φ, h]| h=0 . φ 2 (k) depends on the RG scale k and, if evaluated at the UV cutoff scale, satisfies φ 2 (Λ) = φ 2 . As soon as one lowers the scale k from Λ, the bare term gets renormalized through its coupling to the 1PI vertices. This implies that the flow of the composite operator generates mixing with other operators in the sense that φ 2 (k) = Z 22,k φ 2 + Z 24,k φ 4 + · · · .
It must be noted, however, that scaling operators are not just mere composite operators. A scaling operator is a particular combination of composite operators that diagonalize the flow linearized about the fixed point. As a consequence, in the scaling regime, φ 2 (k) can be expressed as a linear combination of scaling operators. However, by lowering the RG scale k to 0 in such a linear combination, the scaling operators of higher scaling dimension ∆ i are suppressed and only the lowest scaling operator survives. Indeed, numerically solving the flow of e.g. Γ (0,2) k within our approximation scheme, we have been able to check that our procedure reproduces the expected behavior of a scaling operator in the fixed point regime [73].
Let us conclude by noticing that the dependence of the vertices Γ (1,1) k on the background field φ entails an infinite number of 1PI vertices (albeit in a specific momentum configuration). This shows that our ansatz includes non-trivial mixing among field monomials φ 2 , φ 4 and so on.

III. NUMERICAL RESULTS
The flow equations are integrated numerically, see e.g. Refs. [50,51] for details. We work with dimensionless variables,p = p/k andρ = Z k k 2−d ρ. The field dependence of the potential and the vertices is discretized on a finite and evenly spaced gridρ ∈ [0,ρ max ] comprising N ρ points, while the momentum dependence of the vertices is approximated by Chebyschev polynomials of order N p defined on [0,p max ]. The integration of the flow with respect to the RG scale k is done with an adaptive step integration. Convergence of the results with respect to the parameters has been verified; their typical range are N ρ = 40-80,ρ max = 4-8,p max = 4-10 and N p = 20-30 with the precise value depending on d and N .
For each universality class set by d and N and each choice of the cutoff function (25) parameterized by α, the critical point is found by tuning the initial condition of the flow. This enables the computation of G(p), χ s (p) and Γ (2,1) (p, 0) [Eqs. (13) and (31)] at criticality, from which one fits the values of the critical exponents η and ν (or equivalently ∆ 1,2 ) and normalization constants N 1,2 , yielding c 112 through Eq. (21).
A crucial question is that of the regulator dependence. Indeed, while Eq. (30) is exact, any approximation scheme such as BMW introduces a regulator dependence to the results. In order to provide a meaningful prediction for a physical quantity Q(α), a choice of the regulator must be made. The usual rationale is the so-called principle of minimum sensitivity (PMS), according to which the best value of α is that for which the regulator dependence of Q(α) is minimal, i.e., for which ∂ α Q(α) = 0, or failing that for which |∂ α Q(α)| is minimal. However, the PMS for c 112 , shown for the 3d Ising universality (N = 1) class in Fig. 2, does not provide a satisfactory result. Indeed, for a given regulator, c 112 is a monotonous concave function of α, with no extremum or inflection point, varying by about 3% over the range of regulators considered. As a consequence, we choose the regulator that fulfills the PMS for the anomalous dimension η. The value thus obtained for c 112 depends only weakly on the family of regulators considered, with a variation of about 1.15% between the Wetterich and exponential regulators. The regulator dependency is slightly smaller than the 1.2% difference with the conformal bootstrap estimate.
As a side note, we point out a recent proposal for an alternative way to fix the regulator dependence for conformally invariant theories, the principle of maximal conformality (PMC) [74]. Conformal invariance implies a set of (modified) Ward identities associated with scale and special conformal transformations (SCT). While scale invariance is always fulfilled at the fixed point, invariance under SCT is broken within the derivative expansion at high order. PMC suggests to choose the regulator that minimizes the symmetry breaking. While in the present case it is not straightforward to implement the PMC for BMW, because the Ward identities are either trivially fulfilled or involve high-order vertices that cannot be computed using the BMW approximation, its implementation for the derivative expansion shows that the PMS for η and PMC yield very close results, providing a further argument in favor of our regulator choice.
A. Ising university class in dimensions 2 ≤ d ≤ 4 We first consider the OPE coefficients of the Ising university class (N = 1) for dimensions d between the lower and upper critical dimensions d = 2 and 4, for which the results are shown in Fig. 3 and Table I. The FRG results can be compared to the exact values in d = 2 and 4, the conformal bootstrap [14] and Monte Carlo [75] estimates in d = 3 and the expansion up to fourth order [76][77][78][79], where ζ(3) is Apéry's constant and A 4 = −0.158947 . . . is only known numerically [79].
Owing to its asymptotic nature, the expansion does not converge; indeed, for d = 3, its relative error to all estimates (FRG, conformal bootstrap, Monte Carlo) increases from 4% at order O( 2 ) to 6% at order O( 3 ) and 16% at order O( 4 ). In order to make sense of the results, a resummation procedure must be carried out, for instance by approximating c 112 ( ) by a Padé approximant of order (n, m), i.e. a rational fraction of with numerator and denominator of degree n and m, respectively. Following Ref. [79], we pick an approximant of order (3,2), whose coefficients are uniquely determined by imposing the O( 4 ) expansion around = 0 and the exact value c 112 ( = 2) = 1/2. This gives c 112 = 1.0507 for d = 3, in very good agreement with the conformal bootstrap, to be compared to the 0.8% error when the approximant of order (3, 1) is used and the exact result for d = 2 is not imposed. Let us mention that in this case different choices of Padé approximants lead to somewhat different results, which may be used to estimate an average value and its uncertainty: c 112 = 1.048(31), see Table I Compared to the best results (exact in d = 2 and 4, conformal bootstrap and Monte Carlo in 2 < d < 4), the FRG always has an error smaller than 3%, with less than 2% error in d = 3. By contrast with the resummed expansion, which requires additional input in the form of the 2d data to provide accurate results, the FRG and  conformal bootstrap are able to interpolate smoothly between dimensions d = 2 and 4. As the dimension is increased, c 112 increases monotonously, with an almost linear behavior between d = 2 and d = 3, which might partly explain the remarkable agreement of the expansion resummation (when supplemented with the exact 2d result) with conformal bootstrap, Monte Carlo and FRG. However, if the resummation is not supplemented with the 2d result then the -expansion estimate of c 112 is rather poor close to and at d = 2.
In d = 4, the FRG within the BMW approximation scheme gives the exact analytic value of c 112 . The small difference (∼ 0.1%) between the numerical result and the exact value seen in Table I [78]. The black symbols correspond to estimates from Monte Carlo (square) [75,83] and conformal bootstrap (diamond) [14,17,82].
serves as an estimate of this numerical error: in lower dimensions, it is much smaller than the difference to the best estimates.

B. Three-dimensional O(N ) model
We now focus on the three-dimensional O(N ) model. Given that the large-N result is [58,59,81] we consider rather than c 112 the rescaled OPE coefficient √ N c 112 that has a well-defined large-N limit. FRG results and estimates from the expansion [78], conformal bootstrap [14,17,82] and Monte Carlo [75,83] are shown in Fig. 4 and Table II. For N > 1, c 112 is only known up to order O( 3 ) and we resum the expansion using a (2, 1) Padé approximant.

IV. CONCLUSION
We have shown how to extract the OPE coefficients of a conformal theory within the framework of FRG, by determining three-point vertices in specific momentum configurations. We have used our approach to determine the c 112 coefficient, corresponding to the simplest possible OPE coefficient, in the O(N ) universality class for various d and N . This provides the first non-perturbative determination of the OPE coefficients based on field theory, aside from the lattice computations in [75,83] and the conceptually different conformal bootstrap.
While the accuracy of FRG can be sometimes difficult to gauge in the absence of a small expansion parameter, the fact that the results compare extremely well with the values, when available, obtained from Monte Carlo and conformal bootstrap increases confidence in the validity of the method. It is a testament to the versatility of FRG that, in this specific case, tuning such parameters as d or N demands relatively little effort as they only enter the flow equations through their explicit values.
Lastly, we note that the OPE can be used in settings very different from the critical O(N ) models investigated in this work, for instance in theories away from a fixed point or at a non-equilibrium fixed point. In these cases many methods holding for equilibrium critical theories are not available. Our work suggests that the FRG may constitute the right framework to tackle these issues thanks to its aforementioned versatility and to the fact that the FRG equations can be solved without invoking further requirements, such as conformal symmetry.

ACKNOWLEDGMENTS
The authors thank Gonzalo De Polsi, Matthieu Tissier and Nicolás Wschebor for stimulating discussions and are extremely grateful towards Johan Henriksson for pointing out Refs. [59,79,80]. C. P. and N. D. wish to thank the organizers of the meeting FRGIM 2019 where some of the ideas developed in this work were discussed. C. P. thanks Hidenori Sonoda for discussions and collaboration on closely related projects and the Laboratoire de Physique Théorique de la Matière Condenséee in Paris for hospitality.
The expectation value of the field is given by and the effective action is simply We thus obtain where The last expression for Γ (0,2) (p) is obtained for p → 0 and d < 4. From (A5) and (19) we deduce ∆ 1 = d/2 − 1 (in agreement with the vanishing anomalous dimension in the free theory) and On the other hand, comparing (A6) with (20) yields ∆ 2 = d − 2 (in agreement with ν = 1/2) and (A10) From (21) and (A7) we finally obtain (22).
Appendix B: c112 in the large-N limit Following the Appendix of Ref. [51], we introduce the field ρ = ϕ 2 and a Lagrange multiplier λ to write the partition function Z[h] ≡ Z[J = 0, h] of the O(N ) model as Then we split the field ϕ into a σ field and an (N − 1)component field π. Integrating over the π field, we obtain the action S[σ, λ, h] =ˆx − 3N 2u 0 (2h + iλ − r 0 ) 2 + 1 2 (∇σ) 2 + iλσ 2 + N − 1 2 Tr ln g −1 [λ], (B2) where is the inverse propagator of the field π i in the fluctuating λ field. In the limit N → ∞, the action becomes proportional to N (if one rescales the σ field, σ → √ N σ); the saddle point approximation becomes exact for the partition function Z[h] and the Legendre transform of the free energy coincides with the action S [84]. This implies that the effective action is simply equal to S[σ, λ, h]: which is the starting point to compute the vertices Γ (n,m) in the large-N limit.