Perturbative Field-Theoretical Renormalization Group Approach to Driven-Dissipative Bose-Einstein Criticality

The universal critical behavior of the driven-dissipative non-equilibrium Bose-Einstein condensation transition is investigated employing the field-theoretical renormalization group method. Such criticality may be realized in broad ranges of driven open systems on the interface of quantum optics and many-body physics, from exciton-polariton condensates to cold atomic gases. The starting point is a noisy and dissipative Gross-Pitaevski equation corresponding to a complex valued Landau-Ginzburg functional, which captures the near critical non-equilibrium dynamics, and generalizes Model A for classical relaxational dynamics with non-conserved order parameter. We confirm and further develop the physical picture previously established by means of a functional renormalization group study of this system. Complementing this earlier numerical analysis, we analytically compute the static and dynamical critical exponents at the condensation transition to lowest non-trivial order in the dimensional epsilon expansion about the upper critical dimension d_c = 4, and establish the emergence of a novel universal scaling exponent associated with the non-equilibrium drive. We also discuss the corresponding situation for a conserved order parameter field, i.e., (sub-)diffusive Model B with complex coefficients.


I. INTRODUCTION
Experimental systems that are characterized by a strong coupling of light to a large number of matter degrees of freedom 1 hold the potential of developing into laboratories for non-equilibrium statistical mechanics, where phase transitions among stationary states far away from thermodynamic equilibrium could be studied. Instances of such systems have recently been demonstrated in a variety of contexts: In ensembles of ultracold atoms, Bose-Einstein condensates (BEC) placed in optical cavities have allowed to achieve strong lightmatter coupling, and led to the realization of open Dicke models 2,3 . The corresponding phase transition has been studied in real time, including the determination of the associated critical exponent 4 . In systems of trapped ions 5 , Ising models with variable-range interactions of a few hundred quantum spins have been created 6 . Other platforms, which hold the promise of being developed into true many-body systems by scaling up the number of presently existing elementary building blocks in the near future, are provided by arrays of microcavities [7][8][9][10] , and also optomechanical setups [11][12][13] .
Genuine many-body ensembles in the above class are furthermore realized in the context of pumped semiconductor quantum wells in optical cavities 14 . Here, non-equilibrium Bose-Einstein condensation of exciton-polaritons has been achieved 15-17 -the effective bosonic degrees of freedom result from a strong hybridization of cavity light and excitonic matter states 1, 18,19 .
All these systems exhibit the crucial ingredients for non-trivial critical scaling behavior at a continuous nonequilibrium phase transition. This triggers broader theoretical questions on the actual nature and possible universality classes of such non-equilibrium critical points. At first sight, invoking the concept of universality, implying a huge "loss of memory" on details of the microscopic physics, it may seem questionable whether the microscopic non-equilibrium conditions will result in any physically observable consequences at the macroscopic level at all. In particular, for equilibrium dynamical critical behavior there exists a well-developed theoretical framework based on the seminal work of Hohenberg and Halperin (HH) 20 (and other authors), who classified various types of dynamical critical behavior into diverse equilibrium dynamical universality classes, known as Models A to J, depending on the conserved or non-conserved nature of the order parameter itself, and on its dynamical couplings to other slow conserved modes.
However, there are two key ingredients, shared by the systems described above, that place these many-body ensembles apart from equilibrium systems and may in fact cause novel universal physical features. First, they are strongly driven by external fields, such as coherent electromagnetic radiation provided by lasers, and undergo a cascade of internal relaxation mechanisms 1 . This non-equilibrium drive and balancing dissipation adds to the Hamiltonian dynamics, and causes both reversible (coherent) and irreversible (dissipative) dynamics to appear on an equal footing, albeit originating from physically distinct and independent mechanisms. In turn, this induces manifest violations of the detailed-balance conditions characteristic of a many-body system in thermal equilibrium. Indeed, these drive-induced non-equilibrium perturbations transcend mere violations of the Einstein relations (or fluctuation-dissipation theorem) that in equilibrium connect relaxation coefficients with associated thermal noise strengths: Such perturbations have been found to generically become irrelevant in the vicinity of a second-order phase transition (for a concise overview, see Ref. 21). Second, these systems are characterized by the absence of the conservation of particle number. This is due to the admixture of light to the matter constituents, which opens up strong loss channels for the effective hybrid light-matter degrees of freedom, in turn making it necessary to counterpoise these losses by continuous pumping in order to achieve stable stationary states.
Deciding whether or not these ingredients indeed cause universal behavior distinct from equilibrium motivates -and in fact necessitates -a thorough theoretical analysis of the nature of criticality in such non-equilibrium quantum systems.
A key representative of potential non-equilibrium criticality is provided by the driven-dissipative Bose-Einstein condensation transition, relevant to the experiments with excitonpolariton condensates described above. For such systems, indeed, a new independent critical exponent associated with the non-equilibrium drive has recently been identified within a functional renormalization group (RG) approach 22,23 . This exponent describes universal decoherence at long distances, and is observable, e.g., in the momentum-and frequencyresolved single-particle response, as probed in homodyne detections of exciton-polariton systems 24 . Furthermore, an effective thermalization mechanism for the low-frequency distribution function has been found, reflected in an emergent symmetry at the (classical, equilibrium) Wilson-Fisher fixed point.
In this work, we employ the field-theoretical RG [25][26][27][28] in a perturbative dimensional ǫ expansion for a complementary study of driven Bose-Einstein criticality. Here ǫ = 4 − d measures the distance from the upper critical dimension d c = 4, and serves as the effective small parameter in the perturbation series. In this framework, we confirm yet also further develop the physical picture obtained previously within the functional RG approach. Both our perturbative two-loop and the non-perturbative functional RG analysis 22,23 are based on an effective long-wavelength description in terms of a noisy Gross-Pitaevskii equation with complex coefficients [29][30][31][32][33] , which in turn constitutes a variant of the time-dependent complex Ginzburg-Landau equation for a two-component order parameter field 34 . Such complex stochastic differential equations have also found extensive applications in the modeling of spontaneous structure formation in non-equilibrium systems 35,36 . Remarkably, these comprise coupled non-linear oscillators subject to external noise near a Hopf bifurcation instability 37 , and even spatially extended evolutionary game theory and the dynamics of cyclically competing populations 38 .
We present a concise account of the main results of this work in the following Sec. II. The remainder of this paper is organized as follows: In Sec. III, we explain the microscopic model based on a stochastic Gross-Pitaevskii equation with complex coefficients, and introduce the equivalent dynamical response functional integral as appropriate for a subsequent diagrammatic evaluation. We also discuss the relationship of our model with Model E that governs the equilibrium critical dynamics of planar ferromagnets and the normal-to superfluid phase transition. Section IV comprises the bulk of this work. It contains an explanation of the renormalization scheme employed to deal with the emergent ultraviolet (UV) divergences, and details how the critical scaling properties in the infrared (IR) region may subsequently be obtained from solutions of the associated RG flow equations. Finally, Sec. V offers a summary and concluding remarks, and supplementary appendices provide more technical details.

II. KEY RESULTS AND PHYSICAL PICTURE
Generalized dynamic scaling forms. -Near the continuous condensation transition for driven-dissipative boson system, we derive generalized scaling laws for the dynamic response and correlation functions at wavevector q and frequency ω: where a is a non-universal constant, and the correlation length diverges as ξ ∝ |τ| −ν as the critical point is approached, Here, ν, η, and z represent the standard equilibrium static and dynamical critical exponents, while η c and η ′ constitute novel scaling exponents induced by the nonequilibrium drive and associated potential violation of detailed balance. The origin of these new scaling exponents is immediately transparent from the description of the problem in terms of a Janssen-De Dominicis (or Martin-Siggia-Rose) functional integral [39][40][41] : Owing to the competition of coherent and dissipative dynamics in the driven problem, two independent mass scales appear, as compared to a single one in the closely related, purely relaxational Model A in equilibrium critical dynamics. This causes a more complex critical scenario akin to a bicritical point. In addition, the fluctuation-dissipation theorem that relates the dynamic response with the correlation function in thermal equilibrium (for which η ′ = η are hence identical) is in general violated. Indeed, there arise new ultraviolet divergences, specifically for two couplings that are marginal at the Gaussian fixed point.
Asymptotic thermalization. -We establish that the renormalization group flow, already to one-loop order, drives the system towards an effectively equilibrium fixed point, where detailed balance is satisfied. The fluctuation-dissipation theorem then implies that holds exactly for the Fisher exponents that characterize the anomalous algebraic spatial decay of the order parameter dynamic response (1) and correlation (2) functions at criticality. The analysis of the two-loop RG flow equations furthermore yields that the asymptotic fixed point values of all nonequilibrium coupling parameters induced by the external drive vanish. Consequently the static critical exponents are precisely those of the equilibrium O(2)-symmetric Ginzburg-Landau-Wilson Hamiltonian (XY model), namely to lowest non-trivial order in the dimensional ǫ expansion, and the associated dynamic critical exponent is the standard one for the O(2)-symmetric Model A that describes the purely relaxational kinetics of a non-conserved two-component order parameter field 40 , Thus, the system's asymptotic long-wavelength and lowfrequency properties become effectively thermalized; see also Ref. 37. Yet a novel universal scaling exponent appears in the subleading scaling behavior, see Eqs.
(1) and (2), which originates from the driven-dissipative setup 34 . Novel drive exponent. -Despite the fact that the system approaches an effectively equilibrium RG fixed point, one of the marginal couplings induces a novel and independent (but subleading) scaling exponent that captures the fadeout of coherent quantum fluctuations relative to their thermal, dissipative counterparts. We compute the new drive exponent to second order in the loop expansion, i.e., to order ǫ 2 : Eq. (7) is one of the central results of this paper. Together with the asymptotic thermalization, these key findings corroborate the earlier functional RG study of Ref. 22 and 23. They also underscore the well-known remarkable stability of Model A with respect to non-equilibrium perturbations [42][43][44] .
Hierarchical structure of non-equilibrium criticality. -In this way, we confirm the hierarchical structure of the model's critical behavior, in the following sense: The static critical behavior is characterized by the O(2) universality class, described by the rotationally invariant Ginzburg-Landau-Wilson Hamiltonian for a two-component order parameter field. In equilibrium dynamical criticality, the static properties are supplemented -but not modified -by the dynamical critical exponent, e.g. Eq. (6) for Model A. Yet here, the nonequilibrium conditions give rise to the new and independent scaling exponent (7). We thus establish the following pattern: While the non-equilibrium drive modifies neither the universal O(2) static nor even the dynamical critical behavior of Model A, it still adds novel universal scaling features. This situation is reminiscent of (but different from) the emergence of a new critical exponent in Model A associated with the non-equilibrium relaxation following a sudden temperature quench from random initial conditions to the critical point 45,46 .
Relation to equilibrium dynamic criticality. -We elaborate on the relation of the driven Bose-Einstein condensation to its equilibrium counterpart, which is described by Model E in the terminology of HH. The hydrodynamic conservation laws relevant for the latter model have two crucial consequences, which set it apart from our non-equilibrium situation: First, the dynamical critical exponent is modified due to the existence of a new relevant reversible coupling to a diffusive mode, and fixed by rotational invariance in order parameter space to z = d/2 (in the strong dynamic scaling regime), quite distinct therefore from our result for the relaxational dynamical critical exponent coinciding with Model A. Second, these conservation laws exclude the addition of a second mass scale to the problem, and in consequence, there emerges no counterpart of the drive exponent in equilibrium Bose-Einstein condensation criticality.
Intriguingly, the additional drive exponent is absent for a Model B version (in the HH classification) of the complex Ginzburg-Landau equation. Instead of the purely relaxational kinetics of a non-conserved order parameter, in this situation a diffusive relaxation for a conserved order parameter field is implemented. Thermalization along with the conservation law and ensuing structure of the non-linear relaxation vertices now imply the exact scaling laws (3) and 40 In addition, we establish the identity at least to two-loop order. Theoretical approach and renormalization scheme. -The perturbative field-theoretical RG approach 25-28 provides a well-established tool for the quantitative characterization of critical behavior close to the upper critical dimension d c = 4; more precisely, it is perturbatively controlled in the dimensional parameter ǫ = 4 − d. In equilibrium it has moreover been demonstrated that the structure of the RG flow equations and the ensuing universality classes remain robust even in extensions down to three dimensions. In this paper we work at lowest non-trivial order in ǫ, i.e., to O(ǫ) for, e.g., the correlation length exponent ν and the fixed-point value of the non-linear coupling u, but to order ǫ 2 for the Fisher, dynamic, and drive exponents, whose anomalous scaling dimensions require a calculation at the two-loop level. A key advantage of this approach in the context of non-equilibrium criticality is the possibility of a direct and quantitative comparison of our findings with the well-known results for equilibrium dynamical criticality displayed by the phenomenological models of HH 20 , which have not yet been comprehensively studied in a functional RG framework (except for Refs. [47][48][49]. We remark that the field-theoretical RG approach differs conceptually from the functional RG based on Wetterich's equation 50 : The latter constitutes an exact reformulation of a given functional integral in terms of a functional differential equation, in this way at least in principle addressing the full many-body problem [51][52][53][54][55] . Critical behavior can then be studied by a suitable fine-tuning of parameters. In contrast, the field-theoretical RG focuses immediately on the critical surface of the problem, in this way isolating the universal critical behavior from the outset (see, e.g., Refs. 25-28; and for the application to dynamic critical phenomena Refs. 34,40,[56][57][58][59]. While, therefore, non-universal aspects of the problem are projected out, it provides a perhaps more fundamental understanding of the emergence of scaling properties, and moreover allows us to obtain explicit analytical results for the critical exponents, cf. Eqs. (4)- (7).
In short, by means of the field-theoretical RG approach we provide complementary strong evidence that the microscopic non-equilibrium character bears observable consequences up to the largest distance and time scales in driven-dissipative Bose-Einstein condensation.

A. The dissipative Gross-Pitaevskii equation with noise
Driven-dissipative Bose-Einstein condensation in excitonpolariton systems is properly described by a noisy dissipative Gross-Pitaevskii equation with complex coefficients [29][30][31][32][33] , It basically coincides with the time-dependent complex Ginzburg-Landau equation, which has been prominently employed to describe pattern formation in non-equilibrium systems, typically however in the deterministic limit without noise 35,36 . A stochastic variant has been analyzed in the context of coupled anharmonic oscillators 37 . Here, the complex bosonic field ψ describes the polariton degrees of freedom. The complex coefficients have clear physical meanings; in particular, χ = (γ p −γ l )/2 is the net gain, i.e., the balance of the incoherent pump rate γ p and the local single-particle loss rate γ l . The positive parameters κ and λ represent the two-body loss and interaction strength, respectively, while A = 1/2m eff relates to the effective mass of the polaritons. Typically, this equation is not presented with an explicit diffusion coefficient D, but rather with a frequency dependent pump term ∼ η∂ t ψ adding to the left hand side of the equation 60,61 . The form (10) is then recovered upon division by 1 − iη, i.e., with D = Aη and subleading corrections to the other coefficients which are complex to begin with. We emphasize that, due to the freedom of normalizing the time derivative term as above in the equation of motion, this model accurately captures the physics close to the phase transition, since it describes the most general low-frequency dynamics in a systematic derivative expansion that incorporates all relevant couplings (in the sense of the RG) in dimensions d > 2. Finally, the noise described by the fluctuating complex variable ζ is taken to be Gaussian, white, and Markovian, and hence is fully characterized by the correlators Physical stability requires A, D, λ, and κ to be positive. The parameter χ is negative in the disordered phase, where the global U(1) gauge symmetry (or O(2) rotational symmetry in the complex plane) is not spontaneously broken. Our calculations will be carried out in this regime, i.e., we approach the phase transition from the disordered side, in contrast to the non-perturbative RG analysis in Refs. 22 and 23. For increasing pump rate γ p , the gain χ eventually turns positive and the system undergoes a continuous, driven Bose condensation transition: The instability occuring for a state with vanishing polariton field expectation value is cured by the expression of a polariton condensate ψ(x, t) 0. The parameter µ, which effectively assumes the role of a chemical potential, is fixed by the requirement of stationarity, see below. The Langevin equation (10) can be formally derived from a microscopic description in terms of a quantum master equation (see, e.g., Ref. 62) upon employing canonical power counting in the vicinity of the critical point 22,23 .
For the analysis of the critical behavior, it turns out useful to introduce the following ratios (and sign conventions): The relaxation rate D may be small on the microscopic scales of actual experiments, but it plays a key role for the dominantly diffusive dynamics in the vicinity of the phase transition, as will be confirmed in the subsequent calculation. Therefore, it is useful to express all external control parameters relative to D. The quantities r K and r U are dimensionless, giving the ratio of real and imaginary parts of the couplings in Eq. (10), and therefore describing the relative strength of coherent vs. dissipative dynamics. In these units Eq. (10) takes the form with the stochastic noise ξ(x, t) = −iζ(x, t) governed by the same correlations (11) as ζ. Formally, as indicated in Eq. (14), this describes the relaxational kinetics of Model A with a nonconserved order parameter, however with a non-Hermitean effective "Hamiltonian" The complex coefficients in Eq. (15) reflect the presence of the non-equilibrium drive. The theory becomes critical (massless) when r, r ′ → 0 (more precisely, the renormalized counterparts of these parameters τ, τ ′ → 0) simultaneously. We remark that Ref. 37 addressed the distinct physical situation where the uncoupled oscillation frequence r ′ was held fixed. The calculation was then performed in a rotating reference frame, which formally amounts to setting r ′ = 0 in our analysis.

B. Field theory representation
The nonlinear partial differential equation (10) or (14) represents a classical stochastic evolution, which is readily mapped into an equivalent Janssen-De Dominicis (or Martin-Siggia-Rose) functional integral representation [39][40][41] ; for detailed explanations, see, e.g., Refs. 56, 58, 59, and 63. This formulation renders it amenable to straightforward perturbative expansions with respect to the nonlinear coupling u ′ , the use of diagrammatic techniques, and subsequent implementation of the field-theoretical dynamical RG. The Janssen-De Dominicis response functional corresponding to Eq. (14) with noise correlations (11) reads It provides the statistical weight for the stochastic process encoded in the Langevin equation (13) for ψ(x, t). The associated generating function for the dynamic correlation functions and cumulants becomes We note that Z[j = 0, j = 0] = 1 carries no information, in stark contrast to the partition function in thermal equilibrium. The perturbative expansion proceeds around the Gaussian action (u ′ = 0). With the Fourier transform convention and analogously for the response fieldψ(x, t), the Gaussian action reads in frequency-momentum space: with the Hermitean harmonic coupling matrix where R(q) = r+ir ′ +(1 + ir K ) q 2 . Inversion of the 2×2 matrix yields the bare advanced and retarded response propagators as well as the correlation propagator; explicitly, these read: These expressions can be written in scaling form, One may set up the diagrammatic perturbation expansion either with these three propopagators, or equivalently just with the response propagators (21) and the two-point noise vertex Γ 0 ψ * ψ = γ/2, in addition to the nonlinear four-point vertices (1 + ir U ) (computed at symmetrized incoming external wavevectors). The graphical representations for these elements of the perturbation series are depicted and explained in Fig. 1.
The noise correlators (11) imply for η 1 = Re ξ and η 2 = Im ξ: In Eq. (24), the systematic forces in the Langevin equations have been decomposed into the dissipative, relaxational term with the standard O(2)-symmetric Ginzburg-Landau-Wilson Hamiltonian and the reversible contribution with a second Ginzburg-Landau-Wilson Hamiltonian In the mean-field approximation, we merely need to simultaneously minimize both H and H ′ to obtain possible stationary configurations. For the temperature-like control parameter r > 0 (net gain χ < 0), the only homogeneous state is S = 0, or ψ = 0, describing the disordered phase, i.e., the absence of a Bose-Einstein condensate. For r < 0 (χ > 0), on the other hand, we encounter the ordered phase with finite condensate fraction, namely from minimizing H with a constant | S | = √ 6|r|/u ′ = χ/κ = |ψ|. Minimizing the effective Hamiltonian H ′ yields the second condition | S | = √ 6|r ′ |/r U u ′ = µ/λ = |ψ|, whence consistency requires that indeed r ′ = r U r. The chemical potential then adjusts itself to µ = λ|ψ| 2 . In effect, this leaves r as the sole control parameter for the condensation transition.
We may now consider the following two special cases: (i) For parameters r ′ = r K = r U = 0, H ′ vanishes, whence we recover the O(2)-symmetric Model A for purely relaxational critical dynamics towards thermal equilibrium, if we impose Einstein's relation (or rescale the fields appropriately) with an effective temperature T . We shall later employ this parametrization also in a strictly non-equilibrium setting (with Boltzmann's constant set to k B = 1). A distinct scaling behavior of the noise strength γ and the relaxation rate D, and hence the "temperature" T indicate a violation of detailed balance.
(ii) For r ′ = r U r, r K = r U 0, H ′ = r K H, one arrives at an effective equilibrium dynamics with reversible term Its antisymmetry ensures that the associated reversible probability current remains divergence-free in the space of dynamical variables S α . This special situation has been analyzed in Ref. 34. Intriguingly, this kinetics resembles the critical dynamics of Model E for a non-conserved two-component order parameter field (e.g., the in-plane magnetization fluctuations for an XY ferromagnet), reversibly coupled to a conserved scalar field M (corresponding to the z-component of the magnetization in a planar ferromagnet; c.f. Refs. 20 and 59). Here, however, M ∝ Dr K is spatially uniform and stationary. We remark in passing that the dynamic critical exponent of the equilibrium Model E is fixed by the fact that the conserved field M generates rotations in order parameter space. Under the assumption of strong dynamic scaling, i.e., proportional divergent time scales for the critical modes and the conserved quantity, one obtains z = d/2 exactly in dimensions d ≤ 4. Indeed, this constitutes a crucial difference between the equilibrium and the driven models: The uniform magnetization in the driven case does not scale, whereas the slowly varying magnetization field in the equilibrium case does. This distinction causes the dynamic critical exponent in the driven and equilibrium cases to differ markedly.
Adding an external field term to the Hamiltonian, , yields in this special case (ii) the dynamic susceptibilities whereS 1/2 represent the (real) Martin-Siggia-Rose response fields associated with S 1/2 . For the original complex fields, these components combine to since the effective Onsager coefficient is D (1 + ir K ). As the system is in thermal equilibrium, the fluctuation-dissipation theorem relates the dynamic response with the correlation function through or equivalently in Fourier space For both special situations (i) and (ii), and with Eq. (30) the scaling forms (22) for the retarded response propagator and dynamical correlation functions reduce to With Eqs. (33) and (34), i.e., χ(q, ω) = D (1 + ir K ) Gψ * ψ (q, ω) and C(q, ω) = G ψ * ψ (q, ω), Eqs. (37) satisfy the fluctuationdissipation theorem (36). Fourier transformation to the time domain yields explicitly which likewise fulfill Eq. (35).

A. Renormalization scheme for ultraviolet divergences
The perturbation expansion is most conveniently carried out for the one-particle irreducible vertex functions, since redundancies are thus eliminated in the calculations. The generating functional for the vertex functions is related to its counterpart (17) for the connected correlation functions (cumulants) in the standard manner through a Legendre transformation [25][26][27]59 . Here, we merely list the explicit relationships between the two-point vertex functions and cumulants in Fourier space, where the second expression originates from Dyson's equation with the associated self-energy Σ; furthermore and similarly for the four-point functions, etc.
In the field-theoreticalversion of the renormalization group approach to critical phenomena [25][26][27][28]56,57,59 , one sends all ultraviolet (UV) cutoffs originating from the short-distance physics to infinity. At and above the upper critical dimension (here, d c = 4) UV divergences appear in the perturbation expansion that are absorbed into appropriately defined renormalized parameters. In the vicinity of an RG fixed point, where scale invariance ensues, one may then infer the desired infrared (IR) scaling properties of the theory from its UV behavior, which is perturbatively accessible, provided one ensures to work outside the IR-singular critical region (which here is defined by r, r ′ → 0 in the unrenormalized theory along with q, ω → 0). The field theory action (16) entails UV divergences for the vertex functions Γψ ψ * , Γψψ * , and Γψ ψ * ψ * ψ .
The propagator self-energy Σψ * ψ contains quadratic UV divergences (at d c = 4) that first need to be additively renormalized. Physically, this corresponds to a fluctuation-induced downward shift of the critical point (pump rate), τ = r − r c , and of course the subsequent characterization of the critical behavior needs to address the vicinity of the true phase transition point at τ = 0; i.e., r c is determined from the condition Γψ ψ * (q = 0, ω = 0) = 0 at r = r c .
The remaining logarithmic UV divergences are then multiplicatively absorbed into renormalization factors for which we choose the following conventions: We define the renormalized counterpart to the two-point vertex function (39) as The complex renormalization constant Z then follows from the singular part of its frequency derivative, evaluated at the normalization point τ = µ 2 with arbitrary momentum scale µ, but manifestly outside the IR-singular critical regime. Next we introduce dimensionless renormalized counterparts to the parameters defined in Eqs. (12) and (30) via multiplicative renormalization with real Z factors: Here, A d = Γ(3 − d/2)/2 d−1 π d/2 denotes a geometric factor; Γ(x) indicates Euler's Gamma function. Note that independent renormalization constants Z T , Z r K , and Z r U only arise in a genuine non-equilibrium setting; the equilibrium theory can be fully renormalized through a real Z along with Z τ , Z D , and Z u ′ . Indeed, the fluctuation-dissipation theorem (35) or (36) in conjunction with (40) and the definitions (41), (43) implies that must hold in thermal equilibrium. An independent Z T factor means that the effective "temperature" becomes scaledependent in the driven non-equilibrium system. In contrast, the relation (44) reflects the partition invariance of temperature in the renormalization group language: all arbitrary system partitions must be in equilibrium with each other. Its origin can be traced back to a specific equilibrium symmetry of the response functional 23,59 .
The renormalization constants will be determined perturbatively to lowest non-trivial order in the non-linear coupling u ′ , namely first Z D and Z r K from ∂ q 2 Γψ ψ * (q, ω = 0)| sing.

B. Renormalization group equation and RG flow functions
The renormalization group equation exploits the fact that the unrenormalized quantities do not depend on the arbitrary momentum renormalization scale µ. Translated to renormalized correlation or vertex functions, it relates their properties at different momentum (or length, time) scales, and thus provides the desired link between the theory in the ultraviolet, where perturbative computations can safely be carried out, and the physically interesting infrared region governed by non-trivial critical singularities [25][26][27][28]59 . Denoting the set of model parameters as {p} = {D, τ, T, r K , r U }, and introducing u = u ′ T , which turns out to be the proper effective non-linear coupling in the perturbation series, one obtains for example for the two-point vertex function Here we have defined Wilson's flow functions -note that ζ is complex -and the RG beta function The partial differential equation (45) is readily solved by the method of characteristics µ → µℓ, which leads to decoupled first-order ordinary differential flow equations for the running parameters and the coupling The infrared limit is attained as ℓ → 0. Near an infrared-stable RG fixed point given by the zero of the beta function (47), i.e., β u (u * ) = 0 with ∂ u R β u | u * > 0, the model becomes scaleinvariant. The solutions of the flow equations for the running couplings then become simple power lawsp(ℓ) ≈ p R ℓ γ * P , with the anomalous scaling dimensions γ * P = γ P (u * ). Since all perturbative contributions to the vertex and correlations functions consist of integrals over products of Gaussian propagators and vertices, we observe that the renormalized two-point function takes the form Γ R ψψ * (q, ω, {p R }, u R ) = D R µ 2 Γ q 2 (1 + ir KR )/µ 2 , ω/(D R µ 2 ), τ R (1 + ir UR ), u R , c.f. Eq. (22). One thus finally arrives at the asymptotic solution of the RG equation (45) near a fixed point u * , One may similarly proceed for any other vertex function, or, u q q k k q + Figure 2. One-loop propagator renormalization. The loop diagram depicts the first-order correction to the self-energy or two-point vertex function Γψ ψ * (q, ω).
e.g., the dynamical correlation function

C. Renormalization and scaling to one-loop order
We now carry out the explicit perturbational analysis of the fluctuation corrections to first order in u, represented through Feynman diagrams with a single closed propagator loop. For the two-point vertex function Γψ ψ * (q, ω), the corresponding one-particle irreducible graphs are shown in Fig. 2. The loop diagram represents the lowest-order contribution to the associated self-energy −Σψ * ψ (q, ω), see Eq. (39). The ensuing analytic expression is explicitly Upon performing the internal frequency integral via Cauchy's theorem, the integrand simplifies considerably. Setting γ = 4DT (henceforth we employ units where Boltzmann's constant k B = 1), r ′ = r U r and u = u ′ T , one arrives at which is negative and represents a downward shift of the critical point: fluctuations suppress spontaneous long-range order. By means of (A2) in App. A one arrives at a self-consistent equation for |r c | which to this order is solved by . (55) Notice that r c (u) diverges upon approaching the lower critical dimension d lc = 2, correctly indicating that the critical point is driven towards zero and there emerges no truly longrange order with spatially homogeneous condensate (in fact, an isotropic driven-dissipative Bose gas in two dimensions cannot even support quasi long-range order, see Ref. 64). Furthermore, both non-equilibrium parameters r K and r U have dropped out: Eq. (54) just represents the equilibrium critical point shift. Indeed, since the fluctuation loop in Eq. (52) is proportional to the factor 1 + ir U , we may define the true distances from the critical point τ = r − r c and τ ′ = r U (r − r c ) = r U τ, which vanish simultaneously as r → r c (u), even when fluctuation effects are included (to first order), thus preserving the basic mean-field scenario.
In terms of τ, since r c = O(u) one can now rewrite Eq. (52) to this order: Just as in thermal equilibrium, the UV singularities contained in the wavevector integral can be entirely absorbed into a multiplicative renormalization of the parameter τ, see Eq. (43): where we have used (A3) in App. A.1, and in the final step applied the minimal subtraction scheme, wherein only the 1/ǫ pole and its residuum at the upper critical dimension d c are included in the renormalization constant Z τ . Since in addition there exists no one-loop correction to the noise vertex, see Fig. 4 below, i.e., Γψ * ψ(q, ω) = −2DT + O(u 2 ), we infer that to this order The associated Wilson RG flow functions (46) and anomalous dimensions all vanish in the one-loop approximation, and the sole non-trivial RG flow function to first order in u R is In order to compute the beta function (47) for the non-linear coupling u and obtain the RG fixed points, we require the renormalization of the four-point vertex function Γψ ψ * ψ * ψ . Its tree and one-loop contributions are depicted in Fig. 3. Carrying out the internal frequency integrals and combining the three loop contributions, one finally arrives at To this order in u, we may replace r with τ in the integrals, and r ′ = r U r with τ ′ = r U τ, and evaluate at the normalization point τ = µ 2 (τ R = 1, safely outside the IR-singular region), with q = 0, ω = 0. By means of Eqs. (A3) and (A4), we obtain in minimal subtraction: Separating out the real and imaginary parts yields the renormalization constants Notice that consistency with the one-loop analysis of the propagator self-energy that led to Z r U = 1 + O(u 2 ) demands that r * U = r * K at the stable RG fixed point. With Z T = 1 + O(u 2 ), the RG beta function (47) becomes where ∆ R = r UR − r KR . For d > d c = 4 (ǫ < 0), its only stable zero is the Gaussian fixed point u * 0 = 0 (with ∂ u β u | u * 0 =0 = −ǫ), which implies mean-field scaling exponents ν = 1/2, η = 0, z = 2, and also η c = 0 = η ′ . At the upper critical dimension, the RG flow tends to zero only slowly, inducing logarithmic corrections to the mean-field power laws. In dimensions d < 4, a non-trivial fixed point emerges: (provided the denominator is positive). It depends parametrically on r KR and the difference ∆ R , and leads to non-Gaussian critical exponents. We next study the RG beta function associated with ∆ R . Since γ r K = 0, we obtain from Eq. (63) At the Gaussian fixed point u * 0 = 0, any constant values of r KR and r UR are allowed. For d < 4, at the non-trivial, positive, and stable fixed point (65), the only real zero of Eq. (66) is indeed ∆ * = 0, whence, as anticipated, r * U = r * K , and the RG fixed point becomes independent of r K : This is just the equilibrium XY model fixed point for the O(2)symmetric Ginzburg-Landau-Wilson Hamiltonian. We note that the stability matrix eigenvalues at the infrared-stable RG fixed point (67) are ∂ u R β u | ∆ * =0,u * = ǫ and ∂ ∆ R β ∆ | ∆ * =0,u * = ǫ/5. Therefore the RG flow will typically first approach the nonequilibrium fixed line (65), and subsequently tend towards the equilibrium fixed point ∆ R → 0 along this critical surface. At this point, the system has already become effectively thermalized, and is described by the special case (ii) discussed in Sec. IIIC; thus Eq. (44) holds, albeit trivially to one-loop order, see (58). With the one-loop results (59), the solutions (49) and (50) of the RG equations for the inverse response propagator and the dynamical correlation function simplify drastically at the one-loop equilibrium fixed point. According to Eq. (37) and applying the matching condition ℓ = |q|/µ, they reduce to the following scaling laws for the dynamical susceptibility: where we have omitted fixed, constant arguments and identified the inverse correlation length exponent and for the dynamical correlation function   Figure 4. Two-loop renormalization of the noise vertex: There is no fluctuation correction to first order in u; the two-loop graph represents the lowest-order contribution to the vertex function Γψ * ψ(q, ω).

D. Two-loop analysis and renormalization
In order to obtain non-trivial dynamic critical and drive exponents, we need to proceed to the next order in the perturbational and dimensional ǫ expansion. The two-point noise vertex is only renormalized to two-loop order, as shown in Fig. 4. Carrying out both internal frequency integrals associated with the closed propagator loops, one arrives to second order in the non-linear coupling u at Γψ * ψ(q, ω) = −2DT 1 + 2 9 u 2 1 + r 2 Setting r ′ = r U r, and replacing r = τ + O(u) in the integrands, we find Γψ * ψ(0, 0) = −2DT 1+ 2 9 u 2 1 + r 2 U Re I d (τ, r K , r U )+O(u 3 ) , (72) with the nested wave vector integral Evaluating this integral at the normalization point τ = µ 2 , and isolating its UV divergences in the form of 1/ǫ poles then yields the Z factor product (in minimal subtraction) The two-loop Feynman graphs contributing to the retarded response propagator self-energy or vertex function Γψ ψ * (q, ω) are depicted in Fig. 5. The first two closed diagrams (on the left) yield contributions that are independent of the external wavevector q and frequency ω, and combine to . Two-loop one-particle-irreducible Feynman diagrams contributing to the propagator self-energy or two-point vertex function Γψ ψ * (q, ω). The two graphs on the right induce non-classical values for the exponents η, z, η c , and η ′ .
The two graphs to the right, with the wavevectors distributed as indicated in Fig. 5, yield after internal frequency integration −D 2 9 u 2 (1 + ir U ) Symmetrizing with respect to the internal wavevectors k ↔ p ↔ q − k − p simplifies this expression markedly, and the sum of Eqs. (52), (75), and (76) can be written as For vanishing external wavevector and frequency, we obtain with r ′ = r U r: At the stable RG fixed point (67) with r * U = r * K , the righthand side of Eq. (78) reduces to the standard two-loop additive and multiplicative temperature renormalizations for the mass parameter r. The fluctuation-induced T c shift, as well as the renormalization constant Z τ and hence the correlation length exponent ν, remain identical to those of the XY model or O(2)-symmetric Model A in thermal equilibrium to this order. The remaining multiplicative renormalization factors follow from the frequency and wavevector derivatives of Eq. (77) at the normalization point τ = r + O(u) = µ 2 : According to Eq. (42) whence subsequently Z D and Z r K can be determined from the singular contributions to The ultraviolet singularities to be captured in Z r K encode a novel scaling exponent that describes the weight and fadeout of coherent quantum fluctuations relative to their thermal, dissipative counterparts. Appendix A.1 details how the UV-singular part is extracted from this nested wavevector integral in the form of a simple 1/ǫ pole and its residuum, applying dimensional regularization with minimal subtraction. Thus, Eq. (80) yields with Eq. (A9): (82) The evaluation of the integral (73), detailed in Apps. A.1 and A.2, gives where s k = 1 + 4 9 r 2 K and the logarithmic function L(r K ) is given in Eq. (A16). Finally, at the infrared-stable fixed point where r * U = r * K : Carefully separating the real and imaginary parts in Eq. (82) and inserting Eqs. (83), (84) allows us to compute the desired renormalization constants to two-loop order: and at last, Eq. (85) gives with Eq. (87): E. Scaling and critical exponents to order ǫ 2 In order to explore possible fixed points, we compute the RG beta function for the non-equilibrium parameter r K , using Eqs. (43) and (86): At the equilibrium fixed point r * K = 0, we have ∂ r KR β r K | r * K =0 = 2 9 ln 4 3 u 2 R > 0, whence it is infrared-stable. Indeed, as shown in Fig. 6, β r K (r KR ) is a monotonically growing function of the renormalized non-equilibrium parameter r KR , and r * K = 0 is its only zero. Setting r KR = 0, Wilson's flow functions (46), readily derived from Eqs. (83), (84), (86), (87), and (88), simplify drastically: Once the system has reached a thermalized state, which here happens already at the one-loop level, the exact equilibrium relation (44) enforces the identity with a real ζ, which is in fact satisfied by the explicit two-loop results (90). Employing again the matching condition ℓ = |q|/µ to the RG equation solutions (49) and (50), the universal scaling forms (68) and (70) become now generalized to see Eqs.
(1) and (2), that are valid in the vicinity of the critical point. Here we have identified the set of universal scaling exponents as Note that η − η c = γ * r K is determined by the anomalous scaling dimension of the non-equilibrium parameter r KR . Since z − 2, η, η c , and η ′ all vanish to first order in the non-linear coupling, we only need to insert the stable one-loop equilibrium fixed point value (67) to finally recover the standard equilibrium Model A critical exponents (5) and (6) to two-loop order: = 2 + ǫ 2 50 6 ln 4 3 − 1 + O(ǫ 3 ) .
Since the system has reached an equilibrium fixed point and is effectively thermalized, the identity (91) immediately implies the exact scaling relation reflecting the emergence of detailed balance and the ensuing fluctuation-dissipation theorem. This leaves us with a single new non-equilibrium drive scaling exponent If we (daringly) set ǫ = 1, we find the critical exponents ν ≈ 0.625 (to one-loop order), η = η ′ ≈ 0.02, z ≈ 2+0.72609 η ≈ 2.01452, and η c ≈ −0.15073 η ≈ −0.0030146. For comparison, the numerical values found in d = 3 dimensions by means of the non-perturbative RG approach in Refs. 22 and 23 are ν ≈ 0.716, η ≈ 0.039, z ≈ 2.121, and η c ≈ −0.223. The two-loop values thus apparently still underestimate the fluctuation corrections in three dimensions.

F. Observability of the drive exponent
The fact that the drive exponent appears in the scaling form of the single-particle dynamical response makes it accessible to experimental observation 22,23 . In particular, the imaginary part of the dynamical response is probed in radiofrequency spectroscopy in ultracold atoms 65 , or homodyne detection for exciton-polariton systems 24 (in the latter, also the real part is available separately). For probe frequency ω, the scaling form (1) implies for the response at criticality which demonstrates different critical wavevector scaling for the peak position ω 0 ∼ |q| z+η−η c and peak width ∼ |q| z .

G. Complex Model B dynamics
We conclude our field-theoretical approach with a brief RG analysis of a variant of the Langevin equation (14) with conserved order parameter dynamics, i.e., off-critical diffusive relaxation (Model B): with the non-Hermitean "Hamiltonian" (15), and the noise correlations with γ = 4DT such that the standard O(2)-symmetric equilibrium Model B critical dynamics is incorporated as the special case with r ′ = r K = r U = 0. When these parameters are non-zero, we have a conserved complex Ginzburg-Landau equation, or a driven, non-equilibrium version of the noisy Cahn-Hilliard equation with complex coefficients.
As a consequence of the conserved dynamics, both the Onsager relaxation coefficient in Eq. (99) and the noise correlator strength (100) now carry an additional Laplacian operator. This increases the mean-field dynamic critical exponent to z = 4, and generates an external wavevector factor q 2 attached to the outgoing legs in the non-linear relaxation vertices depicted in Fig. 1(c). Therefore, one has to all orders in the perturbation expansion Since these vertex functions carry no UV singularities, thus Z = 1 and Z T = Z −1 D , which imply the exact results Note that these Model B relations represent a special case of the general equilibrium condition (91). In addition, according to Eq. (82) Z = 1 enforces to two-loop order that . For the dynamical susceptibility χ(q, ω) = Dq 2 Γψ ψ * (q, ω) −1 and correlation function C(q, ω) = G ψ * ψ (q, ω), one is led again to the scaling laws (1) and (2); Here, the static critical exponents ν = −1/γ * τ and η = −γ * D are identical to the values (69) and (5) computed in Secs. IV.C and E. In contrast to its Model A counterpart, however, we now have within the two-loop approximation Yet the relations (102) imply for the complex non-equilibrium Model B that to all orders in perturbation theory The dynamic critical exponent thus assumes its standard equilibrium value, as for the driven-dissipative extension of Model A. In addition, η ′ = η is confirmed exactly here, as opposed to the situation in Model A, and to Eq. (105), for which we can establish this identity only up to O(ǫ 3 ) corrections. This shows that unlike the driven-dissipative extension of Model A, no independent critical behavior is found for an analogous extension of Model B, at least to two-loop order.

V. CONCLUSION AND OUTLOOK
We have obtained a detailed picture of the drivendissipative Bose condensation transition using the welldeveloped framework of the perturbative field-theoretical dynamic renormalization group, in this way complementing a previous functional renormalization group study 22,23 . In particular, we traced back the existence of a new, independent critical exponent to additional UV divergences without counterpart in the critical theory for the equilibrium Bose condensation transition, and obtained its analytical value at twoloop order in the dimensional ǫ expansion. This exponent is present in the scaling form of the dynamic single-particle response of the system, and in its two-point correlation functions, and witnesses non-equilibrium conditions at the largest distances in the problem. We furthermore confirmed explicitly the asymptotic thermalization scenario in the renormalization group flow for this system, as in Ref. 37 demonstrating the stability of the equilibrium fixed point, along with the construction of the corresponding dynamical scaling forms and calculation of the associated critical exponents.
The perturbative field theoretic approach offers the possibility for a direct comparison of the results for non-equilibrium systems with well-studied counterparts in thermodynamic equilibrium. This makes it a valuable tool for future investigations of critical driven-dissipative quantum systems, working towards the ultimate goal of a systematic classification of non-equilibrium dynamic criticality to a similar maturity level as has been achieved in thermodynamic equilibrium. In the present context, first steps include the exploration of different internal symmetries, such as general O(n) rotation invariance, or the inclusion of explicit coherent pumping processes. They also comprise an investigation of universal aspects of the dynamics following a parameter quench, such as the determination of the initial slip exponent and critical aging 45,46 in driven-dissipative systems. It remains to be seen whether this approach can also be leveraged to situations where criticality and genuine quantum effects come into play simultaneously. (A7) Therefore one obtains in d = 4 − ǫ dimensions ∂ q 2 D d (τ, q) | q=0 = − Γ(1 + ǫ) Noting that Γ(1 + ǫ)/Γ(1 + ǫ/2) 2 = 1 + O(ǫ 2 ), and that the parameter integrals are regular in the limit ǫ → 0, one finally isolates the 1/ǫ pole and its residuum ∂ q 2 D d µ 2 , q sing.
For the integral (73), we proceed similarly: Feynman's parametrization leads to I d (τ, r K , r U ) = which remarkably comes out independent of the parameter r U .