Self-consistent $O(4)$ model spectral functions from analytically continued FRG flows

In this paper we explore practicable ways for self-consistent calculations of spectral functions from analytically continued functional renormalization group (aFRG) flow equations. As a particularly straightforward one we propose to include parametrizations of self-energies based on explicit analytic one-loop expressions. To exemplify this scheme we calculate the spectral functions of pion and sigma meson of the $O(4)$ model at vanishing temperature in the broken phase. Comparing the results with those from previous aFRG calculations, we explicitly demonstrate how self-consistency at all momenta fixes the tight relation between particle masses and decay thresholds. In addition, the two-point functions from our new semi-analytic FRG scheme have the desired domain of holomorphy built in and can readily be studied in the entire cut-complex frequency plane, on physical as well as other Riemann sheets. This is very illustrative and allows, for example, to trace the flow of the resonance pole of the sigma meson across an unphysical sheet. In order to assess the limitations due to the underlying one-loop structure, we also introduce a fully self-consistent numerical scheme based on spectral representations with scale-dependent spectral functions. The most notable improvement of this numerically involved calculation is that it describes the three-particle resonance decay of an off-shell pion, $\pi^* \to \sigma\pi\to3\pi$. Apart from this further conceptual improvement, overall agreement with the results from the considerably simpler semi-analytic one-loop scheme is very encouraging, however. The latter can therefore provide a sound and practicable basis for self-consistent calculations of spectral functions in more realistic effective theories for warm and dense matter.


I. INTRODUCTION
The properties of strong-interaction matter, i.e. matter subject to the strong force described by Quantum Chromodynamics (QCD), are in the focus of intense experimental and theoretical study worldwide. In particular at finite temperature and high baryon density, where one searches for a QCD critical endpoint [1,2], this is a very challenging task. In the region of the QCD phase diagram relevant for the equation of state in the interior of neutron stars or their merger events neither lattice QCD, perturbation theory, nor chiral perturbation theory can be applied for different reasons without restrictions. In this regime, functional approaches like the Functional Renormalization Group (FRG) or Dyson-Schwinger Equations (DSEs) might help because they can in principle be applied at any temperatures and baryon chemical potential, see e.g. [3,4] for recent compilations of results on the QCD phase diagram. Especially the FRG approach is particularly well suited to describe equilibrium thermodynamics and phase transitions. In order to avoid the peculiarities associated with gauge invariance and confinement in QCD, one thereby often employs universality and structurally simpler replacement theories to study static and dynamic critical phenomena.
However, these functional approaches are originally also formulated in Euclidean spacetime which hampers the access to real-time quantities like spectral functions or transport coefficients due to the analytic continuation problem. In relativistic theories this entails a transition from Euclidean to Minkowski spacetime, for various reconstruction methods to achieve that see for example [5][6][7][8][9][10][11]. This is another challenging problem and, when based on a finite number of discrete data points with errors, even an ill-posed numerical problem.
For ab-initio simulations of the static and dynamic universal critical behavior of spectral functions near continuous thermal phase transitions, one can resort to classical-statistical field theory [12][13][14][15]. Reconstructing spectral functions from Euclidean data can in principle also be avoided in non-perturbative functional methods, e.g. based on the 2-PI effective action [16][17][18][19] or Dyson-Schwinger equations [20,21]. Within the FRG, several techniques have been developed recently that are based on a direct computation of correlation functions in Minkowski spacetime and hence avoid the need for any numerical reconstruction method [22][23][24][25][26][27].
The first FRG calculation of the O(4) model spectral functions for pions and the sigma meson that was based on a direct analytic continuation of the Euclidean FRG flow equations for the corresponding two-point functions was presented in [28]. The calculation was set up in a thermodynamically consistent way, in that the flow of the two-point functions for vanishing external momenta, at fixed field values agreed with the corresponding flow for the curvature of the effective potential at the same values. However, the pion and sigma propagators inside these analytically continued FRG (aFRG) flows were still the simple quasi-particle constituent propagators of the leading order derivative expansion, i.e. the local potential ap-proximation (LPA), that was used to calculate the effective potential. This means that the two-point functions inside the flows only contained the scale-dependent quasiparticle masses but no momentum-dependent self-energy corrections. As a result of this missing self-consistency beyond the zero-momentum limit, a certain mismatch arose between the physical pole masses and the location of decay thresholds in the calculated spectral functions, for example.
Some of these shortcomings were overcome in [29] where a self-consistent approach for the calculation of real-time correlation functions was proposed and applied to the O(4) model in the vacuum. Therein, the analytic continuation was performed as a purely numerical procedure put forward in [30]. This procedure was, however, found to be very computation time intensive and results for the spectral functions were only obtained up to rather low energies. In addition, access to the analytic structure of the real-time correlation functions was limited. The numerical analytic continuation procedure of [30] was extended to finite temperature and finite spatial momentum in [31] where the spectral functions were obtained for a larger range of energies, although not in a fully self-consistent way.
In this work we further develop and study selfconsistent approaches to calculate real-time correlation and spectral functions. A particularly straightforward one is based on the idea [23] of including parametrized self-energies with the required analyticity properties directly in the ansatz for the effective average action and calculate the flows of all parameters. The difficulty then is to find such parametrizations with the correct domain of holomorphy for two-point functions in local quantum field theory which is not entirely straightforward in general. Here, we achieve this by simply using explicit analytic expressions of one-loop form to parametrize the self-energies of pion and sigma meson at zero temperature with broken O(4) symmetry. Because the resulting parametrizations are known analytically, on all Riemann sheets in fact, the analytic continuation from Euclidean to Minkowski spacetime is then straightforward again, and gives immediate access to the spectral functions as well as to the real-time propagators at arbitrary energies. At the same time the analytic expressions for the parametrized self-energies are also used in the propagators inside the flow of the effective potential which thus includes the same non-trivial momentum dependence of the two-point functions as encoded in the spectral functions in a self-consistent manner.
Comparing our results for the O(4) model spectral functions from this new self-consistent scheme with those from a corresponding aFRG calculation that is built on the LPA truncation as in most previous aFRG applications, see e.g. [32][33][34][35][36][37][38][39][40], we explicitly demonstrate that the mismatch between pole masses and decay thresholds is thereby removed. In addition, the effects of the momentum-dependent self energies are now selfconsistently included in the calculation of the effective potential. Moreover, the analytic structure of the propagators is known which gives access to the spectral functions as well as to the propagators in the cut-complex frequency plane, on the physical and other Riemann sheets.
After having established the merits of these selfconsistent flows with parametrized self-energies, called the 'SC1L' flows below, we also assess their limitations. These clearly are inherited from the explicit one-loop forms of the self-energies used here, which limit the structure of the possible imaginary parts from Cutkosky's cutting rules. We therefore also perform fully self-consistent aFRG calculations for comparison, with numerically computed self-energies included. This is possible based on spectral representations with scale-dependent spectral functions that are recomputed with every integration step in the FRG flow as we will explain. As a proof of principle, we demonstrate that we are able for the first time to describe the three particle resonance decay contribution to the pion spectral function, of an offshell pion into three pions via a broad sigma resonance, at the correct threshold of three times the pion mass. This clearly comes from an imaginary part beyond oneloop. At the same time, the overall modifications in the fully self-consistent numerical calculation, as compared to the semi-analytic SC1L flows, are otherwise surprisingly small. This paper is organized as follows: In Sec. II we start with a brief introduction of the general FRG concepts needed here. In Sec. II A we summarize the analytically continued LPA setup which will be used for comparison with the self-consistent SC1L flows introduced in Sec. II B. In Sec. II C we discuss parameter sets, UV conditions, and our numerical implementation. In Sec. III we compare the results for effective potential, particle masses, and the coupling constants from LPA and SC1L flows. The corresponding results for the two-point functions in Euclidean and Minkowski spacetime, from analytic continuation, together with the spectral functions are compared in Sec. IV. In Sec. V we discuss the analytic structure of the SC1L two-point functions in the complex-momentum plane. As a first step to include selfenergies beyond one-loop structure we describe a fully self-consistent numerical framework, based on spectral representations with scale dependent spectral functions, in Sec. VI, where we also compare the fully self-consistent with the SC1L results from the previous sections to assess the overall quality of the latter. We close with our summary and an outlook in Sec. VII. Details on the explicit analytic expressions for the self-energy parametrizations and loop functions are provided together with other technicalities in several appendices.

II. THEORETICAL SETUP
The FRG is a powerful and versatile non-perturbative framework that aims at calculating the effective action by successively integrating out quantum fluctuations in momentum space, for reviews see for example [41][42][43][44][45][46][47][48][49]. In the formulation pioneered by C. Wetterich [50], the RG-scale dependence of the effective average action Γ k , which interpolates between the bare action S Γ Λ at an ultraviolet (UV) scale Λ and the full quantum effective action Γ ≡ Γ k=0 in the infrared (IR), is given by the following flow equation, Therein, the regulator function R k acts as a mass term and suppresses fluctuations of low-momentum modes, p k, while the high-momentum modes, p k, are already integrated out and included in Γ k . For a discussion of how to devise optimized regulators in a particular truncation where this can be quite non-trivial, see [51]. Apart from the regulator function, the Wetterich equation only depends on the second functional field derivative of the effective action, which is denoted as Γ The trace in Eq. (1) represents a summation over internal indices as well as an integration over momentum space, which gives rise to a simple one-loop structure of the Wetterich equation since (Γ In general, the two-point function Γ (2) k [φ] has to be obtained from its respective flow equation, which is given by the second functional derivative of Eq. (1). This flow equation in turn depends on the three-and four-point functions, Γ , for which additional flow equations have to be derived and solved. This leads to an infinite tower of equations which has to be truncated at some point. To achieve this, one may chose a suitable ansatz for the functional form of the effective action. Standard expansion schemes for the effective action are for example given by vertex expansions or a derivative expansion, see e.g. [41]. While vertex expansions allow to take into account the full momentum dependence of n-point correlation functions, the derivative expansion can resolve the full field dependence of the corresponding terms in the effective action already at the lowest order, i.e. in the local potential approximation, see e.g. [52,53].
In this paper we use an ansatz for the effective average action which combines both aspects. It includes two-point functions containing arbitrarily many derivatives which are one-loop exact and self-consistent. All three and higher n-point functions are included as the scale-dependent but structure-less vertices obtained from the LPA. In a first exploratory study we here apply this scheme to the relativistic O(4) model in four space-time dimensions, where our ansatz for the effective average action is then of the form: Here, U k (φ 2 ) is the O(4) symmetric effective potential that includes the mass term and all higher-order local mesonic self-interactions, i.e. the scale dependent n-point vertices with n ≥ 3. The −cσ term is included to account for the explicit symmetry breaking. The model serves as a low-energy effective description of the light meson sector in QCD with N f = 2 flavors. The field φ then encodes the sigma and the pion fields, φ = (σ, π) T , and the term cσ resembles the current quark masses in QCD that explicitly break the corresponding chiral symmetry. The self-energies Σ k (x) for pions and the sigma meson are diagonal in field space, of one-loop structure, and renormalized by zero-momentum subtraction so that they vanish for constant fields in order not to interfere with the effective potential and to avoid double counting. To be precise, this implies that their two lowest moments vanish, To achieve this, the corresponding self-energies in momentum space, referred to as polarization functions Π k below, are renormalized by subtraction at zero momentum so that Π k (0) = 0 with Π k (p) ∼ p 2 for p 2 → 0.
To systematically extend the vertex-expansion scheme, at the next order one might include self-consistent self energies of two-loop structure in combination with a correspondingly soft-momentum subtracted three-point function of one-loop form, supplemented by local n-point vertices for n ≥ 4, and so forth. While straightforward in principle, these higher order expansions will eventually have to rely on numerical methods to evaluate the successively higher loop contributions to the lower n-point functions maintained in the expansion. At the present one-loop order we have the luxury to be able to use explicit and analytic expressions for the self-energies in the entire cut-complex plane of their analytically continued invariant momentum p 2 argument. This makes the analytic continuation of the flow equations for the two-point functions, to calculate the spectral functions in our aFRG framework, particularly transparent, and it therefore provides an illustrative and valuable first step. For constant fields, our ansatz reduces to the local potential approximation by design, and the Wetterich equation reduces to a flow equation for the effective potential U k (φ 2 ) as the only scale-dependent quantity on the left hand side, σ,k (q) + R k (q) where the integration is over four-dimensional momentum space. The two-point functions on the right hand side, on the other hand, will now include the sigma and pion polarization functions Π σ,k (q) and Π π,k (q) as detailed in Sec. II B below. For a better comparison, however, we will first briefly review the standard LPA procedure to calculate pion and sigma spectral functions in the next subsection.
A. Analytically continued LPA flows The LPA setup presented in this section was first proposed in [28] and will be briefly summarized here. The standard LPA form of the two-point functions can be obtained by taking two functional derivatives of the ansatz for the effective action, Eq. (2), here without the subtracted self-energies, which then simply yields where the squares of the (Euclidean or curvature) masses are given by, with ρ = φ 2 . These two-point functions are then inserted into the flow equation for the effective potential, Eq. (5), which leaves to specify the regulator function R k . In order to reproduce the previous LPA setup for the O(4) model spectral functions from the aFRG flows in [28] as our benchmark calculation we will here also use the threemomentum analogue of the LPA-optimized regulator [54] from that previous study, which is given by where Θ(x) is the Heaviside step function. It only regulates spatial momenta but not the energy components at the expense of some breaking of the Euclidean O(4) symmetry. This breaking was assessed and found to be negligible for external momenta well below the UV cutoff scale Λ in the Euclidean two-point functions, with still reasonably small and only quantitative effects in the time-like domain after analytic continuation [28]. While this can be avoided in principle [30,31], the three-dimensional regulators allow to perform the integration over the internal energy component or the corresponding Matsubara sum at finite temperature analytically which tremendously simplifies the analytic continuation procedure. The scale-dependent effective potential is then used as input in the flow equations for the two-point functions. These are obtained from the Wetterich equation by tak- ing two functional field derivatives, leading to see Fig. 1 for a diagrammatic representation. Explicit expressions for the involved loop functions I and J as well as for the four-point vertex functions Γ (4) are given in App. B. As we will use the same three-point vertex functions in the one-loop setup as in the LPA setup, we discuss them in more detail here. They are obtained from appropriate derivatives of the ansatz for the effective action in Eq. (2) which leads to We note that this setup is thermodynamically consistent, in that the flow equations for the three-point functions Γ (3) (p, q) can be expressed in terms of derivatives of the flow equation for the effective potential when evaluated at vanishing external momenta p = q = 0, The analogous consistency relation is known to hold for the two-point functions as well [28].
Since the vertex functions in Eqs. (11)-(12) only carry field and scale dependence, the momentum dependence is generated solely by the loop functions J. As shown in App. B, they carry the momentum dependence of the LPA propagators which in turn depend on the scaledependent curvature masses in Eqs. (8)- (9). The curvature masses therefore also determine the location of thresholds in the ensuing spectral functions which leads to a certain mismatch between the physical pole masses and decay thresholds. We also note that the non-trivial two-point functions obtained from these aFRG flows are not inserted back into the flow equations here. Although these two-point functions with their non-trivial momentum dependence in general provide results beyond any finite order in the derivative expansion, the self-consistency of these flows remains restricted to the lowest order LPA level, valid for the effective potential and the zero-momentum limit of the two-point functions. Fully self-consistent aFRG calculations are possible in principle, but technically and numerically very demanding as we will demonstrate in Sec. VI below. Moreover, the resulting IR two-point functions are then only known numerically which limits access to their structure in the complex plane and on different Riemann sheets. For all these reasons, we here propose a self-consistent but equally simple and illustrative one-loop setup as described in Sec. II B next.
We close this section by a brief discussion of the analytic continuation of the flow equations and the calculation of spectral functions. At zero temperature, the analytic continuation simply involves replacing the Euclidean energy p 0 by a real frequency ω in the following way, where the limit ε → 0 + is implicitly assumed. This prescription yields the flow equations for the retarded twopoint functions at real frequency ω which can then be solved numerically. Finally, the spectral functions are defined as the imaginary parts of the retarded propagators D R k = −D k (−iω + ), with ω + = ω + iε, and can be obtained from the corresponding real-time two-point functions (with retarded self-energies) as follows, In this section we present our new self-consistent oneloop FRG approach. It is based on parametrizing the two-point functions by analytic expressions for the scaledependent self-energies Π k (p) obtained from one-loop formulae given in Appendix A. By construction, these expressions are given by analytic functions in the cutcomplex p 2 -plane and allow to extract the spectral functions from the discontinuities of the corresponding propagators along this cut. With these analytic one-loop expressions, the type of processes that can occur, from Cutkosky's rule, are therefore the same as in the analytically continued LPA flows above. The parameters involved in these one-loop parametrizations are calculated self-consistently from the aFRG flows, however, as described below. In particular, the two-particle thresholds are then determined by the physical pole masses m p of the particles involved. This fixes a slight inconsistency of the aFRG LPA flows where these thresholds are determined by the curvature masses m c , cf. Eqs. (6)- (7), which represent the zero-momentum limit of the twopoint functions and the curvature of the effective potential in the corresponding field direction at the same time. However, these curvature masses will differ, in general, from the physical pole masses, especially for heavier particles where extrapolation from the pole position at −p 2 = (m p ) 2 to p 2 = 0 can induce sizable corrections due to the p 2 dependence of the self-energies.
With these self-energies included, the sigma and pion two-point functions are then given by Here, the self-energies can be split up into the contributions from the different loops, see also Fig. 2, where g σ,k and g σπ,k are the scale-dependent three-point σ and σ-π coupling constants, m p σ,k and m p π,k denote the scale-dependent pole masses. The renormalized selfenergy integrals Π R are defined in App. A and include the symmetry factor of 1/2 in the diagrams, while the factor of 3 in Eq. (22) reflects the three pion fields, and the factor of 2 in Eq. (23) comes from the degeneracy of the two lower loops in Fig. 2. The renormalization of the otherwise UV divergent self-energy integrals is performed by zero-momentum subtraction. Hence, these self-energies vanish at zero momentum, Π R (0) = 0, so that the relation between zero-momentum two-point function and effective potential or curvature mass, Γ (2),1L k (0) = (m c k ) 2 remains unchanged by construction. We also note that the self-energies are analytic with Π R (p 2 ) = O(p 2 ) at p 2 = 0 in presence of a mass gap.
For the scale-dependent coupling constants g σ,k and g σπ,k we will use the same expressions (13) and (14) as in the LPA setup, i.e.
but evaluated at the scale-dependent global minimum of the effective potential ρ 0,k . The effective potential is obtained from Eq. (5), as in the LPA setup, but now evaluated using the one-loop two-point functions, Eqs. (20)- (21). For the regulator function we will here use the four-dimensional Litim regulator, Given the coupling constants g σ,k and g σπ,k , one can immediately obtain the pole masses m p σ,k and m p π,k . For resonances, these are of course complex valued and off the physical Riemann sheet. When the widths are sufficiently small they are close to the physical sheet, and their real parts are then very well approximated by the zero crossing of the real part of the corresponding two-point functions along the timelike real axis. In the following, we therefore define what we mean by the pole masses from the zero crossings of the real parts of the corresponding two-point functions as our proxy for the real parts of the resonance positions, i.e. as the scale-dependent solutions with Eqs. (22) and (23) for the two self-energies. Note that their real parts are continuous along the cut and no ε-prescription is thus required here. Since both two-point functions depend on both pole masses, this coupled set of equations has to be solved simultaneously at each scale k.
Eqs. (26)- (27) are the self-consistency equations for the scale-dependence of the mass parameters which determine that of the resonance positions and decay thresholds at the same time in a consistent manner at the expense of an only moderate increase in numerical costs.
As the analytic form of the two-point functions is fixed by Eqs. (20)- (21), there is no need to solve flow equations for them. In fact, it can be shown that solving the flow equations for the two-point functions, Eqs. (11)-(12), for classical input, i.e. for scale-independent coupling constants and masses, exactly reproduces the one-loop result. The main effects of the LPA aFRG calculation of the two-point functions discussed in Sec. II A are therefore expected to be included in the one-loop setup proposed here, in particular when using scale-dependent parameters obtained from the FRG in the self-consistent one-loop two-point functions.
After solving the self-consistent system of equations for the effective potential and the two-point functions, the spectral functions can be obtained from the real-time two-point functions using Eq. (19). Given the analytic structure of the two-point functions, the spectral functions can be evaluated at arbitrarily high energies. This is not possible in the LPA setup since there the UV cutoff for the RG scale also limits the energy range where the spectral functions give meaningful results, see also the results shown in Sec. IV. Knowing the spectral functions at arbitrary energies also allows us to check their normalization which is given by the energy-weighted sum rule see for example [55] for a more detailed discussion of various sum rules for the spectral functions from aFRG flows.
We have verified numerically that this sum rule is indeed satisfied, independent of the particular values used for coupling constants and masses. Finally, we note that it is also possible to explore the analytic structure of the two-point functions in the complex-momentum plane by simply evaluating them at complex arguments. In particular, it is possible to investigate different Riemann sheets and thus to identify the complex pole of the sigma meson. For more details on how to access different Riemann sheets we refer to App. A while corresponding results are presented in Sec. V.

C. UV conditions and numerical implementation
For both the LPA aFRG and the self-consistent oneloop (SC1L) setup we choose the effective potential in the UV to be of the following form, while the explicit symmetry breaking term, −cσ, is added in the IR since it is a renormalization group invariant term. The parameters b 1 , b 2 and c are chosen such as to obtain phenomenological values for the pion decay constant as well as for the pole masses of the pion and the sigma meson, where, however, in FRG studies of the O(4) model the possible range for the sigma mass is typically below its experimental value, see also [28]. Explicit numerical values for the parameters chosen in the UV and the resulting IR values are given in Tab. I, where we use Λ = 500 MeV for the UV scale and k IR = 40 MeV as the IR scale. In order to obtain results for k < k IR we use an extrapolation. We note that the pion decay constant is identified with the global minimum of the effective potential in the IR, σ 0 , which in the SC1L setup needs to be The effective potentials U k (ρ) obtained from the LPA flow as well as the SC1L flow are shown in the UV and in the IR, both normalized to zero at σ = 0, and with the UV potentials divided by a factor of 2 for better comparison.
renormalized due to the non-trivial wave function renormalization factors, see Sec. III below for a more detailed discussion.
The flow equation for the effective potential, Eq. (5), is solved numerically using the so-called grid method which is based on a discretization of the field variable ρ into equidistant grid points, see for example [56]. The flow equation then reduces to a set of coupled ordinary differential equations which can be solved using standard techniques where we obtain the required first and second derivatives of the effective potential using a finitedifference discretization. In the present work we used different numbers of grid points in the range N ∈ [100, 300] distributed over a ρ-field range of ρ ∈ [0, 140 2 ]MeV 2 . We note that such a finite-difference implementation has certain limitations, see e.g. [57,58]. We therefore checked explicitly that another, more sophisticated numerical technique [59,60], i.e. the Kurganov-Tadmor finite volume technique [61], yields the same numerical results.
Within the LPA setup, the solution for the scaledependent effective potential is then used as input in the flow equations for the two-point functions, Eqs. (11)- (12). In the UV, the LPA real-time two-point functions have the following form, with the curvature masses defined in Eqs. analytically in the SC1L setup.

III. EFFECTIVE POTENTIAL, COUPLINGS, AND MASSES
In this section we present our results for the effective potential, the three-point coupling constants, as well as for the curvature masses and the pole masses. We begin with a discussion of the effective potential U k (ρ) for both the LPA setup as well as the SC1L setup, see Fig. 3. In the UV, the effective potential already shows spontaneous chiral symmetry breaking as its global minimum is at large values of the σ field. This is of course done by construction and necessary since the bosonic fluctuations tend to restore chiral symmetry. Each of the two UV potentials is chosen so as to obtain the same values for the pion decay constant and the pole masses in the IR. We also note that the UV potential in the SC1L setup starts out less strongly broken than in the LPA setup. This can be understood by considering the different propagators that enter the flow equations for the effective potential: In the SC1L setup, the propagators contain the one-loop self-energies which tend to suppress the bosonic fluctuations and hence the resulting symmetry restoration.
From the effective potential we obtain the pion decay constant f π which is identified with the value σ 0 of the σ field at its global minimum. In the case of the SC1L setup, this value has to be renormalized due to the appearance of non-trivial wave function renormalization factors Z π and Z σ . These factors arise from the nontrivial momentum dependence of the SC1L propagators The scale-dependent minimum of the effective potential for the LPA setup, σ 0,k , is shown together with the renormalized minimum obtained for the self-consistent oneloop setupσ 0,k .
and are defined as The renormalized minimum of the effective potential can then be obtained using the pion wave function renormalization factor,σ While in the LPA setup we have Z π,k = Z σ,k = 1 by definition, the wave function renormalization factors for the SC1L setup are shown in Fig. 4 as a function of the RG scale k. We note that Z π,k is almost scale-independent and close to 1 and therefore only has a small effect on the renormalization of the sigma field. The sigma wave function renormalization on the other hand shows a nontrivial scale dependence and is significantly larger, with an IR value of Z σ,k=0 = 1.58 as opposed to Z π,k=0 = 1.06 for the pion. Given the wave function renormalization factors, we can compute the renormalized minimumσ 0,k of the effective potential and compare it to the LPA result σ 0,k , see Fig. 5. As also evident from Fig. 3, the location of the minimum flows from larger values in the UV to smaller values in the IR, with the LPA minimum starting at larger values than the SC1L minimum. In the IR, the LPA minimum σ 0,k and the renormalized SC1L minimumσ 0,k arrive at the same value of f π = 93 MeV since the respective parameter sets have been chosen accordingly.
From the scale-dependent effective potential U k (ρ), we can also determine the coupling constants g σ,k and g σπ,k , cf. (24) with (13)- (14). Their scale dependence is shown in Fig. 6 for the LPA and the SC1L setup. We note that g σ,k is about three times larger than g σπ,k and that the SC1L couplings tend to be smaller than the LPA couplings, see Tab. II for explicit values.
We close this subsection with a discussion of the scale dependence of the various masses obtained within the LPA as well as within the SC1L setup, see Fig. 7.
For the LPA setup we show the Euclidean curvature masses as obtained from the effective potential and the real-time pole masses as obtained from the zero crossings of the real parts of the corresponding two-point functions.
In the UV, the LPA curvature and pole masses agree by construction. In the IR, the pole mass of the sigma meson is significantly smaller than its curvature mass while the pion pole mass stays close to its curvature mass, see also Tab. II. We note that the sigma meson pole mass suddenly jumps to a higher value at k ≈ 110 MeV which can be understood as follows. In the UV, the sigma meson starts out as a stable particle in the LPA setup with its spectral function given by a Dirac-delta function. This delta function then moves to smaller energies, i.e. pole masses, as the RG scale k is decreased. Eventually, the two-pion decay threshold, which also moves to smaller energies due to the decreasing pion curvature mass, overtakes the delta peak which leads to the formation of a broad resonance structure in the spectral function. This Re Im also leads to a sudden increase of the pole mass as visible in Fig. 7.
On the right-hand side of Fig. 7 we show the corresponding masses obtained within the SC1L setup. In this case we also show the renormalized curvature masses which are defined as an thus take the effect of the wave function renormalization factors into account. We note that, already in the UV, the masses are different from one another since the UV two-point functions, which determine the pole masses and the wave function renormalizations, already include the one-loop self energies and thus go beyond LPA. As in the LPA setup, however, the pion masses remain almost constant throughout the flow. The sigma meson masses show a stronger scale dependence and continuously decrease towards the IR, see Tab. II for explicit values.

IV. TWO-POINT FUNCTIONS AND SPECTRAL FUNCTIONS
We now turn to our results for the real-time two-point functions and spectral functions. In Fig. 8 we show the real and the imaginary part of the sigma and pion twopoint functions in the IR as obtained from the LPA aFRG flow as well as the SC1L setup. While the two frameworks give similar results at low external energies ω, one encounters undesired effects from the UV cutoff in We also note that the two-point functions show several kinks which correspond to the location of decay thresholds in the spectral functions. These decay thresholds are determined by the curvature masses in the LPA and by the pole masses in the SC1L setup, respectively. They are more clearly visible in the imaginary parts of the twopoint functions as shown on the right-hand side of Fig. 8. The first threshold in the sigma two-point function is determined by the process i.e. the decay of an off-shell sigma meson with energy ω into two on-shell pions which is kinematically possible when ω is larger than two times the pion mass. The second threshold in the sigma two-point function corresponds to the process i.e. the decay of an off-shell sigma meson into two onshell sigma mesons. The pion two-point function, on the other hand, is only affected by a single process, i.e. the decay of an off-shell pion into a sigma meson and a pion, We now turn to the resulting spectral functions which are shown in Fig. 9. On the left-hand side we show the spectral functions obtained from the LPA aFRG flow and SC1L flow for a small but finite value of ε = 0.1 MeV. Both approaches again give very similar results at lower energies while the LPA approach breaks down at energies close to the UV cutoff scale. The pion spectral function shows a sharp peak at the pion pole mass followed by the π * → σ + π decay threshold and the corresponding continuum part. Note that this decay threshold occurs at different energies in the two different calculations. With the analytically continued LPA flows it is always determined by the curvature-mass parameters, and it therefore here occurs at ω = m c σ + m c π ≈ 575 MeV. With the selfconsistent one-loop flows of the SC1L setup, on the other hand, it is determined by (the real parts of) their pole masses and hence occurs at ω = m p σ + m p π ≈ 460 MeV. This most significant difference is a direct consequence of the self-consistency of the SC1L flows for the two-point functions in the timelike domain.
In contrast, less noticeable differences are observed in the sigma spectral function. Because the difference between curvature mass defined at p 2 = 0 and pole mass at −p 2 = (m p π ) 2 is negligible for the pion in our calculations, as expected for the lightest meson, the σ * → π + π thresholds occur at practically the same energies in both calculations. The small differences in the sigma-meson spectral function near this two-pion threshold arise from the slightly different resonance pole positions of the sigma meson in both calculations, see Table II. Only visible in the self-consistent SC1L calculation, on the other hand, is the σ * → σ + σ threshold at higher energies where the analytically continued LPA flow starts to become cutoff sensitive and less reliable.
In order to exemplify the influence of the small but finite ε > 0 in these results, on the right-hand side of Fig. 9 we also show the corresponding results obtained in the limit ε → 0. In the case of the LPA aFRG flows, this result was obtained by performing the limit ε → 0 analytically in the flow equations for the imaginary parts of the two-point functions, see for example [37] for details. In the case of the SC1L setup this limit can readily be obtained numerically since the two-point functions are known analytically. In contrast to the ε > 0 results on the left, the pion peak is now given by a Dirac delta function and the spectral functions otherwise vanish below the lowest-lying decay thresholds. At energies larger than these lowest-lying thresholds, however, the spectral FIG. 10. The imaginary part of the pion (left) and the sigma (right) two-point function is shown in the complex momentum plane. The pion two-point function exhibits a branch cut corresponding to the π * → σ + π decay process while the sigma two-point function exhibits two branch cuts corresponding to the processes σ * → π + π and σ * → σ + σ. The complex pole of the pion is marked as a black dot.
functions are basically unaffected by this limit. We note that, at finite temperature, a finite value for ε may be used to qualitatively mimic the effects from multi-particle scattering processes in the thermal medium.

V. ANALYTIC STRUCTURE IN THE COMPLEX PLANE
In this section we discuss the analytic structure of the real-time sigma and pion two-point functions in the complex plane of the analytically continued energy variable z, with z = ω along the timelike real axis, as obtained from the self-consistent one-loop aFRG setup. Analyzing this structure is possible since these functions are given in analytic form, cf. Eqs. (20)-(21), while they are only known numerically in the LPA setup. We begin by studying the branch-cut structure of the imaginary part of the two-point functions which is shown in Fig. 10. The pion two-point function shows a branch cut starting at ω = Re z = m p π,k=0 + m p σ,k=0 ≈ 460 MeV corresponding to the process π * → σ + π while the sigma two-point functions shows two branch cuts, one starting at ω = Re z = 2m p π,k=0 ≈ 270 MeV corresponding to the process σ * → π + π and one at ω = Re z = 2m p σ,k=0 ≈ 650 MeV corresponding to the process σ * → σ + σ. The pion pole is represented by a black dot located on the real axis at ω = m p π,k=0 = 135 MeV. The complex pole of the sigma meson is not visible since it is not located on the so-called first Riemann sheet which is shown here.
In Fig. 11 we show contour plots of the real part, the imaginary part, and the absolute value of the pion twopoint function evaluated at complex z as analytic continuation of the energy ω along the real axis. The energy range is chosen such as to encompass the pole of the pion propagator, i.e. the zero-crossing of the two-point function. In this energy range no decay thresholds are present and hence only the first Riemann sheet is shown. We note that the real part is symmetric about the real axis and that the crossing of the zero-line with the real axis defines the pion pole. The imaginary part, on the other hand, is antisymmetric with respect to the real axis and the zero-line coincides with the real axis. Finally, the absolute value of the pion two-point function clearly shows the location of the pion pole as the point in the complex momentum plane where both the real part as well as the imaginary part vanishes.
We now turn to the sigma two-point function, which is shown in Fig. 12, again for complex energy z. This time, we show an energy range between the σ * → π + π threshold and the σ * → σ + σ where two Riemann sheets are of importance. Crossing continuously across the cut along the real axis from above one crosses from the first Riemann sheet in the upper half plane to the second Riemann sheet in the lower half plane, cf. App. A for a discussion of the different sheets. The complex pole of the sigma meson is then found on the second Riemann sheet at We note that the same value can also be obtained using numerical reconstruction techniques such as the Schlessinger point method (SPM), see [11,62,63]. This method can be directly applied to numerical data on the spectral function and thus can in principle also be used for numerical results from FRG flows such as our LPA calculation here. The branch cut structure of the sigma two-point function is more readily visible in a three-dimensional plot, see Fig. 13 where the first and the second Riemann sheet of the imaginary part of the sigma two-point function are shown together with the location of the complex pole of the sigma meson. The two Riemann sheets are continuously connected along the branch cut corresponding to the process σ * → π + π while they are discontinuous along the second branch cut corresponding to the process σ * → σ + σ. We note that there exist additional (infinitely many) Riemann sheets which can be accessed as described in App. A, but which are of no further relevance for our present study.
We close this section by investigating the movement of the poles of the pion and the sigma meson in the complex plane as a function of the RG scale k, as shown in Fig. 14. While the pion pole is always located on the real axis and moves from smaller values in the UV to slightly larger values in the IR, the complex pole of the sigma meson is located deep in the complex plane on the second Riemann sheet and moves from larger values in the UV to smaller values in the IR. We note that having access to the complete analytic structure of the two-point functions, including the location of the complex poles, represents a significant advantage over previous real-time FRG calculations and provides this additional insight. In particular, the techniques developed in this work can be extended to phenomenologically more realistic effective theories, e.g. involving additional (mesonic and baryonic) degrees of freedom at finite temperature and density [64].

VI. TOWARDS FULLY SELF-CONSISTENT SPECTRAL FUNCTIONS
The main virtue of the aFRG flows based on the SC1L ansatz for the effective average action as described in the previous section is its self-consistency beyond the softmomentum limit which is a substantial improvement over previous LPA based aFRG calculations. It is nevertheless straightforward to implement, and the explicit analytic expressions provide a very clear and intuitive understanding of the underlying analytic structures.
These merits of the explicit SC1L ansatz at the same time clearly illustrate its limitations. Based on the explicit analytic expressions for self-energies of one-loop form, the approach is not able to resolve vertex corrections and n-particle thresholds with n ≥ 3 at this lowest order for example. Including two-loop expressions for the self-energies together with one-loop structures in the 3point vertex functions at the next order, on the other hand, the approach will quickly lose its simplicity.
Leaving the systematics of this expansion aside for now, in particular by not attempting to resolve the momentum-dependent substructures of three-point and higher n-point vertex functions, we can test at least the effects of including the higher orders in the self-energies. In order to achieve that, in this section we demonstrate how the numerical self-energies can alternatively be obtained in a fully self-consistent calculation beyond loopexpansion. In a first exploratory study we will use this fully self-consistent framework here to calculate the same O(4) model spectral functions as above for comparison.
The main difference between the SC1L ansatz and the fully self-consistent calculation here will be the feedback of the full sigma two-point function reflecting its character as a broad two-pion resonance. As a result, the two-particle decay threshold π * → σ + π described above will then become a broader three-particle resonance decay π * → σ +π → 3π. Apart from this important further conceptual improvement, however, the overall differences will turn out to be otherwise surprisingly small. Our comparison here might therefore lend further support to continue using the considerably simpler SC1L flows in more realistic effective field theory calculations for technical reasons in the future.

A. Self-consistent Euclidean flows
The fully self-consistent aFRG setup is based on reinserting the full k-dependent two-point functions into the flow of the effective potential. For the Euclidean system of coupled flow equations for effective potential and twopoint functions this is rather straightforward. Based on the Euclidean results as input, the analytically continued flows will then be solved in a second step as described in the next subsection.
First we briefly describe how to solve the Euclidean system consisting of flow equations for Γ (2) π,k , Γ the effective potential U k ; the vertices Γ (3) k and Γ (4) k are extracted momentum independently from the effective potential as before. The resulting system of Euclidean flow equations then is of the following schematic form: As before, the flow of the effective potential is solved on a grid in field direction, whereas the flow equations of the two-point functions Γ (2),E π,k and Γ (2),E σ,k are evaluated at the physical minimum in the IR φ 0 = (σ 0 , π = 0). As a further simplification we assume the field dependence of the two-point functions to be of the following form, i.e. we split the two-point functions into field-independent but momentum-dependent parts plus a fully fielddependent but momentum-independent term which is identified with the Euclidean curvature mass extracted from the effective potential. This separation ansatz seems rather natural in view of the standard LPA or LPA truncation, cf. e.g. [29], but of course contains fully back-coupled objects calculated self-consistently here. The flow equations for the two-point functions contain contributions analogous to those in Eqs.  momentum dependent expressions on every grid-point at every scale k. In order to be ready for a straightforward extension to finite temperature in the future, in our numerical implementation we have introduced discrete explicit Matsubara sums for p 0 n = 2πnT at a small but finite Temperature of T = 10 MeV, and a separate rotationally invariant grid for the magnitude of the spatial momentum variable | p| using an extrapolation for large momenta outside the grid, if necessary. As regulator function we use the standard three-dimensional Litim regulator, cf. Eq. (10). We have checked explicitly that the number of Matsubara modes we included (n max ≈ ±100) are sufficient, and that the specific extrapolation procedure does not affect the results.
For the initial conditions at the UV scale Λ = 500 MeV we use Eq. . And as before, we add the explicit symmetry breaking term cσ with c = 0.014 Λ 3 to the resulting effective potential at the IR cutoff scale. The parameters b 1 , b 2 and c are again chosen to more or less reproduce the same phenomenological values for pion pole mass, with m p π = 135.5 MeV, and decay constant, with f π ≡ σ 0 = 93 MeV. 1 Also the pole mass of the sigma meson then results in the same range as from the previous calculations. Since the numerical effort of further fine-tuning of the UV initial parameters gets unreasonably expensive in the fully self-consistent calculation, we are content here with the infrared values summarized in Tab. III. For the purpose of comparison of the resulting self-consistent spectral functions with the considerably cheaper LPA and SC1L results these IR values are sufficiently close to those in Tab. I.

B. Self-consistent aFRG flows
Using the Euclidean results as input, we solve the analytically continued flow equations for real and imaginary parts of the (retarded) two-point functions of pion and sigma meson self-consistently in the next step. Because we have already calculated the field-dependent but momentum independent constants m 2 E [φ], the remaining system consists only of the momentum-dependent contributions to the flow. The corresponding system of flow equations in the Minkowski spacetime is obtained from its Euclidean counterpart by applying the same analytic continuation procedure as above, cf. Eq. (17), in combination with Källén-Lehmann spectral representations for regulated propagators. The finite-temperature extension of the analytic continuation procedure is then also straightforward and parallels that used for the LPA based aFRG flows [32].
In order to analytically continue the system of FRG flow equations for the retarded two-point functions, we first consider Källén-Lehmann spectral representations for the scale-dependent regulated Euclidean propagators with scale dependent spectral functions, which read with α ∈ {σ, π} for the different fields. The existence of such scale-dependent spectral representations follows from Cauchy's theorem and the use of three-dimensional regulator functions. Inserting these spectral representations of the propagators in the momentum-dependent loop functions J k (p) on the right-hand side of the flow equations, see Eq. (B2), these are of the following form, × ω 1 ω 2 ω 3 ρ α,k (ω 1 , q)ρ α,k (ω 2 , q)ρ β,k (ω 3 , p − q) ω 2 1 + q 2 0,n ω 2 2 + q 2 0,n (ω 2 3 + (p 0 − q 0,n ) 2 ) .
In these expressions we can then perform the Matsubara sum over q 0,n = 2πT n analytically and simply apply the standard analytic continuation procedure for the external p 0 , cf. Eq. (17), after employing the 2πT periodicity of the resulting statistical distribution functions [32]. The loop functions in (44) then obviously still involve the spatial momentum integration. A major simplification thereby occurs, if we assume the scale-dependent spectral functions to be Lorentz invariant as they should be at zero temperature, i.e.
We can then substitute s = p 2 = ω 2 − p 2 which amounts to using the zero-temperature spectral representation in (44). Note that this is strictly speaking only valid at vanishing temperature, and for regulator functions that do not break the corresponding Euclidean O(4) invariance. The latter is of course not the case for the threedimensional regulators considered here, but the effect of the breaking of boost invariance induced by this was assessed already in [28] and found to be negligible as long as the magnitude of the spatial momentum stays well below the UV cutoff scale Λ. Also at finite temperature this will hold at best approximately, if the temperature is sufficiently low. Collisional broadening is not included in the spectral functions in this way [65], for example. We therefore assume Lorentz invariance to hold at least approximately here, in order to reduce the number of momentum arguments that have to be taken into account to resolve the momentum dependence of the two-point function. Eq. (44) then simplifies to J αβ,k (p 0 , p) = q s1,s2,s3 .
Here, we can in principle now also carry out the remaining spatial momentum integration and are hence left with the three spectral integrals over the invariant s i ∈ [0, ∞) variables for every propagator involved in the loops on the right-hand side of the flow equations. These have to be performed numerically in every integration step. The analytic continuation can then proceed in the standard way as described, and the originally Euclidean system of flow equations then turns into aFRG flows for the real and imaginary parts of the retarded two-point functions Γ (2),R α,k (ω) = −Γ (2) α,k (−iω + ) in Minkowski spacetime. We use p = 0 here and hence obtain complex aFRG flow equations of the following form k [φ 0 ]; ρ π,k , ρ σ,k .
The k-dependent spectral function, here with notation ρ k (ω) ≡ ρ k (ω, 0), are thereby recomputed at every integration step of the scale k from for the pion and sigma meson respectively, where the form of the shift in the real part is due to the inclusion of the three-dimensional Litim regulator (10) and the continuation procedure, which throws in a global minus sign from the definition of the two-point functions.
To further reduce the numerical effort, we further consider the product (in momentum space) of the identical propagators in the loop functions, and assume the same form of spectral representation to exist for this product also. Since both propagators are analytic functions in the cut-complex p 2 plane, the product is as well. Positivity of the discontinuity along the cut is not required here. The effective spectral function ρ eff k of the product then reads (again for the three-dimensional sharp regulator and vanishing spatial momentum) In this way we can further reduce the number of spectral integrals in Eq. (44) from three down to two. The effects of using such a spectral representation for products of propagators should certainly be tested more carefully in the future. Our first results here are very promising, however, as we demonstrate below. The form of the loop functions J k (p 0 ) in (48), cf. Eq. (B2), at vanishing spatial momentum and T = 0 for example then simplify as follows, The remaining q 0 integral can in principle also be calculated analytically, but the resulting expression is lengthy and not particularly illuminating. Most importantly, it still is an analytic function in the cut-complex p 2 0 plane and can thus be analytically continued in the external frequency via Eq. (17) as needed for our aFRG approach. The explicit expression of the resulting J R k (ω) ≡ J k (−iω + ) is given in App. C.
In the numerical calculation, the system of equations given by Eqs. (48) is solved on a grid in frequency ω. In practice, it is therefore best to precompute the kernel of the spectral integrals in (52) and store it as a fixed matrix on this frequency grid, from which the two-dimensional spectral integrals of the form (52) are then evaluated numerically in every k step. The necessary cutoff Λ ω in these spectral integrals determines the size of the ω grid and is chosen much larger than the RG cutoff scale Λ. We moreover extrapolate the spectral functions for ω > Λ ω with ρ(ω) ∝ 1/ω 2 and verified explicitly that the results do not depend on the specific value used for this frequency cutoff Λ ω . The initial conditions for the real and imaginary parts of the two-point functions at the UV scale Λ are chosen as in (30) from the effective potential, C. Self-consistent aFRG results and comparison Our first results for the O(4) model spectral functions from the fully self-consistent aFRG flows, taken from [66], are shown in Fig. 15. These are based on the spectral representations as explained in the previous subsection. The discretization of the frequency for the Euclidean input here corresponds to a small residual temperature of T = 10 MeV as mentioned above.
Because distributions are generally difficult to integrate numerically, we have to keep the imaginary parts of the frequencies in all retarded correlations positive with a small but finite ε here. More sophisticated principal value prescriptions are needed to be able to take the limit ε → 0 analytically in spectral-representation based aFRG flows before they are solved numerically [67].
In this first exploratory study we simply stay off the real frequency axis a little bit. Smaller values of ε require finer frequency grids for the analytically continued aFRG flows, and hence increased numerical costs. In order to exemplify the influence of these residual imaginary parts, the results of two such calculations with different values of ε = 1 MeV and ε = 0.5 MeV are shown in Fig. 15.
The vertical red lines in Fig. 15 thereby indicate the location of the pion pole mass (at ω = m p π ≈ 135.5 MeV) in the pion spectral function and multiples thereof: (i) twice that for the two-pion threshold in the sigma spectral function (i.e. at ω = 2m p π ≈ 271 MeV); and, most significantly for the first time here, (ii) also three times that for the three-pion threshold (i.e. at ω = 3m p π ≈ 406.5 MeV) in the pion spectral function again.
While the sigma meson spectral function shows the usual pronounced two-pion threshold at the right location of twice the pion pole mass, as it did already in our SC1L flows above, the figure also nicely illustrates the main effect of the full self-consistency here: towards smaller values of ε the threshold in the pion spectral function nicely settles at ω = 3m p π as it should, because the self-consistent sigma correlations in the flow now for the first time represent the broad two-pion resonance that the sigma should be. As a result, they contribute to the three-particle resonance decay π * → σ + π → 3π which is impossible in any calculation based on one-loop structures as in our SC1L flows above.
This important qualitative improvement comes at the price of an enormous increase in numerical complexity and costs, however. The SC1L aFRG flows of the previous sections are straightforward to implement and numerically very inexpensive compared to the fully selfconsistent flows for the spectral functions, despite the simplifications we have used in this first study for the latter. 2 We therefore close this section with a direct comparison of the fully self-consistent spectral functions with the corresponding ones from the self-consistent SC1L flows in Fig. 16. In order to compare apples with apples instead of oranges we have recomputed the more economic SC1L results of Sec. IV for this comparison once more, with the same value of ε = 0.5 MeV as used also in the fully self-consistent aFRG flows.
Although other parameters such as temperature and sigma pole mass for example do not match perfectly, compare Tabs. II and III, the overall results of both schemes, with vastly different computational complexity, then look surprisingly similar. The particle peaks in the pion spectral functions are almost the same, with an only marginal relative shift due to the choice of parameters which lead to pion pole masses that differ by irrelevant 0.5 MeV.
But even the decay threshold in the pion spectral function is not so dramatically different. While it starts at different values, i.e. at ω = m p π + m p σ ≈ 460 MeV in the SC1L result as compared to ω = 3m p π ≈ 406.5 MeV in the fully self-consistent calculation, this qualitatively important shift is mainly due to a smearing of the unphysically sharp threshold observed in the self-consistent one-loop ansatz for the effective average action. The fully self-consistent aFRG calculation captures the broad nature of the 3-particle resonance decay π → π + σ → 3π certainly in a physically more realistic manner than the SC1L calculation, where the resonance decay in the sec- π -fully self-consistent π -self-consistent 1-loop σ -fully self-consistent σ -self-consistent 1- loop   FIG. 16. The pion and sigma spectral functions are shown vs. external energy ω in the IR as obtained from the selfconsistent one-loop setup and the fully self-consistent framework for ε = 0.5 MeV.
ond step is missing, and which therefore produces an unrealistically sharp threshold for π → π + σ. While it is reassuring to know that it is possible in principle to describe the three-particle resonance decay correctly, the question might arise whether the gain is worth the effort. Especially when one is interested in heavy-ion collision observables such photons and dileptons for the electromagnetic spectral function, for example, which are measured from all stages of a collision and hence convoluted over the various temperatures and densities that occur during the evolution of the fireball, small differences such as those around threshold in the pion spectral function observed here will hardly be observable and will therefore not be relevant from a purely phenomenological perspective.
This kind of smearing of sharp structures is a general trend in the fully self-consistent and hence more realistic calculation. As such, it is also observed in the sigma spectral function, where one can also clearly see that the threshold of the σ * → σ+σ process (which corresponds to four pions in the final state now) gets smeared out completely and is no-longer visible in the fully self-consistent result. The sharp two-pion decay threshold on the other hand is located at the same energies again, owing to the self-consistency of both approaches and the nearly identical values for the pion pole masses in the two calculations. In conclusion, the comparison overall demonstrates how robust our aFRG results for spectral functions are at the present level.

VII. SUMMARY AND OUTLOOK
The main message of this work is the efficiency of the semi-analytic framework proposed here, which is referred to as the SC1L FRG flows above, to calculate real-time correlation functions and spectral functions within the FRG from a self-consistent albeit still economic ansatz that includes parametrized self-energies directly in the effective average action. The explicit expressions for the self-energies of oneloop form in this approach have the correct domain of holomorphy of the two-point functions implemented by construction. The analytic continuation from Euclidean to Minkowski spacetime is therefore straightforward. The scale-dependent parameters in these self-energies are couplings and (real parts of resonance) pole masses. It is only these few parameters that have to be calculated self-consistently during the FRG flow which makes this scheme comparatively economic.
We have compared our results for the spectral functions of pion and sigma meson in the O(4) model obtained from these self-consistent SC1L flows first with our previous analytically continued aFRG flows based on the same local potential approximation (LPA) for the flows of effective potential and two-point functions. Because of this restriction, the self-consistency is then also restricted to the zero-momentum limit of the two-point functions obtained from the corresponding analytically continued aFRG flows. In particular, the non-trivial momentum dependence of these two-point functions is not fed back into the flows for effective potential and two-point functions, which results in a mismatch between the physical pole masses and the location of decay thresholds. This restriction is avoided in our new SC1L flows which are fully self-consistent. As an added bonus, the analytic structure of the two-point functions can be studied not only in the entire complex frequency plane but on all Riemann sheets. This is very illustrative and, in particular, it allows to identify the complex poles on unphysical Riemann sheets associated with resonances which are generally difficult to access numerically otherwise.
In order to overcome the structural limitations due to the underlying one-loop form in the parametrizations of the SC1L self-energies, however, fully self-consistent aFRG calculations are necessary where the complete frequency and momentum dependence of the two-point functions is inserted back into the complete set of Euclidean and analytically continued flows. This is possible in principle by using spectral representations with scaledependent spectral functions that are recomputed selfconsistently in every numerical integration step of the FRG scale k. We demonstrate the feasibility of this procedure here in an exploratory study of fully self-consistent O(4) model spectral functions with several further simplifications for comparison. Our fully self-consistent results demonstrate that it is possible in principle to describe processes, from Cutkosky's cutting rule, that are beyond the one-loop structure of the previous approaches. Most significantly here, this is the three-particle resonance decay π * → σ + π → 3π which we observe at the correct three-pion threshold for the first time in this approach. The fully self-consistent aFRG scheme at this level is basically as good as it gets without including momentum dependent three and higher n-point vertex corrections.
Having said that, our comparison of the fully selfconsistent aFRG with the SC1L results mentioned above also shows that this conceptually important improvement essentially arises by a smearing around the unphysically sharp threshold for π * → σ + π of the latter. It might therefore quantitatively not be so significant in phenomenological applications, when there are other sources of smearing as well. Overall, our comparison can then serve to justify using the numerically considerably more economic and straightforward to implement SC1L framework which turned out to provide otherwise very robust results.
The self-consistent SC1L approach can thus serve as a starting point for phenomenologically more realistic calculations of real-time quantities such as spectral functions and transport coefficients in a strongly interacting warm and dense medium. We therefore plan to extend this SC1L framework to finite temperature and density next, and to then apply it to effective hadronic theories for nuclear and neutron matter. This will eventually include also the electroweak response from calculating (axial-)vector and electroweak spectral functions of strong-interaction matter in the future.
These renormalized expressions are used in Eqs. (22)- (23) for the calculation of sigma and pion self energies.
The analytic continuation is performed by replacing the Euclidean external momentum by a real-time energy ω as discussed in Sec. II A.
Note that Π R (p 2 ) in Eq. (A5) is an analytic function in the cut-complex p 2 -plane as it must, from our zerotemperature one-loop calculation, here. It is therefore very well suited for parametrizing self-energies with the correct analyticity properties. The branch cut for timelike −p 2 = s starts at s ≥ (M + m) 2 where the imaginary parts of both artanh's have the same sign and add up.
In contrast, for 0 ≤ s ≤ (M − m) 2 the two imaginary parts have opposite sign and cancel one another, while the apparent square-root cut for (M − m) 2 ≤ s ≤ (M + m) 2 is removed as long as the same Riemann sheet is used for the square-root in the prefactor as in the argument of artanh (with the identity −i artanh(iz) = arctan(z)).
Finally, the artanh(z) has branch cuts along the real axis for z ∈ (−∞, −1] and [1, ∞) with infinitely many Riemann sheets which are separated by multiples of iπ. We can therefore add or subtract multiples of iπ to the artanh functions in Eqs. (A5)-(A6) in order to access different Riemann sheets. This allows to explore the self energies and hence the two-point functions on different Riemann sheets as well. The second Riemann sheet for example is obtained by adding +iπ to all artanh functions in Eqs. (A5)-(A6), see also [63] for a similar discussion.

Appendix B: LPA loop and vertex functions
In this appendix we provide explicit expressions for the loop and vertex functions used in the flow equations for the two-point functions in the LPA setup, see also [28]. The loop functions are defined as J αβ,k (p) = q ∂ k R k (q)(D α,k (q)) 2 D β,k (q + p) , (B2) with the scale-dependent Euclidean propagator D α,k (q) = Γ (2) αα,k (q) + R k (q) and α ∈ {σ, π}. For the three-dimension Litim regulator, Eq. (10), we then get the following explicit expressions for p = 0, J αβ,k (p) = k 4 12π 2 (E α,k + E β,k ) 2 (2E α,k + E β,k ) + E β,k p 2 0 E 3 α,k E β,k (p 2 0 + (E α,k + E β,k ) 2 ) 2 , The four-point vertex functions can be obtained from the ansatz for the effective potential, Eq. (2), and are given explicitly by Γ (4) σσσσ,k = 12U Γ (4) σσππ,k = 4U Γ (4) ππππ,k = 176/3 U where the last vertex function contains contributions from the case where the internal pion is the same as the external one as well as from the other case where they are different.
The structural similarity with the analytic expression of the self-energy in Eq. (A5) is of course no coincidence. The discussion of the analytic structure follows the same line of arguments. In particular, it also has branch cuts along the real axis for ω 2 ≥ ( √ s 1 + √ s 2 ) 2 .

(C4)
Finally, for k 2 < s 1 , s 2 , a very good further approximation of K J k (s 1 , s 2 , ω 2 ) is obtained from expanding The leading term of this expansion was actually used for simplicity to produce the results presented in Sec. VI C above. The higher-order corrections affect the flow only close to the UV cutoff. These might in turn induce some corrections to the spectral functions of pion and sigma meson for high frequencies ω, of the order of the UV cutoff Λ as well, which are however not the focus of our comparison with the considerably more economic SC1L calculations here.