Dynamical symmetries and crossovers in a three-spin system with collective dissipation

We consider the non-equilibrium dynamics of a simple system consisting of interacting spin-$1/2$ particles subjected to a collective damping. The model is close to situations that can be engineered in hybrid electro/opto-mechanical settings. Making use of large-deviation theory, we find a Gallavotti-Cohen symmetry in the dynamics of the system as well as evidence for the coexistence of two dynamical phases with different activity levels. We show that additional damping processes smoothen out this behavior. Our analytical results are backed up by Monte Carlo simulations that reveal the nature of the trajectories contributing to the different dynamical phases.

Understanding and controlling the dynamical behavior of quantum systems has seen flourishing interest [1][2][3], propelled by theoretical and experimental progress that has made it possible to observe and manipulate such systems with unprecedented accuracy. Much attention has also been devoted recently to the notion of dynamical phase transitions in such systems, relating them to the non-analyticity of, e.g., the Loschmidt echo [4] or the logarithm of a biased partition function in large-deviation (LD) theory [5], which has a natural interpretation in terms of the statistics of rare trajectories observed in experiment. The study of the dynamics of quantum systems through LD methods [6,7] emerged recently both as an extension of the theory as applied to classical systems [8][9][10][11] and as a dynamical complement to the standard analysis of equilibrium phase transitions in many-body systems [12]. Here, non-analyticities in the LD free-energy function of a system, extracted from the equations governing its dynamical behavior, are identified in the literature with dynamical phase transition points [6].
Following Ref. [6], in this paper we are interested in studying the statistical properties of rare quantum-jump trajectories [13] of a system that interacts with a heat bath. We consider the dynamical LD properties of a simple three-spin quantum open model which departs from those recently studied in two respects: first, dissipation is due to non-classical bilinear jump operators; and second, we consider a currentlike dynamical order parameter. The two central results in the paper are: the observation of intermittency between dynamical phases of distinct activity, itself a consequence of the reducibility of the dynamics in an appropriate limit due to the collective jump operators; and the existence of a Gallavotti-Cohen symmetry in LD functions associated to the timeasymmetric order parameter, analogous to that found in driven classical systems [15], which gives rise to a fluctuation theorem [16] relating to the quantum jump rate. Our system therefore provides a minimal but extendible model that uncovers the effects of thermal baths and the nontrivial interplay between local and global decay channels [17] on the nonequilibrium dynamics of a quantum system.
We start with three spins-1/2, which we label j = 1, 2, 3, placed at the vertices of an equilateral triangle (see Fig. 1), in-teracting via an Ising-type interaction, and in a uniform magnetic field. Suppose that spin 1 is in a harmonic trap, allowing it to move on an axis parallel to the line joining the two other spins to the two other spins, whereas spins 2 and 3 are held tightly pinned. Since the spin-spin interaction depends on the distance between each pair of spins, the motion of spin 1 will modulate the interactions between them. Adding damping of the harmonic motion by introducing coupling to an environment, we obtain a coupling of a collective spin degree of freedom to this environment. We explore the trajectories of this collective degree of freedom, finding evidence of coexisting dynamical phases and nontrivial dynamical symmetries. It is worth pointing out that individual spins can be coupled to the motion of a harmonic oscillator through the use of trapped ions [18], by embedding solid-state qubits into mechanically-compliant structures [19], in nanomechanical resonator arrays [20], on graphene layers [21], or on diamond surfaces [22]. Furthermore, this system lends itself well to being extended by adding more spins, thereby changing the symmetries of the model.
We will first derive the master equation for the spin system by adiabatically eliminating the motion of spin 1. The resulting reduced dynamics will be investigated first by means of the quantum version of LD theory, following Ref. [6]. This procedure gives access to the statistics of the trajectories of the system, painting a clear picture of its dynamical behavior. Following this, we complement these anayltical insights through the use of quantum jump Monte Carlo simulations FIG. 1. (Color online) The system we consider. Three spins are arranged on the vertices of an equilateral triangle. Spin 1 is coupled to a harmonic oscillator, represented by the springs, allowing it to move in the horizontal direction. The motion of spin 1 modulates the spin-spin interaction strengths and mediates an interaction between a collective spin degree of freedom and the mechanical bath. [24] that give access to a transparent physical interpretation of the processes occurring. The Hamiltonian describing the three-spin system interacting with a harmonic oscillator as described above can be writtenĤ =Ĥ s +Ĥ m +Ĥ s-m witĥ being the spin-chain Hamiltonian under uniform magnetic fields α and B α,Ĥ m = ω mb †b the harmonic oscillator Hamiltonian, andĤ s-m = −g(b † +b)σ 1 x (σ 2 x −σ 3 x ) the interaction Hamiltonian between the spin chain and the harmonic oscillator; i, j denotes a sum over nearest neighbors. We now move into the interaction picture with respect toĤ 0 =Ĥ s +Ĥ m , settingb →b e −iω m t andσ j x →σ j + e −iBt +σ j − e iBt , whereσ j ± = σ j x ± iσ j y is the spin-flip operator for the j th spin and the tilde distinguishes interaction-picture operators from Schrödingerpicture ones. Assuming that ω m = 2B and performing a rotating-wave approximation allows us to consider only the time-independent terms in the interaction-picture interaction HamiltonianH s-m , resulting inH s-m ≈ −g(bσ + +b †σ − ), where the collective operatorsσ ± =σ 1 ± (σ 2 ± −σ 3 ± ). Assuming that |g| sets the longest time-scale of the dynamics, we can adiabatically eliminate the harmonic-oscillator degree of freedom [24] and obtain a master equation for the reduced interaction-picture density matrix for the spin systemρ(t) with Γ = g 2 /κ. Here we have defined κ to be the damping rate for the harmonic oscillator andn the average number of excitations in the corresponding bath. We also havẽ where D ↑,↓ := e iĤ s tD ↑,↓ e −iĤ s t . We similarly defineσ ± by transformingσ ± and shall drop the label t when it can be understood from the context. By tracing out the motion of the damped harmonic oscillator, we have obtained an effective damping acting on a collective degree of freedom of the spin system through the jump operatorsσ ± . This collective damping is the source of the interesting behavior that we shall explore in the following.
As a first step in exploring the LD behavior of our system, we associate a counting process related to the flow of excitations into or out of the system with the collective jump operatorsσ ± . There are two counting processes K ± associated witĥ σ ± . K + counts excitations emitted into the bath and K − excitations absorbed from it. We then define an overall counting process K, which counts the net number of excitations emitted into the bath due to the collective spin flips K := K + − K − (in contrast to the total activity given by K + + K − [6]). Next, we can unravel the master equation of the reduced density matrix by projecting it onto a particular number of jump events, i.e., ∂ t ρ K = P K W[ρ] where P K is a projector over trajectories with K net jump events, and p K (t) = Tr{ρ K (t)} = Tr{P K ρ(t)} represents the probability to observe such a trajectory [13]. The moment-generating function associated to this probability p K (t) can be written [6]: with ρ s (t) = ∞ K=0 e −sK ρ K (t) the Laplace transform of the density matrix with respect to the net excitation exchanges K between the system and the bath; we call s the 'bias parameter.' The Laplace-transformed density matrix evolves according to the modified master equation In the long-time limit, LD theory applies and we can write Z(t, s) → e tθ(s) , where θ(s) represents the system's dynamical free energy [6,14]. Consequently, we have θ(s) = lim t→∞ ln Tr{ρ s } /t [25]. θ(s) is also given by the eigenvalue of W s := W + V s with the largest real part [6] (which can be shown to be real [26]).
Derivatives of this dynamical free energy with respect to s can be used to obtain the activity (net count rate) k(s) = −∂ s θ(s) of the system. This quantity is represented in Fig. 2 for different values of the bath populationn. As is clearly visible form the upper plot, for the value of the bias parameter s = 0 we have a non-analytic point in θ(s) for any value of n. The lower curves show that this point presents two distinct values of the activity k(s). Unbiased dynamics takes place at s = 0; from this we can conclude the existence of two dynamical phases [6,14]. These two phases have k(0 − ) > 0 and k(0 + ) = 0, i.e., one phase is active in the sense that the net rate of excitations exchanged between the bath and the system is nonzero whereas the other has an exact balance between excitations emitted and absorbed. The two terms composing the dissipative part of the master equation (4) act at different rates; we thus deduce that this balanced phase is inactive and no emission or absorption events occur. Focussing now on the active phase corresponding to k(0 − ) > 0, we also notice that the activity seems not to depend on the thermal population of the bathn, since all the curves converge to the same point as s → 0. This effect, seen for smalln, stems from the weak coupling between the bath and the system (Γ α).
A feature of Fig. 2 that is not seen when considering a "symmetric" dynamical observable such as the total activity K + + K − , as in Refs. [6,14,[27][28][29], is the second point of non-analyticity in θ(s) that occurs for s > 0 whenn > 0. This is related to a Gallavotti-Cohen symmetry [15] due to the driven nature of the system's dynamics [30]. Indeed, for s 0 = ln(1 + 1/n) we find θ(s) = θ(s 0 − s), which yields a fluctuation theorem of the form where we have defined p ∞ K := lim t→∞ p K (t), that relates the infinite-time probability of observing a trajectory with a net count of K to one with a net count of −K. This ratio approaches unity asn → ∞ in which case the rates for the collective jump process balance (s 0 → 0), such that k(0 ± ) = 0 [31]. Conversely, asn → 0, the probability for observing negative K goes to zero and the above ratio diverges.
To make the model more realistic we now add independent damping channels acting on each spin and explore the conse- quences of these channels on the dynamical behavior of the system. For simplicity, the single-spin baths are also taken to have occupationn and be coupled with the rate γ. In Eq. (4) we set W → W : stays unchanged from its definition Eq. (6) when γ/3 Γ and neglecting correlation effects between the various damping channels. D i ↑,↓ are the spin-flip Liouvillians for spin i. As Fig. 3 illustrates, introducing extra damping at the level of the individual spins has a strong effect on the dynamical properties of the system. Concentrating on the dynamical free energy, we see that individual spin damping will smoothen out the non-analyticity at s = 0. The same holds for the second non-analyticity at s > 0. This smoothing effect is also visible in the activity, which becomes well-defined everywhere. Meanwhile, as it is clearly visible in Fig. 3, for small values of the single-spin damping and low thermal population the activity remains approximately constant. Conversely, for high thermal population and strong single-spin damping, the activity can switch sign and become negative. Physically, this corresponds to the case where the single-spin damping channel upsets the balance, making it more likely that excitations enter (K − ) than leave the system (K + ) through the collective channel, leading to a thermally driven system.
The LD approach to dynamical phase transitions yields a transparent physical interpretation based on the statistics of ensembles of trajectories of the system. To explore these statistics, we now conduct an analysis based on a Monte Carlo wave-function (MCWF) simulation (also known as quantum jump Monte Carlo). The MCWF technique is a wellestablished method to simulate open system dynamics, like the one we are interested in, following the ideas set forth in Refs. [24,32]. Using this technique and starting off from a randomly-chosen initial condition, we simulate trajectories of Eq. (8) of jump events related to the operatorsσ ± (t). For each trajectory generated we estimate the activity k(0) by calculating the net rate of jump events. We consider a set of 2000 trajectories, each having approximatively 10 4 jump events, from which we obtain probability distributions for k(0), as represented in Fig. 4. This figure shows the probability distributions for thermal populations ofn = 0 and 5, with (upper panel) our without (lower panel) single-spin damping. Samples of the typical trajectories obtained corresponding to different parts of the distributions are shown in inset; each vertical line corresponds to a jump event, with upper ones representing an emission from the system to the bath (K + ) and lower ones the opposite (K − ). It is clear that when γ = 0 (upper panel), the probability distribution is bimodal, with one peak centred at an activity equal to k(0 + ) = 0 and the other at k(0 − ) > 0. In the former case (cf. inset) the corresponding trajectories have no jump events.
The second peak is centred about the same value for both values ofn, in agreement with what is expected from the middle panel of Fig. 2, where we saw that the activity is independent of the thermal population in the low-temperature limit. Corresponding trajectories are shown in the inset, and as expected demonstrate jumps associated only withσ − (i.e., the system losing excitations to the bath) forn = 0; forn = 5 jumps are observed in both directions. It can be see in Fig. 4 that this second peak of the distribution broadens when the temperature increases. All this yields the interpretation that the Liouville space accessible to the system consists of two disconnected subspaces, one active and one inactive, to which the two peaks in the activity distribution are related. Based on the fraction of K = 0 trajectories in our simulation, we can determine that 25% of the Liouville space is inactive.
Consider now the case with damping on the individual spins, corresponding to the one shown in Fig. 4 in the lower panel. By contrast to the γ = 0 case, we immediately see that the distribution becomes unimodal, with the mean activity decreasing by almost 25% compared to the active trajectories of the γ = 0 situation. This fraction corresponds to the fraction of inactive trajectories observed when γ = 0 (upper panel). The unimodal behavior and the mentioned reduction lead to the interpretation that the single-spin damping channel connects the two previously disconnected parts of the Liouville space. This interpretation is supported by the sample trajectories shown in inset of lower panel of Fig. 4, where one can observe the trajectory 'blinking,' i.e., spontaneously switching between active and inactive behavior. In analogy with what occurs with the total activity Ref. [14], the intermittency in trajectories is a form of "mesoscopic" (i.e., finite time) dynamical phase coexistence, consistent with the fact that the dynamical free-energy is analytic in this case, and the transi-  5. (Color online) Activity as a function of the thermal populationsn beyond the low-temperature limit. The case with (without) single-spin damping is shown in green (black). The solid curve corresponds to the large-deviation analysis and the data line to the Monte Carlo simulations. The error bars correspond to a fifth of the standard deviation. Other thann, the parameters are the same as in Fig. 2. tion therefore becomes a crossover. A simple interpretation of this behaviour can be deduced by looking at the eigenspace of H and the collective spin-flip operators: One can find a subspace of the Liouville space that is both inactive and isolated, in the sense thatĤ does not couple it to the rest of the space, and a subspace that is active. When γ 0, these two partitions are no longer isolated, and the system can switch dynamically between the active and inactive subspaces.
To understand how closely our MCWF results agree with the LD analysis, we present in Fig. 5 a quantitative comparison between the two. We plot the activity of the system at s = 0 as a function of the thermal populations, going beyond the low-temperature limit to which we restricted ourselves for most of the preceding discussion. This plot shows that the MCWF results (data points) agree very well with the results from the LD theory (solid curves). As visible in Fig. 4, increasingn tends to broaden the distribution of the activity. Consequently, the error bars shown in Fig. 5 grow quickly withn. As discussed previously, we clearly see that both curves tend to zero asn → ∞, where the rates of the two counting processes K ± balance such that the net rate of jump events is zero. Note also that the decrease in activity matches the predicted 25%, independently ofn.
We sum up by recalling our main results. We have explored a simple yet intriguing system consisting of three equidistant spins interacting pairwise, one of which moves in a harmonic trap. This motion gives rise to collective spin dynamics, which are dissipated through the mechanical decay channel. Adopting a large-deviation approach to analyse this system, we observe that its dynamics consists of two distinct dynamical phases, one active and one inactive, possessing different emission statistics. We observed that these two phases can be mixed by introducing damping on the individual spins. All our observations were confirmed through Monte Carlo simulations, which lend themselves to a natural interpretation in terms of ensembles of trajectories observed in repeated experimental runs. The system we explore is not overly complex but yields a surprisingly rich behavior.
This work has been supported by the Royal Commission for the Exhibition of 1851, the UK EPSRC (EP/G004579/1, EP/L005026/1 and EP/J009776/1), the John Templeton Foundation (grant ID 43467), ERC Grant ESCQUMA (Grant Agreement 335266), and the EU Collaborative Project Ther-MiQ (Grant Agreement 618074). Part of this work was supported by the COST Action MP1209 "Thermodynamics in the quantum regime."