{'fraction_non_alphanumeric': 0.03652341924023809, 'fraction_numerical': 0.011592757642805946, 'mean_word_length': 4.3635300016714025, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Calculation of the pressure of a hot scalar theory within the Non-Perturbative Renormalization Group

We apply to the calculation of the pressure of a hot scalar field theory a method that has been recently developed to solve the Non-Perturbative Renormalization Group. This method yields an accurate determination of the momentum dependence of n-point functions over the entire momentum range, from the low momentum, possibly critical, region up to the perturbative, high momentum region. It has therefore the potential to account well for the contributions of modes of all wavelengths to the thermodynamical functions, as well as for the effects of the mixing of quasiparticles with multi-particle states. We compare the thermodynamical functions obtained with this method to those of the so-called Local Potential Approximation, and we find extremely small corrections. This result points to the robustness of the quasiparticle picture in this system. It also demonstrates the stability of the overall approximation scheme, and this up to the largest values of the coupling constant that can be used in a scalar theory in 3+1 dimensions. This is in sharp contrast to perturbation theory which shows no sign of convergence, up to the highest orders that have been recently calculated.


I. INTRODUCTION
There has been much effort devoted in the recent years to the development of finite temperature field theory, in particular in the context of Quantum Chromodynamics (QCD) whose thermal properties are central to the understanding of heavy ion collisions at high energy. Equilibrium properties of hot QCD are calculable from lattice gauge theory (for a review, see e.g. [1]), but there is a need to develop semi-analytical tools to understand the results of such calculations, with the hope that such tools may also allow one to approach non-equilibrium situations.
Weak coupling expansions are among such tools. In the case of QCD, the use of perturbation theory is motivated by the asymptotic freedom that leads to a small effective coupling at high temperature. However, strict perturbation theory does not work at finite temperature: it exhibits indeed very poor convergence properties, even in a range of values of the coupling where good results are obtained at T = 0. This difference of behavior of perturbation theory at zero and finite temperature can be understood from the fact that, at finite temperature, the expansion parameter involves both the coupling and the magnitude of thermal fluctuations (for a recent review, see [2]; see also [3]). In that respect, the problem is not specific to QCD: Similar poor convergence behavior appears also in the simpler scalar field theory [4], and has also been observed in the case of large-N φ 4 theory [5]. Similar observations can be made also in the case of QED [6].
In this paper we restrict ourselves to the case of the simple theory of a scalar field ϕ, with a quartic self-coupling g 2 ϕ 4 , for which high order calculations were recently completed. Thus the pressure is known to order g 8 ln g [7], and the screening mass to order g 6 ln g [8]. Such calculations were made by exploiting effective theory techniques, in particular dimensional reduction, that rely on a separation of scales of various degrees of freedom in the hot scalar plasma: hard modes with momenta k ∼ T that contribute dominantly to the pressure and are weakly coupled, and soft modes with momenta k ∼ gT that are more strongly coupled.
In QCD, another scale concerns the ultrasoft modes, with momenta k ∼ g 2 T that remain strongly coupled for any value of the coupling constant. This scale, relevant only in the case where massless modes exist at finite temperature, does not play any role in the scalar field theory. The separation of scales that allows the organization of the calculation using effective field theory disappears when the coupling is not too small: then the various degrees of freedom mix and the situation requires a different type of analysis. The purpose of this paper is to apply to this problem the Non-Perturbative Renormalization Group (NPRG) (for reviews of this method see e.g. [9][10][11]).
To do so, we shall rely on an elaborate approximation scheme that has been developed recently in order to obtain a good determination of the momentum dependence of the npoint functions [12]. This new approximation has been tested on the O(N) models, for which it provides excellent critical exponents, and more generally, an excellent description of the momentum dependence of the 2-point function, from the low momenta of the critical region, all the way up to the large momenta of the perturbative regime [13]. One may then expect this method to capture accurately the contributions to the thermodynamical functions of thermal fluctuations from various momentum ranges, and hence handle properly the mixing between degrees of freedom that takes place as the coupling grows. Since it involves also non trivial momentum dependent self-energies, the method also encompasses effects related to the damping of quasiparticles, or their coupling to complex multi-particle configurations.
The present paper may be viewed, in its spirit and goals, as a follow up of the analysis presented in Ref. [14]. However it departs from it in two ways. The first difference is of a technical nature: the new approximation scheme is better justified when one uses an Euclidean symmetric four dimensional regulator that cuts off the contribution of high Matsubara frequencies. In contrast, in Ref. [14], we used a three dimensional regulator, and performed analytically the (untruncated) sums over the Matsubara frequencies. The second difference is that we use a new, much more accurate, approximation scheme to solve the NPRG equations, as mentioned above: The Local Potential Approximation (LPA) used in Ref. [14] can be viewed as the zeroth order in this new approximation scheme. In physical terms, the LPA corresponds to an approximation where the degrees of freedom of the hot scalar plasma are massive quasiparticles. The new scheme goes beyond that simple picture.
As it turns out, the results obtained are not too different from those of the LPA. This stability of the results against improvements in the approximation suggests that the scheme that we are using to solve the NPRG equations may give already, at the level at which it is implemented here, an accurate representation of the exact pressure, and this over a wide range of coupling constants. It also indicates that for such a system the quasiparticle picture is presumably robust.
The outline of this paper is as follows. In Sec. II, we present a brief introduction to the NPRG and the approximation scheme of Ref. [12], indicating specific features of finite temperature calculations. In Sec. III we integrate the flow equations numerically and discuss the results obtained. The last section summarizes the conclusions. In the Appendices we give details about the numerical integration needed to calculate the flow of the pressure, and we discuss specific features of the exponential regulator used in our calculations at finite temperature.

II. THE NPRG IN THE BMW APPROXIMATION SCHEME
We consider a scalar field theory with the classical (Euclidean) action where T is the temperature. Our goal is to calculate the thermodynamical pressure for this scalar theory, using the NPRG. We follow here Ref. [10], and add to the original action S a where the parameter κ runs continuously from the microscopic scale Λ down to 0. We use the notation with q = (ω n , q), q 2 = ω 2 n + q 2 , and ω n = 2πnT are the Matsubara frequencies. The role of ∆S κ is to suppress the fluctuations with momenta |q| κ, while leaving unaffacted those with |q| κ (|q| = q 2 ). This is achieved with a cut-off function R κ (q) that has the following properties: R κ (q) → κ 2 when |q| ≪ κ, and R κ (q) → 0 when |q| ≫ κ. The precise form of the function R κ (q) used in our calculation will be discussed later.
The effective action Γ κ [φ] associated to the action S + ∆S κ obeys the exact flow equation [15] (with ∂ t ≡ κ∂ κ ): where G κ [q, −q; φ] is the full propagator in the presence of the background field φ: of the scalar field theory is obtained as the solution of Eq. (2.4) for κ → 0, at which point When φ is constant, the functional Γ κ [φ] reduces, to within a volume factor, to the effective potential V κ (φ). The flow equation for V κ follows from that of the effective action Γ κ , Eq. (2.4), when restricted to a constant φ. It reads We used here the simplified notation Γ (2) κ (q, ρ) in place of Γ (2) κ (q, −q, ρ) for the 2−point function in a constant background field, and similarly for G(q, ρ). Also, we have set ρ ≡ φ 2 /2, a notation to be used throughout (when φ is constant). The pressure P is related to the effective potential by (2.8) The equation for the effective potential may be viewed as the equation for the "zeropoint" function in a constant background field. By taking two derivatives with respect to φ and letting φ be constant, one obtains the equation for the 2-point function in a constant background field: (2.9) In this equation, all the n-point functions depend on the constant background field φ. This has not been indicated explicitly in order to alleviate the notation.
The flow equations (2.6) and (2.9) are the first equations of an infinite tower of coupled equations for the n-point functions, whose solution requires some truncation. We shall use the truncation scheme proposed recently by Blaizot, Méndez-Galain and Wschebor (BMW) [12], implemented here at its lowest non-trivial order: in this case we need only consider the equation for the effective potential and that for the 2-point function, that is, Eqs. (2.6) and (2.9), on which we shall perform approximations that are described next. For a more complete discussion we refer to Ref. [16].
The lowest level of the approximation scheme corresponds to a widely used approximation, usually referred to as the Local Potential Approximation (LPA). It consists in assuming that for all values of κ the effective action takes the form [10] which is tantamount to assume that the 2-point function is of the form With this Ansatz, Eq. (2.6) for V κ becomes a closed equation that we write as follows where I 1 is one of the following integrals Note that Eq. (2.12) is formally an exact equation, and it would yield the exact effective potential if I 1 were calculated with the exact propagator. In the LPA, the equation keeps the same form, but the propagator is given by Eq. (2.11), where m 2 κ (ρ) is itself determined by the potential (thereby making Eq. (2.12) a closed, self-consistent equation).
The next order of the approximation scheme, to which we refer to as the leading order of the BMW method [12], or here simply as BMW for briefness, consists in neglecting the loop momentum q in the 3 and 4-point functions in the right hand side of Eq. (2.9). Once this approximation is made, the corresponding n-point functions can be obtained as the derivatives of the 2-point function with respect to the constant background field: 14) The equation for the 2-point function becomes then a closed equation There is however a subtlety: this equation for the 2-point function is coupled to that of the effective potential, because Γ . In order to properly implement this coupling, we treat separately the zero momentum (p = 0) and the non-zero momentum (p = 0) sectors, and define is obtained by solving the equation for the effective potential. The equation for ∆ κ (p, ρ) can be easily deduced from that for Γ (2) , i.e., from Eq. (2.15) by subtracting the corresponding equation that holds for p = 0. It reads where the symbol ′ denotes the derivative with respect to ρ, and we have set u κ (ρ) ≡ ∂m 2 κ (ρ)/∂ρ. This equation (2.17), together with Eq. (2.6) for the effective potential and that for the constitute a closed system of equations for ∆ κ (p, ρ) and V κ (φ). This can be solved with the initial condition Γ Λ (p; ρ) = p 2 + m 2 + uρ, where m 2 and u are essentially (to within small ultraviolet cut-off corrections) the parameters of the action (2.1).
We now specify the regulator R κ (q) that we have used in our calculation. We take it of the generic (Euclidean symmetric) form where Z κ is a function of κ only, to be specified shortly, and the function r(q) is a smooth function of its argument. In most of our calculations, we have used an exponential regulator of the form where α is a free parameter. We have also considered an alternative regulator, given by with parameters β = 1/2 and γ = 1/24 chosen so that the Taylor expansion aroundq = 0 agrees with the regulator (2.20) through order O(q 4 ), and leave only the prefactor α as a free parameter 1 . The regulators (2.20) and (2.21) are suited for the BMW approximation whose justification relies on both momenta and frequencies being cut-off. For comparison, we have also solved the LPA with these exponential regulators. We shall also compare to the LPA results of Ref. [14] where the Litim regulator r(q) = (1 −q 2 )θ(1 −q 2 ) [17] was used in integrals over three momenta, the Matsubara frequencies being integrated over analytically.
In the absence of any approximation, the physical quantities such as the pressure should be strictly independent of the cut-off function, and in particular of the value of the parameter α that we have introduced in the regulators (2.20) and (2.21). In practice, we find a (weak) residual dependence on the parameter α, and a study of this spurious dependence provides an indication of the quality of the approximation. It should be emphasized however that this dependence on α is rather small, as we shall see, and we have made no effort to perform a systematic study of the dependence of our results on its value. Nor did we explore the effects of enlarging the space of possible variations by allowing for instance the parameters β and γ in (2.21) to take arbitrary values.
The factor Z κ in Eq. (2.19) reflects the finite change in normalization of the field between the scale Λ and the scale κ. It is defined by where ρ 0 is a priori arbitrary but chosen here to correspond to the minimum of the effective potential. This factor Z κ enters also the definition of the dimensionless variables that are used in the numerical solution. Thus, for instance we define with m κ ≡ m k (ρ = 0), u κ ≡ u κ (ρ), and K ≡ 1/(32π 2 ). In the next section, when we discuss the results obtained by solving numerically Eqs. (2.17) and (2.6), we shall set g 2 κ ≡ u κ /24.

III. RESULTS
To calculate temperature dependent physical quantities, we need to evaluate the flow of these quantities at both zero temperature and at finite temperature. We follow here the strategy exposed in Refs. [14,18]. Conceptually, one starts with given physical parameters at zero temperature at the scale κ = 0. One then removes quantum fluctuations step by step by integrating the flow equations upwards from κ = 0 to Λ in order to arrive at "bare quantities" at a chosen scale Λ. If Λ is chosen big enough, only renormalizable local operators survive in the effective action, and the system can then be described by a simple set of bare parameters g Λ , m Λ , and Z Λ . Starting from these bare parameters one then follows the flow downwards from κ = Λ to 0, but this time with the temperature T turned on. The physical quantities are then obtained at κ = 0.
In practice, for reasons of numerical stability, the flow is never integrated upward, but always downward, for both zero and finite temperature. For a given bare coupling g Λ (and Z Λ = 1), one adjusts the bare mass m Λ at the scale κ = Λ by bisection, so as to arrive at a vanishing mass m 0 = 0 at the end (κ = 0) of the zero temperature flow, to the desired accuracy. Keeping the same bare parameters, one then turns on the temperature and runs the flow again. Note that we consider specifically here a massless theory in order to be able to compare our results with the known results of high order perturbation theory. But the same analysis could be done for any value of m 0 .
Since the zero-temperature coupling constant vanishes at κ = 0 for any g Λ (the so-called triviality of ϕ 4 theory in d = 4), one adjusts the coupling constant at a finite scale κ > 0.
In accordance with what is commonly done in perturbation theory, we choose to do that at κ = 2πT on the T = 0 flow. This procedure induces a specific scheme dependence, attached to the choice of the regulator, which should be taken into consideration when comparing with results that are obtained in another scheme, for example the minimal subtraction scheme in perturbation theory.
In order to keep the zero-temperature flow under control, we introduce dimensionless variables [10]. Although these are not optimal for the thermal flow, which freezes in dimensionful variables (but diverges in dimensionless variables) as κ → 0, it is still advantageous to use the same scaling of variables as for the zero-temperature flow in order to achieve high-precision cancellations between zero and finite temperature flows for larger κ ≫ T . We have adapted to finite temperature the numerical strategy that has been developed to solve the BMW equation at zero temperature in the context of O(N) models [13,16].
This method puts propagator and potential on a regular two-dimensional grid in |q| and ρ variables and solves the flow equations using the Euler method. At finite temperature, the flow equation (2.6) involves a summation over Matsubara frequencies. This summation is infeasible in practice at large values of κ, but at sufficiently large κ ≫ T , the corresponding vacuum integral constitutes the dominant contribution and any thermal contribution is exponentially suppressed. Practically, a switching temperature T switch ≫ T is introduced such that, at κ = T switch , the numerical code switches from vacuum integration to thermal integration with Matsubara frequencies. At this point, the values of the vacuum propagator G κ (q, ρ) which are stored on a two-dimensional grid (|q|, ρ) have to be interpolated in order to yield values on a three-dimensional grid (ω, |q|, ρ) using cubic splines. Good results have been achieved on a n q ×n ρ = 90×50 grid withq max = 9 andρ = 5 (in dimensionless variables, see Eq. (2.23)) with 30 Matsubara terms at the switching temperature T switch which is varied in the range 4πT ≤ T switch ≤ 6πT . To speed up the calculation, the number of Matsubara terms can be gradually reduced as the flow proceeds, without affecting accuracy. Further details on the integration of the pressure flow are given in Appendix A. Figure 1 shows the flow of various quantities (coupling constant, mass, pressure and Zfactor) in the BMW approximation. (The pressure P 0 is that of the non-interacting scalar plasma, i.e., P 0 = π 2 T 4 /90.) These flows are remarkably similar to the corresponding ones obtained within the LPA approximation in Ref. [14]. At finite temperature, the flow starts to deviate from the vacuum flow between 2πT and T , and stabilizes shortly below the scale gT (with g ≡ g(2πT )). Note that gT is the leading order value at small g of the thermal mass of the excitations (see e.g. Fig. 2 system. This could be confirmed by a complete calculation of the single particle spectral function, which is beyond the scope of this paper. BMW approximations to the NPRG are compared to perturbation theory through order g 6 log g for the mass [8] and g 8 log g for the pressure [7], and also to the result of the 2-loop 2PI resummation from Ref. [14]. The BMW results were all obtained with the exponential regulator (2.20), while the LPA results were obtained with this same exponential regulator, as well as with the 3-dimensional Litim regulator used in Ref. [14]. For a given choice of regulator, one sees in Fig. 2 that the difference between the LPA and the BMW results for both the mass and the pressure is tiny: one does not gain much, for the thermodynamics, in improving the treatment of long wavelength degrees of freedom and incorporating explicit momentum dependence in the self-energy. The stability of the BMW approximation scheme, comparable to that of the 2PI resummation, should be contrasted with the wild oscillatory behavior of the successive orders of perturbation theory.
The plots in Fig. 2 depend on the choice of the prescription for obtaining the coupling at a particular scale (in this case at the scale µ = 2πT ), and hence on the regulator: depending on the regulator, a given value of the coupling constant g(2πT ) will correspond to different bare coupling constants g Λ , and hence to different initial conditions for the temperature flow.
Together with the obvious dependence of the flow itself on the regulator this will affect the final result. The regulator dependence that is visible in the LPA results in Fig. 2 is to be attributed in part to the way the results are plotted, namely as a function of g(2πT ). Most of such a "scheme dependence" can be eliminated by plotting only physical quantities. This is realized in Fig. 3 which shows the pressure as a function of the mass. In this combination, only physical quantities are compared, that do no longer depend on a particular choice of a scale at which one fixes the coupling g, as was the case in Fig. 2. While perturbation theory clearly breaks down above a certain ratio m/T 0.3, RG methods and the 2PI approach give consistent results up to twice this value. As noticed previously [14], the perturbative g 2 contribution seems to be a surprisingly good approximation for the behavior at larger m/T -but only for the plotted ratio (P/P 0 ) versus (m/T ): the g 2 curve representing P/P 0 in Fig. 3 is identical to the corresponding curve in Fig. 2 (right panel), since, as already mentioned, at leading order in the coupling constant, m/T = g. Thus the reason the curves corresponding to the various non-perturbative approximations are brought closer to the g 2 curve in Fig. 3 can be traced back to the fact that the mass increases much less rapidly than g with increasing g, as clearly visible in Fig. 2 (left panel). Note that the dependence of the LPA results on the choice of the regulator remains present, although it is less important than in Fig. 2. We shall return to this issue shortly.
We turn now to more technical aspects of the calculations, namely the dependence of the results displayed above on the choice of the temperature, or on the choice of the regulator.
Consider first the dependence on the temperature, which is measured by its ratio T /Λ to the microscopic scale Λ. As T /Λ increases and becomes close to 1, our numerical calculations loose accuracy for a variety of reasons, the main one being the following: If T is too close to Λ, 2πT may become bigger than Λ and the whole procedure eventually collapses. One indeed assumes that the beginning of the flow is not affected by the temperature (so as to use the 4-dimensional integration procedures), and this assumes Λ > 2πT so that there is room for a 4-dimensional flow between Λ and 2πT , with κ ≈ 2πT the point where thermal fluctuations start to contribute to the flow. This limitation is illustrated in Fig. 4, for the sizable dependence starts to be visible. Note that this is the value of the mass where we observed possible cut-off effects when the temperature is too large. In the present case, the temperature is not to blame. However cut-off effects may show up in the fact that integrals are done with Λ as ultraviolet cut-off. In most cases this is redundant and does not interfere with the regulator (whose derivative with respect to κ also provides an ultraviolet cut-off).
However in the mass region m/T ≈ 0.5 ÷ 0.6, the two cut-offs interfere, and produce a sensitivity of the results to the regulator: indeed a larger α allows larger momenta.
To investigate this further, we have considered the other regulator mentioned in Sect. II, namely, Eq. (2.21). This regulator has the property to cut-off more efficiently the large momenta (because of the presence of the q 4 term). And indeed this is what one sees in the results on the temperature T switch at which one changes from the 4-dimensional to the 3-dimensional flow: essentially no such dependence is observed with the second regulator as long as T switch > ∼ 4π. One should emphasize however that the difference between the two regulators is tiny. Thus for instance, the plots in Fig. 6 show the derivative of the pressure as a function of κ. The oscillatory behavior at the beginning of the flow is generic for regulators that cut-off the sum over the Matsubara frequencies, and the particular oscillations exhibited in Fig. 6 may be understood from the analytical structure of the regulator analyzed in detail in Appendix B. One sees clearly that the larger α the larger are the oscillations. On the other hand, there is hardly any visible difference between the two panels corresponding to the two different regulators. Only if one zooms in a lot, can one see that (2.21) is slightly better in cutting off ∂P κ at larger values of κ than (2.20). The plots show curves for g Λ = 2.8, corresponding to g(2πT ) = 1.1 and T = 10 −2 Λ (that is, the same parameters as the right column of Fig. 1).

IV. CONCLUSIONS
In this paper we have applied a powerful approximation scheme to solve the NPRG equations and we have calculated various thermodynamical quantities for a scalar field theory at finite temperature, in a range of values of the coupling constant covering both weak and strong coupling regimes. Of course, because of the presence of the Landau pole, arbitrarily large values of the coupling cannot be reached. However the range of values of g that can be explored allowed us to demonstrate the stability of the results as successive orders in the approximation scheme are taken into account, in sharp contrast to the strict expansion in terms of the coupling constant. A similar stability also emerges in other non-perturbative schemes, such as "screened perturbation theory" [19], or the 2PI effective action [20] (see also [21]).
In fact, the results that we have obtained using the BMW approximation differ very little from those that we obtained previously using the LPA [14]. This may be attributed to the fact that the thermal mass provides an infrared cut-off that reduces the contributions to the pressure of the long wavelength modes, which are the most strongly coupled. The thermal mass also provides a threshold that hinders the mixing with complex multiparticle configurations, making the quasiparticles well defined. This latter aspect is explicitly verified by the small deviation of the field normalization from unity, as obtained within the BMW approximation. These two effects conspire to give the approximation scheme a remarkable stability, and contribute to the robustness of the quasiparticle picture.
We have also compared the results obtained within the BMW scheme with those of the simple 2-loop 2PI resummation used in [14]: both methods lead to very similar results in the extrapolation to strong coupling. This is not too surprising, given the agreement already observed between the LPA and the 2PI method in Ref. [14]. The stability of the 2PI scheme itself can be assessed from the 3-loop calculation of Ref. [20]. The detailed comparison between the latter calculation and ours is not straightforward however because of the different renormalization schemes used in the two cases. However, the main message is essentially the same: the main qualitative difference between the 2-loop and 3-loop calculations is the presence of momentum-dependent self-energies in the latter, in contrast to simple mass terms in the former. The small difference observed between the results in the two cases, corroborates the conclusions that follow from our NPRG analysis about the robustness of the quasiparticle picture. In fact the stability of the results lead us to conjecture that any corrections to them are presumably very small. It would be interesting to have lattice calculations allowing us to test this conjecture for values of the couplings where perturbation theory breaks down.

It is of the form
, where the vacuum (T = 0) integral with vacuum propagator is subtracted from the thermal sum with thermal propagator: As can be seen in Fig. 7, the two integration domains do not match. For large cutoffs Λ, this mismatch may give rise to spurious contributions to the pressure.
One solution could seem to be to restrict the range of the Matsubara sum to a 4dimensional sphere. This turns out to be problematic for the following reason: Good convergence properties are only obtained if sufficiently many Matsubara terms are summed up, but at large values of q, only few Matsubara terms would contribute if the restriction to a 4-dimensional sphere were imposed.
It is better then to fix the mismatch in an alternative way, i.e., by extending the vacuum integration domain so that it matches the domain of the Matsubara summation. That is, we write for q 0,max = 2πn max T . Practically, only the pressure is sensitive to this correction. At each integration step in t direction, the vacuum integral is calculated twice: Once as in Eq. (A1) to follow the vacuum flow, and additionally as in Eq. (A2) to obtain the contribution for the pressure. The values of G T =0 κ (q, ρ) that are only known on a grid in (|q|, ρ) coordinates have to be interpolated to obtain values on a three-dimensional grid (q 0 , |q|, ρ). Since for each temperature the switching point between 4D vacuum integration and 3D Matsubara summation varies, the vacuum piece for Eq. (A2) has to be calculated for each temperature separately.
Appendix B: Oscillatory behavior of the flow As we have seen in Sect. III, the flow of the pressure exhibits an oscillatory behavior at the beginning of the flow. This generically occurs whenever the sum over the Matsubara frequencies is limited to a finite number of terms. We show in this Appendix that the oscillations seen in Fig. 6 at the beginning of the flow can be understood from the analytic structure of the particular regulator that we are using, and can be simply calculated for small values of κ. The analysis is performed within the LPA with the exponential regulator (2.20) with α = 1. (Note that much wilder oscillations than the one discussed here are observed when one uses a "hard" regulator, such as the Litim regulator [22].) Instead of performing the sum over Matsubara frequencies in Eq. (2.6) explicitly, we can convert it into a contour integration in the following way: The resulting integrals can then be calculated by closing the vertical contours with semi circles of infinite radii, and summing over the encircled residues. It is convenient to group the poles of the integrand in Eq. (2.6) into two categories (see Fig. 8). There are the poles associated with the regulator, which we shall call trivial poles, and those of the propagator, that is the zeros, of G −1 κ (q, ρ) = q 2 + m 2 κ + R κ (q).
It turns out that the non-trivial poles lie to the right of the trivial poles in the complex plane. Due to the statistical factor in Eq. (B1), the corresponding residues are therefore exponentially suppressed. In fact, the first few trivial poles are sufficient to accurately approximate Eq. (2.6) for large values of κ.
The oscillatory behavior follows from the distribution function 1/(e q 0 /T − 1) with complex number q 0 .
The thermal contribution is directly given by the second line of (B1). One can write the contribution of the trivial poles as (instead of adding the residues above and below the real axis, we simply can take two times the real value of one of the residues; another factor 2 comes from f (q 0 ) + f (−q 0 ); the factor 4π = dΩ) 1 e q 0 (n)/T − 1 q 2 −q 2 0 (n) + q 2 q 0 (n) (B9) or explicitly: For large κ this gives a good approximation to the thermal pressure contribution, as can be seen in Fig. 9. Taking only the first term n = 1 already gives a good approximation for large values of κ.