Active energy transport and the role of symmetry breaking in microscopic power grid
Julian Huber Peter Rabl









Vienna Center for Quantum Science and Technology,
Atominstitut, TU Wien, 1040 Vienna, Austria

Vienna Center for Quantum Science and Technology,
Atominstitut, TU Wien, 1040 Vienna, Austria







We study the  transfer of energy through a network of coupled oscillators, which represents a minimal microscopic power grid connecting multiple active quantum machines. We evaluate the resulting energy currents in the macroscopic, the thermal and the quantum regime and describe how transport is affected by the competition between coherent and incoherent processes and nonlinear saturation effects.  Specifically, we show that the transfer of energy through such networks is strongly influenced by a non-equilibrium phase transition between a noise-dominated and a coherent transport regime. This transition is associated with the formation and breaking of spatial symmetries and is identified as a generic feature of active networks.  Therefore, these findings have  important practical consequences for the distribution of energy over coherent microwave, optical or phononic channels, in particular close to or at the quantum limit. 

 
 


Active energy transport and the role of symmetry breaking in microscopic power grids
    Peter Rabl
    December 30, 2023
====================================================================================



§ INTRODUCTION

Motivated by fundamental thermodynamical considerations as well as potential practical implications, there has recently been a growing interest in the performance of microscopic generators, engines or refrigerators, which may even be realized with single quantum systems <cit.>. However, while many theoretical  <cit.> and first experimental <cit.> studies of individual quantum machines have already been performed, there is still little known about interfacing multiple such devices. For example, how can the energy output of a microscopic generator be efficiently delivered to a microscopic engine and how will even larger networks of active quantum machines behave? Compared to conventional transmission lines for electric power, energy can be distributed at the microscale via highly coherent nanophotonic, microwave or micromechanical channels, while at the same time thermal and quantum fluctuations become important. Thus,  the flow of energy through such microscopic networks can differ strongly from the usual Ohm's or Fourier's law.




In this work we analyze the energy transfer between two quantum machines, which are connected through  a `power grid' of coupled harmonic oscillators, as depicted in Fig. <ref>.  Such oscillator chains have been considered in the past to study the passive transfer of heat between thermal reservoirs <cit.>, often with the goal to investigate the transition from ballistic transport to Fourier's law by adding additional dephasing mechanisms. Here we are interested in a different scenario where energy is injected at one end by a microscopic generator and extracted at the other end of the chain by a microscopic engine. 
Such quantum machines are active devices, meaning that they are (i) operated under non-equilibrium conditions and (ii) characterized by a maximal rate at which energy quanta can be emitted or absorbed. These properties make the problem of active energy transfer very distinct from the study of heat or electric transport between large passive reservoirs. By varying a single saturation parameter,  we can tune the degree of microscopicity of our network and investigate the resulting transport phenomena in the macroscopic, in the thermal and deep in the quantum regime.




Despite its conceptual simplicity, this setting already exhibits several surprising features, which will be relevant as well for the operation of more general  networks. Most importantly,  we identify two basic modes of transport, which are separated by a sharp phase transition and differ strongly in their noise characteristics. Such sharp transitions are unfamiliar in heat transport—neither in the ballistic nor in the Fourier regime—where mean currents and fluctuations vary smoothly as a function of the energy injection and extraction rates. The observed transition is accompanied by a breaking of the parity symmetry of the steady-state energy distribution along the channel. This  effect is closely related to the phenomenon of 𝒫𝒯-symmetry breaking <cit.> in systems with exactly balanced gain and loss. Importantly, here we show that the formation and breaking of spatial symmetries plays a much more fundamental  role for energy transport in coherent networks,  even in systems where such symmetries are not reflected in the underlying equations of motion. Therefore, these findings have direct practical consequences for various energy-distribution schemes at the microscopic level, but also reveal an interesting fundamental connection between symmetry-breaking, non-equilibrium phase transitions and the operation of networks of coupled quantum machines. 
 
 



§ MODEL

We consider a chain of N≥ 2 coupled harmonic oscillators, as schematically shown in Fig. <ref>. The oscillators have a frequency ω_0 and they are coupled to their neighbors with strength g. Energy is injected at the first site with a rate Γ_i and extracted at the other end with rate Γ_e. In addition, all oscillators are weakly coupled to a thermal environment at temperature T. In the frame rotating with  ω_0≫ g,Γ_i,e, the whole network is described by a master equation for the density operator ρ, 

    ρ̇=     - i/ħ [H_g,ρ] +Γ_i𝒟[A_1^†]ρ+Γ_e𝒟[A_N]ρ
       + ∑_ℓ=1^Nγ (N_ th+1) 𝒟[a_ℓ]ρ+γ N_ th𝒟[a_ℓ^†]ρ,

where a_ℓ (a_ℓ^†) are the annihilation (creation) operators for each oscillator and 𝒟[a] ρ≡  a ρ a^† - (a^† a ρ- ρ a^† a)/2.
In Eq. (<ref>), H_g= -ħ g/2∑_ℓ=1^N-1(a_ℓ^† a_ℓ+1 + H.c.) describes the coherent exchange of energy along the chain, while the second and the third term model the incoherent pump and dissipation processes, respectively. 

As mentioned in the introduction, in this work we are interested in active energy transport, where source and drain are represented by driven few-level quantum systems. However, to keep our analysis on a general level and to avoid details of specific implementations, we simply mimic the main characteristics of such microscopic generators and engines (providing gain, being saturable) by introducing in Eq. (<ref>) the nonlinear jump operators A_ℓ=1,N=f(a_ℓ^† a_ℓ) a_ℓ. Here the cutoff function f(x), where f(0)=1 and f(x≫ n_0)→0, accounts for the fact that both the injection as well as the extraction of energy saturate above a characteristic occupation number n_0. By changing this saturation parameter, we can tune the degree of microscopicity of the network without changing any other properties of the system. For concreteness, we will focus here on the cutoff function

    f(a^† a)=1/(1+a^† a/n_0),
 
which reproduces the saturation dependence of driven three-level generators and engines <cit.>, as depicted in Fig. <ref>. However, none of the central conclusions of this work depends on this assumption and different shapes of f(x) can be used to model other realizations of quantum machines discussed in the literature <cit.>.

Finally, the second line of Eq. (<ref>) describes the coupling of each oscillator to a local thermal bath, where N_ th=(e^ħω_0/k_BT-1)^-1 is the equilibrium occupation number and γ the damping rate, which we assume to be much smaller than Γ_i and Γ_e. Note that the use of local jump operators in Eq. (<ref>) is justified by the assumption that ω_0 is large compared to both the coherent intra-system coupling g and the dissipation rates <cit.>. Typical systems which can be used to implement this model include coupled nanomechanical resonators <cit.>,  linear chains of trapped ions <cit.> or arrays of coupled LC oscillators <cit.>. For all those platforms various techniques for engineering local gain and loss processes at the quantum level are already experimentally available <cit.>.  
  

For most parts of the following discussion we will be interested in the regime n_0≫1,  where Eq. (<ref>) can be mapped onto a Fokker-Planck equation for the Glauber-Sudarshan P-distribution P({α_ℓ},t) <cit.> (for details see  Appendix <ref>). This distribution can be sampled numerically by integrating the corresponding stochastic Ito equations for the amplitudes α_ℓ <cit.>,  

    α̇_̇1̇   =   Γ_i(α_1)-γ/2α_1 +i g/2α_2+√(D_ th+Γ_i(α_1))ξ_1(t),
    α̇_̇ℓ̇   =   - γ/2α_ℓ+i g/2( α_ℓ-1 + α_ℓ+1) + √(D_ th)ξ_ℓ(t),
    α̇_̇Ṅ   =   -Γ_e(α_N)+γ/2α_N + i g/2α_N-1+ √(D_ th)ξ_N(t).

Here  Γ_i,e(α)= Γ_i,e f^2(|α|^2) and D_ th=γ N_ th is the thermal diffusion rate. The ξ_ℓ(t) are white noise processes, which satisfy ξ_ℓ^*(t)ξ_ℓ'(t')=δ_ℓℓ'δ (t-t'). We are primarily interested in the steady-state energy current ⟨ J_ℓ⟩ = ig/2a_ℓ^† a_ℓ-1-a_ℓ-1^†  a_ℓ= g  Im⟨⟨α_ℓ-1^* α_ℓ⟩⟩, which can be obtained from the longtime average over many trajectories, denoted by ⟨⟨·⟩⟩. In the regime of interest, γ→ 0, the average current is approximately constant throughout the chain and we can drop the index ℓ.









§ ANOMALOUS ENERGY TRANSPORT
 
We first consider the macroscopic regime n_0≫ 1 and N_ th≪ n_0, where both thermal and quantum noise effects in Eqs. (<ref>)-(<ref>) can be neglected.  The steady state  is then described by a set of amplitudes α_ℓ^0 and in Fig. <ref>(a) we plot the corresponding current ⟨ J⟩ for N=10 sites. We see that transport in this system is very different from Ohm's law, but also from the ballistic flow of heat through a coupled chain of harmonic oscillators <cit.>. Overall, we find regimes of normal transport, where for fixed injection rate Γ_i the current increases with increasing extraction rate Γ_e, but also regimes of anomalous transport, where the opposite dependence is observed. 

For Γ_i<g there is a range of rates Γ_e where the current is completely stalled and only re-establishes at higher extraction rates. This counterintuitive behavior <cit.> can be traced back to the fact that within this parameter range all eigenvalues of the linear chain, i.e., when saturation effects are neglected, have a negative real part and the whole network is damped to zero. 
In all other parameter regimes the analysis of the linear chain predicts amplified solutions with a maximal gain rate that simply increases with increasing Γ_i, see Appendix <ref>. This behavior of the linear chain is not at all reflected in the stationary current shown in Fig. <ref>(a), 
which instead has a sharp maximum around Γ_e ≃Γ_i. For Γ_e=Γ_i the current then saturates at  ⟨ J⟩≃ g n_0 above Γ_i=Γ_e=4g as shown in Fig. <ref>(b). Note that for  γ→ 0 the current can exhibit sharp discontinuities near this symmetry line, where it  jumps abruptly within a range δΓ_e∼𝒪(γ).  At high rates, Γ_i/g>4,  also a bistable regime exists, where the current depends on the order in which the rates are switched on. However, in our analysis below we find that these fine-tuned features are washed out in the presence of noise and therefore they are less relevant for understanding transport in the microscopic regime.




§ SYMMETRY-BREAKING


In Fig. <ref>(c) we also plot the occupation numbers |α_ℓ^0|^2, i.e., the stationary distribution of the energy along the channel.  In contrast to conventional transport scenarios, where the energy distribution is flat or monotonically decreasing <cit.>, here the chain exhibits an alternating zig-zag structure. For  γ→ 0 and N even we obtain (see Appendix <ref>)

    |α^0_ℓ|^2=| A sin(k_0 ℓ)+ B cos (k_0 ℓ)|^2,

where k_0=π(N+2)/(2N+2) for  Γ_i>Γ_e and  k_0=π/2 for Γ_i<Γ_e. Eq. (<ref>) shows that the stationary current is carried by a single mode with wavevector k_0≈π/2, which is the mode supporting the highest current. However, since the saturable absorber can only extract a finite amount of energy per unit of time, most of the energy current is reflected at the extraction site and forms a standing wave.

While a standing-wave pattern is observed in all parameter regimes, the boundary conditions depend on the relation between Γ_i and Γ_e. For Γ_i>Γ_e the two ends of the chain have exactly the same amplitude, |α_1^0|^2≃ |α_N^0|^2≃ |A|^2 and |B/A|≪ 1. In contrast, for Γ_i<Γ_e the amplitude of the gain mode is much higher than the amplitude of the loss mode, |α_1^0|^2≫|α_N^0|^2. Therefore, for Γ_i>Γ_e the steady-state energy distribution of this network exhibits a left-right (parity) symmetry, which is broken above the transition point Γ_e≃Γ_i. Exactly at this point we obtain B≃ -i A and the transport becomes unidirectional, α^0_ℓ∼ e^ik_0ℓ. Note that also this behavior of the steady-state amplitudes cannot be derived by looking at the mode function of the most unstable mode of the linear chain. This mode always has the highest amplitude on site ℓ=1, such that gain is maximized. A more detailed derivation and discussion of the steady-state amplitudes is given in Appendix <ref>.

The breaking of a spatial symmetry in systems with gain and loss is reminiscent of the effect of 𝒫𝒯-symmetry breaking <cit.> in systems with equal gain and loss rates. 
Interestingly, in the current system such a symmetry is not present in the underlying equations of motions, since for  Γ_i≠Γ_e the oscillators at the injection and extraction sites evolve with very different rates. Additional numerical and analytical results summarized  in Appendix <ref> show that the emergence of a symmetric stationary phase also does not rely on the specific choice of the cutoff function f(x) and can be found even in situations where the gain and loss processes saturate at different amplitudes, n_0^(1)≠ n_0^(N). For such general cases, the symmetric phase is established as long as there is an amplitude α^0=α_1^0=α_N^0 such that energy conservation 
 
    Γ_i(α^0)-Γ_e(α^0)≃γ N/2
 
 can be satisfied. Only for larger chains this symmetry degrades when either Nγ > g or a finite amount of disorder prevents a free propagation of excitations (see Appendix <ref>).
Therefore, consistent with previous observations in specific two-mode systems <cit.>, we find that the emergence of steady-state symmetries and the breaking thereof is a generic  mechanism in active oscillator networks. In this context, the so-called 𝒫𝒯-symmetric configuration, Γ_i= Γ_e, appears naturally as the phase boundary, along which additional symmetry-breaking transitions can take place <cit.>.    







§ CURRENT FLUCTUATIONS

To understand the  consequences of this symmetry-breaking transition for microscopic transport, we consider now the thermal regime, n_0≫1 and  N_ th∼ n_0. In this case, quantum effects are still small, but noise from the environment can no longer be neglected and induces strong fluctuations of the current,  Δ J=√(⟨ J^2⟩-⟨ J⟩^2). In equilibrium, Δ J=g N_ th/√(2) and for Δ J/⟨ J⟩∼ N_ th/n_0>1 these thermal  fluctuations exceed the average currents discussed above. In Fig. <ref>(a) we consider this high-noise regime and plot ⟨ J⟩ and Δ J for N_ th/n_0=10, Γ_i/g=4 and for varying Γ_e. We see that in the symmetric phase  transport is indeed dominated by fluctuations, which even exceed the thermal level. This behavior changes abruptly after the transition point Γ_e≃Γ_i, beyond which a well-defined current below the thermal noise level is established. 
This transition is also clearly visible in the steady-state distributions of the individual oscillators, P_ℓ(α_ℓ), shown in Fig. <ref>(b). For Γ_i>Γ_e we observe strong fluctuations, but the distributions are still symmetric with respect to the center of the chain, i.e., P_ℓ≃ P_N-ℓ+1.  For Γ_e>Γ_i this symmetry is broken and fluctuations are strongly suppressed.
 


The striking difference in the current noise can be related to an equivalent change in the response of the network. In Fig. <ref>(c) we plot the relaxation time τ_r,  i.e., the time it takes for the amplitude α_1 to relax back to its steady-state value after a small perturbation has been applied to site ℓ=1. For details about the numerical procedure that has been used to determine τ_r, see Appendix <ref>.  In the symmetric phase this time constant is approximately  independent of Γ_i, Γ_e and N. It is essentially determined by the bare damping rate, τ_r∼γ^-1, and diverges in the limit γ→0. In the symmetry-broken phase a much faster response, τ_r∼ O(Γ_e^-1)∼ N^2 is observed. At the transition point the relaxation time diverges as τ_r∼ (Γ_e-Γ_i)^-ξ, where we find ξ≃ 1.3 from numerical simulations. This behavior is very different from a laser or from other non-equilibrium phase transitions, where the relaxation time diverges only at the transition point, but is finite and of similar magnitude in both phases <cit.>. As shown in Fig. <ref>(c), in the current system the relaxation time diverges (in the limit γ→ 0) within the whole symmetric phase. 

To provide a connection between the symmetry of P({α_ℓ}) and the current noise, it is useful to consider the mean damping rate Γ̅=Γ_e(α_N)-Γ_i(α_1) <cit.>, i.e., the average difference between energy injection and extraction rates. Due to the symmetry of the marginal distributions, this rate is vanishing small in the symmetric phase, Γ̅∼ O(γ) [see Fig. <ref>(d)]. By breaking this symmetry, a finite value Γ̅≫γ is established for Γ_e>Γ_i. This then leads—on average—to an efficient cooling of fluctuations and the possibility for subthermal energy transport. Again this behavior shows a close analogy to conventional 𝒫𝒯-symmetric systems <cit.>. In such systems the breaking of the parity symmetry of the eigenstates of a non-Hermitain matrix is accompanied by a transition from real to imaginary eigenvalues, i.e. a transition from a purely oscillatory to an exponentially damped or amplified dynamics <cit.>.  The order parameter Γ̅ generalizes this effect to steady-state distributions of nonlinear gain-loss systems <cit.>, where the conventional definition of 𝒫𝒯-symmetry breaking is no longer meaningful. Note that  the scaling of Γ̅ near the transition point is not related to the branching of eigenvalues near an exceptional point <cit.> and depends solely on the saturation function f(a^† a). This is illustrated by additional numerical examples presented in Appendix <ref>.   Importantly, all these examples show that also the characteristic cancellation of the average dissipation rate is a much more general effect and occurs as well in gain-loss systems where no symmetry is present on a fundamental level.








§ QUANTUM NOISE LIMIT
 

 


From Eq. (<ref>) we see that even for N_ th≈ 0, the network is still affected by quantum noise ∼√(Γ_i(α_1))ξ_1(t). In the regime N_q=Γ_i/γ>N_ th, this noise dominates over thermal fluctuations and represents a fundamental limitation for energy transport deep in the quantum regime, n_0∼ O(1). Fig. <ref>(a) shows that for n_0≫ 1 the sharp transition between a noisy and a coherent transport regime still prevails, even for N_ th=0. As the saturation number n_0 is lowered, the relative level of fluctuations increases, develops a peak at the transition point and becomes much more pronounced also in the symmetry-broken phase. Note that for small n_0≲ 10 the mapping of the master equation onto a Fokker-Planck equation is no longer valid and a full simulation of Eq. (<ref>) must be performed (see Appendix <ref>). Therefore, due to the large Hilbert space and large separation of time scales involved in such simulations, the results in Fig. <ref> are restricted to N=2 oscillators. 







Access to the full density operator also allows us to investigate true non-classical quantities, such as the entanglement established between the injection and extraction sites. As shown in  Fig. <ref>(b) for different n_0=1,2,5, a significant amount of entanglement exists for Γ_e<g, it then vanishes in the rest of the symmetric phase, and peaks again right after the transition point. Therefore, this plot reveals an additional substructure, which is not reflected in the mean current or its fluctuations.  This entanglement between source and drain can  be relevant for thermodynamical considerations, where not only the flow of energy, but also changes in entropy through mutual (quantum) correlations must be taken into account.  Note, however, that for a more detailed study of entanglement it is necessary to go beyond our simply model and explicitly include specific implementations of quantum generators and engines in the dynamics.







§ DISCUSSION AND CONCLUSION



In summary, we have shown that the transfer of  energy between two active quantum machines can be very different from what one would intuitively expect based on Ohm's or Fourier's law.  Most importantly, we have identified the emergence and breaking of spatial symmetries as a generic feature of such networks, which drastically affects their dynamical response and noise properties. This aspect is of particular relevance in the microscopic regime where thermal and quantum fluctuations are unavoidable and play a dominant role. 


The transport effects analyzed here in detail for a single channel will be highly relevant as well for understanding the flow of energy through more complex networks with multiple active sites. To illustrate this point, we consider in Fig. <ref>(a) a small multi-port network where energy is extracted at two sites 7 and 9 with rates  Γ_e^(7) and Γ_e^(9), respectively. Fig. <ref>(b) shows that although Γ_e^(9)≫Γ_e^(7), only a residual thermal current is initially flowing from site 1 to site 9. This counterintuitive behavior can be explained by the fact that a symmetric standing wave is formed between sites 1 and 7, which results in a vanishing amplitude α_4≈ 0 at the crossing site. Once Γ_e^(7) is increased above the value of about Γ_i^(1), the symmetry breaks and α_4≠0 now supports a large current flowing to site 9. As a result, we obtain a transistor-like behavior, where a small increase of losses in one site leads to a sudden increase of the energy current through another part of the network. Simultaneously, there are sharp jumps in the level of current fluctuations, in analogy to what we have found above for the 1D chain.

This brief outlook already shows that the combination of interference, nonlinear symmetry-breaking effects and the prominent role of noise makes the operation of microscopic power grids a very rich and complex problem, which is still little understood. The current analysis reveals the important part in this problem that is played by quantum fluctuations as a fundamental source of noise as well as by the topology of the network, which determines whether these fluctuations are enhanced or suppressed.  



§ ACKNOWLEDGEMENT
 We thank Alexander Carmele, Wolfgang Niedenzu, Stefan Rotter and Henning Schomerus for stimulating discussions. This work was supported by the Austrian Science Fund (FWF) through the SFB FoQuS, Grant No. F40, the START Grant No. Y 591-N16, and the DK CoQuS, Grant No. W 1210.   J.H. is a recipient of a DOC Fellowship of the Austrian Academy of Sciences (ÖAW).






§ FOKKER-PLANCK EQUATION

In the semiclassical regime n_0≫ 1 the cutoff function f(x) in Eq. (<ref>) varies slowly on the scale of individual excitations and the master equation can be mapped onto a Fokker-Planck equation for the Glauber-Sudarshan P-representation <cit.>. This distribution function is defined by

    ρ= ∫∏_ℓ d^2α_ℓ  P({α_ℓ}) |{α_ℓ}⟩⟨{α_ℓ}|,

where |{α_ℓ}⟩ denotes a multi-component coherent state. By using the usual substitution rules <cit.>

    a_ℓρ  →  α_ℓ P,      a^†_ℓρ  →   (α_ℓ^*-∂/∂α_ℓ)  P,
    ρ a^†_ℓ  →  α^*_ℓ P,     ρ a_ℓ  →   (α_ℓ-∂/∂α^*_ℓ)  P,
 
we can convert Eq. (<ref>) for the density operator ρ into a partial differential equation for P({α_ℓ}). We obtain

    ∂ P/∂ t =  .∂ P/∂ t|_ lin + .∂ P/∂ t|_ nl,

where the first term,

    .∂ P/∂ t|_ lin  =1/2[-i g  ( ∑_ℓ=1^N-1∂/∂α_ℓα_ℓ+1.  +∑_ℓ=2^N∂/∂α_ℓα_ℓ-1 )  
     +. γ∑_ℓ=1^N(∂/∂α_ℓα_ℓ + N_ th∂^2/∂α_ℓ∂α_ℓ^*) + c. c. ] P,

describes the linear chain and already has the form of a Fokker-Planck equation. 

The nonlinear dissipative terms in Eq. (<ref>) at the ends of the chain translate into higher order derivatives for the P-distribution and additional approximations are required. To do so we first use the substitution rules from above to translate the action of f(a^†a) on the density operator into a differential operator acting on a coherent state,

    f(a^†a) |α⟩⟨α|  →  ∑_m=0^∞f̅_mα^m (α^* + ∂/∂α)^m |α⟩⟨α|.

Here the coefficients f̅_m follow from an expansion of the  operator f(a^†a) into a normally ordered series

    f(a^†a)=∑_m=0^∞f̅_m (a^†)^m a^m.

By using the binomial theorem

    (α^* + ∂/∂α)^m= ∑_k=0^mmk (α^*)^m-k∂^k/∂α^k
     = ∑_k=0^m1/k!∂^k/∂α^*k (α^*m) ∂^k/∂α^k,

and integrating by parts we obtain the following substitution (omitting the site index) 

    f(a^†a) ρ  →  ∑_m=0^∞∑_k=0^mf̅_mα^m (-1)^k /k!∂^k/∂α^*k (α^*m) ∂^k/∂α^k P(α)
     =∑_k=0^∞(-1)^k /k!∂^k/∂α^*k[f̅(α,α^*)] ∂^k/∂α^k P(α),

where f̅(α,α^*)=αf(a^† a)α. Since f(a^† a) is a function of a/√(n_0) and a^†/√(n_0), the derivatives of f̅(α,α^*) scale as ∂^k/∂α^*kf̅(α,α^*) ∝ n_0^-k/2. Therefore, in the limit n_0→∞, we can neglect all derivatives and approximate 

    f(a^†a) ρ  →  f̅(α,α^*) P(α,α^*)+ 𝒪 (1/√(n_0)).

Note that the definition of f̅(α,α^*) is based on the normally ordered series expansion and in general  f̅(α,α^*)≠ f(|α|^2). Therefore, in our derivation we make a second approximation and neglect this difference, i.e.,

    f̅(α,α^*)=α1/(1+a^† a/n_0)α≈1/(1+α^2/n_0).

To show the validity of this approximation, we compare the function α1/(1+a^† a/n_0)α=e^-α^2 (-α)^-n_0[Γ (n_0,0)-Γ (n_0,-α^2)] <cit.> with  the approximate form 1/(1+α^2/n_0). Here Γ(n,x) denotes the incomplete Gamma function. Even deep in the quantum regime, n_0≈1,  these two expressions agree up to a few percent and become essentially identical for n_0≳ 10. Therefore, we conclude that the main approximation in the derivation of our semiclassical Fokker-Planck equation arises from neglecting higher order derivatives in Eq. (<ref>).  

Based on these considerations we obtain the following approximate substitution rules

    A ρ  →  α/(1+α^2/n_0)  P(α,α^*), 
     A^†ρ  →   (α^*-∂/∂α) 1/(1+α^2/n_0) P(α,α^*),

and analogous relations for ρ A and ρ A^†. All together we then obtain 

    .∂ P/∂ t|_ nl 
    =1/2[ -∂/∂α_1Γ_i(α_1) α_1
    + ∂^2/∂α_1∂α_1^*Γ_i(α_1) 
    . 
    . + ∂/∂α_NΓ_e(α_N) α_N+ c.c.
    ] P,

where Γ_i,e(α)=Γ_i,e f^2(|α|^2). After this approximation, the resulting Fokker-Planck equation (<ref>) can be mapped onto the set of stochastic differential equations (3)-(5) <cit.>. 





§ LINEAR CHAIN


In Fig. <ref>(a) we plot the largest real part of all the eigenvalues obtained from the dynamical matrix of a  linear chain where Γ_i,e(α)=Γ_i,e. As long as all eigenvalues have a negative real part, the chain is damped to zero. This only occurs in the `stalled' phase where Γ_i<g and Γ_i ≤Γ_e<g^2/Γ_i. Otherwise, we see that the structure of the current plotted in Fig. 2(a) is not at all reflected in the eigenvalue structure of the linear chain.







§ STEADY STATE AMPLITUDES


In the limit n_0≫ 1 and N_ th/n_0→ 0, the stochastic terms in Eqs. (<ref>)-(<ref>) can be neglected and we obtain a set of ordinary differential equations with steady-state amplitudes α_ℓ^0. To obtain analytic insights about the steady state of the chain in this regime, we consider in the following the slightly simplified scenario, where only sites ℓ=1 and ℓ=N are affected by the bare decay γ≪ g, while all the other oscillators evolve coherently. 


We are interested in the long-time dynamics of the chain and make the following ansatz for the amplitudes

    α^0_ℓ(t)=√(n_0) e^-i ω t[A sin(k_0ℓ)+B cos (k_0ℓ)],

where A,B ∈ℂ and ω=  g cos(k_0). With this ansatz the current between two sites is

    ⟨ J_ℓ⟩  = g  Im{  (α^0_ℓ-1)^* α^0_ℓ}= g n_0 Im{A B^* }sin(k_0).

To obtain a steady state configuration that maximizes the energy transfer, we look for solutions with k_0 equal or close to π/2. By writing k_0=π/2+δ we obtain the equation

    [Γ_i/2/(1+A cos(δ)-B sin(δ)^2)^2-γ/2]
    × (A cos(δ)-B sin(δ)) - i g/2 B=0,

from the equation of motion for α_1. Similarly, from the equation of motion for α_N we obtain








    [-Γ_e/2/(1+A sin(δ N)+B cos(δ N)^2)^2-γ/2] 
    × (A sin(δ N)+B cos(δ N))
    
    -i g/2 ( A cos(δ (N+1))-B sin(δ(N+1))=0,

for the case where N is even. To proceed with our analysis we must distinguish between the symmetric (Γ_i>Γ_e) and the symmetry-broken regime (Γ_e>Γ_i) and between an even and an odd number of oscillators.  In the following we will only carry out the analysis for an even number of oscillators, however the analysis for N odd can be done in an equivalent manner.







 §.§ Symmetry-broken phase
	
We first consider the regime Γ_e>Γ_i and N even. In this case the choice k_0=π/2 results in ω=0 and a symmetry-broken solution for the amplitudes, |α_1^0|≫ |α_N^0|. 	The remaining parameters A and B are determined by the two coupled equations


    (Γ_i/(1+A^2)^2-γ) A - i g B=0,
     (-Γ_e/(1+B^2)^2-γ) B - i g A=0.

These equations have a solution for Γ_e≥Γ_i, but not for Γ_i>Γ_e. Although these equations can still be solved analytically, the results are already quite involved. However, sufficiently deep in the symmetry-broken phase we can neglect the bare decay γ and approximate Γ_e(B)≈Γ_e. We then obtain

    A^2≃√(Γ_iΓ_e/g^2)-1,      B≃ -i  g/Γ_e  A,

and the current

    ⟨ J⟩≃g^2 n_0/Γ_e  (√(Γ_iΓ_e/g^2)-1 ).







 §.§ Symmetric phase

For N even and  Γ_i>Γ_e, the choice k_0=π/2 would results in an asymmetric steady state and also the resulting equations for A and B do not have a solution for Γ_i>Γ_e. To recover a symmetric solution with a maximal current we choose δ=π/(2(N+1)). In this case the chain undergoes persistent oscillations with frequency ω=g sin(δ). By defining Ã=A cos(δ) and using the approximation Bsin(δ) ≈ 0,  the resulting equations simplify to

    (Γ_i/(1+Ã^2)^2-γ) Ã - i g B=0,
     (- Γ_e/(1+Ã^2)^2-γ) Ã+i g B=0.

Therefore, we obtain the amplitudes

    Ã^2=√(Γ_i-Γ_e/2 γ)-1,      B=-i  γ/gΓ_i+Γ_e/Γ_i-Γ_eÃ,

and, since sin(k_0)=cos(δ), the current

    ⟨ J⟩ =n_0γΓ_i+Γ_e/Γ_i-Γ_e (√(Γ_i-Γ_e/2 γ)-1 ).





 §.§ Symmetry breaking transition
	





Near the transition point we find |A|≃ |B|. More precisely,  from the solution in the symmetric regime we see that B=-iA, at a value of

    Γ_e^*=Γ_i g-γ/g+γ≈Γ_i- 2Γ_i /gγ.

Near this parameter the standing wave turns into a running wave α_ℓ^0∼ e^ik_0ℓ and the current is close to maximum and scales as J_ max∝√(Γ_i).
Although the symmetric solution exists up to  Γ_e^**=Γ_i-2γ, the stability analysis reveals that for Γ_i>g the symmetric solution becomes unstable before, at around Γ≃Γ_e^*. In the regime of interest, γ/g→ 0, these differences become negligible and the transition is simply given by Γ_e= Γ_i. The dependence of the coefficients A and B around the transition point is shown in Fig. <ref>(b). 







 §.§ Damping of all oscillators

The results derived so far for a chain without damping of the oscillators in the middle agree in essence with the results obtained for two coupled oscillators (see also Ref. <cit.>). However, 
while in the symmetry-broken phase the bare damping γ has a negligible effect, it determines the value of the current in the symmetric phase. In this regime it is thus important to analyze the steady state also for the full system, where all oscillators are weakly damped. In this case the equation 

    α̇_ℓ = -γ/2α_ℓ +ig/2(α_ℓ-1+α_ℓ+1),

cannot be fulfilled by the ansatz (<ref>). However, for γ/g≪1 the correction are small and we can still use this ansatz with the same k_0 as above as a first approximation. For simplicity we focus on N odd where k_0=π/2. Then, summing the equations of motion for every other site we obtain

    ∑_ℓ=1^(N-1)/2 (-1)^ℓ+1α̇_2ℓ-1 =
    
       √(n_0) (Γ_i(A) A-γ∑_ℓ=1^(N-1)/2A- i g B )=0

and for the last site,

    α̇_N= √(n_0)[(-Γ_e(A)-γ) A+ i g B ]=0.

From this set of equations we obtain the amplitudes

    A^2=√(2(Γ_i-Γ_e)/γ (N+1))-1,    B=-i γ/2 gΓ_e (N-1)+2 Γ_i/Γ_i-Γ_e A,

and the current

    J=n_0γ/2Γ_e (N-1)+2 Γ_i/Γ_i-Γ_e (√(2(Γ_i-Γ_e)/γ (N+1))-1 ).

Although this result was derived for N odd, it is also a good approximation for N even when N>2.


Note that near Γ_e^* we obtain a single traveling wave. To account first order corrections due to a finite decay γ/g≪1, we can generalize the ansatz to α_ℓ^0∼ e^ik_0ℓe^-κℓ. From Eq. (<ref>) we then obtain κ=γ/(2g). Therefore, all our analytic estimates will remain valid as long as Nγ≪1, although numerical simulations show that most of the qualitative features survive at much larger decay rates.





§ UNIVERSALITY OF THE SYMMETRY-BREAKING TRANSITION
	
For all the results discussed in the main text we have assumed a specific cutoff function and the same saturation occupation number n_0 for the gain and the loss mechanism. While the precise quantitative findings will of course depend on these assumptions, we will now demonstrate with several other examples that the essential qualitative features of the symmetry-breaking transition do not depend on these details.  



 §.§ Different gain/loss mechanisms

As there are many ways to engineer gain and loss, we first show that our findings do not depend on the precise form of the saturation function f(x). In Fig. <ref> we consider the example of a cutoff function f(a^† a)=1/(1+a^† a/n_0)^ν/2, where we have assumed ν=1 to model a system with a weaker saturation dependence. This case corresponds, for example, to the saturation of a regular laser. Again we see the characteristic  structure of the current with a maximum at Γ_e≃Γ_i and that this maximum is associated with a transition between a symmetric and a symmetry-broken energy distribution.  For Γ_i≳  3.4g we obtain a region, where the current does not have a precise value and the whole chain settles into a limit cycle. Such a behavior has previously been predicted for a 𝒫𝒯-symmetric system, Γ_i=Γ_e, where ν=1 has been identified as a special case, where no real symmetry-breaking occurs <cit.>.  However, in the presence of thermal noise [see Fig. <ref>(d)] these limit cycles are no longer visible and qualitatively we obtain the same transition between a noise-dominated and a coherent transport regime as in the main part of the paper. The same behavior is also found for stronger nonlinearities, ν=3, and other saturation functions with different functional dependencies.  
 
 




















 §.§ Different saturation numbers

To further illustrate that the physical effects discussed in this work are very generic, we now return to the cutoff function given in Eq. (<ref>), but consider the case where the gain and the loss oscillator saturate at different amplitudes, i.e., n_0^(1)≠ n_0^(N). The resulting mean currents and fluctuations are shown in Fig. <ref>.  We see that also in this case all the qualitative features of the symmetry-breaking phase transition remain unaffected, except that the transition point is now shifted from Γ_i=Γ_e 

to Γ_i≃Γ_e(n_0^(N)/n_0^(1))^2.






Importantly, this example shows that even when different saturation mechanisms for energy injection and extraction are considered, there is still an emergent symmetric phase, which is characterized by an almost complete cancellation of the average dissipation rate Γ̅=⟨⟨Γ_e(α_N)-Γ_i(α_1)⟩⟩.

















 §.§ Conditions for a symmetric phase

To obtain a more general result for the symmetry-breaking point, we derive a minimal condition under which a symmetric phase can exist. This condition follows from the fact that in the steady state the total energy of the system must be conserved. This means that the absorbed and dissipated energy must be the same, or

    Γ_i(α^0_1) |α_1|^2 = Γ_e(α^0_N) |α^0_N|^2 +  γ∑_ℓ=1^N |α^0_ℓ|^2.

For a symmetric state, where |α_1^0|=|α_N^0|=α^0 and |α^0_ℓ|^2=η_ℓ |α^0|^2, we obtain

    Γ_i(α^0) -  Γ_e(α^0) =  γ𝒩,    𝒩= ∑_ℓ=1^N η_ℓ.

Here 𝒩=2 for N=2, 𝒩≃ (N+1)/2 for N odd and in general 𝒩≈ N/2 for N≫1. For identical cutoff functions this condition can always be satisfied by increasing the value of α^0, as long as Γ_e< Γ_i-γ𝒩. For non-identical saturation parameters, n_0^(1)≠ n_0^(N), and by approximating f(α^0)≃ n_0^2/|α^0|^4, this argument also explains the shift of the transition point discussed above. Thus, Eq. (<ref>) provides a simple minimal condition for the existence of a symmetric phase. Note, however, that for larger systems one find that for Γ_i>g symmetry breaking already  occurs closer to the point where  A=B. For example, for n_0^(1)= n_0^(N) we find the transition point approximately at

    Γ_e^*=2 Γ_i (g-γ)/2g+γ (N-1).

As long as γ N ≪ g, this result does not considerable change by changing the system size and for all results presented in the main text the transition point derived from Eq. (<ref>) is a sufficient approximation.





§ DISORDER




For all our results presented in the main part of the paper we have considered chains of oscillators with identical frequencies ω_ℓ=ω_0. To understand the robustness of the observed effects with respect to small frequency variations, which will be unavoidable in any real system, we numerically simulate the steady state of a chain of N=10 oscillators with frequencies  ω_ℓ=ω_0+Δ_ℓ. Here the random frequency offsets for each site are chosen from a uniform distribution Δ_ℓ∈[-σ_Δ,σ_Δ].

In Fig. <ref> (a) and (b) we plot the disorder-averaged steady-state occupation numbers |α^0_ℓ|^2 for each of the oscillators and the current for a few disorder realizations.  We find that for σ_Δ< 0.1 g, the steady-state amplitudes reproduce almost perfectly the alternating structure predicted for the ideal case, Δ_ℓ=0.  In this regime also the current exhibits the characteristic peak structure for each individual disorder realization and is hardly affected for parameters away from the transition point. This shows that all the effects discussed in the main part of this work are insensitive  to a small amount of disorder. For 0.1 <σ_Δ/g< 0.3, the amplitudes still follow more or less a zig-zag structure, while for σ_Δ/g>0.3 the energy distribution is completely different from the non-detuned case and most of the energy gets localized around the gain mode. 








§ NUMERICAL SIMULATIONS



For the numerical integration of the stochastic equations (<ref>)-(<ref>) we have used the Euler Maruyama method with a time step of Δ t=10^-4/g. For all the main plots, the stochastic equations have been vectorized and n_ traj=50 trajectories have been evolved simultaneously. After a time t=5000/g ≈ 5 τ_r, we have sampled 3000 points per trajectory, separated by  17000 time steps, to get the steady-state distribution.  

The relaxation time in Fig. 3(c) was obtained in the absence of noise, by first determining the steady state amplitudes with high accuracy. This was implemented by a 4-th order Runge-Kutta algorithm with an accuracy of 10^-11.  Then the amplitude of the gain oscillator was changed by an amount δα_1=1/10 and the simulation was continued until the system has relaxed again. From the time difference Δ t_r between the points where the remaining occupation difference of the first oscillator was δα_1^2=10^-5 and δα_1^2=10^-8, we calculated the relaxation time as τ_r=Δ t_r /ln(10^3). Note that in Fig. 3(c) the relaxation rate exhibits a peak in a very small region around the transition point Γ_e=Γ_i, where we find almost no relaxation.  In this regime the numerically extracted values for τ_r depend on very fine details and are no longer meaningful.  


For the cases n_0=1,2,5 in Fig. 4(a) and (b) we have used a stochastic quantum wavefunction method <cit.> to simulate the full master equation (<ref>). The results for n_0=1 were independently verified by calculating directly the steady-state density operator for a system of two coupled oscillators with n_ basis=30 basis states per oscillator. For the quantum trajectory simulations the evolution under the effective non-hermitian Hamiltonian has been implemented by the time evolution operator U=e^-i H_ effΔ t with Δ t= 2× 10^-3g^-1, which is computed once at the beginning of the trajectory. After random times quantum jumps occur and the state gets renormalized. After t=10000/g≈ 10 τ_r, when the system has reached the steady state, we have sampled the state after every 800 time steps for 990000 times to obtain the steady state density matrix. In Fig. 4(a) we used n_ basis=70 and n_ basis=150 states per oscillator for n_0=1 and n_0=2, respectively, while in Fig. 4(b) we used  n_ basis=30,60,100 states per oscillator for n_0=1,2,3. The entanglement negativity was obtained by calculating the 4 n_ basis  lowest eigenvalues after partial transpose.


























































































































































99





 Gelbwaser2015 D. Gelbwaser-Klimovsky, W. Niedenzu, and G. Kurizki, Thermodynamics of quantum systems under dynamical control, Adv. At. Mol. Op. Phys. 64, 329 (2015). 

Vinjanampathy2016 S. Vinjanampathy and J. Anders, Quantum thermodynamics, Contemp. Phys. 57, 1 (2016).

Goold2016 J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, The role of quantum information in thermodynamics — a topical review, J. Phys. A: Math. Theor. 49, 143001 (2016).    
 
Alicki2018 R. Alicki and  R. Kosloff, Introduction to Quantum Thermodynamics: History and Prospects, arXiv:1801.08314 (2018).

Millen2016 J. Millen and  A. Xuereb, Perspective: Quantum Thermodynamics, New J. Phys. 18, 011002 (2016).    

   
Scovil1959 H. E. D. Scovil and E. O. Schulz-DuBois, Three-Level Masers as Heat Engines, Phys. Rev. Lett. 2, 262 (1959).




Geva1994 E. Geva and R. Kosloff, Three-level quantum amplifier as a heat engine: A study in finite-time thermodynamics, Phys. Rev. E 49, 3903 (1994).

Scully2004 M. O. Scully, M. S. Zubairy, G. S. Agarwal, and H. Walther, Extracting Work from a Single Heat Bath via Vanishing Quantum Coherence, Science 299, 862 (2003).

Kieu2004 T. D. Kieu, The Second Law, Maxwell's Demon, and Work Derivable from Quantum Heat Engines, Phys. Rev. Lett. 93, 140403 (2004).

Quan2007 H. T. Quan, Y. Liu, C. P. Sun, and F. Nori, Quantum thermodynamic cycles and quantum heat engines, Phys. Rev. E 76, 031105 (2007).

Linden2010 N. Linden, S. Popescu, and P. Skrzypczyk, How Small Can Thermal Machines Be? The Smallest Possible Refrigerator, Phys. Rev. Lett. 105, 130401 (2010).

Abah2012 O. Abah, J. Ronagel, G. Jacob, S. Deffner, F. Schmidt-Kaler, K. Singer, and E. Lutz, Single-Ion Heat Engine at Maximum Power, Phys. Rev. Lett. 109, 203006 (2012).

Gelbwaser2013 D. Gelbwaser-Klimovsky, R. Alicki, and G. Kurizki, Minimal universal quantum heat machine, Phys. Rev. E 87, 012140 (2013).

Gallego2014 R. Gallego, A. Riera, and J. Eisert, Thermal machines beyond the weak coupling regime, New J. Phys. 16, 125009 (2014).

Zhang2014 K. Zhang, F. Bariani, and P. Meystre, Quantum Optomechanical Heat Engine, Phys. Rev. Lett. 112, 150602 (2014).

Bergenfeldt2014 C. Bergenfeldt, P. Samuelsson, B. Sothmann, C. Flindt, and M. Büttiker, Hybrid Microwave-Cavity Heat Engine, Phys. Rev. Lett. 112, 076803 (2014).

Brunelli2015 M. Brunelli, A. Xuereb, A. Ferraro, G. De Chiara, N. Kiesel, and M. Paternostro, Out-of-equilibrium thermodynamics of quantum optomechanical systems, New J. Phys. 17, 035016 (2015).

Elourad2015 C. Elouard, M. Richard, and A. Auffeves, Reversible work extraction in a hybrid opto-mechanical system, New J. Phys. 17, 055018 (2015).

Dechant2015 A. Dechant, N. Kiesel, and E. Lutz, All-Optical Nanomechanical Heat Engine, Phys. Rev. Lett. 114, 183602 (2015).

Mari2015 A. Mari, A. Farace, and V. Giovannetti, Quantum optomechanical piston engines powered by heat, J. Phys. B 48, 175501 (2015). 

Li2017  S.-W.  Li,  M. B.  Kim,  G. S.  Agarwal,  and  M. O. Scully,  Quantum statistics of a single-atom heat engine, Phys. Rev. A 96, 063806 (2017).


Hugel2002 T. Hugel, N. B. Holland, A. Cattani, L. Moroder, M. Seitz, and H. E. Gaub, Single-molecule optomechanical cycle, Science 296, 1103 (2002).

Steeneken2011 P. G. Steeneken, K. Le Phan, M. J. Goossens, G. E. J. Koops, G. J. A. M. Brom, C. van der Avoort, and J. T. M. van Beek, Piezoresistive heat engine and refrigerator, Nature Phys. 7, 354 (2011).

Blickle2012 V. Blickle and C. Bechinger, Realization of a micrometre-sized stochastic heat engine, Nature Phys. 8, 143 (2012).

Brantut2013 J.-P. Brantut, C. Grenier, J. Meineke, D. Stadler, S. Krinner, C. Kollath, T. Esslinger, and A. Georges, A Thermoelectric Heat Engine with Ultracold Atoms, Science 342, 713 (2013).

Thierschmann2015 H. Thierschmann, R. Sanchez, B. Sothmann, F. Arnold, C. Heyn, W. Hansen, H. Buhmann, and L. W. Molenkamp, Three-terminal energy harvester with coupled quantum dots, Nat. Nanotechnol. 10, 854 (2015).


Rossnagel2016 J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, A single-atom heat engine, Science 352, 325 (2016).

Schmidt2018 F. Schmidt, A. Magazza, A. Callegari, L. Biancofiore, F. Cichos, and G. Volpe, Microscopic Engine Powered by Critical Demixing, Phys. Rev. Lett. 120, 068004 (2018).

Klaers2017 J. Klaers, S. Faelt, A. Imamoglu, and E. Togan, Squeezed thermal reservoirs as a resource for a nano-mechanical engine beyond the Carnot limit, Phys. Rev. X 7, 031044 (2017).



Rieder1967 Z. Rieder, J. L. Lebowitz, and E. Lieb, Heat transport through lattices of quantum harmonic oscillators in arbitrary dimensions, J. Math. Phys. 8, 1073 (1967).

Lin2011 G.-D. Lin and L.-M. Duan, Equilibration and temperature distribution in a driven ion chain, New J. Phys. 13 075015 (2011).

Asadian2013 A. Asadian, D. Manzano, M. Tiersch, and H. J. Briegel, Heat transport through lattices of quantum harmonic oscillators in arbitrary dimensions, Phys. Rev. E 87, 012109 (2013). 

Bermudez2013 A. Bermudez, M. Bruderer, and M. B.  Plenio, Controlling and measuring quantum transport of heat in trapped-ion crystals, Phys. Rev. Lett. 111, 040601 (2013).

Motz2017 T. Motz, J. Ankerhold, and J. T. Stockburger, Currents and fluctuations of quantum heat transport in harmonic chains, New. J. Phys. 19, 053013 (2017).



Bender1998 C. M. Bender and S. Boettcher, Real spectra in non-hermitian hamiltonians having PT symmetry, Phys. Rev. Lett. 80, 5243 (1998).

Hassan2015 A. U. Hassan,  H. Hodaei, M. A. Miri,  M. Khajavikhan, and D. N. Christodoulides,  Nonlinear reversal of the PT-symmetric phase transition in a system of coupled  semiconductor microring resonators, Phys. Rev. A 92, 063807 (2015).
 
Ge2016 L. Ge and R. El-Ganainy, Nonlinear modal interactions in parity-time (PT)  symmetric lasers, Sci. Rep. 6, 24889 (2016).

Kepesidis2016 K.  V.  Kepesidis,  T.  J.  Milburn,  J.  Huber,  K.  G.  Makris,  S. Rotter, and P. Rabl, PT-symmetry breaking in the steady state of microscopic gain-loss systems, New J. Phys. 18, 095003 (2016).

Assawaworrarit2017  S. Assawaworrarit, X. Yu, and S. Fan, Robust wireless power transfer using a nonlinear parity-time-symmetric circuit, Nature 546, 387 (2017).

Ganainy2018 R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Non-Hermitian physics and PT symmetry, Nature Phys. 14, 11 (2018).




Kepesidis2013 K. V. Kepesidis, S. D. Bennett, S. Portolan, M. D. Lukin, and P. Rabl, Phonon Cooling and Lasing with Nitrogen-Vacancy Centers in Diamond, Phys. Rev. B 88, 064105 (2013).

hofer2017  P. P. Hofer, M. Perarnau-Llobet, L D. M Miranda, G. Haack,  R. Silva, J. B. Brask, and N. Brunner,  Markovian master equations for quantum thermal machines: local versus global approach,  New J. Phys. 19, 123037 (2017).
 
Plenio2018 M. T. Mitchison and M. Plenio, Non-additive dissipation in open quantum networks out of equilibrium, New J. Phys. 20, 033005 (2018).

 





Hatanaka2014 D. Hatanaka, I. Mahboob, K. Onomitsu, and H. Yamaguchi, Phonon waveguides for electromechanical circuits, Nat. Nanotechnol. 9, 520 (2014).

Huang2016 P. Huang, L. Zhang, J. Zhou, T. Tian, P. Yin, C. Duan, and J. Du, Nonreciprocal Radio Frequency Transduction in a Parametric Mechanical Artificial Lattice, Phys. Rev. Lett. 117, 017701 (2016).

Cha2018 J. Cha and C. Daraio, Electrical tuning of elastic wave propagation in nanomechanical lattices at
MHz frequencies, Nat. Nanotechnol. 13, 1016 (2018).

Patel2018 R. N. Patel, Z. Wang, W. Jiang, C. J. Sarabalis, J. T. Hill, and A. H. Safavi-Naeini, A single-mode phononic wire, Phys. Rev. Lett. 121, 040501 (2018).








Haeffner2014 M. Ramm, T. Pruttivarasin, and H. Häffner, Energy transport in trapped ion chains, New J. Phys. 16, 063062 (2014).





Underwood2012 D. L. Underwood, W. E. Shanks, J. Koch, and A. A. Houck,
Low-Disorder  Microwave  Cavity  Lattices  for  Quantum Simulation with Photons, Phys. Rev. A 86, 023837 (2012).

Houck2017 M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Observation of a Dissipative Phase Transition in a One-Dimensional Circuit QED Lattice, Phys. Rev. X 7, 011016 (2017).

Mirhosseini2018 M. Mirhosseini, E. Kim, V. S. Ferreira, M. Kalaee, A. Sipahigil, A. J. Keller, and O. Painter, Superconducting metamaterials for waveguide quantum electrodynamics,  Nat. Commun. 9, 3706 (2018).








Aspelmeyer2014 M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. 86, 1391 (2014).

Gu2017 X. Gu, A. Frisk Kockum, A. Miranowicz, Y.-X. Liu, and F. Nori, Microwave photonics with superconducting
quantum circuits, Phys. Rep. 718, 1 (2017).

Leibfried2003 D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Quantum dynamics of single trapped ions, Rev. Mod. Phys. 75, 281 (2003).


WallsMilburn D. F. Walls, and G. J. Milburn, Quantum Optics (Springer, 1994).

GardinerZoller C. W. Gardiner and P. Zoller, Quantum Noise (Springer, 2000).

Gardiner C. Gardiner, Stochastic methods  (Springer, Berlin, 2009).


Peng2014 B. Peng, S. K. Özdemir, S. Rotter, H. Yilmaz, M. Liertzer, F. Moni, C. M. Bender, F. Nori, and L.
Yang, Loss-induced suppression and revival of lasing, Science 346, 328 (2014).

Brandstetter2014 M. Brandstetter, M. Liertzer, C. Deutsch, P. Klang, J. Schöberl, H. E. Türeci, G. Strasser, K. Unterrainer, and S. Rotter, Reversing the pump dependence of a laser at an exceptional point, Nat. Commun. 5, 4034 (2014).




Diehl2010 S. Diehl, A. Tomadin, A. Micheli, R. Fazio, and P. Zoller, Dynamical Phase Transitions and Instabilities in Open Atomic Many-Body Systems,  Phys. Rev. Lett. 105, 015702 (2010).

Nagy2011 D. Nagy, G. Szirmai, and P. Domokos, Critical exponent of a quantum-noise-driven phase transition: The open-system Dicke model, Phys. Rev. A 84, 043637 (2011).

Oztop2012 B. Öztop, M. Bordyuh, Ö. E. Müstecaploglu, and H. E. Türeci, Excitations of optically driven atomic condensate in a cavity: theory of photodetection measurements, New J. Phys. 14, 085011 (2012).

Kessler2012 E. M. Kessler, G. Giedke, A. Imamoglu, S. F. Yelin, M. D. Lukin, and J. I. Cirac, Dissipative phase transition in a central spin system, Phys. Rev. A 86, 012116 (2012).

Casteels2016 W. Casteels, F. Storme, A. LeBoite, and C. Ciuti, Power Laws in the Dynamic Hysteresis of Quantum Nonlinear Photonic Resonators, Phys. Rev. A 93, 033824 (2016).

Hwang2018 M.-J. Hwang, P. Rabl, and M. B. Plenio, Dissipative Phase Transition in the Open Quantum Rabi Model, Phys. Rev. A 97, 013825 (2018). 

Zyczkowski K. Zyczkowski, P. Horodecki, A. Sanpera, and M. Lewenstein, Volume of the set of separable states, Phys. Rev. A 58, 883 (1998).

Vidal G. Vidal and R. F. Werner, Computable measure of entanglement, Phys. Rev. A 65, 032314 (2002).



 
Sudarhan E. C. G. Sudarshan, Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams, Phys. Rev. Lett. 10, 277 (1963).

normalordering J. M. Vargas-Martínez, H. Moya-Cessa, and M. Fernández Guasti,  Normal and anti-normal ordered expressions for annihilation and creation operators, Revista mexicana de fi­sica E 52, 13 (2006).

molmer J. Dalibard, Y. Castin, and K. Mølmer,  Wave-function approach to dissipative processes in quantum optics, Phys. Rev. Lett. 68, 580 (1992).

daley A. J. Daley, Quantum trajectories and open many-body quantum systems, Adv. Phys. 63, 77 (2014).

 































