Stochastic Differential Equations for Quantum Dynamics of Spin-Boson Networks

The quantum dynamics of open many-body systems poses a challenge for computational approaches. Here we develop a stochastic scheme based on the positive P phase-space representation to study the nonequilibrium dynamics of coupled spin-boson networks that are driven and dissipative. Such problems are at the forefront of experimental research in cavity and solid state realizations of quantum optics, as well as cold atom physics, trapped ions and superconducting circuits. We demonstrate and test our method on a driven, dissipative two-site system, each site involving a spin coupled to a photonic mode, with photons hopping between the sites, where we find good agreement with Monte Carlo Wavefunction simulations. In addition to numerically reproducing features recently observed in an experiment [Phys. Rev. X 4, 031043 (2014)], we also predict a novel steady state quantum dynamical phase transition for an asymmetric configuration of drive and dissipation.


I. INTRODUCTION
Our understanding of any physical system inevitably relies on our ability to subdivide the object of study into a "system" and a "bath". In this sense, open quantum many-body systems are ubiquitous in nature and in practical applications. A wide range of theoretical and computational techniques have been developed in condensed matter physics to study many-body systems that equilibriate due to the coupling to their environment and can be described by equilibrium statistical mechanics in the long-time limit. In recent years attention has shifted to nonequilibrium many-body systems. For example their time evolution after a quench [2][3][4][5], where the question of whether and by what mechanism they thermalize [6][7][8][9][10] has become a focus of study. Also of great interest are driven systems; applications range from the dynamics of ultra-cold atoms [11][12][13][14][15][16][17], trapped ions [18], coupled light-matter systems [1,[19][20][21][22][23], transport problems [24][25][26][27], and simulated quantum annealing [28]. The need to account for quantum coherence that may be long-ranged in the presence of external forcing and dissipation is an intriguing problem that calls for the re-evaluation and extension of established techniques developed originally for near-equilibrium correlated systems. In recent years we have seen the further development of a number of such powerful computational techniques: exact diagonalization and density matrix renormalization group methods [29,30], nonequilibrium versions of dynamical mean field theory [31], application of Bethe ansatz techniques to nonequilibrium quantum transport in impurity models [32], and various Quantum Monte Carlo algorithms, among them the continuoustime Monte Carlo algorithm (again for quantum impurity models) [33] [34]. Nonequilibrium quantum dynamics on the other hand has been a sine qua non of quantum optics and laser physics. In particular, an arsenal of pow-erful techniques have been developed in the field of Cavity QED to deal with nonequilibrium dynamics of open quantum systems [35][36][37][38][39][40]. Master equation methods, methods based on Heisenberg-Langevin equations of motion, and Monte Carlo Wavefunction (MCW) approaches are by construction ideally suited to study dynamics of open quantum systems. With the experimental progress in Cavity QED in atomic, semiconductor and superconducting circuit systems the attention has recently been drawn towards exploration of quantum many-body phenomena in extended light-matter systems. In particular lattices or networks of cavity QED systems [19,20,[41][42][43][44][45][46][47] present a challenge to established techniques due to the exponential proliferation of the Hilbert space with the system size. For one dimensional systems DMRGbased approaches [48][49][50][51][52] have been extended to study the dynamics of the density matrix of open quantum systems, but these rely on reduced dimensionality and certain constraints in the generation of entanglement during the evolution of the open system. There is a clear need for advancing computational approaches that are more immune to the exponential growth problem and which scale more favorably with system size. Phase-space representations of quantum mechanics possibly offer such an approach, which we explore in this paper.
Our goal in this paper is to develop a novel stochastic approach to study the dynamics of driven and dissipative systems involving spins and bosons, such as cavity QED systems. We generalize the positive P −representation of quantum mechanics to model the dynamics of an interconnected network of spins and bosons coupled linearly to bosonic quantum baths. Phase-space representations have been employed in the past to study purely bosonic open systems [53,54], but there is little work on spin systems and none that we know of for open spin-boson networks. Barry and Drummond [55] used the positive P −representation for spins to simulate equilibrium ther-modynamic properties of the quantum Ising model. Ng and Sørensen [56] used the mapping to Schwinger Bosons to derive a positive P −representation for a spin 1/2 system. Closest to our approach are Refs. [57,58] which employ spin coherent states to derive a Fokker-Planck equation for the Q−function of the single-site Dicke model. In fact the present work is originally inspired by this work, to go beyond the Fokker-Planck level (which is numerically infeasible to solve) and develop a stochastic description, turning a formal identity into a numerical method. To this end a different representation is needed, as the Q−function does not possess a positive semi-definite diffusion matrix.
As we will present in some detail, we use a combination of bosonic and spin coherent states to map a quantum master equation to a Fokker-Planck equation, and in a second step, onto a stochastic differential equation which can be simulated efficiently. Notably, the latter step is only possible if the corresponding diffusion matrix in the Fokker-Planck equation is positive semi-definite. This is guaranteed in the positive P −representation [54,59,60]. To evaluate the effectiveness and accuracy of our computational approach, we analyze in detail a two-site system -a dimer -each site of which features a photonic mode coupled to a local spin (this system has been recently studied in a circuit QED setup [1,61]). We make a numerical comparison of our approach to the Monte Carlo Wavefunction technique for spin values accessible to the latter. We refer to this system as the Dicke dimer when the spins are taken to be large, a limit which lies beyond the capability of the Monte Carlo Wavefunction approach, but as we demonstrate is accessible in this new approach. We stress that it can be applied to more complicated network geometries.
Our paper is organized as follows: In section II, we specify the quantum model that we use as an example to demonstrate our formalism. The first step in our mapping, the derivation of a Fokker-Planck equation from a quantum master equation proceeds via the introduction of bosonic and spin coherent states, and is presented in III. In a second step, we map the Fokker-Planck equation to a stochastic differential equation in section IV. The method is tested on a physical model in section. V. Finally, we summarize and discuss our results in section VI.

II. SPIN-BOSON NETWORKS
A broad class of models in quantum optics and quantum information theory falls into the class of the following network-type, involving spins and bosons: Parts of the Hamiltonian consists of a sum of local terms, describing e.g. external magnetic fields, on-site ener-gies, coherent drives and interactions between bosons and spins. The kinetic Hamiltonian couples different sites of the network. A more specific class of models neglects the coupling of different spins, J spin ij = 0 , but still allows bosons to hop on the network. We also consider a specific type of on-site interactions: . For spin 1/2 the single site system is known as the Jaynes-Cummings model, and it plays a prominent role in cavity and circuit QED systems, but appears also in other contexts. We consider its generalization to a network, where the spin sizes are arbitrary. We will refer to this model as the Dicke network. It is characterized by a matrix of photon hopping amplitudes J ij , a cavity frequency ω c , a spin frequency ω s , a matter-light coupling of strength g, and coherent drive amplitudes f i at each site. Furthermore, the quantum number s specifies the spin representation. Note that our model does not contain any counter-rotating terms. Hence, the isolated Dicke network (f i = 0) conserves the total excitation number iN i +Ŝ z i , whereN i denotes the photon number on site i, andŜ z i are the corresponding z components of the spin. Models (1) and (2) are tractable with our approach that we present in this paper. We focus on the dynamics of model (2), including possibly a coupling of the system to a bath. We therefore use the Lindblad master equation Dissipation in this master equation is described by the Lindblad superoperators for the spins L S and photons L a , describing weak coupling to a bath (we omit the site indices i to keep the notation simple): Here, the constants κ and γ specify the decay rate of photons from each cavity, and the spontaneous decay rate of each spin.n is the number of photons in the thermal bath and is a measure of temperature (withn = 0 for zero temperature) [37,39,40]. Simulating this master equation numerically becomes intractable for large photon numbers, spins, and/or large networks. As L[ρ] is a linear operator on density matrices, the computational complexity grows quadratically in the Hilbert space dimension, which itself would grow exponentially with network size. We thus aim for a different representation of the problem.

III. FOKKER-PLANCK EQUATION
A. Coherent states and the positive P -representation In this section we extend the positive Prepresentation [37,39,40] to situations involving spin. The positive P -representation makes use of the basis of coherent states, which for bosons are eigenstates of the annihilation operatorâ with eigenvalue α, and |vac denoting the vacuum state for the bosons. The spin coherent states [62][63][64] for the spin s representation of the su(2) spin algebra will be defined as whereŜ + is the raising operator the algebra. Both state labels α and z are complex valued. Note that in our convention, the spin state |z = 0 corresponds to the lowest weight ("spin down") state. We construct the following operators: We similarly define for the spinŝ The normalization ofΛ S comes from the overlap of two spin coherent states [62][63][64]. We now introduce a container variable for the complex numbers that specify sets of operators,Λ a andΛ S , for the network: where n is the size of the network. Combining spin and bosonic degrees of freedom, we then define the following operator which acts on the full many-body Hilbert space of the network: The density operator of the system can now be expanded in our generalized positive P -representation as follows: In the above, we defined the integration measure dα = i d 2 α i d 2 β i d 2 z i d 2 w i . As can be easily verified, normalordered bosonic operator expectation values in the positive P -representation are calculated according to The second line gives an approximation to the expectation value for the case where only N s samples α l (t), β l (t), from the positive P -function are available (in terms of solutions to an equivalent stochastic differential equation to be derived below). An example of an expectation value involving spin is given in appendix D.

B. Mapping to a Fokker-Planck equation
Having introduced coherent states and P -functions, it is worth outlining the general strategy for deriving a Fokker-Planck equation. Consider a general master equation of the formρ where L is an arbitrary Liouvillian or Lindblad operator. As we will show in the next paragraph, the previously introduced operatorsΛ allow us to convert secondquantized operators into differential ones. As a consequence, where D L (α) is a representation of L in terms of derivatives with respect to α. It turns out that this differential operator consists only of first and second order derivatives and can therefore be written as (Einstein summation convention over site indices is implied) The dependence of the vector A and the matrix D on α will be determined shortly. Using this relation (without further specifying D L (α) at this point), we can derive a Fokker-Planck equation from the Lindblad master equation as follows. We first note that The derivatives can be transferred fromΛ to P by partial integration. We use eq. (17) to find from which we conclude that This is the Fokker-Planck equation. The goal for the remainder of this section is to calculate A i (α) and D ij (α). Note that the better known P -representation is simply obtained from the positive P -representation by substituting β → α * , which enforces the variables α and β to be complex conjugates. In the positive P -representation, the presence of quantum noise violates this conjugacy relation; for more details we refer the reader to [37]. Eq. (16) provides a prescription for finding the differential operator D L , by evaluating the action of the Lindblad operator appearing in eq. (3) onΛ: In order to proceed, we need the action of the second-quantized operators appearing in the Hamiltonian onΛ, (see Ref. [37] for a pedagogical introduction to the positive P representation). We list and derive those identities in appendix A. Note that each creation or annihilation operator corresponds to a first order differential operator. Commutation relations reflect themselves in the non-commutativity of these differential operators. Furthermore, the master equation contains only products up to second order in bosonic creation and annihilation operators and spin raising and lowering operators. This guarantees that the resulting partial differential equation is necessarily of second order in α and first order in time, and hence of Fokker-Planck type [65].
We first simplify the contributions arising from the local HamiltoniansĤ i (derivations are presented in the appendices). Note that all complex fields α, β, z, w carry a site index. In order not to overload the notation at this point, we omit these indices and instead indicate the site index on the brackets: We next compute the kinetic term for the network, being the only term coupling the different sites: Finally, we calculate the dissipators, starting with the photons. After a straightforward calculation, making use of (A1), we find Notably, all second-order derivatives are proportional tō n. According to the Feynman-Kac relation, second order derivatives correspond to noise terms in a stochastic description. Asn(T = 0) = 0, there is no quantum noise associated to cavity loss at zero temperature. A much longer calculation, shown in the appendix, results in the following contributions from the spin dissi-pators: +γ n + (n + 1) The double arrow indicates an identical contribution on the last line where the roles of z and w are interchanged. Interestingly, the spin dissipator does contain quantum noise terms that are present even at zero temperature. Our goal will be to cleanly separate quantum from classical dynamics. To this end, we perform a transformation on the following variables: However, in order to keep the notation clean, we omit the tildes in the remainder of the paper. The third transformation is needed as photon densities scale as n ∼ αβ. Intuitively , as the field amplitude gets scaled, the photon density of the bath and the external drive have to be scaled up as well (otherwise the bath temperature would be effectively lowered). After the transformation, the Hamiltonian contribution to the Fokker-Planck equation is independent of s, and hence it has a well-defined limit for s → ∞. In contrast, the second-order differential operators that arise due to the interaction g are proportional to s −1 and therefore vanish in this limit. This means that quantum noise vanishes in the classical limit of large spin, as intuition would suggest [57].

C. Drift vector and diffusion matrix
We are now ready to collect all terms and specify the drift vector A and the diffusion matrix D in the Fokker-Planck equation (19).
The vector A has 4 complex entries for each network site. The 4 entries corresponding to site i are given by As we shall explain below in more detail, these are the right hand sides of the classical equations of motion for the variables α i , β i , z i and w i , respectively. We now move on to the diffusion matrix D. Note that it is block-diagonal in the space of network sites; therefore we will focus on a single site, omitting the site indices. For later convenience, we will separate it in the following way The non-zero entries are themselves 2 × 2 matrices These matrices are proportional to the parameters κ, g and γ, respectively, indicating three distinct sources of noise, namely a quantum noise contribution due to g, a purely thermal noise arising from κn, and a noise contribution from the spin, associated with the spontaneous emission at rate γ. In contrast to the other noise contributions, D (2) couples the photon and spin sectors. Also note that it is proportional to 1/s, and hence vanishes in the classical limit.

IV. STOCHASTIC DIFFERENTIAL EQUATIONS
In order for our approach to provide an efficient basis for numerical simulation, we furthermore map the Fokker-Planck equation we have derived onto a set of stochastic differential equations. A necessary requirement for this step to be possible is that the diffusion matrix be positive semi-definite. If this requirement is not met, then the diffusion matrix possesses contracting directions, and it is impossible to model contracting densities in terms of random walks. This requirement is precisely why we have chosen to work in the positive P -representation, as then the positivity of the diffusion matrix is assured [37,39,40].
It follows from the standard theory of stochastic calculus [39,40] that the Ito stochastic differential equations equivalent to the Fokker-Planck equations (19) are given by Hence, the deterministic evolution is described by the drift vector A, while the equal-time correlator of the noise in the four dimensional complex space is given by the diffusion matrix. We are now faced with the challenge of designing a noise such it satisfies this requirement. For the following discussion, we focus on a single network site and omit the site index (the diffusion matrix is blockdiagonal). An obvious way of creating such a noise term would be to define where 1 4 is the four dimensional identity matrix. The noise is thus decomposed into the product of B(α, t) (the matrix square root of D), with dW , a four dimensional vector of independent Wiener increments. However, as D is a complex 4 × 4 (or real 8 × 8) matrix for each network site, determining this matrix decomposition numerically at each infinitesimal time step would be computationally demanding. Fortunately, we can use a trick to circumvent the need to perform such a time-dependent factorization, making the numerical algorithms more efficient, by using the explicit decomposition of D given in eq. (27). In the following, we will construct three infinitesimal noise vectors dξ 1 (t), dξ 2 (t), and dξ 3 (t), which are taken to be mutually uncorrelated, and which satisfy for each i = 1, 2, 3 (no summation over i) with the individual D (i) given in eq. (27). We then define the total noise dξ(t) = dξ 1 (t) + dξ 2 (t) + dξ 3 (t) .
As the dξ i are mutually uncorrelated, it follows that = D µν δ tt dt , after using (27) and (32). Consequently we have succeeded in generating a random noise vector with the desired correlator (30). It remains still to show how to construct the individual dξ i 's. This is, however, easy. For this purpose, we construct the matrix square roots B (i) for the matrices D (i) such that D (i) = B (i) B (i),T . Due the structural simplicity of the matrices D (i) we find the analytic solutions A similar analytic formula for B (3) can be obtained from diagonalizing a 2 × 2 matrix, but in the following we will set the spontaneous emission rate γ to zero, and so have no need for it. We now define with each dW i (t) being a four dimensional vector of independent Wiener increments. It follows that the individual noises dξ i have the correlators (32). This concludes our construction of the correlated noise. We note here that our noise vector can be multiplied from the right by any complex orthogonal matrix, and would still satisfy eq. (30). This observation is a realization of the stochastic gauge degree of freedom [56,66].

V. NUMERICAL SIMULATIONS Nonequilibrium Dicke-dimer
We shall now apply our stochastic formalism to a physical test case, one involving strong spin-photon interactions, and which has an additional spatial degree of freedom involving a kinetic energy term. For simplicity we will focus on the dynamics of two coupled cavities (a dimer), each of which contains a spin coupled to a single photonic mode. In circuit Quantum Electrodynamics (circuit QED), each site would be realized in terms of a microwave field, coupled to a superconducting qubit, with the sites capacitively coupled to allow photon hopping. As the qubit can be interpreted as a spin s = 1/2, each cavity is described by a Jaynes-Cummings model. A dimer of such cavities has recently been studied in an experiment [1]. We consider a broader class of similar systems, allowing for arbitrary spin s. In particular, we are interested in the scaling limit of large spin and photon numbers and in the quantum to classical crossover. A corresponding dimer could e.g. be realized by using many qubits per cavity that are all coupled to the same photonic mode.
To begin with, let us briefly review the main experimental findings. In the experiment, one of the two cavities (cavity 1, which we will take to be the left cavity) is initially populated with many photons. The system is undriven and is let to evolve in time. As both cavities are dissipative (with photon loss rate κ), the photon number decreases monotonically with time. The following physics is observed in the experiment. At large photon numbers, photons are observed to undergo linear periodic oscillations between the two cavities, while simultaneously exponentially decaying to the outside environment. However, as the photon number drops below a certain critical threshold, the oscillations are seen to cease as the system enters a macroscopic quantum selftrapped state (for details please refer to [1]).
A qualitative explanation of the physics is as follows: There are two competing time scales for the dimer when observing the homodyne signal (equivalently, there is a competition between the on-site interaction energy, with scale set by the cavity-spin coupling g, and the kinetic energy, dictated by the hopping rate J). The Josephson oscillations occur with period t J = 1/2J when the photon number is above the critical threshold and the system is in the delocalized phase. The second time scale is the collapse and revival period associated with the single site Jaynes-Cummings physics, the relevant time scale when the system has localized and the tunneling disappears, wherein the two sites are effectively decoupled. The localization transition is predicted to occur when these two time scales become comparable. This matching argument has been supported by extensive numerical simulations for the spin 1/2 dimer using Monte Carlo Wavefunction simulations [1].
We would like to analyze the quantum transition in the well-controlled semiclassical limit of large spin, going beyond the classical solution and taking into account the impact of quantum fluctuations, as well as the effect of thermal noise. We will also give a theoretical explanation for the super-exponential decay of the homodyne signal that has been observed in [1].
We test our method on the example of a dissipative Dicke dimer. We are interested in two cases. First, we study the undriven lossy dimer, inspired by the experiment, where we prepare an initial state with a fixed number of photons in the left cavity and study the dynamics. Second, we study the corresponding driven system. Here, we are mainly interested in the behavior of the tunneling current between the two sites, in steady state. As we show, the driven system displays a dynamic quantum phase transition visible in the inter-cavity current upon varying the interaction strength.

A. Undriven dissipative Dicke dimer
Zero temperature, infinite spin. We begin by modeling the classical dynamics of the dissipative Dicke dimer at T = 0. Zero temperature (n = 0) in combination with infinite spin implies that all noise terms vanish. The equations (29) are then completely deterministic, and the positive P -representation becomes equivalent to the standard P -representation, allowing for an alternative representation of the stochastic differential equations in terms of the compact angular variables. We have found the corresponding equations (F2) to be more stable at long times in this new representation. Simulation results for a decaying dimer are presented in figure 1. We observe a self-trapping transition setting in at a critical photon number, below which the oscillations die out rapidly. The dynamic equations in this limit are equivalent to the classical Maxwell-Bloch equations that have been studied earlier in this setup [1,61], but whose derivation requires the uncontrolled assumption of the factorization of operator expectation values, whereas our deviation is fully controlled in the scaling limit.
Finite temperature, infinite spin. Next, we explore the impact of thermal noise. To this end, we setn to a finite, positive value, simulating the coupling to an external photon bath with mean occupationn.
The fact that we have taken the spin s → ∞ and the qubit relaxation rate γ → 0 implies that the only remaining noise term in eqs. (35,36) is ξ 1 (t), corresponding to the matrix B (1) . This matrix has an interesting symmetry property: when multiplying on the right any real vector, the resulting complex vector has just two entries which are always complex conjugates. As a consequence, the random thermal noise acting on the variables FIG. 2. Finite temperature, classical simulations for the selftrapping transition (infinite spin). We averaged over 10,000 stochastic trajectories. The plot shows particle number (red) and homodyne signal (blue). At the transition, the homodyne signal decays super-exponentially.
α and β preserves this conjugacy. The drift equations share the same property, and preserve conjugacy, namely A 1 (α) = A * 2 (β) for β = α * , where * denotes complex conjugation. It follows that α(t) = β * (t) for all times and all stochastic trajectories. This mirrors the fact that the positive P -representation is equivalent to the ordinary P -representation in the absence of interaction induced (quantum) noise.
Physically, coupling a thermal bath to our decaying cavity will induce two things: first, coherence will be destroyed over time, and second, the system's equilibrium state (at least for small g) will not be the vacuum state but rather an incoherent photon state at mean photon numbern. To illustrate the loss of coherence, we calculated the homodyne signal h = Î 2 + Q 2 , where the quadratures are defined in terms of the creation and annihilation operators asÎ = (1/2)(â +â † ) andQ = (i/2)(â † −â). This quantity was experimentally measured in Ref. [1] in the closely related setup of a decaying Jaynes-Cummings dimer, i.e. for s = 1/2. Note that for a perfectly coherent system where â †â = â † â , the homodyne signal measures the photon number. In the presence of some incoherence, however, it drops below the photon number. In Ref. [1], the homodyne signal was seen to decay super-exponentially close to the self-trapping transition. Our simulations show a qualitatively similar behavior in figure 2, where the homodyne signal (but not the photon number) is seen to decay super-exponentially. In the experiment, individual photons escape according to a Poisson process. For each single trajectory in the ensemble average, the photon number drops below the critical threshold at a random time, the initial photon number determining the average time at which this occurs. On approaching the transition, the The upper panel shows the photon numbers in the left (red crosses) and right (blue dots) cavities, respectively. The high frequency oscillations are Rabi oscillations, whose frequency depend on the photon number. They are also apparent in the lower panel, showing the z-components of the two spins (same color coding). We also present the system's total excitation number N 1 +N2+Ŝ z 1 +Ŝ z 2 (black dashed line in upper panel), which is a constant of motion for the closed system, but here slowly and smoothly decays due to the cavity loss. Jumps in this line indicate the breakdown of numerical reliability, here at t = 4.3.
oscillations become highly nonlinear, with a diverging period (critical slowing down) [61]. This results in a dephasing of the different trials within an ensemble. Therefore averages of the homodyne signal die out faster than exponentially. Hence, the quantum localization transition in [1] possesses a classical analogue at large spin.
Finite spin In contrast to the semiclassical limit of infinite spin, finite spin simulations require the full machinery of the positive P -representation. In particular, the emergent quantum noise at finite spin violates the conjugacy relation between α and β (and also z and w). Hence, all four complex coordinates evolve according to their individual dynamics, and all are subject to individual (yet correlated) sources of noise. This fact has counterintuitive consequences. Let us consider for example of the photon density in a given cavity. An individual stochastic trajectory in the P -representation necessarily has positive photon numbers n ∼ α(t)α * (t) ∈ R + . In positive P -representation, individual runs have generally complex contributions n ∼ α(t)β(t) ∈ C. It is therefore important to keep in mind that only averaged quantities have a physical meaning. Similarly, the z−component of the spins in the P -representation are given byŜ z ∼ 1−zz * 1+zz * ∈ [−1, 1]. In contrast, the corresponding contribution in positive P -representation readŝ S z ∼ 1−zw 1+zw ∈ C. This behavior is seen for example in the lower panel of figure (3), which shows the real part of this expression for a single stochastic run. Note that the averaged rescaled z-components of the spins are always between −1 and 1.
Studying the lossy cavity, we are faced with typical problems that arise in positive P -representation simulations: individual stochastic trajectories show "spikes", as is apparent in figure (3). Such spikes are a well-known problem in the context of positive P -representation simulations [37,56,66]. They indicate that the underlying P -function is heavy-tailed. When the tails become so heavy that the second moment diverges, stochastic averaging fails to converge beyond the time where spikes proliferate. We analyze the spike statistics in the next section. In extreme cases, the positive P -function can even have non-vanishing mass at infinity, spoiling our derivation of the Fokker-Planck equation, which relied on a partial integration and the dropping of surface terms. In some cases, however, single trajectories can already predict much of the physics. In figure 3, we show such a characteristic stochastic trajectory. We see the onset of self-trapping before the simulations break down. The dashed black curve in figure 3 shows the sum of both z-components of the two spins plus the total number of photons in the two cavities. For the isolated system, this is a conserved quantity, while for the open system this quantity is expected to smoothly decay. This is indeed what can be extracted from the plot, and as long as the dashed black curve is continuous and smooth, the numerical simulations can be trusted. It is interesting to see from this plot that for some time before the simulations break down, single stochastic trajectories undergo Rabi oscillations with the spin amplitude exceeding the classical allowed bound. This is a clear signature of the simulations entering the quantum regime. It was not possible for these parameters to carry out an ensemble average for long times without truncating divergent trajectories. We next consider a driven version of this system, which we can place into a steady state, ameliorating such problems and demonstrating the power of the positive P -simulations.

B. Driven dissipative Dicke dimer
Having studied a strongly interacting quantum system weakly coupled to the environment, we now study the case of a strongly driven, dissipative system. As we shall see, strong drive and dissipation will help to stabilize the positive P -representation simulations. In the steady state of a driven system, autocorrelations quickly decay in time, and the spikes are damped before having the opportunity to grow. As the system will relax into a steady state, we are able to simulate long times and even small spins.
We consider two coupled cavities each supporting an atom with spin s. However, instead of filling the system initially with photons and letting them decay over time, here we start with the cavities in the vacuum state and coherently drive the left cavity, such that a steady state emerges. We choose a hopping rate J = 1 (the other parameters are measured relative to J), κ = 20 for both cavities, and set a coherent drive with amplitude f = 100/ √ 2. In the absence of the second cavity, and for g = 0, this would lead to a steady state photon number of 50. We vary the interaction strength g from 0 to 10. Our observable is the photon current.
Non-interacting limit For g = 0, the stochastic equations become deterministic and reduce to (α's are the expectation values of the annihilation operators in a coherent state)α While an analytic time-dependent solution exists, even more straightforwardly the steady state values can be derived by setting the time derivatives to zero, yielding We will use this result as a reference. In the following we will consider the regime of finite g and s. Classical simulations, finite coupling g In the limit of infinite spin, we simulated the deterministic equations numerically. Figure 4 shows the time-dependent currents for different values of g. Below a critical value of g c ≈ 7J, the current is seen to oscillate around a positive mean value. Note that the current never changes sign. Above g c , the current vanishes. This delocalizationlocalization transition is what we want to simulate for finite spin, using the positive P -representation.
FIG. 5. Steady state currents j as a function of the interaction strength g. Each data point results from an ensemble average over 6, 000 stochastic trajectories and a subsequent time average. While a sharp transition in the current is seen at infinite spin, finite values of s turn this transition into a crossover. Deviations between MCW and PP simulations grow close to the transition for s = 1. As the sampling error shrinks with larger spin size, we would expect better agreement for larger spins, which we cannot test due to the Hilbert space dimensionality constraints of MCW.
Quantum simulations To study the behavior of the asymmetrically driven cavity in the quantum regime, we use the positive P -representation, scanning through all orders of magnitude of the spin s in a range from 1 to 10, 000. We find that in the quantum case (finite s), the currents saturate to steady state values that strongly depend on g. A strict phase transition only exists for s = ∞, but for large spins, the current is strongly suppressed above g c . This behavior is summarized in figure 5. Here, the time averaged current is plotted as a function of g for various spin sizes. Note that close to the transition, the statistical error grows as the system becomes unstable due to the emergence of spikes. Even for the case of s = 1, a strong nonlinear dependence of the intercavity current on the interaction strength g is seen. This effect should be measurable in a circuit QED experiment.
We also compared our method against a numerical simulation based on the Monte Carlo Wavefunction algorithm (MCW) [1,38]. This is an alternative approach based on an unraveling of the master equation, which allows one to simulate reasonable sized systems (the problem of an exponentially growing Hilbert space dimension still exists in this approach). Figure 6 shows the outcome of a comparison of both methods. In this figure, we plot the dynamics of the photon current as a function of time, starting with a "spin down" state and an empty dimer. The common parameters chosen are J = 1, κ = 20, f 1 = 100/ √ 2, f 2 = 0, and we varied g and s. We find good agreement in the time-dependent particle current for values of g that are below the classical critical value of g c ≈ 7J. For larger values of g, small discrepancies appear. A possible explanation for this is the fact that deep in the quantum regime, individual stochastic trajectories may have negative particle numbers and negative currents. This behavior is also shown in the histogram figure 7, which, at a given time t, counts the number of cases where the particle current is found at a given value. Upon normalization, this can be interpreted as a probability distribution for the current. So long as most of the mass sits in the positive range, positive P simulations and MCW simulations agree reasonably well (here for g = J). If much of the probability mass is in the forbidden region of negative currents, deviations become stronger and the positive P simulations lose their validity.
A further comparison was carried out for the spin dynamics of the right and left cavity, as shown in figure 8, where also good agreement between positive P stochastic simulations and MCW is obtained. While the undriven cavity saturates at a negative value for the z component FIG. 7. Probability histogram of finding the current j at a random time t to be X for a stochastic trajectory in steady state (log-scale). The blue histogram in the foreground shows g = 1, where the stochastic fluctuations are much smaller than for g = 7 (background, red). Note that negative currents (flowing from the undriven, lossy cavity to the driven cavity) are unphysical, as are negative photon densities. When the statistical weight of such contributions is too large, the simulations lose their predictive power. The inset contains the same quantities presented on a non-logarithmic scale, showing that the majority of trajectories have positive currents for g = 1, but less so for g = 7.
of the spin, the driven cavity has an S z component that averages to zero. This can be understood as individual stochastic trajectories undergoing Rabi oscillations with different relative phases, which averages out the z component of the spin in the driven cavity.
This concludes our first application of the generalized positive P -representation as a numerical tool for studying spin-boson systems.

VI. SUMMARY AND CONCLUSIONS
We derived stochastic differential equations to model the nonequilibrium dynamics of systems involving bosons and quantum spins. Our approach is based on a generalization of the positive P -representation using spin coherent states. This allows us to map a large class of Lindblad master equations onto Fokker-Planck equations, following in a second step to a set of stochastic differential equations. Our approach can be applied to a variety of systems, including large networks.
Regarding computational efficiency, our approach scales linearly (instead of exponentially) with the number of network sites for nearest-neighbor couplings, and quadratically otherwise. We also note that, in particular for problems involving coherent photons, we arrive at a much lower dimensional representation than in the usual FIG. 8. Spin dynamics of the driven dimer, comparing positive P simulations (PP, averaged over 10000 trajectories) against the Monte Carlo Wavefunction approach (MCW, averaged over 100 trajectories). The z components of the spins in the left (driven) cavity and in the right (undriven) cavity are shown as a function of time (g = 7J, s = 1). Note that the photon number in the undriven cavity is small, as the cavity loss rate κ exceeds the incoming photon flow. Hence, the corresponding spin excitation is small. In contrast, Sz in the left cavity averages to zero due to rapid Rabi oscillations with the photon mode.

Fock state representation.
We also modeled a dimer, each component consisting of a cavity coupled to a spin (of various sizes), as a simple example. Individual stochastic trajectories were found to display heavy-tailed fluctuations, the socalled spikes. Drive and dissipation reduce these fluctuations and bound the sampling variances. We compared our approach against the Monte Carlo Wavefunction method [38], and found good agreement. For the undriven, dissipative dimer, we were able to qualitatively reproduce the super-exponential decay of the homodyne signal that has been observed in a recent circuit QED experiment [1]. We also studied the corresponding driven system where we predicted a new phase transition in the inter-cavity current as a function of the on-site interaction strength.
We plan to study larger systems than the dimer, i.e. large networks of cavities and spins. For these systems, our method will compete very well with other existing methods due to the favorable scaling properties of our approach.
[â †Ŝ − ,Λ] = (∂ α + β)(−z 2 ∂ z + 2s The operators associated with the cavity frequency map according to The spin frequency term yields Operators associated with a coherent drive result in The kinetic energy term results in This concludes the Hamiltonian contributions to the Fokker-Planck equation. a. Photon dissipators First, we will calculate the dissipators of the photon fields, which are given by Similarly, we find for the ingoing term, Hence, Interestingly, note that there is no noise term forn = 0. In the positive P -representation at zero temperature, all noise comes from quantum fluctuations, and its strength depends on g as opposed to κ. b. Spin dissipators In contrast to the dissipators for the photon field, calculating the spin dissipators, Eq. (??), is much more work. Again, we distinguish between "in" dissipators (existing only at finite temperature), and "out" dissipators. Let us calculate them term by term, using Eq. (??): (n + 1) 2(−z 2 ∂ z + 2s z 1 + wz )(−w 2 ∂ w + 2s w 1 + wz ) −(−z 2 ∂ z + 2s z 1 + wz )(∂ z + 2s w 1 + wz ) − (−w 2 ∂ w + 2s w 1 + wz )(∂ w + 2s z 1 + wz ) Λ .
Now, let's consider the "in" term, Again, collecting second and first order differential operators, and doing a similar calculation as above for the latter results in L S (1) in = γ 2n [2(s + 1)z∂ z + 2(s + 1)w∂ w ]Λ.
stochastic differential equations with a better numerical stability. A closely related transformation has been carried out in [57,58], but in contrast to the latter, we keep the variable α. We consider the inverse stereographic projection, mapping the spin field back on the sphere, where φ ∈ R and c ∈ [−1, 1]. We only transform the spin part and leave the equations for the photon field α unchanged. This transformation results in a new stochastic differential equation of the form dα = A 1 (α, φ, c)dt + √n κ (dW 1 + idW 2 )/ √ 2 (F2) dφ = (ω z + g c √ 1 − c 2 (α * e iφ + αe −iφ ))dt The function A 1 is given by Here, we focussed on a single cavity, andᾱ is the field in the other cavity. We used this set of equations when simulating the finite temperature dynamics of the system in the scaling limit of infinite spin.