Hidden time-reversal symmetry, quantum detailed balance and exact solutions of driven-dissipative quantum systems

Driven-dissipative quantum systems generically do not satisfy simple notions of detailed balance based on the time symmetry of correlation functions. We show that such systems can nonetheless exhibit a hidden time-reversal symmetry which most directly manifests itself in a doubled version of the original system prepared in an appropriate entangled thermofield double state. This hidden time-reversal symmetry has a direct operational utility: it provides a general method for finding exact solutions of non-trivial steady states. Special cases of this approach include the coherent quantum absorber and complex-$P$ function methods from quantum optics. We also show that hidden TRS has observable consequences even in single-system experiments, and can be broken by the non-trivial combination of nonlinearity, thermal fluctuations, and driving. To illustrate our ideas, we analyze concrete examples of driven qubits and nonlinear cavities. These systems exhibit hidden time-reversal symmetry but not conventional detailed balance.


I. INTRODUCTION
Time-reversal is a basic symmetry that plays a crucial role in a vast variety of physical systems. For open classical systems subject to dissipation and driving, it manifests itself as detailed balance constraints on transition rates (or equivalently drift and diffusion functions). It also places a strong symmetry constraint on steadystate two-time correlation functions A(t)B(0): they must be invariant when each quantity is replaced by its timereversed version and t → −t. This symmetry is sometimes referred to as Onsager symmetry, as it plays a crucial role in the derivation of Onsager reciprocity relations. In classical systems, this symmetry has a direct operational utility: it provides a simple route for finding steady state probability distributions (i.e. potential conditions that can be used to solve Fokker-Planck equations [1]).
There is a long history of works that extend notions of Onsager symmetry and detailed balance to quantum open systems described by a Markovian master equation in Lindblad form [2][3][4][5][6][7][8][9][10][11][12][13]. The most natural definition requires that steady-state correlation functions in the quantum theory obey an Onsager symmetry analogous to the classical case [2,3,8]; this condition necessarily holds if the microscopic system-bath dynamics obey time-reversal symmetry [3]. Later more formal works considered generalized definitions of quantum detailed balance [9,10], framed in terms of quantities that are not directly measurable and whose physical interpretation is somewhat opaque. The ultimate operational utility of all these quantum definitions of detailed balance are unclear. Unlike classical detailed balance, these quantum symmetries are not known to enable a simple method for finding a non-trivial system's steady state density matrix [14].
In this work, we introduce a powerful, symmetry-based formulation of quantum detailed balance (QDB) that goes beyond the simple definition in Ref. [2], and that directly enables an efficient way for finding non-trivial steady states. Our work builds on Ref. [11], which showed that a particular generalized definition of QDB introduced in Ref. [15] can be formulated using an entangled, thermofield double state [16]. We use this to introduce the notion of "hidden" time reversal symmetry (TRS) in an open quantum system. This anti-unitary symmetry need not reveal itself through some simple invariance of the original master equation, nor through a standard Onsager symmetry of two-time correlation functions. Instead, this symmetry is directly tied to a time symmetry of correlation functions of a doubled version of the original system prepared in an thermofield double state whose form is directly tied to the symmetry operator (see Fig. 1(c)). Crucially, we show that a system can possess hidden TRS even if it fails to have the conventional quantum detailed balance (CQDB) defined in Ref. [2] (though in the limit of infinitely weak dissipation, these notions coincide). Hidden TRS is not just a formal curiosity: it provides a powerful tool for understanding complex non-thermal and non-classical steady states. We show that the existence of hidden TRS directly yields a simple and direct method for analytically finding the steady state density matrix of a Lindblad driven-dissipative quantum system. This method is not limited to situations of weak driving, interactions or dissipation. It represents a generalization of the coherent quantum absorber (CQA) method introduced in Ref. [17], and extended in Ref. [18]. Hidden-TRS is also connected to well-known exact solution methods from quantum optics based on the complex-P phase space quasiprobability [19][20][21][22]: these methods can be viewed as special cases of our more general approach.
While experiments on doubled quantum systems prepared in thermofield double states have recently been performed [23], hidden TRS also has experimental consequences in experiments on just a single system. Unlike CQDB, systems with hidden TRS will not exhibit Onsager time-symmetry of all correlators. However, we show that there are always a class of special correlation functions that are guaranteed to have this time symmetry. This provides a direct means for probing hidden TRS (and its possible breaking) in a variety of experimental platforms. We explore in detail two classes of ubiquitous, experimentally-accessible systems (see Table  I): Rabi-driven qubits subject to dissipation, and drivendissipative nonlinear quantum cavities. These systems exhibit in general no correlation function time-symmetry, and hence do not possess CQDB as defined in [2]. They however do possess hidden-TRS in the low temperature or small nonlinearity limit. This explains the surprising exact solvability of a variety of driven nonlinear cavity models [18][19][20][21][22]. We explore how hidden TRS is broken in these models by the combination of non-zero temperature, driving and nonlinearity. For nonlinear cavities, breaking of detailed balance was extensively studied in the semiclassical limit [24][25][26][27][28][29][30].
The rest of this paper is organized as follows: in Sec. II we review the doubled-system formulation of classical detailed balance and the definition of CQDB [2,3]; we show that CQDB only holds for a limited class of systems that possess trivial steady states. Sec. III introduces the no- tion of hidden TRS, and connects it to the definition of generalized QDB introduced in Ref. [15]. In Secs. IV and V we demonstrate how the existence of hidden TRS enables an extremely direct method for finding exact solutions for steady states. In Sec. VI we discuss how a variety of driven quantum cavity models possess hidden TRS, while in Sec. VII we discuss how thermal fluctuations in some cases can break this symmetry. Sec. VIII shows how hidden-TRS underlies the complex-P function exact solution method. We conclude in Sec. IX.

A. Classical detailed balance
Consider a classical stochastic system with a discrete set of microstates indexed by integers n, whose Markovian dynamics is fully described by a set of transition rates Γ n→m . The time-dependent probability p(n, t) for the system to be in a given state n then obeys: We assume that this equation admits a timeindependent steady-state probability distributionp(n). This steady state is said to have detailed balance if there is a balancing of probability fluxes between any given pair of states and their time-reversed partners. Letting n denote the time-reversed version of the microstate n, the condition is (see e.g. [1]): This definition generalizes directly to systems with a continuous state space: Eq. (1) then becomes a Fokker-Planck equation, and Eq. (2) becomes a constraint on (b) (a) Figure 2. (a) System with discrete microstates, described by a classical master equation. Detailed balance (Eq. (2)) is equivalent to a time-symmetry of stationary correlation functions. (b) Equivalent formulation of classical detailed balance, involving a doubled system prepared in an initial correlated state given in Eq. (6); the auxiliary system B has no dynamics. The detailed balance condition in (a) is equivalent to requiring a time-symmetry of doubled-system correlators.
drift and diffusion matrices (so-called potential conditions). As is well known, these "zero probability flux" conditions give a direct way to find the steady state distribution. In the discrete case, one uses the zero-flux condition to iteratively find the probability of each microstate in terms of the rates. The continuous version of this yields the steady state in terms of a potential function that is determined by drift and diffusion matrices (see e.g. [1]).
One can equivalently define detailed balance by a time symmetry of correlation functions, what we term here an Onsager symmetry. Suppose X and Y are arbitrary functions of the microstate n of our system. Steady state correlation functions can be defined in the usual manner in terms of conditional probabilities p(m, t|n, 0), which can be computed from Eq. (1). For example: The detailed balance condition of Eq. (2) is then equivalent to requiring that the following symmetry hold for all steady state correlation functions: Here, the time-reversed functionX is defined asX(n) ≡ X(ñ).

B. Doubled-system formulation of classical detailed balance
Consider now an alternate but equivalent formulation of classical detailed balance [11]. We imagine making a copy of our original system that has exactly the same set of microstates as the original. This auxiliary system (system B) is completely static, whereas the original system (system A) retains its transition rates and dynamics, see Fig. 2. To be explicit, the doubled-system is described by a probability distribution p AB (n A , n B ; t) which satisfies the master equation: We further assume that the doubled system is initially prepared in a correlated state described by the probability distribution: wherep(n) is the single-system steady state of interest. Note that the correlations in this state partner each microstate n in system A with its time-reversed partnerñ in system B. It is easy to confirm that the marginal probability distribution describing each subsystem alone (A or B) is time independent and equal top(n) for system A, andp(ñ) for system B. In contrast, the total state of the two systems will evolve non-trivially. The result is that while quantities involving only each subsystem are completely static, correlations between the two subsystems will evolve in time.
Given this dynamics and initial state, we now ask whether inter-system correlations at time t > 0 are invariant if we "exchange" the two subsystems. To be concrete, we consider two observable quantities X and Y that we could measure on either system A or B; these are functions of microstates, i.e. X A (n, m) = X(n), X B (n, m) = X(m). We now define a doubled-system Onsager time-symmetry by the requirement: for all possible observables X, Y . We stress that system B is not dynamical. Superficially, this looks very different from the standard Onsager time-symmetry constraint in Eq. (4). We are not exchanging observation times, but rather the subsystem in which each quantity is measured. There is also no explicit time-reversal operation in this equation (it is instead hardwired into the initial correlated state). Despite these differences, one can easily show (see Ref. [11] and App. A) that the above symmetry relation is completely equivalent to standard Onsager time symmetry. While this formulation thus provides nothing new in the classical context, we will see that it motivates an extremely powerful generalized notion for quantum systems, namely the notion of a hidden time-reversal symmetry.

C. Markovian quantum open system: general setting
Onsager-like correlation function time-symmetries can also be considered in the context of open quantum sys-tems [2,3,8]. Our general setting will be an open system whose dynamics is described by a Markovian Lindblad master equation. The dynamics is completely specified by the Hermitian system HamiltonianĤ and a set of jump operatorsĉ l that describe the influence of external dissipative reservoirs. DefiningĤ eff =Ĥ − i 2 lĉ † lĉ l , our general Lindblad equation of motion for the system density matrixρ(t) is: where we have introduced the Liouvillian superoperator L. Given a particular steady state of this equation (i.e. a time-independent solutionρ ss ), steady-state correlators between two system operatorsX andŶ obey time-translation symmetry and can be calculated from the master equation (see, e.g., [31]). To state this compactly, we first introduce the adjoint Liouvillian super-operatorL, determined by (see e.g. Ref. [32]) Correlation functions are then computed as As mentioned in the introduction, a variety of definitions of quantum detailed balance have been formulated for Lindblad Markovian master equations [2,5,11,33]. The simplest and best-known definition of quantum detailed balance involves the same constraint on correlation functions that exists in the classical setting [2,3]. With this definition (which we call conventional quantum detailed balance (CQDB)), detailed balance requires: HereX andŶ are arbitrary system operators, and tilde is used to denote the time-reversed version of an operator, i.e.B ≡TBT −1 , whereT is the anti-unitary timereversal operator. Note that on the RHS of Eq. (11), we Besides being simple to state, the CQDB correlation function symmetry also has a direct connection to a microscopic symmetry: this correlation function time symmetry necessarily holds if the microscopic system-bath dynamics has time-reversal symmetry [3]. This requires both that the entire system-plus-bath Hamiltonian has time-reversal, and that the bath state also relaxes to a time-independent state (not just the system state). This connection to microscopics (and the fact that it involves measurable correlation functions) has led many to consider CQDB the most natural generalization of classical detailed balance to quantum systems [8]. Note that CQDB is also referred to as "GNS detailed balance" [12].

E. CQDB implies a trivial steady state
As has been noted previously [4,9,10], the CQDB condition is extremely restrictive: it only holds for systems with steady states that are diagonal in the energy eigenstate basis (i.e. [Ĥ,ρ ss ] = 0) [34]. Such states are in some sense trivial, as they can be found by solving a classical (Pauli) master equation, obtained by setting all energy-eigenstate coherences in Eq. (8) to zero. The resulting steady state can thus be interpreted classically: the Hamiltonian plays no role, and steady-state probabilities are determined by balancing incoherent transition rates between different eigenstates. A simple proof of this result is presented in App. B.
The upshot is that CQDB is found in an extremely limited class of systems, and is not a useful tool for finding non-trivial quantum steady states; in fact, the presence of CQDB makes it impossible to have such a state. In App. C, we show explicitly that an extremely simple model of a Rabi-driven dissipative qubit fails to have CQDB (though it will possess the hidden time-reversal symmetry we introduce in the next section).

A. Basic formulation
As we have discussed, the simple CQDB condition of Eq. (11) corresponds directly to microscopic timereversal symmetry [3]. We now consider something more general, the notion of hidden time-reversal symmetry, which we will formulate in a doubled version of our original system (in rough analogy to the classical construction in Sec. II B). This unusual symmetry will connect directly to an abstract variant of quantum detailed balance (so called "SQDB-θ") studied in the mathematical physics literature [9,10], and which was recently linked to entanglement [11]. Our work complements these studies by providing a direct physical motivation for this definition and connects it explicitly to symmetry. More importantly, we show that this formulation has great practical utility: it allows us to solve for non-trivial steady states. This connection was not previously known.
Our starting point is again the Lindblad master equation of Eq. (8) and particular steady stateρ ss , which we write in diagonal form aŝ ρ ss = n p n |n n|.
Throughout this work, we will assume thatρ ss is full rank, and further, that the p n have no degeneracies. We next construct a purification of this state: an entangled pure state |ψ T of a doubled version of our system that yieldsρ ss if we trace out the auxiliary system B. We take system B to have the same Hilbert space as the original system A. To construct |ψ T , we first chose an anti-unitary operatorT which will define our hidden time-reversal symmetry. We then use this choice to construct |ψ T in a manner that mimics the classical doubled-system state in Eq. (6): we pair each pointer state |n in the original system with its time-reversed partner in the auxiliary system. We thus have With this definition, |ψ T is invariant under a gauge change of theρ ss eigenkets |n → e iαn |n . The state |ψ T has the form of a so-called "thermofield double state"; such states have been studied in many different contexts [16], though usually without including a time-reversal operation in its definition. We stress that the form of |ψ T is contingent on the choice ofT . We next specify the dynamics of the doubled system: as we did in the classical case, we take subsystem A to evolve as per the master equation in Eq. (8), but take the auxiliary subsystem B to have no dynamics at all. The dynamics of the joint system thus follows the Lindblad master equation: where L is our original single-system Liouvillian, and 1 B denotes the unit superoperator acting on subsystem B. Starting the doubled system at t = 0 in the pure statê ρ AB (0) = |ψ T ψ T |, this dynamics leaves the reduced density matrices of each subsystem invariant. It does not however leave the full joint stateρ AB time invariant. The result is that single-system expectation values are timeindependent, but correlations between them can evolve.
We will now use these evolving inter-system correlations to define the notion of a hidden time reversal symmetry. We define single-subsystem operators in the natural way, i.e.X A ≡X ⊗1,X B ≡1 ⊗X. Starting in the state |ψ T and evolving as per Eq. (14), hidden TRS holds if all inter-system correlations obey the following symmetry: Stated explicitly, the symmetry requires that intersystem correlations at any time t ≥ 0 are invariant if we exchange which subsystem each quantity is measured in. Note that as system B has no dynamics, we can equivalently write this condition in the Heisenberg pic-ture. Defining the condition for having hidden TRS condition then becomes Note that for a general system that does not have hidden TRS, the definition in Eq. (16) implies that not only can C TFD XY (t) be asymmetric, it can even fail to be continuous at t = 0. We also stress again that the hidden-TRS condition above is contingent on the choice of anti-unitaryT . As we will see, there exist physical systems where hidden-TRS is in some sense degenerate: there are a whole family of distinct operatorsT for which Eq. (17) holds.
Despite mirroring the classical doubled-system construction, in the quantum case Eq. (17) gives us something truly new: there are systems that fail to satisfy the Onsager symmetry condition of Eq. (11), but nonetheless satisfy the generalized condition in Eq. (17). In these cases, we describe the particular anti-unitary operatorT used to define |ψ T as a hidden time-reversal symmetry. Note that if a system has hidden TRS, the steady state is invariant under the corresponding time-reversal operation:Tρ ssT −1 =ρ ss (18) This follows directly by assuming hidden TRS (i.e. Eq. (15)), and taking the choiceŶ = 1. The only way the remaining condition can hold for allX is if Eq. (18) holds. To gain intuition for the role of entanglement in our formulation of hidden TRS, it is useful to express the doubled correlation function in terms of the eigenstates ofρ ss , c.f. Eq. (12). Using Eq. (16), we have where For times t < 0, Eqs. (19)- (21) are defined analogously to C TFD XY (t) itself. The above separation has a direct physical significance, as it is the second contribution in Eq. (21) that encodes the contribution from quantum entanglement. To see this explicitly, suppose that we initially had prepared our doubled system in a state that only had classical correlations (but not entanglement) between systems A and B, i.e. ρ AB (0) = n p n |n,ñ n,ñ|. (22) In this case, our entire inter-system correlation functions would be given by Eq. (20): We can thus view C en XY (t) in Eq. (21) as the extra contribution to the TFD correlator that stems from having quantum entanglement (as opposed to purely classical correlations).

B. CQDB as a special case of hidden TRS
At this point, it is natural to ask whether there is any simple relation between our notion of hidden TRS and conventional CQDB (which we stress is directly connected to microscopic time-reversal of the complete system-plus-environment [3]). The first key result here is that CQDB is a special case of hidden-TRS: any system satisfying CQDB automatically has a hidden-TRS, but the converse is not true. Recall from Sec. II E that systems satisfying CQDB necessarily have somewhat trivial steady states (in that they are diagonal in the energy eigenstate basis). In contrast, there are many systems with have hidden TRS but not CQDB, and have steady states with non-zero coherences between energy eigenstates. This is a crucial result: hidden TRS does not preclude having a non-trivial steady state.
To understand the origin of the above result, it is useful to introduce an effective dimensionless Hermitian Hamiltonian defined by the steady-state density matrixρ ss , the so-called modular Hamiltonian: Any system possessing CQDB has an extremely tight constraint on its dynamics [35]: dynamical evolution generated by the full Liouvillian L must commute with evolution generated byĤ ρ (i.e. by the unitary e −iĤρt ). This symmetry allows one to demonstrate that CQDB is a special case of hidden TRS, using the same methods introduced to study a class of generalized QDB conditions in Refs. [9,10]. In Appendix E, we give a non-technical proof of this result: for any system with modular symmetry, CQDB and hidden TRS are equivalent [9], and thus CQDB implies hidden TRS. The second key result connecting hidden TRS and CQDB involves the limit of vanishing incoherent dynamics. Consider a system that possesses hidden-TRS even as the strength of the incoherent dynamics is tuned to zero (i.e. replace the jump operatorsĉ l → λĉ l in the master equation, and let λ → 0). In such systems, hidden TRS reduces to CQDB in the limit of vanishing dissipation. Hence, while in many cases finite dissipation destroys CQDB, hidden TRS can continue to be a symmetry.
This result can also be understood using the modular Hamiltonian. In the limit of vanishing dissipation, the full dynamics is completely generated by the Hermitian HamiltonianĤ, i.e.
Further, in this limit the steady stateρ ss must commute withĤ in order to be stationary. These two facts then necessarily imply that in this asymptotic zero dissipation limit, the full Liouvillian commutes with the modular Hamiltonian. As discussed above, this then implies that our system has CQDB in the zero-dissipation limit. This proves the desired result: systems with hidden TRS always recover CQDB in the weak dissipation limit. At a more physical level, systems with hidden TRS do not necessarily have (single system) correlation functions that are time-symmetric. However, the lack of time-symmetry vanishes in the zero-dissipation limit.
C. Hidden TRS has observable consequences for a single system Hidden TRS might initially seem to be experimentally irrelevant in most settings, as it is defined in terms of TFD correlators (c.f. Eq. (16)) that require one to prepare two copies of the system of interest in an initial entangled state. This is usually extremely challenging (though see the recent trapped ion experiment in Ref. [23]). Formally, one could define this correlator as a single-system quantity involving a state-dependent observable. We first introduce a superoperator J which acts on single-system operators: J is well-defined and unique whenρ ss is full-rank, and is given by the explicit formula (see App. D) The TFD correlator in Eq. (16) can then be written as a single system correlator, e.g. for t > 0: Using this correspondence, the defining symmetry condition of hidden TRS in Eq. (17) can be written in a manner that only involves a single system: The above condition is formally equivalent to one of the many generalized quantum detailed balance conditions discussed in Ref. [9,10] (so-called SQDB-θ). It is also referred to as "KMS detailed balance" [12]. In general, the symmetry condition in Eq. (29) is not helpful as an experimental tool, as it involves an operator whose very form depends on the system stateρ ss , i.e. it is nonlinear inρ ss . This is reminiscent of the problem of measuring Renyi entropies [36][37][38]. Remarkably, all hope is not lost. As we will show in Sec. V, systems with hidden TRS are guaranteed to have their effective HamiltonianĤ eff and jump operatorsĉ k transform very simply under the action of J . For correlation functions involving these operators, Eq. (29) becomes a standard (and often simple) single-system correlator.
We thus have a key result that makes it possible to directly and simply test for hidden TRS in experiment: hidden TRS implies that only a certain restricted class of correlation functions (directly identifiable from the master equation) are guaranteed to have Onsager timesymmetry.

D. Example: Hidden TRS in dissipative
Rabi-driven qubit We show in App. C that a simple Rabi-driven qubit with loss fails to respect CQDB: its correlation functions do not exhibit Onsager time-symmetry, regardless of how one tries to define time-reversal. Here we show that this system nonetheless possesses a hidden TRS. Working in the rotating frame set by the drive, and making a rotating wave approximation, the master equation has the form of Eq. (8) with M = 1 and For simplicity, we consider the resonant-driving case ∆ = 0 in what follows. This system has a unique steady state, and corresponding thermofield doubled states can be constructed according to Eq. (13). This construction depends on the choice of the hidden TRS operatorT ; as discussed in App. C, the requirement that the steady state be TRS invariant (c.f. Eq. (18)) constrainsT to the form: where b = Ω/κ,K z is the complex conjugation operator acting in theσ z basis, and ψ is at this stage an arbitrary real parameter.
To determine whether our system has hidden-TRS, we must find aT such that Eq. (17) holds (i.e. intra-system TFD correlators have a time symmetry). We thus compute TFD correlators between different pairs of Pauli operators. Here we will consider only the correlation function C TFD yz (t), and we leave the other two correlation functions to App. G.
Using Eq. (19) we can decompose this into the classical correlation C cl yz (t) and the entanglement correction C en yz (t) as: C TFD yz (t) = C cl yz (t) + C en yz (t). The classical correlation is independent of ψ, and its time asymmetry is nonzero irrespective of howT is chosen: Here, α = √ 16b 2 − 1. We thus see that were we to neglect the entanglement correction, the system could never have hidden TRS. Now, we look at the time asymmetry of the entanglement correction, which is dependent on ψ: Comparing Eqs. (32) and (33), we see that for the TRS with ψ = π, the entanglement correction to C TFD yz (t) modifies the classical correlation in just the right way to cancel the net time asymmetry. We see just how stark the effect is in Fig. 3 which compares the full TFD correlation function with the classical correlation terms for b = 1 at the TRS ψ = π. For reference, the single qubit correlation function C yz (t) for b = 1 at ψ = π is also included. This result shows the importance of the entanglement correction to restoring a notion of detailed balance to the Rabi-driven qubit and highlights the fact that the notion of hidden TRS has a distinctly quantum nature.
From the above, we conclude that our model does have a unique hidden TRS, described by the anti-unitary op-eratorT hT In App. G we confirm that the remaining two correlation functions are time symmetric for this TRS.
It is interesting to consider the form ofT h in various limits. For weak Rabi driving (i.e. b → 0), Eq. (34) reduces toT h =K z . In this limit the qubit system in fact satisfies SDQB withT =K z (i.e. all correlation functions have standard Onsager time-symmetry). In the strong drive limit b 1,T h → −iσ xKz . Up to a phase, this just complex conjugation in theσ y basis. To make sense of this, consider the steady state Eq. (C2) in this limit. To first order in small b −1 1, the steady state reduces toρ ss =1 2 + b −1 2σ y . The form of the hidden-TRS operator in this limit thus directly reflects the eigenvectors of ρ ss . Furthermore, one can show that for any b, the hidden TRST h corresponds to complex conjugation in the steady state eigenbasis. Figure 3. Correlation functions and hidden-TRS in a driven qubit. Stationary, connected σy(t)σz(0) correlation functions for the dissipative Rabi-driven qubit system in Eq. (30), for a drive Ω equal to the decay rate κ. Blue: the standard single-system correlation function Cyz(t) is asymmetric as a function of time, reflecting the fact that this system does not satisfy conventional quantum detailed balance. Red: Two-qubit correlator for a system prepared in a TFD state corresponding to the hidden-TRS operatorT defined in Eq. (34). All TFD correlators symmetric in time, reflecting the presence of hidden TRS. Green: "classical" part of the TFD correlator (c.f. Eq. (19)), which has no time symmetry. The lack of symmetry shows that the importance of entanglement in the definition of hidden-TRS.

IV. HIDDEN TIME REVERSAL SYMMETRY AND DYNAMICAL CONSTRAINTS
We have now introduced our notion of hidden TRS (c.f. Eqs. (15), (17)), and demonstrated that this symmetry can hold even when the more standard CQDB symmetry is broken. It still however may seem that hidden-TRS is nothing more than a formal curiosity. We show here that this is not the case: hidden-TRS is a symmetry that has direct operational utility in helping us understand complex phenomena, as it enables the exact solution of steady-states of non-trivial systems. In particular, it is the symmetry condition that enables the surprising but powerful coherent quantum absorber method introduced in Ref. [17] and extended in Ref. [18].

A. Equivalent subsystem dynamics and hidden TRS as a self-dual condition
We start by demonstrating that the hidden TRS condition can also expressed as a kind of dynamical equivalence between the two subsystems in our TFD state. Consider a general system and a TFD state which does not necessarily satisfy the hidden TRS condition of Eq. (17). We stress that the TFD state is defined completely by the steady state of interestρ ss and choice of anti-unitaryT . We will takeρ ss to be full rank in what follows, and will consider intra-system correlations in this TFD state as defining a bilinear form: whereX,Ŷ denote arbitrary single-system operators. This bilinear form can then be used to define the dual E * of any given single-system superoperator E via Of particular interest is the case where E is the adjoint evolution operator E t = exp[Lt] defined in Eq. (10) in terms of the adjoint LiouvillianL (c.f. Eq. (9)). The LHS of Eq. (36) then describes the correlation of a subsystem-A operator at time t and a subsystem B operator at time zero. In this case, the dual (E t ) * has a direct physical interpretation: it represents an alternate and equivalent time evolution of subsystem B that would result in the same inter-system correlation. This dual time evolution can be written as (E t ) * = exp L * t . Thus, for a given subsystem-A dynamicsL, we have a corresponding "mirrored" dynamicsL * for subsystem-B, defined by the constraint that it yield identical inter-system correlations, i.e.
These notions now give an extremely transparent way to rephrase the hidden-time reversal condition of Eq. (17): the original system-A dynamics and its mirrored version must be identical, that isL is self-dual, To see this, note first that if a system satisfies the hidden TRS condition of Eq. (17), then the bilnear form in Eq. (35) must be symmetric, i.e. X ,Ŷ T = Ŷ ,X T ; this follows from the t = 0 limit of Eq. (17). This in turn implies that the steady stateρ ss of the original master equation must be invariant under our hidden timereversal operatorT (i.e. consider the caseŶ =1). Combining these two conditions lets us express the hidden-TRS condition of Eq. (17) as: where for either subsystem, . This now looks more like a standard Onsager-type correlation function symmetry, except that the two operators are measured on different subsystems. Finally, comparing this equation against Eq. (36) directly yields the self-duality condition in Eq. (38) [39].

B. Hidden TRS as a symmetry of the Liouvillian
We now show that the hidden TRS condition can be viewed as a dynamical symmetry that directly constrains the system's adjoint LiouvillianL. To do this, we step back and consider a general system andT , such that hidden TRS is not necessarily satisfied. We then explicitly construct the dual LiouvillianL * that generates the mirrored-system dynamics, by considering each term in Eq. (9). Our construction will explicitly make use of the exchange superoperator J introduced in Eqs. (26) and (27); recall that this superoperator lets us convert subsystem-A into corresponding subsystem-B operators (and vice-versa) such that TFD expectation values are preserved.
The exchange superoperator allows us to efficiently express the desired dual of the adjoint LiouvillianL. This is done using the following relations, that follow directly from the definition of J in Eq. (26): We can thus obtain an explicit expression for the desired dualL * as RecallĤ eff is the effective non-Hermitian Hamiltonian, andĉ l the jump operators in our original Lindblad master equation Eq. (8). We thus see that the properties of the system-B "mirrored dynamics" are encoded in the exchange superoperator J .
We now ask what constraints ensue when we insist that the hidden TRS condition holds, and henceL * =L. For two Liouvillians (each defined with traceless jump operators) to be equivalent, the effective Hamiltonians for each must be identical (up to an additive real constant), and the jump operators related by a unitary mixing matrix (see e.g. Ref. [40]).
Hence, insisting that our system has hidden TRS leads to the following constraint equations: E is a real number, and U lk is a M × M unitary matrix. The last constraint on U (i.e. that it is involutory) follows from the fact that if hidden TRS holds, then the steady state is itself invariant underT (c.f. Eq. (18)). This resulting additional symmetry of the TFD state then implies (via Eq. (26)) that two exchanges yield the identify superoperator: J 2 = 1. This immediately constrains the unitary matrix U lk to have purely real eigenvalues [41].
Eqs. (43) represent necessary and sufficient conditions for our system to have a hidden TRS. They are however somewhat unwieldy, as they directly involve the exchange superoperator, which is itself a function ofT and the steady stateρ ss . We can eliminate the explicit appearance of J by using the fact that sinceρ ss is full rank, the TFD state is separating (see e.g. [42]): if two subsystem-A operators have the same action on the TFD state, then they must be identical operators. Stated explicitly: Using this, we can eliminate J from each equation in Eqs. (43) by taking each side of each equation to be a system-A operator, and applying it to the TFD state |ψ T . Using the definition of the exchange operator on the resulting state, this gives us an equivalent but more useful set of constraint equations: Eqs. (45)- (47) are the main result of this subsection: they express the existence of a hidden TRS symmetry directly as a constraint on the Hamiltonian and jump operators that define our open system dynamics. Heuristically, these conditions imply that the action ofĤ eff and the jump operators are "almost" the same whether they act on subsystem A or B. Hidden TRS requires that these equations must hold for some pure state |ψ T of the doubled system, some constant E and some involutory M × M unitary matrix U . One can view this as a generalization of the classical detailed balance condition in Eq. (2). While the classical condition only involves transition rates, our quantum conditions above constrain both the incoherent dynamics generated by theĉ l operators, as well as the coherent system HamiltonianĤ.
We stress that the above equations are equivalent to those derived in Ref. [9] when considering the abstract "SQDB-θ" version of quantum detailed balance. By phrasing these conditions directly in terms of the thermofield double state, we will be able to directly exploit them as a means for efficiently finding unknown steady states (something that was not considered previously).
Finally, we note that Eq. (45) -Eq. (46) tells us that if hidden-TRS holds, than the action of the exchange superoperator is extremely simple when acting onĤ eff or the jump operatorsĉ l . This means that for doubled-system TFD correlation functions involving these operators can be directly converted to single-system correlation functions. As discussed extensively in Sec. III C, this gives a direct means for experimentally testing for hidden-TRS in a single system: hidden-TRS ensures that a certain reduced class of correlation functions will obey Onsagerlike time symmetry (c.f. Eq. (29)).

A. Basic idea
Classical detailed balance has a profound operational utility: it provides an extremely efficient method for finding the steady state of a given dynamical model (i.e. so called potential solutions of Fokker-Planck equations [1]). It is thus natural to ask whether something similar is possible using our notion of hidden TRS. If a system satisfies this symmetry, does this directly provide a method for solving for the steady state? As we now show, the answer is a resounding yes. The existence of hidden TRS places a strong constraint on the form of our dynamics via Eqs. (45) - (47). These equations also provide an efficient method for finding an unknown steady state. To see this, we change perspective, and view |ψ T in these equations as an unknown pure state of a doubled version of the original system. The goal is then to find a pure state |ψ T , constant E, and unitary matrix U such that Eqs. (45) -(47) are satisfied. If we are able to do this, then as we will show, our system has hidden TRS, and the desired system-A steady stateρ ss is obtained by tracing out system-B from |ψ T . Conversely, if we cannot do this, then our system does not have hidden TRS, and there is no generic simple route to finding the steady state.
We stress that the above procedure for finding the steady state is simpler than a direct brute-force approach. Suppose our original system has a Hilbert space dimension d. Without assuming hidden TRS, solving for the steady state of Eq. (8) reduces to the problem of solving for the null space of a matrix with dimensions d 2 × d 2 . Without additional assumptions, this matrix does not have any obvious sparseness properties. In contrast, with the assumption of hidden TRS, we need to solve Eqs. (45) and (46). Each of these M + 1 equations also involves a d 2 ×d 2 matrix. However, each of these matrices has a simplified structure as there are no terms corresponding to an interaction between the two subsystems. As a result, there can be at most O[d 3 ] non-zero matrix elements. In addition, our constraint equations decouple the effective Hamiltonian physics (Eq. (45)) from the incoherent "jump" physics (Eq. (46)). This effective non-interacting property provides a considerable simplification, as we will exploit more fully in the next section.

B. Connection to perfect quantum absorbers
As we show below, the presence of hidden TRS guarantees that we can construct a simple mirrored system that perfectly absorbs everything emitted by the main system into its environment. Such absorbing systems have been studied previously as a method for deriving exact solutions of certain Lindblad master equations [17,18]. Our discussion here will provide a generalization of this "coherent quantum absorber" method to systems with multiple jump operators, and also show that its success is indeed intimately connected to hidden time-reversal symmetry.
To establish this connection, we again consider a general system described by the master equation Eq. (8) with a steady stateρ ss . We also construct a doubled system as in Sec. III A with a TFD state given by Eq. (13). To start, we do not assume that the system has hidden TRS. As discussed in Sec. IV, for a given subsystem A dynamics (generated byL), we can always construct a corresponding "mirrored" dynamics on subsystem B (generated bȳ L * ), such that either evolution generates the same timedependent inter-system correlations, c.f. Eq. (37).
Somewhat remarkably, this mirrored dynamics is also exactly what is needed make subsystem-B a "perfect absorber" of energy and information emitted by subsystem-A into its environment (Fig. 4). This can be established by using the exchange superoperator J introduced in Eq. (26), which converts the action of a subsystem-A operator acting on the TFD state to a subsystem-B operator (and vice-versa). From the definition of J we have: whereĤ eff is the effective Hamiltonian in our master equation, andĉ l are the jump operators. As shown in App. H, these equations can be re-written as:Ĥ Here, the Hermitian HamiltonianĤ AB describes an interaction between the two subsystems in our doubled system:Ĥ withĤ We denote the Hermitian part of an operatorÂ as Re[Â]. Note thatĤ is nothing but the Hermitian Hamiltonian in our original master equation. Eq. (50) has an extremely suggestive form: it tells us that |ψ T is necessarily a zero-energy eigenstate of a Hermitian Hamiltonian describing a doubled system with an inter-system coupling, and that it is also annihilated by particular combinations of jump operators. Together, these conditions imply that |ψ T is a zero energy pure-state, steady-state of the cascaded doubled system sketched in Fig. 4c. In this cascaded system [43,44], there is an independent chiral (directional) waveguide associ-  ated with each jump operatorĉ l in our original master equation. These channels mediate a directional coupling between systems A and B, with B downstream from A.
Using the standard theory of cascaded quantum systems [32,43], the full master equation for this system is: Hereρ AB is the density matrix of the doubled system, and D[ẑ]ρ =ẑρẑ † − {ẑ †ẑ ,ρ}/2 is the standard Lindblad dissipation operator. One can easily verify that if |ψ T satisfies Eq. (50), then it is a steady state of Eq. (53). We thus have established the desired connection: the same formal construction that gives us a correlationconserving mirrored dynamics on subsystem-B also tells us the precise dynamics that is needed for subsystem-B to be a perfect absorber for subsystem-A. We stress that for each possible choice of candidate time-reversal operatorT , we have a different TFD state, a different mirrored-dynamics (i.e. HamiltonianĤ , jump operators J [ĉ l ]), and hence a different possible coherent quantum absorber.

C. Hidden TRS and simple absorbing dynamics
The cascaded master equation in Eq. (53) in principle provides a route for finding the steady state of the physical system A. If one could find the steady state of this master equation, then tracing out system B necessarily yields a steady state of the original single-system master equation. One could simplify this procedure by trying to find a pure state solution to Eq. (53). Of course, there is an obvious problem to this approach: the construction of Eq. (53) is contingent on already knowing the steady stateρ ss , as this is needed to construct the exchange superoperator J .
Things simplify considerably though in the case where our system possesses a hidden TRS. In this case, we can use Eqs. (45)- (47) to dramatically simplify the cascaded master equation for our system. The system-B jump operators and Hamiltonian are then given by for some involutory unitary M × M matrix U and real constant E, and the Hamiltonian of the coupled system becomeŝ where E is now implicitly absorbed into an energy-shift of the dark state.
We can now view this as a method for finding an unknown steady state of our original system-A master equation in Eq. (8). If we assume the existence of hidden TRS, finding this steady state is equivalent to finding a involutory M × M unitary matrix U and energy E, such that the cascaded master equation in Eq. (53) (with the simplifications of Eqs. (54) and (55)) yields a pure-state, steady state. This pure state then gives us the desired system-A steady state by just tracing over system B.
The technique detailed above is a generalized version of the CQA exact solution method introduced in Ref. [17] for solving master equations with a single jump operator. Our extension to systems with multiple jump operators involves a new object, the involutory unitary matrix U . We have shown that this solution technique is thus intimately connected to the notion of a hidden TRS, and thus to the generalized SQDB-θ quantum detailed balance conditions introduced earlier on purely formal grounds [9,10]. As far as we know, this is the first example of this notion of quantum detailed balance having an operational utility.

VI. HIDDEN-TRS IN NONLINEAR DRIVEN-DISSIPATIVE QUANTUM CAVITIES
At this stage, we have established the basic notion of hidden TRS. This symmetry can hold even if the more conventional CQDB condition (Sec. II D) is broken. Moreover, it directly enables a simple but powerful method for finding non-trivial steady states (Sec. V). We have illustrated these ideas by explicitly considering hidden TRS in a model of a dissipative Rabi-driven qubit (see App. C and III D). In this section, we turn to a more complex class of models. These describe a bosonic mode (canonical annihilation operatorâ) with a Kerr or Hubbard-U type interaction, subject to both one and two particle coherent driving, as well as one and two-particle losses. The system is described by a Lindblad master equation Eq. (8) with a coherent Hamiltonian: and with jump operatorŝ This model describes a dissipative cavity mode driven both with linear and parametric drives that have commensurate frequencies (working within rotating-wave approximations, and in a rotating frame that eliminates the time dependence of the drives). It is a ubiquitous system, having both been studied extensively in quantum optics, and more recently in the field of superconducting quantum circuits as a route to error-corrected quantum memories [45][46][47][48][49].
It is well known that the steady state of this class of systems can be found analytically using quantum-optical phase space methods [19,21,50]; more recent work has shown that these exact solutions can be derived more directly (and even extended) using the coherent quantum absorber method (CQA) [17,18]. An underlying explanation however for why these models are solvable has been lacking. We now have such an explanation: this class of models possess hidden-TRS, which directly leads to their solvability. In what follows, we discuss the nature of hidden-TRS in these systems, showing that hidden-TRS is present even though (generically) CQDB does not hold. Crucially, we show that the presence of hidden-TRS has observable consequences even in experiments on a single system: while generic correlation functions do not exhibit a time-symmetry, there is a special class of correlators that do. In Sec. VIII, we will show how hidden-TRS directly enables the required symmetry exploited in the complex-P function phase space methods that were first used to solve these systems [19,50].

A. Multiple non-trivial hidden-TRS symmetries
We start by noting that our general driven-dissipative resonator problem does not satisfy CQDB, and thus its correlation functions do not all exhibit a simple timesymmetry. An example of such a lack of correlation function symmetry is shown in Fig. 6. More generally, as discussed in Sec. II E, CQDB can only hold if the system's steady state commutes withĤ. This condition is violated except in the vanishing dissipation limit κ 1 , κ 2 → 0 + .
Despite the lack of CQDB, these systems always possess hidden-TRS, which explains their solvability. The specific nature however of the symmetry operator (or operators) depends on the particular version of the model. Consider first the most common case, where there is no two-photon loss, κ 2 = 0. To determine whether our system has hidden-TRS, we consider a doubled two-cavity system and a two-cavity state |ψ T . The question is whether this state could represent a TFD state constructed using an anti-unitary operatorT which describes a hidden-TRS (c.f. Eq. (13)). From Eqs. (45)-(47), such a state must satisfy: for some real energy E and constant u = ±1. Here (as always)Ĥ eff =Ĥ −iκ 1â †â /2 is the effective non-Hermitian Hamiltonian in our master equation. If we can find a twocavity state |ψ T satisfying the above equations, then we are guaranteed both to have hidden TRS, and to be able to solve for our system using the CQA method of Sec. V.
If we have a non-zero single-photon drive Λ 1 , one can only solve Eqs. (58)-(59) if u = 1 and E = 0. With these choices, there is a unique solution for the two-cavity state |ψ T . This was explicitly found and expressed as a confluent hypergeometric function in Ref. [18], which demonstrated that this model can be solved using CQA. Hence, the system has a unique hidden-TRS operatorT in this case. As we will see in the next section, this gives us more than just a way to understand the solvability of the model: it also directly lets us predict a surprising correlation function symmetry.
It is also interesting to consider the special case where there is no single-photon drive, Λ 1 → 0. Because of the single photon loss, the system still has a unique steady state. However, there are now two distinct hidden-TRS symmetriesT ± , each corresponding to distinct TFD states |ψ ± T : We stress that both these states each yield the sameρ ss when the auxiliary second cavity is traced out. Formally, the two TFD states (and correspondingT ± ) are found by solving Eqs. (58)-(59) for E = 0, u = 1 and E = 0, u = −1. The explicit states can be found analytically in terms of Bessel functions [18].
We thus have our first example of a physical system with multiple, distinct hidden TRS symmetries; other examples are listed in Table I. Using the explicit forms of the TFD states, we can explicitly compute the action of the hidden TRS symmetry operatorsT + andT − . In general, their action is highly non-trivial (as can be seen in Fig. 5, where we show their action in phase space on an initial coherent state).
In the limit of vanishing nonlinearity K → 0, the hidden-TRS operatorsT σ take a simple form. In this case, the two TFD states limit to simple two-mode squeezed states: where ∆ eff ≡ ∆ + iκ 1 /2 and |0, 0 is the two-cavity vacuum state. Expanding out the exponential allows us to pick out the corresponding hidden time-reversal operators, which correspond to simple phase-space reflections about the axes θ = ± arg(Λ 2 /∆ eff ) in phase space: where here,K denotes complex-conjugation in the Fock basis. For non-zero Kerr, the corresponding time-reversal operations become highly nontrivial and non-Gaussian, and must be extracted via a numerical Schmidt decomposition. In Fig. 5, we show the action ofT ± for both K = 5 × 10 −4 κ 1 (weak nonlinearity) and K = κ 1 (strong nonlinearity).

B. Experimental consequences of hidden-TRS
Our finding that driven-dissipative nonlinear cavities possess a hidden-TRS does more than simply explain why these systems are exactly solvable: it also lets us predict observable phenomena that are accessible in a standard single-system experiment. Recall our discussion in Sec. III C: while hidden-TRS (by definition) guarantees a symmetry of doubled-system TFD correlation functions, for certain operators, this directly implies timesymmetry of standard, single-system correlators. In particular, these special operators are ones that transform simply under the exchange operator J . By virtue of Eqs. (45)- (47), the effective HamiltonianĤ eff and jump operatorsĉ k are guaranteed to be such special operators.
As a specific example, consider the following steady-  Figure 6. Time symmetry of special correlation functions in a driven Kerr resonator. Real part of the connected, steady-state correlation function C a 3 ,a (t) ≡ â 3 (t)â (c.f. Eq. (63)) for a parametrically driven nonlinear cavity with Λ2 = K, κ1 = 0.4K and κ2 = Λ1 = ∆ = 0. This correlation function is symmetric in t, something that is guaranteed by the existence of hidden TRS. We also plot another quartic correlation function C X 2 P 2 (t) (whereX,P are canonical quadrature operators). This correlator is clearly asymmetric as a function of time. Hidden-TRS only ensures that a certain restricted class of correlators are time symmetric (in contrast to the more commonly studied CQDB which guarantees all correlators exhibit a form of time-symmetry).
state, single-system correlation function: If we set κ 2 = 0, Eqs. (45)-(47) ensure that J [â] = ±â. From the definition of the exchange operator, it follows that J [â m ] = (−1) mâm . As a result, hidden-TRS guarantees (via Eq. (29)) the above correlator has an Onsagerlike time symmetry: C a 3 ,a (t) = C a 3 ,a (−t). We stress that this correlation function symmetry is special: unlike the case with CQDB, most correlation functions will not exhibit any time-symmetry. This behavior is shown explicitly in Fig. 6, where we contrast the correlator C a 3 ,a (t) (time-symmetric) with a more generic correlator involving quadrature operatorsX = (â +â † )/ √ 2, Hidden-TRS does not enforce any special symmetry of this latter correlator; hence, as expected, it is manifestly not symmetric in time. We stress that even though J [â 2 ] is simple, this does not imply that J [â †â ] is simple.
We thus have a clear experimental test for confirming the existence of hidden-TRS in this class of systems.

VII. BREAKING OF HIDDEN-TRS BY THERMAL FLUCTUATIONS AND INTERACTIONS
We have now demonstrated that hidden-TRS holds in two very different zero-temperature dissipative models: a Rabi-driven qubit with loss (Sec. III D), and a driven nonlinear cavity with one and possibly two photon loss processes (Sec. VI). Within the setting of our Lindblad master equations, zero temperature corresponds to only having dissipators that remove (and not add) excitations.
The natural next question is to ask what happens to hidden-TRS if we introduce a non-zero temperature to the above systems. This corresponds to adding dissipative processes that can add excitations. We show that in a generic setting where there is both coherent (Hamiltonian) driving as well as nonlinearity, introducing such thermal dissipators can break hidden-TRS. The only exceptions to this are the case of no driving (where the system is effectively classical), or the case of no nonlinearity (where the steady state is Gaussian). Our results here suggest that for a generic nonlinear drivendissipative system, hidden-TRS is a symmetry associated with vacuum fluctuations, and hence only emerges as one approaches the zero-temperature limit.
Our work here is inspired by and complements seminal studies from Dykman and co-workers of related phenomena in driven nonlinear oscillators [24][25][26][27][28][29][30]. These works studied the basic nonlinear resonator model of Eqs. (56)-(57) in the limit of weak dissipation, where the quantum master equation can be reduced to a simpler Pauli master equation (i.e. one can drop off-diagonal elements of the density matrix in the energy eigenstate basis). The resulting classical master equation was found to satisfy the classical detailed balance condition of Eq. (2) at zero-temperature; in a semiclassical limit, this could be shown analytically. Further, it was shown that this classical detailed balance failed to hold at non-zero temperatures, and that in the semiclassical limit, the corresponding transition temperature became exponentially small. Our work extends these results: by formulating detailed balance in completely quantum manner using hidden TRS, we are not limited to weak-damping or semiclassical regimes. We also discuss how the breaking of hidden-TRS by thermal fluctuations is contingent on having driving and nonlinearity; without both these ingredients, there is no symmetry breaking. Finally, we discuss how this symmetry breaking could be directly probed in experiment by measuring the time-symmetry of correlation functions.
wheren th represents the bath thermal occupancy at the qubit frequency.

Thermal dissipation with no drive
In the absence of a Rabi drive (i.e. Ω = 0), the unique steady state of our master equation has the thermal equilibrium form: where |g , |e denoteσ z eigenstates. This steady state commutes withĤ, and it is easy to confirm that the system has CQDB. Due to the lack of coherences, the problem is analogous to a classical two-state system; hence, the presence of detailed balance is not surprising. Formally, the system still possesses a set of hidden-TRS symmetries; this symmetry is however not unique. There is a one parameter family of hidden TRS operatorŝ T = (e iψ |e e| + |g g|)K z .
The presence of hidden TRS in the finite-temperature, undriven qubit system implies that it can be solved using the coherent absorber method of Sec. V B. The qubit-plus-absorber system has the cascaded Hamiltonian Eq. (55) whereĤ A is the qubit HamiltonianĤ = 1 2 ω 0σz acting on the physical qubit A andĤ B is the Hamiltonian acting on the auxiliary qubit B (the absorber). The cascaded system also has the collective jump operatorsĈ The pure state which is simultaneously dark with respect toĤ AB ,Ĉ 1 , andĈ 2 is After tracing out the absorber system, the single site steady state density matrix is precisely Eq. (65).

Thermal dissipation with a non-zero drive
We now ask what happens to our thermal qubit when a non-zero drive is added (Ω = 0). For simplicity we take ∆ = 0 (resonant driving), and define the dimensionless driving b ≡ Ω/κ(1 + 2n th ). With this definition, the steady stateρ ss,T of the driven qubit with thermal dissipation is given by the zero temperature result in Eq. (C2) with the simple substitution b → b .
Furthermore the eigensystem of the Liouvillian at finite temperature is obtained from the zero temperature results in Eqs. (C5)-(C8) by replacementsρ ss →ρ ss,T , b → b , and κ → κ(1 + 2n th ). Finally, the permissible TRS of the finite temperature system are given by Eq. (31) We consider the TFD correlation function C TFD yz (t) defined in Eq. (16) for Pauli operatorsσ y ,σ z . As in the zero-temperature case, we decompose this into classical correlation and the entanglement correction using Eq. (19), and we look at the time asymmetry of each. At finite temperature, the classical correlation asymmetry picks up new temperature dependent terms: where we have defined In contrast, the ψ-dependent entanglement correction is which also gains temperature-dependent terms. At zero temperature the above expressions reduce to Eqs. (32)-(33) so for ψ = π, C TFD yz (t) has time symmetry, and our system has a hidden-TRS. However, as soon asn th is non-zero, hidden-TRS is broken. For non-zero temperature, there is choice of ψ for which the timeasymmetry of the correlation function. To see this explicitly, suppose we set the time asymmetry of C TFD yz (t) to zero and attempt to solve for cos ψ. We obtain cos ψ For any finite temperature η > 0 so the right-hand side has a magnitude greater than 1, thus there is no solution for ψ. This shows explicitly that for finite temperature and non-zero driving, the driven-dissipative qubit problem has no hidden-TRS. This breaking of correlation function time symmetry is shown explicitly in Fig. 7. At a heuristic level, this symmetry breaking is a result of the classical contribution growing faster withn th than the entanglement contribution, breaking the cancellation that occurs atn th = 0.

B. Parametrically-driven nonlinear cavity at finite temperature
The qubit example above corresponds to a system where the strength of the nonlinearity is essentially infinite. We now consider a driven-dissipative system where the strength of the nonlinearity is tuneable: the parametrically-driven damped nonlinear cavity of Sec. VI, but now with thermal dissipation. For simplicity, we consider where there is only a parametric (twophoton) drive, and there are only single-photon dissipation processes. We then have a Lindblad master equation with a Hamiltonian given by Eq. (56) (with Λ 1 = ∆ = 0),  80)), which is guaranteed to be symmetric if hidden-TRS holds. There is a sudden onset of asymmetry above a threshold temperature, indicating a sharp temperature at which hidden-TRS is broken. In contrast, we also plot the total time asymmetry of a correlation function whose behavior is not constrained by hidden-TRS, function C X 2 ,P 2 (t) (dashed red curve); hereX andP are standard quadrature operators. This correlator is asymmetric already at zero temperature, and shows no strong temperature dependence. Inset: real-part of the correlation function asymmetry of C(t) ≡ C a 2 ,H eff (t) forn th = 0 (solid green), andn th = 0.2 (dashed green). For all plots we take Λ2 = 3K, κ1 = 0.01K, ∆ = Λ1 = κ2 = 0. and dissipatorŝ As is standard in the derivation of quantum optics master equations, the thermal occupancyn th corresponds to a Bose-Einstein factor evaluated at the bath temperature T and cavity resonance frequency ω cav :n th = 1/ (exp[ ω cav /k B T ] − 1). For simplicity, we ignore twophoton dissipative processes (i.e. κ 2 = 0). At non-zero temperature, hidden TRS is broken unless one or both of the nonlinearity K or drive Λ 2 are zero. To see this, note that for hidden-TRS to hold, Eq. (46) requires that for some two-mode state |ψ T and 2 × 2 involutory unitary matrix U the jump operators must satisfy The equations imply that |ψ T is annihilated by two independent canonical annihilation operators. As such, this state must be Gaussian, which in turn implies thatρ ss must be Gaussian. However, this steady state is Gaussian only if one or both of Λ 2 and K are zero. We thus have an important conclusion: the combination of thermal fluctuations, driving and nonlinearity together can break hidden-TRS. Note that for a more explicit proof that hidden-TRS does not hold, one can explicitly try to solve both Eqs. (45) and (46); for both Λ 2 and K nonzero, one can confirm that it is impossible to solve these equations.
It is interesting to consider the simple case of an undriven, linear thermal cavity (i.e. K = 0, Λ 2 = 0, n th > 0). In this case, the steady state is essentially classical (no Fock-state coherences), and it is well known that CQDB holds [2]. Formally, our system also has hidden-TRS, implying that this system can be solved using the CQA method. This can be explicitly shown by solving Eqs. (45)- (47). We find that solutions are possible when E = 0 and when U is chosen to have the form Here, θ is a real parameter. We thus have a continuous family of distinct hidden-TRS operatorsT θ (see Tab. I), in contrast to the pair of hidden time-reversal operatorŝ T σ seen for nonzero parametric driving and zero temperature (i.e. K = 0, Λ 1 = 0,n th = 0).
For each possible hidden-TRS operatorT θ , we have a corresponding thermofield double state. These always have the form of Gaussian, two-mode squeezed states: Returning now to the more interesting case where we add both parametric driving and nonlinearity, we can study how thermal fluctuations break the hidden-TRS that is present at zero temperature. We focus on an experimentally-accessible quantity that shows this symmetry breaking: the time-asymmetry of the steady-state correlation function where as always,Ĥ eff is the non-Hermitian effective Hamiltonian associated with our master equation. As discussed above, at zero temperature hidden-TRS guarantees that this special correlator has time-symmetry. This time-symmetry is lost asn th is increased.
To see this explicitly, we plot the total time-asymmetry Potential conditions (U must be unity) Table II. Dictionary connecting hidden TRS and an effective classical notion of TRS in the effective phase space used in the complex-P function method. We list objects/conditions commonly appearing in the solution of quantum master equations via quantum detailed balance (i.e. hidden TRS), and their counterparts in the language of classical detailed balance in the complex-P representation. This correspondence only exists for multimode bosonic systems coupled to local zero-temperature dissipation.
vs. temperature, which we define as: As shown in Fig. 8, the total asymmetry remains zero as long as the temperature T is small, but then suddenly jumps at a critical "transition" temperature, consistent with a sudden breaking of detailed balance.
This temperature at this transition can be understood heuristically as corresponding to having the thermal excitation rate κn th be comparable to the dissipative gap of the zero-temperature system. This dissipative gap Γ (i.e. slow relaxation rate) corresponds to switching between two coherent states |±α with |α| = (Λ 2 2 − κ 2 /4)/K 2 [48]. One finds that Γ is exponentially small due to the small overlap of these coherent states [45]: Setting this rate equal to κn th using the parameters in Fig. 8 yields a temperature k B T / ω c ≈ 0.2; this is consistent with the temperature scale for hidden-TRS breaking.
We stress that even at zero temperature, most system correlation functions do not exhibit any time-symmetry. Such correlation functions do not show any dramatic behavior as a function of temperature. As an example, we plot the asymmetry of the correlation function C X 2 ,P 2 (t), defined as m X 2 P 2 (T ) ≡ˆ∞ 0 |C X 2 ,P 2 (t) − C X 2 ,P 2 (−t)|dt, (83) in Fig. 8; here,X,P are standard quadrature operators.

VIII. HIDDEN TRS AND PHASE-SPACE METHODS: A QUANTUM-CLASSICAL CORRESPONDENCE
In this final section, we turn to driven-dissipative systems comprised of one or more bosonic modes, and connect our notion of hidden TRS to phase space methods that are well known in quantum optics, and have been used to solve non-trivial problems using an effective Fokker-Planck equation in an expanded phase space. The focus is on Lindblad master equations of the form Here, [â j ,â † k ] = δ jk is a standard set of independent bosonic modes, κ j represents a decay rate for each mode, andĤ is an arbitrary bosonic quantum many-body Hamiltonian. We will establish that for this restricted class of models, the fully quantum notion of hidden-TRS (described by an anti-unitary operatorT ) coincides with a classical notion of time-reversal in an expanded phase space, i.e. an involution of the form (x, p) → (x,p) (where x, p are classical phase space coordinates). Hence, the effective detailed balance properties (and potential conditions) of a complex-P Fokker-Planck equation is directly tied to hidden-TRS. This allows us to directly understand the success of the complex-P method in solving several non-trivial driven cavity problems [19]. Interestingly, we show that for extended models, this correspondence no longer necessarily holds. For example, by simply adding two-photon loss processes, there exist hidden-TRS operators that have no correspondence to a simple operation in an extended phase space.
The context of our discussion will be the complex-P phase-space representation of the general bosonic master equation in Eq. (84). This is a particular non-diagonal expansion of the system's density matrix in terms of coherent states that can be used to convert the master equation into a well-behaved, Fokker-Planck-like equation (see [51] for a pedagogical introduction). Consider the single-mode case first for simplicity. We consider a doubled phase space described by complex coordinates (α, β), and chose appropriate integration contours C, C for each of these variables. This lets us express the density matrix aŝ where P (α, β) is the complex-P quasi-distribution function. Using standard techniques [19], one can often convert the Lindblad master equation forρ into a Fokker-Planck equation for this function, which is required to be nonsingular on the integration surface C × C defined by the contours. The resulting equation has the standard form: Here, C µ (α, β) represents a generalized drift vector, and D µν (α, β) a generalized diffusion tensor. The above derivatives are holomorphic derivatives [19], and Einstein summation notation is implied.
We can now state our quantum-classical correspondence principle: if the quantum master equation Eq. (84) has a hidden quantum time-reversal symmetryT and corresponds to a well-defined Fokker-Planck equation in the complex-P representation, then this associated Fokker-Planck equation has a well-defined classical TRS corresponding toT . In the case where this classical TRS operation is trivial (i.e. the identity operation), this symmetry then correspond to a standard detailed balance condition, meaning that the drift and diffusion matrices satisfy potential conditions [1]. This directly enables an efficient solution for the steady-state distribution function P (α, β), and hence the steady-state density matrix.
The fact that the complex-P Fokker-Planck equations satisfy potential conditions is precisely the property that enabled exact solutions of a variety of nonlinear driven cavity models [20]. Our result shows that this surprising property is directly tied to a more general, and fully quantum symmetry: hidden-TRS. In what follows, we will describe precisely how to construct the classical timereversal operator corresponding to a hidden TRST , and then show how this correspondence can be broken by considering higher-order Markovian loss channels.

A. Detailed balance in generalized P -representations
We start by defining a notion of time-reversal symmetry that is meaningful for the complex-P distribution function P (α, β). We stress the well-known fact that this distribution function is in general complex valued, and thus does not represent a meaningful probability distribution. Nonetheless, we can formally use it to define quantities analogous to expectation values and correlation functions.
Note first that the expectation value of a holomorphic function X(α, β) defined on our complex phase space is defined as: We can also define a time-evolved function X(α, β; t) de-fined by the solution of the dual Fokker-Planck equation With these definitions in hand, we define time-reversal symmetry in our doubled, complex classical phase space as the existence of a phase-space involution such that all two-time correlation functions (calculated as defined above) are time-symmetric: Here X, Y are any two holomorphic functions, and the time-reversed functions are given as where (α, β) is the time-reversed counterpart to (α, β), i.e. another point in the integration surface on which the complex Fokker-Planck evolution is taking place. All we require is that this time-reversal operation squares to the identity, namely, time-reversing a point twice recovers the original point on C × C .
For complex-P distributions, one can establish a generalization of a standard result in classical probability theory, which we rigorously establish in Appendix J: in the limit where the time-reversal operation is just the identity, the classical detailed balance condition Eq. (90) is equivalent to the potential conditions on the Fokker-Planck equation [1,52]. Recall that these conditions correspond to the having the (pseudo)probability current vanish in the steady state at every point in phase space, where the pseudoprobability current J µ (α, β, t) is defined by rewriting the Fokker-Planck equation as a continuity equation: This constraint allows a direct method for solving for the steady state in terms of a potential function.
Note that in the case where the time-reversal operation (α, β) → (α, β) is not the identity, there is no simple potential-condition method for solving Fokker-Planck equations, unless the time-reversal symmetry is known beforehand (see e.g. Ref. [52].)

B. Constructing the classical TRS corresponding to a hidden TRST
We briefly outline the correspondence here in the single-mode case, and leave the discussion of the multimode case to App. F. For these systems, hidden TRS implies (among other constraints) the constraint where J is the exchange superoperator as always, and u = ±1. In App. F we show that, under the assumption that Eq. (93) holds, hidden TRS is equivalent to classical detailed balance in the complex-P representation with respect to the following classical time-reversal operation: Note that the fact that u squares to one, an intrinsic property of the exchange superoperator, ensures that this is a valid classical time-reversal operation in the effective phase space used for the complex-P representation. This surprising correspondence is the result of TFD correlation functions of normally-ordered operators coinciding with complex-P correlations. More precisely, if J [â] = uâ with u = ±1, then we have (see App. F): where X, Y are the classical representatives ofX,Ŷ in the complex-P representation. Explicitly, without loss of generality,X,Ŷ can be written aŝ In terms of the normally-ordered expressions above, the classical representatives X, Y have the following form: Finally, the classical time-reversal operation used to define the reversibility of the Fokker-Planck equation is given in Eq. (94). Therefore, in this situation, hidden (quantum) TRS is equivalent to classical TRS in the complex-P representation.
C. Breakdown of the correspondence principle: going beyond phase-space methods The simplest situations in which the above correspondence principle breaks down is in systems with higher-order loss dissipators, e.g. a system of the form An example of this is the driven cavity problem considered in Sec. VI, in the regime where all single-photon terms are set to zero: Λ 1 = κ 1 = 0. We are left with a model with an interaction, detuning, two-photon drive and two-photon loss. The full master equation in this case conserved photon number parity, and thus does not have a unique steady state. For this model, we still have hidden-TRS for each of the steady states. The full set of hidden-TRS compatible TFD states, i.e. obtained by solving Eqs. (45)- (47), has the form: where the individual terms ψ ± T in the superposition correspond to the two simple thermofield doubled states encountered in Sec. VI, and which, upon tracing-out the auxiliary cavity B, correspond to a single steady state of the master equation Eq. (100): According to the quantum-classical correspondence outlined in this section, this stationary state corresponds to a stationary complex-P distribution with both a trivial TRS (corresponding to the potential conditions) and an inversion TRS (corresponding to something more complicated): In light of the above observation, one might interpret the more exotic thermofield doubled state |ψ(γ, ν) as describing a quantum TRS which corresponds to a "superposition" of both the trivial and inversion TRS, and thus this hidden time-reversal symmetry no longer has a classical analogue. Indeed, it is well-known that the stationary state obtained via solution of the potential conditions is not the only stationary state of the quantum master equation Eq. (100) [21,47]. Tracing out the auxiliary cavity B for arbitrary parameters γ, ν reveals a family of quantum steady states:ρ While anyT and TFD must yield an exchange superoperator J that acts simply onâ 2 (via Eq. (46)), the action onâ need not be simple. In fact, we only get a simple action when ν = 0 (γ = 0), in which case J [â] =â ( J [â] = −â). For the more general case, the identity Eq.
(93) is broken. The more complicated nature of the hidden TRS operator and the corresponding J implies that the steady states corresponding to |ψ T (γ, ν) cannot be easily found using the complex-P phase space solution method.

IX. SUMMARY & OUTLOOK
In this work, we have introduced a new symmetry that can exist in driven-dissipative systems described by a Lindblad master equation: hidden time-reversal symmetry. We have shown explicitly how this goes beyond the conventional definition of quantum detailed balance (CQDB) introduced by Agarwal [2]; crucially, hidden-TRS can exist in systems whose steady states have nonzero energy-eigenstate coherences, something that makes it impossible to have CQDB. While hidden-TRS is most naturally formulated in terms of a doubled system prepared in a thermofield double state, we demonstrated that it has a direct observable consequence: a certain class of single-system correlation functions are guaranteed to be time-symmetric. This is in contrast to CQDB, which requires all correlation functions to obey a timesymmetry. To illustrate our ideas, we have analyzed how several ubiquitous driven quantum systems (qubit and nonlinear cavity models) can have hidden-TRS despite not having CQDB.
Perhaps most importantly, we established how hidden TRS provides a powerful means to derive analytic solutions for non-trivial steady states of quantum master equations. In particular, hidden-TRS underlies both the coherent quantum absorber exact solution method [17,18], as well as phase space methods based on the complex-P function [19].
We hope that our results will lay the groundwork for many further fruitful studies exploiting hidden-TRS as a means to understand even more complex systems. This symmetry could provide an interesting means for finding non-trivial, exactly solvable many-body driven dissipative systems, both of bosons, qubits and possibly of fermions. It could also lead to novel perturbative techniques for studying systems that weakly break hidden-TRS. Finally, it would also be extremely interesting to rephrase this symmetry fully in terms of a dissipative field theory describing the system of interest (i.e. in terms of a Keldysh action [53][54][55]). This could yield further insights, and also perhaps enable an extension of these ideas into non-Markovian regimes. from the Simons Foundation via a Simons Investigator Award.
Appendix A: Doubled-system classical detailed balance In this section, we prove that classical detailed balance may be equivalently stated as the following symmetry condition: The reason why the above condition is equivalent to the standard definition of detailed balance is that the doubled-system correlation function X A (t)Y B (0) is actually a single-system correlation function in disguise: Therefore, the doubled-system correlation function Now, making the replacement Y →Ỹ in this single-site symmetry condition yields the definition of classical detailed balance used in the main text. Therefore, Eq. (A4) is equivalent to classical detailed balance: Note that we have implicitly used the fact that Y → Y is a bijection of the algebra of random variables. In summary, we have shown that the doubled definition Eq. (A1) of classical detailed balance is completely equivalent to the standard definition Eq. (A5).

Appendix B: CQDB rules out stationary coherences between energy eigenstates
We now demonstrate that systems with CQDB have steady state density matrices that are always guaranteed to be diagonal in the energy eigenbasis. There are many references that show this explicitly [4,5,35]. However, here we will assume an intermediate result, namely that CQDB implies modular symmetry [35], that is, a symmetry of the driven-dissipative dynamics with respect to the unitary dynamics generated by the modular Hamiltonian The reason for taking this symmetry-based perspective is that it informs most of the central results in the theory of quantum detailed balance [9,10]. Indeed, once the above symmetry is established, the proof that steady states with CQDB are diagonal in the energy eigenbasis is very easy. We provide here a simple argument that works in the finite-dimensional case. We begin with Lindblad's original expression for the effective Hamiltonian as a classical average [56]: where here, U is a Haar-random unitary, andL is the Heisenberg-picture Lindbladian, which generates timeevolution of observables. We now observe how the effective Hamiltonian evolves under the modular (dynamical) The above identity holds for any Lindbladian, and thus contains no hidden assumptions about the system in question. Now suppose, however, thatL satisfies CQDB, and thus has modular symmetry. Then, in particular, we haveρ Substituting the above identity into the Haar-random average defining the effective non-hermitian Hamiltonian yields a simple result -Ĥ eff is time-independent with respect to the unitary dynamics generated by the modular Hamiltonian: Therefore, the effective Hamiltonian commutes with the modular Hamiltonian, and thus the effective Hamiltonian commutes with the steady state itself: [Ĥ eff ,ρ ss ] = 0. By taking the anti-hermitian part of the commutator, we immediately have thatĤ commutes withρ ss .
Appendix C: Example of broken CQDB: dissipative Rabi-driven qubit To make the ideas of sections II D and II E more concrete, we consider simple but ubiquitous system which does not satisfy CQDB, does have hidden TRS.

Violation of CQDB via correlation function asymmetry
Our example is a qubit (with Pauli operatorsσ x,y,z ) subject to a coherent Rabi drive in the presence of loss. Working in the rotating frame set by the drive, and making a rotating wave approximation, the master equation Here ∆ is the detuning of the drive from the qubit splitting frequency, Ω is the Rabi frequency, and κ is the decay rate of the qubit excited state. The steady state for this system is easy to find and given in many textbooks, see e.g. Ref. [32]: Given the external driving, this steady state does not correspond to thermal equilibrium, and hence a priori there is no reason to expect that it will satisfy CQDB. While this may seem obvious, we will now show that CQDB is broken explicitly, by directly uncovering correlation function time-asymmetry in this system. We stress that the CQDB condition in Eq. (11) is contingent on the choice of time-reversal operatorT . We will take a general approach here (and throughout this paper): we do not preselect the definition ofT using on additional knowledge of our system, but rather ask where there is any possi-ble anti-unitaryT which would give rise to a symmetry. Hence to truly rule out CQDB, one must check Eq. (11) for all permissible choices ofT . We will thus show that CQDB does not hold no matter what choice is made for T .
In what follows, for simplicity we assume a resonant drive (i.e. ∆ = 0), and introduce the dimensionless Rabi frequency b ≡ Ω/κ; CQDB is violated even for nonresonant drives, see Sec. C 2 below. The first step is to parameterize all possible TRS operators. Since CQDB only holds ifρ ss is itself invariant under TRS, this constrains the form of TRS. The only permissible TRS operators are then parameterized by a single phase ψ and have the form: whereK z is the complex conjugation operator acting in theσ z basis.
To show that CQDB cannot hold in this system, it is sufficient to show that at least one correlation function violates Eq. (11) for each TRS angle ψ. As our main object of study, we introduce the correlation function defined for positive and negative times: whereσ j =Tσ jT −1 . CQDB holds, then Eq. (11) implies time symmetry: C yz (t) = C yz (−t). In what follows we show that the properties of the TRS and the eigenmodes of the Liouvillian do not allow for the time symmetry of C yz (t) and other correlation functions.
The Liouvillian of the driven qubit system is readily diagonalized. Letting λ j denote its eigenvalues andr j the corresponding right-eigenvectors, we find: where b ≡ Ω/κ, α = √ 16b 2 − 1 is the dimensionless damped Rabi frequency. For b < 1/4 the qubit is overdamped and the frequency is imaginary: The essential feature of the eigensystem is that theσ x coherence behaves differently from theσ y coherence or the classical population (σ z ). Theσ x coherence decays exponentially with rate (−λ 1 ) = κ/2. Thê σ y coherence and the classical populations decay with (−Re λ 2,3 ) > κ/2 for any finite b. This implies that σ x (t)σ k (0) has a time dependence that is always different from σ y (t)σ k (0) or σ z (t)σ k (0) for anyσ k .
With this in mind, we turn to the time-reversed operatorσ z which, in general is a linear combination of all three Pauli operators. In particular, theσ x component isσ We can see that for any ψ = 0, π, for which sin ψ = 0, the expression C yz (t < 0) will have terms measuring the decay ofσ x coherence. Thus the time dependence at t < 0 must be qualitatively different from that at t > 0. Even at ψ = 0, π the correlation function is generically not time symmetric. As an example, we show a generic plot of the time asymmetry of C yz (t) in Fig. 9 for b = 1 computed for various ψ. Although the argument breaks down at ψ = 0, π, we can show definitively that CQDB cannot hold by considering C xz (t), which is defined analogously with C yz (t). Since at ψ = 0, π there are noσ x terms iñ σ z , the time dependence of C xz (t < 0) must be qualitatively different from the time dependence of C xz (t > 0). Therefore we conclude that the Rabi-driven qubit does not satisfy CQDB.

Violation of CQDB for any detuning
In the preceding section we restricted our analysis to the resonantly driven qubit for which ∆ = 0 primarily because the diagonalization of the Liouvillian becomes unwieldy for ∆ = 0. Here we show by an alternate route that for any detuning, the Rabi driven system violates CQDB.
Recall from Appendix B that a system which satisfies CQDB must have a steady state that is diagonal in the energy eigenbasis. The commutator [Ĥ,ρ ss ] is equivalent to taking the cross product between the Hamiltonian "vector" (Ω/2, 0, ∆) and the traceless part ofρ ss . Imposing the constraint that the commutator is zero requires the form ofρ ss to bê where α is a real constant of proportionality. Crucially, the above expression is linear in ∆ and Ω, whereas the true steady state is a quadradic rational function of these parameters and the decay rate κ, c.f. Eq. (C2). Therefore we conclude that for any detuning, the Rabi-driven qubit does not satisfy CQDB. Furthermore, one can numerically diagonalize the Liouvillian at any drive detuning and verify that there does not exist any TRS for which the steady state is invariant and for which all correlation functions are time symmetric.

Permissible TRS
Here we show how the permissible TRS of Eq. (C3) are determined from the steady state Eq. (C2). We require only thatT leaves the steady state invariant as a necessary condition of CQDB. All possible TRS take the form T =VK ρ for unitaryV and complex conjugationK ρ in the steady state eigenbasis such thatTρ ssT −1 =ρ ss . This condition implies that the action ofT on the eigenstates ofρ ss (i.e. pointer states) is restricted to bê up to a global phase. Thus in the steady state eigenbasis the permissible TRS take the formT = (|1 1| + e iψ |2 2|)K ρ for any ψ, whereK ρ is complex conjugation in the steady state eigenbasis and the pointer states are assumed time-reversal invariant:K ρ |n = |n .
As the final step, we representT in the familiarσ z basis. Given the change of basis unitaryÛ that diagonalizes the steady state,Ûρ ssÛ † = diag(p 1 , p 2 ), the permissible TRS are given in theσ z basis aŝ For completenessÛ is given in terms of b ≡ Ω/κ and s ≡ √ 4b 2 + 1 aŝ Substituting this into the expression forT above recovers Eq. (C3) (which is also Eq. (31) in the main text).
Appendix D: Explicit construction of exchange superoperator J By definition, the exchange superoperator is supposed to act on a single site observableÔ and produce a new single-site observable, J [Ô], which, upon acting on site B of the thermofield double, produces the same state as one would obtain by acting on site A with the observablê O: One may interpret the above equation as the vectorization of an operator equation, according to the ruleŝ where |ψ K is the thermofield doubled state where the time-reversal operationT ≡K is complex-conjugation in the eigenbasis of the steady state. Explicitly: Furthermore,Ô →Ô T denotes matrix transposition in the eigenbasis of the steady state. An arbitrary time-reversal operation can be decomposed asT ≡VK, whereV commutes withρ ss . Explicitly, an arbitrary thermofield doubled state corresponding to a hidden TRS must always have the expression Under the above rules, the thermofield doubled state with an arbitrary time-reversal operation for a unitaryV satisfies the following identities: We can then proceed to verify Eq. (D1) directly: Taking advantage of the fact thatV commutes withρ ss , we havê which is just the operator representation ofÔ A |ψ T according to the rules Eqs. (D7-D8).

Appendix E: From CQDB to hidden TRS
In this appendix, we will show that, for systems with modular symmetry, CQDB and hidden TRS are equivalent. Since CQDB already implies modular symmetry by itself [5,35], this will demonstrate that CQDB necessarily implies hidden TRS. Therefore, CQDB is a strict subphenomenon of hidden TRS.
Indeed, consider an arbitrary driven-dissipative system described by a Lindblad master equation. The definition of the exchange superoperator tells us that The above correlation function, for negative times, is harder to rephrase as a two-site quantity: However, the above expression simplifies considerably in systems with modular symmetry. Indeed, let us now assume that the driven-dissipative system in question is symmetric with respect to the modular HamiltonianĤ ρ , as is always the case with systems that have CQDB [5]. Then, we can write (ρ ss . Substituting this identity into Eq. (E2), we get Therefore, for any driven-dissipative system with modular symmetry, Therefore, since J is a bijection of the observable algebra, CQDB and hidden TRS are equivalent for this class of systems, as time-symmetry of one set of correlation functions implies time-symmetry of the other.
steady state), is expressible using the same bilinear form: where the relevant bilinear form (which depends on a particular choice of time-reversal operatorT ) is defined in the main text as follows: Clearly, the condition of hidden TRS is equivalent to identifying the left-hand side of Eq. (F3) with the lefthand side of Eq. (F4), which, upon examining the righthand sides of said equations, is equivalent to the symmetry of the Heisenberg-picture time-evolution superoperator E t ≡ e −tL with respect to the bilinear form ·, · T .
One can form an analogous construction to express classical correlations in the complex-P distribution in terms of a bilinear form. Indeed, for any pair of normallyordered quantum observablesX,Ŷ and a classical timereversal operation X →X, one can write and for negative times (with a stationary complex-P distribution which is time-reversal invariant), we can also write: Where we have introduced a new symmetric bilinear form where X, Y are the normally-ordered symbols ofX,Ŷ , i.e. the classical representatives of the observablesX,Ŷ in the complex-P representation. Explicitly, ifX = IJ c IJ (â † 1 ) i1 · · · (â † n ) inâ j1 1 · · ·â jn n , then the normallyordered symbol X is where I ≡ {i 1 , · · · , i n } and J ≡ {j 1 , · · · , j n } are multiindices. The condition of classical detailed balance in the complex-P representation, i.e. the time-symmetry of Eq. (F1), is equivalent to identifying the left-hand side of Eq. (F6) with the left-hand side of Eq. (F7), which is equivalent to the symmetry of the Heisenberg-picture time-evolution superoperator E t ≡ e −tL with respect to the bilinear form ·, · P .
With these definitions in hand, we return to the general problem of interest: a Markovian multi-mode bosonic system where each mode is subject to loss. For simplicity, we start by assuming each mode has the same loss rate, implying a master equation of the form If this system has hidden TRS with respect to a particular quantum time-reversal operationT , then as discussed in the main text (c.f. Eq. (43)) there must exist a N × N unitary matrix U such that: Here (as always) J is the exchange superoperator. What we will show in this section is a remarkable coincidence between the two families of bilinear forms, under the above jump operator constraint. That is, Eq. (F11) implies that we can identify both bilinear forms, i.e.
where the classical time-reversal operation on the complex-P side is none other than the change-of-Kraus representation given in Eq.
In what follows, we first prove this result in the single mode case, then extend to the multi-mode case where each mode has an identical damping rate κ. Finally, we extend the result to the more general case where each mode has a different damping rate κ j .
a. Single-mode case Consider Eq. (F10) in the single mode limit N = 1. If this system has hidden TRS with respect to a particular quantum time-reversal operationT , then the jump operator constraint reduces to J [â] = uâ, with u = ±1 a scalar quantity. We show in what follows that this in turn implies the following identity: where the bilinear form on the right-hand side is defined using the classical time-reversal operation (α, β) = (uα, uβ). To see this, note that both the left-and righthand sides are bilinear with respect toX,Ŷ , and so it suffices to verify the above identity in a basis of normallyordered monomials. That is, without loss of generality, we may assume that We begin by computing, e.g.
Substituting-in the definition of the exchange superoperator, we get Finally, by direct computation, we also have that J [Ô † ] =ρ ss J [Ô] †ρ−1 ss for genericÔ, and so (via the cyclic nature of the trace) Immediately, we recognize here the normally-ordered symbols ofX,Ŷ , which are, explicitly: With the above observation, we have thus proved Eq. (F13), which establishes the equivalence of generalized quantum detailed balance, that is, hidden TRS, with classical detailed balance in the complex-P representation. Note that the potential conditions, as well as the original CQA method, as originally formulated in [17], both correspond in this context to the special case of a trivial TRS, i.e. U = 1.

b. Many-body case
Now, consider the multi-mode master equation Eq. (F10). If this system has hidden TRS with respect to a particular quantum time-reversal operationT , then we have the constraint Eq. (F11) on the jump operators, which we write as What we will demonstrate in this section is that this jump operator constraint implies the following identity: where the bilinear form on the right-hand side is defined using the classical time-reversal operation ( α, β) = (U α, U β). To see this, note that both the left-and righthand sides are bilinear with respect toX,Ŷ , and so it suffices to verify the above identity in a basis of normallyordered monomials. That is, without loss of generality, we may assume that X = (â † 1 ) k1 · · · (â † n ) knâl1 1 · · ·â ln n , (F22) Y = (â † 1 ) p1 · · · (â † n ) pnâq1 1 · · ·â qn n .
We begin by computing, e.g.
With the above observation, we have thus proved Eq. (F21), which establishes the equivalence of generalized quantum detailed balance, that is, hidden TRS, with classical detailed balance in the complex-P representation.
Finally, consider the general case where each mode has a different damping rate κ j . In this case, the corresponding classical TRS is modified to be which can be proven by trivially repeating the steps above.
Appendix G: Doubled system correlation function symmetry for TRS ψ = π In the main text we showed that the TFD correlation function C TFD yz (t) is time symmetric for the TRS ψ = π. We must still explicitly verify that the other two correlation functions, C TFD xy (t) and C TFD xz (t), are time symmetric for the same TRS.
The expansion of the TFD correlation functions as in Eq. (19) is practically useful for computing correlation functions if the time-dependent operators are easily found because it reduces the problem to computing matrix elements. For the qubit system, the time-dependent operators are readily written in terms of the Liouvillian eigenmodes asσ With these in hand we proceed to compute the matrix elements.
We thus arrive at an interesting result already: the doubled system classically correlated state has identically zero correlation functions C cl xy (t) = 0 = C cl xz (t). We therefore need to compute only the off-diagonal elements ofσ y (t) andσ z (t). Given the TRS ψ = π and the form of the pointer states, it is straightforward to compute these in the time-reversed pointer state basis. The relevant result is that they are real and thus − |σ y |+ = + |σ y |− , (G10) − |σ z |+ = + |σ z |− , as required by their hermiticity. To see why this must be true without explicitly computing the matrix elements, note that in the time-reversed pointer state basis,σ x (t) is off-diagonal and imaginary and thus plays the role of an effective "σ y " in this basis. Thereforeσ y andσ z must be linear combinations of the effective "σ x " and "σ z " and hence their off-diagonal matrix elements must be real and equal. Furthermore, we can draw the same conclusion about the matrix elements of the time-dependentσ y (t) and σ z (t): −|σ y (t)|+ = +|σ y (t)|− , (G12) −|σ z (t)|+ = +|σ z (t)|− .
(G13) Againσ x is off-diagonal and imaginary and becausê l 2 andl 3 have onlyσ y andσ y components, the timedependent operators must remain as linear combinations of the effective "σ x " and "σ z " which have real and equal off-diagonal matrix elements. Finally putting everything together, the TFD correlation functions are whereĤ eff is the effective Hamiltonian in our master equation, andĉ l are the jump operators. We denote the Hermitian (anti-Hermitian) parts of an operatorÂ as Re Â (iIm Â ). One can then tautologically rewrite Eq. (H1) as ishes.
Here, Σ is a 2n-dimensional closed integration surface (the many-body analogue of the pair of integration contours C, C mentioned in the main text). Furthermore, as the diffusion tensor D µν ( α, β) is generically complex, it may fail to be positive-definite. Despite these differences, in this appendix we will nonetheless prove that the standard classical result in [57] still applies (however, here there will be no assumption on the diffusion tensor D µν ( α, β)): the potential conditions are equivalent to detailed balance with respect to a trivial time-reversal operation: ( α, β) ≡ ( α, β). (J6) This result is critical for bosonic many-body systems with local onsite loss, as, with a nontrivial U parameter, hidden TRS corresponds to a nontrivial classical TRS in the complex-P representation, and thus potentially corresponds to a richer symmetry than that encapsulated by the potential conditions. We first establish the forward implication: if the stationary pseudoprobability current Eq. (J4) vanishes, then a straightforward calculation shows that L(XP ss ) = P ss L * (X), where L is the Liouvillian for the pseudo Fokker-Planck equation Eq. (J3), and L * is the adjoint Liouvillian, obtained via integration by parts [19]. The proof almost exactly follows [57], except while carrying out the calculation we notice that positive-definiteness of the diffusion tensor D µν ( α, β) is not needed in the proof, and thus the assumption of positive-definiteness may be relaxed. From Eq. (J7), detailed balance is immediate, as, by integrating by parts on the surface Σ (which is valid, as X, Y are holomorphic), the dual Liouvillian L * can then be shown to be symmetric with respect to the bilinear form X, Y ≡ˆΣ d n αd n β P ss XY. (J8) From this and the expressions Y (t)X(0) = X, e −tL * (Y ) , we have detailed balance. For the converse direction, the proof follows analogously to the classical case as well: one starts by proving the following (assuming P ss is nonvanishing on Σ): The above is a generalization of the lemma in Section 4.6 of [57], to pseudo Fokker-Planck equations, and may be proven by invoking the holomorphicity of X, Y . With the above, one can compute the asymmetry of the Liouvillian: Now, if detailed balance holds with respect to a trivial TRS, then the Liouvillian L * is symmetric, and thus the right-hand side must vanish for all pairs of holomorphic functions X, Y . The only way that this can happen is if the stationary current vanishes everywhere on the integration surface, i.e.
(J13) However, the above is precisely the statement of the potential conditions.