nPI Resummation in 3D SU(N) Higgs Theory

We test the utility of the nPI formalism for solving nonperturbative dynamics of gauge theories by applying it to study the phase diagram of SU(N) Higgs theory in 3 Euclidean spacetime dimensions. Solutions reveal standard signatures of a first order phase transition with a critical endpoint leading to a crossover regime, in qualitative agreement with lattice studies. The location of the critical endpoint, x sim 0.14 for SU(2) with a fundamental Higgs, is in rough but not tight quantitative agreement with the lattice. We end by commenting on the overall effectiveness and limitations of an nPI effective action based study. In particular, we have been unable to find an nPI gauge-fixing procedure which can simultaneously display the right phase structure and correctly handle the large-VEV Higgs region. We explain why doing so appears to be a serious challenge.


Introduction
Thermal or off-equilibrium dynamics of hot nonabelian gauge theories have applications in heavy ion physics (see e.g. [1] and references therein) and in early Universe cosmology [2,3]. An important feature of nonabelian gauge theory is asymptotic freedom, which makes the coupling smaller at shorter length (or higher energy) scales. Naively this means that perturbation theory should work better at higher temperatures, where the relevant energy scale T is large. In fact this is only partly the case. As shown already in 1980 [4,5], the behavior on scales > 1/T is in fact that of a 3-dimensional theory, which goes rapidly to strong coupling at longer distances. Therefore the long-distance T behavior of nonabelian gauge theory is strongly coupled at any temperature.
There is a dearth of tools for computing real-time dynamics of nonabelian gauge theories when nonperturbative physics is involved. This frustrates efforts to understand dynamics, both at strong coupling and at weak coupling. A method which has shown much promise in scalar and Yukawa theories is the n-particle irreducible (nPI) method [6][7][8][9][10][11][12]. While the motivation for developing such methods lies largely in the hopes that they can be applied to nonabelian gauge theories [13,14], almost no work in that direction has occurred yet. There have been some results in abelian theories [15][16][17][18], and some arguments as well as a computation demonstrating that a 3-particle-irreducible treatment of QCD would automatically capture the leading perturbative effects relevant for transport and equilibriation [19]. But no one has made a concerted effort to apply the nPI approach to nonabelian gauge dynamics.
In a previous paper [20] we took a first step in this direction, by applying the 3-particle irreducible (3PI) method to the study of Yang-Mills theory in three Euclidean dimensions. The main motivation was to test out the methodology in a context where we do have other computational tools at our disposal, so we can appraise whether it is successful before undertaking the more challenging problem of applying 3PI methods to dynamics. But if successful, the 3PI method could still have real utility as a potentially faster or more efficient method of studying 3D theories, which are in fact of intrinsic interest. In particular, the 3D theory captures the large-distance nonperturbative physics of weakly-coupled hot gauge theories alluded to above, at least at the thermodynamical level.
Unfortunately, in that work we were only able to solve the 3PI problem for pure gauge theory in 3 dimensions. There are few long-distance sensitive observables in that theory, so we lacked gauge-invariant measurements to match to (lattice) nonperturbative calculations in 3D Yang-Mills theory and did not find effective ways to test whether the method "works." In the present paper we intend to address this by extending our nPI treatment to 3dimensional Yang-Mills Higgs theory, a theory which has a nontrivial phase structure. We will study whether the nPI approach can successfully reproduce the phase structure of the theory, a nontrivial and nonperturbative test of the technique. We emphasize that our purpose is as a test of whether the nPI approach in gauge theory can reproduce nonperturbative phenomena. The goal is not to improve our understanding of 3D Yang-Mills Higgs theory, which has been thoroughly studied using lattice gauge theory techniques [21][22][23]. In the remainder of the introduction, we will explain both 3D Yang-Mills Higgs theory and the nPI approach in a little more detail.
The study of three dimensional nonabelian gauge theory is motivated by its relationship to electroweak theory and QCD via dimensional reduction [24][25][26][27][28]. At high temperature T Λ QCD , QCD exhibits a natural separation of scales g 2 T gT T so that non-zero bosonic and all fermionic Matsubara modes become heavy compared to the soft scales of the theory. These modes can be integrated out to obtain an effective 3 dimensional description of the soft physics, which is precisely SU(3) Yang Mills coupled to an adjoint scalar with gauge coupling g 2 3D = g 2 4D T and mass m 2 A = g 2 (N/3 + N f /6)T 2 (identified with the 0-mode of the A 0 component of the 4D gauge field). If one is only interested in physics at the supersoft scales, this can be taken one step further by integrating out the A 0 field to yield pure 3D Yang-Mills.
Yang-Mills theory, QCD and electroweak theory are known to undergo a phase transition [29][30][31][32] over certain ranges of the model parameters. Naturally, for physical values of these parameters, one would ask whether we are in a first order, second order or crossover regime. 3D effective models could potentially shed some light on this matter, except that for QCD, where the effective 3D description is an SU(3)+adjoint Higgs theory, T c ∼ Λ QCD . Thus, in the vicinity of the QCD phase transition (or crossover) the effective description breaks down, since the underlying assumption of weak 4D coupling and a separation of scales is not valid. A 3D effective model is still useful for studying the nonperturbative infrared dynamics of hot QCD, just not at temperatures in the vicinity of the scale Λ QCD . However, the situation is different for electroweak theory near its phase transition.
The effective 3D description of electroweak theory resulting from dimensional reduction is an SU(2) × U(1) gauge theory coupled to both fundamental and adjoint scalars. The adjoint scalars arise via the dimensional reduction, while the fundamental scalar is identified with the 4D Higgs field. In practice, an accurate study of the 4D theory does not require such elaborate field content; rather, quantitative predictions can be made by considering the much simpler SU(2) + fundamental case [21-23, 33, 34]. Then, as a further refinement, one may study the effects due to the inclusion of an adjoint field [35,36], as well as a U(1) gauge field [37]. Or, in the context of GUTs, the model with an SU(5) or SU(3) × SU(2) gauge group may be of interest [38,39]. These models have received a fair amount of attention in the past due to the significance of a phase transition on electroweak baryogenesis [2,3]. For the models considered, a first order phase transition at physical values of the Higgs mass has been ruled out.
We show a cartoon of the phase diagram of SU(N ) fundamental-Higgs theory in Fig. 1. It is parametrized by two dimensionless variables x and y, describing the scalar self-coupling and renormalized mass respectively, normalized to the appropriate power of g 2 . In terms of the 4D thermal theory, x is predominantly set by m H /m W and y is predominantly set by the temperature. Therefore an evolution in temperature in the 4D theory corresponds to a nearly vertical line on the phase diagram; moving horizontally is changing the vacuum parameters of the theory. The diagram has a first order line which ends at a second-order 3D Ising universality [23] endpoint at (x c , y c ); if y is varied holding x < x c fixed, the system encounters a first order phase transition. But the transition is not a symmetrybreaking phase transition, and no order parameter distinguishes one phase from the other, similar to the liquid-gas phase transition (which is in the same universality class). We have also indicated upper and lower metastability lines, which show how deeply the system can "superheat" or "supercool" before encountering spinodal instability. The locations of these lines cannot be rigorously defined; but they will enter in our analysis so we indicate them for completeness.
At small x the phase transition can be studied perturbatively by computing the oneloop effective potential for the Higgs vacuum expectation value (VEV). Indeed one finds that in this region, the phase transition is first order. However, the perturbative treatment then goes on to predict a first order phase transition for all values of x! Higher-order computations [40,41] show that the perturbative expansion parameter is actually x, indicating that perturbation theory fails for large values of x, which turns out to mean x > ∼ 1/10. Therefore, the end point and crossover must be resolved nonperturbatively. Since the lattice has already provided us with a very accurate determination of the phase diagram, we are able to use these results to test the reliability and accuracy of an alternative nonperturbative approach to the lattice, namely that of nPI resummation in a gauge theory setting. In this paper we will study the application of the nPI (specifically 2PI) formalism to SU(N ) Higgs theory.
In the context of a hot gauge theory, the use of an nPI based resummation scheme is primarily motivated by the extremely poor convergence of a weak-coupling expansion [42], since it provides a systematic procedure for reorganizing a perturbation series. Our approach here is along a trajectory which differs from many of the previous works on this subject mentioned earlier. That is, by applying the nPI formalism to SU(N ) Higgs theory our goal is to directly solve the resulting self-consistent, Schwinger-Dyson (SD) resemblant 1 integral equations in a manner reminiscent of [20], and then subsequently derive gaugeinvariant quantities from the solutions.
An nPI effective action Γ[φ, G, ...] generates equations of motion for n-point resummed vertices by variation with respect to these n-point functions. We will specifically consider a three-loop truncation of the case n = 2, which in terms of diagrams can be interpreted as resumming one-and two-loop self-energy topologies. (In [20] we treated the pure-gauge theory at the n = 3 level, that is, including as well a self-consistent one-loop resummation of three-point vertices. The result established that the corrections to these vertices are small, so we have avoided this technical complication.) By solving the resulting "SD" equations, we can subsequently compute the gauge-invariant scalar condensate φ † φ as a function of the parameters x and y on the phase diagram. Then, at a specific value x, from the behavior of φ † φ over a range of y we can infer whether or not we are in the crossover or first order phase transition region. This will allow us to bracket and locate the critical end point.
We will present the technical details of the computation for a single complex scalar field in representation R of SU(N ). Results will be given for N = 2 (fundamental repre-sentation) in Landau and Feynman gauges, however it should be noted that the method straightforwardly generalizes to the inclusion of additional scalar fields, higher representations, and larger gauge groups.
The text is organized as follows: in Section 2 we will present the three-loop truncated 2PI effective action for SU(N ) Higgs theory, as well as the SD equations that it generates. Additionally, some details pertinent to regularization and renormalization will be reviewed here. In Section 3 we will present certain extensions to the algorithm described in [20] which are needed to solve the 2PI equations of motion numerically. In Section 4 we will give an overview of the results, as well as derived quantities such as the scalar condensate and the location of the critical end point. Finally, throughout Sections 3 and 4 we will also discuss the properties of the effective action, and comment on the overall effectiveness of the method.
2 SU(N ) Yang-Mills + Higgs theory in the nPI formalism

General remarks and notation
It is useful to begin by reviewing a number of the basic conventions that are used throughout. It should be assumed that T a R is a generator of some representation R of SU(N ). The fundamental and adjoint representations are denoted by F and A respectively, and d R is and additional group theory identities needed in this computation can be found in [43]. Following gauge fixing, the Lagrangian can be divided into a Yang-Mills component and a Higgs component, so that L = L YM + L φ (in general covariant gauge as written). Note that we have not fixed to R ξ gauge; our gauge fixing only acts on the gauge degrees of freedom. We will discuss this more in Subsection 2.3. We define the dimensionless ratios which is the same as the definition introduced in Ref. [21] and commonly used throughout the literature. In Eq. (2.4) an additive counterterm has been explicitly included to cancel the divergent two-loop self-energy graphs (its value is given in Appendix B). This leads to a scale dependence in m 2 and y accordingly; for a fundamental SU(2) Higgs we have The mass renormalization scale is fixed at µ = g 2 throughout. For phenomenological applications, the relation between x, y and 4D theory parameters are [21] x = −0.00550 + 0.12622 m H 80.6GeV 2 (2.7) assuming a value of g 4D = 2/3 for the 4D gauge coupling.

The three-loop 2PI effective action
The 2PI effective effective action Γ[G ij ] is formally defined as the Legendre transform of the generating function of connected diagrams W [K ij ] with respect to a two particle source [7]. Using the generic label Φ i for fields, Even correlation functions can be obtained by differentiation with respect to K ij . For instance, yields the two-point function G ij . For the two-point functions of SU(N ) Higgs theory, we can assume that G ij is proportional to the color identity of the corresponding species, and hence so is K ij . Then, for a rotationally symmetric Lagrangian i.e., the presence of K ij does not alter the global rotational invariance of the original Lagrangian. So in fact, Eq. (2.10) generates the connected two-point function. The consequences of this statement in the context of a spontaneously broken gauge theory will be discussed towards the end of this section, but for now we can proceed with the Legendre transform (2.12) In setting K ij = 0, equations of motion for G ij are obtained from the stationarity condition The solutions we seek correspond to extrema of Γ[G ij ]. Specializing now to the field content of SU(N ) Higgs theory, we can write Γ = Γ YM + Γ φ and explicitly state the loop expansion, which we will truncate at three loops. Γ YM is defined so that it contains those diagrams encountered in the pure Yang-Mills problem while Γ φ contains the additional diagrams which arise when a single arbitrary representation Higgs field is included. Using the diagrammatic notation G µν = (2.14) without loss of generality we write the two-point functions as .
All vertices appearing in the 2PI effective action are at tree level. These are drawn as with the corresponding expressions given in Appendix A. Finally, the Higgs mass renormalizes at two loops; it is necessary to subtract the divergence with an additive counterterm of the form m 2 = m 2 φ + δm 2 , with the corresponding vertex Explicitly factoring out minus signs due to ghost loops, we have For an n-loop pure Yang-Mills planar diagram, the tracing over internal color indices generically results in an overall color factor of (N 2 − 1)N n−1 . Furthermore, since an n-loop vacuum bubble is also proportional to g 2(n−1) , one finds as earlier that factors of g 2 always appear in the form of the 't Hooft coupling g 2 N . Hence, for the pure gauge problem, the natural scale is not g 2 , but rather g 2 N . The Higgs contribution is These diagrams have a somewhat more complicated dependence on N (the associated color factors are stated in Table 1). We will nevertheless continue to use g 2 N as the standard mass scale, but for clarity, units of g 2 N will be explicitly stated throughout. The power of the 2PI formalism becomes apparent when we perform the variation of Γ with respect to G T , G L , ∆ and D. For instance, from δΓ/δD = 0, we have with (omitting charge arrows) Π φ (p) = + 2 + 1 2 Equations of the type Eq.  truncation of the effective action. By solving this equation self-consistently in a three-loop truncation, we fully resum one-and two-loop self-energy topologies to all orders.
Eq. (2.25) contains terms that are linearly and logarithmically divergent; in dimensional regularization, only the logarithmic divergences appear explicitly as 1/ 's, and these are subtracted by the counterterm. This implies that the entire computation must be performed in MS, which requires the analytic continuation of these integrals to D dimensions. The regularization procedure which we adopt is described at length in [20]; to quickly recap the key points, consider the tadpole graph is an arbitrary function of q, this integral would need to be performed numerically; in doing so we must set D → 3. To implement dimensional regularization, we adopt a procedure of "addition and subtraction," as follows: The rightmost term is simple enough that it can be computed analytically (in MS its value is zero), and the leftmost term is now only logarithmically divergent. Thus, we have removed the linear divergence by subtracting 1/q 2 , and now the next step is to remove the logarithmic one. At large momenta, and near 3 dimensions, D(q) can be expanded as where we have been careful to keep O( ) corrections in the 1/q 3 term. Now, we can add and subtract the subleading term, The first line of Eq. (2.29) is finite, so we can set D = 3 and perform the integral numerically. What we have effectively done is shuffled all of the dependence into terms which can be integrated analytically. Thus the regularized expression for A has the form We can then subtract the 1/ divergence with the counterterm, and take the limit → 0. Note here that the subleading term appears with a mass ω. Its value is arbitrary, but it must be included, otherwise one would introduce an IR divergence where originally there was none. For simplicity, we can set ω = g 2 N noting that the final results of the calculation are ω independent. Though it is certainly permitted, it is not a requirement that ω be set to the scalar mass m (and our reasoning for not doing so is explained in Appendix B).
Other diagrams which appear in Γ are regularized in much the same fashion. In the end we need to compute all of the one-and two-loop gluon and Higgs self-energy diagrams which appear in perturbation theory (ensuring that IR divergent diagrams are not introduced inadvertently); the results of this computation are contained in Appendix B.

Remarks on gauge fixing
Since our approach is diagrammatic and is founded on trying to determine correlation functions of gauge dependent fields, we are obliged to perform some sort of gauge fixing. We have chosen covariant gauge with gauge-fixing functional ∂ µ A µ , that is, a gauge choice which does not make reference to the scalar field one-point function or vacuum value. This choice differs from what is usually done in perturbative treatments of the Higgs phase, and requires some explanation. First we will argue that one can gauge fix as we do here; then we will explain why we believe it is preferable.
That the gauge-fixing approach we have adopted is possible, has already been explained clearly by Buchmüller, Fodor, and Hebecker [44] in the context of the electroweak phase transition. Perturbatively, in the broken symmetry phase we expect the typical contribution to the path integral to have a nonvanishing scalar VEV Φ i ; but since the gauge fixing does not eliminate the integration over the global gauge rotations, there are equal contributions from every field direction choice, and the ensemble average of Φ i is zero. However the existence of a VEV will still appear as a delta-function contribution to the scalar two-point function, so the approach will still capture that physics. Nonperturbatively, while infrared gauge fields are suppressed, we do not expect them to vanish, and they can still destroy any infinite-range order in the scalar field. If this is the case then the scalar two-point function will not in fact have a strict delta-function contribution. Instead it will have a very sharp structure near zero momentum, corresponding to long (but not infinite) distance correlations in the scalar field. Indeed, we expect this must be the correct behavior, since the symmetric and Higgs phases are analytically connected and are not distinct in the sense of being distinguished by a true order parameter. But the existence of infinite-range Higgs-field correlations in part but not all of the phase diagram would constitute an order parameter and would forbid an analytic connection between the symmetric and Higgs sides of the transition line. Now let us consider the alternative approach. It is to include explicitly a one-point source for the scalar field, The value of W [0, 0] is gauge-invariant [45], but the inclusion of nonzero J i explicitly breaks gauge invariance. Naturally we are then only interested in the J i → 0 limit. Depending on our gauge-fixing procedure, this limit may or may not be smooth. That is, we can interpret the gauge-fixed one-point function as a directional derivative where 0(ϑ) means "zero is approached along a direction ϑ on the manifold of SU(N ) rotations." Perturbatively, we expect the symmetric phase to display smooth behavior, so an approach from any direction will yield the same result, consistent with a zero VEV. In the broken phase, W is expected to develop a conical singularity, so the direction has a significance. Nonperturbatively, and still working in covariant gauge, we have just argued that we do not in fact expect infinite-range correlations in the scalar field when J i is taken to zero; so the behavior of W near J i = 0 should always be differentiable, albeit the derivatives can become very large. Therefore, in covariant gauge we expect that there should be no VEV, even if J i is taken to zero along a particular direction. Therefore the introduction of J i changes nothing; the J i → 0 limit is nonsingular and the VEV of the field vanishes.
Alternatively, we can change our gauge-fixing procedure so that it makes explicit reference to the VEV Φ i -that is, we can use R ξ gauge. The gauge-fixing choice introduces into the action the gauge-fixing term where Φ i is the VEV and φ j is the field. This is balanced as usual by the appropriate Fadeev-Popov determinant. Physically, the role of covariant gauge fixing can be understood as using up the gauge freedom to force ∂ µ A µ to be as small as possible, which minimizes the total size of fluctuations in the gauge fields. R ξ gauge is instead a compromise, in which the gauge fixing is used to try to minimize ∂ µ A µ (gauge field fluctuations), but also to minimize fluctuations in the components of the scalar field not in the direction of the VEV Φ i (pseudo-Goldstone fluctuations). The limit ξ → 0, Landau gauge, is when all gauge freedom is used to control gauge field fluctuations. The opposite limit, ξ → ∞, Unitary gauge, is when all freedom is used to align the scalar in the direction of its VEV. One challenge with this approach is that in the current context the value of the VEV Φ i must be determined self-consistently as part of the procedure. In general the VEV will depend on ξ [46], growing larger at large ξ as more fluctuations are forced into the VEV. Perturbatively this effect is suppressed by g 2 in 4 dimensions. Here it will be suppressed by x. Since we are interested in a regime where perturbation theory requires resummation, the ξ dependence can be large.
The first problem with this approach is that whether Φ i vanishes or not would constitute an order parameter, but we know that there should not be an order parameter for this system. Second, it is possible that there are multiple self-consistent solutions for Φ i , in which case it is not clear which to use. Finally, the existence of a VEV Φ i for a given x, y value can and almost surely will depend on ξ, as large ξ biases the gauge fixing towards the development of a VEV. Therefore we anticipate that the details of the transition will have no stability as a function of ξ, in other words, the methodology would not be reliable. (Similar issues are discussed in Appendix A of Ref. [47].) On the other hand, we emphasize that using covariant gauge will present considerable challenges when the transition is strong or when the value of y places us deep in the Higgs phase. In this case we will be keeping track of very long-distance correlations in the φ field via the corresponding very low-momentum structure in the two-point function. As we will explain in the following, this proves numerically challenging, but it can also be a problem from the point of view of the convergence of the loop expansion and the stability of the solutions we find within the space of possible G, D, ∆ choices. It is also not clear what ξ dependence the phase diagram will display in covariant gauge; if the dependence is strong, it indicates a problem with the method's reliability.

Variational Ansätze
In [20], we extensively described an algorithm which can be used to extremize the effective action when only gauge fields are present. Now we have to address the additional complications which arise due to the presence of a Higgs field. In the symmetric phase, the presence of the Higgs does not really change much at the technical level, and obtaining self-consistent solutions for the gauge field and Higgs propagators proceeds much as earlier.
To begin, we will review the details of the functions which enter into the problem. Since we have assumed a general covariant gauge, we are attempting to solve self-consistently for the following 4 functions: G T (p), G L (p), ∆(p) and D(p), which are respectively the transverse and longitudinal gauge field propagators, the ghost propagator and the Higgs propagator. We can opt for the most part to simplify the problem further by working in Landau gauge, where G L falls out of the picture; however, computations in Feynman gauge do require a treatment of G L .
To realize the extremization, we will specify Ansätze for these functions in terms of a finite set C = {c i } of variational coefficients, such that the variational equations become Writing the functions in terms of a finite number of parameters in this way replaces the infinite-dimensional functional space with a finite-dimensional subspace; and the problem becomes finding the extremum in this subspace. By increasing the size of C we enlarge the space of allowed functions, and the true extremum should be more closely approached.
Our choice is to fit the self-energies as rational functions (Padé approximants), since this gives a very flexible class of smooth functions. Specifically, for some {a i } ∪ {b j } ⊂ C, we will define Here we have incorporated the one-loop large-p behavior exactly and allowed the rest of the self-energy to be determined by extremization. Previous work [20] shows that third order Padé approximants are sufficient and we will use them here. The resulting SD equations have the simple form Π Ansatz is the gluonic analogue of Eq. (2.25)) and similarly for Π L , Σ and Π φ .
However, we anticipate that the scalar propagator may display a very narrow structure near zero momentum. Therefore we will add to the scalar propagator Ansatz an additional term: .
The added term is designed to allow for a sharp structure at small p; its form has been chosen phenomenologically. Technically the functional form allows for D(p) ∝ p −2 smallp behavior, whereas we expect that lim p→0 D(p) should be a constant. 3 However, the extremization procedure is not sensitive to a turnover at very small p, so this functional form near p = 0 is not very important. Generally, in the symmetric phase extremization setting R 0 to zero produces essentially the same extremum as allowing the term to be nonzero, whereas deep in the Higgs region the inclusion of R 0 is essential to getting a good solution to the SD equation, which is modified to −D −1 (p) + D (0)−1 (p) = Π 2PI φ (p). We can also define the renormalized scalar "condensate" as The twice-subtracted integral is both UV and IR finite; as indicated it equals the expectation value of the field squared when renormalized usingμ = g 2 N . This condensate will be useful in distinguishing between coexisting phases. In Eq.
as well as the variation of Γ with respect to c Since convergence to the perturbative limit is only necessary at large p, looking back at Eq. (3.5) and Eq. (3.6), we opted to include the one-loop contributions with an additional IR suppression factor of p/(p + ω). At sufficiently small momenta, a linear term in the denominator of a propagator can lead to the formation of a pole. The gauge fields dynamically generate a mass that is sufficiently large to prevent this sort of thing from happening so a suppression factor of this sort is not required. For the Higgs and ghost this is not the case. In solving the problem we have set ω = g 2 N . This arbitrary choice does not affect the final result, since a different choice of ω together with an appropriate shift in the Ansatz parameters leaves the self-energy unchanged.

Initial conditions and root finding
To extremize Γ[G T , G L , ∆, D], we employ and algorithm based on conjugate gradient descent specialized to the problem at hand. We can visualize the root-finding algorithm as a dynamical system where we choose an initial value for the coefficients c i and subsequently follow a flow through the gradient field ∂Γ/∂c i until we reach an attracting fixed point, which corresponds to a solution.
For x below the critical end point, there is a region of metastability where two attractors coexist. We will denote the solution with larger condensate D (and smaller small-p behavior in G T ) as the Higgs solution and the solution with smaller condensate as the symmetric phase solution. We write the propagator (condensate) in the Higgs solution as D − (D − ) and that in the symmetric phase as D + (D + ). In the space of initial guesses for the Ansatz parameters, each solution has a basin of attraction, which we write as γ + and γ − . These basins of attraction do not cover the full space of initial guesses; because of the saddle-like nature of Γ, there is also a set of initial conditions c i ∈ γ 0 which evolve towards divergent values of D. This is to be discussed in greater detail in Section 4.
The number of iterations of gradient descent in the extremization procedure (denoted by N ) can be thought of as time evolution, and we are interested in the results at late times. We can observe convergence of the algorithm by plotting the evolution of the LHS and RHS of the SD equations with N ; this is shown generically in Fig. 2. From this figure it is apparent that convergence is attained, and that the variational Ansatz has captured a choice for the self-energy where the SD equations are quite accurately solved.
To establish two distinct metastable solutions, it is also important not to be fooled by slow convergence to an extremal solution. We illustrate this idea in Fig. 3, which shows algorithm convergence for two cases. On the left, we see convergence from two different starting configurations to two distinct final solutions. On the right, we see slow evolution to a single solution. To distinguish these cases, it is important to fit the N dependence of D (or some other measure of the solution), to test convergence. We find that the fit form gives a good description.

Issues of stability
Is the extremization of Γ a minimization/maximization or a saddlepoint-seeking procedure? It is a saddlepoint-seeking procedure, as can be seen by considering the one-loop value, which is solved trivially by D = D 0 , G = G 0 , ∆ = ∆ 0 . For the sign choice above, this extremum is clearly a maximum for G, but a minimum for ∆, since the ghost enters with the opposite sign. At least in the ultraviolet this property is not affected by additional diagrams, since the tree terms dominate in the UV. This does not cause a problem in practice, since we can alternately extremize with respect to G, D holding ∆ fixed and with respect to ∆ holding G, D fixed. The former involves maximization, the latter involves minimization. This procedure shows rapid convergence, and was used in our previous work [20]. The interpretation of this saddle behavior is benign; it arises because of the peculiarities of gauge fixing and the presence of the "fermionic" ghost species which it introduces.
But we are in trouble if the bosonic part of Γ at fixed ∆ is unbounded from above and below. When we move from one to three loops and we consider the possibility of large infrared contributions in D, we find that precisely this problem arises. The pure scalar diagrams, omitting group theoretic factors, are of the form where we have flipped the overall sign so the one-loop piece opens upwards. With the three-loop term present, Γ is unbounded from above and below. This unbounded behavior becomes important whenever D becomes sufficiently large in some narrow momentum range. For instance, when a small p range around zero supports a finite value of the integral p D(p), then the small p, k, q contribution to the three-loop (basketball) term becomes large, and it diverges as the phase space region supporting p D ∼ 1 goes to zero. Unfortunately, the expected behavior deep in the Higgs phase is precisely that p D(p) should receive a finite contribution from a very narrow momentum range near p = 0. Therefore, the extremum we seek is at best a local maximum as a function of G, D; and in particular, we can expect trouble deep in the Higgs phase. The origin of this problem is that, when the field develops large long-distance correlations, the loopwise expansion of the 2PI functional breaks down. For instance, at the four-loop level we will encounter which diverges still more strongly. The sequence of such divergent graphs is resummed by including the one-point function in the procedure, and stripping the square of the one-point function from the two-point correlator. However, as we have emphasized, any procedure for including the one-point function appears to damage the properties which ensure the possibility of a phase transition endpoint.
Here we work in terms of the two-point function only, which will mean that we are unable to study cases where the solutions become strongly Higgs-like. We anticipate that, when we seek solutions which show strong Higgs-like behavior, we will instead find runaway behavior in our extremization algorithm.
Naively it appears that this problem is less severe at small x where the high-loop Higgs diagrams are suppressed by more explicit powers of x. But the Higgs-phase value of the condensate p D(p) grows as 1/x 2 , so in fact the problem is more severe, not milder, at small x. Therefore it will not be possible to make contact with the perturbative part of the phase diagram.

Analysis and results
We will concentrate on the analysis in Landau gauge (which eliminates the longitudinal gluon propagator), and set N = 2 with the scalar field in the fundamental representation. A comparison with the results in Feynman gauge appears towards the very end, in Section 4.2.
Solutions for the Higgs, gauge and ghost propagator are shown in Fig. 4 (Higgs) and Fig. 5 (gauge/ghost). These plots are generated for the specific value of x = 0.125 (and a range of y); however, at generic values of (x, y) solutions (when they exist) will take on either of these forms. On the right panel in Fig. 4, we can distinguish between the peak-like and massive behavior of Higgs (D − ) and symmetric (D + ) phase solutions. Furthermore, at x = 0.125, we observe that Γ simultaneously admits two solutions over a range of y, which is evidence of metastability. As shown and explained in Fig. 5, both solutions display gauge correlators with "massive" behavior (in the sense that lim p→0 G T (p) is finite; we are not claiming that the propagator has a pole at imaginary p and we have not investigated the behavior of spatial Wilson loops). In the symmetric phase this is due primarily to pure-glue loops; in the Higgs phase the mass is larger, due to additional Higgs-loop contributions.
On the left panel of Fig. 4 we see that the symmetric solution D + terminates. This indicates that the symmetric phase has lost its metastability and become spinodally unstable; so we identify the y value where this occurs as y − (x). The Higgs solution also terminates, and the possibility of metastability ceases to occur, at a larger y value which we interpret as y + (x). There is a third special value of y, where the Higgs solution becomes unstable to runaway behavior as described in Subsection 3.3. We will label this value y end . It does not have a physical interpretation in terms of the phase diagram; it is simply the point where the solution becomes so Higgs-like that our three-loop truncation encounters uncontrolled stability issues when we try to analyze the Higgs branch.
We can map out the region of the (x, y) plane between the y + and y − curves by finding those regions where two (meta)stable solutions for the propagators exist. The "symmetric" (small-D) solution is obtained by seeding the gradient solver with a configuration found at larger y, while the "broken" solution is found by starting at a smaller value of y with an initial guess for the scalar propagator with strong small-p behavior. The critical value x c is the largest x such that metastability is observed.
Plots of D(y) for x = 0.125 and x = 0.150 are shown in Fig. 6. In both cases we see a branch of symmetric phase solutions which terminates at y − . But while x = 0.125  supports a Higgs branch, x = 0.15 does not; so x c must occur between these two values. To determine where, we carry out a scan of the phase structure for several values of x, as shown in Fig. 7. The figure displays two-branch behavior at x = 0.14 but not x = 0.15, so we conclude that 0.14 < x c < 0.15.

Comparison with the lattice
Our original purpose in applying the 2PI formalism to SU(N ) Higgs theory was not specifically to determine the phase diagram (which is already known), but rather, to test the accuracy with which nPI resummation is able to make predictions about the nonperturbative sector of a nonabelian gauge theory. The nPI method relies on approximating the effective action by its truncation at a finite loop order, which results in a selective resummation to all orders of a certain class of topologies. In a gauge theory, this induces gauge-fixing dependence [48,49], since at least perturbatively, one should include all diagrams at every loop order. This effect could potentially be very mild, but a priori it is not clear that accurate results can be obtained from this method anywhere on the phase diagram. The only way to test the reliability of the approximation is to directly compute gauge-invariant observable quantities.
Here we will attempt a direct comparison between lattice and 2PI determined values of x c , y − (x) and y + (x). An overview of many of the pertinent results from 3D lattice studies of SU(2) Higgs theory can be found in [50], which incorporates the original studies [21][22][23]33]. The most relevant quantity to compare is the location of the critical endpoint. We find (x c , y c ) (0.145, 0.118). The accepted nonperturbative lattice value is (x c , y c ) = (0.0983±.0015, −.0175±.0013). There is qualitative agreement, but quantitatively the 2PI method has ∼ 50% relative errors in x c (establishing relative errors in y c is harder since it depended on an arbitrary renormalization point prescription).
We could also try to compare the spacing y + − y − to the lattice, at a comparable distance below x c . Unfortunately, the locations of the upper and lower metastability lines lack a clean nonperturbative definition. Technically, at any (x, y) value there is only one possible phase, and the transition line is where the is an abrupt change in that phase's properties, such as D. In practice, for systems near the transition line there are very longlived metastable states, and the transition from the metastable to the stable state involves an extremely rare and spatially inhomogeneous configuration. The spatial inhomogeneity is the reason that our 2PI approach cannot explore such states, allowing us to explore the supercooled or superheated phases. The lattice avoids this problem by sampling over all such states, typically using reweighting to make it more likely to sample the inhomogeneous states which carry us between metastable and stable phases. Nevertheless, Ref. [21] presented a definition of the metastability limits. For x = .0645 they find y + = .0009 and y − = −.0086, for a range of (y + − y − ) = .0095. This is comparable to the ranges we see in Fig. 7 for x = 0.10. So there is at least qualitative agreement here.
We are not able to compare the discontinuity between condensates at y c as a function of x c − x to the lattice, because we have not implemented a procedure to find the Γ difference between the two phases and thereby determine the transition value y c .

Comparison between Landau and Feynman gauges
Up to now, we have argued diagrammatically that critical values of y are expected to exhibit dependence on the gauge parameter ξ. However, since it is difficult to quantify this effect without an explicit computation, we will now briefly present a comparison between Landau and Feynman gauges. The results in Feynman gauge are best summarized by a ξ = 1 analogue of Fig. 7, shown in Fig. 8.
In setting ξ = 1 and resolving the SD equations (following the usual procedure), we observe that qualitatively very little has changed. Feynman gauge solutions exhibit similar features to those in Landau gauge, and once again we observe a disappearance of a stable Higgs branch somewhere between x = 0.125 and x = 0.150. This is consistent with the observation that dependence on x enters primarily through diagrams without gauge field lines. The biggest change though is the observed shift in the critical range of y, from around y ∼ 0.120 to y ∼ 0.250. This is interpreted as a contribution to the scalar mass from G L propagators, which does not fully cancel between diagrams. For instance, the mass contributions (at vanishing external momentum) of G L in the two self-energy corrections only cancel if the scalar propagator takes the free massless value 1/p 2 ; otherwise the tadpole contribution is larger and leads to a positive mass contribution which is proportional to the gauge parameter ξ.

Concluding remarks
We directly solved the three-loop 2PI effective action for 3D SU(N ) Higgs theory and obtained resummed correlators which correspond to both the symmetric and Higgs phases of the theory. We found that these solutions coexist over a region of the phase diagram, indicative of metastability and a first order phase transition. Subsequently, we have also observed that there is a point x where the metastability ceases to be observed, which we identified with the critical end point of the theory, x c . Concerning the numerical accuracy of the predictions made in Landau gauge, the location of the critical end point we inferred differs from the lattice value with a relative error of ∼ 50%. We also found that the critical y c depends surprisingly strongly on the gauge parameter ξ.
The most promising finding regarding the applicability of the nPI formalism to a nonabelian gauge theory is the apparent qualitative evidence for a critical end point (x c , y c ) located relatively close to its known nonperturbative value. In this sense the 2PI approach has successfully seen nonperturbative behavior in the phase diagram. However, the method has shown serious weaknesses as well. The quantitative level of agreement with the lattice is not very impressive, and the strong ξ dependence in y c is also worrying. More urgently, the method has failed completely to resolve the behavior of the Higgs phase when the scalar condensate is large. The most straightforward way to fix this problem, via the introduction of a scalar one-point function and the use of R ξ gauge, would introduce new problems. As we have argued, the ξ dependence should be significant where the transition is weak, and the gauge-fixing procedure may destroy the existence of a critical endpoint and analytic connection between the two "phases." Thus, the study of SU(N ) Higgs theory has therefore revealed several limitations to the nPI method in the context of a nonabelian gauge theory. In addition to the described ambiguities in physical observables, the application of the formalism is difficult numerically, especially if one wishes to consider higher-loop truncations or higher n-particle-irreducibility. If qualitative predictions can be made at best, then it may be hard to justify the numerical expense. However, this work does not preclude the possibility that further refinements may be possible with the goal of obtaining quantitatively accurate answers to nonperturbative and gauge-invariant questions. This matter is left open for a future investigation.

Acknowledgments
We would like to thank Meg Carrington and Marcus Tassler for useful comments. This work was supported in part by the Natural Science and Engineering Research Council of Canada (NSERC).

A Feynman rules for SU(N ) Higgs theory
The Feynman rules for covariant-gauge perturbative calculations in SU(N ) Higgs theory are derived from the Lagrangian Gauge field, scalar and ghost propagators are denoted by the symbols G, D and ∆. In Euclidean space at tree level these are where the gauge-field propagator is specified by the transverse and longitudinal projectors With all momenta assumed to be flowing outwards, the bare Yang-Mills vertices are where (a 1 , p 1 ) are the color indices and momentum of the outgoing ghost in V. The presence of a complex scalar results in the following additional vertices, λV (0)a 1 a 2 a 3 a 4 = −λ(δ a 1 a 2 δ a 3 a 4 + δ a 1 a 4 δ a 2 a 3 ). (A. 13) where the outgoing scalar(s) are indexed by (a 1 , p 1 ) (Eq. (A.11) and Eq. (A.12)) and a 1 , a 3 (Eq. (A.13)).

B Self-energies computed in dimensional regularization
In regularizing the 2PI effective action, one makes use of one-and two-loop self-energy corrections computed in perturbation theory. In pure Yang-Mills, all one-and two-loop integrals are massless from the onset. However, the inclusion of a Higgs field now in principle adds massive propagators to many of the diagrams. But, since we really only need to know the UV limit of these diagrams, it actually suffices to compute them with a massless scalar field. In this Appendix, though some results are valid for arbitrary D, should be treated as a small parameter, i.e., it is assumed that we are working at or near 3 dimensions, D = D 0 + 2 , with D 0 = 3. Finally, since the Higgs mass renormalizes at the two-loop level in three dimensions, it is useful to define the MS scaleμ 2 = µ 2 e γ /4π. The master one-loop topology is (B.2)

B.1 One-loop gluon self-energy
The presence of a scalar field adds two additional diagrams to the one-loop gluon self-energy relative to the the pure Yang-Mills expression, The result is strictly transverse; we will separate the Yang-Mills and Higgs contributions as follows, The terms which appear in the limit D → 3 are and it is also useful to take the m → 0 limit and keep terms O( ),

B.2 One-loop Higgs self-energy
The calculation of the one-loop correction of the Higgs self-energy proceeds forward in much the same manner, For the D → 3 and massless limits we have in terms of the dimensionless quartic coupling x = λ/g 2 .

B.3 One-loop ghost self-energy
The one-loop ghost self-energy is constructed out of a a single diagram, for which in D = 3 + 2 , ξ dependence only appears at O( ),

B.4 Two-loop topologies
The massless two-loop master topology is using a notation where the subscript zero refers to the mass of the scalars in the loops being set to m 2 = 0. As mentioned at the start of this appendix, for the purpose of regularizing this calculation genuinely two-loop topologies can be computed in the massless limit. However, recursively one-loop diagrams (labeled with the superscript IR2) will exhibit IR divergences without the inclusion of a regulator mass ω. We have (retaining the superscript IR to indicate that the full expression involves the specifically IR regulated diagrams)  Now, regarding the notation: at this point there are two quantities which can be regarded as masses, m 2 and ω 2 . m 2 refers to the Higgs mass which enters the problem via the scalar propagator, which we have already set to zero. Whereas, ω 2 is an unphysical regulator mass introduced to regulate IR divergences in some two-loop diagrams. So, for instance, the diagrams which comprise π (IR2, ) YM;0;µν are calculated using finite ω 2 , but setting m 2 = 0. One may ask why we do not simply regulate the IR divergences by keeping the scalar field massive from the onset? There are two reasons. First, a number of divergences arise from a 1/p 3 gauge field propagator, so this would not solve the problem entirely. Second, in general these diagrams are introduced to regularize the UV divergences in the problem. To compute the leading order UV behavior, it is sufficient to set m 2 = 0, which drastically simplifies the majority of the diagrams which must be calculated. Then, the IR divergences which would arise in the bare perturbation theory are handled with ω 2 , of which the final results will be independent regardless.