Nonequilibrium Fixed Points of Coupled Ising Models

Driven-dissipative systems are expected to give rise to nonequilibrium phenomena that are absent in their equilibrium counterparts. However, phase transitions in these systems generically exhibit an effectively classical equilibrium behavior in spite of their nonequilibrium origin. In this paper, we show that multicritical points in such systems lead to a rich and genuinely nonequilibrium behavior. Specifically, we investigate a driven-dissipative model of interacting bosons that possesses two distinct phase transitions: one from a high- to a low-density phase—reminiscent of a liquid-gas transition—and another to an antiferromagnetic phase. Each phase transition is described by the Ising universality class characterized by an (emergent or microscopic) ℤ2 symmetry. However, they coalesce at a multicritical point, giving rise to a nonequilibrium model of coupled Ising-like order parameters described by a ℤ2×ℤ2 symmetry. Using a dynamical renormalization-group approach, we show that a pair of nonequilibrium fixed points (NEFPs) emerge that govern the long-distance critical behavior of the system. We elucidate various exotic features of these NEFPs. In particular, we show that a generic continuous scale invariance at criticality is reduced to a discrete scale invariance. This further results in complex-valued critical exponents and spiraling phase boundaries, and it is also accompanied by a complex Liouvillian gap even close to the phase transition. As direct evidence of the nonequilibrium nature of the NEFPs, we show that the fluctuation-dissipation relation is violated at all scales, leading to an effective temperature that becomes “hotter” and “hotter” at longer and longer wavelengths. Finally, we argue that this nonequilibrium behavior can be observed in cavity arrays with cross-Kerr nonlinearities.


I. INTRODUCTION
The increasing control over synthetic quantum systems has provided new avenues into studying many-body physics that are not accessible in conventional condensed matter systems. In particular, driven-dissipative systems, defined by the competition between a coherent drive and dissipation due to the coupling to the environment, have emerged as a versatile setting to investigate nonequilibrium physics [1]. They are very naturally realized by a variety of emerging quantum simulation platforms ranging from exciton-polariton fluids [2][3][4][5][6][7], to trapped ions [8,9], to Rydberg gases [10][11][12][13], to circuit-QED platforms [14,15]. At long times, these systems approach a nonequilibrium steady state due to the interplay of drive and dissipation. The steady states can potentially harbor novel phases and exhibit exotic dynamics. Situated far from equilibrium, understanding the properties of these steady states requires methods beyond those suitable in or near equilibrium. The quest to realize and characterize macroscopic phases of these nonequilibrium systems has sparked a flurry of theoretical and experimental investigations.
Given their nonequilibrium dynamics, driven-dissipative systems are expected to exhibit universal, critical properties different from their equilibrium counterparts. In spite of this, it has become increasingly clear that an effective temperature, and thus an effectively classical equilibrium behavior, emerges in a large class of driven-dissipative phase transitions [16,17]. In particular, the equilibrium Ising universality class and, more generally, the model A dynamics of the Hohenberg-Halperin classification-describing the critical behavior of a nonconserved order parameter in or near equilibrium-have emerged in a variety of drivendissipative phase transitions; these include bosonic or photonic Bose-Hubbard systems [18][19][20][21], various driven-dissipative spin models near an Ising [12,[22][23][24][25][26][27], antiferromagnetic [24][25][26][27][28][29][30][31][32], or limit-cycle phase [25][26][27][28], as well as driven-dissipative condensates consisting of polaritons [33,34]. A possible exception is a two-dimensional driven-dissipative condensate, where it has been argued that the nonequilibrium Kardar-Parisi-Zhang universality class governs the long-wavelength dynamics [35,36]. While existing experiments are consistent with an effective thermal behavior [37,38], the Kardar-Parisi-Zhang dynamics is expected to emerge under different experimental conditions. In general, an important goal of the field is to identify whether generic driven-dissipative systems can escape the pull of an effective equilibrium behavior and instead give rise to new nonequilibrium universality classes. In particular, it has proved difficult to identify nonequilibrium universal behavior, which is genuinely of a quantum nature; see Refs. [39][40][41] for recent proposals to achieve nonequilibrium quantum criticality and Ref. [42] for numerical evidence of an equilibrium quantum critical point in a driven-dissipative system.
An effective equilibrium behavior is not special to driven-dissipative quantum systems. In driven-diffusive classical systems too, where the drive, as well as the dynamics, is entirely classical, effective equilibrium seems to be remarkably robust. For instance, an Ising-type dynamics governing a nonconserved order parameter is argued to be stable against all dynamical, nonequilibrium perturbations [43]. More generally, the universal dynamics of various models in the Hohenberg-Halperin classification [44] are shown to be robust against a class of nonequilibrium perturbations that violate detailed balance [45][46][47][48][49][50][51][52][53][54]; truly nonequilibrium behavior emerges under more constrained dynamics, for example, in the presence of a conserved order parameter in an anisotropic medium [55][56][57][58][59][60][61]. In much of the previous work, situations have been considered where the phase transition is governed by a single order parameter. Because of the restriction that this places on the dynamics, a description based on an effective Hamiltonian often becomes available, hence the emergence of an effective equilibrium behavior.
In this work, we consider a driven-dissipative model that gives rise to multicritical points defined by the joint transition of two order parameters. In particular, we investigate the interplay of two phase transitions, each of which has been studied extensively in drivendissipative settings: One is the many-body analog to optical bistability; in the other, a sublattice symmetry is spontaneously broken, leading to an antiferromagnetic ordering. A schematic illustration of this combination is shown in Fig. 1. With two order parameters, the nonequilibrium dynamics is much less constrained than that of equilibrium, and an immediate identification of an effective Hamiltonian is no longer possible. Remarkably, we show that a new, genuinely nonequilibrium universal behavior emerges at the multicritical point, giving rise to exotic critical behavior and dynamics. Our proposal to observe nonequilibrium critical behavior relies on tuning the system parameters (such as drive and detuning, which are easy to control) to a multicritical point. In fact, the driven-dissipative setting of our model can be experimentally realized using cross-Kerr nonlinearities in cavity arrays [30,31,62].
In order to determine the critical behavior, we will employ the Keldysh-Schwinger and functional integral formalism suited for the nonequilibrium setting of driven-dissipative systems [20,24,[33][34][35][36]39,40,[63][64][65][66][67]. While the presence of two order parameters prevents an immediate Hamiltonian description, the long-wavelength universal behavior-and whether or not the macroscopic behavior escapes an equilibrium fixed point-is determined by investigating how the parameters evolve under a dynamical version of renormalization-group (RG) techniques [68].
The remainder of this paper is organized as follows. In Sec. II, we present the main results of this paper and a summary of the nonequilibrium critical properties emerging in our drivendissipative model. In Sec. III, we discuss the phase diagram of the model and identify the multicritical points where two distinct phase transitions meet. In Sec. IV, we present the RG analysis and show that a pair of new classical nonequilibrium fixed points (NEFPs) emerge that exhibit a variety of novel critical behaviors. In Sec. V, we discuss an experimental setting based on cavity arrays to realize the multicritical points of our model. Finally, in Sec. VI, we conclude our paper with a discussion of possible future directions, which are motivated by the results of our work. In the Appendixes, we present technical details omitted from the main text.

II. MAIN RESULTS
In this section, we present the main results of this paper. We consider a driven-dissipative model that displays two distinct phase transitions, each of which arises generically in various settings. The first one is a many-body version of bistability where two stable solutions arise with a low or high population of photons (or excitation of spins). In the thermodynamic limit, the bistable region is reduced to a line of first-order phase transitions that terminates at a critical point, reminiscent of a liquid-gas phase transition. The second type of phase transition is one to an antiferromagnetic phase where the population takes different values on the two sublattices (a or b) of the system. We consider a model where these phase transitions coalesce at a multicritical point and investigate the exotic dynamics that arise due to the interplay of the respective order parameters. These features are provided, for example, in a driven-dissipative model of weakly interacting bosons with nearest-neighbor density-density interactions on a d-dimensional cubic lattice. The coherent dynamics of the model is where Δ = ω D − ω 0 is the detuning of the drive (ω 0 is the frequency of the bosons and ω D the drive frequency), Ω the drive strength, J the hopping strength, and V the strength of the nearest-neighbor interactions. The incoherent dynamics is due to loss of bosons, characterized by the Lindblad operators L i = Γa i , where Γ defines the loss rate. The (mixed) state of the system ρ evolves under the quantum master equation until it approaches a nonequilibrium steady state at long times where ρ = 0. The interplay of the coherent drive (the linear term in the Hamiltonian) and dissipation, together with the interactions, tends to give rise to bistability, while the nearest-neighbor interactions can lead to an antiferromagnetic phase. We stress that our general results should hold beyond the specific model considered here; for example, the addition of on-site interactions or densitydependent hopping terms to our model also gives rise to multicritical points whose universal properties should be independent of the microscopic model considered. More generally, the relevant features of our bosonic model also arise in a variety of driven-dissipative systems including spin models [20,23,24,26,28,30,31]. We have chosen this particular model as a minimal driven-dissipative setting that gives rise to multicritical points, although our general conclusions should apply to a large class of models.
Each phase transition in our model is characterized by an Ising-like order parameter (low or high density in bistability and sublattice a or b in the antiferromagnetic transition). The simple structure of the order parameter, together with the incoherent nature of the dynamics, puts a strong constraint on the universal properties of the phase transition. Thus, it may be expected that each phase transition alone is described by the Ising universality class that also governs the Ising-type transitions in equilibrium. It can be argued, on more formal grounds, that this is indeed the case. Associating the order parameter ϕ 1 with bistability and ϕ 2 with antiferromagnetic ordering, their long-wavelength properties in the steady state are governed by a thermal distribution but with respect to the effective (classical) Hamiltonians, with D i characterizing the stiffness, g i the interaction strength, r i the distance from the critical point (which will shift once fluctuations are taken into account), and h an effective magnetic field. Note that due to sublattice symmetry, there is no magnetic field associated with the antiferromagnetic phase. Furthermore, the incoherent nature of the model leads to stochastic Langevin-type dynamics of the order parameters as [68] where ζ i is a "friction" coefficient and ξ i describes Gaussian white noise with correlations where T i is the effective temperature of the system. Near equilibrium, the "friction" coefficients ζ i control the rate at which the system relaxes to a thermal state via dissipating energy and thus is a purely dynamical quantity. The noise itself is related to the dissipation (i.e., friction) through temperature in what is known as the Einstein relation, which itself is a consequence of the fluctuation-dissipation theorem [68]. In the nonequilibrium context of driven-dissipative models, where there is no intrinsic temperature, the ratio of the noise level to the dissipation can be used to define an effective temperature at long wavelengths. Even a nonequilibrium system that is effectively (i.e., at large scales) governed by the Hamiltonian dynamics [as in Eq. (4)] is effectively in thermal equilibrium. This condition is often satisfied for a single Ising-like order parameter [58], although notable examples exist where this is not the case [55][56][57][58][59][60][61]. Notice that, with the appropriate scaling of the fields ϕ i , the effective temperatures can be set to T i = 1. Therefore, as long as the two order parameters are not coupled, their distribution in the steady state is given by e −ℋ 1 − ℋ 2 .
The situation is entirely different in the vicinity of multicritical points where the two order parameters are generally coupled. Given the underlying symmetries, the dynamics can always be brought to the form ζ 2 ∂ t ϕ 2 = − δℋ 2 δϕ 2 − g 21 ϕ 2 ϕ 1 2 + ξ 2 . (6b) Notice that the new terms that couple the two fields respect the underlying Ising symmetry of both order parameters. The noise correlations are given by Eq. (5); we again exploit our freedom in scaling the fields to set T 1 = T 2 = 1. With the two order parameters coupled, the condition for an effective equilibrium description becomes much more restrictive. A thermal description requires the entire dynamics to be described by a single Hamiltonian, which only occurs when g 12 = g 21 , leading to the effective Hamiltonian ℋ = ℋ 1 + ℋ 2 + g 12 2 ∫ x ϕ 1 2 ϕ 2 2 , (7) in which case the steady-state distribution is given by e −ℋ . However, this will generally not be the case, so we must consider how various parameters flow under RG. While the microscopic (though coarse-grained) dynamics is not immediately described by a thermal state, it could very well be the case that the RG flow pulls the system into a thermal fixed point where g 12 * = g 21 * at macroscopic scales. Indeed, we show that this is the case roughly when the microscopic values of the coupling constants are both positive, i.e., when g 12 , g 21 > 0. It is rather remarkable that equilibrium restores itself under RG, showcasing another instance in which equilibrium proves to be a robust fixed point even when the system is driven far from equilibrium. However, this is not the end of the story: We show that, when the microscopic couplings have opposite sign (g 12 g 21 < 0), a pair of two NEFPs emerge where g 12 * = − g 21 * . (Both equilibrium fixed points and NEFPs are shown in Fig. 2.) Furthermore, we argue that for the model under consideration, the critical behavior will be governed by one of the NEFPs. These fixed points give rise to a new nonequilibrium universality class exhibiting a variety of exotic features that generically do not, or even cannot, arise in any equilibrium setting. A summary of the most interesting features, including the critical behavior, of the new NEFPs is illustrated in Fig. 3. In the following subsections, we discuss these features in detail.

A. Scaling phenomena
In the vicinity of a RG fixed point governing a phase transition, the system exhibits universal scaling behavior characterized by critical exponents, regardless of the microscopic model. The scaling behavior of the correlation and response functions at a NEFP or in its vicinity takes, respectively, the form where ℱ denotes the Fourier transform in both space (x) and time (t), with q the momentum and ω the frequency, the curly brackets denote the anticommutator, and the subscript c indicates the connected part of the correlation function. Here, r = r 1 2 + r 2 2 is the distance from the multicritical point, P is a 2π-periodic function, and the functions C and χ are scaling functions. While, in principle, the scaling behavior could be different for the two order parameters (ϕ 1 and ϕ 2 ), we argue, based on a systematic RG analysis, that a stronger notion of scaling emerges where the critical (static and dynamic) behavior and exponents characterizing the two order parameters become identical. Therefore, we can express either the correlation or the response function via a single scaling function (and not one for each order parameter) with the same set of exponents. The exponents η and η′ define the anomalous dimensions corresponding to correlation and response functions, respectively, and z is the dynamical critical exponent characterizing the relative scaling of time with respect to spatial coordinates. The correlation length ξ is described by the critical exponent ν′ via ξ ∝ r −ν′ . Typically, it is the exponent ν, associated with the scaling behavior of r 1 and r 2 , that describes the scaling of the correlation length. However, the latter exponent becomes complex valued at the NEFPs, ν −1 = ν′ −1 + iν″ −1 , with the real part determining the scaling of the correlation length and the imaginary part leading to a discrete scale invariance, as we shall discuss shortly. Altogether, there are five independent critical exponents of interest: ν′, ν″, η, η′, z.
Critical points are generically associated with a continuous scale invariance where the system becomes self-similar upon an arbitrary rescaling of the momentum and frequency. However, because of the "log-periodic" function in the scaling functions, the correlations are self-similar upon the rescaling q → b*q and ω b * z ω for a particular scaling factor b * = e 2πν″ , (9) or any integer powers thereof. Rather than a typical continuous scale invariance, this behavior is indicative of a discrete scale invariance, which is reminiscent of fractals, shapes that are self-similar under particular choices of scaling [69]. A schematic depiction of the correlation functions with discrete scale invariance is shown in Fig. 3(a). Additionally, since the origin of the discrete scale invariance is the scaling behavior of r i that characterizes the distance from the critical point, the phase boundaries themselves also exhibit a form of discrete scale invariance in r i ; see Fig. 4 and the discussion in the next subsection.
The critical exponents η and η′ characterize the anomalous dimensions corresponding to fluctuations and dissipation, respectively. In an equilibrium setting, the fluctuationdissipation theorem dictates that the correlation and response functions are related as [68] C(q, ω) = 2T ω Imχ(q, ω) . (10) (We have assumed the classical limit of the fluctuation-dissipation theorem at low frequencies and at a finite temperature.) In an equilibrium setting, the temperature T is just a constant set by an external bath and thus is scale invariant. Therefore, the overall scaling behavior of the correlation and response functions is identical apart from the dynamical scaling (due to ω −1 on the rhs of the above equation) set by the critical exponent z. This in turn puts a constraint on the critical exponents as η = η′ for effectively equilibrium phase transitions. However, we find η ≠ η′ at the NEFPs, indicating the violation of the fluctuation-dissipation theorem and resulting in a new exponent characterizing the nonequilibrium nature of the fixed point. This in turn results in an effective temperature that remains scale dependent at all scales. Inspired by the fluctuation-dissipation theorem, we define an effective "temperature" as Env[C(q, ω)] = 2T eff (q, ω) ω Env[Im[χ(q, ω)]] . (11) To factor out the log-periodic nature of the correlation and response functions, we have made a convenient choice by postulating a fluctuation-dissipation relation between the envelope (Env) functions of the correlation and the response functions. This relation can be defined via either the upper or lower envelope functions. We can then identify the scaling behavior of the effective temperature at the NEFP, which we find to be T eff ∼ | q| η − η′ at long wavelengths and fixed ω/|q|z. Interestingly, we find that η′ > η, so the system gets "hotter" and "hotter" at longer and longer scales. Of course, the divergence of the effective temperature at long wavelengths does not imply an infinitely energetic state; rather, it reflects the fact that, at longer wavelengths, the correlation function is increasingly larger compared to the response function than one would expect in an equilibrium setting based on the fluctuation-dissipation theorem. This behavior is illustrated in Fig. 3(b) individually for the two effective temperatures corresponding to the two order parameters. At long wavelengths, these effective temperatures become identical to each other and to T eff (q, ω) defined above. Finally, the values of the critical exponents at the NEFPs are provided to the lowest nontrivial order in ϵ = 4 − d in Fig. 3(d).

B. Phase diagram
The critical point described by the new fixed points is a tetracritical point. In the vicinity of the tetracritical point (with h = 0), there are four different phases where none, one, or both of the order parameters undergo a continuous phase transition. A particularly exotic feature of the phase diagram is that it exhibits spiraling phase boundaries. This feature leads to the discrete scale invariance of the phase diagram itself, a property that follows from the same feature of the scaling functions in Eq. (8). In contrast, depending on the microscopic model, the equilibrium fixed points can give rise to either a bicritical point-in which case there will not be a phase where both order parameters undergo a continuous phase transition-or a tetracritical point; neither of these will exhibit spiraling phase boundaries. Note that since the ℤ 2 symmetry associated with the bistability transition (ϕ 1 → −ϕ 1 ) is only an emergent one (when h = 0), the full phase diagram (including h ≠ 0) can be better described as a threedimensional plot that also includes the first-order phase transitions characteristic of bistability; see Fig. 7. The contrast between the equilibrium fixed points and NEFPs can further provide a route to experimentally identify the new fixed points. An overview of the properties of bicritical and tetracritical points in equilibrium systems can be found in Refs. [70][71][72][73][74][75][76].

C. Spectral properties
The NEFP can be further distinguished by its particular dynamics that governs the relaxation of the system to the steady state. In the nonequilibrium setting of our model, the dynamics is described by the Liouvillian ℒ via [cf. Eq. (2)] rather than a Hamiltonian. However, in analogy with the ground state that is described by the smallest eigenvalue of the Hamiltonian, the steady state(s) is given by the 0 eigenvalue(s) of the Liouvillian; all the other eigenvalues of the Liouvillian have a negative real part characterizing the decay into the steady state. Furthermore, the spectral gap of the Hamiltonian is naturally generalized to the eigenvalue of the Liouvillian with the smallest (in magnitude) nonzero real part. We denote this eigenvalue by Γ ℒ . For a continuous phase transition, just like the spectral gap, the closing of the Liouvillian gap results in the divergence of a timescale associated with a slow or soft mode of the dynamics. The fashion that the latter gap closes reveals characteristic information about the nature of the phase transition. In equilibrium phase transitions at finite temperature, this gap becomes real (i.e., purely dissipative) as the critical point is approached. Even when the microscopic dynamics is far from equilibrium, the Liouvillian gap may (and typically does) become real, leading to effectively thermal equilibrium. In contrast, the dynamics near the NEFPs can close away from the real axis. This indeed occurs in the doubly ordered phase; let M i = ⟨ϕ i ⟩ ≠ 0 define the nonzero order parameters, and redefine the fields as ϕ i → ϕ i + M i . We then find the linearized equations of motion as where, at the NEFPs, g 12 * = − g 21 * and g 1 * = g 2 * , while we can choose ζ 1 * /ζ 2 * = 1; noise, gradient, and higher-order terms have been dropped. Because of the opposite signs of g 12 and g 21 in the two equations, we find a spiral relaxation to the steady state. This relaxation in turn is characterized by a complex Liouvillian gap-defined by a conjugate pair of complex eigenvalues-which exhibits both a dissipative (real) and a "coherent" (imaginary) part depending on the values of M 1 and M 2 . We find that when |M 1 | = |M 2 |, the angle of this complex gap relative to the real line achieves its maximum value of π/3. This behavior is illustrated in Fig. 3(c). The corresponding mean-field relaxational dynamics is illustrated in Fig. 5.

III. MODEL
The representative model we have focused on is a driven-dissipative system of weakly interacting bosons defined in Eqs. (1) and (2). In order to understand how this model gives rise to bistability and antiferromagnetic ordering, we begin this section with a detailed discussion of mean-field theory and corrections, or fluctuations, on top of the mean-field solutions. Along the way, we identify the soft modes of the dynamics that ultimately describe the critical behavior of the multicritical point. Finally, we conclude this section by presenting a mapping of our nonequilibrium model to a model of coupled Ising-like order parameters with a ℤ 2 × ℤ 2 symmetry, corresponding to the sublattice symmetry as well as the emergent Ising symmetry due to bistability.

A. Mean-field theory
In order to analyze the phase diagram of our model, we begin with a mean-field analysis, in which we assume different sites are uncorrelated; that is, for any two operators A i and B j on neighboring sites, we have ⟨A i B j ⟩ = ⟨A i ⟩ ⟨B j ⟩ [77,78]. Additionally, we assume that individual sites are described by coherent states. While the latter assumption follows from the former in our model, this will generally not be the case, for example, when on-site Hubbard interactions are present. However, a systematic path-integral formalism (adopted in subsequent sections) beyond mean-field theory is perfectly suited to analyzing the latter type of interaction. Finally, in anticipation of the antiferromagnetic phase transition, we separate the system into two sublattices, a and b, and assume each to be described by a single coherent state.
Following these assumptions and using the fact that ∂ t O = Tr(ρO) for an arbitrary operator O, the resulting mean-field equations of motion are given by where ψ i corresponds to the coherent state ⟨a⟩ on sublattice i ∈ {a; b} and z is the coordination number; from here on, we absorb z in the microscopic parameters via z J J and z V V . It is clear from these equations that the density-density interaction behaves as an effective detuning that depends on the density of the other sublattice. This behavior results in physics that is similar to Rydberg excitations in stationary atoms, in which the presence or absence of a Rydberg excitation on one site can either prevent (blockade) [79][80][81] or facilitate (antiblockade) [82,83] a Rydberg excitation on a neighboring site by shifting it away or towards the effective resonance.
Setting ψ a = ψ b , one can immediately see that the mean-field equations become identical to those describing bistability; cf. Ref. [20] where the nonlinearity due to the Hubbard interaction should be replaced by the density-density interactions in this context. The emergence of bistability can be understood in simple terms: Away from resonance, there is a low population on each site. However, once a sufficient number of sites are highly excited, they begin to facilitate the excitation of neighboring sites, resulting in a high-population steady state. This process occurs when the shift in detuning due to interactions is comparable to the detuning. This condition is satisfied approximately when [Ω2/[Γ2/4 + (Δ + J)2)]V ≈ Δ + J, where J behaves like an effective detuning while the product of the interaction strength V and the noninteracting steady-state population [Ω2/[Γ2/4 + (Δ + J)2)] gives the interaction-induced shift of the detuning. For Γ ≳ Δ + J, this reasoning becomes blurred as the drive is effectively always on resonance due to the larger linewidth. As a result, a finite region of bistability emerges with low-and high-population steady states. Beyond meanfield theory, the bistable region is replaced by a line of first-order phase transitions that terminates at a critical point.
The presence of antiferromagnetic ordering in this system can be understood by inspecting the role of the density-density interactions. Since the interaction affects neighboring sites only, the blockade effects occur between sublattices but not within each sublattice. For example, if one sublattice has a high population, it can prevent further excitations in the other sublattice. Similar to the case of bistability, the phase boundary occurs approximately when the shift in detuning due to interactions takes the system out of resonance. This approximately occurs when |[Ω2/[Γ2/4 + (Δ + J)2)]V − Δ -J| ≳ Γ, i.e., when one sublattice is effectively more than a linewidth out of resonance due to interactions. Unlike bistability, this process does not break down as Γ and Ω are increased. As the decay Γ is increased, the drive strength Ω can be further increased so that the interaction-induced shift in the detuning compensates for the increase of the linewidth.
In order to better understand the mean-field structure of the model, it is convenient to introduce a new set of fields corresponding to the two order parameters as The field ψ B captures the effects of bistability while ψ AF describes the antiferromagnetic ordering. The mean-field equations can in turn be cast in terms of these fields as as well as Ψ c = 2/3V e −iπ/3 as the steady-state value of ψ B at the critical point (by virtue of symmetry, ψ AF = 0 The two fields ψ B/AF are complex valued, thus comprising four real (scalar) fields. However, given the Ising nature of each transition, we must anticipate that two scalar fields would be sufficient to describe the critical behavior of both types of ordering. Indeed, we find that, at the multicritical point, two massless fields emerge-defined by appropriate components of the original fields-corresponding to the soft (or slow) modes ϕ i , while the other components ϕ i ′ remain massive and are therefore noncritical (or fast). We then adiabatically eliminate the two noncritical modes by setting φ i ′ = 0 and solving for ϕ i ′ in terms of ϕ i . Upon substituting our solutions for the massive fields into φ i , we find an effective description in terms of the soft modes. We closely follow Refs. [20,24] to identify these modes. For the bistability order parameter, we can identify ψ B = Ψ c + e iπ/3 ϕ 1 + ϕ 1 ′, (18) with the real fields ϕ 1 and ϕ 1 ′ characterizing the slow and fast modes, respectively. A similar identification has been made in Refs. [20,24]; see also Refs. [84,85] for a similar reasoning, although the slow and fast modes identified there make an angle of π=2. For the antiferromagnetic field, the massless and massive components depend on the choice of the multicritical point in Eq. (17) as Again, the unprimed fields are massless while the primed fields are massive. The slow and fast modes of the fields are illustrated pictorially in Fig. 6.
Next, we adiabatically eliminate the massive modes to find an effective description in terms of the soft modes. Including the gradient terms-describing the coupling between neighboring sites-as well as the noise terms due to the coupling to the environment, we find the Langevin equations ζ 1 φ 1 = ℎ − r 1 ϕ 1 + D 1 ∇ 2 ϕ 1 + ξ 1 + A 20 ϕ 1 2 + A 02 ϕ 2 2 + A 12 ϕ 1 ϕ 2 2 + A 30 ϕ 1 3 , (20a) with Gaussian noise Higher-order terms that are irrelevant in the sense of RG have been neglected. We have expressed the noise coefficients in a convenient notation that mimics the dissipative dynamics in thermal equilibrium, in spite of the underlying nonequilibrium dynamics. Finally, the details of the adiabatic elimination together with the explicit values of all the coefficients (h, r, D, A, B, ζ, and T) in terms of microscopic parameters of the model are provided in Appendix A.
It turns out that, at the level of mean-field analysis, A 20 = 0 in the vicinity of the multicritical point. The resulting mean-field dynamics (neglecting the gradient and noise terms) of the two soft modes is then described by a cusp-Hopf bifurcation; a detailed analysis of this type of bifurcation can be found in Ref. [86]. However, since A 20 is not protected by any symmetries, the corresponding term can be generated in the course of RG and become of the order of the other quadratic terms. While we focus on the multicritical points, further details about the full mean-field phase diagram of our model and slight variations on it can be found in Refs. [30,31].

B. Nonequilibrium Ising model for two fields
Before proceeding with our perturbative RG analysis, it is important to identify what are known as redundant operators. These are terms in the action that are generated under suitably local symmetry-preserving transformations of the fields. Since such a transformation should not change the long-distance behavior of the system, this redundancy can be used to simplify our analysis. This is an important step for perturbative RG and to identify the upper critical dimension; see Ref. [87] for a discussion of redundant operators in an equilibrium setting.
As a simple illustrative example, consider the generic Hamiltonian of a scalar field ϕ in the absence of a ℤ 2 symmetry: Shifting the field by a constant as ϕ → ϕ + ϕ 0 , the Hamiltonian is given by the same expression (up to an unimportant additive constant) with possibly different coefficients. This underscores a redundancy in Hamiltonians that describe the same physical system. The change of the Hamiltonian Δℋ (or, rather, the integrand) due to a constant shift of the field defines a redundant operator. In particular, the cubic term transforms as By choosing the value of ϕ 0 properly, the ϕ 3 term can be dropped from the Hamiltonian while shifting the coefficients of the terms ϕ and ϕ 2 [88]. Similar to the above example, we should first identify the redundant operators in the nonequilibrium setting of the two coupled scalar fields ϕ 1 and ϕ 2 . In this case, we allow for a more general, nonlinear transformation that is suitably local and retains the underlying symmetries. We find that the set of redundant operators in our model is sufficient to remove all the quadratic terms in the Langevin equation (or, equivalently, the cubic terms in the action, similar to the Hamiltonian in the above example); the details of this analysis are presented in Appendix B. In particular, we find that, under this transformation, the new value for the ratio A 12 /B 21 is given by 2A 02 /B 11 ; therefore, the relative sign of the quadratic terms (prior to the transformation) determines the relative sign of the cubic terms in the final equations of motion. In the model that we have considered here, the two quadratic terms have opposite signs (see Appendix A). This fact will be important in determining the fixed point of the RG flow. In fact, we show that the above sign difference leads the system to one of the NEFPs.
With the above considerations, the Langevin equations can finally be brought into a canonical form as with Gaussian noise Therefore, the dynamics exhibits a ℤ 2 × ℤ 2 symmetry when h = 0, corresponding to the emergent symmetry ϕ 1 → −ϕ 1 in addition to the sublattice symmetry ϕ 2 → −ϕ 2 . Such emergent symmetry has previously been identified in the bistability transition [20,23,24]; our analysis shows that such symmetry emerges even in the vicinity of a multicritical point where bistability and antiferromagnetic transitions coalesce. We must point out that, even in the absence of the latter symmetry, the sublattice symmetry alone prevents any mixing of the gradient and mass terms between the two fields as well as the noise terms, a property that should hold to all orders of perturbation theory.

IV. RENORMALIZATION GROUP ANALYSIS
In this section, we derive the perturbative RG equations to the two-loop order (for reasons that will be explained shortly), identify the fixed points, and characterize the critical exponents that determine the scaling properties of correlations near the multicritical point.

A. RG equations
The Langevin-type equations can be cast in terms of the response-function formalism. This method allows us to study our nonequilibrium model by extending the standard techniques of the RG analysis to a dynamical setting; see, e.g., Ref. [68] for more details. The nonequilibrium partition function is defined by where the functional integral measure, as well as the "action" A, involves both fields ϕ i , with i = 1, 2, and their corresponding "response" fields ϕ i . In the language of Keldysh field theory, ϕ corresponds to the classical field, while ϕ/2i corresponds to the quantum field. The statistical weight of ϕ i (t; x) can be obtained by integrating out both response fields as While the partition function Z = 1 by construction, the expectation value of any quantity-the fields themselves or their correlations-can be determined by computing a weighted average in the partition function. For our model defined by Eqs. (24) and (25), we write the action as the sum of quadratic and nonlinear (beyond quadratic) terms with the quadratic action given by and the nonlinear interaction terms A int [ϕ i , ϕ i ] = ∫ t, x g 1 ϕ 1 3 ϕ 1 + g 2 ϕ 2 3 ϕ 2 + g 12 ϕ 1 ϕ 2 2 ϕ 1 + g 21 ϕ 2 ϕ 1 2 ϕ 2 . (26c) Our goal is to determine the RG flow of various parameters in the action and, specifically, of the coefficients g of the interaction terms.
We begin by considering the subspace defined by g 12 g 21 = 0 when either g 12 = 0 or g 21 = 0.
This subspace is special in that it is closed under renormalization to all orders. The reason is that when g 12 = 0 or g 21 = 0, one of the two fields is not affected by the other at the microscopic level, a property that should hold at all scales. This result can also be understood perturbatively in a diagrammatic scheme: If, say, g 12 = 0, then all diagrams that could generate g 12 involve a causality violation; hence, it should remain zero to all orders. An important consequence of this fact is that the relative sign of g 12 and g 21 cannot change, as this would require passing through the closed subspace.
Before performing the RG analysis, we first use our freedom in rescaling the fields to cast the action in a more convenient form. In Sec. II, we used this freedom to set both temperatures to unity; here, for the convenience of the RG analysis, we make a different choice. Note that rescaling ϕ 2 → cϕ 2 and ϕ 2 ϕ 2 /c maps g 2 → c 2 g 2 , g 12 → c 2 g 12 , and T 2 → T 2 /c 2 . Exploiting this freedom, we can set the rescaled value of g 12 to be identical to g 21 up to a sign. In doing so, we have effectively shifted the renormalization of g 12 /g 21 onto T 1 / T 2 , simplifying the RG analysis later. Note, however, since g 12 is rescaled by a factor c 2 , this transformation cannot change the relative sign of g 12 and g 21 . This behavior is indeed consistent with the closure of the g 12 g 21 = 0 subspace discussed above. Having rescaled the fields appropriately, we can write the action as (the quadratic action is repeated for completeness) A int [ϕ i , ϕ i ] = ∫ t, x u 1 ϕ 1 3 ϕ 1 + u 2 ϕ 2 3 ϕ 2 + u 12 ϕ 1 ϕ 2 (ϕ 2 ϕ 1 + σϕ 1 ϕ 2 ), where σ = ±1 indicates the relative sign of g 12 and g 21 and the coefficients u 1 , u 2 , and u 12 define the rescaled values of the interaction strengths (in an abuse of notation, we use the same notation for the other rescaled parameters of the model as well as the rescaled fields).
Let us first briefly consider σ = 1, in which case the action can be written in a suggestive form as where the function ℋ is given by Put in this form, Eq. (28) bears close resemblance to an equilibrium setting where the dynamics is governed by a Hamiltonian (in this case, ℋ). However, with each field at a different temperature, their coupled dynamics does not generally satisfy fluctuationdissipation relations, and thus an (effective) equilibrium behavior cannot be established, at least at the microscopic level. (Note that unlike Sec. II, we have already used the scaling freedom in redefining the interaction parameters, which in turn fixes the ratio T 1 /T 2 .) One then should resort to a RG analysis to determine whether or not effective equilibrium is restored at long wavelengths, that is, if T 1 /T 2 → 1 under RG. We shall see shortly that equilibrium proves to be a robust fixed point even when T 1 /T 2 ≠ 1 at the microscopic level.
In contrast, a Hamiltonian dynamics [similar to Eqs. (28) and (29)] is not possible when σ = −1 since a term proportional to ϕ 1 2 ϕ 2 2 in the Hamiltonian leads to equations of motion that couple the two fields with the same coefficient and hence the same sign. Therefore, in this case, the dynamics cannot flow to an equilibrium fixed point even when T 1 = T 2 , with the exception of a decoupled fixed point where u 12 = 0 (or g 12 = g 21 = 0). Indeed, we shall argue that a pair of genuinely nonequilibrium fixed points emerge in this case.
At a technical level, a RG analysis would be complicated as we need to consider diagrams up to two loops because, at one loop, no renormalization occurs for the temperatures (due to causality) as well as the diffusion constants and friction terms (owing to their momentum and frequency dependence). In contrast, the interaction terms (u 1 , u 2 , and u 12 ) are all renormalized already at one loop. This observation-besides aesthetic reasons-has motivated the representation adopted here; in the original description in terms of g 12 and g 21 , the ratio g 12 =g 21 would not be renormalized at one loop.
To perform the RG analysis, we first define the renormalized parameters as where is the Euler's Gamma function, μ is an arbitrary small momentum scale (compared to the lattice spacing), and ϵ = 4 − d defines the small parameter of the epsilon expansion. The effect of renormalization is captured in the Z factors that contain the divergences according to the minimal subtraction procedure. We determine these factors perturbatively to the lowest nontrivial order in ϵ or loops (the details are provided in Appendix C). The lowest-order corrections to Z r and Z u occur at one loop (~ϵ), while those of Z ζ , Z T , Z D appear at two loops (~ϵ2). These perturbative corrections, while having some similarities with their equilibrium counterparts, are more complicated due to their nonequilibrium nature.
Using the above Z factors, we determine the RG flow and beta functions via where p ∈ {r i ; ζ i , D i , T i } and u a ∈ {u 1 , u 2 , u 12 }. These functions describe the flow of various parameters in the action under the change of the momentum scale μ. In particular, the beta functions identify the fixed points of the interaction coefficients via β u a = 0. At any such fixed point, the scaling behavior of the remaining parameters is governed by power laws whose exponents depend on γ p . Here, we report the beta functions for the interaction parameters u a (the details are provided in Appendix C): where we have introduced D i R ≡ D i R /ζ i R . These equations exhibit a number of important features. First, for u 12 R = 0, we can absorb a factor of T i R /D i R 2 into u i R , leaving the two beta functions for u i independent of T i , D i , ζ i . We thus immediately recover a pair of uncoupled equilibrium Ising phase transitions, as one should expect. Second, under equilibrium conditions where σ = 1 and T 1 R = T 2 R ≡ T R , we recover the standard beta functions in equilibrium. In a similar fashion, we can absorb the factors of T R /D i R 2 into u i R and T R / D 1 R D 2 R into u 12 R , again leaving the beta functions dependent only on the coupling coefficients. This observation underscores the important fact that, in equilibrium, static properties are entirely decoupled from the dynamics. On the other hand, in the setting of our nonequilibrium model, statics and dynamics are inherently intertwined. Indeed, no redefinition of the coupling terms can lead to beta functions that would be independent of T i R and D i R . This is not the case for ζ i R as they can always be absorbed in other parameters; for example, we can still absorb 1/ζ i R 2 into u i R and 1/ ζ 1 R ζ 2 R into u 12 R in the beta functions.
This reflects the fact that, through an appropriate rescaling of the fields, one can always rescale ζ i arbitrarily without changing T i , D i , or the overall structure of the action.
To set up the full RG equations, let us define the parameters With these definitions, the beta functions for the new interaction parameters u a depend only on the renormalized parameters v R and w R . To obtain the full RG equations, we further need to determine the RG evolution of the latter parameters. As we shall see, their RG equations are also closed in the (five) parameters defined in Eq. (33). To see why, first notice that there are ten marginal parameters in the original action at the upper critical dimension (ζ i , D i , T i , g i , g 12/21 ), which can define the basin of attraction for the RG flow. Since all four fields and time can be rescaled relative to an overall momentum scale, a total of five parameters are needed to define the fixed point. The remaining parameters (r i , h) define relevant directions of the RG flow and thus must be tuned to their critical values. In order to determine the RG equations for the parameters v and w, we use the identity We now report the full set of beta functions of the parameters of our model (with r i and h set to zero at the fixed point), where we have defined C = 9 log(4/3) -3/2 and the functions The functions F, G, H always appear in the RG equations in pairs, with one taking w R and the other w R −1 as an argument. This is because the diagrams that contribute to the beta functions come in pairs, corresponding to one from the renormalization of the terms involving ϕ 1 only and the other from those that involve ϕ 2 only. Similarly, under the mapping u 1 R u 2 R , u 12 R σv R u 12 R , v R v R −1 , and w R w R −1 , the beta functions are left unchanged. This invariance reflects the fact that we can switch the role of ϕ 1 and ϕ 2 without changing the physics. As a result, if either σ = −1, v R ≠ 1, or w R ≠ 1 at a given fixed point, there will always be a second fixed point paired with it.
The above equations determine the full RG equations of our nonequilibrium model, but it is instructive to first consider the RG equations under equilibrium conditions where the temperatures are equal, i.e., v R = 1 and σ = 1.We then immediately find that the temperature ratio does not flow, β v = 0; hence, the two temperatures remain identical at all scales.
Furthermore, the temperature itself-and not just the ratio-remains scale invariant (γ T = 0), indicating (effective) thermal equilibrium. Finally, as remarked earlier, the RG equations for the interaction terms become independent of w R under equilibrium conditions, highlighting once again the fact that, in equilibrium, the statics is decoupled from the dynamics.
There are two distinct scenarios with respect to the beta function β w . The first scenario is that the beta function vanishes when γ D 1 = γ D 2 . Since the dynamical critical exponents are related to the flow of D as z i = 2 + γ D i , we find that z 1 = z 2 under this scenario. This result means that both fields are governed by the same dynamical critical exponent, giving rise to a "strong dynamic scaling." The second scenario occurs when γ D 1 ≠ γ D 2 , which would lead to the fixed point w R = 0 or w R = ∞ depending on the sign of γ D 1 − γ D 2 . This behavior is then described by a "weak dynamic scaling" where the two fields exhibit different dynamical scaling properties and exponents [89][90][91]; see also Ref. [68]. Similarly, one can consider the beta function β v characterizing the RG flow of the ratio of the temperatures. In this case, too, there are two scenarios: Either the beta function vanishes for a fixed temperature ratio, or rather, depending on the sign of γ T 1 − γ T 2 , the RG flow leads to either v R = 0 or v R = ∞, which both correspond to the g 12 g 21 = 0 subspace. However, this subspace does not appear to be amenable to perturbative RG. In this sector, we find that w flows to either 0 or ∞, indicating weak dynamic scaling where the two fields are governed by distinct dynamical universality classes. However, in both cases, the fixed-point values of the coupling terms diverge, resulting in a nonperturbative regime that is not accessible within the perturbative RG analysis. This result indicates that an alternative approach from our present analysis should be considered in this scenario. In this work, we restrict ourselves to the case where v R and w R are both finite and nonzero.

B. Fixed points of RG flow
With the RG beta functions, we can now identify the resultant fixed points. In the σ = 1 sector, the only fixed points of the RG equations are those where w R * = 1, exhibiting a strong dynamic scaling, as well as v R * = 1, indicating that the two temperatures become identical at the fixed point. Indeed, aside from the case of u 12 R = 0, the only possible fixed point value of v R at this order is 1. This result can be seen by noting that the only other root of Eq. (35d) is −F w R −1 /F w R , which is always negative and thus unphysical. Similarly, noting that the beta functions for u 1 R , u 2 R are identical at this order, all coupled fixed points in this sector will satisfy u 1 R = u 2 R . In light of this, we immediately identify w R = 1 as the only possible solution of Eq. (35e). Remarkably, an effective equilibrium behavior emerges in this sector despite the underlying nonequilibrium nature of the dynamics. In particular, we recover the familiar equilibrium O(2) and biconical fixed points as well as various decoupled fixed points involving combinations of Gaussian and Ising fixed points. However, there are no additional NEFPs in this sector (possibly with the exception of a kind of weak dynamical scaling in the g 12 g 21 = 0 subspace). Note that the emergent equilibrium is not achieved by a simple rescaling of the terms in the action to mimic an effective Hamiltonian, but it is truly the result of a nontrivial two-loop RG analysis.
In the σ = −1 sector, any nontrivial fixed point is truly nonequilibrium as it cannot be described by effective Hamiltonian dynamics that defines equilibrium. Therefore, we should first determine if there exists any nontrivial fixed point in this sector or, alternatively, if the RG evolution flows to a trivial (decoupled) fixed point. Interestingly enough, the former is the case; we find a pair of genuinely nonequilibrium fixed points as v R * = 1, w R * = 1, These fixed points also exhibit a strong dynamic scaling since w R * = 1, so the two fields are governed by the same dynamical scaling. Furthermore, we find v R * = 1, implying that the two temperatures are equal, which might suggest an equilibrium behavior; however, the latter temperatures only characterize the strength of the noise (more precisely, ζ i T i defines the noise), while a true equilibrium description (and a genuine notion of temperature) requires Hamiltonian dynamics [similar to Eq. (28)], which is inherently impossible in this sector.
While we have identified a new pair of NEFPs, this does not guarantee that they would govern the critical behavior near the multicritical point. If these fixed points are unstable under RG, further fine-tuning would be necessary to access them. Even if they are stable, depending on the initial microscopic parameters, the system could still flow to an equilibrium fixed point under renormalization. Nevertheless, we argue that the multicritical point is indeed governed by the new NEFPs.
To determine the stability of the fixed points, we need to consider the stability matrix Λ ab = ∂β a ∂s bR , (38) where s b denotes the set of parameters that enter the RG beta functions. A fixed point is stable if all of the eigenvalues of Λ are positive. Although we have determined the lowestorder corrections to all five parameters, we can only determine these eigenvalues up to O(ϵ).
This is because in order to fully determine Λ to O ϵ 2 , we need to consider the two-loop corrections to the coupling terms u and the three-loop corrections to v, w. However, when a fixed point possesses a higher symmetry than the underlying field theory, then it is possible to determine some of the O ϵ 2 eigenvalues without including higher-order corrections. This possibility is a consequence of the fact that a symmetry-preserving perturbation will not generate a symmetry-violating term, so Λ finds a block-triangular form and the two sectors can be diagonalized separately. In the equilibrium limit of our model (in the sector σ = 1 when v = 1), a similar situation occurs with respect to statics and dynamics, where perturbations in the dynamics (w) cannot affect the behavior of the statics (u). This structure makes it possible to inspect the stability of w up to O ϵ 2 at the same order of the RG calculations. In the full nonequilibrium model, the statics is not decoupled from the dynamics, so the stability in w cannot be determined using such an approach. However, for the equilibrium fixed points, equilibrium plays a role similar to a higher symmetry because equilibrium perturbations do not generate nonequilibrium terms. As a result, it is possible to determine the stability in v for the two coupled equilibrium fixed points, and we find both to be stable in v. However, for all of the coupled fixed points, w remains marginal. In short, to the lowest order in our perturbative expansion, the system flows to the NEFP (equilibrium fixed point) in the σ = −1 (σ = 1) sector. While, in principle, nonperturbative effects or higher-order terms in ϵ could modify this behavior, this is a generic feature of perturbative RG and not specific to our model. A qualitative sketch of the expected RG flow is illustrated in Fig. 2 in terms of the original g 12 , g 21 couplings.
Finally, we remark that in the case of the original microscopic model, A 02 and B 11 have opposite signs, which thus carries over to the relative sign of g 12 and g 21 . Thus, it is plausible to expect a critical behavior governed by the NEFPs.

C. Universal scaling behavior
Any fixed point-equilibrium or not-exhibits critical behavior and exponents characterizing correlations and dynamics among other properties of the system. In particular, we consider the anomalous dimensions η and η′ of the original and response fields, the dynamical critical exponent z, as well as the exponent ν characterizing the divergence of the correlation length as the critical point is approached. These exponents describe the scaling behavior of the correlation and response functions at or near criticality as where C i , χ i are general scaling functions. We have dropped the subscript i from η, η′, z due to the strong dynamic scaling and in anticipation of the same spatial scaling dimensions for the two fields; however, we have kept the subscript in r j for j = 1, 2 since the RG equations couple them in a nontrivial way.
The exponents at the fixed point can be extracted via what is known as the method of characteristics (see Appendix D for details). Noting that, for fixed bare (microscopic) parameters, the correlation and response functions are not affected by changing the RG momentum scale μ, we can relate these critical exponents to the flow functions as The renormalization of the parameters r j and the corresponding exponent ν j requires a special treatment and will be discussed later in this section. At the nonequilibrium critical point, we find [cf. Eqs. (31) and (37) Interestingly, we see that in contrast to an equilibrium fixed point where the temperature becomes scale invariant, the effective temperature at the NEFPs changes with the scale. In particular, the system becomes "hotter" at longer length scales since γ T < 0. Using Eq. While, in equilibrium, η = η′ as a consequence of the fluctuation-dissipation theorem, we have η ≠ η′ since the temperature itself is scale dependent, γ T ≠ 0, at the NEFP. Note also that the critical exponents z, η, η′ are the same for both fields. While strong dynamic scaling already guarantees the same dynamical critical exponent, the anomalous dimensions are also identical, owing to the emergent symmetry of the fixed point where u 1 R = u 2 R and v R = w R = 1. However, the latter do not reflect any actual symmetry of the model and could be modified at higher orders in the epsilon expansion.
An interesting feature of the NEFPs is that η < 0 to the first nontrivial order in the epsilon expansion. This result is in contrast with equilibrium where η > 0, a fact that can even be proved on general grounds (e.g., unitarity in a related quantum field theory) [92]. If this feature (η < 0) extends beyond perturbation theory to, say, two dimensions, it would indicate that the correlation function (C(r) ∝ |r| −d+2−η ) diverges at large distances. However, this would invalidate the starting point of our field-theoretical treatment based on an expansion in field powers since large-scale fluctuations grow without bound. However, it might also indicate the absence of ordering in low dimensions. This possibility seems particularly natural in light of the effective temperature increasing at larger scales, which in turn tends to disallow ordering in low dimensions. While this may be an artifact of perturbative RG, it indicates that the behavior of the NEFPs in low dimensions is governed by different principles than their equilibrium counterparts. Finally, we note that, to the lowest nontrivial order considered, the dynamical critical exponent z is related to η′ in an identical fashion as in equilibrium.
Next we consider the renormalization of the mass terms. This case requires special care as their renormalization is intertwined. Similar to our redefinition of u, we should instead consider the renormalization of r i /D i so that we only need to consider two flow equations, which is consistent with the scaling analysis in Appendix D. In a slight abuse of notation, we simply replace r → D i r i . Defining μ(l) = μl and the flowing parameters r i (l) with r i (1) = r i R , we find the flow equations [cf. Eqs. (31) and (37) together with the Z factors in Appendix C] l dr 1 (l) dl = γ r 1 r 1 = −2 + ϵ 2 r 1 (l) ± ϵ 2 3 r 2 (l), l dr 2 (l) dl = γ r 2 r 2 = −2 + ϵ 2 r 2 (l) ∓ ϵ 2 3 r 1 (l), where the ± refer to the two NEFPs with opposite signs of u 12 R . The flow equations can be solved as r 1 (l) = l −1/ν′ r 1 R cos logl ν″ + r 2 R sin logl ν″ , r 2 (l) = l −1/ν′ r 2 R cos logl ν″ − r 1 R sin logl ν″ , where These equations can be cast in a more compact notation as r 1 (l) + ir 2 (l) = l −1/ν′ − i/ν″ r 1 R + ir 2 R .
Defining r ≡ r 1 + ir 2 , we can recast this equation as r(l) = l −1/ν r R , where the critical exponent ν emerges as Interestingly, the exponent ν becomes complex valued at the NEFP. We can then express the scaling functions in Eq. (39) as while r R = r 1 R + ir 2 R = r 1 R 2 + r 2 R 2 , while ∠r R denotes the polar angle in the r 1 R − r 2 R plane. Additionally, P is a 2π-periodic function. To obtain these equations, we have used the fact that r/l 1/ν′ + i/ν″ can instead be written as a function of |r | /l 1/ν′ and e i(logl)/ν″ − i∠r . The former expression often appears in scaling functions of this type and characterizes the scaling of the correlation length; however, the latter gives rise to a log-periodic function as a change of log l → log l + 2πν′′ leaves the exponential invariant.
The appearance of log-periodic functions has important consequences for the critical nature of the fixed points. They lead to a discrete scale invariance rather than the characteristic continuous scale invariance at a typical critical point [69]. Rather than a self-similar behavior at all length scales, a preferred scaling factor emerges as b * = e 2πν″ , (49) rescaling by which, or any multiple integer thereof, leaves the system scale invariant. In this sense, discrete scale invariance mimics a fractal-like structure, in which rescaling the system by a particular factor leaves the system self-similar. Note, however, that the discrete scale invariance and the fractal-like structure only emerges at long length scales (in the continuum) as opposed to the microscopic structure of a fractal (in discrete space). Additionally, if we were to consider, e.g., the effect of a physical momentum cutoff Λ, this would enter the periodic function as a phase shift, thus determining the phase of the oscillations.
Similar phenomena appear to arise in stock markets [93], earthquakes [94], equilibrium models on fractals [95], and several other systems [69]. Log-periodic functions and the emergence of a preferred scale have been identified in the early developments of renormalization group theory [96][97][98], but they have been rejected as artifacts of positionspace RG. On the other hand, their recent surge in diverse contexts from earthquakes to stock markets has instead relied on simple dynamical systems (with one or few variables) where the dynamics involves a discrete map itself [69]. This phenomenon has also emerged in recent works in the context of driven-dissipative quantum criticality [39,99] as well as the dynamics of strongly interacting nonequilibrium systems [100]. A particularly well-known example of RG limit cycles is the behavior of Efimov states [101,102], whose binding energies form a geometric progression similar to discrete scale invariance. These quantum RG limit cycles have been noted to be closely related to Berezinskii-Kosterlitz-Thouless phase transitions [103,104]. Disordered systems provide another context where complexvalued exponents and discrete scale invariance have been noted in both classical [105][106][107][108][109] and quantum [110] settings. The discrete scale invariance reported in this work, however, appears to be unique as it has emerged in an effectively classical yet nonequilibrium model in the absence of disorder.
The discrete scale invariance approaches a continuous one as the upper critical dimension, d c = 4, is approached. In three dimensions, perturbative values at the NEFP (with ϵ = 1) yield a very large scaling factor (b * ~ 10 9 ); however, with the exponential dependence on the critical exponents, the scaling factor is sensitive even to small corrections of the exponent beyond the lowest-order perturbation theory. Nevertheless, our results should be viewed as a proof of principle for the emergence of discrete scale invariance in macroscopic nonequilibrium systems. Additionally, higher harmonics in the periodic function P can be significant, which then should be observed over smaller variations in the physical scale.
Finally, we elaborate on a possible connection between the log-periodic behavior and limit cycles. Indeed, the microscopic mean-field phase diagram near the multicritical point also includes a limit-cycle phase that displays persistent oscillations. For a rapidly oscillating limit cycle, the corresponding phase transition can be described from the viewpoint of a rotating frame (defined by the oscillation frequency) and by making the rotating-wave approximation. With this mapping, a limit-cycle phase transition is no different from a dissipative phase transition with an emergent U (1) symmetry [26]. Near our multicritical point, however, the frequency of oscillations becomes small, and thus no such mapping is possible. On the other hand, the discrete scale invariance discussed above also leads to an oscillatory behavior (in both space and time), albeit one that is log-periodic. Nevertheless, a natural possibility is that, at some intermediate regime away from the multicritical point, the discrete scale invariance merges into a limit-cycle solution. Moreover, we see that in the doubly ordered phase, the Liouvillian gap becomes complex valued (see Sec. IV F). This result furthers the possible connection to the limit-cycle phase as a nonzero imaginary part implies that the system undergoes oscillations-which, however, decay-as the steady state is approached.
In Table I, we summarize all of the fixed points (aside from those involving the trivial Gaussian fixed point) and their critical exponents to the lowest order.

D. Phase diagram
The phase diagram itself is distinct in the vicinity of the NEFPs. In contrast to their equilibrium counterparts, these fixed points only give rise to a tetracritical point. With the effective magnetic field set to zero, h = 0, four different phases emerge: A disordered phase with ϕ 1 = ϕ 2 = 0; two phases with either ϕ 1 ≠ 0 corresponding to bistability or ϕ 2 ≠ 0 leading to antiferromagnetic ordering; and, finally, a doubly ordered phase where both fields become ordered, ϕ 1 ≠ 0 ≠ ϕ 2 . While the first three phases also emerge in the mean-field theory of the microscopic model, the doubly ordered phase only arises in the course of RG when the A 20 term is generated.
The phase boundaries are governed by the scaling behavior of r i . Let us set the effective magnetic field to zero, h = 0, and consider the scaling functions characterizing the correlation and response functions in Eq. (39). To determine the phase boundary, it suffices to take the limit ω, q → 0. In this case, the scaling functions are solely determined as a 2πperiodic function of (ν′/ν′′) log (|r R |) − ∠r R ; this is achieved by eliminating the momentum scale in Eq. (48) in favor of r R . Since the correlation functions only depend on the mass terms through the above combination, the phase boundary itself-characterized by the divergence of correlations-arises at a fixed value of this quantity (up to integer multiples of 2π). Therefore, the shape of the phase boundary is given by ν′ ν″ log r R − ∠r R = const, (50) which is a spiral, leading to the phase diagram illustrated in Fig. 4. Similar to our discussion of the discrete scale invariance, the perturbative values of the exponents at ϵ = 1 require very large scales to observe a full spiral. But again these scales are highly sensitive to corrections to perturbative RG. Moreover, partial spirals can still be observed for reasonable scales.
Since the spiraling boundaries all spiral in the same direction, distinguishing them from equilibrium critical points, the effects of this may be seen even for very weak spiraling. Additionally, the two NEFPs are distinguished from each other by the direction of spiraling since each has a different sign of ν′′.
When the effective magnetic field is nonzero, there is no such symmetry as ϕ 1 → −ϕ 1 , leading to a surface of first-order phase transitions where ϕ 1 undergoes spontaneous symmetry breaking. This first-order transition occurs in both the uniform phase (defined by the B phase at h = 0) and the antiferromagnetic phase (defined by the B AF phase at h = 0).
In both cases, the average population changes discontinuously while the sublattice population difference changes continuously. Finally, we take into account the magnetic field in determining the phase boundaries by considering an effective mass as r R + ℎ R 1/βδ , where β and 1/δ describe the scaling behavior of the order parameter with r and h, respectively, within the ordered phase. The effect of the magnetic field h is to "unravel" the spirals for small r i . The corresponding phase diagram for the nonzero effective magnetic field is illustrated in Fig. 7.

E. Hyperscaling relations
As η ≠ η′ and the exponent ν is complex, the standard hyperscaling relations will be modified. For example, consider the exponents characterizing the order parameter and magnetic susceptibility, In the standard equilibrium setting, these exponents are related to ν, η via To see how the above hyperscaling relations are modified near the NEFPs, first note ∂ ϕ ∂ℎ = lim q, ω 0 χ q, ω, r j , (53b) and that in the ordered phases, 〈ϕ 2 〉 and 〈ϕ〉 2 will have the same scaling behavior. Similar to our phase-diagram analysis, we take the limit ω, q → 0 or t, |x| → ∞ and allow |r R | to define the momentum scale in place of |q|. This method amounts to replacing the |q| prefactors in Eq. (39) with r R −ν′ . We then find the scaling of the order parameter and magnetic susceptibility to be where P C , P χ are 2π-periodic functions. Given this scaling behavior, there are two approaches to identifying the exponents β, γ. The first is that, like ν, the exponents β and γ assume complex values, corresponding to the log-periodic behavior. However, this is simply a reflection of the structure of the phase diagram: Approaching the critical point directly, the system crosses in and out of all four phases, giving rise to discrete scale invariance. The second and perhaps more natural approach is to tune the system to the critical point along paths that fix the argument of the periodic functions. In doing so, the angle relative to the phase boundaries is fixed. According to this second approach, which confines the discrete scale invariance to ν, we have the new hyperscaling relations β = ν′(d − 2 + η)/2, γ = ν′ 2 − η′ , (55) where only the real part ν′ enters.

F. Liouvillian gap closure
In this section, we investigate how the Liouvillian gap closes upon approaching the multicritical point. In contrast to the equilibrium fixed points where the gap always closes along the real axis (hence, purely dissipative or relaxational dynamics), the NEFPs exhibit a qualitatively different behavior with the gap closing along a complex path, indicating an interplay between reversible and irreversible relaxation in this phase.
We consider the system in the doubly ordered phase where both fields take nonzero expectation values M i ; for notational convenience, we make the change of variables ϕ i → ϕ i + M i , where the fields ϕ i now represent the fluctuations around the order parameter. In addition to the original action, this transformation introduces new quadratic and linear terms as (including the original r 1 and r 2 terms, too) ∫ x, t r 1 + 3u 1 M 1 2 + u 12 M 2 2 ϕ 1 ϕ 1 + r 2 + 3u 2 M 2 2 + σu 12 M 1 2 ϕ 2 ϕ 2 + 2u 12 M 1 M 2 ϕ 2 ϕ 1 + 2σu 12 M 1 M 2 ϕ 1 ϕ 2 + M 1 r 1 + u 1 M 1 2 + u 12 M 2 2 ϕ 1 + M 2 r 2 + u 2 M 2 2 + σu 12 M 1 2 ϕ 2 . (56) In addition, several cubic terms are also introduced, which are not reported for simplicity. We set the vertices ϕ 1 and ϕ 2 to zero since, by definition, ϕ i solely represent the fluctuations.
where the coupling terms have been replaced with their renormalized values due to the inclusion of fluctuations in the form of counterterms. The other parameters are not renormalized at this order, but if we were to include higher-order fluctuations, they too would be replaced with their renormalized values since the ordered phases do not give rise to any new Z factors.
Putting these terms together with the quadratic part of the action, we find S 0 = ∫ x, t ∑ i ϕ i ζ ∂ t − D ∇ 2 + R i ϕ i − ζT ϕ i 2 + R 12 (ϕ 2 ϕ 1 + σϕ 1 ϕ 2 ), (58) where R i = 2u i R M i 2 and R 12 = 2u 12 R M 1 M 2 . The poles of the corresponding propagators are then obtained as 0 = σR 12 2 − Dk 2 + R 1 + iζω Dk 2 + R 2 + iζω (59) or, more explicitly, as From this equation, we can find a simple condition for when the poles do not lie along the imaginary axis(corresponding to the negative real eigenvalues of the Liouvillian) as Indeed, in equilibrium, where σ = 1, this condition cannot be satisfied. This implies that the relaxation is purely relaxational in equilibrium as expected (in this case, for model A). However, at the NEFPs where σ = −1, the above condition can be satisfied. To see why this is the case, let us cast the above condition for σ = −1 in terms of u and M i as Recalling that u 1 R = u 2 R at the NEFPs, at least to the lowest order in the epsilon expansion, the above condition is trivially satisfied whenever |M 1 | = |M 2 |. In this case, the pole with the lowest nonzero decay rate takes the form (with |M 1 | = |M 2 |≡ M and k → 0) −iζω = 2M 2 u 1 R ± iu 12 R . (63) In fact, with |M 1 | = |M 2 |, the Liouvillian gap achieves its largest imaginary value relative to its real part. We can identify the ratio of the imaginary part to the real part of the gap as u 12 R /u 1 R = 3, which corresponds to the Liouvillian gap closing at the angle π/3 with respect to the real axis. Again, higher orders in the epsilon expansion may modify the value of this angle. A generalization of the model considered here to the O(N) × O(N) model of Ncomponent vectorlike order parameters leads to a similar behavior. In fact, we find that, for N = 2 and N = 3, the corresponding Liouvillian gap closes at the angles π/4 and π/6, respectively. We should then conclude that different nonequilibrium universality classes of our model and its generalization give rise to different angles of the Liouvillian gap closure in the complex plane. Further details on these generalized models will be presented in followup papers.
The scaling of (the magnitude of) the gap itself as a function of the distance from the critical point can be directly obtained by observing that the gap defines an inverse timescale, which itself is associated with the exponent z. Thus, the gap scales as |r| zν′ , where the exponent ν′ is due to the scaling of r itself. With the order parameters M 1 and M 2 scaling similarly, the angle that defines the gap closure in the complex plane only depends on (the absolute value of) their ratio. As remarked earlier, this angle achieves its maximum when |M 1 | = |M 2 |. We further note that the gap is purely real (relaxational) near phase boundaries where only one of the order parameters undergoes a transition since the lhs of Eq. (62) would be suppressed compared to the rhs.
Finally, much like the discrete scale invariance in the previous section, the complex Liouvillian gap is somewhat reminiscent of limit cycles, although a true limit-cycle phase is characterized by purely imaginary eigenvalues that characterize the steady state itself.

V. EXPERIMENTAL REALIZATION
An ideal avenue for realizing these multicritical points is via the use of cavity or circuit quantum electrodynamics (QED). Individual cavities and circuits have been studied experimentally in great depth due to their potential applications in quantum computation [14,111,112]. Furthermore, both cavity QED and circuit QED have been proposed as platforms for realizing many-body states of light via nearest-neighbor coupling arrays of cavities or circuits [113][114][115][116]. Generally, these cavities and circuits have non-negligible loss due to dissipation. While dissipation is detrimental when it comes to realizing the quantum ground state of a given system, it is a crucial ingredient in realizing driven-dissipative phase transitions. There have been a variety of theoretical proposals to realize different drivendissipative models in cavity-and circuit-QED systems [27,30,31,65,117,118]. Recent experiments have even identified a driven-dissipative many-body phase transition [15].
For the model considered in this work (see Sec. II), many-body experimental platforms already exist that include drive and hopping, as well as dissipation. The remaining ingredient is then the nearest-neighbor interaction [the quartic term in Eq. (1)] to be contrasted with a Hubbard term that characterizes on-site interaction. Both types of interaction are generally known as Kerr nonlinearities; we are interested in what is known as a cross-Kerr nonlinearity, which has been utilized experimentally in several few-mode systems [119][120][121]. A more general version of our model has been considered in Refs. [30,31], along with a discussion on how the nonlinear interaction terms can be tuned experimentally via Josephson nano-circuits. In a recent theoretical proposal, a setting consisting of a capacitor in parallel with a superconducting quantum interference device is put forth as an alternative means of achieving tunable Kerr nonlinearities [62].
While generic experimental settings introduce other nonlinear terms (e.g., Hubbard interactions and correlated hopping) in addition to the density-density interactions, we do not expect them to dramatically affect the results of this paper. While such terms can change the location of the multicritical point [30,31], the universal properties of the latter should not be affected by the details of the microscopic model.
We close this section with a discussion of the sign of various terms (e.g., the negative cross-Kerr nonlinearity) arising in the proposals of Refs. [30,31,62].While a negative interaction term could lead to unbounded energy spectra, it would not pose a problem in the context of driven-dissipative systems where the steady state is not concerned with a minimum-energy ground state. Furthermore, one can change the sign of various terms in the Hamiltonian of a driven-dissipative system with a proper mapping [122]. For example, by sending Ω → −Ω and a → −a on one of the two sublattices, the sign of J can be changed while leaving the remaining terms fixed. Similarly, one can also map H → −H by taking the complex conjugate of the master equation, which, together with the previous mapping, allows an appropriate choice for the sign of J, V. Finally, the overall phase of Ω is unimportant, while the parameter Δ can be easily tuned to a desired sign.

VI. CONCLUSION AND OUTLOOK
In this work, we have considered an experimentally relevant driven-dissipative system where two distinct order parameters emerge that characterize a liquid-gas-type transition (associated with the average density) as well as an antiferromagnetic transition (associated with the difference in the sublattice density). The two phase transitions coalesce and form a multicritical point where both transitions occur at the same time. We have investigated the nontrivial interplay of two order parameters at the multicritical point. Using a fieldtheoretical approach-appropriate in the vicinity of the phase transition-we have shown that the critical behavior at this point can be mapped to a nonequilibrium stochastic model described by a ℤ 2 × ℤ 2 symmetry. Using perturbative renormalization group techniques, we have determined the RG flow equations of the model and identified a pair of new classical nonequilibrium fixed points that exhibit several exotic properties. First, we obtain two different exponents for the critical scaling of fluctuations and dissipation at the critical point, underscoring the violation of the fluctuation-dissipation relation at all scales and resulting in a behavior where the system becomes hotter and hotter at larger and larger scales. Furthermore, these NEFPs are distinguished by the emergence of discrete scale invariance and a complex Liouvillian gap even close to the critical point. Additionally, the phase diagram near these multicritical points displays spiraling phase boundaries. The latter properties could be particularly useful in identifying these NEFPs in experiments.
While generic driven-dissipative phase transitions tend to have effective equilibrium dynamics, we have shown that the interplay between several order parameters (in this case, two) could very well lead to exotic nonequilibrium behavior. This perspective opens a new avenue to investigate and experimentally realize nonequilibrium phases and phase transitions in the context of driven-dissipative systems without relying on the engineering of complicated nonlocal or non-Markovian dissipation. models involving vectorlike order parameters [70][71][72][73][74][75][76]. While a driven-dissipative condensate of polaritons has been investigated theoretically in detail [33,34], recent experimental studies into condensate supersolids [123][124][125][126][127] can provide excellent platforms for probing any emergent NEFPs. In addition to the U(1) symmetry of the condensate, the two coupled optical cavities can provide either an additional ℤ 2 × ℤ 2 symmetry (corresponding to a lattice supersolid) or an approximate O(2) symmetry (corresponding to a continuous supersolid).

APPENDIX A: LANGEVIN EQUATIONS NEAR THE MULTICRITICAL POINTS
In this Appendix, we present the details of the derivation of the Langevin equations in the main text. To this end, we follow the procedure detailed in Ref. [20]. We begin by constructing the Keldysh path integral, then identify a semiclassical limit, and derive a pair of complex Langevin equations that describe the dynamics near the steady state. Finally, we identify a pair of two massless real fields (i.e., soft modes) and two massive real fields (i.e., fast modes). We adiabatically eliminate the massive fields to obtain a pair of Langevin equations presented in the main text.
We first ignore the sublattice symmetry for simplicity; this would not affect the analysis presented here. We return to the latter symmetry once we identify the semiclassical limit and corresponding Langevin equations. We cast our model in terms of a Keldysh path integral as where the action S K is defined as S K = ∫ x, t ψ q * ∂ t ψ cl + ψ q ∂ t ψ cl * − ∫ x, t H n ψ cl + ψ q − H n ψ cl − ψ q + iΓ ∫ x, t ψ q 2 − ψ cl ψ q */2 − ψ cl *ψ q /2 , with ψ cl/q the classical/quantum fields and H n /ψ the normal-ordered form of the Hamiltonian. The third line corresponds to the particular case of the Lindblad operator Γa, although this approach may easily be extended to more general Lindbladians [36]. With the Hamiltonian in Eq. (1), the action in the continuum (with the nearest-neighbor interactions expanded in powers of the gradient) is given by S K = ∫ x, t ψ q * i ∂ t + J ∇ 2 + Δ + zJ + i Γ 2 ψ cl + c . c .
− ∫ x, t V ψ cl 2 ∇ 2 + z ψ cl ψ q * − 2Ωψ q * + c . c . + ∫ x, t iΓ ψ q 2 − V ψ q 2 ∇ 2 + z ψ cl ψ q * + c . c . , where z = 2d is the coordination number. This expression bears a close resemblance to the action of Eq. (12) in Ref. [20], but they differ in the form of their interactions (which involve gradient terms here) and because of our use of normal ordering rather than the Weyl ordering of Ref. [20]. This motivates a similar rescaling of the parameters as The parameter N effectively describes a density scale for the microscopic model via Since varying the density scale also modifies the interaction energy per particle, the interaction strength should be reduced correspondingly such that V|ψ cl | 2 =v|Ψ cl | 2 ; similarly, the drive should be increased so that Ωψ q = ΩΨ q . We can then rewrite the action as S K = ∫ x, t Ψ q * i ∂ t + J ∇ 2 + Δ + zJ + i Γ 2 Ψ cl + c . c .
In the limit of large N, the last term (the second term in the last line) can be dropped, leading to an action that is at most quadratic in Ψ q . This result is simply because a large population N corresponds to the semiclassical limit represented by a large classical field ψ cl and small fluctuations due to the quantum field ψ q . Using this fact, we can map the action to a Langevin equation as [68] i ∂ t Ψ = − J ∇ 2 + Δ + zJ + iΓ/2 ξ*(x, t)ξ x′, t′ = Γ N δ x − x′ δ t − t′ , 〈ξ〉 = 0. (A6b) Note that the noise level is further suppressed at larger N as should be expected from our semiclassical treatment.
Next, we include the sublattice symmetry by defining Ψ 1 as the sublattice average and Ψ 2 as the sublattice difference of the field Ψ; see Eq. (15), which differs by a factor of N due to our semiclassical limit. Our new Langevin equations are now Notice that we have dropped all the gradient terms as they do not play a role in identifying the massive fields and their adiabatic elimination. We also follow our convention in the main text to set zJ J.
Our model exhibits two multicritical points where two modes (each a component of one of the two fields Ψ i ) become critical. Because of the sublattice symmetry, Ψ 2 = 0 at the multicritical points up to fluctuations. Working in units where Δ + J = 1, the two multicritical points occur at with Ψ 1 = Ψ c = 2/3ve −iπ/3 .
Next, we expand the Langevin equations in the vicinity of the two multicritical points as The soft and gapped modes can be determined as linear combinations of the real and imaginary parts of the two fields ψ 1 , ψ 2 . As in the case of bistability, the first pair is identified as [20] ψ 1 = ϕ 1 ′ + e iπ/3 ϕ 1 , where ϕ 1 ′ is massive and relaxes quickly while ϕ 1 defines the slow field. Identifying the massive or massless components of the field ψ 2 critical point as Δ c = 1/3: ψ 2 = 1 3 ϕ 2 e −iπ/6 + ϕ 2 ′e iπ/6 , Δ c = 2/3: ψ 2 = 1 3 ϕ 2 + ϕ 2 ′e iπ/3 , where again the primed (unprimed) field indicates the massive (massless) field. Note that these differ from the main text by a factor of V , which is done to simplify the resulting parameters in the Langevin equations by moving all the V dependence to the noise term.
Upon adiabatically eliminating the massive fields and restoring the gradient terms, we arrive at the Langevin equations φ 1 = ℎ − r 1 ϕ 1 + D 1 ∇ 2 ϕ 1 + ξ 1 + A 20 ϕ 1 2 + A 02 ϕ 2 2 + A 12 ϕ 1 ϕ 2 2 + A 30 with Gaussian noise ξ i (x, t)ξ j x′, t′ = 2 κ i N δ ij δ x − x′ δ t − t′ . (A13) The various numerical factors are summarized in Table II for the two multicritical points under consideration. Note that ζ i = 1 and T i = κ i /N. We can see at this point that the opposite signs of A 02 and B 11 indicate that no Hamiltonian description is possible. Indeed, as will be discussed in the following section, this will carry over to the signs of g 12 and g 21 , leading to the critical behavior defined by the NEFPs.

APPENDIX B: REDUNDANT OPERATORS
In this Appendix, we identify the redundant operators in the Langevin equations (A12). In general, this can be done at the level of the Schwinger-Keldysh action; however, we focus on the equivalent description in terms of the Langevin equations. This perspective is particularly suitable in dealing with the (Itô) regularization that is required to properly define the stochastic equations.
Consider a pair of Langevin equations that define an Itô process in the differential form [128] dϕ 1 = f 1 ϕ 1 , ϕ 2 dt + κ 1 dW 1 , dϕ 2 = f 2 ϕ 1 , ϕ 2 dt + κ 2 dW 2 , where dW i is the stochastic noise that obeys the Itô rules: In order to identify the redundant operators, we should examine the Langevin equations under a general change of the field variables. We should then find the dynamics in terms of new variables defined as Φ 1 = g 1 (ϕ 1 , ϕ 2 ) and Φ 2 = g 2 (ϕ 1 , ϕ 2 ); the functions g i are general (but local) nonlinear maps that are invertible in a neighborhood around the multicritical point and that preserve the sublattice symmetry ϕ 2 → −ϕ 2 . We have used the Itô rules to derive the above equation, which is known as Itô's formula or Itô's lemma [128]. A similar stochastic equation can be derived for Φ 2 by switching 1 ↔ 2.
Note that the terms on the rhs should be expressed in terms of Φ i through the inverse functions g i −1 . One notices that there are new contributions to the deterministic dynamics due to the noise. However, since we are working under the assumption that the noise is parametrically small compared to the deterministic terms (with a strength proportional to 1/N), we can ignore such terms. Additionally, the noise terms are no longer additive but are instead multiplicative, which introduces new terms in the action. However, since κ i ≠ 0, the latter are irrelevant in the sense of RG and can be neglected as well. The same holds for nonlinear terms (beyond quadratic terms) that involve gradients.
Without loss of generality, we assume that ∂g i ∂ϕ j ϕ i , ϕ j = 0 = δ ij ; (B5) rescaling the fields by a constant factor does not allow us any additional freedom, while rotations obscure the symmetry ϕ 2 → −ϕ 2 . Additionally, we do not consider a constant shift in the field ϕ 1 (a shift in ϕ 2 is disallowed due to symmetry) for now, but we discuss it separately later in this Appendix. Based on the structure of Eq. (B4), we notice that the quadratic terms in f i and g i result in additional cubic terms. All other new terms in the deterministic part of the dynamics involve fourth-or higher-order terms that are irrelevant under RG. Expressing a general nonlinear transformation as g 1 ϕ 1 , ϕ 2 = ϕ 1 + c 20 ϕ 1 2 + c 02 ϕ 2 2 , (B6a) g 2 ϕ 1 , ϕ 2 = ϕ 2 + c 11 ϕ 1 ϕ 2 , the modification of the cubic terms to lowest order in the coefficients and interaction terms is given by The effective mass and magnetic field terms also change, but this simply shifts the location of the critical point. Having exhausted the three redundant operators to remove the three quadratic terms, there is no further freedom in tuning other terms, and, specifically, all the cubic terms are fixed. While we could, in principle, include cubic or higher-order terms in the nonlinear transformation [Eq. (B6)], these would only modify the fourth-or higher-order terms, which are irrelevant under RG due to the presence of the cubic terms. We also see that the relative sign of A 12 and B 21 is indeed directly determined by the relative sign of A 02 and B 11 , leading to criticality described by the NEFPs. Finally, we note that the coefficient A 20 appears in the denominator of the above transformations. We assume that this term is generated under coarse graining and thus should pose no problem in making the above transformations. However, if there is a mechanism where this coefficient could be tuned to zero, the above transformations are no longer valid and the two nonzero cubic terms should be kept.

APPENDIX C: PERTURBATIVE RG
In this Appendix, we discuss the details of the calculations in our perturbative RG analysis.
In the first part, we introduce the diagrammatic techniques we have used in the main text. In the second part, we compute the one-loop diagrams, while, in the third part, we compute the two-loop diagrams for terms that are unrenormalized at one loop.

Diagrammatic techniques
To define the Gaussian propagators, we start with the Gaussian model with the corresponding action The Gaussian response and correlation functions are then given by where ℱ denotes the Fourier transform in both space and time. These propagators can be expressed in a diagrammatic representation as shown in Fig. 8.
The four interaction vertices from Eq. (27b) are illustrated in Fig. 9. Because of the structure of the action, one can find the corresponding Z factors for ϕ 2 by switching the subscripts 1 ↔ 2 and multiplying u 12 by a factor of σ. Finally, in order to determine the Z factors, we employ dimensional regularization using the minimal subtraction procedure. This means that only the ultraviolet divergences, in the form of powers of 1/ϵ, are incorporated into the Z factors. For simplicity, we only present these divergences in the evaluation of integrals in the following sections; nondivergent terms are dropped. Technical details for evaluating integrals of the type presented here may be found in, e.g., Ref. [68].

One-loop diagrams a. Mass terms
In this section, we consider corrections to r 1 . The corrections to r 2 can be obtained by switching the roles of ϕ 1 and ϕ 2 . There are two diagrams that provide corrections, which are illustrated in Figs. 10(a) and 10(b). The combinatorial and interaction factors are (a) 3u 1 and (b) u 12 . Thus, the one-loop contribution to the renormalization of r 1 is so we should evaluate the integral of C 0 i . The latter can be written as where the last line follows once we integrate over frequency. The phase transition occurs where the renormalized mass term vanishes. The critical value of the mass parameter r i c is then determined by and similarly for r 2 c ; note that the factors of r in the denominator have been dropped as they introduce O(u 2 ) corrections. Defining an additive renormalized mass term r i = r i − r i c , we can determine the Z factors for r i . To this end, we need to compute integrals of the form One-loop corrections to (a,b) r 1 ϕ 1 ϕ 1 , (c,d) u 1 ϕ 1 3 ϕ 1 , and (e)-(h) u 12 ϕ 2 2 ϕ 1 ϕ 1 . Analogous diagrams for r 2 ϕ 2 ϕ 2 , u 2 ϕ 2 3 ϕ 2 and σu 12 ϕ 1 2 ϕ 2 ϕ 2 can be obtained by switching thin black and thick cyan lines.
where D i = ζ i D i and r i = ζ i r i . Noting that r i = μ 2 r i R + O(u) and r i R is a finite constant, then according to the minimal subtraction procedure, r j + r i x 2 − d/2 ≈ μ −ε r j R + r i R x where ϵ has been set to 0 in the nondivergent part to extract the residue of the pole.

FIG. 11.
Two-loop corrections to (a)-(c) ζ 1 and D 1 as well as (d,e) ζ 1 T 1 . Analogous diagrams for ζ 2 , D 2 , and ζ 2 T 2 can be obtained by switching thin black and thick cyan lines.
where i corresponds to the field with one propagator and j to the field with two.

APPENDIX D: METHOD OF CHARACTERISTICS
In this section, we employ the method of characteristics in order to derive the scaling behavior of the correlation and response functions at or near a given fixed point. Since the correlation and response functions do not depend on the renormalized parameters, they are independent of the momentum scale of renormalization μ. Thus, for the response functions, we have μ d dμ χ i q, ω, p R , u R = 0, where {p} = {r i , ζ i , D i , T i } and {u} = {u 1 , u 2 , u 12 } are the interaction strengths.
Additionally, we define a scaling function via χ i = μ −2 χ i , where the scaling factors are due to the scaling dimensions of the fields as well as the delta functions-which are factored out -corresponding to momentum and energy conservation. We can rewrite the total derivative μ(d/dμ) in terms of the partial derivatives with respect to other parameters as Schematic illustration of the physical setup. The contrast of the two checkerboard (light yellow and dark blue) sublattices defines the antiferromagnetic order parameter. On each sublattice, there is a low-population (low n ph ; solid curve) and a high-population (high n ph ; dashed curve) steady state, corresponding to the bistability order parameter. The large arrow indicates the drive, while the wavy arrows represent the dissipation. J and V denote the hopping and nearest-neighbor interactions, respectively.

FIG. 2.
A schematic RG flow diagram projected to the g 12 − g 21 plane. (The full RG flow requires a five-dimensional space, which precludes a more complete sketch of the RG flow in this twodimensional space; see Sec. IV.) In addition to the equilibrium fixed points where g 12 = g 21 (green circles), a pair of stable NEFPs (orange diamonds) emerge in the sector defined by the opposite signs of g 12 and g 21 . These new fixed points exhibit exotic critical behavior reflecting their truly nonequilibrium nature. Filled (black) arrows represent the stability, while partial (gray) arrows indicate the expected stability of the various fixed points in different directions. Stability is known to lowest order in ϵ = 4 − d only in directions that preserve the ratio g 12 /g 21 . The RG flow cannot cross the g 12 , g 21 axes, which is represented by double lines. The stable equilibrium fixed point is characterized by an emergent O (2) symmetry, while the unstable equilibrium fixed points include a biconical fixed point as well as various decoupled fixed points (which all lie at the origin in this diagram) corresponding to combinations of Ising and Gaussian fixed points.

FIG. 3.
Summary of the main features of the NEFPs contrasted with effective equilibrium fixed points. (a) Schematic correlation functions. A generic continuous scale invariance characteristic of criticality is reduced to a discrete scale invariance at the NEFP. (b) Effective temperatures T i representing the two order parameters as a function of the length scale (q −1 ).
The two temperatures become identical at long length scales, but while they approach a constant at the equilibrium fixed point, they diverge at large scales at the NEFPs. (c) Gap closure upon approaching the critical point. Here, Γ ℒ denotes the Liouvillian gap with the real part describing the relaxation rate (a.k.a. the dissipative gap) and the imaginary part characterizing the "coherence gap." For the equilibrium fixed point, the gap can close only along the real line, indicated by the arrow. In contrast, the gap for the NEFP can take complex values and close along any path lying in the shaded region, making a maximum angle of π/3 with the real line. (d) Critical exponents to lowest nontrivial order in ϵ = 4 − d.
The exponent ν, typically associated with the divergence of the correlation length, becomes complex valued at the NEFPs, with its imaginary part characterizing the discrete scale invariance [cf. part (a)]. Here, η and η′ are anomalous dimensions characterizing fluctuations and dissipation with η ≠ η′ at the NEFPs indicating the violation of the fluctuation-dissipation theorem. Note that z is the dynamical critical exponent.

FIG. 4.
Phase diagram associated with the NEFPs of the nonequilibrium Ising model of two coupled fields for h = 0. The white region indicates the disordered phase, B (red vertical shading) corresponds to the phase where the bistability order parameter undergoes spontaneous symmetry breaking, AF (blue horizontal shading) denotes antiferromagnetic ordering, and B + AF (purple square shading) corresponds to the phase where both order parameters are nonzero. The solid black lines denote second-order phase transitions. The NEFP phase diagram exhibits logarithmic spirals in the phase boundaries. The other NEFP is described by an analogous diagram upon switching the roles of the two order parameters (B ↔ AF, r 1 ↔ r 2 ).

FIG. 5.
Mean-field dynamics near the NEFP within the doubly ordered phase with |M 1 | = |M 2 |. The arrows denote how the fields ϕ i evolve in time, with four possible steady states. At each steady state, there is a dissipative relaxation process as well as a "coherent" rotation, resulting in a spiraling relaxation to the steady state. Two of the steady states spiral clockwise, while two spiral counterclockwise.

FIG. 6.
Dynamics of gapped (fast or massive) and soft (slow or massless) modes with arrows indicating the linear (in ϕ i , ϕ i ′) relaxation of the two modes. Near the critical point, the gapped mode quickly decays to a straight line defined by the slow direction of the soft mode.
(a) Relaxation of the field ψ B with the soft and gapped modes lying along the angles π/3 and 0, respectively. (b) Relaxation of the field ψ AF for Δ c = 1/3 with the soft and gapped modes lying along the angles -π/6 and π/6, respectively. (c) Relaxation of the field ψ AF for Δ c = 2/3 with the soft and gapped modes lying along the angles 0 and π/3, respectively. (We have adopted units where Δ + J = 1.)

FIG. 7.
Phase diagram for (one of the two) nonequilibrium tetracritical points as a function of effective mass terms and magnetic field. The transparent boundary indicates the location of an antiferromagnetic phase transition from a uniform phase. The horizontal (red) surface denotes a surface of first-order phase transitions. In both AF and uniform phases, the latter indicates a transition from a low-to a high-population phase; the sublattice population difference changes continuously across this surface. The (black) square boundary in the middle indicates the plane of h = 0. The diagram for the second nonequilibrium tetracritical point will spiral in the opposite direction.