{'fraction_non_alphanumeric': 0.03816079037317975, 'fraction_numerical': 0.011362394843220802, 'mean_word_length': 4.385291027526791, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-equilibrium Functional Renormalization for Driven-Dissipative Bose-Einstein Condensation

We present a comprehensive analysis of critical behavior in the driven-dissipative Bose condensation transition in three spatial dimensions. Starting point is a microscopic description of the system in terms of a many-body quantum master equation, where coherent and driven-dissipative dynamics occur on an equal footing. An equivalent Keldysh real time functional integral reformulation opens up the problem to a practical evaluation using the tools of quantum field theory. In particular, we develop a functional renormalization group approach to quantitatively explore the universality class of this stationary non-equilibrium system. Key results comprise the emergence of an asymptotic thermalization of the distribution function, while manifest non-equilibrium properties are witnessed in the response properties in terms of a new, independent critical exponent. Thus the driven-dissipative microscopic nature is seen to bear observable consequences on the largest length scales. The absence of two symmetries present in closed equilibrium systems - underlying particle number conservation and detailed balance, respectively - is identified as the root of this new non-equilibrium critical behavior. Our results are relevant for broad ranges of open quantum systems on the interface of quantum optics and many-body physics, from exciton-polariton condensates to cold atomic gases.


I. INTRODUCTION
In recent years, there has been tremendous progress in realizing systems with many degrees of freedom, in which matter is strongly coupled to light. 1 This concerns vastly different experimental platforms: In ensembles of ultracold atoms, the immersion of a Bose-Einstein condensate (BEC) into an optical cavity has allowed to achieve strong matter-light coupling, and lead to the realization of open Dicke models; 2-4 In the context of semiconductor quantum wells in optical cavities, non-equilibrium Bose condensation has been achieved 5-7 -here the effective degrees of freedom, the exciton-polaritons, result from a strong hybridization of cavity light and excitonic matter degrees of freedom. 1, 8,9 Further promising platforms, which are at the verge of the transition to true many-body systems, are arrays of microcavities [10][11][12][13] or trapped ions, 14,15 as well as optomechanical setups. [16][17][18] Those systems have three key properties in common. First, they are strongly driven by external fields, such as lasers, placing them far away from thermodynamic equilibrium even under stationary conditions. Equilibrium detailed balance relations therefore are not generically present. Second, they exhibit the characteristics of quantum optical setups, in that coherent and dissipative dynamics occur on an equal footing, but at the same time are also genuine many-body systems. Finally, a third characteristic is the absence of the conservation of particle number. In particular, the admixture of light opens up strong loss channels for the effective hybrid light-matter degrees of freedom, and it becomes necessary to counterpoise these losses by continuous pumping mechanisms in order to achieve stable stationary flux equilibrium states. The pumping mechanisms can be either coherent or incoherent. In the latter case, e.g., single particle pumping directly counteracts the incoherent single particle loss; once it starts to dominate over the losses, a second order phase transition results on the mean-field level, in close analogy to a laser threshold.
At this point a clear difference between the quantum optical single mode problem of a laser and a driven-dissipative many-body problem becomes apparent: While the inclusion of fluctuations in the treatment of a laser smears out the meanfield transition, in a system with a continuum of spatial degrees of freedom a genuine out-of-equilibrium second order phase transition with true universal critical behavior can be expected. The theoretical challenge is then to understand the universal phenomena that can emerge due to the many-body complexity in a driven non-equilibrium setting.
In this work we address this challenge, focusing on a key representative that shows all the above characteristics: The driven-dissipative Bose condensation transition, relevant to experiments with exciton-polariton condensates, or more generally to any driven-dissipative system equipped with a U(1) symmetry of global phase rotations tuned to its critical point. We provide a comprehensive characterization of the resulting non-equilibrium critical behavior in three dimensions, extending and corroborating results presented recently. 19 A key finding concerns the existence of an additional, independent critical exponent associated with the non-equilibrium drive. It describes universal decoherence at long distances, and is observable, e.g., in the single particle response, as probed in homodyne detection of exciton-polariton systems. 20 This entails evidence that the microscopic non-equilibrium character bears observable consequences up to the largest distances in driven Bose condensation. Furthermore an asymptotic thermalization mechanism for the low frequency distribution function is found. Such a phenomenon has been observed previously in other contexts. [21][22][23][24][25][26][27][28][29] Here it is reflected in a symmetry that is emergent in the critical system on the longest scales.
By contrast, in systems at true thermal equilibrium, this symmetry is present at all scales as a microscopic symmetry. It then places severe restrictions on the relations between the noise and the coherent and dissipative dynamics in the system, 30,31 leading to fluctuation dissipation relations valid at all frequencies and wavelengths. This is the case, e.g., in the models for dynamical universality classes established by Hohenberg and Halperin (HH). 32 Non-equilibrium perturbations to these models that have been discussed in the literature concern, e.g., modifications of the noise term by spatial anisotropies, violating fluctuation-dissipation relations on a microscopic scale. For models without conserved order parameter, such as model A (MA) of HH, it has been shown that this does not lead to the existence of new universal critical behavior, but rather to a modification of nonuniversal amplitude ratios. 33,34 Genuinely non-equilibrium universal critical behavior has been found in several classical, driven systems with different microscopic origins. Examples include models with conserved order parameter with spatially anisotropic temperature, 33 the driven-diffusive lattice gas, 35 reaction-diffusion systems, [36][37][38][39] the problems of directed percolation 40,41 and self-organized criticality, 42 or kinetic roughening phenomena such as described by the Kardar-Parisi-Zhang equation. [43][44][45][46] At the technical level, the purpose of this paper is to lay out a general framework for addressing universal critical phenomena in open markovian many-body quantum systems. This framework may be further generalized and applied to a large variety of non-equilibrium situations, such as driven or driven-dissipative systems with different symmetries, [47][48][49][50] driven-dissipative systems with disorder, 51 and even superfluid turbulence. [52][53][54][55] We start from a microscopic, second quantized description of the system in terms of quantum master equations, and show how to translate the master equation into a Keldysh real-time functional integral, which opens up the toolbox of well-established techniques of quantum field theory. Next, we develop a functional renormalization group (FRG) approach based on the Wetterich equation, 56 which allows us to compute both the dynamical critical behavior as well as certain non-universal aspects of the problem. For example, in addition to determining critical exponents we can also extract a Ginzburg scale which marks the extent of the critical fluctuation regime.
The paper is organized as follows. In the next section we present our key results and sketch the resulting physical picture. Section III introduces to our model and provides the mapping of the master equation to an equivalent Keldysh functional integral. Using this framework, in Sec. IV we reproduce the results from mean-field and Bogoliubov theory, and show how the physics of a semi-classical drivendissipative Gross-Pitaevskii equation emerges naturally as a low frequency limit of the full quantum master equation. We highlight the additional challenges which arise from the need to treat a continuum of spatial degrees of freedom in order to capture critical behavior, and show in Sec. V how they are properly addressed by means of the FRG approach. The precise manifestation of the non-equilibrium character of the problem is worked out in Sec. VI. A detailed comparison of our non-equilibrium versus more conventional equilibrium models highlights a symmetry which is only present in thermal equilibrium and expresses detailed balance. We summarize the computation of the flow equations in Sec. VII and explain the hierarchical structure of the universal critical behavior implied by the flow. In Sec. IX we discuss the numerical analysis of the flow equations. We conclude in Sec. X.
At this point, we remark that the physical picture described in this work, and summarized in the following section, has been fully confirmed and further developed in a recent complementary perturbative field theoretical study presented in Ref. 57. There, in particular, analytical estimates for the critical exponents are provided.

II. KEY RESULTS AND PHYSICAL PICTURE
Driven-dissipative Bose condensation transition -Driven open quantum systems are commonly modeled microscopically by means of quantum master equations or in terms of Keldysh functional integrals as shown below. Starting from such a microscopic model of a driven Bose condensate we derive in Sec. IV C an effective long-wavelength description of the critical dynamics. The result, after dropping all irrelevant terms in the sense of renormalization group (RG), is a stochastic equation of motion for the order parameter, which may be cast in Langevin form, Such a dissipative stochastic Gross-Pitaevskii equation has been used as a model for exciton-polariton condensates. [58][59][60][61][62][63] This equation includes terms describing coherent dynamics, as well as ones capturing the dissipative processes and the drive. The coherent terms are the inverse mass A = 1/ (2m), the chemical potential µ and the elastic two-body interaction λ, whereas dissipative contributions include a kinetic coefficient D, the effective single-particle loss rate κ 1 as the difference between single-particle loss and pump rates, as well as two-body loss κ. The loss and gain processes induce noise, which is taken into account by the Gaussian white noise source ξ of strength γ with zero mean, ξ(t, x) = 0, and correlations Unlike the models of critical dynamics classified in Ref. 32, the coherent and dissipative terms in a driven condensate stem from completely independent physical processes. In particular, this implies that the steady state of the Langevin equation (1) is not characterized by a thermal (Gibbs) distribution of the fields, and this leads to the distinct critical behavior analyzed in this paper. Equation (1) admits a time-independent homogeneous mean-field solution |ψ 0 | 2 = −κ 1 / (2κ) if the single-particle pump rate exceeds the corresponding loss rate, i.e., the effective single-particle loss rate κ 1 becomes negative, and the chemical potential is set to be µ = 2λ |ψ 0 | 2 . Thus at the meanfield level a continuous transition is tuned by varying the single particle pump rate: ψ 0 vanishes as κ 1 goes from negative values to zero. Mean-field theory, however, breaks down in the vicinity of the phase transition as the inclusion of fluctuations may induce non-trivial scaling behavior or even render the transition first order. 64,65 In this paper we verify that the system described by (1) in three spatial dimensions indeed has a critical point characterized by universal dynamics. We argue that this dynamics is governed by a "Wilson-Fisher" like fixed point, but with another layer of dynamical critical behavior that is not found in non-driven systems.
Universality and extent of the critical domain -Our main technical tool for the analysis is a functional RG carried out for the dynamical problem. Emergence of a universal critical point is evident from the flow of the coupling constants to a fixed point independent of the initial conditions, as long as the system is tuned to the phase transition (cf. also Sec. IX). This is demonstrated in Fig. 1 (a), showing the flow of the real and imaginary parts of the complex interaction parameterũ 2 =λ + iκ (see Sec. VIII A). At the fixed point the real parts of all couplings vanish, which implies that the effective long-wavelength dynamics is purely dissipative. Integrating out fast fluctuations in the course of the RG flow, therefore, leads to a loss of coherence.
An important aspect of the phase transition which we analyze in detail here concerns the extent of the critical domain, which is delimited by the Ginzburg momentum scale k G . Knowledge of this non-universal scale is important for assessing the requirements from experiments aimed at measuring the critical phenomena. We find it to be given by (cf. Sec. IX) where C ≈ 14.8 and γ Λ , κ Λ , and D Λ are, respectively, the noise-strength, two-body loss rate and dissipative kinetic coefficient appearing in the description of the system at a mesoscopic scale Λ (see Sec. V B). Here we confirm this behavior quantitatively by a full numerical solution of the flow equations outside the critical domain, highlighting the capability of the FRG approach to compute universal and non-universal physics in a single framework. Asymptotic thermalization of the distribution function -An interesting result of the RG analysis is that the distribution function of the order parameter field at the critical point effectively thermalizes at long wavelengths and low frequencies. The effective thermalization is manifest as an emergent symmetry of the equations of motion at the fixed point that is not present at the mesoscopic level, cf. Sec. VI. For this reason the dynamical critical exponent z is the same as that of MA of the equilibrium classification. The presence of this symmetry implies a fluctuation-dissipation theorem (FDT), or, more physically speaking, a detailed balance condition valid at asymptotically long wavelengths.
In order to better understand this aspect, consider an equilibrium problem with detailed balance. This means that all subparts of the system are in equilibrium with each other. In other words, temperature is invariant under the system's partition in such a state. This statement is easily translated into a RG language: Natural system partitions are the momentum shells. Partition invariance of the temperature thus becomes a renormalized two-body couplingũ 2 =λ + iκ (see Sec. VIII A) for various initial valuesũ 2Λ . As a result of finetuning the initial values w Λ of the dimensionless mass parameter close to criticality, all flow trajectories approach the Wilson-Fisher fixed pointũ 2 * = i5.308 (indicated by the black dot) before eventually bending away. Shown are numerical solutions to the flow equations for r KΛ = 10, r u 3 Λ = 1,κ 3Λ = 0.01, and values ofũ 2Λ lying on a rectangle with sidesλ ∈ [0, 10],κ = 2, 10 andλ = 10,κ ∈ [2, 10]. (See Sec. VIII A for definitions of the parameters.) (b) Flow ofκ as a function of the dimensionless infrared cutoff t = ln (k/Λ) for various starting valuesκ Λ . Dots on the lines indicate the extent of the critical domain, which is set by the Ginzburg scale Eq. (3). Initial values are the same as in (a), apart fromκ Λ = 0.1, 1, 2, . . . , 10 and r u 2 Λ = 10. scale invariance of temperature under renormalization, which successively integrates out high momentum shells. The "equilibrium symmetry" expresses precisely this physical intuition.
In a non-equilibrium problem such as the driven condensate we discuss, this symmetry is in general absent at arbitrary momentum scales. In order to demonstrate how it emerges at long scales, we compute the scale dependence of an effective temperature, entering the (non-equilibrium) FDT, cf. Sec. VI. Indeed, we find scale dependent behavior at high momenta, which becomes universal and scale independent within the critical region delimited by the Ginzburg scale, cf. Fig. 2.
We note that, in principle, it is conceivable that the system might allow for different stationary scaling solutions far from equilibrium with different universal scaling behavior, not captured in the present approach. Indeed, in two dimensions, such a scenario is realized. 66 In three dimensions, however, such a behavior could be present only beyond a threshold value for the microscopic strength of violation of detailed balance.
Hierarchical shell structure of non-equilibrium criticality -A key result of the RG analysis is the hierarchical organization of the non-equilibrium criticality. This structure consists of three shells of critical exponents. The innermost shell in this hierarchy contains the two independent exponents ν, η describing the static (spatial) critical behavior of the classical O(2) model. 67 We find that the static exponents coincide with those of an ab initio computation of the classical O(2) exponents at the same level of approximation. Thus the nonequilibrium conditions do not modify the static critical behavior.
The intermediate shell contains the so-called dynamical exponent z which describes the dynamical (temporal) critical behavior. This intermediate shell is already present in models for equilibrium dynamical criticality. Crucially, it extends the static critical behavior but does not modify it. In fact there is a certain dynamical fine structure: The same static universality class splits up into various dynamical universality classes, classified in models A to J by HH. 32 Again, we find the dynamic exponents to coincide with the one of an ab initio computation for one of HH's models (MA) -the non-equilibrium conditions do not modify the dynamical critical behavior either. A stronger physical consequence of this finding is discussed in the next subsection.
The unique element found only in the driven system is the outer shell of the aforementioned hierarchy. The related exponent η r identified in Ref. 19, which we refer to as the "drive exponent", physically describes universal decoherence of the long-wavelength dynamics as explained above. Crucially, η r relates to the dynamical MA in the same way as MA relates to the classical O(2) model: It adds a new shell, but does not "feed back" or modify the inner shells of the hierarchy. In Sec. VI we argue that this exponent manifestly witnesses nonequilibrium conditions.
Independence of the drive exponent and maximality of the extension -It is important to demonstrate the independence of the drive exponent: At a second order phase transition, many critical exponents can be defined, each characterizing a different observable. However, only few of them are independent, i.e., cannot be expressed in terms of a smaller set by means of scaling relations.
The independence of the four critical exponents identified with our FRG approach is manifest in the block diagonal structure of the linearized RG flow in the vicinity of the Wilson-Fisher fixed point, cf. Sec. VIII: There are two blocks, and the lowest eigenvalue of each of them determines an independent critical exponent. In addition we have the independent anomalous dimension η and the dynamical exponent z.
A general way to determine the number of independent exponents and thereby see the need for one and only one additional exponent in this system (as compared to equilibrium MA dynamics) comes from the UV limit of the problem. Any independent critical exponent must be related to a shortdistance mass scale in the problem. 68 For example, this can be seen in the case of the anomalous dimension associated with the spatial two-point correlation function. An anomalous dimension η implies decay of the correlation function as φ * (x)φ(0) ∼ |x| 2−d+η . Since the physical units of this correlation are [L] 2−d we require a microscopic scale a, to fix the units so that φ * (x)φ(0) ∼ a −η |x| 2−d+η ∼ [L] 2−d . In the same way any non-trivial independent exponent requires such a microscopic scale.
To determine the number of independent critical exponents in our problem we therefore need to count the microscopic mass scales in the bare action. The corresponding quadratic part of the action reads with real parameters µ, κ 1 , γ, f . κ 1 and f , which describe the tuning parameter of the phase transition and an external ordering field respectively, have direct counterparts in the equilibrium O(2) model. They give rise to the two critical exponents ν, which characterizes the divergence of the correlation length, and η, the anomalous dimension of the static two-point function. γ is introduced in the theory of dynamical critical phenomena and is associated to the dynamical exponent in the purely relaxational MA of HH. 32 In the full non-equilibrium problem however, there is yet another mass scale µ. This scale is at the origin of the additional independent exponent identified in Ref. 19. From this discussion we conclude that the extension of the critical behavior at the condensation transition is maximal, i.e., no more independent exponents can exist. This is due to general requirements on the mass matrix above: the offdiagonal elements must be hermitian conjugates; the lower diagonal must be anti-hermitian; and the upper diagonal must be zero due to the conservation of probability.
It is worth noting how this analysis would change if the critical point in question involved breaking of a Z 2 symmetry rather than a continuous O(N) symmetry as we discuss here. Such an Ising transition in a driven system is relevant for the formation of a super-solid due to interaction of a BEC with the modes of an optical cavity. 2,3 In this case the reality of the Ising fields rules out an imaginary mass term (κ 1 = 0). Hence the maximal number of independent critical exponents is 3, which implies that there can be no modification of MA dynamics.
Interpretation and observability of the drive exponent -The drive critical exponent describes the universal flow behavior of all possible ratios of coherent vs. dissipative couplings (real vs. imaginary parts, see Sec. VIII A) to zero upon moving to larger and larger distances. In the competition of coherent and dissipative dynamics, loosely speaking dissipation always wins. Physically, this should be interpreted as a universal mechanism of decoherence. The drive exponent therefore is subleading and not observable in the correlation functions of the system. However, it is directly observable in the single particle dynamical response (single particle retarded Green's function).
The dynamical response can be measured with any probe that couples directly to the field operatorψ(x), i.e., any probe that out-couples single particles from the system. This is the case, e.g., in the angle-resolved detection of leakage photons in exciton-polariton systems, 20 or in angle-resolved radiofrequency spectroscopy in ultra-cold atoms. 69 As we argue in Sec. VIII C, the excitation spectrum close to the critical point is given by (q = |q|) where A 0 and D 0 are non-universal constants. This excitation spectrum leads to a broadened signal in the experiment, cf. Fig. 3. The drive exponent η r can be observed due to the different scaling with momentum of the location (∼ q z−η r ) and the width (∼ q z ) of the measured peak. We note here, however, that technical noise and other uncertainties in the measurement setup will unavoidably also lead to a broadening of the spectrum. The small value of η r = −0.101 thus challenges experiments to verify this prediction.

III. THE MODEL
In this section we introduce a generic microscopic description of driven-dissipative Bose systems, written in terms of a second quantized master equation. We then show how to translate this model into the Keldysh functional integral framework, which provides a convenient starting point for obtaining the long wavelength universal properties of the system. Moreover we introduce the concept of the effective action, which generalizes the action principle to include all quantum and statistical fluctuations and is the key object for the formulation of the FRG. energy momentum absorption Figure 3. (Color online) Illustration of the observability of the drive exponent: The absorption peak for a measurement that observes the single particle dynamical response. The drive exponent η r reveals itself in a different scaling of the peak location and peak-width as a function of the momentum.

A. Quantum Master equation
Our model with particle loss and pumping is described microscopically by a many-body master equation that determines the time evolution of the system density operator (units are chosen such that = 1), This equation incorporates both coherent dynamics generated by the HamiltonianĤ and dissipation that is subsumed in the action of the Liouville operator L. The HamiltonianĤ describes interacting bosonic degrees of freedom of mass m and is given by (we use the shorthand x = d d x) whereψ(x) are bosonic field operators. Note that we do not explicitly introduce any system chemical potential, as the density of the system will be fixed by the balance of pumping and losses. Two-body interactions are described by a densitydensity interaction with coupling constant g. In the following we shall be interested in dynamically stable systems which are characterized by a positive coupling constant g > 0. This modeling of interactions is valid on length scales which are not sufficient to resolve details of the microscopic interaction potential. In our model, dissipative dynamics comes in the form of one-body pumping (p) and losses (l) as well as two-body losses (t). Accordingly, the Liouville operator can be decomposed into the sum of three terms L = α L α with α = p, l, t which have the common Lindblad structure with local Lindblad or quantum jump operatorsL α (x) that create (p) and destroy (l) single particles; for α = t two particles are destroyed at the same instant in time, i.e., the quantum jump operators are given bŷ These processes occur at rates γ p , γ l , and γ t , respectively. The net effect of single-particle pumping and losses is determined by the relative size of the respective rates: For γ p > γ l , there is an effective gain of single particles. Nevertheless, Eq. (6) leads (in a suitably chosen rotating frame, as we will show below) to a stationary stateρ ss in which the gain of single particles is balanced by two-body losses. In this situation, a finite condensate amplitude builds up, ψ (x) ss = tr ψ (x)ρ ss = ψ 0 0, γ p > γ l .
That is, in stationary state the system is in a condensed phase in which the symmetry of the dynamics described by Eq. (6) under global U(1) transformations of the field operatorsψ(x) →ψ(x)e iφ is broken. When the loss rate γ l exceeds the pumping rate γ p , on the other hand, no condensate emerges in stationary state, and the expectation value of the bosonic field operator is zero, Equations (10) and (11) can be derived from the master equation (6) in mean-field approximation by making the ansatz of a coherent stationary stateρ ψ = |ψ ψ|, where we assume that the amplitude in |ψ = 1 N exp ψ xψ † (x) |0 is spatially homogeneous but possibly time-dependent. Proper normalization of the coherent state is ensured by the choice N = e V|ψ| 2 with the system volume V. The time-dependence of the condensate amplitude is determined by taking the time derivative on both sides of the equality ψ = tr ψ (x)ρ ψ and using the master equation (6), which results in For γ p > γ l this equation allows for a solution of the form ψ = ψ 0 e −iµt , where the condensate density is determined by the imaginary part of the term in brackets on the right-hand side (RHS) as The parameter µ is then given by µ = g |ψ 0 | 2 . We obtain the steady state density matrix of Eq. (10) by means of a transformation to a rotating frame with the unitary opera-torÛ = exp iµNt , where the particle number operator iŝ N = xψ † (x)ψ(x): We haveρ ss =Ûρ ψÛ † , which is indeed time-independent, and recover Eqs. (10) and (11). Under the transformation to this rotating frame, the Hamiltonian acquires a contribution −µN, whereas the Liouvillian L remains invariant. In the following we will always be working in the rotating frame.
In summary, the steady state phase diagram of our model exhibits two phases: A symmetric one characterized by Eq. (11) and an ordered one where the global U(1) symmetry is broken by a finite condensate amplitude Eq. (10) with definite phase. These two phases are separated by a continuous phase transition with order parameter ψ 0 . The transition is crossed by tuning the single-particle pumping rate from γ p < γ l in the "symmetric" to γ p > γ l in the "symmetrybroken" or "ordered" phase.
In the following we shall be interested in the critical behavior that is induced by tuning γ p − γ l to zero. Powerful tools for investigating critical phenomena at a second order phase transition are provided by a multitude of variants of the RG. The particular flavor we employ here is the FRG in the formulation of Wetterich 56 (for reviews see Refs. [70][71][72][73][74][75], which builds upon the use of functional integrals. Therefore, as a first step towards implementing a FRG investigation of our model, we will reformulate the physics that is encoded in the quantum master equation (6) in terms of Keldysh functional integrals. 76,77

B. Keldysh functional integral
The Keldysh approach provides a means to tackle general non-equilibrium problems in the language of functional integrals. For the model at hand, the dynamics described by the master equation (6) can be represented equivalently as a Keldysh partition function (see App. A): By Ψ σ = ψ σ , ψ * σ T for σ = +, − we denote Nambu spinors of fields on the forward-and backward-branch of the closed time contour, respectively. Then, collecting time and space in a single variable X = (t, x) and using the abbreviation X = dt d d x, the Keldysh partition function reads The fields J σ = j σ , j * σ T are external sources inserted here for the purpose of calculating correlation functions of the bosonic fields in the usual manner by means of functional differentiation. When they are set to zero, J + = J − = 0, the partition function reduces to unity, 76,77 i.e., we have the normalization Z[0, 0] = 1. While the Keldysh approach can in principle be utilized to study time evolution, here we are assuming translational invariance in time, as appropriate for the investigation of steady state properties.
In complete analogy to the separation of coherent and dissipative contributions to the time evolution of the density operator in Eq. (6), the action S in the functional integral Eq. (14) can be decomposed as S = S H + S D into a Hamiltonian part S H and a part S D corresponding to the dissipative Liouvillian L in the master equation. The former is given by (from now on we will use units such that 2m = 1) As a general rule (see App. A), normally ordered operators in Eq. (6) acting on the density matrixρ from the left (right) result in corresponding fields on the σ = + (σ = −) contour. Consequently, the commutator with the Hamiltonian in Eq. (6) is transferred into the two contributions to Eq. (15) with a relative minus sign. The same rule applies to the dissipative part in the master equation (6). Passing from the Liouvillian L on to a dissipative action S D , quantum jump operatorsL α are replaced by corresponding jump fields L α,σ on the σ = + (σ = −) contour. (In App. A we will discuss regularization issues related to normal ordering of Lindblad operators.) As above we have the three contributions S D = α S α that are due to singleparticle pumping (p) and losses (l) as well as two-body losses (t). The form of the jump fields can directly be inferred from Eq. (9) as Then, for the dissipative parts of the action we find the expression As we can see, the transition from a description of a specific problem in terms of a master equation to one in terms of Keldysh functional integrals reduces to the application of simple rules. For our model, Eqs. (14), (15) and (17) provide us with a convenient starting point for the investigation of the steady state phase transition described in the previous section. While the translation rules from the master equation to the Keldysh functional integral are most simply applied in a basis of fields ψ ± that can be ascribed to the forward and backward branches of the Keldysh contour, subsequently we will find it advantageous to introduce so-called classical and quantum fields, given by the symmetric and anti-symmetric combinations Condensation is described by a time-independent, homogeneous expectation value of the fields on the σ = ± contours, ψ + (X) = ψ − (X) = ψ 0 , cf. Eq. (10). In the basis of classical and quantum fields, this is expressed as φ c (X) = φ 0 = √ 2ψ 0 , φ q (X) = 0, i.e., only φ c can condense (and, therefore, become a "classical" variable), whereas φ q is a purely fluctuating field with zero expectation value by construction. By means of the transformation Eq. (18), the inverse propagator, determined by the quadratic part of the action, is cast in the characteristic causality structure 76,77 with retarded, advanced, and Keldysh components P R , P A , and P K , respectively (in the following we will denote the two-body coupling constant and loss rate by, respectively, λ = g/2 and κ = γ t /2), The inverse retarded and advanced single-particle Green's functions are given by P R = P A † = i∂ t + ∆ + µ + iκ 1 where κ 1 = γ l − γ p /2. For the Keldysh component of the inverse propagator we have P K = iγ, where γ = γ l + γ p is the sum of single-particle pumping and loss rates -both of them increase the noise level in the system. The spectrum of single-particle excitations is encoded in the poles of the retarded propagator in frequency-momentum space or, equivalently, in the zeros of the inverse propagator. Solving P R (Q) = 0 for ω, where Q = (ω, q) collects the frequency and spatial momentum, we obtain the dispersion relation For κ 1 > 0 (i.e., γ p < γ l ) the pole is located in the lower complex half-plane, and the effective loss rate κ 1 takes the role of an inverse lifetime. One has single-particle excitations that decay exponentially in time, a situation that is well-known from the general theory of the analytic structure of correlation functions. 78 As κ 1 is tuned to negative values (i.e., as we cross the phase transition), however, the pole Eq. (20) is shifted into the upper complex half-plane, signaling an instability. After crossing this threshold, the system develops a condensate, and the proper analytical structure of the retarded propagator is restored only by taking the tree-level shifts due to the condensate into account. We will discuss the corresponding modifications of the dispersion relation Eq. (20) below in Sec. IV A. Inversion of the 2 × 2 matrix in Eq. (19) yields the propagator with retarded, advanced, and Keldysh components, The components along with their respective usual diagrammatic representation 76,77 are given by which shows that the poles of G(Q) are determined solely by the zeros of P R (Q) and P A (Q). The Keldysh component P K of the inverse propagator enters the expression for G K (Q) multiplicatively. Therefore, even in a situation where P K is a polynomial in frequency and/or momentum, it can not give rise to further poles in the propagator G(Q).
In the Keldysh formalism elastic two-body collisions and two-body losses are treated on an equal footing: Both appear in the action Eq. (19) as quartic vertices, however, with a real coupling constant λ in the case of elastic collisions and a purely imaginary coupling constant iκ for two-body losses. The vertices in Eq. (19) can further be distinguished by the number of quantum fields they contain: We have the so-called classical vertex − X (λ + iκ) φ * 2 c φ c φ q + c.c. which contains only one quantum field and three classical fields, and two quantum vertices: The first one i4κ X φ * c φ c φ * q φ q containing two and the second one − X (λ + iκ) φ * 2 q φ c φ q + c.c. containing three quantum fields. Diagrammatically, these vertices are depicted as The fact that there are no vertices consisting only of classical fields is a manifestation of causality or conservation of probability in the Keldysh framework. 76,77 Below we will find that only the classical vertex is relevant (in the sense of the RG) once the system is tuned close to the phase transition.

C. Effective action
Having established a description of our model in terms of a Keldysh functional integral, we proceed by introducing the concept of the effective action, [79][80][81] which is central to the FRG. It is also a convenient starting point for a discussion of the phase transition on the mean-field level (see Sec. IV A).
In equilibrium statistical physics the effective action Γ is related to the free energy as a functional of a space-dependent order parameter, and the equilibrium state is determined as the order parameter configuration that minimizes Γ. However, in the present context of non-equilibrium statistical physics we do not have a sensible notion of a free energy. In fact, already the Keldysh partition function Eq. (14) reduces for vanishing external sources to a representation of unity Z[0, 0] = 1, independently of the parameters that characterize the action. 76,77 Still, the Keldysh effective action, defined analogously to its equilibrium counterpart as the Legendre transform of the generating functional for connected correlation functions, is a very useful object. From Γ we can derive, e.g., field equations that determine the stationary configurations of classical and quantum fields Φ ν = φ ν , φ * ν T , ν = c, q. On a more formal level, Γ is the generating functional of one-particle irreducible vertices. [79][80][81] Most importantly for our model, however, the FRG provides us with a means of calculating critical exponents for the phase transition by studying the RG flow of Γ as a function of an infrared cutoff k.
Our starting point for introducing the effective action is the generating functional Eq. (14) for correlation functions, expressed in the basis of classical and quantum fields Φ ν , i.e., the action is given by Eq. (19), and we introduce classical and quantum sources J ν = j ν , j * ν T with ν = c, q according to the Keldysh rotation where the matrix M is defined in Eq. (18). For the generating functional W for connected correlation functions and Z we have the relation The idea is now to express W, which is a functional of the external sources J ν , in terms of the corresponding field expectation valuesΦ ν = Φ ν | J c ,J q = δW/δJ ν where ν = q for ν = c and vice versa. Introducing these as new variables is accomplished by means of a Legendre transform: The difference between the in this way defined effective action Γ and the action S consists in the inclusion of both statistical and quantum fluctuations in the former. This becomes apparent in the representation of Γ as a functional integral, 70 which holds for the equilibrium states that obey δΓ/δΦ c = δΓ/δΦ q = 0 at vanishing external sources J c = J q = 0. The most straightforward way of evaluating the functional integral Eq. (27) approximately is by performing a perturbative expansion around the configuration that minimizes the action S. To zeroth order this corresponds to mean-field theory, an approach we will discuss in the following section. In the FRG, the fluctuations δΦ ν are included stepwise by introducing an infrared regulator which suppresses fluctuations with momenta less than an infrared cutoff scale k. A short review of this method, adapted to the Keldysh framework, is provided in Sec. V A. We will apply it to our model in Sec. VII.

IV. PREPARATORY ANALYSIS
Here we carry out a basic analysis of the model in preparation for setting up a full functional RG calculation used to obtain the critical properties at the phase transition. We summarize the mean-field theory for the effective action and discuss the generic emergence of infrared divergences near a critical point. Furthermore, using dimensional analysis we identify the important terms in the action which are potentially relevant at the critical point. These terms are then included in the ansatz of the effective action used to carry out the FRG calculation. Finally we contrast this ansatz with the equilibrium dynamical models of HH. 32

A. Mean-field theory
In Sec. III A we identified the precise balance between single-particle losses and pumping as the transition point, cf. Eqs. (10) and (11). Here we will derive this result from the Keldysh functional integral Eq. (27), again employing a mean-field approximation. We will then proceed by calculating the excitation spectrum above the stationary mean-field by treating quadratic fluctuations in a Bogoliubov (tree-level) expansion. While this issue, as well as going beyond the meanfield approximation by perturbative methods, can equally well be addressed in the master equation formalism of Sec. III A, 82 in performing a perturbative expansion at and below the critical point we encounter infrared divergences. Proper treatment of these requires RG methods, which are well-established and elegantly formulated in terms of functional integrals.
Mean-field theory corresponds to a saddle-point approximation of the functional integral in Eq. (27) in which fluctuations around the classical fields are completely neglected. In the present context, by classical fields we mean spatially homogeneous solutions to the classical field equations As already mentioned above, there are no terms in the action Eq. (19) that have zero power of both φ * q and φ q , and the same is obviously true for δS/δφ * c . Therefore, the first equation (28) is solved by φ q = 0. Inserting this condition in the second equation (28), we have The solution φ c = φ 0 is determined by the imaginary part of Eq. (29): For κ 1 ≥ 0, in the symmetric phase, the classical field is zero, ρ 0 = |φ 0 | 2 = 0, whereas for κ 1 < 0 we have a finite condensate density ρ 0 = −κ 1 /κ. In a second step, the parameter µ is determined by the real part of Eq. (29) as µ = −λκ 1 /κ. Quadratic fluctuations around the mean-field order parameter can be investigated in a Bogoliubov or tree-level expansion: We set φ c = φ 0 + δφ c , φ q = δφ q in the action Eq. (19) and expand the resulting expression to second order in the fluctuations δφ ν . The poles of the retarded propagator (which is now a 2 × 2 matrix in the space of Nambu spinors Real and imaginary parts of both branches are shown in Fig. 4 in panels (a) and (b), respectively. Due to the tree-level shifts ∝ ρ 0 the instability of Eq. (20) for κ 1 < 0 is lifted: Both poles are located in the lower complex half-plane, indicating a physically stable situation with decaying single-particle excitations. For κ = 0, Eq. (30) reduces to the usual Bogoliubov result, 83 where for q → 0 the dispersion is phononic, ω R 1,2 = ±cq, with speed of sound c = 2λρ 0 whereas particlelike behavior ω R 1,2 ∼ q 2 is recovered at high momenta. Here, due to the presence of two-body loss κ 0, the dispersion is strongly modified: While at high momenta the dominant behavior is still given by ω R 1,2 ∼ q 2 , at low momenta we obtain purely dissipative non-propagating modes ω R 1 ∼ −i λ κ q 2 and ω R 2 ∼ −i2κρ 0 . In particular, for q = 0 we have ω R 1 = 0: This is a dissipative Goldstone mode, 27,58,60 associated with the spontaneous breaking of the global U(1) symmetry in the ordered phase. The existence of such a mode is not bound to the meanfield approximation but rather an exact property of the theory guaranteed by the U(1) invariance of the effective action, even in the present case of a driven-dissipative condensate.

B. Infrared divergences near criticality
The discussion of our model on the mean-field level has illustrated some of the benefits of the Keldysh approach: Not only have we gained a simple physical picture of the phase transition as a condensation instability in the retarded and advanced propagators, but we were able to investigate excitations in both the symmetric and ordered phases quite straightforwardly. Mean-field theory, however, while providing us with a good qualitative understanding of the stationary state physics of our model far away from the phase transition, has major shortcomings when it comes to the discussion of critical phenomena. In particular, the critical exponents that can be extracted from an analysis of quadratic fluctuations around the mean-field configuration are not indicative of the universality class of the phase transition, as they correspond to the RG flow in the vicinity of a non-interacting (or Gaussian) fixed point. Critical behavior at the phase transition, however, is encoded in the RG flow in the vicinity of an interacting (or Wilson-Fisher) fixed point.
In a many-body system, excitations and their interactions get dressed due to scattering from other particles. The meanfield results of this section can be taken as the starting point for a calculation of the effective dressed parameters in a perturbative expansion. In the functional integral Eq. (27), diagrammatically this amounts to an expansion in the number of loops around the mean-field configuration. To lowest (one-loop) order, the correction ∆λ to the real part of the bare classical vertex (the first diagram in Eq. (23)) reads where the elements appearing in the diagram are defined in Eqs. (31) and (22) (here, however, lines correspond to propagators of fluctuations δΦ ν and acquire an additional 2 × 2 matrix structure in Nambu space), and the ellipsis indicates that all diagrams with four external legs and one closed loop corresponding to a single internal momentum integration have to be included. In the integrand we have only kept the dominant contribution for q → 0, and we have introduced an infrared cutoff q IR in order to regularize the divergence at low momenta. Such infrared divergences, however, appear not only in our specific example of the loop correction to λ, but rather are characteristic of perturbative expansions in symmetry broken phases. They are due to the presence of a massless Goldstone mode, which results in a pole of the retarded and advanced propagators at ω = q = 0. This problem is even enhanced as we approach the phase transition: Then both modes become degenerate, with also the second mode ω R 2 ∼ −i2κρ 0 for q → 0 becoming massless. A method that allows us to go beyond mean-field theory, therefore, has to provide for a proper treatment of infrared divergences. In the FRG, this is achieved by effectively introducing a mass term ∝ k 2 in the inverse propagators by hand. In consequence, the integrand in Eq. (31) is replaced by ∞ q IR dqq d−1 / q 2 + k 2 2 and we may safely set q IR to zero since the effective mass k 2 acts as an in-frared cutoff. The resulting loop-corrected coupling is a function of this cutoff, λ = λ(k), and we obtain the fully dressed or renormalized coupling by following the RG flow of the running coupling λ(k) for k → 0. This procedure can be implemented efficiently by introducing the cutoff in the functional integral Eq. (27). We will discuss how this is done in practice for the present non-equilibrium problem [84][85][86][87][88][89][90][91] in the following section. Critical exponents can then be extracted from the flow of the critical system, i.e., when κ 1 is fine-tuned to zero.
So far we have discussed only corrections to the bare interaction vertices due to the inclusion of loop diagrams. However, also the propagators appearing in these diagrams are themselves renormalized. In particular, the inverse propagator can be written as P(Q) − Σ(Q), i.e., as the sum of the bare inverse propagator P(Q) and a self-energy correction Σ(Q). [79][80][81] The self-energy contribution at one-loop order to the retarded propagator is represented diagrammatically as where effective cubic couplings, which are obtained upon expanding the interaction vertex around the field expectation value, appear in the second diagram. Lines beginning and terminating in crosses indicate that particles are scattered out of and into the condensate, respectively. Due to momentum conservation, the first diagram does not depend on the external momentum Q = (ω, q) and gives a correction to the constant part of the inverse propagator, i.e., the so-called mass terms. Since the coupling λ + iκ associated with the vertex appearing in this diagram is complex, both the real and imaginary masses, µ and κ 1 , are affected by the loop correction. The second diagram in Eq. (32) gives a frequency-and momentum-dependent contribution to the selfenergy. Symmetry under spatial rotations implies that it depends only on the modulus of the momentum and we may write Σ R (Q) = Σ R (ω, q 2 ). For small ω and q 2 we can expand Σ R (ω, q 2 ) ≈ Σ R (0, 0)+ω∂ ω Σ R (0, 0)+q 2 ∂ q 2 Σ R (0, 0). Transforming back to the time domain and real space, the derivatives of the self-energy with respect to frequency and momentum give corrections to the coefficients of ∂ t and ∆ in the inverse propagators, which are again complex valued. An imaginary part of the coefficient of the Laplacian corresponds to an effective dissipative kinetic coefficient due to the interaction with other particles; A complex prefactor of the time derivative, on the other hand, has significant consequences for the physical interpretation of all other couplings, as we will discuss in detail in later sections.

C. Canonical power counting
While the proper theoretical approach to critical phenomena has to cope efficiently with the infrared divergences discussed above, such systems also exhibit an important ordering principle, which is provided by the classification of couplings according to their canonical scaling dimension. In the following we will briefly review this procedure, often referred to as canonical power counting or dimensional analysis. It lays the basis for a suitable choice of ansatz for the effective action that will contain only couplings which are relevant or marginal according to this counting scheme. [79][80][81] At second order phase transitions, physical quantities exhibit scaling behavior, which means that they depend on the distance from the phase transition (in our case this distance is measured by κ 1 ) in a power-law fashion ∼ κ τ 1 , with a generally non-integer exponent τ. In order to study critical behavior in the RG, we investigate the RG flow starting from the action fine-tuned to criticality, i.e., with κ 1 = 0, and approach the critical point by lowering k. Then, scaling behavior of a physical quantity g shows up as power-law dependence g ∼ k θ on k for k → 0 with a critical exponent θ. In other words, phase transitions are associated to scaling solutions of the RG flow (not all scaling solutions correspond to phase transitions 92 ), or -equivalently -fixed points of the flow of rescaled couplingsg = k −θ g. The dominant contribution to the exponent θ associated to a coupling g is determined by its physical dimension measured in units of momentum k, i.e., the canonical scaling dimension or engineering dimension [g] (we have [k] = 1). Anticipating that deviations from canonical scaling will be small (see Sec. VIII), let us study the flow of the dimensionless two-body elastic collision couplingλ = λ/k (we will see below thatλ is indeed dimensionless). In Sec. IV A we saw that the flow of λ is generated by the loop diagrams Eq. (31). Then, to the flow of the dimensionless variableλ we have an additional contribution due to the engineering dimension, where we are taking the derivative with respect to the dimensionless logarithmic scale t = ln(k/Λ) which is zero for k = Λ and goes to −∞ for k → 0. The loop contribution to the flow ofλ is of orderλ 2 ,λκ,κ 2 and higher in the dimensionless twobody couplingsλ,κ. We find, therefore, a trivial fixed point ∂ tλ = 0 forλ * =κ * = 0. The flow for smallλ in the vicinity of this Gaussian fixed point is determined by the canonical scaling contribution on the RHS of Eq. (33) and is directed towards higher values ofλ, i.e., the couplingλ is relevant at the Gaussian fixed point. For increasingλ, the loop contributions become important and balance canonical scaling at a second fixed point. This non-trivial Wilson-Fisher fixed point at finiteλ * ,κ * corresponds to the phase transition in the interacting system, and for small deviationsλ −λ * the flow is attracted toλ * . The described scenario changes drastically for a coupling with negative canonical scaling dimension, i.e., when instead of the prefactor −1 for the first term on the RHS in Eq. (33) we had a positive integer. Such a coupling is irrelevant at the Gaussian fixed point, which means that its flow is attracted to that fixed point. We can, therefore, as a starting point for a systematic expansion in the relevance of couplings, set all irrelevant couplings to zero. Unlike perturbative expansions, the inclusion of irrelevant couplings in higher orders in the expansion in canonical scaling dimensions results not only in enhanced quantitative accuracy, but rather refines our picture of the phase transition, as it involves higher order ver-tices and a refined treatment of the momentum dependence of propagators. 93 We proceed by determining the canonical dimensions of the couplings appearing in the action Eq. (19). They are not uniquely fixed by the requirement that the action is dimensionless, [S] = 0: Still we have the freedom of assigning different scaling dimensions to the classical φ c and quantum fields φ q . We exploit this freedom in order to impose a scaling dimension upon the Keldysh component of the inverse propagator in Eq. (19) that is the same as in finite-temperature thermodynamic equilibrium, 76,77 i.e., we require [γ] = 0. While this choice yields a consistent picture of the driven-dissipative Bose condensation transition as detailed below, it is inappropriate for the investigation of stationary transport solutions that define genuine nonequilibrium states with nonvanishing flux which might be contained in our model. As already pointed out in Sec. II, in two dimensions, such a scenario indeed has been recently identified in Ref. 66, showing that the Kardar-Parisi-Zhang non-equilibrium fixed point 43 governs the long wavelength behavior. In three dimensions, a similar scenario is conceivable in principle, however only beyond a certain threshold value for the strength of violation of detailed balance.
Denoting the dynamical exponent by [∂ t ] = z we find, from the quadratic part of the action and in d dimensions, The different scaling dimensions of classical and quantum fields result in different behavior of the complex couplings associated with the classical and quantum vertices Eq. (23) under renormalization, even though their values at k = Λ are the same. In particular, for a local vertex that contains n c classical and n q quantum fields, the canonical scaling dimension of the corresponding coupling is We observe that all couplings λ n c ,n q with n q > 2 (n q ≥ 1 is required by causality [76][77][78] or n c > 5 are irrelevant. The coupling λ 3,1 associated with the classical quartic vertex has canonical dimension 4 − d, i.e., its upper critical dimension is d = 4 and, in the case of interest d = 3, it is relevant with canonical scaling dimension equal to unity. All other quartic couplings are irrelevant, as are sextic couplings with n q > 1. The classical three-body coupling λ 5,1 is marginal and we will include it (with both real and imaginary parts) in our ansatz for the running effective action below, even though it is not present in the action S. Higher order couplings λ n c ,n q with n c + n q > 6 are irrelevant and we will discard them.

D. Equilibrium symmetry
According to the canonical power counting scheme outlined in the previous section, in order to describe critical properties at the driven-dissipative Bose condensation transition we may disregard quantum vertices in the action Eq. (19) -this corresponds to a semiclassical approximation, [76][77][78]94 and the resulting simplified "mesoscopic" action has the same structure as the classical dynamical models considered in Ref. 32 inasmuch as it is linear in the quantum fields apart from the noise term which is quadratic. Therefore, like in the classical dynamical models, the functional integral with the mesoscopic action is equivalent to a Langevin equation for the classical field. This is just the stochastic dissipative Gross-Pitaevskii equation Eq. (1) for a single non-conserved complex field ψ and bears close resemblance to the equation of motion of MA of HH with N = 2 real components. There are, however, two key differences: First the dynamics in MA is purely relaxational whereas Eq. (1) contains both coherent and dissipative contributions. Second, and more importantly, dropping all coherent contributions on the RHS of Eq. (1) we find that it is invariant under the transformation of the fields: 91,95,96 This symmetry of the dynamics implies a FDT for the retarded response and correlation functions. Its absence in the driven-dissipative model (DDM), therefore, may be seen as indicating non-equilibrium conditions. In Sec. VI below we discuss a generalized version of the symmetry transformation Eq. (36) and we are led to consider an extension of MA by coherent dynamics that then differs from the DDM precisely in the obedience to this generalized symmetry. With regard to critical phenomena, the difference in symmetries between equilibrium and non-equilibrium situations renders it possible that novel universal behavior may be found in the latter case. We proceed to perform an FRG analysis of the critical properties of both models in the following sections.

V. FUNCTIONAL RENORMALIZATION GROUP
A. FRG approach for the Keldysh effective action The transition from the action S to the effective action Γ consists in the inclusion of both statistical and quantum fluctuations in the latter (see Eq. (27)). In the FRG, the functional integral over fluctuations is carried out stepwise by introducing an infrared regulator which suppresses fluctuations with momenta less than an infrared cutoff scale k. 70 This is achieved by adding to the action in Eq. (14) a term with a cutoff function R k,K which will be specified below in Sec. V B. We denote the resulting cutoff-dependent Keldysh partition function and generating functional for connected correlation functions by, respectively, Z k and W k . The effective running action Γ k is then defined as the modified Legendre Here the subtraction of ∆S k on the RHS guarantees that the only difference between the functional integral representations for Γ and Γ k is the inclusion of the cutoff term in the latter, Physically, Γ k can be viewed as the effective action for averages of fields over a coarse-graining volume with size ∼ k −d .
We choose the form of the cutoff term ∆S k such that it modifies the inverse retarded and advanced propagators: Comparing Eqs. (19) and (37), we see that associated with the action S + ∆S k are the regularized retarded and advanced inverse propagators P R (Q) + R * k,K (q 2 ) and P A (Q) + R k,K (q 2 ) respectively, whereas the Keldysh part P K of the inverse propagator remains unchanged. In other words, by introducing the cutoff ∆S k we manipulate the spectrum of single-particle excitations, which is encoded in the zeros of the inverse propagators P R/A (Q) or, equivalently, in the poles of the propagators Eq. (22). At the transition, these poles are determined by Eq. (20) with κ 1 = 0, i.e., we have a pole at ω = q = 0, and as we have pointed out in the paragraph following Eq. (31), this leads to infrared divergences that drive critical behavior. For the regularized propagators, on the other hand, we have To regulate infrared divergences, it is sufficient to introduce the cutoff function in the retarded and advanced inverse propagators, as becomes clear from the discussion following Eq. (22). We have seen that the effective action Γ k defined by Eq. (39) has an infrared-finite loop expansion. Its main usefulness, however, lies in the fact that it interpolates between the action S for k → Λ where Λ is an ultraviolet cutoff scale, and the full effective action Γ for k → 0. This is ensured by the requirements on the cutoff function 85 where under the condition that Λ exceeds all energy scales in the action by far, for k → Λ we may evaluate the functional integral Eq. (39) in a stationary phase approximation. Then, to leading order we find Γ Λ ∼ S. The evolution of Γ k from this starting point in the ultraviolet to the full effective action in the infrared for k → 0 is described by the exact Wetterich flow equation 56,70 where Γ (2) k andR k denote, respectively the second variations of the effective action and the cutoff ∆S k and will be specified in Sec. V B; Tr denotes summation over internal field degrees of freedom as well as integration over frequencies and momenta. The flow equation provides us with an alternative but fully equivalent formulation of the functional integral Eq. (39) as a functional differential equation. Like the functional integral, the flow equation can not be solved exactly. It is, however, amenable to various systematic approximation strategies. Here we perform an expansion of the effective action Γ k in canonical scaling dimensions as outlined above in Sec. IV C, keeping only those couplings which are -in the sense of the RG -relevant or marginal at the phase transition.

B. Truncation
In three dimensional classical O(N)-symmetric models, already the inclusion of non-irrelevant couplings gives a satisfactory description of critical phenomena. 70 As we will show below, static critical properties of our non-equilibrium phase transition are described by such a model with N = 2. Therefore, in the following, we will as well restrict ourselves to the inclusion of relevant and marginal couplings in the ansatz for the effective action, i.e., we choose a truncation of the form (Here all couplings depend on the infrared cutoff scale k. However, for the sake of keeping the notation simple, we will not state this dependence explicitly.) All terms involving derivatives are contained inD R = iZ * ∂ t +K * ∆ andD A =D R † . In contrast to the action Eq. (19), however, we allow for complex coefficients Z = Z R + iZ I andK =Ā + iD: Due to the presence of complex couplings λ + iκ in the classical action, imaginary parts of Z andK will be generated in the RG flow as indicated at the end of Sec. IV A, even though they are zero initially at k → Λ.
A complex prefactor Z of the time derivative -often referred to as wave-function renormalization -obscures the physical interpretation of the other complex couplings: The field equation δΓ k /δφ * q = 0 contains iZ * ∂ tφc = −K * ∆φ c + · · · . The physical meaning of the gradient coefficientK becomes clear only after division by Z * , i.e., in the form i∂ tφc = − (A − iD) ∆φ c + · · · where we introduced the decomposition K =K/Z = A + iD into real and imaginary parts. In this form, the interpretation of A and D as encoding coherent propagation and diffusive behavior of particles is apparent. Similar considerations hold for the other couplings in Eq. (43), and we will elaborate on this point in Sec. VI D.
In our truncation containing only non-irrelevant contributions, the only momentum-independent couplings we keep are the Keldysh and spectral masses,γ andū 1 = −μ + iκ 1 , as well as the classical quartic and sextic couplings (i.e., those vertices containing only one quantum field but three and five classical field variables respectively). These are included in the part in Eq. (43) that involves the potentialŪ, which is a function of the U(1) invariantρ c = φ c 2 and given bȳ where bothū 2 =λ + iκ andū 3 =λ 3 + iκ 3 are complex. In the symmetric phase, we keepū 1 0 as a running coupling and setρ 0 = 0, whereas in the ordered phase we set the masses to zero,ū 1 = 0, and regard the condensate amplitude as a running coupling,ρ 0 0. Then, the parameterization Eq. (44) corresponds to an expansion of the potential around its minimum in both the symmetric and ordered phases. It ensures that the field equations δΓ k /δφ * c = 0, δΓ k /δφ * q = 0 are solved byρ c = 0 andρ c =ρ 0 in the symmetric and ordered phases respectively (in both cases we requireφ q =φ * q = 0) for all values of k.
In what follows we will find it advantageous to introduce renormalized fields φ c =φ c , φ q = Zφ q (the various symbols for bare/renormalized fields etc. are summarized in Tab. I).
With this choice the complex wave-function renormalization Z that multiplies the time derivative in Eq. (43) is absorbed in the field variables and we can write the effective action in the form (σ z denotes the Pauli matrix) The renormalized Keldysh mass is γ =γ/ |Z| 2 . For the variational derivatives with respect to the classical fields we are using the notation δ/δΦ * c = δ/δφ * c , δ/δφ c T , and the renormalized potential functionals that encode unitary and dissipative terms respectively, read where A and D are the real and imaginary parts of the renormalized gradient coefficient K =K/Z = A+iD. Primes denote derivatives with respect to ρ c = |φ c | 2 of the real and imaginary parts of the renormalized potential U =Ū/Z = U H + iU D , which is given by with renormalized couplings u 1 =ū 1 /Z = −µ + iκ 1 , u 2 = u 2 /Z = λ + iκ, and u 3 =ū/Z = λ 3 + iκ 3 . The inclusion of the classical three-body coupling u 3 adds the vertex to the building blocks Eqs. (22) and (31) dimensionless couplings VIII A Table I. Summary of notation. The columns are symbols, their meaning, and the section in which they are introduced.
As we have already indicated, the first variational derivative of the effective action yields field equations that determine the stationary state values of the classical and quantum fields. In the ordered phase, these are constant in space and time and read φ c | ss = φ * c | ss = √ ρ 0 (our choice of a real condensate amplitude does not cause a loss of generality) and φ q | ss = φ * q | ss = 0. Then, the scale-dependent inverse connected propagator is given by the second variational derivative of the effective action, 79-81 evaluated in stationary state. We will carry out this variational derivative in a basis of real fields, which we introduce by decomposing the classical and quantum fields into real and imaginary parts according to φ ν = 1 √ 2 χ ν,1 + iχ ν,2 for ν = c, q. The inverse propagator at the scale k is then given by Here the indices i, j enumerate the four components of the field vector Analogous to the inverse propagator in the action Eq. (19), the inverse propagator at the scale k is structured into retarded, advanced, and Keldysh blocks, However, here these blocks are themselves 2 × 2 matrices. (This additional Nambu structure emerges in the ordered phase.) We have explicitly These expressions can be used to deduce the dispersion relation for single-particle excitations. It is determined by solving det P(Q) = det P R (Q) det P A (Q) = 0 for ω. Due to the second relation Eq. (52), two of the four solutions to Eq. (53) are complex conjugate. The zeros of the determinant of the retarded inverse propagator encode the two branches ω R 1,2 = −iDq 2 − iκρ 0 ± Aq 2 Aq 2 + 2λρ 0 − (κρ 0 ) 2 , (54) which differ from the mean-field expression Eq. (30) by the contribution −iDq 2 due to the explicit inclusion of a dissipative kinetic term in our truncation, and by the appearance of the scale dependent gradient coefficient A. The dissipative Goldstone mode is now characterized by the low-momentum behavior ω R 1 ∼ −i D + A λ κ q 2 , whereas for the "massive" (the mass is purely imaginary) mode we reproduce the form of the mean-field expression ω R 2 ∼ −i2κρ 0 -however, in a scaledependent version with all couplings running in the course of the RG. In this way, structural properties such as Goldstone's theorem are preserved during the flow. The dispersion relation Eq. (54) is depicted in Fig. 4.
We proceed by specifying the cutoff function R k,K which appears in Eq. (37). We will work with an optimized cutoff 97 which obviously meets the requirements Eqs. (40) and (41). The regularized propagator, which appears in the loop diagrams that generate the RG flow, reads where the 4 × 4 matrix R k (Q) is defined in analogy to the inverse propagator Eq. (49) as the second variational derivative of the cutoff Eq. (37) with respect to the real fields Eq. (50), Due to the cutoff R k (Q) in the denominator in Eq. (56), the poles of G k (Q) are given by Eq. (54), however, with Aq 2 and Dq 2 replaced by p A (q 2 ) and p D (q 2 ) respectively, where the function p a (q 2 ) for a = A, D reads p a (q 2 ) = aq 2 − R k,a (q 2 ) =        ak 2 for q 2 < k 2 , aq 2 for q 2 ≥ k 2 .
The thus modified dispersion relations are finite for q → 0, i.e., infrared divergences of loop diagrams are regularized.
In panel (d) in Fig. 4 the regularized dispersion relations are shown as dashed-dotted lines. In Sec. V we introduced most of the ingredients for a FRG investigation of the steady state driven-dissipative Bose condensation transition. Before we present the explicit flow equations in Sec. VII, we will now provide a detailed discussion of the relation between our non-equilibrium model and the classical equilibrium dynamical MA of HH. 32

VI. RELATION TO EQUILIBRIUM DYNAMICAL MODELS
Here we extend the discussion of Sec. IV D and work out the precise relation of the DDM to MA with N = 2 com- ponents. We reemphasize that these considerations rely on the power counting introduced in Sec. IV C, which implies that we may omit quantum vertices from an effective longwavelength description close to criticality; The resulting action Eq. (45) is equivalent to a Langevin equation of the form of Eq. (1). 76,77,94,98 Originally, MA was formulated in terms of such a Langevin equation for a non-conserved, coarse-grained order parameter. It provides for a phenomenological description of the relaxational dynamics of the order parameter subject to stochastic fluctuations, which are introduced necessarily as a consequence of the coarse-graining over a volume of extent k −d cg : The effects of fluctuations with momenta q greater than the coarse-graining scale k cg are included by introducing random noise sources in the evolution equation.
For our model, coarse-graining amounts to integrating out fluctuations with momenta q greater than k cg in the functional integral Eq. (39), 70 which results in an effective action Γ cg that can be regarded as the starting point of a phenomenological description in the spirit of HH, i.e., we may interpret it as the action S cg = Γ cg for slow modes with momenta q < k cg .
The equation of motion of MA is constructed such that its stationary state is thermodynamic equilibrium, which manifests itself in a FDT 32 relating the order parameter retarded response and correlation functions. The FDT can be derived as a consequence of a specific equilibrium symmetry of the dynamics which is related to time reversal and expresses detailed balance. 91,95,96 This symmetry, however, does not restrict the dynamics to be purely relaxational as is the case in MA. In fact, one can conceive an extension of MA by reversible mode couplings (MAR) which differs from the DDM only in the obedience of the symmetry. (Note that the DDM generically features both coherent and dissipative contributions). As universality classes are fully characterized by the spatial dimensionality and symmetries of a system, however, this opens up the possibility of novel critical behavior in the DDM.
In the remainder of this section we illuminate the consequences of the equilibrium symmetry through a detailed comparison between MAR in which it is present at the outset and the DDM model, where it is only emergent at long scales. We give a simple geometric interpretation of the restriction that the symmetry imposes on the couplings that parameterize the effective action and specify the submanifolds in the coupling space for the DDM that correspond to MA and MAR.
While these considerations demonstrate formally the nonequilibrium character of the DDM, the equilibrium MAR constructed in the above way may seem a bit academic. In fact, as we will see in Sec. VI D it amounts to an unrealistic finetuning of the ratios of all coherent vs. dissipative couplings. The physically relevant model which the DDM should be compared to is model E, which describes the equilibrium Bose condensation transition. An important difference between the DDM and model E is the presence of an exact particle number conservation in the latter case which can be seen to rule out a finite κ 1 mass term. 99 Therefore, according to the arguments given in Sec. II, the standard equilibrium Bose condensation transition exhibits only three independent exponents (as opposed to four in the DDM) and, in particular, no counterpart to η r . Moreover, as a consequence of the exact particle number conservation an additional slow mode occurs at criticality and modifies the dynamical exponent.

A. Model A with N = 2 and reversible mode couplings (MAR)
We specify the equilibrium symmetry in terms of fieldsΦ ν which are related to the bare fieldsΦ ν of Eq. (43) viã where Z R,cg and Z I,cg denote the real and imaginary parts of the wave-function renormalization at the coarse-graining scale k cg andr is a real parameter, the physical meaning of which will become clear in the following. The symmetry transformation is denoted by T and reads 91,95,96 TΦ c (t, x) = σ xΦc (−t, x), cf. the implementation in the Langevin formulation Eq. (36). It includes complex conjugation (in the form of multiplication with the Pauli matrix σ x ) and time reversal; T is the temperature. As outlined above, we now construct the action for MAR as follows: We identify the effective action Eq. (43) at the coarse-graining scale k cg with the action for low-momentum modes, S cg = Γ cg , and enforce thermodynamic equilibrium by requiring invariance of S cg under the transformation T , which results in (See App. B for details of the derivation.) The action S MAR cg contains coherent dynamics in the form ofŪ H,cg =rŪ D,cg , i.e., the parameterr plays the role of the common fixed ratio between coherent and dissipative couplings. This relation ensures compatibility of coherent dynamics with the equilibrium symmetry. We note that here, crucially, both the irreversible and the reversible dynamics have the same physical origin, being generated by the same functionalŪ D,cg . This is motivated in the frame of a phenomenological, effective model for relaxation dynamics in the absence of explicit drive. However, not only the values of the couplings encoding coherent dynamics are restricted by the symmetry, but also the Keldysh massγ cg is determined by the temperature that appears in the symmetry transformation as Finally we note that Eq. (61) includes MA with effectively purely dissipative dynamics as a special case: Indeed we can derive the action for MA in the same way as we derived the action for MAR from the truncation for the DDM, i.e., by enforcing an additional symmetry. Requiring invariance of S MAR cg under complex conjugation of the fields, we find the additional constraintr = −Z I,cg /Z R,cg (see App. B), reducing the number of independent parameters further. Then, after rescaling the quantum fields with Z cg it becomes apparent that this model describes purely dissipative dynamics as we will show in Sec. VI D.

B. Truncation for MAR
We proceed by specifying the truncation for a FRG analysis of MAR. Here it is crucial to note that the transformation T Eq. (60) not only leaves the action Eq. (61) invariant, but is actually a symmetry of the full theory, 96 i.e., of the effective action. Then, if the cutoff ∆S k in Eq. (39) is T -invariant as well (this is indeed the case for the choice Eq. (37)), also the scale-dependent effective action Γ MAR k must obey the symmetry. This requirement implies restrictions on the RG flow: Invariance of the effective action on all scales is guaranteed by the ansatz which follows by enforcing the symmetry on the truncation Eq. (43) (see App. B for details). We note in particular that compatibility of coherent and dissipative dynamics is conserved in the RG flow. In contrast to the DDM, here the Keldysh mass is not an independent running coupling, as it is linked to the wave-function renormalization Z = Z R + iZ I by the Ward identity of the symmetry Eq. (60), In comparison to the DDM, therefore, MAR is described by a reduced number of couplings: Our truncation Eq. (43) for the DDM is parameterized by a vector of couplings where Z,K,ū,ū 3 are complex whereasρ 0 ,γ are positive real numbers. In MAR, the real parts of the complex couplings in the functionalŪ are determined by imaginary ones and the ratior which appears as a fixed parameter in the action at the coarse-graining scale k cg . Additionally the Keldysh mass is related to the wave-function renormalization via Eq. (65), so that a reduced set of running couplings, is sufficient to fully specify the truncation Eq. (64). In the purely dissipative MA, finally, the symmetry Eq. (63) determines the ratio of imaginary to reals parts of the wavefunction renormalization Z asr = −Z I /Z R (see App. B), so that Z can be parametrized in terms of a single real running coupling. The truncation for MA, therefore, is described by the couplings:ḡ MA = Z R ,D,ρ 0 ,κ 1 ,κ,κ 3 T .

C. Fluctuation-dissipation theorem
In the following we will show that the symmetry Eq. (60) implies a classical FDT for MAR. 91,95,96 If we regard the full propagators of the theory as the k → 0 limits of the RG flow of scale-dependent propagators, we may say that the FDT holds for MAR (and, a fortiori, for MA) for all 0 < k < k cg . In addition we will see that this is not the case for the drivendissipative system we consider. There the equilibrium symmetry is not present at mesoscopic scales but rather emergent for the system at criticality in the infrared for k → 0. As a result, thermalization sets in only at low frequencies and long wavelengths.
As indicated at the beginning of the preceding section, the transformation T Eq. (60) is a symmetry of the full theory. In particular, for two-point correlation functions we have and corresponding relations hold for higher correlation functions.
Here expectation values are defined as The relation Eq. (69) implies a FDT: For the particular choice of correlations between quantum fields ν = ν = q which vanish by construction of the Keldysh functional integral, [76][77][78] we find Inserting here explicit expressions for the T -transformed fields and performing a Fourier transformation, we obtain the classical FDT Such a relation is in general not valid in the DDM. It is, however, emergent for the critical system in the longwavelength limit: In the basisφ c =φ c ,φ q = i (Z/ |Z|)φ q we have for the inverse propagators at the scale k (for convenience we are working here in the symmetric phase; the scale-dependent inverse propagators are determined by the quadratic part of the effective action Eq. (43))P R (ω, q) = i |Z| (ω − ξ * (q)) =P A (ω, q) † where ξ(q) = Kq 2 + u 1 (note that here the renormalized quantities appear) andP K =P K . With these inverse propagators we form the ratio which would equal the temperature if a FDT were valid. 100 As we will see in Sec. VIII C, the effective action for the critical system becomes purely dissipative for k → 0. In particular we have Re ξ(q) → 0 so that Eq. (73) indeed reduces to an FDT with an effective temperature Note that for purely dissipative dynamics Eq. (65) implies that the ratioγ/ |Z| is a constant of the RG flow. For the DDM the emergence of an FDT with T eff manifests itself in the relation Eq. (114) between the anomalous dimensions ofγ and Z valid at the fixed point. The flow ofγ/ (4 |Z|) is shown in Fig. 2.

D. Geometric interpretation of the equilibrium symmetry
For our truncation of the effective action Γ MAR k , the relation U H =rŪ D between the real and imaginary parts of the func-tionalŪ =Ū H + iŪ D implies that the couplings parameteriz-ingŪ H andŪ D share a common ratior of real to imaginary partsr The same applies to the renormalized couplings, however, with a different value r: With z = −Z I /Z R we have This can be visualized conveniently in the complex plane, where the ratio of real to imaginary parts contains the same information as the argument of a complex number (the argument is tan(1/r)): The renormalization of a complex couplingḡ with Z corresponds to a rescaling |g| = |ḡ| / |Z| of the modulus and a rotation of the phase by the argument of Z, arg g = argḡ − arg Z. The condition Eq. (75) corresponding to MAR is depicted in Fig. 5 (b): All bare 101 couplings lie on a single ray. In the purely dissipative case with r = 0 andr = z, which is shown in Fig. 5 (a), this ray is perpendicular to Z. As a result, in this case the renormalized couplings are purely imaginary. Generally, only the renormalized quantities allow for an immediate physical interpretation: A and D describe propagation and diffusive behavior of particles, respectively, while λ (λ 3 ) and κ (κ 3 ) are two-body (three-body) elastic collisions and loss. In the generic driven-dissipative case, we have no a priori relations between these couplings because they are due to different physical mechanisms: Dissipative couplings describe local incoherent single particle pump and loss, as well as local two-body loss. On the other hand, unitary dynamics is given by coherent propagation and elastic collisions. Geometrically, the physical couplings point in different directions in the first quadrant of the complex plane (see Fig. 5 (c)), the latter restriction being due to the physical stability of the system (see Sec. III A). This concludes our discussion of the relation of the DDM to dynamical equilibrium models. In the following section we will proceed to derive explicit flow equations for the couplings Eq. (66).

VII. NON-EQUILIBRIUM FRG FLOW EQUATIONS
In the following we discuss how the functional differential equation Eq. (42) for the effective action is reduced to a set of ordinary differential equations by virtue of the ansatz Eq. (43) for Γ k . First we derive the flow equation for the effective potential, i.e., the part of the effective action that involves all momentum-independent couplings. Then we proceed to specify the flow of the inverse propagator which determines flow equations for the wave-function renormalization Z and the gradient coefficientK. In the FRG, we approach the critical point from the ordered (symmetry-broken) side of the transition. This allows us to capture the leading divergences of two-loop effects in a calculation that is formally one-loop 70 by means of diagrams like the second one in Eq. (32) in the spirit of the background field method in gauge theories. 102 We denote the truncation Eq. (43), evaluated for homogeneous, i.e., space-and time-independent "background fields" by (the subscript cq indicates that we have classical and quantum background fields) where Ω is the quantization volume and the U(1) invariant combinations of fields areρ cq =φ * cφq =ρ * qc and ρ q = φ q 2 . This representation of Γ k,cq implies that the flow equation for the potentialŪ can be obtained from Eq. (42) by taking the derivative with respect toρ cq and setting the quantum background fields to their stationary value (which is zero) afterwards, where the dimensionless RG flow parameter t is related to the cutoff scale k via t = ln(k/Λ). The flow equation for the renormalized potential follows straightforwardly by taking the scale derivative of the relationŪ = ZU, which results in where we introduced the anomalous dimension of the wavefunction renormalization, Then, using ∂ρ cq = Z∂ ρ cq , the flow equation for the renormalized potential can be written as We proceed by specifying the projection prescriptions that allow us to derive the flow of the couplings u n in the ordered phase from the flow equation (81). Taking the scale derivatives of the relation u n = U (n) (ρ 0 ) we find Based on the power-counting arguments of Sec. IV C, our truncation includes terms up to cubic order in the U(1) invariants, i.e., for derivatives of the effective potential of the order of n ≥ 4 we have U (n) = 0. The flow equations for the quartic and sextic couplings are then given by (the RHS of these equations determine the so-called β-functions) where according to Eq. (82) in ζ we specify the classical background field ρ c it to its stationary value ρ c | ss = ρ 0 . As we have seen above (cf. Secs. III A and IV A, the latter is determined by the dissipative part of the field equation, i.e., by the condition Im U (ρ 0 ) = 0. Taking here the derivative with respect to the RG parameter t, we find Having thus specified the flow equations for the couplings that parameterize the potential U, we proceed to the Keldysh mass γ, which is the coefficient of the term that is proportional to the quantum U(1) invariantρ q in Eq. (77). We obtain the flow equation forγ as For the renormalized Keldysh mass, which is related to the bare one via γ =γ/ |Z| 2 , we have (the transformation from bare to renormalized fields implies ∂ρ q = |Z| 2 ∂ ρ q ) While the flow of Γ k,cq (i.e., the flow equation evaluated at homogeneous background fields) yields flow equations for all momentum-independent couplings, we have to consider the flow of the inverse propagator in order to derive flow equations for the wave-function renormalization Z and the gradient coefficientK. The retarded component of the inverse propagator in the presence of real stationary background fieldsφ c =φ * c =φ 0 reads Then, for the kinetic coefficientK we choose from the flow equation (88) the elements of the inverse propagator that do not have mass-like contributions 70 2λρ 0 and 2κρ 0 , The flow equation for the wave-function renormalization Z as specified below, on the other hand, mixes massive and massless components symmetrically This choice allows for the locking of the flows of the Keldysh mass and Z as implied by the emergence of the symmetry Eq. (60) in the purely dissipative IR regime (see Sec. VIII). Finally, the flow equation for the renormalized coefficient K follows by straightforward differentiation of its definition K =K/Z in terms of bare quantities. We find The truncation Eq. (43) is parameterized in terms of the couplings Eq. (66). Renormalization of the fields with Z leads to a description in terms of g = (K, ρ 0 , u 2 , u 3 , γ) T (where we omit the mass u 1 : as indicated above we approach the critical point from the ordered phase, i.e., we parameterize the effective action in terms of the stationary condensate density ρ 0 instead of the mass u 1 ). In this section we derived the βfunctions for these renormalized couplings, i.e., we have specified a closed set of flow equations ∂ t g = β g (g) from which Z can be completely eliminated (the anomalous dimension η Z entering the β-functions can again be expressed in terms of the couplings g alone). More explicit expressions for the βfunctions are provided in App. C 3).

VIII. SCALING SOLUTIONS
As one considers an effective description of a system at a continuous phase transition at longer and longer scales (which is equivalent to following the RG flow to smaller values of k), physical observables and the couplings that describe the system exhibit scaling behavior. The search for such scaling solutions to the flow equations is facilitated by introducing rescaled dimensionless (in the sense of the canonical power counting introduced in Sec. IV C) couplings which remain constant, i.e., by searching for a fixed point of the flow equations of these rescaled couplings instead. In the following section we introduce such rescaled couplings and derive the corresponding flow equations.

A. Scaling form of the flow equations
As a first step we trade the real parts of K, u 2 , and u 3 for the ratios of real to imaginary parts which measure the relative strength of coherent and dissipative dynamics. As we will show below, at criticality all these ratios flow to zero signaling decoherence. Their flow is given by (The β-functions for the real and imaginary parts of K, u 2 , and u 3 are specified in App. C 2, see Eq. (C34).) We proceed by introducing a dimensionless mass term the flow equation of which mixes contributions from the βfunctions of ρ 0 , κ, and D, and reads where the anomalous dimension of D is defined as Finally we replace the quartic and sextic couplings by dimensionless ones. For a momentum-independent n-body coupling u n we can construct a corresponding dimensionless coupling by means of the relatioñ The flow equations for the imaginary partsκ andκ 3 of the dimensionless quartic and sextic couplings, therefore, are given by and include contributions from the anomalous dimension Thus we are left with six dimensionless running couplings, which we collect in vectors r = r K , r u 2 , r u 3 T and s = (w,κ,κ 3 ) T . Their flow equations form a closed set, The β-functions on the RHS of these equations contain the anomalous dimensions η Z , η D , and η γ , which in turn can be expressed as functions of the running couplings r and s alone. We note in passing that according to the discussion of Sec. VI D, the equilibrium model MAR is described by r K = r u 2 = r u 3 = r, i.e., (MA is realized for the special case r = 0). Inserting the same value r for all three ratios in the respective β-functions we find β r K = β r u 2 = β r u 3 , which shows that for MAR the common ratio is preserved by the flow as it should be. Our analysis of the flow equations (104) will proceed in two steps: First we will search for fixed points r * and s * , which are solutions to the algebraic equations β r (r * , s * ) = β s (r * , s * ) = 0.
In Sec. VIII B we briefly discuss the trivial Gaussian fixed point and then turn to the Wilson-Fisher fixed point that describes the critical system in VIII C). Second we will solve the full flow equations numerically and provide our results in Sec. IX. While already the linearized flow equations in the vicinity of the Wilson-Fisher fixed point encode universal physics at the phase transition and determine the asymptotic flow of the system for k → 0 (or t → −∞), the numerical integration of the full flow equations provides us with information on non-universal aspects such as the extent of the scaling regime.

B. Gaussian fixed point
All β-functions vanish on the manifold of Gaussian fixed points which is parameterized by s * = 0 and r * ∈ R 3 . We note that the combination of vanishing imaginary partsκ * andκ 3 * of the quartic and sextic couplings and arbitrary finite ratios of real to imaginary parts implies that also the real parts ofũ 2 * andũ 3 * are zero on this fixed point manifold. In a linearization of the flow equations around s * = 0, the fluctuation contributions vanish and the scaling behavior is determined solely by the canonical scaling dimensions, implying in particular that the Gaussian fixed point is unstable (for small values s s * the flow is directed away from the fixed point) and, therefore, physically not relevant. Non-trivial scaling behavior at criticality is governed by the Wilson-Fisher fixed point which we will discuss in the next section.

C. Wilson-Fisher fixed point: critical behavior
As discussed in Sec. VI, our driven-dissipative model reduces to MA when we set the real parts of all renormalized couplings to zero, cf. Fig. 5, i.e., for r = 0. It is well-known that MA exhibits a non-trivial Wilson-Fisher fixed point, 32 and indeed we find this fixed point at r * = (r K * , r u * , r u * ) = 0, s * = w * ,κ * ,κ * = (0.475, 5.308, 51.383) .
The values of the coupling s * are identical to those obtained in an equilibrium classical O(2) model from functional RG calculations at the same level of truncation. 70 We note that this fixed point is also contained in the subspace of couplings corresponding to MAR, which is characterized by Eq. (105), i.e., the phase transitions in both the equilibrium and nonequilibrium models are described by the same fixed point.  Table II. Results for the correlation length exponent ν, the anomalous dimension η, the dynamical critical exponent z, and the decoherence exponent η r in our truncation.
Critical behavior, however, is determined by the RG flow in the vicinity of the fixed point. Here the non-equilibrium setting adds two more independent directions, thereby opening up the possibility for deviations from equilibrium criticality as we will now show.
The asymptotic flow for k → 0 of the critical system is determined by a linearization of the flow equations in the deviations δs = s − s * , δr = r from the fixed point. In the linear regime the two sectors corresponding to s and r decouple as described by the block diagonal stability matrix where the 3 × 3 submatrices S and N are given by The matrix N would be identically zero in the absence of anomalous additions to the canonical scaling dimensions (note that the ratios r have canonical scaling dimension zero), or even if coherent and dissipative couplings would exhibit identical anomalous scaling. The non-vanishing of this block thus indicates a different universal behavior of these two types of couplings. Due to the decoupling of the flows of r and s we may discuss the linearized flow of each set of couplings separately.
In the matrix S we find one negative eigenvalue s 1 corresponding to the correlation length exponent ν = −1/s 1 = 0.716 (our findings for critical exponents are summarized in Tab. II). Considering that we are restricting ourselves to relevant and marginal terms in our truncation, the agreement of the numerical value of ν with results from more sophisticated calculations 103 is reasonable. Furthermore there are two complex conjugate eigenvalues s 2,3 = 1.124 ± i0.622 with positive real parts (indicating that these directions are stable). The imaginary parts are known artifacts of this level of truncation for the O(2) model and vanish upon inclusion of higher order terms in the effective potential. 104 The scaling behavior of the couplings Z, D, and γ is determined by the values of the respective anomalous dimensions at the fixed point. In addition we define the anomalous dimension for the bare kinetic coefficientK as (111) where the representation in terms ofη A andη D follows from the definition of these quantities in Eq. (C32). At the fixed point η takes the value which is again the result for the anomalous dimension of the classical O(2) model in d = 3 dimensions at the same level of truncation 70 and agrees well with results from more accurate calculations. 103 In summary, the static critical behavior coincides precisely with the one of the classical O(2) model, implying that the dynamical anomalous dimension η Z effectively does not enter the corresponding equations. This can be seen as follows: Inserting r = 0 in the expressions for the anomalous dimensions, we find (We note that this holds for all values of the static couplings s, i.e., it is always realized in MA.) These relations ensure that η ZR and η γ compensate each other in all flow equations. 105 Moreover they imply that the ratioγ/ |Z| appearing on the RHS of the fluctuation-dissipation relation Eq. (73) approaches a constant value at the fixed point: According to the definition of the anomalous dimensions Eqs. (80) and (103), close to the fixed point the flow of Z and γ is described by Z ∼ k −η Z (note that η Z is real so that this behavior indeed describes algebraic scaling and does not contain oscillatory parts) and γ ∼ k −η γ with η Z and η γ evaluated at r * and s * . Thus we findγ/ |Z| = |Z| γ ∼ k −η Z −η γ = const., i.e., the symmetry Eq. (60), which manifests itself in this quantity approaching a constant value (cf. Eq. (74)), emerges in the IR without imposing it in the microscopic model. In other words, the driven-dissipative system obeys a classical FDT in the long-wavelength limit (see Fig. 2). At the fixed point we find the value Let us now consider the upper left block N of the stability matrix. It has three positive eigenvalues, which indicates that the ratios r are attracted to their fixed point value zero. The corresponding eigenvectors are The smallest of the eigenvalues determines the scaling behavior of r in the deep IR. In order to see this let us expand r in the basis of eigenvectors of the matrix N, The coefficients in this expansion are given by c i = v i ·r, where v i are the left eigenvectors of N (the latter is not symmetric and its left and right eigenvectors, therefore, are not equal), normalized such that u i · v j = δ i j . The asymptotic behavior of the flow of the so-called scaling fields 78 c i is given by c i ∼ e n i t = k n i , which implies that for r we indeed find with only subdominant contributions in the directions of u 2 and u 3 . This leads us to identify the decoherence exponent From the scaling behavior of the ratios r we may infer the one of the coherent couplings. For the coefficient of coherent propagation A, in particular, we have Then, with the anomalous dimension of the dissipative kinetic coefficient D at the fixed point, we obtain the value Let us discuss the consequences of this result for the effective dispersion relation of long-wavelength excitations, which is encoded in the running inverse propagator Eq. (49). Once the cutoff scale k becomes smaller than the external momentum q, the effective infrared cutoff is given by q instead of k. 106 Then, in the dispersion relation Eq. (54), which we rewrite here in terms of the scaling variables introduced in Sec. VIII A as ω R 1,2 = Dq 2 −i (1 + w/2) ± r 2 K + r K r u 2 w − (w/2) 2 , (123) we may insert the scaling forms w ∼ w * , r K ∼ r K0 q −η r , r u 2 ∼ r u 2 0 q −η r , and D ∼ D 0 q −η D . For q → 0 both modes are purely diffusive with ω R 1 ∼ −iD 0 q 2−η D and ω R 2 ∼ −iD 0 q 2−η D (1 + w * ), and for the dynamical critical exponent z which is defined via the relation ω R ∼ −iq z we find the value Above the purely diffusive IR regime, when w r K , r u 2 , the dispersion relation simplifies to i.e., coherent propagation and diffusive contributions scale differently with the momentum q. In the symmetric phase the branch ω R 2 is absent and the bare retarded response function is dominated by a single pole at ω = ω R 1 , i.e., we have (the bare scale-dependent propagators are determined by the quadratic part of the effective action Eq. (43)) As explained at the end of Sec. II this quantity and, in particular, the spectral density which is related to its imaginary part, 78 are direct experimental observables. For ω ≈ ω R 1 the spectral density has the shape of a Lorentzian centered at Re ω R 1 and with width determined by Im ω R 1 , Inserting here the scaling forms Z ∼ Z 0 k −η Z and Eq. (125) for ω R 1 with different scaling behavior of real and imaginary parts, the structure sketched in Fig. 3 emerges. For the specific values of the anomalous dimensions obtained in this section the spectral density is shown in Fig. 6, where a pronounced feature is clearly visible.
Before moving on to a numerical integration of the flow equations in the next section, we briefly contrast our findings for the DDM with the equilibrium case of MAR. There, analyzing the stability of the fixed point Eq. (107) we have to take into account only one direction r = r K = r u 2 = r u 3 , and we find (as r * = 0 we have δr = r) where the matrix S is the same as above and the element R is given by the "middle" eigenvalue Eq. (115) of the stability matrix N in the non-equilibrium problem, i.e., also in the equilibrium setting we find decoherence at the longest scales, however, with a value of the decoherence exponent that is different from the one in non-equilibrium. Let us finally remark that in the linearized regime, the fact that MAR is contained as a special case in the non-equilibrium problem, becomes visible in the form of the second eigenvector Eq. (116) which realizes the constraint Eq. (105). We finally comment on the relation of the critical exponents obtained here with other approaches. The static critical exponent η shows very good agreement with sophisticated high order perturbative calculations. For the dynamical critical exponent z, to the best of our knowledge, currently there are no high-precision results for MA with N = 2 components available. (The situation is different for the Ising-symmetric case with N = 1, where the dynamical critical exponent has been calculated with high accuracy, see Ref. 107 and references therein.) Thus the value z = 2.121 obtained here has to be compared to z = 2.026 which corresponds to the third order in -expansion 81,108 . The decoherence exponent (η r = −0.101 here) has been computed in a recent complementary perturbative field theoretical study to second order in -expansion, where it takes the value η r = −0.003, see Ref. 57. The discrepancy between these values can only be resolved by extending the truncation advocated here, by including higher order corrections in pertrubative field theoretical approaches, or by means of large-scale computer simulations.

IX. NUMERICAL INTEGRATION OF FLOW EQUATIONS
In the previous section we have seen that the flow equations Eq. (104) entail non-trivial critical behavior governed by the Wilson-Fisher fixed point Eq. (107). While these results were based on an analysis of the linearized flow equations in the vicinity of the fixed point, we will now turn to a numerical integration of the full non-linear equations. One the one hand, this serves to illustrate the concept of universality: Independently from the initial values r Λ ,κ Λ , andκ 3Λ at the mesoscopic starting point of the RG flow, critical behavior can be induced by a proper fine-tuning of w Λ and becomes apparent in the approach of the RG flow to the scaling solution. Apart from that, the availability of the full flow in the framework of the FRG allows us to extract non-universal aspects. In particular, we will give an estimate of the Ginzburg scale, i.e., the scale that separates the region of non-universal flow from the universal scaling regime and thus is important for determining experimental requirements on the necessary frequency resolution.
Our approach for finding numerical solutions to the flow equations that exhibit critical behavior is as follows: We choose initial values r Λ ,κ Λ , andκ 3Λ at the mesoscopic scale k = Λ (t = 0), which are appropriate for the description of the model introduced in Sec. III. This model contains two-body elastic interactions and loss, while three-body terms are contained only in an effective low-momentum description, implyingκ Λ ≈ 1 andκ 3Λ 1. The dissipative kinetic coefficient D is very small in the microscopic description, so that r KΛ 1 initially, while for the two-body terms we have r u 2 Λ ≈ 1. The latter generate the three-body couplings and we assume that r u 3 Λ ≈ 1 as well. For such a choice of initial values, there is a critical value w Λ = w c so that the resulting RG trajectories r(t) and s(t) approach the scaling solution, i.e., the fixed point, for k → 0 (t → −∞). Any solution obtained by numerically integrating the flow equations with w Λ fine-tuned to w c , however, eventually always flows away from the fixed point, as due to limited accuracy the solution develops a non-zero component in the unstable direction of the fixed point at some stage. For all solutions shown in the figures we choose w Λ slightly below w c , so that the trajectory at large RG "times" t flows to the symmetric phase with w = 0.
When such a near-critical trajectory approaches the scaling solution, the couplings s flow towards their fixed point values s * on a scale 1/ Re s 2,3 ≈ 1 determined by the eigenvalues s 2,3 of the stability matrix S , cf. Fig. 7, and stay there for a long "time" t s . Depending on how close w Λ is to w c , this duration is typically t s = 10 to 20 which corresponds to several orders of magnitude in k/Λ. During t s the ratios r decay according to Eq. (117), i.e., as the sum of three exponentials, with decay rates given by the eigenvalues Eq. (115) of the stability matrix N. In order to extract these eigenvalues from the numerical solution, we consider the flow of the coefficients c i ∼ e n i t in the expansion of r in the basis of eigenvectors of N Eq. (117). Figure 7 shows c 1,2 along with exponential fits, which reproduce the eigenvalues n 1,2 to satisfactory accuracy.
An important result of the previous section is the scaling relation Eq. (114) between the anomalous dimensions η Z and η γ of the wave-function renormalization and the Keldysh mass respectively, which implies that when evaluated along a critical trajectory, the value of −η γ approaches the one of the real part η ZR of η Z , while the imaginary part η ZI flows to zero. This prediction -physically implying asymptotic thermalizationis verified numerically in Fig. 8.
As the anomalous dimensions η a of a = Z, D, and γ are functions of the renormalized dimensionless couplings r and s alone and not the quantities a themselves, we get the solutions to the flow equations ∂ t a = −η a a simply by exponentiating the integrals of the anomalous dimensions along RG trajectories r(t) and s(t), i.e., In this way we obtain the trajectories of K shown in Fig. 9 and the flow of the effective temperature T eff =γ/ (4 |Z|) = γ |Z| /4 which according to the discussion in Sec. VI C at low frequencies saturates to a constant value as illustrated in Fig. 2. While this asymptotic value depends on the initial values ofγ and Z at the scale Λ and is, therefore, non-universal, the manner in which it is approached is universal as it is determined by the exponent η r : According to Eq. (131) the flow of T eff is given by Close to the fixed point we may expand the anomalous dimensions in the exponential in powers of δr = r and δs. As both η ZR and η γ are even functions of r there is no linear term in the expansion and we may write for a = ZR and γ (here we are indicating the anomalous dimension evaluated at the fixed point explicitly as η a * ): |δs| ∼ e Re s 2,3 t , where s 2,3 are eigenvalues with positive real parts Re s 2 = Re s 3 of the stability matrix S Eq. (109) and determine the leading corrections to scaling in the static sector. Therefore, neglecting exponentially small corrections, we have Note that the quantities η a depend on the precise prefactor in the scaling form Eq. (118) of r, i.e., they depend on microscopic parameters and are thus non-universal. Then, using Eq. (134) and keeping in mind that η γ * = −η ZR * , we find the asymptotic behavior where in the last line we are again dropping exponentially small corrections. This form confirms the physical intuition that long-wavelength thermalization of the DDM is governed by the exponent that is unique to this model and manifestly witnesses the microscopic non-equilibrium nature of this model. We finally note that the effective temperature defined in Eq. (74) is not the one that enters the FDT Eq. (72) for MAR. The latter can be established by means of the basis transformation Eq. (59) which involves the parameterr. This parameter, however, is characteristic the of MAR and has no counterpart in the DDM. The near-critical trajectories we consider in this section illustrate the concept of universality in that they show how details of the microscopic model, which determine the initial conditions of the RG flow, are lost as we lower k → 0, where all of these trajectories converge towards the scaling solution, cf. Fig. 1. However, a distinctly non-universal feature of these trajectories is the point where the crossover to the universal regime takes place, which is known as the Ginzburg scale. [79][80][81]  Physically, the Ginzburg scale marks the breakdown of meanfield theory as we approach the fluctuation-dominated critical region. In a perturbative estimate in the symmetric phase, we compare the bare distance from the phase transition κ 1 to the corresponding one-loop correction. Demanding these quantities to be of the same order of magnitude yields 19 where C is a numerical constant (we find C = 2π if we set the bare value κ 1G exactly equal to its one-loop correction). Expressing κ 1G through a momentum scale as κ 1G = D Λ k 2 G we find Eq. (3), and for the dimensionless RG "time" t G = ln (k G /Λ), in terms of the dimensionless two-body loss rateκ introduced in Sec. VIII A, we have Fitting this logarithmic dependence to numerically obtained trajectories in Fig. 1, we find C ≈ 14.8. The Ginzburg scale delimits also the region where the driven-dissipative system obeys a FDT and the ratio T eff =γ/ (4 |Z|) saturates to a constant value as shown in Fig. 2.

X. CONCLUSIONS
We have studied the nature of Bose criticality in driven open systems. To this end, starting from a description of the microscopic physics in terms of a many-body quantum master equation, we have developed and put into practice a FRG approach based on a Keldysh functional integral reformulation of the quantum master equation for the quantitative determination of the universality class. The absence of both an exact particle number conservation and the detailed balance condition were seen to underly the existence of a new and independent critical exponent governing universal decoherence, while the distribution function shows asymptotic thermalization despite the microscopic driven nature of the system. This work is just a first step in the exploration of nonequilibrium critical behavior. Key questions for future studies concern the status of critical points in lower dimensionality as, e.g., relevant for current exciton-polariton systems. In particular, in Ref. 66 it has been shown that the thermal fixed point is unstable in two dimensions, and instead is replaced by the non-equilibrium Kardar-Parisi-Zhang 43 fixed point. It is also a key issue to investigate different symmetries beyond the O(2) case. For example, Heisenberg models realized with ensembles of trapped ions may exhibit O(3) symmetry. 109 Furthermore, given the fact that many light-matter systems are pumped coherently as opposed to the incoherent pump considered here, it will be important to understand the impact of the coherent drive on potential criticality in these classes of systems. Finally, it is an intriguing question whether driven open systems which realize non-equilibrium counterparts of quantum criticality can be identified. In the long run, it remains to be seen whether a classification of non-equilibrium criticality with similarly clear structure as familiar from equilibrium dynamical criticality 32 can be reached. action. We start from a master equation governing the time evolution of a system density matrix, Here,Ĥ s is a system Hamiltonian generating the unitary evolution andL is a Lindblad operator making up the dissipative part of the Liouvillian. For simplicity we consider only a single dissipative channel. The generalization to several channels as in Eq. (6), realized through the coupling to several baths, is straightforward. Equation (A1) results from a more general system-bath setting,Ĥ tot =Ĥ s +Ĥ b +Ĥ sb (Ĥ b and H sb are a quadratic bath Hamiltonian with a continuum of frequencies and a system-bath Hamiltonian linear in the bath operators, respectively) under the following three assumptions: (i) the system-bath coupling γ(ω) is weak compared to a typical scale ω 0 in the system (or, by energy conservation, in the bath) indicating, e.g., the level spacing in an atom (Born approximation γ(ω)/ω 0 1) (ii) the frequency dependence of the system-bath coupling is negligible over the bandwidth ϑ of the bath centered around ω 0 , implying δ-correlations in the time domain (Markov approximation γ(ω) ≈ const.), and (iii) the system is driven with an external field with frequency ν to bridge the large energy separation of the levels, (ν − ω 0 ) / (ν + ω 0 ) 1. This makes it possible to work in the rotating wave approximation, in which only the detuning ∆ = ν − ω 0 occurs as a physical scale, while all fast terms involving ν + ω 0 are dropped. From this consideration, it is clear that the master equation is an accurate description of strongly driven systems coupled to an environment. A typical realization in quantum optics is an atom with two relevant levels separated by ω 0 , connected by an external laser drive with frequency ν, which is detuned from resonance by ∆ = ν − ω 0 . Only the laser drive makes the excited level accessible and gives rise to two-level dynamics such as Rabi oscillations, with frequency determined by the laser intensity. The excited level is unstable and can undergo spontaneous emission by coupling to the radiation field, providing for the reservoir -this mechanism is physically completely independent of the coherent dynamics. Alternatively but fully equivalent to the operator formalism, the above approximations can be performed in a Keldysh path integral setting (see below). In this way, the physics of a given quantum master equation becomes amenable to quantum field theoretical approaches, which is particularly useful for bosonic and fermionic drivendissipative many-body systems. Here the starting point is the Keldysh partition function which results from a "Trotterization" of the Hamiltonian dynamics (after normal ordering) acting on the density matrix in the integrated form of the von Neumann equation in the basis of coherent states; in this process, the second quantized system and bath field operators,â i (the index i denotes both position and internal indices, such as different particle species) andb µ (µ labels the bath modes and will be chosen a continuous index below) respectively, are mapped to time-dependent complex valued fields in the action where H tot,σ (t) is a quasilocal polynomial of these fields. The relative minus sign for the evolution on the forward (+) and backward (-) contours clearly reflects the commutator structure in the von Neumann equation of motion for the systembath density operator above. We have omitted an imaginary regularization term ensuring convergence of the functional integral [76][77][78] for simplicity, as it does not affect any of the next steps. Integrating out the harmonic bath variables using approximations (i) -(iii) and considering for the moment Lindblad operatorsL which are linear in the system field operators, we arrive at the following effective Markovian dissipative action: While the relative minus sign for the system Hamiltonian H s on the + and -contours preserve the commutator structure, the dissipative terms clearly reflect the temporally local Lindblad structure of Eq. (A1). We thus arrive at a simple translation rule for bosonic 110 master equations into the corresponding Keldysh functional integral: (i) the temporal derivative terms can be read off from the last equation; (ii) for all (normal ordered) operators on the right (left) of the density matrix, introduce a contour index + (-) and write down the Markovian dissipative action. The linear Lindblad operators we consider here are not affected by normal ordering. For the more general case of Lindblad operators that are quasilocal polynomials in the system field operators, operator ordering can be tracked by a suitable temporal regularization procedure as elaborated in the next section.

Derivation in the Keldysh setting
Here we present a derivation of the Markovian dissipative action in the ± basis for arbitrary (non-linear) Lindblad jump operators, which allows for the most direct comparison with the master equation. In particular, we pay special attention to the question how the operator ordering in the master equation is reflected in the path integral formulation. We leave the system action unspecified, requiring only the property that after proper rotating frame transformation the evolution of the system is much slower than the correlation time of the bath τ c = 1/ϑ (broadband bath). The action of the bath is, in the ± basis, The Green's functions for the oscillators of the bath are assumed to be in thermal equilibrium and read The linear coupling between system and the bath is (note that the case of several dissipative channels and local baths as in Eq. (6) can be implemented by adding appropriate indices to the L σ and b µ,σ and summing over these indices) where L ± correspond to the quantum jump operators which are typically quasilocal polynomials of the system's creation and annihilation operators. To be consistent with the derivation of the path integral, we require the jump operators to have been normal ordered before the Trotter decomposition giving rise to the path integral. The partition function is of the general form Now we integrate out the bath via completion of the square which results in an effective action S eff for the system degrees of freedom. The contribution S eff,µ of the µth mode to the effective action reads The signs for the operators on the − contour comes from the backward integration in time. Thus the mixed terms will occur with an overall − sign, while the ++ and −− terms come with an overall +. Summing over all the modes µ we obtain the effective action for the field variables of the subsystem due to the coupling to the bath. We now take the continuum limit of densely lying bath modes, centered around some central frequency ω 0 and with bandwidth ϑ. That is, we substitute the sum over the modes with an integral in the energy Ω weighted by a (phenomenologically introduced) density of states ν(Ω) of the bath µ γ µ ∞ 0 dΩγ(Ω)ν(Ω), and obtain where in addition we have used the translation invariance of the bath Green's function, G αβ Ω (t, t ) = G αβ Ω (t − t ) to suitably shift the integration variables. We consider the various terms separately. In doing the Markov approximation, we use (a) that by assumption it is possible to choose a rotating frame in which the evolution of the system is slow compared to the scales in the bath, ω sys ω 0 , ϑ. In this case, a zeroth order temporal derivative approximation for the jump operators is appropriate. This gives rise to a temporally local form of the markovian dissipative action. However, for the evaluation of tadpole diagrams for this action, ambiguities due to a temporally local vertex arises. In these diagrams -and only in these -it is then important to specify the proper regularization of the system's Green's function at equal time arguments. To keep track of this, we indicate the sign of the next time step in the approximated jump operators by t ±δ = t ± δt. In step (b) below, we assume that the density of states and the coupling of the system to bath are well approximated as constant over the relevant reservoir width, where we have shifted the frequency integration domain by −ω 0 and taken the limit ϑ → ∞, as well as κ = γν andn =n(ω 0 ). Further note the relation to the operator formalism . For the terms on the forward contour, we obtain In the last line we have used and dτθ(−τ) Importantly, note the sign change in the regularization of the time argument upon reversal of integration direction. This gives a hint which operator "comes first" in the coarse grained evolution where the bath has been integrated out, and reflects the fact that in the corresponding master equation, the "cooling" dissipation terms ∼ (n + 1) are normal ordered in the jump operators (∼L †L ), while the "heating" terms ∼n are anti-normal ordered (∼LL † ). Similarly, we obtain on the backward contour, where the changes in the signs relative to the forward term emerge from the reverse signs in the θ-functions. In summary, we obtain the following dissipative contribution to the action: In addition, there is a "Lamb shift" which reads This gives a contribution to the coherent dynamics which has the same physical origin as the dissipative dynamics. However, typically there is a dominant independent Hamiltonian contribution, such that the effective Hamiltonian parameters after the Lamb shift renormalization are properly regarded as independent of the Liouvillian ones. by D(Q) = P(Q) − P(0), and the effective potential V that contains all momentum-independent couplings is given by The second functional derivative of the effective action can then be expressed as the sum of two contributions, where the second term is just the functional derivative of the effective potential, which can be reduced to ordinary (i.e., not functional) partial derivatives with respect to the fields in the time domain and real space, Setting the fluctuation components of the fields to zero in Eq. (C9) we obtain the inverse propagator in the presence of homogeneous classical and quantum background fields, (C12) Note that the difference between P cq (Q) and the inverse propagator Eq. (49) is that in the latter the background fields are set to their stationary values while in the former they remain unspecified. The background fields are all contained in the second contribution V (2) cq which we split into 2 × 2 blocks according to While the upper left block V (2)H cq is linear in the quantum fields (and, therefore, vanishes when we set these to zero, giving rise to the causality structure of the inverse propagator Eq. (49)), V (2)H cq,11 = ρ cq + ρ qc U (3) H + i ρ cq − ρ qc U (3) D χ 2 c,1 + ρ cq + ρ qc + 2χ c,1 χ q,1 U H + i ρ cq − ρ qc + i2χ c,1 χ q,2 U D , V (2)H cq,12 = V (2)H cq,21 = χ c,2 χ q,1 + χ c,1 χ q,2 U H + χ c,1 χ q,1 − χ c,2 χ q,2 U D + χ c,1 χ c,2 ρ cq + ρ qc U (3) H + i ρ cq − ρ qc U (3) D , V (2)H cq,22 = ρ cq + ρ qc U (3) H + i ρ cq − ρ qc U (3) D χ 2 c,2 + ρ cq + ρ qc + 2χ c,2 χ q,2 U H + i ρ cq − ρ qc − i2χ c,2 χ q,1 U D , where the matrix F is given by the sum F = F 1 +F 2 with F 1,2 being of first and second order in δχ. The explicit dependence of these matrices on the fluctuations reads Inserting the decomposition Eq. (C17) in Eq. (C6) and expanding the logarithm in the fluctuations δχ yields ∂ t Γ k = i 2 Tr∂ t ln P k,cq − 1 2∂ t Tr G k,cq F 1 2 , where G k,cq (Q) = P k,cq (Q) −1 is the propagator in the presence of classical and quantum background fields. Note that the appearance of G 2 k,cq makes the trace in the last term UVconvergent and thereby allowed us to commute∂ t with Tr. In the expansion Eq. (C21) we are keeping only terms of zeroth and second order, as these determine, respectively, the flow of the effective potential and the frequency-and momentum-dependent contributions to the inverse propagator. We also omit a term∂ t Tr G k,cq F 2 which in our truncation with momentum-independent vertices does not contribute to the flow of Z andK. where det cq (ω, q 2 ) = det P k,cq (Q) denotes the determinant of the regularized inverse propagator Eq. (C16) in the presence of classical and quantum background fields. Since our model is symmetric under simultaneous phase rotations φ ν → e iα φ ν of the classical and quantum fields, the determinant det cq (ω, q 2 ) can be expressed as a function of the U(1)invariant field combinations ρ c , ρ cq , ρ qc , and ρ q . It can not be written as a function of these invariants without ambiguity though, as can be seen by noting that the product of four fields φ * c φ * q φ c φ q equals both ρ c ρ q and ρ cq ρ qc . However, the form of the field-dependent contribution Eq. (C13) to the inverse propagator implies that det cq (ω, q 2 ) contains terms that are at most quadratic in the quantum fields and that there is no contribution that contains φ * q φ q but no classical fields. All contributions containing quantum and classical fields can be expressed in powers of ρ c , ρ cq , and ρ qc . Therefore, in the following we will consider det cq (ω, q 2 ) to be a function of this reduced set of invariants. Then, inserting Eq. (C22) in the definition of ζ in Eq. (81) we find where det c (ω, q 2 ) = det P k,c (Q) is the determinant of the regularized propagator with only classical background fields, Recalling the definition Eq. (C4) of the differential operator∂ t according to which it effectively acts as a scale derivative of the bare cutoff, we find ∂ t R k,A (x) = Re ∂ t R k,K (x)/Z , ∂ t R k,D (x) = Im ∂ t R k,K (x)/Z .
(C33) The further evaluation of this expression is most conveniently performed on the computer using Mathematica.

Flow equation for the inverse propagator
The second term on the RHS of Eq. (C21) determines the flow of both the wave-function renormalization and the gradient coefficient. It is quadratic in the fluctuations δχ, hence we can write it as where we set the fields to their stationary values. Σ(Q) can be visualized as consisting of one-loop diagrams with four external legs two of which are attached to the condensate (cf. the second diagram on the RHS of Eq. (32)) and is given by where G k (Q) = P k (Q) −1 with the inverse propagator given by Eqs. (51) and (52) to which the cutoff R k (q 2 ) has to be added. For φ c = φ * c = φ 0 and φ q = φ * q = 0 the matrices V (3) i have the structure (C41) Inserting this expression in Eq. (C40) above and taking the causality structure of the propagator into account, we can rewrite the integrand in the form (P + = P + Q) Then the second and third equalities in Eq. (C41) imply that Σ(Q) has the same causality structure as the inverse propagator. For the retarded block we find where now the indices i and j take the values 1, 2, and the Keldysh component is given by The frequency integrals appearing in Eqs. (C43) and (C44) can be evaluated by straightforward application of the resiude theorem: G R k (Q) has simple poles ω R 1,2 given by Eq. (54) with Aq 2 and Dq 2 replaced by p A (q 2 ) and p D (q 2 ) respectively. While the poles of the advanced propagator ω A 1,2 are complex conjugate to the poles of the retarded propagator, G K k (Q) has poles at both ω R 1,2 and ω A 1,2 . We omit the lengthy expression for Σ(Q) after frequency integration.
Combining Eqs. (88) and (C21), the flow equation for frequency-and momentum-dependent part of the the bare inverse propagator can be written as with the matrix z defined in Eq. (C2). Inserting this expression in the flow equations for the wave-function renormalization Z and the gradient coefficientK, Eqs. (91) and (90) respectively, we find after some algebra, The real and imaginary parts of the anomalous dimension η Z , which appear in the flow equations (C34) of the real quartic and sextic couplings, are then given by Here we used the relation Σ R (Q) = Σ R (−Q) * which implies ∂ ω Σ R (0) = −∂ ω Σ R (0) * . To further evaluate η ZR and η ZI we switch to Mathematica. The derivatives with respect to the frequency can be carried out without any difficulty and∂ t can be calculated as in Eq. (C28) above. Again the integral over spatial momenta is facilitated by the θ-function contained iñ ∂ t R k,a (x) and can be carried out analytically. Finally, for the real and imaginary parts of the renormalized kinetic coefficient K =K/Z = A + iD we have where using ∂ q 2 Σ R (0) = ∂ q 2 Σ R (0) * (note that Σ(Q) depends only on the norm squared q 2 of the spatial momentum) we may express the quantitiesη A andη D defined in Eq. (C32) as We will proceed with the evaluation of these expressions in the next section.