Non-equilibrium fixed points of coupled Ising models
Jeremy T. Young, Alexey V. Gorshkov, Michael Foss-Feig, Mohammad F., Maghrebi
TL;DR
This paper uncovers non-equilibrium fixed points in coupled Ising models at multicritical points, revealing exotic phenomena like discrete scale invariance, complex critical exponents, and violation of fluctuation-dissipation relations, with potential experimental realization.
Contribution
It introduces a new class of non-equilibrium fixed points in coupled Ising models at multicriticality, characterized by complex exponents and discrete scale invariance, using a dynamical renormalization-group approach.
Findings
Emergence of non-equilibrium fixed points governing critical behavior.
Complex-valued critical exponents and spiraling phase boundaries.
Violation of fluctuation-dissipation relation at all scales.
Abstract
Driven-dissipative systems are expected to give rise to non-equilibrium 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 non-equilibrium origin. In this paper, we show that multicritical points in such systems lead to a rich and genuinely non-equilibrium 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) symmetry. They, however, coalesce at a multicritical point, giving rise to a non-equilibrium model of coupled…
| Fixed Point | |||||||
| NEFP | |||||||
| Biconical | |||||||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Non-equilibrium fixed points of coupled Ising models
Jeremy T. Young
Joint Quantum Institute, NIST/University of Maryland, College Park, Maryland 20742 USA
Alexey V. Gorshkov
Joint Quantum Institute, NIST/University of Maryland, College Park, Maryland 20742 USA
Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, Maryland 20742 USA
Michael Foss-Feig
United States Army Research Laboratory, Adelphi, Maryland 20783, USA
Mohammad F. Maghrebi
Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824 USA
Abstract
Driven-dissipative systems are expected to give rise to non-equilibrium 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 non-equilibrium origin. In this paper, we show that multicritical points in such systems lead to a rich and genuinely non-equilibrium 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) symmetry. They, however, coalesce at a multicritical point, giving rise to a non-equilibrium model of coupled Ising-like order parameters described by a symmetry. Using a dynamical renormalization-group approach, we show that a pair of non-equilibrium 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 non-equilibrium 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 non-equilibrium 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 non-equilibrium physics Diehl et al. (2008). They are very naturally realized by a variety of emerging quantum simulation platforms ranging from exciton-polariton fluids Deng et al. (2002); Kasprzak et al. (2006); Szymańska et al. (2006); Byrnes et al. (2014); Rodriguez et al. (2016, 2017), to trapped ions Bohnet et al. (2016); Schindler et al. (2013), to Rydberg gases Peyronel et al. (2012); Firstenberg et al. (2013); Carr et al. (2013); Malossi et al. (2014), to circuit-QED platforms Houck et al. (2012); Fitzpatrick et al. (2017). At long times, these systems approach a non-equilibrium 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 non-equilibrium systems has sparked a flurry of theoretical and experimental investigations.
Given their non-equilibrium 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 Mitra et al. (2006); Wouters and Carusotto (2006). 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 non-conserved order parameter in or near equilibrium—have emerged in a variety of driven-dissipative phase transitions; these include bosonic or photonic Bose-Hubbard systems Le Boité et al. (2013); Weimer (2015); Foss-Feig et al. (2017); Vicentini et al. (2018), various driven-dissipative spin models near an Ising Carr et al. (2013); Overbeck et al. (2017); Marcuzzi et al. (2014); Maghrebi and Gorshkov (2016); Owen et al. (2018); Chan et al. (2015); Wilson et al. (2016), antiferromagnetic Lee et al. (2011, 2013); Jin et al. (2013, 2014); Hoening et al. (2014); Maghrebi and Gorshkov (2016); Owen et al. (2018); Wilson et al. (2016); Chan et al. (2015), or limit-cycle phase Lee et al. (2011); Owen et al. (2018); Chan et al. (2015); Wilson et al. (2016), as well as driven-dissipative condensates consisting of polaritons Sieberer et al. (2013); Täuber and Diehl (2014). A possible exception is a two-dimensional driven-dissipative condensate, where it has been argued that the non-equilibrium Kardar-Parisi-Zhang universality class governs the long-wavelength dynamics Altman et al. (2015); Sieberer et al. (2016). While existing experiments are consistent with an effective thermal behavior Roumpos et al. (2012); Nitsche et al. (2014), 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 non-equilibrium universality classes. In particular, it has proved difficult to identify non-equilibrium universal behavior which is genuinely of a quantum nature; see Refs. Marino and Diehl (2016a); Dalla Torre et al. (2010); Cheung et al. (2018) for recent proposals to achieve non-equilibrium quantum criticality and Ref. Rota et al. (2019) 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 non-conserved order parameter is argued to be stable against all dynamical, non-equilibrium perturbations Bassler and Schmittmann (1994). More generally, the universal dynamics of various models in the Hohenberg-Halperin classification Hohenberg and Halperin (1977) are shown to be robust against a class of non-equilibrium perturbations which violate detailed balance Garrido and Marro (1989); Wang and Lebowitz (1988); Marques (1989, 1990); Tome et al. (1991); de Oliveira et al. (1993); Achahbar et al. (1996); Godoy and Figueiredo (2002); Risler et al. (2004, 2005); truly non-equilibrium behavior emerges under more constrained dynamics, for example, in the presence of a conserved order parameter in an anisotropic medium Grinstein et al. (1990); Garrido et al. (1990); Cheng et al. (1991); Täuber and Rácz (1997); Täuber et al. (2002); Katz et al. (1984); Zia (2010). In much of the previous work, situations have been considered where the phase transition is governed by a single order parameter. Due to 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 driven-dissipative 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 non-equilibrium 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 non-equilibrium universal behavior emerges at the multicritical point, giving rise to exotic critical behavior and dynamics. Our proposal to observe non-equilibrium 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 Jin et al. (2013, 2014); Kounalakis et al. (2018).
In order to determine the critical behavior, we will employ the Keldysh-Schwinger and functional integral formalism suited for the non-equilibrium setting of driven-dissipative systems Dalla Torre et al. (2010); Gopalakrishnan et al. (2010); Dalla Torre et al. (2012, 2013); Sieberer et al. (2013, 2014); Täuber and Diehl (2014); Altman et al. (2015); Sieberer et al. (2016); Maghrebi and Gorshkov (2016); Marino and Diehl (2016a); Foss-Feig et al. (2017); Maghrebi (2017). 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 Täuber (2014).
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 non-equilibrium critical properties emerging in our driven-dissipative 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 non-equilibrium 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 Appendices, 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 which 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 shall 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 -dimensional cubic lattice. The coherent dynamics of the model is governed by the Hamiltonian
[TABLE]
where is the detuning of the drive ( is the frequency of the bosons and the drive frequency), the drive strength, the hopping strength, and the strength of the nearest-neighbor interactions. The incoherent dynamics is due to loss of bosons, characterized by the Lindblad operators , where defines the loss rate. The (mixed) state of the system evolves under the quantum master equation
[TABLE]
until it approaches a non-equilibrium steady state at long times where . 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 density-dependent 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 Lee et al. (2011); Jin et al. (2013, 2014); Chan et al. (2015); Maghrebi and Gorshkov (2016); Foss-Feig et al. (2017); Marcuzzi et al. (2014). 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/high density in bistability and sublattice a/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 with bistability and 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
[TABLE]
with characterizing the stiffness, the interaction strength, the distance from the critical point (which will shift once fluctuations are taken into account), and 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 Täuber (2014)
[TABLE]
where is a “friction” coefficient and describes Gaussian white noise with correlations
[TABLE]
where is the effective temperature of the system. Near equilibrium, the “friction” coefficients 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 Täuber (2014). In the non-equilibrium 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 non-equilibrium 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 Täuber and Rácz (1997), although notable examples exist where this is not the case Grinstein et al. (1990); Garrido et al. (1990); Cheng et al. (1991); Täuber and Rácz (1997); Täuber et al. (2002); Katz et al. (1984); Zia (2010). Notice that, with the appropriate scaling of the fields , the effective temperatures can be set to . Therefore, as long as the two order parameters are not coupled, their distribution in the steady state is given by .
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
[TABLE]
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 shall again exploit our freedom in scaling the fields to set . 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. This only occurs when , leading to the effective Hamiltonian
[TABLE]
in which case the steady-state distribution is given by . However, this will not generally 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 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 , . 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 (), a pair of two NEFPs emerge where . (Both equilibrium fixed points and NEFPs are shown in Fig. 2.) Furthermore, we will 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 non-equilibrium 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.
II.1 Scaling phenomena
In the vicinity of an 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
[TABLE]
where denotes the Fourier transform in both space () and time () with the momentum and the frequency, the curly brackets denote the anti-commutator, and the subscript indicates the connected part of the correlation function. Here, is the distance from the multicritical point, is a -periodic function, and the functions and are scaling functions. While in principle the scaling behavior could be different for the two order parameters ( and ), 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 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 . Typically, it is the exponent , associated with the scaling behavior of and , that describes the scaling of the correlation length. However, the latter exponent becomes complex-valued at the NEFPs, , 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: .
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, due to the “log-periodic” function in the scaling functions, the correlations are self-similar upon the rescaling and for a particular scaling factor
[TABLE]
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 Sornette (1998). 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 that characterize the distance from the critical point, the phase boundaries themselves also exhibit a form of discrete scale invariance in ; see Fig. 4 and the discussion in the next subsection titled “Phase diagram”.
The critical exponents and characterize the anomalous dimensions corresponding to fluctuations and dissipation, respectively. In an equilibrium setting, the fluctuation-dissipation theorem dictates that the correlation and response functions are related as Täuber (2014)
[TABLE]
(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 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 on the rhs of the above equation) set by the critical exponent . 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 non-equilibrium 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
[TABLE]
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 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 at long wavelengths and fixed . 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 defined above. Finally, the values of the critical exponents at the NEFPs are provided to the lowest nontrivial order in in Fig. 3(d).
II.2 Phase diagram
The critical point described by the new fixed points is a tetracritical point. In the vicinity of the tetracritical point (with ), 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 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 symmetry associated with the bistability transition () is only an emergent one (when ), the full phase diagram (including ) can be better described as a three-dimensional 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. Fisher and Nelson (1974); Nelson et al. (1974); Bruce and Aharony (1975); Kosterlitz et al. (1976); Folk et al. (2008a, b); Eichhorn et al. (2013).
II.3 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 non-equilibrium setting of our model, the dynamics is described by the Liouvillian via [cf. Eq. (2)]
[TABLE]
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 time scale 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 define the nonzero order parameters and redefine the fields as . We then find the linearized equations of motion as
[TABLE]
where, at the NEFPs, and while we can choose ; noise, gradient, and higher-order terms have been dropped. Due to the opposite signs of and 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 and . We find that when , the angle of this complex gap relative to the real line achieves its maximum value of . This 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,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 will 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 non-equilibrium model to a model of coupled Ising-like order parameters with a symmetry, corresponding to the sublattice symmetry as well as the emergent Ising symmetry due to bistability.
III.1 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 and on neighboring sites, we have Rokhsar and Kotliar (1991); Diehl et al. (2010). 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/b and assume each to be described by a single coherent state.
Following these assumptions and using the fact that for an arbitrary operator , the resulting mean field equations of motion are given by
[TABLE]
where corresponds to the coherent state on sublattice and is the coordination number; from here on, we absorb in the microscopic parameters via and . 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 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) Jaksch et al. (2000); Lukin et al. (2001); Urban et al. (2009) or facilitate (antiblockade) Ates et al. (2007); Amthor et al. (2010) a Rydberg excitation on a neighboring site by shifting it away or towards the effective resonance.
Setting , one can immediately see that the mean-field equations become identical to those describing bistability; cf. Ref. Foss-Feig et al. (2017) 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 where behaves like an effective detuning while the product of the interaction strength and the non-interacting steady-state population [] gives the interaction-induced shift of the detuning. For , 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 mean field 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 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
[TABLE]
The field captures the effects of bistability while describes the antiferromagnetic ordering. The mean-field equations can in turn be cast in terms of these fields as
[TABLE]
At the multicritical point, and the equation governing the steady-state value of is no different than if we had not considered antiferromagnetic ordering. Thus, the critical values of , , as well as the steady-state value of are determined according to the critical point associated with bistability only. This leaves a single free parameter in the equation of motion for : . By properly tuning the latter parameter, the antiferromagnetic phase boundary can be manipulated so that it intersects the critical point associated with bistability. Working in units of , two multicritical points occur at
[TABLE]
as well as as the steady-state value of at the critical point (by virtue of symmetry, there).
The two fields 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 , while the other components remain massive and are therefore noncritical (or fast). We then adiabatically eliminate the two noncritical modes by setting and solving for in terms of . Upon substituting our solutions for the massive fields into , we find an effective description in terms of the soft modes. We closely follow Refs. Maghrebi and Gorshkov (2016); Foss-Feig et al. (2017) to identify these modes. For the bistability order parameter, we can identify
[TABLE]
with the real fields and characterizing the slow and fast modes, respectively. A similar identification has been made in Refs. Maghrebi and Gorshkov (2016); Foss-Feig et al. (2017); see also Refs. Drummond and Dechoum (2005); Dechoum et al. (2016) for a similar reasoning, although the slow and fast modes identified there make an angle of . For the antiferromagnetic field, the massless and massive components depend on the choice of the multicritical point in Eq. (17) as
[TABLE]
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
[TABLE]
with Gaussian noise
[TABLE]
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 non-equilibrium dynamics. Finally, the details of the adiabatic elimination together with the explicit values of all the coefficients (, s, s, s, s, s, and s) in terms of microscopic parameters of the model are provided in Appendix A.
It turns out that, at the level of mean field analysis, 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. Harlim and Langford (2007). However, since 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. Jin et al. (2013, 2014).
III.2 Non-equilibrium 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 which 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. Wegner (1974) 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 symmetry:
[TABLE]
Shifting the field by a constant as , 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
[TABLE]
By choosing the value of properly, the term can be dropped from the Hamiltonian while shifting the coefficients of the terms and Cardy (2000).
Similar to the above example, we should first identify the redundant operators in the non-equilibrium setting of the two coupled scalar fields and . In this case, we allow for a more general, nonlinear transformation which 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 is given by ; 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
[TABLE]
with Gaussian noise
[TABLE]
Therefore, the dynamics exhibits a symmetry when , corresponding to the emergent symmetry in addition to the sublattice symmetry . Such emergent symmetry has previously been identified in the bistability transition Marcuzzi et al. (2014); Foss-Feig et al. (2017); Maghrebi and Gorshkov (2016); 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.
IV.1 RG equations
The Langevin-type equations can be cast in terms of the response-function formalism. This allows us to study our non-equilibrium model by extending the standard techniques of the RG analysis to a dynamical setting; see, for example, Ref. Täuber (2014) for more details. The non-equilibrium partition function is defined by , where the functional integral measure and the “action” involve both fields with and their corresponding “response” fields . In the language of Keldysh field theory, corresponds to the classical field while corresponds to the quantum field. The statistical weight of can be obtained by integrating out both response fields as . While the partition function 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,25), we write the action as the sum of quadratic and nonlinear (beyond quadratic) terms
[TABLE]
Our goal is to determine the RG flow of various parameters in the action and, specifically, of the coefficients of the interaction terms.
We begin by considering the subspace defined by when either or . This subspace is special in that it is closed under renormalization to all orders. The reason is that when or , one of the two fields is not affected by the other at the microscopic level, a property that should hold at all scales. This can also be understood perturbatively in a diagrammatic scheme: If, say, , then all diagrams that could generate involve a causality violation; hence it should remain zero to all orders. An important consequence of this fact is that the relative sign of and 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 and maps , , and . Exploiting this freedom, we can set the rescaled value of to be identical to up to a sign. In doing so, we have effectively shifted the renormalization of onto , simplifying the RG analysis later. Note, however, since is rescaled by a factor , this transformation cannot change the relative sign of and . This is indeed consistent with the closure of the subspace discussed above. Having rescaled the fields appropriately, we can write the action as (the quadratic action is repeated for completeness)
[TABLE]
where indicates the relative sign of and and the coefficients , , and 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 , in which case the action can be written in a suggestive form as
[TABLE]
where the function is given by
[TABLE]
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 fluctuation-dissipation 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 .) One then should resort to an RG analysis to determine whether or not effective equilibrium is restored at long wavelengths, that is, if under RG. We shall see shortly that equilibrium proves to be a robust fixed point even when at the microscopic level.
In contrast, a Hamiltonian dynamics [similar to Eqs. (28,29)] is not possible when since a term proportional to 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 , with the exception of a decoupled fixed point where (or ). Indeed, we shall argue that a pair of genuinely non-equilibrium fixed points emerge in this case.
At a technical level, an RG analysis would be complicated as we need to consider diagrams up to two loops. This is 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). This is while the interaction terms (, , and ) 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 and , the ratio would not be renormalized at one loop.
To perform the RG analysis, we first define the renormalized parameters as
[TABLE]
where is a geometrical factor, is the Euler’s Gamma function, is an arbitrary small momentum scale (compared to the lattice spacing), and defines the small parameter of the epsilon expansion. The effect of renormalization is captured in the 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 and occur at one loop (), while those of appear at two loops (). These perturbative corrections, while having some similarities with their equilibrium counterparts, are more complicated due to their non-equilibrium nature.
Using the above factors, we determine the RG flow and beta functions via
[TABLE]
where and . 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 . At any such fixed point, the scaling behavior of the remaining parameters is governed by power laws whose exponents depend on . Here, we report the beta functions for the interaction parameters (the details are provided in Appendix C):
[TABLE]
where we have introduced . These equations exhibit a number of important features. First, for , we can absorb a factor of into , leaving the two beta functions for independent of . We thus immediately recover a pair of uncoupled equilibrium Ising phase transitions, as one should expect. Second, under equilibrium conditions where and , we recover the standard beta functions in equilibrium. In a similar fashion, we can absorb the factors of into and into , 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 non-equilibrium model, statics and dynamics are inherently intertwined. Indeed, no redefinition of the coupling terms can lead to beta functions that would be independent of and . This is not the case for as they can always be absorbed in other parameters; for example, we can still absorb into and into in the beta functions. This reflects the fact that, through an appropriate rescaling of the fields, one can always rescale arbitrarily without changing , or the overall structure of the action.
To set up the full RG equations, let us define the parameters
[TABLE]
With these definitions, the beta functions for the new interaction parameters depend only on the renormalized parameters , and . 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 () 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, this leaves a total of five parameters needed to define the fixed point. The remaining parameters () 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 and , we use the identity
[TABLE]
We now report the full set of beta functions of the parameters of our model (with and set to zero at the fixed point)
[TABLE]
where we have defined and the functions
[TABLE]
The functions always appear in the RG equations in pairs, with one taking and the other 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 only and the other from those that involve only. Similarly, under the mapping , and , the beta functions are left unchanged. This reflects the fact that we can switch the role of and without changing the physics. As a result, if either , , or 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 non-equilibrium model, but it is instructive to first consider the RG equations under equilibrium conditions where the temperatures are equal, i.e., , and . We then immediately find that the temperature ratio does not flow, ; hence the two temperatures remain identical at all scales. Furthermore, the temperature itself—and not just the ratio—remains scale invariant (), indicating (effective) thermal equilibrium. Finally, as remarked earlier, the RG equations for the interaction terms become independent of 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 . The first scenario is that the beta function vanishes when . Since the dynamical critical exponents are related to the flow of as , we find that under this scenario. This means that both fields are governed by the same dynamical critical exponent, giving rise to a “strong dynamic scaling”. The second scenario occurs when , which would lead to the fixed point or depending on the sign of . This behavior is then described by a “weak dynamic scaling” where the two fields exhibit different dynamical scaling properties and exponents Halperin and Hohenberg (1969); De Dominicis and Peliti (1977); Dohm and Janssen (1977); see also Täuber (2014). Similarly, one can consider the beta function 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 , the RG flow leads to either or , which both correspond to the subspace. However, this subspace does not appear to be amenable to perturbative RG. In this sector, we find that 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 indicates that an alternative approach from our present analysis should be considered in this scenario. In this work, we shall restrict ourselves to the case where and are both finite and nonzero.
IV.2 Fixed points of RG flow
With the RG beta functions, we can now identify the resultant fixed points. In the sector, the only fixed points of the RG equations are those where , exhibiting a strong dynamic scaling, as well as , indicating that the two temperatures become identical at the fixed point. Indeed, aside from the case of = 0, the only possible fixed-point value of at this order is 1. This can be seen by noting that the only other root of Eq. (35d) is , which is always negative and thus unphysical. Similarly, noting that the beta functions for are identical at this order, all coupled fixed points in this sector will satisfy . In light of this, we immediately identify as the only possible solution of Eq. (35e). Remarkably, an effective equilibrium behavior emerges in this sector despite the underlying non-equilibrium nature of the dynamics. In particular, we recover the familiar equilibrium 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 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 is truly the result of a nontrivial two-loop RG analysis.
In the sector, any nontrivial fixed point is truly non-equilibrium 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 non-equilibrium fixed points as
[TABLE]
These fixed points also exhibit a strong dynamic scaling since , so the two fields are governed by the same dynamical scaling. Furthermore, we find , 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, 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 shall 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
[TABLE]
where 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 lowest-order corrections to all five parameters, we can only determine these eigenvalues up to . This is because in order to fully determine to , we need to consider the two-loop corrections to the coupling terms and the three-loop corrections to . However, when a fixed point possesses a higher symmetry than the underlying field theory, then it is possible to determine some of the eigenvalues without including higher-order corrections. This 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 when ), a similar situation occurs with respect to statics and dynamics, where perturbations in the dynamics () cannot affect the behavior of the statics (). This makes it possible to inspect the stability of up to at the same order of the RG calculations. In the full non-equilibrium model, the statics is not decoupled from the dynamics, so the stability in 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 non-equilibrium terms. As a result, it is possible to determine the stability in for the two coupled equilibrium fixed points, and we find both to be stable in . However, for all of the coupled fixed points, remains marginal. In short, to the lowest order in our perturbative expansion, the system flows to the NEFP (equilibrium fixed point) in the () sector. While, in principle, non-perturbative 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 couplings.
Finally, we remark that in the case of the original microscopic model, and have opposite signs, which thus carries over to the relative sign of and . Thus it is plausible to expect a critical behavior governed by the NEFPs.
IV.3 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 , 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
[TABLE]
where are general scaling functions. We have dropped the subscript from 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 for 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
[TABLE]
The renormalization of the parameters and the corresponding exponent requires a special treatment and will be discussed later in this section. At the non-equilibrium critical point, we find [cf. Eqs. (31a,37) together with the factors in Appendix C]
[TABLE]
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 . Using Eq. (40), the critical exponents at the NEFPs are given by
[TABLE]
While, in equilibrium, as a consequence of the fluctuation-dissipation theorem, we have since the temperature itself is scale dependent, , at the NEFP. Note also that the critical exponents 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 and . 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 to the first nontrivial order in the epsilon expansion. This is in contrast with equilibrium where , a fact that can eve ben proved on general grounds (e.g., unitarity in a related quantum field theory) Zinn-Justin (1996). If this feature () extends beyond perturbation theory to, say, two dimensions, it would indicate that the correlation function () 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 is related to in an identical fashion as in equilibrium.
Next we consider the renormalization of the mass terms. This requires special care as their renormalization is intertwined. Similar to our redefinition of , we should instead consider the renormalization of 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 . Defining and the flowing parameters with , we find the flow equations [cf. Eqs. (31a,37) together with the factors in Appendix C]
[TABLE]
where the refer to the two NEFPs with opposite signs of . The flow equations can be solved as
[TABLE]
where
[TABLE]
These equations can be cast in a more compact notation as
[TABLE]
Defining , we can recast this equation as where the critical exponent emerges as
[TABLE]
Interestingly, the exponent becomes complex-valued at the NEFP. We can then express the scaling functions in Eq. (39) as
[TABLE]
where while denotes the polar angle in the - plane. Additionally, is a -periodic function. To obtain these equations, we have used the fact that can instead be written as a function of and . 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 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 Sornette (1998). Rather than a self-similar behavior at all length scales, a preferred scaling factor emerges as
[TABLE]
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 Feigenbaum and Freund (1996), earthquakes Sornette and Sammis (1995), equilibrium models on fractals Karevski and Turban (1996), and several other systems Sornette (1998). Log-periodic functions and the emergence of a preferred scale have been identified in the early developments of renormalization group theory Jona-Lasinio (1975); Nauenberg (1975); Niemeijer and Van Leeuwen (1976), but they have been rejected as artifacts of position-space 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 Sornette (1998). This phenomenon has also emerged in recent works in the context of driven-dissipative quantum criticality Marino and Diehl (2016a, b) as well as the dynamics of strongly interacting non-equilibrium systems Maki et al. (2018). A particularly well-known example of RG limit cycles is the behavior of Efimov states Braaten and Hammer (2006); Efimov (1970), 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 LeClair et al. (2003); Kaplan et al. (2009). Disordered systems provide another context where complex-valued exponents and discrete scale invariance have been noted in both classical Aharony (1975); Chen and Lubensky (1977); Khmelnitskii (1978); Weinrib and Halperin (1983); Boyanovsky and Cardy (1982) and quantum Hartnoll et al. (2016) settings. The discrete scale invariance reported in this work, however, appears to be unique as it has emerged in an effectively classical yet non-equilibrium model in the absence of disorder.
The discrete scale invariance approaches a continuous one as the upper critical dimension, , is approached. In three dimensions, perturbative values at the NEFP (with ) yield a very large scaling factor (); 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 non-equilibrium systems. Additionally, higher harmonics in the periodic function 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 symmetry Chan et al. (2015). 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 find that in the doubly-ordered phase, the Liouvillian gap becomes complex-valued (see Sec. IV.6). This 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 1, we summarize all of the fixed points (aside from those involving the trivial Gaussian fixed point) and their critical exponents to the lowest order.
IV.4 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, , four different phases emerge: A disordered phase with ; two phases with either corresponding to bistability or leading to antiferromagnetic ordering; and, finally, a doubly-ordered phase where both fields become ordered, . 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 term is generated.
The phase boundaries are governed by the scaling behavior of . Let us set the effective magnetic field to zero, , 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 . In this case, the scaling functions are solely determined as a -periodic function of ; this is achieved by eliminating the momentum scale in Eq. (48) in favor of . 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 ). Therefore, the shape of the phase boundary is given by
[TABLE]
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 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 , leading to a surface of first-order phase transitions where undergoes spontaneous symmetry breaking. This first-order transition occurs in both the uniform phase (defined by the B phase at ) and the antiferromagnetic phase (defined by the B+AF phase at ). 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 , where and describe the scaling behavior of the order parameter with and , respectively, within the ordered phase. The effect of the magnetic field is to “unravel” the spirals for small . The corresponding phase diagram for nonzero effective magnetic field is illustrated in Fig. 7.
IV.5 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,
[TABLE]
In the standard equilibrium setting, these exponents are related to via
[TABLE]
To see how the above hyperscaling relations are modified near the NEFPs, first note
[TABLE]
and that in the ordered phases, and will have the same scaling behavior. Similar to our phase-diagram analysis, we take the limit or and allow to define the momentum scale in place of . This amounts to replacing the prefactors in Eq. (39) with . We then find the scaling of the order parameter and magnetic susceptibility to be
[TABLE]
where are -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 which 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
[TABLE]
where only the real part enters.
IV.6 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 ; for notational convenience, we make the change of variables where the fields 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 and terms too)
[TABLE]
In addition, several cubic terms are also introduced, which are not reported for simplicity. We set the vertices and to zero since, by definition, solely represent the fluctuations. This in turns sets and . Upon including the effect of fluctuations to , the remaining quadratic vertices are then given by
[TABLE]
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 factors.
Putting these terms together with the quadratic part of the action, we find
[TABLE]
where and . The poles of the corresponding propagators are then obtained as
[TABLE]
or, more explicitly, as
[TABLE]
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
[TABLE]
Indeed, in equilibrium, where , 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 , the above condition can be satisfied. To see this, let us cast the above condition for in terms of and as
[TABLE]
Recalling that at the NEFPs, at least to the lowest order in the epsilon expansion, the above condition is trivially satisfied whenever . In this case, the pole with the lowest nonzero decay rate takes the form (with and )
[TABLE]
In fact, with , 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 , which corresponds to the Liouvillian gap closing at the angle 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 model of -component vector-like order parameters leads to similar behavior. In fact, we find that, for and , the corresponding Liouvillian gap closes at the angles and , respectively. We should then conclude that different non-equilibrium 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 follow-up 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 time scale which itself is associated with the exponent . Thus the gap scales as , where the exponent is due to the scaling of itself. With the order parameters and 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 . 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 Blais et al. (2007); Houck et al. (2012); Devoret and Schoelkopf (2013). 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 Hartmann et al. (2008); Tomadin and Fazio (2010); Carusotto and Ciuti (2013); Peropadre et al. (2013). 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 driven-dissipative models in cavity- and circuit-QED systems Dalla Torre et al. (2013); Jin et al. (2013, 2014); Wilson et al. (2016); Biondi et al. (2017); Iemini et al. (2018). Recent experiments have even identified a driven-dissipative many-body phase transition Fitzpatrick et al. (2017).
For the model considered in this work (see Sec. II), many-body experimental platforms already exist that include drive, 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 Kumar and DiVincenzo (2010); Hoi et al. (2013); Holland et al. (2015). A more general version of our model has been considered in Refs. Jin et al. (2013, 2014), along with a discussion on how the nonlinear interaction terms can be tuned experimentally via Josephson nanocircuits. 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 Kounalakis et al. (2018).
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 Jin et al. (2013, 2014), 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. Jin et al. (2013, 2014); Kounalakis et al. (2018). 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 Li and Koch (2017). For example, by sending and on one of the two sublattices, the sign of can be changed while leaving the remaining terms fixed. Similarly, one can also map by taking the complex conjugate of the master equation, which, together with the previous mapping, allows an appropriate choice for the sign of . 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 field-theoretical approach—appropriate in the vicinity of the phase transition—we have shown that the critical behavior at this point can be mapped to a non-equilibrium stochastic model described by a symmetry. Using perturbative renormalization group techniques, we have determined the RG flow equations of the model and identified a pair of new classical non-equilibrium 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 non-equilibrium behavior. This perspective opens a new avenue to investigate and experimentally realize non-equilibrium phases and phase transitions in the context of driven-dissipative systems without relying on the engineering of complicated nonlocal or non-Markovian dissipation.
Future experimental and numerical studies into the NEFPs discussed in this work are crucial to develop a more complete understanding of their properties. Characterizing the discrete scale invariance, either in the dynamics or the form of the phase boundary, defines a particularly important direction. Investigating the possible emergence and the critical behavior of such non-equilibrium phase transitions in low dimensions is worthwhile. It would be particularly interesting to identify low-dimensional ordering and phase transitions that are not otherwise possible in equilibrium settings. Another question that remains open is the fate of the subspace , namely if it contains new NEFPs. Beyond these non-equilibrium generalizations of model A systems, one can further consider similar non-equilibrium versions of other equilibrium universality classes. While we focused on the particular case of an experimentally relevant model with only two scalar order parameters, our analysis indicates that a large class of new non-equilibrium multicritical points are yet to be discovered. A natural extension of our work is to identify possibly new NEFPs in models involving vector-like order parameters Fisher and Nelson (1974); Nelson et al. (1974); Bruce and Aharony (1975); Kosterlitz et al. (1976); Folk et al. (2008a, b); Eichhorn et al. (2013). While a driven-dissipative condensate of polaritons has been investigated theoretically in detail Sieberer et al. (2013); Täuber and Diehl (2014), recent experimental studies into condensate supersolids Baumann et al. (2010); Mottl et al. (2012); Landig et al. (2016); Léonard et al. (2017a, b) can provide excellent platforms for probing any emergent NEFPs. In addition to the symmetry of the condensate, the two coupled optical cavities can provide either an additional symmetry (corresponding to a lattice supersolid) or an approximate symmetry (corresponding to a continuous supersolid).
Acknowledgements.
We thank M. Kardar, S. Diehl, R. Fazio, J. Marino, I. Boettcher, S. Mathey, R. Matthew, S.-K. Chu, J. Curtis, and N. Grabon for helpful discussions. J. T. Y. and A. V. G. acknowledge support from the NSF PFC at JQI, the DOE BES QIS program (Award No. DESC0019449), NSF PFCQC program, AFOSR, ARL CDQI, ARO MURI, and the DOE ASCR Quantum Testbed Pathfinder program (Award No. DE-SC0019040). M. M. acknowledges support from NSF under Grant No. DMR1912799 and start-up funding from Michigan State University.
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. Foss-Feig et al. (2017). We begin by constructing the Keldysh path integral, then identify a semi-classical 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 shall return to the latter symmetry once we identify the semi-classical limit and corresponding Langevin equations. We cast our model in terms of a Keldysh path integral as
[TABLE]
where the action is defined as
[TABLE]
with the classical/quantum fields and the normal-ordered form of the Hamiltonian. The third line corresponds to the particular case of the Lindblad operator , although this approach may easily be extended to more general Lindbladians Sieberer et al. (2016). 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
[TABLE]
where is the coordination number. This expression bears a close resemblance to the action of Eq. (12) in Ref. Foss-Feig et al. (2017), but they differ in the form of their interactions (which involve gradient terms here) and due to our use of normal ordering rather than the Weyl ordering of Ref. Foss-Feig et al. (2017). This motivates a similar rescaling of the parameters as
[TABLE]
The parameter effectively describes a density scale for the microscopic model via , where is . Since varying the density scale also modifies the interaction energy per particle, the interaction strength should be reduced correspondingly such that ; similarly, the drive should be increased so that . We can then rewrite the action as
[TABLE]
In the limit of large , the last term (the second term in the last line) can be dropped, leading to an action that is at most quadratic in . This is simply because a large population corresponds to the semi-classical limit represented by a large classical field and small fluctuations due to the quantum field . Using this fact, we can map the action to a Langevin equation as Täuber (2014)
[TABLE]
Note that the noise level is further suppressed at larger as should be expected from our semi-classical treatment.
Next we include the sublattice symmetry by defining as the sublattice average and as the sublattice difference of the field ; see Eq. (15), which differs by a factor of due to our semi-classical limit. Our new Langevin equations are now
[TABLE]
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 .
Our model exhibits two multicritical points where two modes (each a component of one of the two fields ) become critical. Due to the sublattice symmetry, at the multicritical points up to fluctuations. Working in units where , the two multicritical points occur at
[TABLE]
with .
Next, we expand the Langevin equations in the vicinity of the two multicritical points as
[TABLE]
The soft and gapped modes can be determined as linear combinations of the real and imaginary parts of the two fields . As in the case of bistability, the first pair is identified as Foss-Feig et al. (2017)
[TABLE]
where is massive and relaxes quickly while defines the slow field. Identifying the massive and massless components of the field depends on the corresponding multicritical point as
[TABLE]
where again the primed (unprimed) field indicates the massive (massless) field. Note that these differ from the main text by a factor of , which is done to simplify the resulting parameters in the Langevin equations by moving all the dependence to the noise term.
Upon adiabatically eliminating the massive fields and restoring the gradient terms, we arrive at the Langevin equations
[TABLE]
with Gaussian noise
[TABLE]
The various numerical factors are summarized in Table 2 for the two multicritical points under consideration. Note that and . We can see at this point that the opposite signs of and indicate that no Hamiltonian description is possible. Indeed, as will be discussed in the following section, this will carry over to the signs of and , 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 (75). In general, this can be done at the level of the Schwinger-Keldysh action; however, we shall 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 Gardiner (2009)
[TABLE]
where is the stochastic noise that obeys the Itô rules:
[TABLE]
At the multicritical point, where the effective masses and the magnetic field are set to zero, the Langevin equations (75) can be written in the form of Eq. (77) with
[TABLE]
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 and ; the functions are general (but local) nonlinear maps that are invertible in a neighborhood around the multicritical point and preserve the sublattice symmetry . The equations governing the dynamics of the new variables take the form
[TABLE]
We have used the Itô rules to derive the above equation, which is known as Itô’s formula or Itô’s lemma Gardiner (2009). A similar stochastic equation can be derived for by switching . Note that the terms on the rhs should be expressed in terms of through the inverse functions . 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 ), 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 , 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.
We shall assume without loss of generality that
[TABLE]
rescaling the fields by a constant factor does not allow us any additional freedom while rotations obscure the symmetry . Additionally, we do not consider a constant shift in the field (a shift in is disallowed due to symmetry) for now, but discuss it separately later in this appendix. Based on the structure of Eq. (80), we notice that the quadratic terms in and 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
[TABLE]
the modification of the cubic terms to lowest order in the coefficients and interaction terms is given by
[TABLE]
where the matrix is given by
[TABLE]
Notice that the coefficient is left unchanged. The rank of matrix is 2, which then determines the number of corresponding redundant operators.
In addition to the two redundant operators above, a third one emerges due to a constant shift . Under this transformation, the quadratic terms transform as
[TABLE]
The effective mass and magnetic field terms also change, but this simply shifts the location of the critical point.
The three redundant operators derived here can be used to always set the terms to zero. This can be understood by noting that the transformation corresponding to allows one to adjust the ratios of and relative to without changing the quadratic terms. By properly using this redundancy, these ratios can be tuned until the constant shift in Eq. (85) shifts all three quadratic terms to zero. At the same time, the cubic terms transform as
[TABLE]
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. (82)], 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 and is indeed directly determined by the relative sign of and , leading to criticality described by the NEFPs. Finally, we note that the coefficient appears in the denominator of the above transformations. We assumed 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.
C.1 Diagrammatic techniques
To define the Gaussian propagators, we start with the Gaussian model with the corresponding action
[TABLE]
The Gaussian response and correlation functions are then given by
[TABLE]
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. Due to the structure of the action, one can find the corresponding factors for by switching the subscripts and multiplying by a factor of .
Finally, in order to determine the factors, we employ dimensional regularization using the minimal subtraction procedure. This means that only the ultraviolet divergences, in the form of powers of , are incorporated into the factors. For simplicity, we only present these divergences in the evaluation of integrals in the following sections; non-divergent terms are dropped. Technical details for evaluating integrals of the type presented here may be found in, e.g., Ref. Täuber (2014).
C.2 One-loop diagrams
C.2.1 Mass terms
In this section, we consider corrections to . The corrections to can be obtained by switching the roles of and . There are two diagrams that provide corrections, which are illustrated in Fig. 10(a,b). The combinatorial and interaction factors are (a) and (b) . Thus the one-loop contribution to the renormalization of is
[TABLE]
so we should evaluate the integral of . The latter can be written as
[TABLE]
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 is then determined by
[TABLE]
and similarly for ; note that the factors of in the denominator have been dropped as they introduce corrections. Defining an additive renormalized mass term , we can determine the factors for . To this end, we need to compute integrals of the form
[TABLE]
where we have included the geometrical factor . We can then compute the corresponding factors as
[TABLE]
From this point on, we simply write as .
C.2.2 Coupling terms
We first consider one-loop corrections to . We need to consider two diagrams as illustrated in Fig. 10(c,d). The combinatorial and interaction factors are (c) and (d) . Thus the one-loop contribution to the renormalization of is
[TABLE]
whose evaluation we put aside for the moment.
The renormalization of involves eight diagrams, four of which are illustrated in Fig. 10(e-h) and renormalize . The combinatorial and interaction factors are (e) , (f) , (g) , and (h) . The remaining four diagrams, corresponding to the renormalization of , can be simply obtained by interchanging the two fields. The resulting set of internal diagrams and combinatorial factors are the same up to a factor of , so and are renormalized in the same way at this order. Indeed, this reflects the fact fact that is not renormalized at the one-loop order.
Thus the one-loop contribution to the renormalization of is
[TABLE]
This expression involves a nontrivial integral of the form
[TABLE]
which can be reduced to
[TABLE]
where and . Noting that and is a finite constant, then according to the minimal subtraction procedure, , where has been set to 0 in the non-divergent part to extract the residue of the pole. Similarly, we expand to extract the pole. Thus the above integral becomes
[TABLE]
where we have included the geometrical factor . This result in combination with the diagrams considered above results in the following factors
[TABLE]
C.3 Two-loop diagrams
First, we consider the two-loop corrections that arise from the terms. There are three two-loop diagrams that renormalize and as shown in Fig. 11(a-c). The corresponding combinatorial and interaction factors are (a) , (b) , and (c) .
As in the case of the coupling terms, the internal diagrams are all of the same form, so we consider the generic internal diagram, which corresponds to the integral
[TABLE]
where correspond to the top, bottom, and middle propagators, respectively, in the diagrams. This may be reduced to the form
[TABLE]
In order to determine corrections to and , we consider and , respectively, in the limit . Additionally, noting and taking while extracting a factor , we have
[TABLE]
where we have included the geometrical factor . We evaluate both of these integrals exactly to find the resulting corrections
[TABLE]
from which we identify the following factors
[TABLE]
where the three corrections correspond to diagrams (a), (b), and (c) in Fig. 11, respectively.
Next, we consider the two-loop corrections to the terms. There are two such diagrams, which are illustrated in Fig. 11(d,e). The combinatorial and interactions factors are (d) and (e) . Note the lack of minus sign due to the sign difference in . Again, the internal diagrams are all of the same form, so we instead consider the generic internal diagram, which corresponds to the integral
[TABLE]
where correspond to the the three propagators in the diagrams. This integral can be reduced to the form
[TABLE]
We note that . Taking the limit and extracting a factor , we rewrite the latter integral as
[TABLE]
where we have included the geometrical factor . Taking advantage of the fact that at least two of the must be the same, this integral can be evaluated analytically as
[TABLE]
where corresponds to the field with one propagator and to the field with two.
Thus we identify the corresponding factors
[TABLE]
where the two corrections correspond to diagrams (d) and (e) in Fig. 11, respectively. Note that the factors are half of their combinatorial factors. This is because the zeroth-order vertex is rather than .
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
[TABLE]
where and are the interaction strengths. Additionally, we define a scaling function via , 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 in terms of the partial derivatives with respect to other parameters as
[TABLE]
Next we employ the method of characteristics and define . We then introduce the flowing dimensionless parameters , via
[TABLE]
At the fixed point, for all and the parameters and flow functions assume their fixed-point values and . This allows us to easily solve the flowing parameters as . Additionally, since all the terms in the perturbation series involve integrals of Gaussian propagators, we can reduce the scaling functions to [cf. Eq. (88)],
[TABLE]
where we have utilized the fact that the scaling behavior of is the same for both fields. Upon applying the matching condition and carrying out a similar procedure for the correlation function, we can express the response and correlation functions as
[TABLE]
where we have further simplified the arguments of the scaling functions by dropping factors of and and excluding the argument in a slight abuse of notation. Comparing these scaling functions against those in Eq. (8), we identify the critical exponents
[TABLE]
Similarly, we can identify , although the subtleties of a complex-valued exponent at the NEFPs are discussed in detail in the main text.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Diehl et al. (2008) S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller, Quantum states and phases in driven open quantum systems with cold atoms, Nat. Phys. 4 , 878 (2008) . · doi ↗
- 2Deng et al. (2002) H. Deng, G. Weihs, C. Santori, J. Bloch, and Y. Yamamoto, Condensation of Semiconductor Microcavity Exciton Polaritons, Science 298 , 199 (2002) . · doi ↗
- 3Kasprzak et al. (2006) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Bose-Einstein condensation of exciton polaritons, Nature 443 , 409 (2006) . · doi ↗
- 4Szymańska et al. (2006) M. H. Szymańska, J. Keeling, and P. B. Littlewood, Nonequilibrium Quantum Condensation in an Incoherently Pumped Dissipative System, Phys. Rev. Lett. 96 , 230602 (2006) . · doi ↗
- 5Byrnes et al. (2014) T. Byrnes, N. Y. Kim, and Y. Yamamoto, Exciton-polariton condensates, Nat. Phys. 10 , 803 (2014) . · doi ↗
- 6Rodriguez et al. (2016) S. R. K. Rodriguez, A. Amo, I. Sagnes, L. Le Gratiet, E. Galopin, A. Lemaître, and J. Bloch, Interaction-induced hopping phase in driven-dissipative coupled photonic microcavities, Nat. Commun. 7 , 11887 (2016) . · doi ↗
- 7Rodriguez et al. (2017) S. R. K. Rodriguez, W. Casteels, F. Storme, N. Carlon Zambon, I. Sagnes, L. Le Gratiet, E. Galopin, A. Lemaître, A. Amo, C. Ciuti, and J. Bloch, Probing a Dissipative Phase Transition via Dynamical Optical Hysteresis, Phys. Rev. Lett. 118 , 247402 (2017) . · doi ↗
- 8Bohnet et al. (2016) J. G. Bohnet, B. C. Sawyer, J. W. Britton, M. L. Wall, A. M. Rey, M. Foss-Feig, and J. J. Bollinger, Quantum spin dynamics and entanglement generation with hundreds of trapped ions, Science 352 , 1297 (2016) . · doi ↗