The 2PI effective theory at next-to-leading order using the functional renormalization group

We consider a symmetric scalar theory with quartic coupling in 4-dimensions. We show that the 4 loop 2PI calculation can be done using a renormalization group method. The calculation involves one bare coupling constant which is introduced at the level of the Lagrangian and is therefore conceptually simpler than a standard 2PI calculation, which requires multiple counterterms. We explain how our method can be used to do the corresponding calculation at the 4PI level, which can not be done using any known method by introducing counterterms.


I. INTRODUCTION
There are many interesting systems for which there is no small expansion parameter that could be used to implement a perturbative approach, and for this reason much work has been done in recent years on the development non-perturbative methods. Lattice calculations are valuable in situations where the underlying microscopic theory is known, but issues with the continuum and finite volume limit arise. Various forms of reorganized/improved hard-thermal-loop resummations have been formulated and applied to the calculation of thermodynamic quantities [1][2][3][4][5][6][7][8][9]. Schwinger-Dyson equations are another popular and familar approach to non-perturbative problems (for a classic introduction see [10], a more recent review can be found in [11]). One significant issue with the Schwinger-Dyson approach is that the hierarchy of coupled equations must be truncated by introducing some external prescription. Various methods to construct a truncation that preserves gauge invariance have been proposed [12][13][14].
Another powerful technique is the renormalization group [15], which is traditionally used to study systems where scale-dependent behaviour is important. Its functional formulation can be cast into the form of an exact flow equation for the scale-dependent effective action.
A hierarchy of coupled flow equations for the n-point functions of the theory can be obtained from the action flow equation, but this hierarchy must again be truncated using some prescription [16,17]. Some useful reviews include [18][19][20][21][22][23][24].
The n-particle-irreducible effective action is another method to include non-perturbative effects. In the context of non-relativistic statistical mechanics, the original formalism can be found in Refs. [25][26][27]. In its modern form, the method involves writing the action as a functional of dressed vertex functions, which are determined self consistently using the variational principle [28,29]. The technique has been used to study the thermodynamics of quantum fields [30][31][32], transport coefficients [33][34][35][36], and non-equilibrium quantum dynamics [37][38][39][40][41][42][43][44]. An advantage of the nPI method is that it provides a systemmatic expansion for which the truncation occurs at the level of the action. One major disadvantages is a violation of gauge invariance [45,46]. A method to mimimize gauge dependence has been proposed [47], and some issues with applying the method are discussed in [48][49][50]. Another significant difficulty with the nPI formalism is renormalization. The renormalization of the symmetric 2PI effective theory using a counterterm approach was developed by a number of authors over a period of several years [51][52][53][54]. The renormalization of non-symmetric theories is more subtle, but important for the study of phase transitions and Bose Einstein condensation. Phase transitions are studied in a scalar φ 4 theory in [55][56][57][58][59], and in the SU (N ) Higgs model in 3 dimensions (where vertex divergences are absent) in [60]. Bose Einstein condensation has been studied in Refs. [61,62].
All of the 2PI calculations cited above have used a counterterm method to perform the renormalization. The introduction of multiple sets of vertex counterterms is required, and the complexity of the procedure is such that it is unknown how to extend it to the 4PI theory. However, since the introduction of higher order variational vertices is numerically very difficult, one might be tempted to ignore vertex corrections and try to improve previous 2PI calculations by increasing the order of the truncation (typically the loop order). In calculations where infrared divergences play an important role, such 2PI calculations at higher loop order can be useful [47,59,[61][62][63]. However, it is known that nPI formulations with n > 2 are necessary in some situations. Transport coefficients in gauge theories cannot be calculated, even at leading order, using a 2PI formulation [35]. Numerical calculations using a symmetric scalar φ 4 theory have shown the importance of 4PI vertex corrections in 3 dimensions [64], and the breakdown of the 2PI approximation at the 4 loop level in 4 dimensions [32].
There is a general hierarchial relationship between the order of the truncation and the number of variational vertices that can be included [65]. If the effective action is truncated at L loops in the skeleton expansion, the corresponding nPI effective actions are identical for n ≥ L (equivalently, one necessarily works with L ≥ n). In this sense, a 3 loop calculation done within the 3PI formalism, a 4 loop calculation done within the 4PI formalism, etc, is complete. There is additional evidence that an L loop calculation in the nPI formalism should, in general, be done with L = n. In gauge theories, the n loop nPI effective action respects gauge invariance to the order of the truncation [45,46], and one therefore expects that a 3 loop 2PI calculation will have stronger gauge dependence than a 3 loop 3PI one. In QED a 2 loop 2PI calculation (which is complete at 2 loop order according to the hierarchial relationship discussed above) found weak dependence on the gauge parameter [66], while a recent 3 loop 2PI calculation in SU (N ) Higgs theory [60] has found strong dependence on the gauge parameter.
There is evidence therefore that higher order nPI calculations are important and worth-while to persue. Higher order effective actions have been derived using different methods [65,[67][68][69], but little progress has been made in solving the resulting variational equations.
As mentioned above, one major problem is that the renormalization of such theories in 4 dimensions cannot be done (using any known method) by introducing counterterms. This is, in part, the reason that a method has been developed to apply the FRG to a 2PI theory [70][71][72][73]. Similar techniques have been used in a condensed matter context in [74][75][76].
Another significant technical problem is the size of the phase space involved in self consistent calculations of vertex functions. Because of limitations of memory and computation time, very few calculations that include variational vertices have been done, and typically various ansätze are introduced for the vertex functions. In Ref. [77], the authors study Yang Mills QCD in the 3PI approximation, but they work in 3 dimensions and use Padé ansätze for the variational functions. A set of self consistent vertex equations obtained from a 3 dimensional 3PI Yang Mills theory were solved in [78], but the full structure of the vertices was replaced with comparatively simple ansätze. Probably the most complete calculation to date was done in Ref. [79] where the authors study QCD at the 3 loop 3PI level. They use a clever technique to exploit the symmetry of the vertices and simplify the variational equations, but they do not actually solve the fully self consistent integral equations, but rather obtain the ghost and gluon propagators using a separate truncation, and input these results into the vertex calculations.
The ultimate goal of our research program is to do a 4 loop 4PI calculation. We work (so far) with a symmetric scalar theory, in order to avoid the complications associated with the Lorentz and Dirac structures of fields in gauge theories. As described above, there are two main obstacles: a conceptual one (renormalizability) and a technical one (memory and computation time contraints encountered because of the large phase space associated with self consistent vertices). In this paper we develop a method to resolve both of these issues, and test it by performing a 4 loop 2PI calculation. We renormalize the theory using the FRG method that was introduced in [73] at the 3 loop 2PI level. Using this method, no counterterms are introduced, and all divergences are absorbed into the bare parameters of the Lagrangian, the structure of which is fixed and completely independent of the order of the approximation. The RG method should therefore work at any loop order, and at any order in the nPI approximation. In [73] we tested our RG method by applying it to a symmetric 3 loop 2PI calculation. However, at 3 loops the traditional calculation requires only one vertex counterterm, and in this sense does not involve the full complexity of the 2PI counterterm renormalization procedure. One could therefore suspect that the agreement of our 3 loop RG calculation with the standard counterterm calculation is an artifact of the approximation. One of the motivations for the calculation in this paper is to verify that this agreement holds at the 4 loop level. The results of this paper show that a calculation that requires two vertex counterterms using the traditional method, can be done using the RG method by appropriately defining the one bare coupling constant that appears in the original Lagrangian. The success of this calculation is therefore strong evidence that the RG method will also work on the technically more difficult 4PI calculation.
The success of our approach is verfied by comparing results with our previous calculation [32] which used Cartesian coordinates; performed all integrals using fast Fourier transforms (which implement periodic boundary conditions); and used counterterm renormalization.
Cartesian coordinates are required when using fast Fourier transforms, but beyond the 2PI level they are impractical because of the size of the phase space involved when vertex functions are represented on a Cartesian grid. In order to reduce the size of the vertex phase space to a numerically manageable level, we use spherical coordinates. However, switching to spherical coordinates means giving up the speed obtained from fast Fourier transforms. Adequate speed is obtained by exploiting the symmetries of the vertex function, and developing efficient interpolation and integration methods. This paper is organized as follows. In section II we present our notation and the setup of the calculation. In section III we describe our method and derive the flow equations that we will solve. In section IV we give some details of our numerical procedure. Our results are presented in section V, and some further discussion and conclusions are give in section VI.

A. Notation
In most equations in this paper we suppress the arguments that denote the space-time dependence of functions. As an example of this notation, the quadratic term in the action is written: The notation G no·int indicates the bare propagator. The classical action is For notational convenience we use a scaled version of the physical coupling constant (λ phys = iλ). The extra factor of i will be removed when rotating to Eucledian space to do numerical calculations.
Using the functional renormalization group method, we add to the action in (2) a nonlocal regulator term The parameter κ has dimensions of momentum. The regulator function is chosen to have the following properties: lim Q κRκ (Q) ∼ κ 2 and lim Q≥κRκ (Q) → 0. The effect is therefore that for Q κ the regulator acts like a large mass term which suppresses quantum fluctuations with wavelengths 1/Q 1/κ, but fluctuations with Q κ and wavelengths 1/Q 1/κ are not affected by the presence of the regulator. The 2PI generating functionals are calculated from the regulated action (3): The 2PI effective action is obtained by taking the double Legendre transform of the generating functional W κ [J, J 2 ] with respect to the sources J and J 2 and taking φ and G as the independent variables: After performing the Legendre transform, the functional arguments of the effective action φ and G are formally independent of the regulator function and the parameter κ, but the non-interacting propagator does depend on κ. We define Using this notation the effective actionΓ κ [φ, G] can be written where Γ 2 contains all 2PI graphs with two and more loops. We define an effective action that corresponds to the original classical action at the scale µ: Throughout this paper we use the notation Γ = −iΦ where both Γ and Φ carry the same subscripts or superscripts. For example, the interacting part of the action is written To make the equations look nicer we also define an imaginary regulator function R κ = −iR κ (the extra factor i will be removed when we rotate to Eucledian space). Using this notation equations (8 -10) are rewritten We extremize the effective action by solving the variational equations of motion for the self consistent 1 and 2 point functions. These self consistent solutions will depend on the parameter κ and are therefore denoted φ κ and G κ . We will work throughout with the symmetric theory by setting φ κ = 0. In figure 1 we show the diagrams that contribute to Φ int to 4 loop order, and the names we will use for these diagrams.

C. Kernels
We define a set of n-point kernels by taking functional derivatives of the effective action Since we work with the symmetric theory, we consider only kernels with n = 0 and we suppress the corresponding 0 index in the superscript. Substituting the self consistent solutions into the definition of the kernels (14) we obtain We introduce specific names for the kernels we will write repeatedly: The kernels Σ, Λ and Υ have two, four and six legs, respectively. The Fourier transformed functions are written Σ(P ), Λ(P, K) and Υ(P, K, Q). The diagrams that contribute to Λ and Υ in the 4 loop approximation are shown in Fig. 2. We note that the 4-kernel contains 2 loop diagrams that involve nastly overlapping divergences. When we use the RG method, any kernel that contains subdivergences is calculated from a flow equation. This is explained in detail in the next section. produced it (see Fig. 1).

A. Flow equations
Using the chain rule, and the fact that Φ int and therefore Λ (m) does not explicitly depend on κ, we find In momentum space (restoring arguments) this equation has the form Equation (18) is an infinite heirarchy of coupled integral equations for the n-point kernels.
This structure is typical of continuum non-perturbative methods, for example Schwinger-Dyson equations and traditional (1PI) RG calculations. In our formalism however, the hierarchy truncates at the level of the action. This can be seen immediately from equation (18) since it is clear that when the effective action is truncated at any order in the skeleton expansion, the kernel on the right side of (18) is zero when the largest number of propagators that appears in any diagram in the effective action is less than m + 1. As will be explained in section III C, the hierarchy can be truncated at an even earlier level.
We show below how to rewrite the flow equations (18) in a more convenient form. The stationary condition is and using equations (11 -14, 16) the variation produces a Dyson equation for the nonperturbative 2-point function in terms of the 2-kernel Σ: Using (20) we have The first two equations in the hierarchy (18) can now be rewritten using (16,21) as These equations are shown in Fig. 3. The grey boxes denote the kernels Σ, Λ and Υ (the specific kernel is indicated by the number of legs attached to the box). The crosses indicate Finally we note that the flow equation for the 2-kernel Σ can be rewritten in terms of the Bethe-Salpeter (BS) vertex. Iterating equation (22) we obtain with In the 2PI formalism, we can also define non-perturbative vertices in terms of the change in the effective action with respect to variations in the field evaluated at the stationary point. These are usually called 'physical' vertices. The physical 4 point function, which we call V , can be written in terms of the BS vertex as We note that the vertex V involves a resummation in all three (s, t and u) channels, but using our shorthand notation which suppresses indices, the three channels combine to produce the factor (3) in equation (26). Details of the derivation of (26) are given in Refs. [31,72].
In order to do numerical calculations, we rotate to Eucledian space by defining the Eucledian variables: The extra factors of i in the definitions of λ E and R E remove the factors that were introduced in the definitions λ phys = iλ (under equation (2)) andR = iR (under equation (10)). From this point forward we suppress the subscripts which indicate Eucledian space quantities. The flow equations (22,23) and the BS equation (25) have the same form in Eucledian space.
The Eucledian space Dyson equation is and the physical vertex becomes The regulator function we use has the Eucledian momentum space form

B. Tuning
Physical considerations dictate the renormalization conditions which are enforced on the non-perturbative quantum n-point functions that are obtained at the κ = 0 end of the flow.
We use standard renormalization conditions: The quantum n-point functions in (32) are obtained by solving the integro-differential flow equations (22,23), starting from some set of initial conditions at κ = µ. The regulator function R κ (equation (31)) is chosen so that at the scale κ = µ the theory is described by Since RG procedure involves initializing the flow equations at an ultraviolet scale κ = µ, and also enforcing renormalization conditions at the κ = 0 end of the flow, we must address the question of whether or not the initial and renormalization conditions can be defined consistently.
Consider an arbitrary n-point kernel of the form whereΛ (m) (P 1 , P 2 . . . ) is the result obtained from integrating the corresponding flow equation, and C(P 1 , P 2 . . . ) is a κ independent integration constant. In the limit κ → µ we require that the vertex function approaches a momentum independent constant: Comparing equations (33, 34) we have so that (33) becomes Now we look at the κ → 0 end of the flow and determine the conditions under which we can enforce the renormalization condition Λ 0 (0, 0, . . . ) = −λ. We start by adding and subtracting two different terms to the original vertex, and grouping into square brackets the differences we will consider below. This gives The renormalization condition is satisfied if the square brackets go to zero in the limit that κ and the momentum arguments go to zero. Setting κ = 0 and using equation ( The second term in (38) is zero in the limit µ P , since this is (by construction) the limit in which loop contributions are suppressed and the momentum dependence of the vertex disappears. We have therefore shown that the renormalization condition will be consistent with the initial condition if In the next section we will show that this condition is satisfied if the truncation of the hierarchy in (18) is performed correctly.
We comment on the fact that the discussion above is misleading in one important way. It appears that the bare vertex λ µ cancels exactly when we go from equation (36) to equation (37). If this were true, the initial condition and the renormalization condition would be unconnected to each other, and tuning would not be possible. The apparent cancellation occurs because the self consistent nature of the set of coupled equations for kernels with different numbers of legs is not evident when we consider one kernel in isolation.
Finally, we note that renormalization conditions are typically defined on vertices constructed from resummed kernels (see equation (32)) and not the kernels themselves. However, the condition derived above (39) is still sufficient to guarantee the consistency of the procedure. The crucial point is that the kernels are 2-particle-irreducible (see Fig. 2) and therefore when they are chained together to form a resummed BS vertex, no additional subdivergences are produced. It is easy to see how this works in our calculation, where we only need to enforce a renormalization condition on the 4 point function (as will be explained in the next section). When we define the renormalization condition on the BS vertex instead of the kernel Λ, the result is only a shift in the final value of the bare coupling λ µ that is produced by the tuning procedure.

D. Truncation
Now we return to the issue of the truncation of the hierarchy of flow equations in equation (18). The kernels obtained from direct functional differentiation of the action using (15) will automatically satisfy the correct flow equations (18). As explained in section III A, this is just an obvious application of the chain rule. Next we observe that if a given kernel obtained from functional differentiation satisfies the condition (39), then using the analysis of the previous section we know it will also satisfy Λ In order to identify the terminal kernel, we need to know under what circumstances the condition (39) will be satisfied by a given 2PI kernel Λ (m) κ . We start by looking at an example where it will not. We consider the self energy diagram on the right side of Figure 4 which gives Σ 0 (P ) = − λ 2 6 dQ dL G 0 (L)G 0 (L + Q)G 0 (P + Q) .
This sunset contribution to the self energy is produced by the BBALL diagram in the effective action (see Fig. 1). The quantity Z in equation (39) now takes the form where in the last line we have expanded around P 2 = 0. The prime denotes differentiation with respect to Q 2 and the dots represent terms that are higher order in P 2 /Q 2 . It is clear that the divergent L integral is unaffected by the subtraction, and therefore we cannot conclude that Z → 0 when P approaches zero.
Another example is the contribution to the 4-kernel Λ from the last diagram in the first line of Fig. 2. If we route the momenta as shown in the right side of Fig. 4 and rescale momenta as described above, the divergence in the L 2 integration is unaffected by the subtraction, and therefore the condition (39) is not satisfied.
In general, any loop that does not necessarily carry one of the external momenta is a 'bad' loop. If a kernel does not have any bad loops, it satisfies (39) and its flow equation does not have to be solved. The smallest of these kernels is the terminal kernel. In the following two subsections we explain how to apply these ideas to the 2PI calculation at the 3 and 4 loop level.

3 loop
If we truncate the effective action at the 3 loop (BBALL) level, the self energy includes the sunset diagram which has a bad loop, as explained above. The kernel Λ has no bad loops, since the 1 loop BBALL contributions in Fig. 2 always carry external momenta. The vertex Λ is therefore the terminal kernel and can be substituted directly into the Σ flow equation. In order to satisfy the initial condition, we replace the tree vertex (the EIGHT contribution in Fig. 2) with the bare parameter −λ µ . The Λ flow equation (23) is not affected by any constant shift of the 4-kernel. Thus we see that we can obtain, directly from the action, an expression for the 4-kernel that obeys the initial condition and satisfies the correct flow equation. We have also explicitly checked that the results are the same as those obtained from solving the coupled set of Σ and Λ flow equations.

4 loop
At the 4 loop level, the kernel Λ contains a bad loop (as illustrated in the right side of Fig. 4) and does not satisfy (39). We therefore cannot not use the explicit expression for the 4 kernel Λ shown in Fig. 2 directly in the flow equation for the 2 kernel (equation (22)), as we did at the 3 loop level. The kernel Υ (the bottom line of Fig. 2) contains only 1 loop diagrams that always carry external momenta, and therefore satisfies (39) and is the terminal kernel. Since there is no bare 6-vertex in the Lagrangian, we know the integration constant should be set to zero. We therefore substitute the result for Υ shown in Fig. 2 directly into the Λ flow equation. The Σ and Λ flow equations must then be solved self consistently.

A. Procedure
We initialize the flow of the 2 and 4 kernels and we take the propagator in the ultraviolet limit as G −1 µ (P ) = P 2 + m 2 µ . Schematically the numerical procedure can be described as follows: 1. Choose values for the physical mass m and coupling λ. We use always m = 1, which means we give all dimensionful quantities in mass units.
3. Compare the chosen and found values for the mass and coupling, adjust the bare values up or down accordingly, and repeat all steps until the renormalization conditions are satisfied to the desired accuracy.

B. Parameters
The differential equations are solved using a logarithmic scale by defining the variable t = ln κ/µ, in order to increase sensitivity to the small κ region where we approach the quantum theory. We use κ max = µ = 100, κ min = 10 −2 and N κ = 50. We have checked that our results are insensitive to these choices. In addition, we have checked for possible dependence on the form of the regulator function by using a generalization of (31): The original expression (31) corresponds to the choice of exponent z = 1. Using z = 1/2 and z = 1/4 produces results that are virtually identical.
The 4-dimensional momentum integrals are written in the imaginary time formalism as with m t = 2πT . Numerically we take N t terms in the summation with β = 1 T = N t a t where the parameter a t is the lattice spacing in the temporal direction.
The integrals over the 3-momenta are done in spherical coordinates and using Gauss-Legendre integration. We use typically N x = N φ = 8 points for the integrations over the cosine of the polar angle and the azimuthal angle. The dependence on these angles is weak, and results are very stable when N x and/or N φ are increased. To calculate the integral over the magnitude of the 3-momenta, we define a spatial length scale analogous to the inverse temperature L = a s N s where a s is the spatial lattice spacing and N s is the number of lattice points.

C. Restrictions
The numerical method replaces a continuous integration variable with infinite limits by a discrete sum over a finite number of terms. For numerical accuracy, we need generically that the upper limit of the sum is big and the step size is small. This means we require Finally, we note that it is well known that scalar φ 4 theory in 4-dimensions is noninteracting if it is considered as a fundamental theory valid for arbitrarily high momentum scales (quantum triviality), but the renormalized coupling is non-zero if the theory has an ultraviolet cutoff and an infrared regulator. In our calculation the mass m regulates the infrared and the lattice spacing parameter provides an ultraviolet cutoff.

V. RESULTS AND DISCUSSION
We use a t = 1/10, a s = 1/5 and N s = 18 and vary the temperature by changing the number of lattice points in the temporal direction. The renormalization is done at N t = 37 and the highest temperature we consider is T = 2 which corresponds to N t = 5. In Fig. 5 we show the BS and physical vertices as functions of the temperature. The graph agrees well with the results of our previous calculation [32] which used two separately defined counterterms. Some differences between the two calculations are expected, due to the different boundary conditions that must be implemented using the spherical and cartesian/fft methods. -M(0) and -V(0) In Fig. 6 we show the dependence on the renormalization scale. The two curves correspond to the physical vertex V (0) versus temperature with the renormalization done at two different temperatures. The scale dependence of the calculation is indicated by the shaded grey region between the two curves, and is very small. In Fig. 7  cross at N s = 18, in order to provide the best means of comparison. The curves are shifted so that they cross at N s = 18.
In order to test the renormalization, we check that the results are unchanged when p max is increased while ∆p ∼ 1/L is held fixed. This is done by increasing N s while adjusting a s so that L = a s N s is constant. In Fig. 8

VI. CONCLUSIONS
In this paper we have done a 4 loop 2PI calculation in a symmetric φ 4 theory. We have renormalized the theory using the FRG method that was introduced in [73] at the 3 loop 2PI level. Using this method, no counterterms are introduced, and all divergences are absorbed into the bare parameters of the Lagrangian, the structure of which is fixed and completely independent of the order of the approximation. We therefore expect that our RG method will work at any order in the nPI approximation. The next step in our program is to apply our method to a 4 loop 4PI calculation. The structures of the flow and Bethe-Salpeter equations are the same [80], but there is now a variational 4 vertex that must be calculated self-consistently. In spherical coordinates, the phase space for this vertex is comparable with that of the 3 dimensional self consistent 4PI vertex function that was calculated in [64]. This calculation is currently in progress.