Absorbing state phase transition with competing quantum and classical fluctuations

Stochastic processes with absorbing states feature remarkable examples of non-equilibrium universal phenomena. While a broad understanding has been progressively established in the classical regime, relatively little is known about the behavior of these non-equilibrium systems in the presence of quantum fluctuations. Here we theoretically address such a scenario in an open quantum spin model which in its classical limit undergoes a directed percolation phase transition. By mapping the problem to a non-equilibrium field theory, we show that the introduction of quantum fluctuations stemming from coherent, rather than statistical, spin-flips alters the nature of the transition such that it becomes first-order. In the intermediate regime, where classical and quantum dynamics compete on equal terms, we highlight the presence of a bicritical point with universal features different from the directed percolation class in low dimension. We finally propose how this physics could be explored within gases of interacting atoms excited to Rydberg states.


I. INTRODUCTION
Non-equilibrium phenomena are ubiquitous in nature, ranging from the microscopic scales of chemical reactions to the macroscopic ones of disease-spreading. Remarkably, analogously to the equilibrium case, nonequilibrium ensembles can show the emergence of universal behavior, signaling the irrelevance of the microscopic details of the dynamics for macroscopic observables. This occurs when such out-of-equilibrium systems start to act collectively [1][2][3][4]. On a fundamental level, a distinction arises depending on the presence or absence of detailed balance [5][6][7][8], between systems which evolve towards a stationary equilibrium state (e.g., quenched systems coupled to thermal baths [9]) or that preserve their nonequilibrium character even in the long-time limit, representing flux equilibrium states. The universal dynamical features of purely classical systems have been extensively studied and classified both for unbroken [10] and broken [11][12][13][14] detailed balance, i.e., genuine non-equilibrium systems. Recently, experiments in various platforms have started to systematically probe driven open quantum systems. The spectrum includes light-driven semiconductor heterostructures [15], arrays of driven microcavities [16,17], cold atoms in optical lattices [18], cavities [19,20] and microtraps [21][22][23]. Several among these instances employ excitation of the atoms to high-lying Rydberg orbitals [24][25][26] in order to achieve strong interatomic interactions and to study cooperative effects [27][28][29]. In all these systems, the driving/dissipation introduces coherence loss and explicitly violates the equilibrium conditions at the microscopic level [7,30]. It is thus a challenge to identify to what extent the non-equilibrium and the quantum nature of the dynamics impact on the macroscopic phase diagram and phase transition properties. Oftentimes, upon coarse graining such systems lose their quantum character and equilibrium conditions are effectively restored [31][32][33][34][35][36]. But there are instances where non-equilibrium [37,38] and quantum [39,40] aspects persist even at asymptotically large wavelength.
Directed percolation (DP) represents an instance of a classical, but intrinsically non-equilibrium system (for a review, see [11]). Despite its robustness, its experimental observation has so far been elusive [41], with a single remarkable exception [42,43]. It was suggested recently to realize and explore DP dynamics in cold gases of atoms excited to high-lying Rydberg states [44]. In that case, the non-equilibrium nature of dynamics persists macroscopically, but the impact of quantum dynamics fades out completely under coarse graining.
In this work, we harness the opportunities that result from the fact that such Rydberg gases indeed represent driven open quantum systems to go beyond the realm of classical physics, and establish a novel absorbing state phase transition characterized by the interplay of classical and quantum terms on equal footing. This transition does not fall into the DP universality class, and its origin can be unambiguously traced back to the presence of coherent dynamics. More precisely, the latter introduces a strong coupling first-order non-equilibrium phase transition without counterpart in the purely classical DP problem. Remarkably, this discontinuous phase transition terminates in a novel bicritical point which even asymptotically at large distances and in dimensions d < 2, does not feature the symmetries underlying DP, or any equilibrium problem.

II. MODEL
In the following we reproduce a quantum variant of the contact process (for an introduction we refer to Ref. [11]). Its defining property is the following: in a lattice of "active" and "inactive" sites, the former can spontaneously decay to inactive, whereas activation can only occur in the proximity of already active sites. Thus, the fullyinactive state is absorbing, i.e., once reached it cannot be left. Specifically, we consider a lattice of quantum two-level systems with spacing r. On every site k we define the basis |a k (active) and |i k (inactive), the den-arXiv:1601.07305v1 [cond-mat.stat-mech] 27 Jan 2016 sity of active sites n k = |a k a k | and the ladder operators σ + k = |a k i k | and σ − k = |i k a k |. Under the ac- Figure 1. (Color online) (a) Fundamental processes. We consider a lattice whose sites admit two states: active (red) and inactive (green). Active sites decay at a rate γ and become inactive. Proliferation of active sites is possible through classical (rate κ) and quantum (strength Ω) branching. (b) Phase diagram constructed from the effective action (5) in saddlepoint approximation (color code corresponds to density of active sites). In the classical limit (Ω = 0) the system exhibits a continuous (2nd order) directed percolation phase transition between an absorbing state and one with finite density. This transition extends into the quantum regime (thick red line) until the critical point α is reached. In the quantum limit (κ = 0) a first-order transition is found which also extends into the classical regime (dashed yellow line) up to point α.
In the neighborhood of this line, a narrow region of coexistence of two attractive stationary solutions is present, which is not resolved here. The high values of the density reached in the active phase stem from neglecting higher orders in n in the action, which would otherwise prevent it from exceeding 1/2.
is the quantum Hamiltonian, σ x k = σ + k + σ − k , and "nn k" denotes nearest neighbors (nn) of site k; D[X]ρ = XρX † − (X † Xρ + ρX † X)/2 is the dissipator and L a,k are the so-called jump operators, having the index a labeling the process type, and k the lattice site. These jump operators are chosen to define a modified contact process [11], which is known to feature a DP transition, and include decay L d,k = √ γσ − k (|a k → |i k ) and -for every neighbor j of k -branching L b,j,k = √ κ n j σ + k (an active site can activate a neighboring one |a j i k → |a j a k ) and coagulation L c,j,k = √ κ n j σ − k (the inverse process |a j a k → |a j i k ). The operator Π k in H represents the simplest choice which effectively reproduces the requirement of an active site nearby to flip a spin; this makes H the "minimal quantum equivalent" of the noisy branching/coagulation above. Similar "constrained" Hamiltonians have been studied in the past with a focus on many-body localization [47,48].

III. EQUATIONS OF MOTION AND DENSITY PATH INTEGRAL
We infer here the properties of the phase diagram by exploiting an effective path integral description for the density variable n k alone. We start by deriving the Heisenberg-Langevin equations of motion (EOM) [49] for the single-site operators n k , σ x k and σ y k = −iσ + k + iσ − k . For convenience we introduce the coordination number z (number of nearest neighbors per lattice site), the shorthand Σ x/y k = σ x/y k j nn k σ x j , rescale time by t → τ = γt and the rates accordingly, i.e. χ = κ/γ and ω = Ω/γ: The quantum noise termsξ α k consider the fluctuations of the bath and depend on the structure of the jump operators. They show vanishing averages but non-trivial, Markovian correlations, which for the present setup are (in rescaled units) ξx = −δ k,k σ + k and ξn kξ y k = iδ k,k σ + k up to leading order in the density (see Appendix A).
In the following, we work in the continuum limit (k, t) → ( x, t) ≡ X and derive an effective path integral for the density field n X via a Martin-Siggia-Rose (MSR) construction [3,[50][51][52], presented in Appendix B. Crucially, the σ x,y -fields are gapped, and thus can be integrated out perturbatively. The resulting long wavelength field theory depends on the density variable n alone, and is obtained by additionally performing a derivative expansion of the action. It reads where D = r 2 χ represents a diffusion constant (lattice spacing r) and are the microscopic coupling constants. The response fieldñ encodes the linear response properties of n under small perturbations. At this point we emphasize two key properties of the action (5): First, the absence of a density independent Markovian noise level ∼ Tñ 2 X (which is necessarily present in any classical system in thermal equilibrium). This is characteristic of DP dynamics, which feature the absence of density fluctuations in the absorbing state n X = 0 and consequently a multiplicative kernel ∼ n X . An additive noise introduced by the dissipative terms L d = √ γσ − only occurs in the eliminated spin variables σ x,y . Second, the presence of a non-zero coherent coupling ω = 0 -i.e. the intrinsic quantum effect -leads to the appearance of non-zero couplings u 4 and µ 4 as well as a negative contribution to u 3 . This new "quantum" scale ω breaks a fundamental symmetry of the DP class (specified below) and strongly modifies the phase diagram compared to the purely dissipative model [see Fig. 1].

IV. EFFECTIVE POTENTIAL AND MEAN-FIELD PHASE DIAGRAM
The discussion of the various phases and phase transitions of the system is considerably simplified by realizing that the deterministic contribution to the action S is a local effective potential. In the absence of fluctuations Γ characterizes the mean-field phases, which are determined by the properties around its minima. The corresponding phase diagram is shown in Fig. 1(b). The active phase is identified by ∆ < 0, u 4 ≥ 0 and u 3 > 0, which leads to a single minimum of the effective potential at finite density. On the other hand, when both ∆ and u 4 are positive, there is a local minimum of Γ at n = 0. For negative and sufficiently strong cubic coupling u 3 < −2 √ u 4 ∆, there exists a second local minimum at finite density n > 0. In this regime, the mean field evolution features two attractive fixed points and the thermodynamic phase is determined within the optimal path approximation in phase space [53].
Three different types of phase transitions from the active to the inactive state can be thus identified, their nature depending on the specific choice of parameters and the dimensionality. When the gap ∆ vanishes with both u 3 , u 4 > 0 the system undergoes a second order phase transition [see Fig. 2(a)], corresponding to a diverging correlation length ξ = 1/ |∆| → ∞. Numerical evidence for this transition is presented in panel (b) of Fig. 2, which displays the stationary density of active sites obtained for ω = 0 in a chain of 200 sites. For ∆ > 0 and u 3 ≤ −2 √ u 4 ∆, the transition from the active to the inactive phase takes place instead at finite (c) Effective potential Γ(n) (dashed lines) and corresponding "optimal-path" potential W (n) (solid lines), see Eq. (7), across the first-order transition. At the transition point, W (n1 = 0) = W (n2) = 0. (d) Steady-state histogram of the density in the quantum limit χ = 0 (12 spins) obtained via a quantum-jump Monte Carlo (QJMC) method, indicating a first-order transition (ωc ≈ 2) as ω increases. Two stable stationary solutions, one with zero and one with finite density, emerge. The inset displays a section of the histogram taken at ω = 8.
correlation length ξ = 1/ |∆| < ∞. The form of the effective potential Γ(n) suggests a first-order transition line in this regime featuring the coexistence of the zero and finite-density solutions. This case, however, requires additional care due to the specific form of the noise, as detailed further below.

V. FLUCTUATIONS AT THE CONTINUOUS TRANSITION
The competition between quantum and classical dynamics strongly affects the nature of the active-toinactive transition. In the absence of the coherent coupling, u 4 , µ 4 = 0, the action (5) is equivalent to the socalled Reggeon field theory for classical DP [54]. It features -upon rescaling the fields -the characteristic rapidity inversion symmetry, which leaves the system invariant under the transformation n ↔ −ñ and t → −t [3,7,53]. For u 4 > 0, this symmetry is broken by the microscopic action. The implications depend on the dimension d: For d > 2, u 4 is RG irrelevant and can be discarded in the infrared-dominated dynamics close to the second order transition. Consequently, in d > 2, rapidityinversion is restored and the line of continuous transitions displays universal scaling behavior corresponding to classical DP.
At the α point [white dot in Fig. 1(b)], u 3 = 0 and the leading-order coupling becomes u 4 . For d > 2, the second order transition at this point is governed by mean-field scaling behavior, since u 4 is RG-irrelevant and cannot introduce infrared divergent corrections to the vanishing couplings u 3 , ∆. On the other hand, for d < 2, u 4 becomes RG relevant and generates a non-trivial RG flow of ∆ and u 3 on the entire second order transition line. This leads to a violation of rapidity-inversion which persists at long wavelength, and thus drives the system away from the DP critical point to a new non-equilibrium universality class, without specific symmetries. In d < 2, therefore, only the isolated point χ = 1/z, ω = 0 lies in the DP class, while the presence of quantum fluctuations imprints a new universal scaling behavior on the entire line, including the α point. For marginal u 4 in d = 2, the scaling of the fluctuation corrections to u 4 determines whether this coupling becomes relevant, making the scenario equivalent to d < 2, or irrelevant, which has to be determined by an RG analysis of the problem.

VI. NON-EQUILIBRIUM DISCONTINUOUS TRANSITION
For (∆ > 0, u 3 < −2 √ ∆u 4 ) the effective potential Γ displays two distinct minima, n 1 = 0 and n 2 = |u3| 2u4 + ( , suggesting a first-order phase transition. The actual transition line, however, lies where the finite-density minimum becomes statistically preferred. In equilibrium, this would be the point at which the minima of Γ are at the same height. However, the present non-equilibrium noise shows more pronounced fluctuations at larger densities and thus favors the absorbing minimum n 1 with respect to n 2 . To estimate the steady state distribution function P (n), we apply the optimal path approximation to the action [3,53]; this involves treating the coefficient Ξ(n) = 1 2 n + µ 4 n 2 ofñ 2 as a kind of mean-field and density-dependent temperature. It yields (see Appendix D) with volume V and normalization Z. Both potentials W (n) and Γ(n) vanish in n 1 and share the finite-density minimum n 2 . In the thermodynamic limit V → ∞, P (n) → δ(n − n l ), where l = 1, 2 depending on which one is the global minimum of W , accounting for the physical constraint n ≥ 0. The transition takes place when W (n 2 ) = 0, which identifies the non-equilibrium firstorder line [dashed line in Fig. 1(b)]. Due to the nonequilibrium nature of the fluctuations, this does not coincide with the naive prediction Γ(n 2 ) = 0, as shown in Fig. 2(c). In Fig. 2(d) we report the full-counting statistics of the density n obtained via QJMC techniques [55] for a chain of 12 spins. Despite the presence of strong finite-size effects, a bimodal structure is still highlighted for large values of ω. This implies that trajectories bunch together around two possible values, the absorbing one and a finite-density one and is a signature of the aforementioned coexistence.

VII. REALIZATION WITH RYDBERG ATOMS
Instances of this crucial competition between classical and quantum processes can be implemented with cold atoms excited to Rydberg states [18,[56][57][58][59][60]. They are represented with two internal states, the ground state |GS ≡ |i (inactive site) and the excited one |Ryd ≡ |a (active site). Rydberg gases feature strong van-der-Waals interactions in state |a [24][25][26], which rapidly decay as r −6 with the interparticle distance r. For the sake of simplicity, we approximate it here as a nearest-neighbor interaction of strength V nn in a one-dimensional configuration.
Quantum branching/coagulation is realized via coherent driving by a laser field of Rabi frequency Ω and detuning ∆ L with respect to the atomic transition frequency; fixing ∆ L = −V nn enables an "anti-blockade" [58,61,62] mechanism which favors the excitation of a Rydberg atom next to an already excited one, e.g. |iai → |iaa . Differently from the idealized model above, the constraint requires here a single excitation nearby, and processes such as |aia → |aaa are highly suppressed. The Hamiltonian is therefore approximately given by H ryd = Ω k Π k σ x k where Π k = n k−1 + n k+1 − 2n k−1 n k+1 . To generate the incoherent branching/coagulation the atoms are coupled (with coupling g) to a second equallydetuned light field with strong phase noise (dephasing rate λ g) [63]; for a correlation length shorter than the interatomic distance, the bath is modeled as independent bosonic modes b k , b † k acting on each lattice site. The effective equation of motion for the atoms is obtained by performing second order perturbation theory in the small parameter g/λ [44,60,64]. The resulting master equation for the reduced atomic density matrix ρ iṡ For sufficiently high ( b † k b k 1) and homogeneous intensity, one can identify κ = (4g 2 b † k b k )/λ, leading to the branching/coagulation jump operators: The final process is radiative decay of an atom from its Rydberg state to the ground state, modeled by the jump operator L ryd d,k = √ γ σ − k [26]. Although the microscopic formulation of the dynamics is slightly different from the previously-discussed model, the resulting phase structure is similar, as the EOMs only differ from Eqs. (2)(3)(4) by RG irrelevant higher order density terms. In particular, they leave the universal properties near the continuous transition points unchanged.

VIII. OUTLOOK
We have investigated the effects of quantum dynamical processes on a prototypical absorbing-state phase transition. We highlighted the emergence of a richer structure in the phase diagram, which includes both a discontinuous and a continuous non-equilibrium transition. In low dimension d < 2 the presence of a quantum coherent process leads to a breaking of the only fundamental symmetry of DP in a way that persists at long wavelengths, and thus leads to a phase transition of a different nature. In equilibrium, the interplay between classical (thermal) and quantum fluctuations typically leads to a dimensional crossover [2,65]. The present work shows that out of equilibrium the picture is not as straightforward and opens the path for further investigations in this field, including the quantitative characterization of the new universality class.
where D [X]O = X † OX − (X † XO + OX † X)/2 andξ O is the quantum noise term for the operator O. The noise-less equation of motion is only formally correct on the level of single-operator expectation values, while the noise contributes by preserving the (anti-)commutation relations of the operators during the evolution. [49]. For a linear coupling of the system to the bath, the noise is typically Gaussian, with zero mean, but nonvanishing time-and space-local correlations. This prescription leads to Eqs. (2)-(4); below we provide the main conceptual steps. There exist several canonical (and equivalent) strategies to determine the properties of the noise operatorŝ ξ n,x,y , as for instance outlined in Ref. [49]. Here, we follow a path relying on the unitary Heisenberg equations of motion for system operators in the presence of a bath. As a simplifying assumption, we imagine the spatial correlations of this bath to be shorter than the typical interparticle distance in the system. This allows us to describe every spin as coupled to its own bath. We can therefore focus on a single spin as a representative and we model the spontaneous emission dynamics via the simple Hamiltonian where the b l s represent a set of bosonic bath operators. We further assume that this bosonic reservoir is kept at zero temperature and that the number of modes is sufficiently large to allow a continuum description with a density of states D(ω) = q δ(ω − ω q ). Taking the von Neumann equation for the global (spin plus bath) density matrix and eliminating the bath degrees of freedom leads then to the jump operator L d = √ γσ − with γ = 2π[λ(0)] 2 D(0) being proportional to the bath density of states D(0) and the couplings λ(0) evaluated at zero frequency (see e.g. Chapter 8 of [49]). The Heisenberg equations of motion for the operators are thereforė Formally solving Eq. (12) yields Inserting this solution into Eqs. (10), (11) and performing the Born-Markov approximation leads tȯ Definingξ − (t) = (ξ + (t)) † and taking the bath to be in the vacuum state (corresponding to spontaneous emission), we find therefore the noise properties in the Born Markov approximation to be By noticing thatξ n ≡ −ξ + σ − − σ +ξ− one gets all the remaining non-vanishing correlations Rotating into the (x, y, n) basis and introducingξ x = ξ + +ξ − andξ y = −iξ + + iξ − (and analogously σ x and The dependence of the ξ n noise on the density keeps track of the fact that the absorbing configuration n = 0 represents a fluctuationless state in the entire parameter regime, which forbids a density independent contribution to ξ n in Eq. (17). The Markovian noise level introduced by the decay terms L d = γσ − only appears as an additive noise in the σ x,y variables. Including classical coagulation and branching processes yields additional noise terms. However, due to the presence of the term j n j these contributions will always be higher-order in the density and are therefore subleading with respect to the ones derived above in the absorbing phase. Extending the system from a single spin to a lattice of individual spins, an equivalent computation shows To leading order in the density, this yields the same noise terms reported above. Since the coherent branching and coagulation does not produce an additional noise, this concludes the derivation of the Heisenberg-Langevin equations. retarded regularization [3,53]. Introducing three sets of imaginary response fieldsñ X ,σ x X ,σ y X , exploiting the Fourier transform δ(f (n)) = Dñ exp(−ñf (n)) and integrating over the noise variables ξ n/x/y , Z can be cast into a path-integral form where, up to leading order in a (spatial) derivative expansion, the action reads Since zχ+1 2 ≥ 1/2 throughout the physical parameter region χ ≥ 0, both the σ x and the σ y fields remain gapped and one can therefore neglect the subleading derivative and fluctuating terms within square brackets. This yields an action which is separately quadratic in (σ x X ,σ x X ) and (σ y X ,σ y X ). These modes can be actually integrated out exactly and, up to RG-irrelevant terms, one obtains action (5). The couplings correspond to the following combinations of miscroscopic parameters: µ 4 = 2z 2 ω 2 (zχ + 1) 2 + 128z 4 ω 4 (zχ + 1) 6 . (34d) Our procedure differs conceptually from the approach advocated in [66]. We found it necessary to accurately capture the short distance physics of the problem.

C. ADDITIONAL DETAILS ON THE NATURE OF THE OBSERVED PHASE TRANSITIONS
Except for the α point, in the proximity of the transitions u 3 = 0 and one can rescale the fields according to n → Kn,ñ →ñ/K with the choice K = 1/ √ 2u 3 . Thus, one finds with κ 3 = Ku 3 = 1/2K = u 3 /2. In the absence of the u 4 coupling, i.e. for ω = 0, this is the action describing classical directed percolation. It features the characteristic rapidity inversion symmetry under n ↔ −ñ, t → −t.
For ω > 0, the relevance of the u 4 coupling has to be considered, which is determined by the scaling behavior of the fields n,ñ. In the absence of a thermal fluctuation dissipation relation, both fields n,ñ typically have the same scaling dimension [3,53,67]. This leads to an upper critical dimension of d c = 4 for the cubic coupling u 3 and an upper critical dimension of d c = 2 for the quartic coupling u 4 . In dimensions d > 2, u 4 renormalizes to zero in the RG flow and the rapidity inversion symmetry is restored in the infrared regime. Hence, the effective low frequency theory, and therefore the long time dynamics, is again described by the directed percolation class. On the other hand, for d < 2, u 4 is relevant in the renormalization group sense and the absence of rapidity inversion introduces a different non-equilibrium dynamics at the phase transition, which is not captured by the DP universality class. In d = 2, the quartic couplings are marginal and whether they become relevant or irrelevant in the RG flow has to be determined by a renormalization group analysis of the problem.
At the point α, u 3 vanishes microscopically and the rescaling leading to the action (35) is not defined. In dimensions d > 2 this point features a second order phase transition in the absence of the rapidity inversion symmetry. Since the leading order term in the effective potential Γ (Eq. (6)) is RG irrelevant in d > 2, one expects meanfield scaling behavior at this point. On the other hand, in d < 2 the effective theory at the α-point corresponds to the new universality class in the absence of rapidity inversion.
In any experiments with cold atoms, the presence of small fluctuation-inducing terms ∼ ∆ x σ x or ∆ y σ y is hardly avoidable at the microscopic level. These will generate fluctuations on top of the absorbing state and lead to a temperature type term ∼ñ 2 X T in the action (35), with T ≈ ∆ x , ∆ y . The present discussion of the non-equilibrium phase transitions is then valid on length scales l −1 ≥ √ T .

D. DETAILS ON OPTIMAL PATH APPROXIMATION
In the present setting, the noise Ξ X ≡ 1 2 n X + µ 4 n 2 X increases monotonically with the density. As a consequence, it favors the (fluctuationless) zero density solution over the finite-density one. In order to determine the distribution function for the density variable in the vicinity of the active-to-inactive transition, we apply the optimal path approximation [3,53] to the partition function. Note that the system remains gapped for ∆ > 0 (where the first-order transition is expected to take place) and therefore we can -as a first approximation -neglect spatial fluctuations and approximate n X ,ñ X by spatially homogeneous but temporally fluctuating fields n t ,ñ t . This yields the action with the shorthand Γ (n) = δΓ/δn. The optimal path for the configurations n t ,ñ t corresponds to the configurations for which the non-fluctuating part of the action vanishes, i.e. for whichñ t Γ (n t ) −ñ 2 t Ξ(n t ) = 0. This equation shows the two trivial solutionsñ t = 0 and n t arbitrary as well as n t = 0 andñ t arbitrary. Apart from this, there exists the non-trivial solutioñ n op t = Γ (n t ) Ξ(n t ) .
Considering only configurations which correspond to the optimal path the action becomes with the change of variable ∂ t n t dt → dm and where n 0 is the initial condition -whose specific value is irrelevant -and n the steady-state value of the density. The corresponding density distribution function is with Z = n e −V W (n) , if the integral exists. − ∆ u4 and coincide with the ones of Γ(n). The two functionals, however, may differ significantly. In particular, the global minimum of W does not necessarily coincide with the global minimum of Γ, such that the presence of a non-equilibrium noise term can strongly modify the phase boundary as a function of the noise strength. Since W (0) = 0 for all parameters, the first order transition line separating the active from the inactive phase for ∆ > 0 is determined by the equation W (n 2 ) = 0.