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

David Roberts1,2, Andrew Lingenfelter1,2, A. A. Clerk  Pritzker School of
Molecular Engineering, University of Chicago, Chicago, IL, USA
2Department of Physics, University of Chicago, Chicago, IL, USA

(August 27, 2024)

###### Abstract

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 steady-state two-time correlation functions
$\overline{A(t)B(0)}$: they must be invariant when each quantity is replaced
by its time-reversed version and $t\rightarrow-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 Gardiner (2009)).

Figure 1: (a) The concept of hidden time-reversal symmetry helps unify
disparate-seeming topics in the study of driven-dissipative quantum systems.
(b) A driven dissipative system described by a Lindblad master equation. The
conventional definition of quantum detailed balance is formulated as a time
symmetry of steady-state correlation functions. (c) Hidden TRS is a symmetry
ensuring that correlation functions of a doubled version of the original
system, prepared in a thermofield double state, are time symmetric. It is
defined by an anti-unitary operator $\hat{T}$. Hidden TRS can hold even when
the correlation function symmetry in (b) fails; it also enables a powerful
method for finding exact solutions.

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 Agarwal (1973); Carmichael and Walls (1976); Alicki
(1976); Kossakowski  _et al._ (1977); Majewski (1984); Majewski and Streater
(1999); Denisov  _et al._ (2002); Fagnola and Umanita (2007); Fagnola and
Umanità (2010); Duvenhage and Snyman (2018); Carlen and Maas (2017); Ramezani
_et al._ (2018). The most natural definition requires that steady-state
correlation functions in the quantum theory obey an Onsager symmetry analogous
to the classical case Agarwal (1973); Carmichael and Walls (1976); Denisov
_et al._ (2002); this condition necessarily holds if the microscopic system-
bath dynamics obey time-reversal symmetry Carmichael and Walls (1976). Later
more formal works considered generalized definitions of quantum detailed
balance Fagnola and Umanita (2007); Fagnola and Umanità (2010), 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 matrix111 The
only attempt at such connections in the past were limited to systems that were
easily solvable by other means (e.g. linear bosonic systems, or systems that
could be reduced to a classical master equation) Agarwal (1973)..

In this work, we introduce a powerful, symmetry-based formulation of quantum
detailed balance (QDB) that goes beyond the simple definition in Ref. Agarwal
(1973), and that directly enables an efficient way for finding non-trivial
steady states. Our work builds on Ref. Duvenhage and Snyman (2018), which
showed that a particular generalized definition of QDB introduced in Ref.
Goldstein and Lindsay (1995) can be formulated using an entangled, thermofield
double state Takahashi and Umezawa (1996). 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. Agarwal (1973) (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. Stannigel  _et al._
(2012), and extended in Ref. Roberts and Clerk (2020). Hidden-TRS is also
connected to well-known exact solution methods from quantum optics based on
the complex-$P$ phase space quasiprobability Drummond and Gardiner (1980);
Drummond and Walls (1980); Bartolo  _et al._ (2016); Elliott and Ginossar
(2016): 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 Zhu and Monroe (2020), 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 1): Rabi-driven
qubits subject to dissipation, and driven-dissipative nonlinear quantum
cavities. These systems exhibit in general no correlation function time-
symmetry, and hence do not possess CQDB as defined in Agarwal (1973). 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 Drummond and Gardiner (1980); Drummond and Walls
(1980); Bartolo  _et al._ (2016); Elliott and Ginossar (2016); Roberts and
Clerk (2020). 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 Dykman and Krivoglaz (1979); Dykman and Smelyanskii
(1988); Marthaler and Dykman (2006); Dykman (2012); Guo  _et al._ (2013); Guo
(2013); Zhang and Dykman (2019).

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 Agarwal (1973); Carmichael and Walls (1976); we show that CQDB only holds
for a limited class of systems that possess trivial steady states. Sec. III
introduces the notion of hidden TRS, and connects it to the definition of
generalized QDB introduced in Ref. Goldstein and Lindsay (1995). 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.

|  | Hidden | Hidden TRS:
---|---|---|---
System | CQDB | TRS | unique?
Thermal qubit | Yes | Yes | $\hat{T}_{\theta},\leavevmode\nobreak\ e^{i\theta}\equiv 1$
Thermal linear cavity | Yes | Yes | $\hat{T}_{\theta},\leavevmode\nobreak\ e^{i\theta}\equiv 1$
Kerr cavity, 1-ph. drive | No | Yes | Yes
Kerr cavity, 2-ph. drive | No | Yes | No; $\hat{T}_{+},\hat{T}_{-}$
with nonzero temp. | No | No | N/A
Driven qubit | No | Yes | Yes
with non-zero temp. | No | No | N/A
Table 1: Common driven-dissipative quantum systems and their status both with
respect to conventional quantum detailed balance (CQDB) (c.f. Sec. II.4), and
our new notion of hidden-TRS. Italics indicate systems with hidden-TRS that do
not have CQDB. Some of these systems possess multiple distinct hidden TRS
(right-most column).

##  II Classical detailed balance and conventional quantum detailed balance

###  II.1 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 $\Gamma_{n\to m}$. The time-dependent probability $p(n,t)$
for the system to be in a given state $n$ then obeys:

$\displaystyle\frac{d}{dt}p(n,t)=\sum_{m}p(m,t)\Gamma_{m\to
n}-\sum_{m}p(n,t)\Gamma_{n\to m}.$ (1)

We assume that this equation admits a time-independent steady-state
probability distribution $\bar{p}(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 $\tilde{n}$
denote the time-reversed version of the microstate $n$, the condition is (see
e.g. Gardiner (2009)):

$\displaystyle\bar{p}(n)\Gamma_{n\to
m}=\bar{p}(\tilde{m})\Gamma_{\tilde{m}\to\tilde{n}}.$ (2)

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 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. Gardiner (2009)).

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:

$\overline{X(t)Y(0)}\equiv\sum_{m,n}X(m)p(m,t|n,0)Y(n)\bar{p}(n)$ (3)

The detailed balance condition of Eq. (2) is then equivalent to requiring that
the following symmetry hold for all steady state correlation functions:

$\displaystyle\overline{X(t)Y(0)}=\overline{\tilde{Y}(t)\tilde{X}(0)}.$ (4)

Here, the time-reversed function $\tilde{X}$ is defined as $\tilde{X}(n)\equiv
X(\tilde{n})$.

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.

###  II.2 Doubled-system formulation of classical detailed balance

Consider now an alternate but equivalent formulation of classical detailed
balance Duvenhage and Snyman (2018). 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:

$\displaystyle\frac{d}{dt}p_{AB}(n,n^{\prime};t)$
$\displaystyle=\sum_{m}p_{AB}(m,n^{\prime};t)\Gamma_{m\to n}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ -\sum_{m}p_{AB}(n,n^{\prime};t)\Gamma_{n\to m}.$ (5)

We further assume that the doubled system is initially prepared in a
correlated state described by the probability distribution:

$\displaystyle p_{AB}(n,m;0)=\bar{p}(n)\delta_{m,\tilde{n}},$ (6)

where $\bar{p}(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 $\tilde{n}$ 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 to $\bar{p}(n)$ for system $A$, and
$\bar{p}(\tilde{n})$ 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:

$\displaystyle\overline{X_{A}(t)Y_{B}(0)}=\overline{Y_{A}(t)X_{B}(0)},$ (7)

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. Duvenhage and Snyman (2018) 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.

###  II.3 Markovian quantum open system: general setting

Onsager-like correlation function time-symmetries can also be considered in
the context of open quantum systems Agarwal (1973); Carmichael and Walls
(1976); Denisov  _et al._ (2002). 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 $\hat{H}$
and a set of jump operators $\hat{c}_{l}$ that describe the influence of
external dissipative reservoirs. Defining
$\hat{H}_{\text{eff}}=\hat{H}-\frac{i}{2}\sum_{l}\hat{c}_{l}^{\dagger}\hat{c}_{l}$,
our general Lindblad equation of motion for the system density matrix
$\hat{\rho}(t)$ is:

$\displaystyle\frac{d}{dt}\hat{\rho}=-i(\hat{H}_{\text{eff}}\hat{\rho}-\hat{\rho}\hat{H}_{\text{eff}}^{\dagger})+\sum_{l=1}^{M}\hat{c}_{l}\hat{\rho}\hat{c}_{l}^{\dagger}\equiv\mathcal{L}[\hat{\rho}],$
(8)

where we have introduced the Liouvillian superoperator $\mathcal{L}$. Given a
particular steady state of this equation (i.e. a time-independent solution
$\hat{\rho}_{\rm ss}$), steady-state correlators between two system operators
$\hat{X}$ and $\hat{Y}$ obey time-translation symmetry and can be calculated
from the master equation (see, e.g., Breuer and Petruccione (2002)). To state
this compactly, we first introduce the adjoint Liouvillian superoperator
$\bar{\mathcal{L}}$, determined by (see e.g. Ref. Gardiner and Zoller (2000))

$\displaystyle\bar{\mathcal{L}}[\hat{A}]\equiv
i(\hat{H}_{\text{eff}}^{\dagger}\hat{A}-\hat{A}\hat{H}_{\text{eff}})+\sum_{l=1}^{M}\hat{c}_{l}^{\dagger}\hat{A}\hat{c}_{l}.$
(9)

Correlation functions are then computed as

$\displaystyle\langle\hat{X}(t)\hat{Y}(0)\rangle\equiv\textrm{Tr}\left(\mathcal{E}_{t}[\hat{X}]\hat{Y}\hat{\rho}_{\rm
ss}\right),\,\,\,\,\,\mathcal{E}_{t}\equiv\exp\left(\bar{\mathcal{L}}t\right).$
(10)

###  II.4 Conventional quantum detailed balance (CQDB)

As mentioned in the introduction, a variety of definitions of quantum detailed
balance have been formulated for Lindblad Markovian master equations Duvenhage
and Snyman (2018); Agarwal (1973); Cipriani (1997); Kossakowski  _et al._
(1977). The simplest and best-known definition of quantum detailed balance
involves the same constraint on correlation functions that exists in the
classical setting Agarwal (1973); Carmichael and Walls (1976). With this
definition (which we call conventional quantum detailed balance (CQDB)),
detailed balance requires:

$\displaystyle\langle\hat{X}(t)\hat{Y}(0)\rangle=\langle\tilde{Y}^{\dagger}(t)\tilde{X}^{\dagger}(0)\rangle.$
(11)

Here $\hat{X}$ and $\hat{Y}$ are arbitrary system operators, and tilde is used
to denote the time-reversed version of an operator, i.e.
$\tilde{B}\equiv\hat{T}\hat{B}\hat{T}^{-1}$, where $\hat{T}$ is the anti-
unitary time-reversal operator. Note that on the RHS of Eq. (11), we have
$\tilde{Y}^{\dagger}(t)\equiv\mathcal{E}_{t}\left[\tilde{Y}^{\dagger}\right]$.

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 Carmichael and Walls (1976). 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 Denisov  _et
al._ (2002). Note that CQDB is also referred to as “GNS detailed balance”
Carlen and Maas (2017).

###  II.5 CQDB implies a trivial steady state

As has been noted previously Alicki (1976); Fagnola and Umanita (2007);
Fagnola and Umanità (2010), the CQDB condition is extremely restrictive: it
only holds for systems with steady states that are diagonal in the energy
eigenstate basis (i.e. $[\hat{H},\hat{\rho}_{\rm ss}]=0$)222The Hamiltonian in
a Lindblad master equation is in general not unique. Our condition holds for
the case where jump operators are chosen to be traceless.. 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).

##  III Hidden time-reversal symmetry and generalized detailed balance

###  III.1 Basic formulation

As we have discussed, the simple CQDB condition of Eq. (11) corresponds
directly to microscopic time-reversal symmetry Carmichael and Walls (1976). 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.2). This unusual
symmetry will connect directly to an abstract variant of quantum detailed
balance (so called “SQDB-$\theta$”) studied in the mathematical physics
literature Fagnola and Umanita (2007); Fagnola and Umanità (2010), and which
was recently linked to entanglement Duvenhage and Snyman (2018). 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 $\hat{\rho}_{\rm ss}$, which we write in diagonal form
as

$\displaystyle\hat{\rho}_{\rm ss}$
$\displaystyle=\sum_{n}p_{n}|n\rangle\langle n|.$ (12)

Throughout this work, we will assume that $\hat{\rho}_{\rm 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
$\left|\psi_{T}\right>$ of a doubled version of our system that yields
$\hat{\rho}_{\rm 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
$\left|\psi_{T}\right>$, we first chose an anti-unitary operator $\hat{T}$
which will define our hidden time-reversal symmetry. We then use this choice
to construct $\left|\psi_{T}\right>$ in a manner that mimics the classical
doubled-system state in Eq. (6): we pair each pointer state $\left|n\right>$
in the original system with its time-reversed partner in the auxiliary system.
We thus have

$\displaystyle|\psi_{T}\rangle\equiv\sum_{n}\sqrt{p_{n}}\left|n\right>_{A}\left(\hat{T}\left|n\right>_{B}\right)\equiv\sum_{n}\sqrt{p_{n}}\left|n\right>_{A}\left|\tilde{n}\right>_{B}$
(13)

With this definition, $\left|\psi_{T}\right>$ is invariant under a gauge
change of the $\hat{\rho}_{\rm ss}$ eigenkets $\left|n\right>\rightarrow
e^{i\alpha_{n}}\left|n\right>$. The state $\left|\psi_{T}\right>$ has the form
of a so-called “thermofield double state”; such states have been studied in
many different contexts Takahashi and Umezawa (1996), though usually without
including a time-reversal operation in its definition. We stress that the form
of $|\psi_{T}\rangle$ is contingent on the choice of $\hat{T}$.

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:

$\displaystyle\frac{d}{dt}\hat{\rho}_{AB}(t)=\left[\mathcal{L}\otimes\mathds{1}_{B}\right]\hat{\rho}_{AB}(t)$
(14)

where $\mathcal{L}$ is our original single-system Liouvillian, and
$\mathds{1}_{B}$ denotes the unit superoperator acting on subsystem $B$.
Starting the doubled system at $t=0$ in the pure state
$\hat{\rho}_{AB}(0)=\left|\psi_{T}\right>\left<\psi_{T}\right|$, this dynamics
leaves the reduced density matrices of each subsystem invariant. It does not
however leave the full joint state $\hat{\rho}_{AB}$ time invariant. The
result is that single-system expectation values are time-independent, 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. $\hat{X}_{A}\equiv\hat{X}\otimes\hat{1}$,
$\hat{X}_{B}\equiv\hat{1}\otimes\hat{X}$. Starting in the state
$\left|\psi_{T}\right>$ and evolving as per Eq. (14), hidden TRS holds if all
inter-system correlations obey the following symmetry:

$\mathrm{Tr}\left[\hat{X}_{A}\hat{Y}_{B}\,\hat{\rho}_{AB}(t)\right]=\mathrm{Tr}\left[\hat{Y}_{A}\hat{X}_{B}\,\hat{\rho}_{AB}(t)\right]$
(15)

Stated explicitly, the symmetry requires that inter-system correlations at any
time $t\geq 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 picture. Defining

$C_{XY}^{\text{TFD}}(t)\equiv\begin{cases}\langle\psi_{T}|\hat{X}_{A}(t)\hat{Y}_{B}|\psi_{T}\rangle&t\geq
0\\\
\langle\psi_{T}|\hat{Y}_{A}(-t)\hat{X}_{B}|\psi_{T}\rangle&t<0\end{cases}$
(16)

the condition for having hidden TRS condition then becomes

$\displaystyle C_{XY}^{\rm TFD}(t)$ $\displaystyle=C_{XY}^{\rm TFD}(-t)$ (17)

Note that for a general system that does not have hidden TRS, the definition
in Eq. (16) implies that not only can $C_{XY}^{\rm TFD}(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-unitary $\hat{T}$. As
we will see, there exist physical systems where hidden-TRS is in some sense
degenerate: there are a whole family of distinct operators $\hat{T}$ 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 operator $\hat{T}$ used to define
$\left|\psi_{T}\right>$ 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:

$\hat{T}\hat{\rho}_{\rm ss}\hat{T}^{-1}=\hat{\rho}_{\rm ss}$ (18)

This follows directly by assuming hidden TRS (i.e. Eq. (15)), and taking the
choice $\hat{Y}=1$. The only way the remaining condition can hold for all
$\hat{X}$ 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 $\hat{\rho}_{\rm ss}$, c.f. Eq. (12). Using Eq. (16), we have

$\displaystyle C_{XY}^{\text{TFD}}(t>0)$
$\displaystyle=\sum_{nm}\sqrt{p_{n}p_{m}}\langle
n|\hat{X}(t)|m\rangle\langle\tilde{n}|\hat{Y}|\tilde{m}\rangle$
$\displaystyle=C_{XY}^{\text{cl}}(t)+C_{XY}^{\text{en}}(t),$ (19)

where

$\displaystyle C_{XY}^{\text{cl}}(t)$ $\displaystyle\equiv\sum_{n}p_{n}\langle
n|\hat{X}(t)|n\rangle\langle\tilde{n}|\hat{Y}|\tilde{n}\rangle$ (20)
$\displaystyle C_{XY}^{\text{en}}(t)$ $\displaystyle\equiv\sum_{n\neq
m}\sqrt{p_{n}p_{m}}\langle
n|\hat{X}(t)|m\rangle\langle\tilde{n}|\hat{Y}|\tilde{m}\rangle.$ (21)

For times $t<0$, Eqs. (19)-(21) are defined analogously to $C_{XY}^{\rm
TFD}(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.

$\displaystyle\hat{\rho}_{AB}(0)=\sum_{n}p_{n}|n,\tilde{n}\rangle\langle
n,\tilde{n}|.$ (22)

In this case, our entire inter-system correlation functions would be given by
Eq. (20):

$\displaystyle\mathrm{Tr}\left[\hat{X}_{A}(t)\hat{Y}_{B}\,\hat{\rho}_{AB}\right]$
$\displaystyle=\sum_{n}p_{n}\langle
n|\hat{X}(t)|n\rangle\langle\tilde{n}|\hat{Y}|\tilde{n}\rangle$
$\displaystyle=C_{XY}^{\text{cl}}(t).$ (23)

We can thus view $C_{XY}^{\text{en}}(t)$ in Eq. (21) as the extra contribution
to the TFD correlator that stems from having quantum entanglement (as opposed
to purely classical correlations).

###  III.2 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 Carmichael and Walls (1976)). 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.5 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 $\hat{\rho}_{\rm ss}$, the so-called modular Hamiltonian:

$\displaystyle\hat{H}_{\rho}\equiv-\log\hat{\rho}_{\rm ss}.$ (24)

Any system possessing CQDB has an extremely tight constraint on its dynamics
Bratteli and Haagerup (1978): dynamical evolution generated by the full
Liouvillian $\mathcal{L}$ must commute with evolution generated by
$\hat{H}_{\rho}$ (i.e. by the unitary $e^{-i\hat{H}_{\rho}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. Fagnola and Umanita (2007); Fagnola and Umanità (2010). In Appendix E,
we give a non-technical proof of this result: for any system with modular
symmetry, CQDB and hidden TRS are equivalent Fagnola and Umanita (2007), 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 $\hat{c}_{l}\rightarrow\lambda\hat{c}_{l}$ in the master
equation, and let $\lambda\rightarrow 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 $\hat{H}$, i.e.

$\displaystyle\mathcal{L}[\hat{\rho}]\sim-i[\hat{H},\hat{\rho}].$ (25)

Further, in this limit the steady state $\hat{\rho}_{\rm ss}$ must commute
with $\hat{H}$ 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.

###  III.3 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. Zhu and Monroe (2020)). Formally, one could
define this correlator as a single-system quantity involving a state-dependent
observable. We first introduce a superoperator $\mathcal{J}$ which acts on
single-system operators:

$\displaystyle\mathcal{J}[\hat{X}]_{A}|\psi_{T}\rangle$
$\displaystyle\equiv\hat{X}_{B}|\psi_{T}\rangle$ (26)

$\mathcal{J}$ is well-defined and unique when $\hat{\rho}_{ss}$ is full-rank,
and is given by the explicit formula (see App. D)

$\displaystyle\mathcal{J}[\hat{X}]$ $\displaystyle=\hat{\rho}_{\rm
ss}^{1/2}\tilde{X}^{\dagger}\hat{\rho}_{\rm ss}^{-1/2}.$ (27)

The TFD correlator in Eq. (16) can then be written as a single system
correlator, e.g. for $t>0$:

$\displaystyle
C^{\text{TFD}}_{XY}(t)=\langle\hat{X}(t)\mathcal{J}[\hat{Y}]\rangle=\textrm{tr}\left(\hat{X}(t)\mathcal{J}[\hat{Y}]\hat{\rho}_{\rm
ss}\right)$ (28)

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:

$\displaystyle\textrm{tr}\left(\hat{X}(t)\mathcal{J}[\hat{Y}]\hat{\rho}_{\rm
ss}\right)=\textrm{tr}\left(\hat{Y}(t)\mathcal{J}[\hat{X}]\hat{\rho}_{\rm
ss}\right)$ (29)

The above condition is formally equivalent to one of the many generalized
quantum detailed balance conditions discussed in Ref. Fagnola and Umanita
(2007); Fagnola and Umanità (2010) (so-called SQDB-$\theta$). It is also
referred to as “KMS detailed balance” Carlen and Maas (2017).

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 $\hat{\rho}_{\rm ss}$, i.e. it is nonlinear in $\hat{\rho}_{\rm
ss}$. This is reminiscent of the problem of measuring Renyi entropies Ekert
_et al._ (2002); Islam  _et al._ (2015); Kaufman  _et al._ (2016). Remarkably,
all hope is not lost. As we will show in Sec. V, systems with hidden TRS are
guaranteed to have their effective Hamiltonian $\hat{H}_{\rm eff}$ and jump
operators $\hat{c}_{k}$ transform very simply under the action of
$\mathcal{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 time-symmetry.

###  III.4 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

$\hat{H}=\Delta\hat{\sigma}_{z}+\frac{\Omega}{2}\hat{\sigma}_{x},\,\,\,\hat{c}_{1}=\sqrt{\kappa}\hat{\sigma}_{-}.$
(30)

For simplicity, we consider the resonant-driving case $\Delta=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 operator $\hat{T}$; as discussed in App. C, the
requirement that the steady state be TRS invariant (c.f. Eq. (18)) constrains
$\hat{T}$ to the form:

$\hat{T}=\left[\frac{\sin(\psi/2)}{\sqrt{4b^{2}+1}}\left(\hat{1}-2ib\hat{\sigma}_{x}\right)+i\cos(\psi/2)\hat{\sigma}_{z}\right]\hat{K}_{z},$
(31)

where $b=\Omega/\kappa$, $\hat{K}_{z}$ is the complex conjugation operator
acting in the $\hat{\sigma}_{z}$ basis, and $\psi$ is at this stage an
arbitrary real parameter.

To determine whether our system has hidden-TRS, we must find a $\hat{T}$ 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_{yz}^{\rm TFD}(t)$, and
we leave the other two correlation functions to App. G.

Using Eq. (19) we can decompose this into the classical correlation
$C_{yz}^{\rm cl}(t)$ and the entanglement correction $C_{yz}^{\rm en}(t)$ as:
$C_{yz}^{\rm TFD}(t)=C_{yz}^{\rm cl}(t)+C_{yz}^{\rm en}(t)$. The classical
correlation is independent of $\psi$, and its time asymmetry is nonzero
irrespective of how $\hat{T}$ is chosen:

$C_{yz}^{\rm cl}(t)-C_{yz}^{\rm
cl}(-t)=-\frac{32b^{5}\sin\left(\frac{\alpha}{4}\kappa
t\right)e^{-\frac{3}{4}\kappa t}}{\alpha(2b^{2}+1)(4b^{2}+1)}.$ (32)

Here, $\alpha=\sqrt{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 $\psi$:

$C_{yz}^{\rm en}(t)-C_{yz}^{\rm
en}(-t)=-\frac{32b^{5}\sin\left(\frac{\alpha}{4}\kappa
t\right)e^{-\frac{3}{4}\kappa t}}{\alpha(2b^{2}+1)(4b^{2}+1)}\cos\psi.$ (33)

Comparing Eqs. (32) and (33), we see that for the TRS with $\psi=\pi$, the
entanglement correction to $C_{yz}^{\rm TFD}(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 $\psi=\pi$.
For reference, the single qubit correlation function $C_{yz}(t)$ for $b=1$ at
$\psi=\pi$ 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 operator $\hat{T}_{h}$

$\hat{T}_{\rm
h}=\frac{1}{\sqrt{4b^{2}+1}}(\hat{1}-2ib\hat{\sigma}_{x})\hat{K}_{z}.$ (34)

In App. G we confirm that the remaining two correlation functions are time
symmetric for this TRS.

It is interesting to consider the form of $\hat{T}_{h}$ in various limits. For
weak Rabi driving (i.e. $b\to 0$), Eq. (34) reduces to
$\hat{T}_{h}=\hat{K}_{z}$. In this limit the qubit system in fact satisfies
SDQB with $\hat{T}=\hat{K}_{z}$ (i.e. all correlation functions have standard
Onsager time-symmetry). In the strong drive limit $b\gg 1$,
$\hat{T}_{h}\rightarrow-i\hat{\sigma}_{x}\hat{K}_{z}$. Up to a phase, this
just complex conjugation in the $\hat{\sigma}_{y}$ basis. To make sense of
this, consider the steady state Eq. (117) in this limit. To first order in
small $b^{-1}\ll 1$, the steady state reduces to $\hat{\rho}_{\rm
ss}=\frac{\hat{1}}{2}+\frac{b^{-1}}{2}\hat{\sigma}_{y}$. The form of the
hidden-TRS operator in this limit thus directly reflects the eigenvectors of
$\hat{\rho}_{\rm ss}$. Furthermore, one can show that for any $b$, the hidden
TRS $\hat{T}_{\rm h}$ corresponds to complex conjugation in the steady state
eigenbasis.

Figure 3: Correlation functions and hidden-TRS in a driven qubit. Stationary,
connected $\langle\sigma_{y}(t)\sigma_{z}(0)\rangle$ correlation functions for
the dissipative Rabi-driven qubit system in Eq. (30), for a drive $\Omega$
equal to the decay rate $\kappa$. Blue: the standard single-system correlation
function $C_{yz}(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 operator $\hat{T}$ 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.
Stannigel  _et al._ (2012) and extended in Ref. Roberts and Clerk (2020).

###  IV.1 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 $\hat{\rho}_{\rm ss}$ and choice of
anti-unitary $\hat{T}$. We will take $\hat{\rho}_{\rm ss}$ to be full rank in
what follows, and will consider intra-system correlations in this TFD state as
defining a bilinear form:

$\displaystyle\langle\langle\hat{X},\hat{Y}\rangle\rangle_{\mathcal{T}}\equiv\langle\psi_{T}|\hat{X}_{A}\hat{Y}_{B}|\psi_{T}\rangle.$
(35)

where $\hat{X}$, $\hat{Y}$ denote arbitrary single-system operators. This
bilinear form can then be used to define the dual $\mathcal{E}^{*}$ of any
given single-system superoperator $\mathcal{E}$ via

$\displaystyle\langle\langle\mathcal{E}[\hat{X}],\hat{Y}\rangle\rangle_{\mathcal{T}}\equiv\langle\langle\hat{X},\mathcal{E}^{*}[\hat{Y}]\rangle\rangle_{\mathcal{T}}.$
(36)

Of particular interest is the case where $\mathcal{E}$ is the adjoint
evolution operator $\mathcal{E}_{t}=\exp[\bar{\mathcal{L}}t]$ defined in Eq.
(10) in terms of the adjoint Liouvillian $\bar{\mathcal{L}}$ (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
$\left(\mathcal{E}_{t}\right)^{*}$ 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
$\left(\mathcal{E}_{t}\right)^{*}=\exp\left[\bar{\mathcal{L}}^{*}t\right]$.
Thus, for a given subsystem-$A$ dynamics $\bar{\mathcal{L}}$, we have a
corresponding “mirrored” dynamics $\bar{\mathcal{L}}^{*}$ for subsystem-$B$,
defined by the constraint that it yield identical inter-system correlations,
i.e.

$\displaystyle\left\langle\left(\exp\left[\bar{\mathcal{L}}t\right]\hat{X}\right)_{A}\hat{Y}_{B}\right\rangle_{T}=\left\langle\hat{X}_{A}\left(\exp\left[\bar{\mathcal{L}}^{*}t\right]\hat{Y}\right)_{B}\right\rangle_{T}$
(37)

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 is $\bar{\mathcal{L}}$ is self-dual,

$\bar{\mathcal{L}}^{*}=\bar{\mathcal{L}}.$ (38)

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.
$\langle\langle\hat{X},\hat{Y}\rangle\rangle_{\mathcal{T}}=\langle\langle\hat{Y},\hat{X}\rangle\rangle_{\mathcal{T}}$;
this follows from the $t=0$ limit of Eq. (17). This in turn implies that the
steady state $\hat{\rho}_{\rm ss}$ of the original master equation must be
invariant under our hidden time-reversal operator $\hat{T}$ (i.e. consider the
case $\hat{Y}=\hat{1}$). Combining these two conditions lets us express the
hidden-TRS condition of Eq. (17) as:

$\displaystyle\langle\psi_{T}|\hat{X}_{A}(t)\hat{Y}_{B}(0)|\psi_{T}\rangle=\langle\psi_{T}|\hat{X}_{A}(0)\hat{Y}_{B}(t)|\psi_{T}\rangle.$
(39)

where for either subsystem, $\hat{O}(t)=\mathcal{E}_{t}[\hat{O}]$. 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) 333Note that we are implicitly using the fact that our
bilinear form is non-degenerate (as $\hat{\rho}_{\rm ss}$ is full rank),
something which guarantees the uniqueness of the dual..

###  IV.2 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 Liouvillian
$\bar{\mathcal{L}}$. To do this, we step back and consider a general system
and $\hat{T}$, such that hidden TRS is not necessarily satisfied. We then
explicitly construct the dual Liouvillian $\bar{\mathcal{L}}^{*}$ that
generates the mirrored-system dynamics, by considering each term in Eq. (9).
Our construction will explicitly make use of the exchange superoperator
$\mathcal{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 Liouvillian $\bar{\mathcal{L}}$. This is done using the
following relations, that follow directly from the definition of $\mathcal{J}$
in Eq. (26):

$\displaystyle\langle\langle\hat{X}\hat{O},\hat{Y}\rangle\rangle_{\mathcal{T}}$
$\displaystyle=\langle\langle\hat{X},\hat{Y}\mathcal{J}[\hat{O}]\rangle\rangle_{\mathcal{T}},\leavevmode\nobreak\
$ (40)
$\displaystyle\langle\langle\hat{O}^{\dagger}\hat{X},\hat{Y}\rangle\rangle_{\mathcal{T}}$
$\displaystyle=\langle\langle\hat{X},\mathcal{J}[\hat{O}]^{\dagger}\hat{Y}\rangle\rangle_{\mathcal{T}}.$
(41)

We can thus obtain an explicit expression for the desired dual
$\bar{\mathcal{L}}^{*}$ as

$\displaystyle\bar{\mathcal{L}}^{*}[\hat{A}]$ $\displaystyle\equiv
i(\mathcal{J}[\hat{H}_{\text{eff}}]^{\dagger}\hat{A}-\hat{A}\mathcal{J}[\hat{H}_{\text{eff}}])$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\sum_{l}\mathcal{J}[\hat{c}_{l}]^{\dagger}\hat{A}\mathcal{J}[\hat{c}_{l}].$
(42)

Recall $\hat{H}_{\rm eff}$ is the effective non-Hermitian Hamiltonian, and
$\hat{c}_{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 $\mathcal{J}$.

We now ask what constraints ensue when we insist that the hidden TRS condition
holds, and hence $\bar{\mathcal{L}}^{*}=\bar{\mathcal{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. Parthasarathy (1992)).

Hence, insisting that our system has hidden TRS leads to the following
constraint equations:

$\displaystyle\mathcal{J}[\hat{H}_{\text{eff}}]=\hat{H}_{\text{eff}}+E,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\mathcal{J}[\hat{c}_{l}]=\sum_{k=1}^{M}U_{lk}\hat{c}_{k},\,\,\,\,\,U^{2}=1.$
(43)

$E$ is a real number, and $U_{lk}$ is a $M\times 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 under $\hat{T}$
(c.f. Eq. (18)). This resulting additional symmetry of the TFD state then
implies (via Eq. (26)) that two exchanges yield the identify superoperator:
$\mathcal{J}^{2}=\mathds{1}$. This immediately constrains the unitary matrix
$U_{lk}$ to have purely real eigenvalues 444Note that the condition in Eq.
(43) on the jump operators implies the exchange superoperator acts as a
unitary on the subspace spanned by the Lindblad operators $c_{l}$. Since
$\mathcal{J}$ squares to one, the unitary operation in Eq. (43), which is the
restriction of $\mathcal{J}$ to the subspace spanned by the Lindblad operators
$c_{l}$, also squares to one. It follows that its eigenvalues can only be $1$
or $-1$. .

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 of $\hat{T}$ and the steady
state $\hat{\rho}_{\rm ss}$. We can eliminate the explicit appearance of
$\mathcal{J}$ by using the fact that since $\hat{\rho}_{\rm ss}$ is full rank,
the TFD state is separating (see e.g. Bratteli and Robinson (1997)): if two
subsystem-$A$ operators have the same action on the TFD state, then they must
be identical operators. Stated explicitly:

$\displaystyle\left(\hat{X}_{A}|\psi\rangle=\hat{Y}_{A}|\psi\rangle\right)\Leftrightarrow\left(\hat{X}=\hat{Y}\right).$
(44)

Using this, we can eliminate $\mathcal{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 $\left|\psi_{T}\right>$. Using the definition of the exchange
operator on the resulting state, this gives us an equivalent but more useful
set of constraint equations:

$\displaystyle\hat{H}_{{\rm eff},B}\left|\psi_{T}\right>$
$\displaystyle=(\hat{H}_{{\rm eff},A}+E)\left|\psi_{T}\right>$ (45)
$\displaystyle\hat{c}_{l,B}\left|\psi_{T}\right>$
$\displaystyle=\left(\sum_{k=1}^{M}U_{lk}\hat{c}_{k,A}\right)\left|\psi_{T}\right>$
(46) $\displaystyle U^{2}$ $\displaystyle=UU^{\dagger}=1$ (47)

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 $\hat{H}_{\rm 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
$\left|\psi_{T}\right>$ of the doubled system, some constant $E$ and some
involutory $M\times 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 $\hat{c}_{l}$
operators, as well as the coherent system Hamiltonian $\hat{H}$.

We stress that the above equations are equivalent to those derived in Ref.
Fagnola and Umanita (2007) when considering the abstract “SQDB-$\theta$”
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 $\hat{H}_{\rm eff}$ or the jump operators $\hat{c}_{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.3, 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 Onsager-like time symmetry
(c.f. Eq. (29)).

##  V Hidden TRS as a route to exact solutions

###  V.1 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 Gardiner
(2009)). 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 $\left|\psi_{T}\right>$ 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 $\left|\psi_{T}\right>$,
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 $\hat{\rho}_{\rm ss}$ is obtained by
tracing out system-$B$ from $\left|\psi_{T}\right>$. 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}\times 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}\times 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 $\mathcal{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.

###  V.2 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 Stannigel  _et al._ (2012); Roberts and Clerk (2020). 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 $\hat{\rho}_{\rm ss}$. We also
construct a doubled system as in Sec. III.1 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 by
$\bar{\mathcal{L}}$), we can always construct a corresponding “mirrored”
dynamics on subsystem $B$ (generated by $\bar{\mathcal{L}}^{*}$), such that
either evolution generates the same time-dependent inter-system correlations,
c.f. Eq. (37).

Figure 4:  Hidden TRS and perfect quantum absorbers. (a) A Markovian quantum
system evolves according to a master equation in Lindblad form, with
Hamiltonian $\hat{H}$ and jump operators $\hat{c}_{l}$. (b) A particular
realization of the environment as a collection of unidirectional waveguides.
(c) The dual Lindbladian $\overline{\mathcal{L}}^{*}$ always formally solves
the "perfect absorber" problem for the Lindbladian $\overline{\mathcal{L}}$
depicted in panel (a): when a system described by $\overline{\mathcal{L}}^{*}$
placed downstream, it absorbs all of the output radiation (red squiggly
arrows) emitted by the original $A$ systems. As a result, the two quantum
systems $A$ and $B$ relax into a pure entangled state (which has the general
form of a thermofield double state). In general, the Hamiltonian
$\hat{H}^{\prime}$ and jump operators $\mathcal{J}[\hat{c}_{l}]$ of the $B$
system are extremely complex and difficult to find. (d) If the master equation
in panel (a) has hidden TRS, then it is extremely easy to construct the
Hamiltonian and jump operators of the absorber $B$ system.

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 $\mathcal{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 $\mathcal{J}$
we have:

$\displaystyle\hat{H}_{{\rm eff},A}\left|\psi_{T}\right>$
$\displaystyle=\mathcal{J}[\hat{H}_{{\rm eff}}]_{B}\left|\psi_{T}\right>$ (48)
$\displaystyle\hat{c}_{l,A}\left|\psi_{T}\right>$
$\displaystyle=\mathcal{J}[\hat{c}_{l}]_{B}\left|\psi_{T}\right>$ (49)

where $\hat{H}_{\rm eff}$ is the effective Hamiltonian in our master equation,
and $\hat{c}_{l}$ are the jump operators.

As shown in App. H, these equations can be re-written as:

$\displaystyle\hat{H}_{AB}\left|\psi_{T}\right>=0,\,\,\,\,\,\,\,\,(\mathcal{J}[\hat{c}_{l}]_{B}-\hat{c}_{l,A})|\psi_{T}\rangle$
$\displaystyle=0.$ (50)

Here, the Hermitian Hamiltonian $\hat{H}_{AB}$ describes an interaction
between the two subsystems in our doubled system:

$\displaystyle\hat{H}_{AB}=\hat{H}_{A}-\hat{H}^{\prime}_{B}-\frac{i}{2}\sum_{l=1}^{M}(\hat{c}_{l,A}^{\dagger}\mathcal{J}[\hat{c}_{l}]_{B}-h.c.)$
(51)

with

$\displaystyle\hat{H}^{\prime}\equiv\textrm{Re}\left[\mathcal{J}[\hat{H}_{\text{eff}}]\right],\,\,\,\,\,\hat{H}\equiv\textrm{Re}\left[\hat{H}_{\text{eff}}\right]$
(52)

We denote the Hermitian part of an operator $\hat{A}$ as
$\textrm{Re}[\hat{A}]$. Note that $\hat{H}$ is nothing but the Hermitian
Hamiltonian in our original master equation.

Eq. (50) has an extremely suggestive form: it tells us that
$\left|\psi_{T}\right>$ 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 $\left|\psi_{T}\right>$ is a zero energy
pure-state, steady-state of the cascaded doubled system sketched in Fig. 4c.
In this cascaded system Gardiner (1993); Carmichael (1993), there is an
independent chiral (directional) waveguide associated with each jump operator
$\hat{c}_{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 Gardiner (1993);
Gardiner and Zoller (2000), the full master equation for this system is:

$\displaystyle\partial_{t}\hat{\rho}_{AB}$
$\displaystyle=-i[\hat{H}_{AB},\hat{\rho}_{AB}]+\sum_{l=1}^{M}\mathcal{D}\Big{[}\hat{c}_{l,A}-\mathcal{J}[\hat{c}_{l}]_{B}\Big{]}\hat{\rho}_{AB}.$
(53)

Here $\hat{\rho}_{AB}$ is the density matrix of the doubled system, and
$\mathcal{D}[\hat{z}]\hat{\rho}=\hat{z}\hat{\rho}\hat{z}^{\dagger}-\\{\hat{z}^{\dagger}\hat{z},\hat{\rho}\\}/2$
is the standard Lindblad dissipation operator. One can easily verify that if
$\left|\psi_{T}\right>$ 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 correlation-conserving 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 operator $\hat{T}$, we have a different TFD state, a
different mirrored-dynamics (i.e. Hamiltonian $\hat{H}^{\prime}$, jump
operators $\mathcal{J}[\hat{c}_{l}]$), and hence a different possible coherent
quantum absorber.

###  V.3 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 $\hat{\rho}_{\rm
ss}$, as this is needed to construct the exchange superoperator $\mathcal{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

$\mathcal{J}[\hat{c}_{l}]\rightarrow\left(\hat{d}_{l}=\sum_{m=1}^{M}U_{lm}\hat{c}_{m}\right),\leavevmode\nobreak\
\textrm{Re}\left[\mathcal{J}[\hat{H}_{\text{eff}}]\right]\to\hat{H}+E$ (54)

for some involutory unitary $M\times M$ matrix $U$ and real constant $E$, and
the Hamiltonian of the coupled system becomes

$\displaystyle\hat{H}_{AB}=\hat{H}_{A}-\hat{H}_{B}-\frac{i}{2}\sum_{l=1}^{M}\left(\hat{c}_{l,A}^{\dagger}\hat{d}_{l,B}-h.c.\right),$
(55)

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\times 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. Stannigel  _et al._ (2012) 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-$\theta$
quantum detailed balance conditions introduced earlier on purely formal
grounds Fagnola and Umanita (2007); Fagnola and Umanità (2010). 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

Figure 5:  Hidden time-reversal symmetry operations in a parametrically-driven
Kerr cavity. For vanishing single-photon drive $\Lambda_{1}$, the driven-
dissipative cavity model in Eqs. (56)-(57) has two distinct hidden time-
reversal symmetries, corresponding to anti-unitary operators $\hat{T}_{\pm}$.
Here, we plot the Wigner functions of the states
$\hat{T}_{\sigma}\left|\alpha\right>$, where $\left|\alpha\right>$ is a
coherent state (amplitude $\alpha=\sqrt{2}i$, black dots). (a) For weak
nonlinearity, $K=5\times 10^{-4}\kappa_{1}$, and two-photon drive
$\Lambda_{2}=6.25\times 10^{-5}\kappa_{1}$, $\hat{T}_{\pm}$ are simple phase-
space reflections about the axes $\theta=\pm\arg(\Lambda_{2}/i\kappa)$
(indicated by dashed black lines). (b) For strong nonlinearity,
$K=\kappa_{1},\Lambda_{2}=\kappa/8$, hidden-TRS operations become highly non-
Gaussian, as indicated by the presence of significant Wigner negativity in the
final states.

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.4) 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.4). In this section, we turn to a more complex class
of models. These describe a bosonic mode (canonical annihilation operator
$\hat{a}$) 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:

$\displaystyle\hat{H}$
$\displaystyle=\frac{K}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}\hat{a}+\Delta\hat{a}^{\dagger}\hat{a}$
(56) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\bigg{(}\Lambda_{1}\hat{a}^{\dagger}+\frac{\Lambda_{2}}{2}\hat{a}^{\dagger}\hat{a}^{\dagger}+h.c.\bigg{)}$

and with jump operators

$\displaystyle\hat{c}_{1}=\sqrt{\kappa_{1}}\hat{a},\,\,\,\,\,\hat{c}_{2}=\sqrt{\kappa_{2}}\hat{a}^{2}$
(57)

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 Mirrahimi  _et al._ (2014); Leghtas  _et al._ (2015); Wang  _et al._
(2016); Puri  _et al._ (2017); Grimm  _et al._ (2020).

It is well known that the steady state of this class of systems can be found
analytically using quantum-optical phase space methods Drummond and Gardiner
(1980); Wolinsky and Carmichael (1988); Bartolo  _et al._ (2016); more recent
work has shown that these exact solutions can be derived more directly (and
even extended) using the coherent quantum absorber method (CQA) Stannigel  _et
al._ (2012); Roberts and Clerk (2020). 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 Drummond and Gardiner (1980); Wolinsky and
Carmichael (1988).

###  VI.1 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 time-symmetry. An example of such a lack of correlation function
symmetry is shown in Fig. 6. More generally, as discussed in Sec. II.5, CQDB
can only hold if the system’s steady state commutes with $\hat{H}$. This
condition is violated except in the vanishing dissipation limit
$\kappa_{1},\kappa_{2}\rightarrow 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,
$\kappa_{2}=0$. To determine whether our system has hidden-TRS, we consider a
doubled two-cavity system and a two-cavity state $\left|\psi_{T}\right>$. The
question is whether this state could represent a TFD state constructed using
an anti-unitary operator $\hat{T}$ which describes a hidden-TRS (c.f. Eq.
(13)). From Eqs. (45)-(47), such a state must satisfy:

$\displaystyle\hat{a}|\psi_{T}\rangle$
$\displaystyle=u\hat{b}|\psi_{T}\rangle$ (58)
$\displaystyle\hat{H}_{\text{eff},A}|\psi_{T}\rangle$
$\displaystyle=(\hat{H}_{\text{eff},B}+E)|\psi_{T}\rangle$ (59)

for some real energy $E$ and constant $u=\pm 1$. Here (as always)
$\hat{H}_{\rm eff}=\hat{H}-i\kappa_{1}\hat{a}^{\dagger}\hat{a}/2$ is the
effective non-Hermitian Hamiltonian in our master equation. If we can find a
two-cavity state $\left|\psi_{T}\right>$ 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 $\Lambda_{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 $\left|\psi_{T}\right>$. This was explicitly
found and expressed as a confluent hypergeometric function in Ref. Roberts and
Clerk (2020), which demonstrated that this model can be solved using CQA.
Hence, the system has a unique hidden-TRS operator $\hat{T}$ 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, $\Lambda_{1}\to 0$. Because of the single photon loss, the
system still has a unique steady state. However, there are now two distinct
hidden-TRS symmetries $\hat{T}_{\pm}$, each corresponding to distinct TFD
states $|\psi_{T}^{\pm}\rangle$:

$\displaystyle|\psi_{T}^{\pm}\rangle=\sum_{n}\sqrt{p_{n}}|n\rangle_{A}\hat{T}_{\pm}|n\rangle_{B}.$
(60)

We stress that both these states each yield the same $\hat{\rho}_{\rm ss}$
when the auxiliary second cavity is traced out. Formally, the two TFD states
(and corresponding $\hat{T}_{\pm}$) 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 Roberts and Clerk (2020).

We thus have our first example of a physical system with multiple, distinct
hidden TRS symmetries; other examples are listed in Table 1. Using the
explicit forms of the TFD states, we can explicitly compute the action of the
hidden TRS symmetry operators $\hat{T}_{+}$ and $\hat{T}_{-}$. 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\rightarrow 0$, the hidden-TRS
operators $\hat{T}_{\sigma}$ take a simple form. In this case, the two TFD
states limit to simple two-mode squeezed states:

$\displaystyle|\psi_{T}^{\pm}\rangle\underset{K\to
0}{\sim}e^{\frac{\Lambda_{2}}{2\Delta_{\text{eff}}}(\hat{a}^{\dagger}\pm\hat{b}^{\dagger})^{2}}|0,0\rangle,$
(61)

where $\Delta_{\text{eff}}\equiv\Delta+i\kappa_{1}/2$ and $\left|0,0\right>$
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
$\theta=\pm\arg(\Lambda_{2}/\Delta_{\text{eff}})$ in phase space:

$\displaystyle\hat{T}_{\pm}\underset{K\to
0}{\sim}e^{\arg(\pm\Lambda_{2}/\Delta_{\text{eff}})\hat{a}^{\dagger}\hat{a}}\hat{K},$
(62)

where here, $\hat{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 of $\hat{T}_{\pm}$ for both $K=5\times
10^{-4}\kappa_{1}$ (weak nonlinearity) and $K=\kappa_{1}$ (strong
nonlinearity).

###  VI.2 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.3: while hidden-TRS (by
definition) guarantees a symmetry of doubled-system TFD correlation functions,
for certain operators, this directly implies time-symmetry of standard,
single-system correlators. In particular, these special operators are ones
that transform simply under the exchange operator $\mathcal{J}$. By virtue of
Eqs. (45)-(47), the effective Hamiltonian $\hat{H}_{\rm eff}$ and jump
operators $\hat{c}_{k}$ are guaranteed to be such special operators.

As a specific example, consider the following steady-state, single-system
correlation function:

$\displaystyle C_{a^{3},a}(t)$
$\displaystyle\equiv\begin{cases}\langle\hat{a}^{3}(t)\hat{a}(0)\rangle&t\geq
0\\\ \langle\hat{a}(-t)\hat{a}^{3}(0)\rangle&t<0.\end{cases}$ (63)

If we set $\kappa_{2}=0$, Eqs. (45)-(47) ensure that
$\mathcal{J}[\hat{a}]=\pm\hat{a}$. From the definition of the exchange
operator, it follows that $\mathcal{J}[\hat{a}^{m}]=(-1)^{m}\hat{a}^{m}$. As a
result, hidden-TRS guarantees (via Eq. (29)) the above correlator has an
Onsager-like 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 operators
$\hat{X}=(\hat{a}+\hat{a}^{\dagger})/\sqrt{2}$,
$\hat{P}=-i(\hat{a}-\hat{a}^{\dagger})/\sqrt{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
$\mathcal{J}[\hat{a}^{2}]$ is simple, this does not imply that
$\mathcal{J}[\hat{a}^{\dagger}\hat{a}]$ is simple.

We thus have a clear experimental test for confirming the existence of hidden-
TRS in this class of systems.

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)\equiv\langle\hat{a}^{3}(t)\hat{a}\rangle$ (c.f. Eq. (63)) for
a parametrically driven nonlinear cavity with $\Lambda_{2}=K$,
$\kappa_{1}=0.4K$ and $\kappa_{2}=\Lambda_{1}=\Delta=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)$ (where $\hat{X}$, $\hat{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).

##  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.4),
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 driven-
dissipative 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 Dykman and
Krivoglaz (1979); Dykman and Smelyanskii (1988); Marthaler and Dykman (2006);
Dykman (2012); Guo  _et al._ (2013); Guo (2013); Zhang and Dykman (2019).
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.

###  VII.1 Rabi-driven qubit subject to thermal dissipation

A driven-dissipative qubit is a simple example to illustrate the breaking of
hidden TRS due to thermal fluctuations. While this model has hidden-TRS at
zero temperature (c.f. Sec. III.4), this symmetry is broken in the presence of
both coherent driving and thermal fluctuations. The master equation is of the
form Eq. (8) but now with $M=2$ and

$\displaystyle\hat{H}$
$\displaystyle=\Delta\hat{\sigma}_{z}+\frac{\Omega}{2}\hat{\sigma}_{x},$
$\displaystyle\hat{c}_{1}$ $\displaystyle=\sqrt{\kappa(1+\bar{n}_{\rm
th})}\hat{\sigma}_{-},\,\,\,\hat{c}_{2}=\sqrt{\kappa\bar{n}_{\rm
th}}\hat{\sigma}_{+}.$ (64)

where $\bar{n}_{\rm th}$ represents the bath thermal occupancy at the qubit
frequency.

####  VII.1.1 Thermal dissipation with no drive

In the absence of a Rabi drive (i.e. $\Omega=0$), the unique steady state of
our master equation has the thermal equilibrium form:

$\hat{\rho}=\frac{1+\bar{n}_{\rm th}}{1+2\bar{n}_{\rm
th}}\left|g\right>\left<g\right|+\frac{\bar{n}_{\rm th}}{1+2\bar{n}_{\rm
th}}\left|e\right>\left<e\right|.$ (65)

where $\left|g\right>,\left|e\right>$ denote $\hat{\sigma}_{z}$ eigenstates.
This steady state commutes with $\hat{H}$, 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
operators

$\hat{T}=(e^{i\psi}|e\rangle\langle e|+|g\rangle\langle g|)\hat{K}_{z}.$ (66)

For each $\psi$ there is a corresponding $U$ matrix (c.f. Eq. (43))

$U_{\psi}=\begin{pmatrix}0&e^{-i\psi}\\\ e^{i\psi}&0\end{pmatrix}$ (67)

for which the dynamical constraints Eqs. (45)-(47) are satisfied.

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.2.
The qubit-plus-absorber system has the cascaded Hamiltonian Eq. (55) where
$\hat{H}_{A}$ is the qubit Hamiltonian
$\hat{H}=\frac{1}{2}\omega_{0}\hat{\sigma}_{z}$ acting on the physical qubit
$A$ and $\hat{H}_{B}$ is the Hamiltonian acting on the auxiliary qubit $B$
(the absorber). The cascaded system also has the collective jump operators

$\displaystyle\hat{C}_{1}$
$\displaystyle=\hat{c}_{1,A}-e^{-i\psi}\hat{c}_{2,B},$ (68)
$\displaystyle\hat{C}_{2}$
$\displaystyle=\hat{c}_{2,A}-e^{i\psi}\hat{c}_{1,B}.$ (69)

The pure state which is simultaneously dark with respect to $\hat{H}_{AB}$,
$\hat{C}_{1}$, and $\hat{C}_{2}$ is

$|\psi_{0}\rangle=\frac{1}{\sqrt{1+2\bar{n}_{\rm
th}}}\left(\sqrt{1+\bar{n}_{\rm th}}|gg\rangle+e^{i\psi}\sqrt{\bar{n}_{\rm
th}}|ee\rangle\right).$ (70)

After tracing out the absorber system, the single site steady state density
matrix is precisely Eq. (65).

####  VII.1.2 Thermal dissipation with a non-zero drive

We now ask what happens to our thermal qubit when a non-zero drive is added
($\Omega\neq 0$). For simplicity we take $\Delta=0$ (resonant driving), and
define the dimensionless driving
$b^{\prime}\equiv\Omega/\kappa(1+2\bar{n}_{\rm th})$. With this definition,
the steady state $\hat{\rho}_{ss,T}$ of the driven qubit with thermal
dissipation is given by the zero temperature result in Eq. (117) with the
simple substitution $b\rightarrow b^{\prime}$.

Furthermore the eigensystem of the Liouvillian at finite temperature is
obtained from the zero temperature results in Eqs. (120)-(123) by replacements
$\hat{\rho}_{\rm ss}\to\hat{\rho}_{ss,T}$, $b\to b^{\prime}$, and
$\kappa\to\kappa(1+2\bar{n}_{\rm th})$. Finally, the permissible TRS of the
finite temperature system are given by Eq. (31) with $b\to b^{\prime}$.

We consider the TFD correlation function $C_{yz}^{\rm TFD}(t)$ defined in Eq.
(16) for Pauli operators $\hat{\sigma}_{y},\hat{\sigma}_{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:

$\displaystyle C_{yz}^{\rm cl}(t)-C_{yz}^{\rm cl}(-t)=$
$\displaystyle-\frac{8b^{\prime}\sin\left(\frac{\alpha}{4}\kappa
t\right)e^{-\frac{3}{4}\kappa t}}{\alpha(2b^{\prime 2}+1)(4b^{\prime
2}+1)}\left[4b^{\prime 4}+(4b^{\prime 2}+1)\eta_{\rm th}\right]$ (71)

where we have defined

$\eta_{\rm th}=1-\frac{1}{(1+2\bar{n}_{\rm th})^{2}}.$ (72)

In contrast, the $\psi$-dependent entanglement correction is

$\displaystyle C_{yz}^{\rm en}(t)-C_{yz}^{\rm en}(-t)=$
$\displaystyle-\frac{8b^{\prime}\sin\left(\frac{\alpha}{4}\kappa
t\right)e^{-\frac{3}{4}\kappa t}}{\alpha(2b^{\prime 2}+1)(4b^{\prime
2}+1)}\left[2b^{\prime 2}\sqrt{4b^{\prime 4}+(4b^{\prime 2}+1)\eta_{\rm
th}}\right]\cos\psi,$ (73)

which also gains temperature-dependent terms.

At zero temperature the above expressions reduce to Eqs. (32)-(33) so for
$\psi=\pi$, $C_{yz}^{\rm TFD}(t)$ has time symmetry, and our system has a
hidden-TRS. However, as soon as $\bar{n}_{\rm th}$ is non-zero, hidden-TRS is
broken. For non-zero temperature, there is choice of $\psi$ for which the
time-asymmetry of the correlation function. To see this explicitly, suppose we
set the time asymmetry of $C_{yz}^{\rm TFD}(t)$ to zero and attempt to solve
for $\cos\psi$. We obtain

$\cos\psi\overset{?}{=}-\sqrt{1+\eta\frac{4b^{\prime 2}+1}{4b^{\prime 4}}}$
(74)

For any finite temperature $\eta>0$ so the right-hand side has a magnitude
greater than 1, thus there is no solution for $\psi$. 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 with $\bar{n}_{\rm th}$
than the entanglement contribution, breaking the cancellation that occurs at
$\bar{n}_{\rm th}=0$.

Figure 7: Hidden-TRS breaking in a driven qubit. Inset: The doubled-system TFD
$\sigma_{y}$-$\sigma_{z}$ connected correlation function $C_{yz}^{\rm TFD}(t)$
(c.f. Eq. (16)) as a function of time, for a Rabi-driven dissipative qubit for
a non-zero temperature corresponding to $\bar{n}_{\rm th}=1$. The TFD state is
defined by the hidden-TRS operator $\hat{T}$ (c.f. Eq. (34)). Main plot: The
time asymmetry of the TFD correlation function $C_{yz}^{\rm th}(t)-C_{yz}^{\rm
th}(-t)$ versus time $t$ for various temperatures. The values of $\bar{n}_{\rm
th}$ are in order from top to bottom: 0, 0.01, 0.03, 0.1, 0.3, and 1. The
$\bar{n}_{\rm th}=0$ trace is identically zero which reflects the presence
hidden TRS at zero temperature. The onset of asymmetry heralds the breaking of
hidden-TRS with the introduction of thermal fluctuations. All functions are
computed for a resonant Rabi drive with amplitude
$\Omega=\kappa(1+2\bar{n}_{\rm th})$, where $\kappa$ is the loss rate.

###  VII.2 Parametrically-driven nonlinear cavity at finite temperature

Figure 8: Hidden-TRS breaking in a driven nonlinear cavity. Main plot: total
correlation function time-asymmetry $m(T)$ vs. temperature for a
parametrically-driven Kerr resonator, c.f. Eq. (81). Solid green: integrated
asymmetry for the special correlation function $C_{a^{2},H_{\text{eff}}}(t)$
(c.f. Eq. (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);
here $\hat{X}$ and $\hat{P}$ 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)\equiv C_{a^{2},H_{\text{eff}}}(t)$ for $\bar{n}_{\text{th}}=0$ (solid
green), and $\bar{n}_{\text{th}}=0.2$ (dashed green). For all plots we take
$\Lambda_{2}=3K$, $\kappa_{1}=0.01K$, $\Delta=\Lambda_{1}=\kappa_{2}=0$.

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 (two-photon)
drive, and there are only single-photon dissipation processes. We then have a
Lindblad master equation with a Hamiltonian given by Eq. (56) (with
$\Lambda_{1}=\Delta=0$), and dissipators

$\displaystyle\hat{c}_{1}$ $\displaystyle=\sqrt{\kappa(1+\bar{n}_{\rm
th})}\hat{a},\,\,\,\hat{c}_{2}=\sqrt{\kappa\bar{n}_{\rm
th}}\hat{a}^{\dagger}.$ (75)

As is standard in the derivation of quantum optics master equations, the
thermal occupancy $\bar{n}_{\rm th}$ corresponds to a Bose-Einstein factor
evaluated at the bath temperature $T$ and cavity resonance frequency
$\omega_{\rm cav}$: $\bar{n}_{\rm th}=1/\left(\exp[\hbar\omega_{\rm
cav}/k_{B}T]-1\right)$. For simplicity, we ignore two-photon dissipative
processes (i.e. $\kappa_{2}=0$).

At non-zero temperature, hidden TRS is broken unless one or both of the
nonlinearity $K$ or drive $\Lambda_{2}$ are zero. To see this, note that for
hidden-TRS to hold, Eq. (46) requires that for some two-mode state
$\left|\psi_{T}\right>$ and $2\times 2$ involutory unitary matrix $U$ the jump
operators must satisfy

$\displaystyle\bigg{[}\hat{c}_{1,A}-(U_{11}\hat{c}_{1,B}+U_{12}\hat{c}_{2,B})\bigg{]}|\psi_{T}\rangle$
$\displaystyle=0,$ (76)
$\displaystyle\bigg{[}\hat{c}_{2,A}-(U_{21}\hat{c}_{1,B}+U_{22}\hat{c}_{2,B})\bigg{]}|\psi_{T}\rangle$
$\displaystyle=0.$ (77)

The equations imply that $\left|\psi_{T}\right>$ is annihilated by two
independent canonical annihilation operators. As such, this state must be
Gaussian, which in turn implies that $\hat{\rho}_{\rm ss}$ must be Gaussian.
However, this steady state is Gaussian only if one or both of $\Lambda_{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 $\Lambda_{2}$ and
$K$ non-zero, 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$, $\Lambda_{2}=0$, $\bar{n}_{\rm th}>0$). In this case, the
steady state is essentially classical (no Fock-state coherences), and it is
well known that CQDB holds Agarwal (1973). 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

$\displaystyle U=\begin{pmatrix}0&e^{i\theta}\\\ e^{-i\theta}&0\end{pmatrix}.$
(78)

Here, $\theta$ is a real parameter. We thus have a continuous family of
distinct hidden-TRS operators $\hat{T}_{\theta}$ (see Tab. 1), in contrast to
the pair of hidden time-reversal operators $\hat{T}_{\sigma}$ seen for nonzero
parametric driving and zero temperature (i.e.
$K=0,\Lambda_{1}=0,\bar{n}_{\text{th}}=0$).

For each possible hidden-TRS operator $\hat{T}_{\theta}$, we have a
corresponding thermofield double state. These always have the form of
Gaussian, two-mode squeezed states:

$\displaystyle|\psi_{T}(\theta)\rangle\equiv\sum_{n=0}^{\infty}e^{-n(\beta/2+i\theta)}|n,n\rangle.$
(79)

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

$\displaystyle C_{a^{2},H_{\text{eff}}}(t)$
$\displaystyle\equiv\begin{cases}\langle\hat{a}^{2}(t)\hat{H}_{\text{eff}}(0)\rangle&t\geq
0\\\ \langle\hat{H}_{\text{eff}}(-t)\hat{a}^{2}(0)\rangle&t<0\end{cases},$
(80)

where as always, $\hat{H}_{\rm 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 as $\bar{n}_{\rm th}$ is increased.

To see this explicitly, we plot the total time-asymmetry vs. temperature,
which we define as:

$\displaystyle m(T)$
$\displaystyle\equiv\int_{0}^{\infty}|C_{a^{2},H_{\text{eff}}}(t)-C_{a^{2},H_{\text{eff}}}(-t)|dt.$
(81)

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 $\kappa\bar{n}_{\rm th}$
be comparable to the dissipative gap of the zero-temperature system. This
dissipative gap $\Gamma$ (i.e. slow relaxation rate) corresponds to switching
between two coherent states $\left|\pm\alpha\right>$ with
$|\alpha|=\sqrt{(\Lambda_{2}^{2}-\kappa^{2}/4)/K^{2}}$ Puri  _et al._ (2017).
One finds that $\Gamma$ is exponentially small due to the small overlap of
these coherent states Mirrahimi  _et al._ (2014):

$\Gamma=\kappa|\alpha|^{2}e^{-2|\alpha|^{2}}$ (82)

Setting this rate equal to $\kappa\bar{n}_{\rm th}$ using the parameters in
Fig. 8 yields a temperature $k_{B}T/\hbar\omega_{c}\approx 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

$\displaystyle m_{X^{2}P^{2}}(T)$
$\displaystyle\equiv\int_{0}^{\infty}|C_{X^{2},P^{2}}(t)-C_{X^{2},P^{2}}(-t)|dt,$
(83)

in Fig. 8; here, $\hat{X}$, $\hat{P}$ are standard quadrature operators.

##  VIII Hidden TRS and phase-space methods: a quantum-classical
correspondence

Hidden TRS |  Classical TRS
---|---
$|\tilde{\psi}\rangle\equiv\hat{T}|\psi\rangle$ |  $(\tilde{\alpha},\tilde{\beta})\equiv(U\alpha,U\beta)$
$C^{\text{TFD}}_{X,Y}(t)\equiv\langle\hat{X}_{A}(t)\hat{Y}_{B}(0)\rangle$ |  $C^{\mathcal{P}}_{X,Y}(t)\equiv\overline{X(t)\tilde{Y}(0)}$
$\mathcal{J}[\hat{a}]\equiv U\hat{a},$
$\mathcal{J}[\hat{H}_{\text{eff}}]\equiv\hat{H}_{\text{eff}}+E.$ |  Potential conditions
($U$ must be unity)
Table 2: 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.

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

$\displaystyle\partial_{t}\hat{\rho}=-i[\hat{H},\hat{\rho}]+\sum_{j}\kappa_{j}\mathcal{D}[\hat{a}_{j}]$
(84)

Here, $[\hat{a}_{j},\hat{a}_{k}^{\dagger}]=\delta_{jk}$ is a standard set of
independent bosonic modes, $\kappa_{j}$ represents a decay rate for each mode,
and $\hat{H}$ 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 operator $\hat{T}$) coincides with
a classical notion of time-reversal in an expanded phase space, i.e. an
involution of the form $(x,p)\to(\tilde{x},\tilde{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 Drummond and Gardiner (1980). 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 Walls and Milburn (2008) for a
pedagogical introduction). Consider the single-mode case first for simplicity.
We consider a doubled phase space described by complex coordinates
$(\alpha,\beta)$, and chose appropriate integration contours
$\mathcal{C},\mathcal{C}^{\prime}$ for each of these variables. This lets us
express the density matrix as

$\displaystyle\hat{\rho}=\int_{\mathcal{C}}d\alpha\int_{\mathcal{C}^{\prime}}d\beta\frac{|\alpha\rangle\langle\beta^{*}|}{\langle\alpha|\beta^{*}\rangle}P(\alpha,\beta).$
(85)

where $P(\alpha,\beta)$ is the complex-$P$ quasi-distribution function. Using
standard techniques Drummond and Gardiner (1980), one can often convert the
Lindblad master equation for $\hat{\rho}$ into a Fokker-Planck equation for
this function, which is required to be nonsingular on the integration surface
$\mathcal{C}\times\mathcal{C}^{\prime}$ defined by the contours. The resulting
equation has the standard form:

$\displaystyle\partial_{t}P(\alpha,\beta)$
$\displaystyle=\partial_{\mu}\bigg{[}C^{\mu}(\alpha,\beta)P(\alpha,\beta)\bigg{]}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\partial_{\mu}\partial_{\nu}\bigg{[}D^{\mu\nu}(\alpha,\beta)P(\alpha,\beta)\bigg{]}$
(86)

Here, $C^{\mu}(\alpha,\beta)$ represents a generalized drift vector, and
$D^{\mu\nu}(\alpha,\beta)$ a generalized diffusion tensor. The above
derivatives are holomorphic derivatives Drummond and Gardiner (1980), 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 symmetry
$\hat{T}$ 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 to $\hat{T}$. 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 Gardiner
(2009). This directly enables an efficient solution for the steady-state
distribution function $P(\alpha,\beta)$, 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 Drummond and Walls (1980). 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 time-reversal operator corresponding to a hidden
TRS $\hat{T}$, and then show how this correspondence can be broken by
considering higher-order Markovian loss channels.

###  VIII.1 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(\alpha,\beta)$. 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(\alpha,\beta)$ defined on our complex phase space is defined as:

$\displaystyle\overline{X}\equiv\int_{\mathcal{C}}d\alpha\int_{\mathcal{C}^{\prime}}d\beta
P_{ss}(\alpha,\beta)X(\alpha,\beta).$ (87)

We can also define a time-evolved function $X(\alpha,\beta;t)$ defined by the
solution of the dual Fokker-Planck equation

$\displaystyle\partial_{t}X(\alpha,\beta)$
$\displaystyle=-\partial_{\mu}\bigg{[}C^{\mu}(\alpha,\beta)X(\alpha,\beta)\bigg{]}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\partial_{\mu}\partial_{\nu}\bigg{[}D^{\mu\nu}(\alpha,\beta)X(\alpha,\beta)\bigg{]}.$
(88)

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

$\displaystyle(\alpha,\beta)\to\widetilde{(\alpha,\beta)},$ (89)

such that all two-time correlation functions (calculated as defined above) are
time-symmetric:

$\displaystyle
C^{\mathcal{P}}_{X,Y}(t)\equiv\begin{cases}\overline{X(t)Y(0)}&t\geq 0\\\
\overline{\tilde{Y}(-t)\tilde{X}(0)}&t<0\end{cases},$ (90)

Here $X,Y$ are any two holomorphic functions, and the time-reversed functions
are given as

$\displaystyle\tilde{X}(\alpha,\beta)\equiv X(\widetilde{\alpha,\beta})$ (91)

where $\widetilde{(\alpha,\beta)}$ is the time-reversed counterpart to
$(\alpha,\beta)$, 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
$\mathcal{C}\times\mathcal{C}^{\prime}$.

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 Risken and Frank
(1996); Gardiner (2009). 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^{\mu}(\alpha,\beta,t)$ is
defined by rewriting the Fokker-Planck equation as a continuity equation:

$\displaystyle\partial_{t}P(\alpha,\beta;t)\equiv\partial_{\mu}J^{\mu}(\alpha,\beta;t).$
(92)

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
$(\alpha,\beta)\to\widetilde{(\alpha,\beta)}$ 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. Risken and Frank
(1996).)

###  VIII.2 Constructing the classical TRS corresponding to a hidden TRS
$\hat{T}$

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

$\displaystyle\mathcal{J}[\hat{a}]=u\hat{a},$ (93)

where $\mathcal{J}$ is the exchange superoperator as always, and $u=\pm 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:

$\displaystyle\widetilde{(\alpha,\beta)}=(u\alpha,u\beta)$ (94)

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 $\mathcal{J}[\hat{a}]=u\hat{a}$ with $u=\pm 1$, then we have
(see App. F):

$\displaystyle
C_{X,\tilde{Y}}^{\mathcal{P}}(t)=C^{\text{TFD}}_{\hat{X},\hat{Y}}(t)$ (95)

where $X,Y$ are the classical representatives of $\hat{X},\hat{Y}$ in the
complex-$P$ representation. Explicitly, without loss of generality,
$\hat{X},\hat{Y}$ can be written as

$\displaystyle\hat{X}$
$\displaystyle=\sum_{n,m}\chi_{n,m}(\hat{a}^{\dagger})^{n}\hat{a}^{m},$ (96)
$\displaystyle\hat{Y}$
$\displaystyle=\sum_{n,m}\lambda_{n,m}(\hat{a}^{\dagger})^{n}\hat{a}^{m}.$
(97)

In terms of the normally-ordered expressions above, the classical
representatives $X,Y$ have the following form:

$\displaystyle X(\alpha,\beta)\equiv\sum_{n,m}\chi_{n,m}\beta^{n}\alpha^{m},$
(98) $\displaystyle
Y(\alpha,\beta)\equiv\sum_{n,m}\lambda_{n,m}\beta^{n}\alpha^{m}.$ (99)

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.

###  VIII.3 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

$\displaystyle\frac{d}{dt}\hat{\rho}=-i[\hat{H},\hat{\rho}]+\kappa\mathcal{D}[\hat{a}^{2}]\hat{\rho}.$
(100)

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:
$\Lambda_{1}=\kappa_{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:

$\displaystyle|\psi_{T}(\gamma,\nu)\rangle\equiv\gamma|\psi_{T}^{+}\rangle+\nu|\psi_{T}^{-}\rangle.$
(101)

where the individual terms $\left|\psi_{T}^{\pm}\right>$ 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):

$\displaystyle\hat{\rho}_{\rm ss}^{(0)}$
$\displaystyle=\text{Tr}_{B}[|\psi^{+}_{T}\rangle\langle\psi_{T}^{+}|]=\text{Tr}_{B}[|\psi^{-}_{T}\rangle\langle\psi_{T}^{-}|].$
(102)

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):

$\displaystyle\hat{\rho}_{\rm ss}^{(0)}$
$\displaystyle=\int_{\mathcal{C}}d\alpha\int_{\mathcal{C}^{\prime}}d\beta\frac{|\alpha\rangle\langle\beta^{*}|}{\langle\alpha|\beta^{*}\rangle}P_{ss}(\alpha,\beta)$
(103)
$\displaystyle=\int_{\mathcal{C}}d\alpha\int_{\mathcal{C}^{\prime}}d\beta\frac{|\alpha\rangle\langle\beta^{*}|}{\langle\alpha|\beta^{*}\rangle}P_{ss}(-\alpha,-\beta)$
(104)

In light of the above observation, one might interpret the more exotic
thermofield doubled state $|\psi(\gamma,\nu)\rangle$ 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) Bartolo  _et al._ (2016); Wang  _et al._ (2016).
Tracing out the auxiliary cavity $B$ for arbitrary parameters $\gamma,\nu$
reveals a family of quantum steady states:

$\displaystyle\hat{\rho}_{\rm ss}(\gamma,\nu)$
$\displaystyle\equiv\text{Tr}_{B}[|\psi_{T}(\gamma,\nu)\rangle\langle\psi_{T}(\gamma,\nu)|].$
(105)

While any $\hat{T}$ and TFD must yield an exchange superoperator $\mathcal{J}$
that acts simply on $\hat{a}^{2}$ (via Eq. (46)), the action on $\hat{a}$ need
not be simple. In fact, we only get a simple action when $\nu=0$ ($\gamma=0$),
in which case $\mathcal{J}[\hat{a}]=\hat{a}$ (
$\mathcal{J}[\hat{a}]=-\hat{a}$). For the more general case, the identity Eq.
(93) is broken. The more complicated nature of the hidden TRS operator and the
corresponding $\mathcal{J}$ implies that the steady states corresponding to
$|\psi_{T}(\gamma,\nu)\rangle$ 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 Agarwal (1973); crucially, hidden-TRS can exist in systems whose
steady states have non-zero 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
time-symmetry. 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 Stannigel  _et al._ (2012); Roberts and
Clerk (2020), as well as phase space methods based on the complex-$P$ function
Drummond and Gardiner (1980).

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 Kamenev (2011); Sieberer  _et al._ (2015,
2016)). This could yield further insights, and also perhaps enable an
extension of these ideas into non-Markovian regimes.

## Acknowledgements

This work is supported by the Air Force Office of Scientific Research MURI
program under Grant No. FA9550-19-1-0399, and by the Army Research Office
under Grant No. W911NF-19-1-0380. AC also acknowledges support 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:

$\displaystyle\overline{X_{A}(t)Y_{B}(0)}=\overline{Y_{A}(t)X_{B}(0)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall X,Y.$
(106)

The reason why the above condition is equivalent to the standard definition of
detailed balance is that the doubled-system correlation function
$\overline{X_{A}(t)Y_{B}(0)}$ is actually a single-system correlation function
in disguise:

$\displaystyle\overline{X_{A}(t)Y_{B}(0)}$
$\displaystyle=\sum_{n,m}\bar{p}(n)\delta_{m,\tilde{n}}X(t,n)Y(0,m)$ (107)
$\displaystyle=\sum_{n}\bar{p}(n)X(t,n)Y(0,\tilde{n})=\overline{X(t)\tilde{Y}(0)},$
(108)

Therefore, the doubled-system correlation function
$\overline{X_{A}(t)Y_{B}(0)}$ is time-symmetric for all random variables $X,Y$
if and only if

$\displaystyle\overline{X(t)\tilde{Y}(0)}=\overline{Y(t)\tilde{X}(0)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall X,Y.$
(109)

Now, making the replacement $Y\to\tilde{Y}$ in this single-site symmetry
condition yields the definition of classical detailed balance used in the main
text. Therefore, Eq. (109) is equivalent to classical detailed balance:

$\displaystyle\overline{X(t)Y(0)}=\overline{\tilde{Y}(t)\tilde{X}(0)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \forall X,Y.$
(110)

Note that we have implicitly used the fact that $Y\to\tilde{Y}$ is a bijection
of the algebra of random variables. In summary, we have shown that the doubled
definition Eq. (106) of classical detailed balance is completely equivalent to
the standard definition Eq. (110).

##  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 Kossakowski  _et al._ (1977);
Bratteli and Haagerup (1978); Alicki (1976). However, here we will assume an
intermediate result, namely that CQDB implies modular symmetry Bratteli and
Haagerup (1978), that is, a symmetry of the driven-dissipative dynamics with
respect to the unitary dynamics generated by the modular Hamiltonian

$\displaystyle\hat{H}_{\rho}\equiv-\log\hat{\rho}_{\rm ss}.$ (111)

The reason for taking this symmetry-based perspective is that it informs most
of the central results in the theory of quantum detailed balance Fagnola and
Umanita (2007); Fagnola and Umanità (2010).

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 Lindblad (1976):

$\displaystyle
i\hat{H}_{\text{eff}}=\overline{\bar{\mathcal{L}}[U^{\dagger}]\cdot U}$ (112)

where here, $U$ is a Haar-random unitary, and $\bar{\mathcal{L}}$ is the
Heisenberg-picture Lindbladian, which generates time-evolution of observables.
We now observe how the effective Hamiltonian evolves under the modular
(dynamical) group $\hat{O}(t)\equiv\hat{\rho}_{\rm
ss}^{it}\hat{O}\hat{\rho}_{\rm ss}^{-it}$:

$\displaystyle i\hat{H}_{\text{eff}}(t)=\overline{\hat{\rho}_{\rm
ss}^{it}\bar{\mathcal{L}}[\hat{U}^{\dagger}]\cdot\hat{U}\hat{\rho}_{\rm
ss}^{-it}}.$ (113)

The above identity holds for any Lindbladian, and thus contains no hidden
assumptions about the system in question. Now suppose, however, that
$\bar{\mathcal{L}}$ satisfies CQDB, and thus has modular symmetry. Then, in
particular, we have

$\displaystyle\hat{\rho}_{\rm
ss}^{it}\bar{\mathcal{L}}[\hat{U}^{\dagger}]=\bar{\mathcal{L}}[\hat{\rho}_{\rm
ss}^{it}\hat{U}^{\dagger}\hat{\rho}_{\rm ss}^{-it}]\hat{\rho}_{\rm ss}^{it}.$
(114)

Substituting the above identity into the Haar-random average defining the
effective non-hermitian Hamiltonian yields a simple result -
$\hat{H}_{\text{eff}}$ is time-independent with respect to the unitary
dynamics generated by the modular Hamiltonian:

$\displaystyle
i\hat{H}_{\text{eff}}(t)=\overline{\bar{\mathcal{L}}[\hat{U}^{\dagger}(t)]\cdot\hat{U}(t)}=\overline{\bar{\mathcal{L}}[\hat{U}^{\dagger}]\cdot\hat{U}}=i\hat{H}_{\text{eff}}(0).$
(115)

Therefore, the effective Hamiltonian commutes with the modular Hamiltonian,
and thus the effective Hamiltonian commutes with the steady state itself:
$[\hat{H}_{\text{eff}},\hat{\rho}_{\rm ss}]=0$. By taking the anti-hermitian
part of the commutator, we immediately have that $\hat{H}$ commutes with
$\hat{\rho}_{\rm ss}$.

##  Appendix C Example of broken CQDB: dissipative Rabi-driven qubit

Figure 9: The Rabi-driven qubit violates conventional quantum detailed
balance. Inset: The correlation function $C_{yz}(t)$ is shown as a function of
time for TRS $\psi=0$. Although time symmetry is not explicitly ruled out for
$\psi=0$, the lack of time symmetry is clear. We show only the connected part
which decays to zero for $|t|\gg\kappa$. Main plot: The time asymmetry of the
symmetrized correlation function ${\rm Re}(C_{yz}(t)-C_{yz}(-t))$ as a
function of time $t/\kappa$ for different permissible TRS. The correlation
function is complex in general for $\psi\neq 0,\pi$; however, the time
asymmetry is manifest even in the real part alone. The values of $\psi$ in
order from top to bottom are: $\pi$, $5\pi/6$, $3\pi/4$, $5\pi/8$, $\pi/2$,
and $0$. All functions are computed for Rabi drive strength $b=1$.

To make the ideas of sections II.4 and II.5 more concrete, we consider simple
but ubiquitous system which does not satisfy CQDB, does have hidden TRS.

###  C.1 Violation of CQDB via correlation function asymmetry

Our example is a qubit (with Pauli operators $\hat{\sigma}_{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 has the form of Eq. (8) with $M=1$ and

$\hat{H}=\Delta\hat{\sigma}_{z}+\frac{\Omega}{2}\hat{\sigma}_{x},\,\,\,\hat{c}_{1}=\sqrt{\kappa}\hat{\sigma}_{-}.$
(116)

Here $\Delta$ is the detuning of the drive from the qubit splitting frequency,
$\Omega$ is the Rabi frequency, and $\kappa$ 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. Gardiner and Zoller (2000):

$\displaystyle\hat{\rho}_{\rm ss}=\frac{\hat{1}}{2}$
$\displaystyle-\frac{4\Delta\Omega}{16\Delta^{2}+2\Omega^{2}+\kappa^{2}}\hat{\sigma}_{x}+\frac{\Omega\kappa}{16\Delta^{2}+2\Omega^{2}+\kappa^{2}}\hat{\sigma}_{y}$
$\displaystyle-\frac{16\Delta^{2}+\kappa^{2}}{2(16\Delta^{2}+2\Omega^{2}+\kappa^{2})}\hat{\sigma}_{z}$
(117)

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 operator $\hat{T}$. We will take a general
approach here (and throughout this paper): we do not pre-select the definition
of $\hat{T}$ using on additional knowledge of our system, but rather ask where
there is any possible anti-unitary $\hat{T}$ which would give rise to a
symmetry. Hence to truly rule out CQDB, one must check Eq. (11) for all
permissible choices of $\hat{T}$. We will thus show that CQDB does not hold no
matter what choice is made for $\hat{T}$.

In what follows, for simplicity we assume a resonant drive (i.e. $\Delta=0$),
and introduce the dimensionless Rabi frequency $b\equiv\Omega/\kappa$; CQDB is
violated even for non-resonant drives, see Sec. C.2 below. The first step is
to parameterize all possible TRS operators. Since CQDB only holds if
$\hat{\rho}_{\rm ss}$ is itself invariant under TRS, this constrains the form
of TRS. The only permissible TRS operators are then parameterized by a single
phase $\psi$ and have the form:

$\hat{T}=\left[\frac{\sin(\psi/2)}{\sqrt{4b^{2}+1}}\left(\hat{1}-2ib\hat{\sigma}_{x}\right)+i\cos(\psi/2)\hat{\sigma}_{z}\right]\hat{K}_{z}$
(118)

where $\hat{K}_{z}$ is the complex conjugation operator acting in the
$\hat{\sigma}_{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 $\psi$. As
our main object of study, we introduce the correlation function defined for
positive and negative times:

$C_{yz}(t)=\begin{cases}\langle\hat{\sigma}_{y}(t)\hat{\sigma}_{z}(0)\rangle&t\geq
0\\\ \langle\tilde{\sigma}_{z}(-t)\tilde{\sigma}_{y}(0)\rangle&t<0\end{cases}$
(119)

where $\tilde{\sigma}_{j}=\hat{T}\hat{\sigma}_{j}\hat{T}^{-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
$\lambda_{j}$ denote its eigenvalues and $\hat{r}_{j}$ the corresponding
right-eigenvectors, we find:

$\displaystyle\lambda_{0}=0$ $\displaystyle\qquad\hat{r}_{0}=\hat{\rho}_{\rm
ss}$ (120) $\displaystyle\lambda_{1}=-\frac{\kappa}{2}$
$\displaystyle\qquad\hat{r}_{1}=\hat{\sigma}_{x}$ (121)
$\displaystyle\lambda_{2}=-\frac{\kappa}{4}\left(3+i\alpha\right)$
$\displaystyle\qquad\hat{r}_{2}=\hat{\sigma}_{y}+\frac{(1+i\alpha)}{4b}\hat{\sigma}_{z}$
(122) $\displaystyle\lambda_{3}=-\frac{\kappa}{4}\left(3-i\alpha\right)$
$\displaystyle\qquad\hat{r}_{3}=\frac{(1+i\alpha)}{4b}\hat{\sigma}_{y}+\hat{\sigma}_{z}$
(123)

where $b\equiv\Omega/\kappa$, $\alpha=\sqrt{16b^{2}-1}$ is the dimensionless
damped Rabi frequency. For $b<1/4$ the qubit is overdamped and the frequency
is imaginary: $\alpha\to i\sqrt{1-16b^{2}}$. The essential feature of the
eigensystem is that the $\hat{\sigma}_{x}$ coherence behaves differently from
the $\hat{\sigma}_{y}$ coherence or the classical population
($\hat{\sigma}_{z}$). The $\hat{\sigma}_{x}$ coherence decays exponentially
with rate $(-\lambda_{1})=\kappa/2$. The $\hat{\sigma}_{y}$ coherence and the
classical populations decay with $(-{\rm Re}\,\lambda_{2,3})>\kappa/2$ for any
finite $b$. This implies that
$\langle\hat{\sigma}_{x}(t)\hat{\sigma}_{k}(0)\rangle$ has a time dependence
that is always different from
$\langle\hat{\sigma}_{y}(t)\hat{\sigma}_{k}(0)\rangle$ or
$\langle\hat{\sigma}_{z}(t)\hat{\sigma}_{k}(0)\rangle$ for any
$\hat{\sigma}_{k}$.

With this in mind, we turn to the time-reversed operator $\tilde{\sigma}_{z}$
which, in general is a linear combination of all three Pauli operators. In
particular, the $\hat{\sigma}_{x}$ component is

$\tilde{\sigma}_{z}=-(2b/\sqrt{4b^{2}+1})\sin\psi\hat{\sigma}_{x}+\cdots$
(124)

We can see that for any $\psi\neq 0,\pi$, for which $\sin\psi\neq 0$, the
expression $C_{yz}(t<0)$ will have terms measuring the decay of
$\hat{\sigma}_{x}$ coherence. Thus the time dependence at $t<0$ _must be_
qualitatively different from that at $t>0$. Even at $\psi=0,\pi$ 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 $\psi$. Although the argument breaks down at
$\psi=0,\pi$, we can show definitively that CQDB cannot hold by considering
$C_{xz}(t)$, which is defined analogously with $C_{yz}(t)$. Since at
$\psi=0,\pi$ there are no $\hat{\sigma}_{x}$ terms in $\tilde{\sigma}_{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.

###  C.2 Violation of CQDB for any detuning

In the preceding section we restricted our analysis to the resonantly driven
qubit for which $\Delta=0$ primarily because the diagonalization of the
Liouvillian becomes unwieldy for $\Delta\neq 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
$[\hat{H},\hat{\rho}_{\rm ss}]$ is equivalent to taking the cross product
between the Hamiltonian “vector” $(\Omega/2,0,\Delta)$ and the traceless part
of $\hat{\rho}_{\rm ss}$. Imposing the constraint that the commutator is zero
requires the form of $\hat{\rho}_{\rm ss}$ to be

$\hat{\rho}_{\rm ss}=\frac{1}{2}\left(\hat{1}+\alpha\hat{H}\right)$ (125)

where $\alpha$ is a real constant of proportionality. Crucially, the above
expression is linear in $\Delta$ and $\Omega$, whereas the true steady state
is a quadradic rational function of these parameters _and_ the decay rate
$\kappa$, c.f. Eq. (117). 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.

###  C.3 Permissible TRS

Here we show how the permissible TRS of Eq. (118) are determined from the
steady state Eq. (117). We require only that $\hat{T}$ leaves the steady state
invariant as a necessary condition of CQDB. All possible TRS take the form
$\hat{T}=\hat{V}\hat{K}_{\rho}$ for unitary $\hat{V}$ and complex conjugation
$\hat{K}_{\rho}$ in the steady state eigenbasis such that
$\hat{T}\hat{\rho}_{\rm ss}\hat{T}^{-1}=\hat{\rho}_{\rm ss}$. This condition
implies that the action of $\hat{T}$ on the eigenstates of $\hat{\rho}_{\rm
ss}$ (i.e. pointer states) is restricted to be

$\displaystyle\hat{T}|1\rangle=$ $\displaystyle|1\rangle$ (126)
$\displaystyle\hat{T}|2\rangle=$ $\displaystyle e^{i\psi}|2\rangle$ (127)

up to a global phase. Thus in the steady state eigenbasis the permissible TRS
take the form $\hat{T}=(|1\rangle\langle 1|+e^{i\psi}|2\rangle\langle
2|)\hat{K}_{\rho}$ for any $\psi$, where $\hat{K}_{\rho}$ is complex
conjugation in the steady state eigenbasis and the pointer states are assumed
time-reversal invariant: $\hat{K}_{\rho}|n\rangle=|n\rangle$.

As the final step, we represent $\hat{T}$ in the familiar $\hat{\sigma}_{z}$
basis. Given the change of basis unitary $\hat{U}$ that diagonalizes the
steady state, $\hat{U}\hat{\rho}_{\rm
ss}\hat{U}^{\dagger}=\mathrm{diag}(p_{1},p_{2})$, the permissible TRS are
given in the $\hat{\sigma}_{z}$ basis as

$\hat{T}=\hat{U}^{\dagger}\hat{V}\hat{K}_{\rho}\hat{U}=\hat{U}^{\dagger}\hat{V}\hat{U}^{*}\hat{K}_{z}.$
(128)

For completeness $\hat{U}$ is given in terms of $b\equiv\Omega/\kappa$ and
$s\equiv\sqrt{4b^{2}+1}$ as

$\hat{U}=\frac{1}{\sqrt{2s}}\left(\begin{array}[]{cc}i\sqrt{s-1}&\sqrt{s+1}\\\
i\sqrt{s+1}&-\sqrt{s-1}\end{array}\right).$ (129)

Substituting this into the expression for $\hat{T}$ above recovers Eq. (118)
(which is also Eq. (31) in the main text).

##  Appendix D Explicit construction of exchange superoperator $\mathcal{J}$

By definition, the exchange superoperator is supposed to act on a single site
observable $\hat{O}$ and produce a new single-site observable,
$\mathcal{J}[\hat{O}]$, 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 observable $\hat{O}$:

$\displaystyle\mathcal{J}[\hat{O}]_{B}|\psi_{T}\rangle\equiv\hat{O}_{A}|\psi_{T}\rangle$
(130)

One may interpret the above equation as the vectorization of an operator
equation, according to the rules

$\displaystyle\hat{O}_{A}|\psi_{K}\rangle$
$\displaystyle\to\hat{O}\hat{\rho}^{1/2}_{ss},$ (131)
$\displaystyle\hat{O}_{B}|\psi_{K}\rangle$
$\displaystyle\to\hat{\rho}^{1/2}_{ss}\tilde{O}^{T},$ (132)

where $|\psi_{K}\rangle$ is the thermofield doubled state where the time-