Coupled activity-current fluctuations in open quantum systems under strong symmetries

The study of the statistical and thermodynamical properties of open quantum systems is one of the most fundamental problems nowadays in modern theoretical physics [1, 2]. As the size of technological devices reduces, the understanding of quantum effects becomes crucial for the development of new solutions. Indeed, quantum effects can be engineered to increase the performance of microscopic thermal machines. Examples abound, e.g. quantum refrigerators [3, 4], engines [5–7], batteries [8–10], and switches [11, 12]. These devices can be implemented using different nanotechnologies that are already available, including trapped ions [13, 14], cold atoms [15, 16], and molecular spins [17, 18]. In most situations of interest, these systems are externally driven and operate under out-of-equilibrium conditions, making the study of nonequilibrium quantum thermodynamics crucial for the development of this emerging field. The natural framework to study this set of problems is the theory of open quantum systems [19, 20]. Armed with this toolbox we can study the behavior of an open system in contact with one or several leads that drive it far from equilibrium. Such nonequilibrium systems typically evolve to a steady-state characterized by a finite current (of energy, excitations, etc) and a well-defined stationary activity. In this way, the last years have witnessed the appearance of a number of interesting results concerning quantum nonequilibrium systems, ranging from detailed analyses of Fourier's law in one [21–24] and several dimensions [25–27] to the study of quantum transport in photosynthetic compounds [28–31] as well as in condensed-matter systems [14, 32, 33].

Fluctuations in small quantum systems play a key role as they crucially affect both their function and response to external driving. Moreover, large fluctuations (though rare) can result in drastic changes of behavior in the system of interest, and therefore the investigation of their statistics as well as the typical paths leading to them has been the focus of an intense research effort in recent years, both in the classical [34–48] and quantum [49–66] realms. The mathematical framework to analyze the physics of fluctuations is large deviation theory [51, 67]. The central objects of the theory are the large deviation functions (LDFs) of the different observables of interest, which allow the calculation of probabilities related to typical and not-so-typical fluctuations. The relevant observables are usually the currents (of energy, excitations, particles, etc) characterizing nonequilibrium behavior, which comprise the time-antisymmetric sector, and the dynamical activity featuring the time-symmetric sector. LDFs for the current or the activity are of fundamental importance in nonequilibrium statistical physics as they play a role equivalent to the equilibrium free energy and related potentials, and govern macroscopic behavior out of equilibrium [34–48]. Interesting results along this research line include the existence of dynamical phase transitions (DPTs) in the fluctuations of driven systems [36, 37, 45, 53, 66, 68–94], or the formulation of different fluctuation theorems for currents based on microscopic time reversibility [48, 95–101].

Interestingly, for open quantum systems governed by a Lindblad-type master equation [19, 20], it has been recently shown that the existence of symmetries leads to different invariant subspaces and multiple (degenerate) steady-states [102, 103]. From a quantum jump (trajectory) perspective, the symmetry-induced emergence of multiple steady states is related to an underlying DPT in the current statistics [108, 109] that leads to a dynamical coexistence of different transport channels classified by symmetry. Such DPT manifests as a non-analyticity of the associated current LDF, a phenomenon that has been confirmed in different setups including quantum networks [108] and optical switches [12]; see also [109]. Symmetries, and the associated dynamic phase transitions in current statistics, are very sensitive to external perturbations such as dephasing noise. However, if the noise is sufficiently weak, one can show that the existence of symmetries in the noise-free case can be inferred from the time-dependent current behavior of the (noisy) system of interest [110]. Moreover, symmetries can be also manipulated by the presence of magnetic fields, resulting in a detailed control of nonequilibrium currents [111].

Up to now, research has been focused on understanding the effects of symmetries on the statistics of a single relevant observable, typically the current. A natural question hence concerns the effect of symmetries in the statistics of the dynamical activity, a time-symmetric observable of direct experimental relevance which may constraint the range of current fluctuations. Moreover, it is important to understand how symmetries are reflected in the joint large deviation statistics of these two key observables, the current and the activity, which characterize respectively the time-antisymmetric and the time-symmetric sectors of the dynamics. In this paper we address this issue and investigate how symmetries affect both the transport properties and the activity patterns of a particular class of Markovian open quantum system, a three-qubit model under the action of a magnetic field and in contact with a thermal bath. As expected, we find that activity constraints current fluctuations and viceversa, giving rise in particular to an activity-driven current lockdown phase induced by symmetry for this model. Interestingly, the DPT predicted at the steady state and its Gallavotti–Cohen dual are extended into lines of first order DPTs in the current-activity plane, with a nontrivial structure which depends on the transport and activity properties of each of the symmetry phases. The average current and activity of the different symmetry subspaces are analyzed in detail, as well as their dependence on the bath temperature and the external magnetic field. In addition, conditional averages such as the average current for a given value of the activity and the average activity conditioned to a given current are also explored. Finally, we also study the effect of a symmetry-breaking, ergodicity-restoring dephasing channel on the joint activity-current statistics for this model.

We structure the paper as follows. In section 2 we introduce the model of interest in detail, as well as its open dynamics in terms of a Lindblad master equation for the density matrix. Section 3 is devoted to a symmetry analysis of the resulting dynamical equations and the associated degenerate steady states, while section 4 includes quantum Monte Carlo simulations of individual quantum trajectories. These allow us to better understand how the existence of a strong symmetry constraints the system evolution, leading to a remarkable dissipative freezing behavior, as well as to study the effect of a symmetry-breaking dephasing channel on the dynamics of qubits and the restoration of ergodicity. In section 5 we introduce the counting statistics or large deviation approach to investigate the thermodynamics of quantum trajectories biased over both the current and the activity. Section 6 is then devoted to analyze the spectral consequences of a strong symmetry in the system, and how this leads to DPTs in univariate LDFs, which affect both current and activity statistics. In section 7 we investigate the joint activity-current statistics and the symmetry-induced lines of DPTs appearing in the current-activity plane. Finally, section 8 presents our conclusions and outlook for future investigation.

Zoom In Zoom Out Reset image size

Our system consists in three spins with an XX interaction term, forming an equilateral triangle as displayed in figure 1 . Together with the XX coupling, there is a magnetic field acting along the Z direction on each spin. This kind of spin models have been broadly studied in literature, as they are relevant for a number of interesting problems ranging from quantum heat conduction and Fourier's law [22, 27, 112] to noise-assisted transport [29, 31, 113] or phase transitions [108, 114, 115], to mention just a few. Interestingly, spin systems like this one can be experimentally realized using different current technologies, including spins in semi-conductors [116] and ion traps [117].

Zoom In Zoom Out Reset image size

The total dimension of the system is that of three-qubits. We name the pure states Hilbert space as $\mathcal{H}$, having a dimension d = 2³. Mixed states are defined by density matrices ρ that are positive trace-one operators in the space of the bounded operators $\mathcal{B}(\mathcal{H})$. The Hamiltonian controlling the system coherent dynamics is

$$\begin{equation}H=\sum\limits _{i,j=0}^{2}{\sigma }_{i}^{x}{\sigma }_{j}^{x}+{B}_{z}\sum\limits _{i=0}^{2}{\sigma }_{i}^{z},\end{equation} \tag{ 1 }$$

where ${\sigma }_{i}^{x},{\sigma }_{i}^{z}$ are Pauli matrices, and B_z is the strength of the external magnetic field along the Z direction. In addition, the system is driven by the action of a bosonic thermal bath that interacts locally with spin 0. This bath can trigger incoherent jumps in qubit 0 and it is modeled by the Lindblad jump operators [19, 20]

$$\begin{align}{L}_{0}^{+}& =\sqrt{{\Gamma}n}\enspace {\sigma }_{0}^{+}\\ {L}_{0}^{-}& =\sqrt{{\Gamma}(n+1)}\enspace {\sigma }_{0}^{-},\end{align} \tag{ 2 }$$

were Γ is the coupling strength to the bath, ${\sigma }_{0}^{{\pm}}={\sigma }_{0}^{x}{\pm}{\sigma }_{0}^{y}$ are raising/lowering operators acting on spin 0, and $n=1/[{\text{e}}^{(\hslash {B}_{z})/({k}_{\mathrm{B}}T)}-1]$ is the average number of excitations in the bath at the resonance frequency, given by a Bose–Einstein distribution at temperature T. For simplicity, we use units such that ℏ = k_B = 1 throughout the paper.

In addition to the previous ingredients, we also consider a dephasing channel acting locally and independently on all three spins. This allows us to analyze the quantum-to-classical transition at the level of trajectories, as well as the role of dephasing on symmetry-breaking and the restoration of ergodicity (see below). This dephasing channel is modeled by the jump operators

$$\begin{equation}{L}_{i}^{\mathit{\text{D}}}=\sqrt{\gamma }{\sigma }_{i}^{+}{\sigma }_{i}^{-},\quad i=0,1,2\end{equation} \tag{ 3 }$$

with γ being the dephasing strength. The main effect of this channel is to reduce the coherences between different spins without affecting the populations in the site basis. Dephasing in nonequilibrium quantum system is known to have important effects, such as e.g. noise-enhanced transport [29, 31, 113, 118], current suppression [119] and emergence of diffusive heat conduction [22, 25]. In the specific topic of symmetries and invariant subspaces, this kind of channel is often responsible of a noise-induced symmetry breaking that can collapse the multiple invariant subspaces of a Liouvillian into a single, unique steady state [102, 108–110]. The global dynamics of the system is thus given by a Lindblad (or Lindblad–Gorini–Kossakowski–Sudarshan) master equation [19, 20] of the form

$$\begin{equation}\dot {\rho }=-i\left[H,\rho \right]+\sum\limits _{k={\pm}}{L}_{0}^{k}\rho {L}_{0}^{k{\dagger}}-\frac{1}{2}\left\{\rho ,{L}_{0}^{k{\dagger}}{L}_{0}^{k}\right\}+\sum\limits _{i=0}^{2}{L}_{i}^{\mathit{\text{D}}}\rho {L}_{i}^{\mathit{\text{D}}{\dagger}}-\frac{1}{2}\left\{\rho ,{L}_{i}^{\mathit{\text{D}}{\dagger}}{L}_{i}^{\mathit{\text{D}}}\right\}\equiv \mathcal{L}\rho \end{equation} \tag{ 4 }$$

with [A, B] = AB − BA the commutator of two operators A and B, {A, B} = AB + BA the anti-commutator, and $\mathcal{L}$ the Liouvillian superoperator. This superoperator can be expressed in the Fock–Liouville space as a d² × d² complex matrix. In this way, if we have an initial state described by the density matrix ρ(0), its time evolution is formally given by $\rho (t)=\mathrm{exp}(\mathcal{L}t)\rho (0)$. According to Evan's theorem [120], a bounded system like the one discussed here should have at least one steady state, meaning that the superoperator $\mathcal{L}$ defined in equation (4) should have at least one eigenvalue with zero real part [102, 121]. The eigenoperators associated with these null eigenvalues then correspond to the steady-state solutions of the master equation.

Master equations with local couplings have been extensively used in several fields including quantum transport [22, 25] and quantum thermodynamics [104, 105]. Note however that this family of equations does not arise naturally from a microscopic derivation for quantum systems locally coupled to a bath, but they can always be engineered by careful dissipation control. Indeed, in the case of a master equation microscopically derived by tracing over a bath locally coupled to a three-spin system, one would expect to obtain also global jump operators considering the collective modes of the system [106, 107]. This interesting case can also exhibit symmetries leading to a phenomenology similar to the one presented in this paper, but its analysis goes beyond the scope of this work.

In the absence of the dephasing channel (i.e. γ = 0), our system presents an obvious topological symmetry given by the exchange of spins 1 and 2, see figure 1 and equation (4). Using the language of reference [102], when γ = 0 we have a strong symmetry in the system. This is given by a unitary (and Hermitian) operator, ${\pi }_{12}=\frac{1}{2}\left({\sigma }_{1}^{x}{\sigma }_{2}^{x}+{\sigma }_{1}^{y}{\sigma }_{2}^{y}+{\sigma }_{1}^{z}{\sigma }_{2}^{z}+\mathbb{1}\right)$, being $\mathbb{1}$ the identity operator in the three-spins Hilbert space. This operator commutes with all the generators of the system dynamics when there is no dephasing, i.e.

$$\begin{equation}\left[{\pi }_{12},H\right]=\left[{\pi }_{12},{L}_{0}^{{\pm}}\right]=0,\end{equation} \tag{ 5 }$$

thus defining a strong symmetry of the dynamics [102]. As the operator π₁₂ has two different eigenvalues $\left\{-1,+1\right\}$, we can spectrally-decompose the system's Hilbert space in symmetric and antisymmetric subspaces with respect to this symmetry operator, $\mathcal{H}={\mathcal{H}}^{A}\oplus {\mathcal{H}}^{S}$. The symmetric subspace ${\mathcal{H}}^{S}$ has dimension d_S = 6 and it can be spanned by the basis $\left\{{S}_{i},i\in [1,{d}_{S}]\right\}\equiv \left\{\vert {0\rangle }_{0},\vert {1\rangle }_{0}\right\}\otimes \left\{\vert {00\rangle }_{12},\vert {+\rangle }_{12}\equiv \frac{1}{\sqrt{2}}\left(\vert {01\rangle }_{12}+\vert {10\rangle }_{12}\right),\vert {11\rangle }_{12}\right\}$. As expected, the basis of the subsystem formed by spins 1 and 2 in this ${\mathcal{H}}^{S}$ subspace is given by the triplet states due to its symmetric nature. On the other hand, the antisymmetric subspace ${\mathcal{H}}^{A}$ has dimension d_A = 2, and we can define its basis as $\left\{{A}_{i},i\in [1,{d}_{A}]\right\}\equiv \left\{\vert {0\rangle }_{0},\vert {1\rangle }_{0}\right\}\otimes \left\{\vert {-\rangle }_{12}\equiv \frac{1}{\sqrt{2}}\left(\vert {01\rangle }_{12}-\vert {10\rangle }_{12}\right)\right\}$, i.e. in terms of the singlet state for the spin 1 and 2 subsystem. Note that the exchange property of the symmetry operator can be made explicit in the computational basis, i.e. ${\pi }_{12}={\mathbb{1}}_{0}\otimes \left(\vert 10\rangle \langle 01{\vert }_{12}+\vert 01\rangle \langle 10{\vert }_{12}+\vert 11\rangle \langle 11{\vert }_{12}+\vert 00\rangle \langle 00{\vert }_{12}\right)$, with ${\mathbb{1}}_{0}$ the identity operator acting on spin 0. Furthermore, the Hilbert space $\mathcal{B}(\mathcal{H})$ of bounded operators acting on $\mathcal{H}$ can be also decomposed using the symmetry π₁₂ in the form $\mathcal{B}(\mathcal{H})={\mathcal{B}}_{AA}\otimes {\mathcal{B}}_{AS}\otimes {\mathcal{B}}_{SA}\otimes \enspace {\mathcal{B}}_{SS}$, with ${\mathcal{B}}_{\alpha \beta }=\text{span}\left\{\vert {\alpha }_{i}\rangle \langle {\beta }_{j}\vert :i\in [1,{d}_{\alpha }],j\in [1,{d}_{\beta }]\right\}$, $\alpha ,\beta \in \left\{S,A\right\}$, and S_i , A_i represent the elements of the basis of ${\mathcal{H}}^{S}$ and ${\mathcal{H}}^{A}$ respectively. Note that $\mathcal{B}(\mathcal{H})$ is equipped with the Hilbert–Schmidt inner product [19],

$$\begin{equation}\langle \langle \sigma \vert \rho \rangle \rangle =\mathrm{Tr}({\sigma }^{{\dagger}}\rho ),\quad \forall \sigma ,\rho \in \mathcal{B}(\mathcal{H}),\end{equation} \tag{ 6 }$$

where Tr(ω) is the trace of the operator $\omega \in \mathcal{B}(\mathcal{H})$.

Due to the symmetry, the subspaces ${\mathcal{B}}_{\alpha \beta }$ remain invariant under the action of the Liouvillian, meaning that $\mathcal{L}{B}_{\alpha \beta }\subset {B}_{\alpha \beta }$ [102, 103, 108, 109], so $\mathcal{L}$ can be block-decomposed into 2² = 4 invariant subspaces. This can be easily proved by defining the left and right superoperators ${{\Pi}}_{12}^{l,r}$ such that

$$\begin{equation}{{\Pi}}_{12}^{l}\rho ={\pi }_{12}\rho ,\hspace{25.0pt}{{\Pi}}_{12}^{r}\rho =\rho {\pi }_{12}^{{\dagger}},\end{equation} \tag{ 7 }$$

and noticing that (a) the subspaces ${\mathcal{B}}_{\alpha \beta }$ are the joint eigenspaces of both ${{\Pi}}_{12}^{l}$ and ${{\Pi}}_{12}^{r}$, and (b) $\left[{{\Pi}}_{12}^{r},\mathcal{L}\right]=0=\left[{{\Pi}}_{12}^{l},\mathcal{L}\right]$ due to the commutation relations (5). In this way we obtain that, if ${\rho }_{\alpha \beta }\in {\mathcal{B}}_{\alpha \beta }$, then $\mathcal{L}{\rho }_{\alpha \beta }$ is still an eigenoperator of both ${{\Pi}}_{12}^{l,r}$ with the same eigenvalues, so that $\mathcal{L}{\rho }_{\alpha \beta }\in {B}_{\alpha \beta }$. As our system is bounded, we can use now Evan's theorem [120] to prove the existence of at least two fixed points of the dynamics, corresponding to the two null eigenoperators of $\mathcal{L}$ in the diagonal subspaces ${\mathcal{B}}_{\alpha \alpha }$, with α = A or S, which are the only ones to contain unit trace (physical) density matrices. Note that the subspaces ${\mathcal{B}}_{AS}$ and ${\mathcal{B}}_{SA}$ have no physical fixed points as they contain only zero-trace density matrices [111]. In this way, there are two orthogonal steady states; if we initialize the system with a normalized (unit trace) density matrix ${\rho }_{\alpha }(0)\in {\mathcal{B}}_{\alpha \alpha }$, with α = A or S, it will evolve in the long-time limit to a steady state

$$\begin{equation}{\rho }_{\alpha }^{SS}={\mathrm{lim}}_{t\to \infty }\enspace {\text{e}}^{\mathcal{L}t}{\rho }_{\alpha }(0)\in {\mathcal{B}}_{\alpha \alpha },\quad (\alpha =A,S).\end{equation} \tag{ 8 }$$

The existence of two steady states with typically different transport properties can be understood at the quantum trajectory level [51] as a consequence of an underlying DPT of first-order type in the current statistics [108, 109], that leads to a dynamical coexistence of different transport channels classified by symmetry (as observed above). Such dynamical coexistence has been reported in a variety of systems [108, 109], including a three-spin model similar to the one described here [115]. Note also that the steady state degeneracy of open quantum systems with non-abelian symmetries has been recently addressed [122].

When restricted to the antisymmetric subspace, spins 1 and 2 stay frozen into the singlet state, i.e. they fall into a dark (decoherence-free) state and the system dynamics is exclusively due to spin 0, connected to the bath. This is equivalent to effectively removing spins 1 and 2 from the total system. In this way, when restricted to the antisymmetric subspace, the system Hamiltonian can be simply written as

$$\begin{equation}{H}_{A}=-({\mathbb{1}}_{0}\otimes \vert -\rangle \langle -{\vert }_{12})+{B}_{z}{\sigma }_{0}^{z}.\end{equation} \tag{ 9 }$$

For the sake of clarity, in what follows we will not make explicit identity operators (like ${\mathbb{1}}_{0}$ above) acting on subspaces when they can be inferred from the context. In this antisymmetric case, the steady-state density matrix is separable in the partition between qubit 0 and qubits 1, 2. It takes the form ${\rho }_{A}^{SS}={\rho }_{0}^{A}\otimes \vert -\rangle \langle -{\vert }_{12}$, with

$$\begin{equation}{\rho }_{0}^{A}=\frac{1}{1+2n}\left(\begin{matrix}\hfill n\hfill & \hfill \hfill & \hfill 0\hfill \\ \hfill \hfill & \hfill \hfill \\ \hfill 0\hfill & \hfill \hfill & \hfill 1+n\hfill \end{matrix}\right),\end{equation} \tag{ 10 }$$

being the thermal density matrix for one qubit in contact with a bosonic thermal bath with mean number of excitations n. On the other hand, the symmetric subspace interaction is described by the Hamiltonian

$$\begin{align}{H}_{S}& =\sqrt{2}\enspace {\sigma }_{0}^{x}\left(\vert 00\rangle \langle +{\vert }_{12}+\vert +\rangle \langle 00{\vert }_{12}+\vert +\rangle \langle 11{\vert }_{12}+\vert 11\rangle \langle +{\vert }_{12}\right)\\ & \quad +\left(\vert +\rangle \langle +{\vert }_{12}+\vert 00\rangle \langle 11{\vert }_{12}+\vert 11\rangle \langle 00{\vert }_{12}\right)-{B}_{z}\left({\sigma }_{0}^{z}+\vert 00\rangle \langle 00{\vert }_{12}-\vert 11\rangle \langle 11{\vert }_{12}\right).\end{align} \tag{ 11 }$$

Before studying the large deviation statistics of the activity and the current in our three-qubits system (next section), we focus momentarily our interest in understanding how symmetry affects individual quantum jump trajectories. In the presence of a dephasing channel, i.e. when the dephasing rate γ ≠ 0, see equation (3), the exchange symmetry of the three-qubits system is broken as $\left[{L}_{1,2}^{D},{\pi }_{12}\right]\ne 0$. This symmetry-breaking channel restores ergodicity and leads to a unique steady-state independently of the initial state. However, for small dephasing rate the effects of the symmetry subspaces may still be observable in the transient (relaxation) behavior of the system [110]. The purpose of this section is thus to analyze the role of symmetry and its breaking in the stochastic quantum jump dynamics.

To study the dynamical interplay between the different symmetry sectors as a function of dephasing, we define now a symmetry parameter ξ(t) at time t in terms of the projector to the antisymmetric subspace ${P}_{-}={\mathbb{1}}_{0}\otimes \vert -\rangle \langle -\vert$. For a given pure state |ψ(t)⟩ we thus define $\xi (t)\equiv {\left\vert {P}_{-}\vert \psi (t)\rangle \right\vert }^{2}$, while for mixed states we can define its ensemble average ⟨ξ (t)⟩ = Tr(P₋ ρ(t)), with ρ(t) the density matrix at time t. In this way the symmetry parameter ξ(t) captures how individual quantum jump trajectories are projected onto one of the symmetry sectors (in this case the antisymmetric one) as a function of time, quantifying the amount of symmetry selection in the dynamical evolution. We hence generate quantum jump trajectories using a quantum Monte Carlo simulation [123] of the Lindblad Equation (4) for our three-qubits model. Figure 2 shows the symmetry parameter ξ(t) as measured for different quantum trajectories generated in this way, as well as its ensemble average ⟨ξ(t)⟩, as a function of time and for different values of the dephasing rate. The initial state for all quantum trajectories is taken as

$$\begin{equation}\vert {\psi }_{0}\rangle =\frac{1}{2}\left(\vert {0\rangle }_{0}+\vert {1\rangle }_{0}\right)\otimes \left(\vert {00\rangle }_{12}+\vert {-\rangle }_{12}\right),\end{equation} \tag{ 12 }$$

that corresponds to ξ(0) = 1/2, i.e. a fair quantum superposition of both the symmetric and antisymmetric sectors.

As expected, in the no-dephasing limit γ = 0 (left panel in figure 2) both subspaces remain unmixed at the ensemble level so that the value of ⟨ξ(t)⟩ remains constant and equal to ξ(0). Interestingly, however, each individual stochastic trajectory selects randomly one of the symmetry sectors, collapsing in a finite time to the corresponding subspace (making ξ either 0 or 1) and remaining there from that time on. This remarkable behavior, known as dissipative freezing, is a particular instance of a general observation put forward in [124], and implies a breakdown at the individual trajectory level of a conservation law associated to the symmetry operator at the ensemble level. Individual quantum trajectories have equal chances to decay into either symmetry subspaces, restoring the unmixing of the symmetry sectors at the ensemble level and preserving the value of ⟨ξ(t)⟩ = ξ(0).

When dephasing noise is switched on (γ > 0) the dynamics is more complex as there appears mixing between the symmetry subspaces. On one hand we find that, at the ensemble level, the average symmetry parameter ⟨ξ(t)⟩ decreases in time to reach a non-trivial steady-state value below ξ(0) = 1/2, see center and right panels in figure 2, meaning that the noisy dynamics favors trajectories to collapse into the symmetric subspace. This happens because this subspace is of higher dimension (d_S = 6 vs d_A = 2), and hence it is entropically favored in the evolution. Interestingly, the value of γ does not affect appreciably the steady-state value of ⟨ξ⟩, see the red cross in the center and right panels of figure 2. In comparison, the time required for ⟨ξ(t)⟩ to relax to its steady-state value increases as the dephasing rate γ decreases (essentially as 1/γ). The story at the level of individual quantum jump trajectories is rather different. Remarkably, for weak dephasing (γ ≪ 1) the system dynamics is characterized by an intermittent, punctuated evolution, see right panel in figure 2, with long periods of time where the state is trapped in one of the symmetry sectors followed by quick jumps between different sectors. Such intermittent behavior is a dynamical signature of the underlying exchange symmetry [110], and remains observable as far as dephasing noise is weak. As the dephasing strength increases, this intermittent behavior tends to disappear in favor of a rapid succession of jumps between the symmetry sectors, see central panel in figure 2, although the overall picture is similar if time is properly rescaled by the associated relaxation timescale.

In section 3 we have seen how the existence of a symmetry leads to multiple invariant subspaces and degenerate steady states in in our open quantum system governed by a Lindblad master equation (4). The purpose of this section is to set the stage for trajectory statistics in open quantum systems in order to understand in subsequent sections how such symmetry affects the joint statistical properties of two key dynamical observables, the current of excitations and the activity, for the particular case of the three-qubits model of interest in this paper. In order to do so, we use tools from large deviation theory and full-counting statistics [47, 51, 67, 108, 109] to study the thermodynamics of quantum jump trajectories conditioned to a given total current and total activity. The current is the key measure of transport out of equilibrium and a main token of the time-antisymmetric sector of dynamics, while the activity is a direct measure of the open quantum dynamics in the time-symmetric sector, readily accessible in experiments or simulations. Understanding their joints large-deviations statistics thus opens the door to a full characterization of the quantum jump dynamics [51, 67].

For a given quantum jump trajectory (see [109] for a precise definition), the total current Q after a time t corresponds to the net exchange of excitons with the thermal bath. It is defined as

$$\begin{equation}Q\equiv {K}^{+}-{K}^{-},\end{equation} \tag{ 13 }$$

where K⁺ (K⁻) is the total number of quanta absorbed by (emitted from) the system from (to) the bath in a time interval t. On the other hand, the activity is just the total number of quantum jumps during such time interval,

$$\begin{equation}A\equiv {K}^{+}+{K}^{-}.\end{equation} \tag{ 14 }$$

Clearly, these two magnitudes are independent but correlated. In particular, the value of the activity restricts the possible values of the current since −A ⩽ Q ⩽ A, with A ⩾ 0. In addition, constraining the current to take a certain value Q is expected to affect the probability distribution of the activity A, and viceversa. Such interplay between current and activity fluctuations is captured by their joint probability distribution [51, 67, 109]. To describe this joint statistics, we first introduce the reduced density matrix ρ_Q,A (t) which is the projection of the full density matrix to the subspace defined by particular, fixed values for the total current Q and activity A. This reduced density matrix is the solution of a current- and activity-resolved quantum master equation which can be derived from the unraveling of the Liouvillian superoperator $\mathcal{L}$ in equation (4) [109, 125]. The joint probability of observing certain values for Q and A after a time t is thus given by ${P}_{t}(Q,A)=\mathrm{Tr}\left[{\rho }_{Q,A}(t)\right]$, and scales in a large-deviation form

$$\begin{equation}{P}_{t}(Q,A)\asymp \mathrm{exp}[+t\enspace G(q,a)]\end{equation} \tag{ 15 }$$

in the long-time limit, with q = Q/t and a = A/t the time-averaged (intensive) associated quantities. The symbol '≍' means asymptotic logarithmic equality, i.e. ${\mathrm{lim}}_{t\to \infty }\enspace \frac{1}{t}\enspace \mathrm{ln}\enspace {P}_{t}(Q,A)=G(q,a)$. The function G(q, a) ⩽ 0 is the joint LDF for the current and the activity, and contains all the information of the coupled fluctuations of these two central observables.

As usual in statistical physics, working with global constraints (in this case on Q and A) makes the problem cumbersome from a mathematical point of view (think for instance on the microcanonical ensemble in equilibrium) [47, 126, 127]. In order to understand the joint fluctuations it is therefore appropriate to perform a change of ensemble by introducing the Laplace transform of the reduced density matrix,

$$\begin{equation}{\rho }_{\lambda ,{\epsilon}}(t)=\sum\limits _{Q}\sum\limits _{A}{\rho }_{Q,A}(t){\text{e}}^{-\lambda Q-{\epsilon}A},\end{equation} \tag{ 16 }$$

with λ and different counting fields conjugated to the current and activity, respectively, and controlling their averages. By applying this transformation to the current- and activity-resolved master equation obtained from the unraveling of equation (4), we obtain a closed master equation for ρ_λ, [109],

$$\begin{align}{\dot {\rho }}_{\lambda ,{\epsilon}}(t)=& -i\left[H,{\rho }_{\lambda ,{\epsilon}}\right]+\sum\limits _{i=1}^{3}{L}_{i}^{\mathit{\text{D}}}{\rho }_{\lambda ,{\epsilon}}{L}_{i}^{\mathit{\text{D}}{\dagger}}-\frac{1}{2}\left\{{\rho }_{\lambda ,{\epsilon}},{L}_{i}^{\mathit{\text{D}}{\dagger}}{L}_{i}^{\mathit{\text{D}}}\right\}\\ & +{e}^{-\lambda -{\epsilon}}{L}_{0}^{+}{\rho }_{\lambda ,{\epsilon}}{L}_{0}^{+{\dagger}}-\frac{1}{2}\left\{{\rho }_{\lambda ,{\epsilon}},{L}_{0}^{+{\dagger}}{L}_{0}^{+}\right\}\\ & +\enspace {e}^{+\lambda -{\epsilon}}{L}_{0}^{-}{\rho }_{\lambda ,{\epsilon}}{L}_{0}^{-{\dagger}}-\frac{1}{2}\left\{{\rho }_{\lambda ,{\epsilon}},{L}_{0}^{-{\dagger}}{L}_{0}^{-}\right\}\equiv {\mathcal{L}}_{\lambda ,{\epsilon}}{\rho }_{\lambda ,{\epsilon}},\end{align} \tag{ 17 }$$

which defines ${\mathcal{L}}_{\lambda ,{\epsilon}}$, the tilted (or deformed) Liouvillian superoperator for the three-qubits dynamics, that no longer preserves the trace during the time evolution [51, 108, 109]. Interestingly, the moment generating function of the activity-current statistics is given by

$$\begin{equation}{Z}_{\lambda ,{\epsilon}}(t)\equiv \mathrm{Tr}[{\rho }_{\lambda ,{\epsilon}}(t)],\end{equation} \tag{ 18 }$$

which for long times also obeys a large deviation principle of the form Z_λ, (t) ≍ exp[+tμ(λ, )]. The function μ(λ, ) is nothing but the scaled cumulant generating function of the activity-current probability density function, and defines an additional LDF corresponding to the Legendre transform of G(q, a), i.e.

$$\begin{equation}\mu (\lambda ,{\epsilon})={\mathrm{max}}_{q,a}\left[G(q,a)-\lambda \enspace q-{\epsilon}\enspace a\right],\end{equation} \tag{ 19 }$$

a relation equivalent to the Legendre duality between different thermodynamic potentials [47, 126, 127]. It can be shown [108, 109], see also below, that μ(λ, ) is directly related to the spectral properties of the tilted superoperator ${\mathcal{L}}_{\lambda ,{\epsilon}}$ and the symmetry decomposition of the initial state. Note also that if we make λ = 0 = we recover the canonical (trace-preserving) Lindblad master equation (4). By inverting the Legendre transform we can conversely obtain the joint current-activity LDF G(q, a) –or at least its convex envelope (see below)–from the LDF μ(λ, ), i.e. $G(q,a)={\mathrm{max}}_{\lambda ,{\epsilon}}\left[\mu (\lambda ,{\epsilon})+\lambda \enspace q+{\epsilon}\enspace a\right]$.

On the other hand, if we make λ = 0 (or = 0) we recover the tilted Liouvillian superoperator for the activity (or current) statistics alone. In particular, let P_t (Q) be the probability of observing a total exciton current Q in a time t. This probability obeys for long times another large deviation principle P_t (Q) ≍ exp[+tF(q)] which defines the current LDF F(q)⩽0 and an associated scaled cumulant generating function for the current θ(λ) = max_q [F(q) − λq] = F(q_λ ) − λq_λ , with q_λ the current associated to a given λ, solution of the equation F'(q_λ ) = λ. Similarly, if P_t (A) is the probability of observing a total dynamical activity A in a time t, it can be shown to scale as P_t (A) ≍ exp[+tI(a)], with I(a)⩽0 the activity LDF such that ζ() = max_a [I(a) − a] = I(a ) − a is the scaled cumulant generating function for the activity. Here a is the activity for a given , solution of I'(a ) = . It is now easy to show that

$$\begin{equation}\theta (\lambda )=\mu (\lambda ,0),\hspace{25.0pt}\zeta ({\epsilon})=\mu (0,{\epsilon}).\end{equation} \tag{ 20 }$$

In this way, the kth-order cumulants of the current and the activity are given by

$$\begin{equation}{\langle {q}^{k}\rangle }_{c}={\left.-\frac{{\partial }^{k}\theta (\lambda )}{\partial {\lambda }^{k}}\right\vert }_{\lambda \to 0}={\left.-\frac{{\partial }^{k}\mu (\lambda ,0)}{\partial {\lambda }^{k}}\right\vert }_{\lambda \to 0}\enspace ,\end{equation} \tag{ 21 }$$

$$\begin{equation}{\langle {a}^{k}\rangle }_{c}={\left.-\frac{{\partial }^{k}\zeta ({\epsilon})}{\partial {{\epsilon}}^{k}}\right\vert }_{{\epsilon}\to 0}={\left.-\frac{{\partial }^{k}\mu (0,{\epsilon})}{\partial {{\epsilon}}^{k}}\right\vert }_{{\epsilon}\to 0}\end{equation} \tag{ 22 }$$

which correspond to the central moments of the associated distributions up to k = 3. Moreover, using Bayes theorem we can now define the conditional probability P_t (Q|A) = P_t (Q, A)/P_t (A) to observe a total exciton current Q given that the total activity is A, or similarly the conditional probability P_t (A|Q) = P_t (Q, A)/P_t (Q) of measuring a total activity A given a fixed total current Q. These conditional probabilities scale in the long-time limit in a large deviation form

$$\begin{equation}{P}_{t}(Q\vert A)\asymp \mathrm{exp}[+t{G}_{Q}(q\vert a)],\hspace{25.0pt}{P}_{t}(A\vert Q)\asymp \mathrm{exp}[+t{G}_{A}(a\vert q)],\end{equation} \tag{ 23 }$$

with

$$\begin{equation}{G}_{Q}(q\vert a)=G(q,a)-I(a),\hspace{25.0pt}{G}_{A}(a\vert q)=G(q,a)-F(q)\end{equation} \tag{ 24 }$$

the associated conditional LDFs.

In order to analyze the role of symmetry in the joint activity-current statistics of the three-qubits system, note first that the existence of the exchange symmetry operator π₁₂ implies that the associated symmetry superoperators ${{\Pi}}_{12}^{l,r}$ defined in equation (7) commute with the tilted Liouville superoperator ${\mathcal{L}}_{\lambda ,{\epsilon}}$ of equation (17), i.e. $[{{\Pi}}_{12}^{l,r},{\mathcal{L}}_{\lambda ,{\epsilon}}]=0$. Therefore there exists a complete biorthogonal basis of common left (${~{\omega }}_{\alpha \beta \nu }(\lambda ,{\epsilon})$) and right (ω_αβν (λ, )) eigenoperators in $\mathcal{B}(\mathcal{H})$, connecting eigenvalues of ${\mathcal{L}}_{\lambda ,{\epsilon}}$ to particular symmetry subspaces, such that

$$\begin{align}{{\Pi}}_{12}^{l}{\omega }_{\alpha \beta \nu }(\lambda ,{\epsilon})& ={\text{e}}^{\text{i}{\phi }_{\alpha }}{\omega }_{\alpha \beta \nu }(\lambda ,{\epsilon}),\\ {{\Pi}}_{12}^{r}{\omega }_{\alpha \beta \nu }(\lambda ,{\epsilon})& ={\text{e}}^{-\text{i}{\phi }_{\beta }}{\omega }_{\alpha \beta \nu }(\lambda ,{\epsilon}),\\ {\mathcal{L}}_{\lambda ,{\epsilon}}{\omega }_{\alpha \beta \nu }(\lambda ,{\epsilon})& ={\mu }_{\nu }(\lambda ,{\epsilon}){\omega }_{\alpha \beta \nu }(\lambda ,{\epsilon})\end{align} \tag{ 25 }$$

with ϕ_α = 0, π (and similarly for left eigenfunctions). Note that, due to the orthogonality of symmetry eigenspaces, Tr[ω_αβν (λ, )] ∝ δ_αβ , and we introduce in what follows the normalization Tr[ω_ααν (λ, )] = 1 for simplicity. The solution to equation (17) can be formally written as ${\rho }_{\lambda ,{\epsilon}}(t)=\mathrm{exp}(+t{\mathcal{L}}_{\lambda ,{\epsilon}})\rho (0)$, so a spectral decomposition of the initial density matrix in terms of the common basis yields ${Z}_{\lambda ,{\epsilon}}(t)={\sum }_{\alpha \nu }{\text{e}}^{+t{\mu }_{\nu }(\lambda ,{\epsilon})}\langle \langle {~{\omega }}_{\alpha \alpha \nu }(\lambda ,{\epsilon})\vert \rho (0)\rangle \rangle$ for the moment generating function of the activity-current statistics. For long times we thus have

$$\begin{equation}{Z}_{\lambda ,{\epsilon}}(t)t\to \infty {\to }{\text{e}}^{+t{\mu }_{0}^{({\alpha }_{0})}(\lambda ,{\epsilon})}\langle \langle {~{\omega }}_{{\alpha }_{0}{\alpha }_{0}0}(\lambda ,{\epsilon})\vert \rho (0)\rangle \rangle .\end{equation} \tag{ 26 }$$

Here ${\mu }_{0}^{({\alpha }_{0})}(\lambda ,{\epsilon})$ is the eigenvalue of ${\mathcal{L}}_{\lambda ,{\epsilon}}$ with largest real part and symmetry index α₀ among all symmetry diagonal eigenspaces ${\mathcal{B}}_{\alpha \alpha }$ (symmetric or antisymmetric, i.e. with α = A, S) with nonzero projection on the initial density matrix ρ(0). In this way, this eigenvalue defines the Legendre transform of the joint activity-current LDF, i.e. $\mu (\lambda ,{\epsilon})\equiv {\mu }_{0}^{({\alpha }_{0})}(\lambda ,{\epsilon})$, see equation (19) above. In other words, if ${\mu }_{0}^{(\alpha )}(\lambda ,{\epsilon})$ is the leading eigenvalue of ${\mathcal{L}}_{\lambda ,{\epsilon}}$ with symmetry index α, then

$$\begin{equation}\mu (\lambda ,{\epsilon})={\mathrm{max}}_{\alpha }[{\mu }_{0}^{(\alpha )}(\lambda ,{\epsilon})]={\mu }_{0}^{({\alpha }_{0})}(\lambda ,{\epsilon}),\end{equation} \tag{ 27 }$$

with the maximum taken over all symmetry subspaces with non-zero overlap with ρ(0), i.e. such that $\langle \langle {~{\omega }}_{\alpha \alpha 0}(\lambda ,{\epsilon})\vert \rho (0)\rangle \rangle \ne 0$. It is interesting to note that the long time limit in equation (26) selects a particular symmetry eigenspace α₀, effectively breaking at the fluctuating level the original symmetry of the three-qubits system. Remarkably, as shown in [108, 109], distinct symmetry eigenspaces may dominate different fluctuation regimes, separated by first-order-type DPTs. Moreover, the symmetry projections of the initial mixed state ρ(0) can be harnessed to control both the average transport properties of the three-qubit system and its joint activity-current statistics. As an example of this mechanism applied for currents, a symmetry-controlled quantum thermal switch was introduced in [108] that allows the control of heat flow using initial state preparation and symmetry tools, see also [12, 109].

We next show that the existence of a symmetry such as π₁₂ implies non-analyticities in the univariate LDFs θ(λ) and ζ() associated to the current and the activity, respectively, as well as in the joint LDF μ(λ, ) (see next section). These non-analyticities signal DPTs separating fluctuation regimes where the original symmetry is broken in different ways. For simplicity, we start with the current LDF θ(λ), see equation (20). We hence proceed by noting that ${\theta }_{0}^{(\alpha )}(\lambda )$, the leading eigenvalue of ${\mathcal{L}}_{\lambda ,0}$ with symmetry index α, can be expanded to first order for λ → 0 as

$$\begin{equation}{\theta }_{0}^{(\alpha )}(\lambda )\approx {\theta }_{0}^{(\alpha )}(0)+\lambda {\partial }_{\lambda }{\theta }_{0}^{(\alpha )}(\lambda ){\vert }_{\lambda =0}=-\lambda \langle {q}_{\alpha }\rangle ,\end{equation} \tag{ 28 }$$

where we have used that ${\theta }_{0}^{(\alpha )}(0)=0\enspace \forall \alpha$ due to the existence of a well-defined steady state ${\rho }_{\alpha }^{SS}$ in each symmetry sector α, see equation (8). Moreover, $\langle {q}_{\alpha }\rangle =-{\partial }_{\lambda }{\theta }_{0}^{(\alpha )}(\lambda ){\vert }_{\lambda =0}$ is the average current for the steady state ${\rho }_{\alpha }^{SS}$. Using now that the cumulant generating function for the current can be written as $\theta (\lambda )={\mathrm{max}}_{\alpha }[{\theta }_{0}^{(\alpha )}(\lambda )]$ [108, 109], we thus find

$$\begin{equation}\theta \enspace (\lambda )\underset{\vert \lambda \vert \to 0}{=}\begin{cases}+\vert \lambda \vert \langle {q}_{{\alpha }_{\mathrm{max}}}\rangle \quad & \quad \text{for}\enspace \lambda \lesssim 0\\ -\vert \lambda \vert \langle {q}_{{\alpha }_{\mathrm{min}}}\rangle \quad & \quad \text{for}\enspace \lambda \gtrsim 0\end{cases},\end{equation} \tag{ 29 }$$

where α_max (α_min) denotes the symmetry sector with maximal (minimal) average current $\langle {q}_{{\alpha }_{\mathrm{max}}}\rangle$ ($\langle {q}_{{\alpha }_{\mathrm{min}}}\rangle$) among those with nonzero overlap with ρ(0). Therefore the LDF θ(λ) will exhibit a kink at λ = 0 whenever $\langle {q}_{{\alpha \enspace }_{\mathrm{max}}}\rangle \ne \langle {q}_{{\alpha }_{\mathrm{min}}}\rangle$, characterized by a finite, discontinuous jump in the dynamic order parameter q_λ ≡ −θ'(λ) at λ = 0 of magnitude ${\Delta}{q}_{0}=\langle {q}_{{\alpha }_{\mathrm{max}}}\rangle -\langle {q}_{{\alpha }_{\mathrm{min}}}\rangle$, a behavior reminiscent of first order phase transitions [51, 108, 109]. The top-left panel in figure 3 shows the LDF θ(λ) measured for our three-qubits system for a general initial state (with nonzero projections on both the symmetric and antisymmetric sectors) and a particular set of parameters (B_z = 0.5, Γ = 0.1, γ = 0, and n = 0.1), and the presence of a kink at λ = 0 is apparent, as predicted. The symmetry sector corresponding to the minimal current phase is in this case the symmetric subspace, $\langle {q}_{{\alpha }_{\mathrm{min}}}\rangle =\langle {q}_{S}\rangle {< }0$, i.e. when qubits 1 and 2 are restricted to stay in any mixed state based on triplet states, see section 3 above. The reason is that, interestingly, and despite the presence of a single thermal reservoir, there is a net average current of excitons from the system to the bath in the symmetric steady state ${\rho }_{S}^{SS}$, with $\langle {q}_{S}\rangle =-{\theta }^{\prime }(\lambda ){\vert }_{\lambda \to {0}^{+}}{< }0$. This results from the nontrivial interplay between the Hamiltonian XX-interaction, the magnetic field along the z-direction, which induces rotation of the spins, and the thermal bath, which projects qubit 0 at a constant rate. On the other hand, as discussed in section 3, in the antisymmetric subspace spins 1 and 2 stay frozen into the singlet state, a dark (decoherence-free) state of the dynamics, effectively decoupling from the system evolution. In this symmetry sector the system thus behaves as a single qubit connected to a thermal reservoir, and the corresponding exciton current in the antisymmetric steady state ${\rho }_{A}^{SS}$ is hence $\langle {q}_{A}\rangle =-{\theta }^{\prime }(\lambda ){\vert }_{\lambda \to {0}^{-}}=0$, as deduced from the flat part of the kink at λ = 0⁻ of top-left panel in figure 3, so that ⟨q_A ⟩ > ⟨q_S ⟩.

Zoom In Zoom Out Reset image size

As a result of microscopic time reversibility (since the governing Lindblad superoperator obeys a local detailed balance condition [98–100]), the probability of every quantum jump trajectory is related to the probability of its time-reversed trajectory. Both trajectories share the same value of the activity (it is a time-symmetric observable), but the current sign is reversed as it falls into the time-antisymmetric sector. Consequently, the system will obey a Gallavotti–Cohen-type fluctuation theorem for the current statistics [95–97] (and the joint activity-current fluctuations, see below), linking the probability of a current fluctuation with its time-reversal event. This fluctuation theorem implies that θ(λ) = θ(κ − λ) for the cumulant generating function of the current, with κ a constant related to the rate of entropy production in the system. For our three-qubits system in contact with a thermal bath characterized by an average excitation number n, see section 2, we have that κ = ln[n/(n + 1)]. In this way, the kink in θ(λ) predicted at λ = 0 has a specular image at λ = κ, where a twin DPT emerges in current statistics. This twin kink is confirmed in the top-right panel of figure 3.

The behavior of the activity LDF ζ() can be analyzed in similar terms to the current. In particular, reasoning along the same lines it is easy to show that ζ() will show another kink around = 0, i.e.

$$\begin{equation}\zeta ({\epsilon})\underset{\vert {\epsilon}\vert \to 0}{=}\begin{cases}+\vert {\epsilon}\vert \langle {a}_{{\beta }_{\mathrm{max}}}\rangle \quad & \quad \text{for}\enspace {\epsilon}\lesssim 0\\ -\vert {\epsilon}\vert \langle {a}_{{\beta }_{\mathrm{min}}}\rangle \quad & \quad \text{for}\enspace {\epsilon}\gtrsim 0\end{cases},\end{equation} \tag{ 30 }$$

where now β_max (β_min) denotes the symmetry sector with maximal (minimal) average activity $\langle {a}_{{\beta }_{\mathrm{max}}}\rangle$ ($\langle {a}_{{\beta }_{\mathrm{min}}}\rangle$) among those with nonzero overlap with ρ(0). For the particular case of the three-qubits system with exchange symmetry studied in this paper, the maximum current and maximum activity subspaces do not coincide, as the maximum current subspace corresponds to the antisymmetric one while the maximum activity subspace corresponds to the symmetric sector, i.e. $\langle {a}_{{\beta }_{\mathrm{max}}}\rangle =\langle {a}_{S}\rangle$ and $\langle {a}_{{\beta }_{\mathrm{min}}}\rangle =\langle {a}_{A}\rangle {< }\langle {a}_{S}\rangle$. We stress however that this is not a necessary condition in general cases (for other open quantum systems with different symmetries, and/or different joint observables other than the current and the activity). Top-right panel in figure 3 shows the measured ζ() for the three-qubits system and the same particular set of parameters, confirming the presence of this additional kink. Notice also that, as opposed to the current, the dynamical activity is a time-symmetric observable which does not change sign upon time-reversal (indeed the activity is a positive-definite observable). Therefore no twin kink in ζ() is expected in this case, as confirmed in the top-right panel of figure 3.

As an interesting corollary, note that the kinks in the current LDF θ(λ) can only happen out of equilibrium, disappearing in equilibrium. In particular, the average currents for the multiple steady states are by definition zero in equilibrium, ⟨q_α ⟩ = 0∀α, so no symmetry-induced first-order DPT appears in θ(λ) at λ = 0 in equilibrium. On the other hand, the average activities of the different symmetry subspaces can be still different even when the system is in equilibrium, so the kink in the activity LDF ζ() and the associated activity DPT may still be present in equilibrium.

We may now obtain the current LDF F(q) from the inverse Legendre transform of θ(λ), i.e. F(q) = max_λ [θ(λ) + qλ]. It is then easy to show [67] that the twin kinks in θ(λ) correspond to two different current intervals, |q| ∈ (0, |⟨q_S ⟩|] (related by time-reversibility, q ↔ −q) where F(q) is affine or nonconvex [108], as corresponds to a multimodal current distribution P_t (Q) reflecting the dynamical coexistence of multiple transport channels (or steady states) classified by symmetry. Bottom-left panel in figure 3 shows F(q) as obtained by numerical inverse Legendre transform of θ(λ) in the top-left panel of the same figure. Note that the maximum current regime reduces to a single point in q-space as $\langle {q}_{{\alpha }_{\mathrm{max}}}\rangle =\langle {q}_{A}\rangle =0$. In this way, in order to sustain a given current fluctuation q such that |q| ∉ (0, |⟨q_S ⟩|], the open quantum system breaks the original symmetry and selects the particular symmetry sector that maximally facilitates a given current fluctuation: the statistics during a current fluctuation with |q| ⩾ |⟨q_S ⟩| is dominated by the symmetric subspace, whereas for zero current q = 0 the antisymmetric subspace prevails. Moreover, for currents |q| ∈ (0, |⟨q_S ⟩|] the dominant quantum jump trajectory spends some time t₀ = pt (p < 1) in the symmetric sector and a complementary time t − t₀ = t(1 − p) in the antisymmetric subspace, with p = |q/⟨q_S ⟩|, a sort of dynamical Maxwell-like construction [128]. Equivalent arguments hold for the activity LDF I(a), shown in the bottom-right panel of figure 3, which exhibits an affine or nonconvex regime for activities $a\in [\langle {a}_{{\beta }_{\mathrm{min}}}\rangle ,\langle {a}_{{\beta }_{\mathrm{max}}}\rangle ]=[\langle {a}_{A}\rangle ,\langle {a}_{S}\rangle ]$ as a result of the kink in ζ() at = 0, with ⟨a_A ⟩ < ⟨a_S ⟩. Note however that, due to the time-symmetric character of the activity, no twin DPT is expected in this case (note also that a ⩾ 0 in all cases).

Next we study how the average current and activity of the different symmetry-classified steady states behave as a function of the average number of excitations n in the thermal bath, see top-right and bottom-right panels in figure 4. We first note that, while the average activity in both the symmetric and antisymmetric sectors increases with n (bottom-right panel in figure 4), the absolute value of the average current in the symmetric sector decreases instead with n (the current in the antisymmetric sector vanishes ∀n as explained above). In this way the average current in the symmetric sector recedes from its activity-related bound, −⟨a_S ⟩ ⩽ ⟨q_S ⟩ ⩽ ⟨a_S ⟩, as n increases. The activity of both the symmetric and antisymmetric steady states grows linearly with n for large enough n, although the activity of the symmetric sector is always (slightly) above the one of the antisymmetric sector. Note however that activity differences between both symmetry sectors are only apparent for low values of the bath average excitation number n, see top inset in the bottom-right panel of figure 4. Moreover, the differences in the transport and activity patterns between the symmetric and antisymmetric sectors of the dynamics tends to vanish as n increases, thus reducing the range of controllability of the quantum current and activity, possible by tuning the symmetry projections of the initial state [12, 108, 109]. We also studied the dependence with n of the typical current q_λ = −θ'(λ) for a bias parameter λ, and the typical activity a = −ζ'() for bias , see top left and bottom-left panels in figure 4, respectively. The discontinuities in both q_λ and a around the kinks in their respective LDFs are apparent. Moreover, the width of the discontinuous gap in q_λ , denoted here as δ and associated to the regime in λ-space dominated by the antisymmetric sector, quickly decreases with λ. This gap is related to the λ-distance between a given event and its time-reversal, and is simply given by δ = −κ = ln[(n + 1)/n], see inset in top-left panel of figure 4. In comparison, the size of the discontinuous jumps in q_λ , related to the difference in average currents between the symmetric and antisymmetric sectors, decreases at a much slower pace with n. For the activity jump, on the other hand, the decrease with n is much faster.

Zoom In Zoom Out Reset image size

Figure 5 explores the dependence of the same observables with the strength of the external magnetic field B_z . Interestingly, both the average current and the average activity of the symmetric steady state depend non-monotonously on the strength of the magnetic field, see top-right and bottom-right panels in figure 5, while B_z does not affect the values of ⟨q_A ⟩ and ⟨a_A ⟩. In particular, there exists a particular (common) value of B_z where the differences ⟨q_S ⟩ − ⟨q_A ⟩ and ⟨a_S ⟩ − ⟨a_A ⟩ are maximal. This may be used to optimize the range of controllability of the transport and activity properties of the three-qubits systems [12, 108, 109]. Due to the activity constraint on the current, −⟨a_S ⟩ ⩽ ⟨q_S ⟩ ⩽ ⟨a_S ⟩, the decrease of ⟨a_S ⟩ with B_z forces ⟨q_S ⟩ to diminish also as the magnetic field increases, a sort of enslaved behavior. Note also that, contrary to what happens when n is changed, see figure 4, the width of the discontinuous gap δ in q_λ does not depend on the magnetic field intensity B_z , see top-left panel in figure 5, indicating that B_z does not affect the rate of entropy production in the system.

Zoom In Zoom Out Reset image size

We next turn to investigate the joint LDF μ(λ, ), a bivariate LDF for which the arguments are similar to those discussed in previous section but some care is needed. Proceeding as in section 6, see e.g. equation (28), we first note that the leading eigenvalue of the tilted superoperator ${\mathcal{L}}_{\lambda ,{\epsilon}}$ with symmetry index α, denoted as ${\mu }_{0}^{(\alpha )}(\lambda ,{\epsilon})$, can be expanded to first order for λ, → 0 as

$$\begin{equation}{\mu }_{0}^{(\alpha )}(\lambda ,{\epsilon})\underset{\vert \lambda \vert ,\vert {\epsilon}\vert \ll 1}{\approx }{\mu }_{0}^{(\alpha )}(0,0)+\lambda {\partial }_{\lambda }{\mu }_{0}^{(\alpha )}(\lambda ,0){\vert }_{\lambda =0}+{\epsilon}{\partial }_{{\epsilon}}{\mu }_{0}^{(\alpha )}(0,{\epsilon}){\vert }_{{\epsilon}=0}=-\lambda \langle {q}_{\alpha }\rangle -{\epsilon}\langle {a}_{\alpha }\rangle ,\end{equation} \tag{ 31 }$$

where we have used in the second equality that ${\mu }_{0}^{(\alpha )}(0,0)=0\enspace \forall \alpha$ due to stationary conditions in each symmetry subspace, as above. Using now that the joint LDF $\mu (\lambda ,{\epsilon})={\mathrm{max}}_{\alpha }[{\mu }_{0}^{(\alpha )}(\lambda ,{\epsilon})]$, we thus find that $\mu (\lambda ,{\epsilon})\underset{\lambda ,{\epsilon}\to 0}{=}{\mathrm{max}}_{\alpha }[-\lambda \langle {q}_{\alpha }\rangle -{\epsilon}\langle {a}_{\alpha }\rangle ]$. Since we are working in the limit where λ, → 0, both at comparable rates (otherwise the faster-decaying biasing field would be effectively zero in this limit, thus reducing the discussion to that in the previous section for univariate LDFs), we can write = uλ with u some proportionality constant so

$$\begin{equation}\mu (\lambda ,{\epsilon})\underset{\lambda ,{\epsilon}\to 0}{=}{\mathrm{max}}_{\alpha }[-\lambda {M}_{\alpha }(u)]=\begin{cases}-\lambda \enspace {\mathrm{min}}_{\alpha }\enspace {M}_{\alpha }(u)\quad & \quad (\lambda { >}0)\\ \vert \lambda \vert \enspace {\mathrm{max}}_{\alpha }\enspace {M}_{\alpha }(u)\quad & \quad (\lambda {< }0)\end{cases},\end{equation} \tag{ 32 }$$

where we have defined a linear function

$$\begin{equation}{M}_{\alpha }(u)=\langle {q}_{\alpha }\rangle +u\langle {a}_{\alpha }\rangle \end{equation} \tag{ 33 }$$

of slope ⟨a_α ⟩ and intercept ⟨q_α ⟩, characteristic of each symmetry subspace with index α. Figure 6 shows a sketch of M_α (u) for the three-qubits system studied here and the particular set of parameters used e.g. in figure 3 (B_z = 0.5, Γ = 0.1, γ = 0, n = 0.1), with α = S (symmetric subspace) and α = A (antisymmetric sector). In general, there exists a threshold value

$$\begin{equation}{u}_{0}=-\frac{\langle {q}_{S}\rangle -\langle {q}_{A}\rangle }{\langle {a}_{S}\rangle -\langle {a}_{A}\rangle },\end{equation} \tag{ 34 }$$

defined by the equality M_S (u₀) = M_A (u₀), such that M_A (u) < M_S (u) for u > u₀ while M_S (u) < M_A (u) for u < u₀, see figure 6. In this way, the symmetry sector dominating the joint activity-current fluctuations near the steady state will depend on how we approach the origin in (λ, )-space, i.e. on the particular value of parameter u = /λ. For λ > 0 we have that $\mu (\lambda ,{\epsilon})\underset{\lambda ,{\epsilon}\to 0}{=}-\lambda \;{\mathrm{min}}_{\alpha }{M}_{\alpha }(u)$, so for u > u₀ the antisymmetric sector dominates activity-current fluctuations, while for u < u₀ the symmetric subspace is responsible of activity-current fluctuations, see figure 6. Conversely, for λ < 0 we have that $\mu (\lambda ,{\epsilon})\underset{\lambda ,{\epsilon}\to 0}{=}\vert \lambda \vert \;{\mathrm{max}}_{\alpha }{M}_{\alpha }(u)$, so for u > u₀ the symmetric sector prevails, while the antisymmetric subspace dominates for u < u₀. Therefore μ(λ, ) exhibits a kink line (i.e. with discontinuous derivative) near (λ, ) → 0, with local slope u₀ near the origin, such that the symmetric sector dominates below this line while the antisymmetric subspace prevails above it.

Zoom In Zoom Out Reset image size

Figure 7 (left panel) shows the measured LDF μ(λ, ) for the three-qubits model as a function of λ and , for parameters B_z = 0.5, Γ = 0.1, γ = 0 and n = 0.1, as in previous plots. Interestingly, the infinitesimal kink segment predicted around λ = 0 = of local slope u₀ is confirmed, and extends over the entire (λ, )-plane into a line of first-order DPTs along which dynamical coexistence of the different (symmetric and antisymmetric) transport channels appears. Furthermore, due to the microscopic time-reversibility of the open quantum dynamics [98–100], the joint activity-current LDF G(q, a) obeys a Gallavotti–Cohen-type fluctuation theorem along the current (time-antisymmetric) axis [95–97], which for the Legendre-dual LDF can be simply written as μ(λ, ) = μ(κ − λ, ), where κ = ln[n/(n + 1)] has been defined above. This immediately implies a twin kink branch in the (λ, )-plane signaling a twin line of DPTs, as confirmed in the left panel of figure 7. In this way, the (λ, )-plane is divided into two regions, an inner zone where the associated joint activity-current fluctuations are dominated by the antisymmetric subspace, and an outer region where the symmetric sector prevails. The presence of a double kink along the λ-axis and a single kink along the -axis can be confirmed by representing constant- and constant-λ slices of the LDF μ(λ, ), see respectively top-right and bottom-right panels in figure 7.

Zoom In Zoom Out Reset image size

As discussed in section 4, a noisy dephasing channel acting on all qubits (γ ≠ 0, see equation (3) and section 2) breaks the exchange symmetry of the original system, restoring global ergodicity and leading to a unique steady state [102, 120]. The presence of this symmetry-breaking dephasing noise then immediately implies the disappearance of the symmetry-induced DPTs and the kinks in μ(λ, ) described above, as well as the kinks in the univariate LDFs θ(λ) and ζ(). This is illustrated in figure 8, where constant- and constant-λ slices of the LDF μ(λ, ) are represented, both in the absence (γ = 0) and in the presence (γ ≠ 0) of the dephasing channel. The existence of the kink is apparent in the case with no dephasing (solid lines), but it disappears when γ ≠ 0 (dashed lines). It is also remarkable that, in the symmetric-subspace phase, μ(λ, ) is very similar with and without dephasing, while in the antisymmetric-subspace regime the difference is appreciable, e.g. μ(λ, ) vs λ for constant- is flat in this antisymmetric regime when γ = 0, but it seems to continue analytically the form of μ(λ, ) in the symmetric case when γ ≠ 0, see left panel in figure 8. This happens because the action of the dephasing channel changes dramatically the symmetry properties of the system, but not its transport properties. The observed behavior is another indication that, once the symmetry is broken by the dephasing channel, the dominant subspace in the dynamical evolution is the symmetric one (due to entropic reasons), as already discussed in section 4.

Zoom In Zoom Out Reset image size

We can now obtain the joint current-activity LDF G(q, a) by inverse Legendre transforming μ(λ, ), see left panel in figure 9. First, it is important to note that the activity constraint on the current, −a⩽ q ⩽ a, has significant consequences in the joint activity-current statistics. For instance, the constraint implies that the absolute value of the current cannot be larger than the activity under any circumstances, so G(q, a) → −∞ for any |q| > a and therefore G(q, a) takes finite values only in a triangle defined by the constraint |q| ⩽ a. On the other hand, since both the symmetric and antisymmetric subspaces have definite values of the average current, ⟨q_A ⟩ = 0 and ⟨q_S ⟩ < 0 respectively, the constraint −a ⩽ q ⩽ a implies that there exists a critical value for the activity a_c = |⟨q_S ⟩| such that the current for any activity a < a_c can only be q = ⟨q_A ⟩ = 0. This immediately implies that G(q = 0, a < a_c) = 0 while G(q ≠ 0, a < a_c) = −∞, see main plot in figure 9. Remarkably, this clear-cut observation opens up a new route to control quantum transport: by biasing the activity of our three-qubits system below the critical activity a_c, one is able to shut down completely the exciton current in the system, since in this fluctuation regime the antisymmetric sector prevails. This novel activity-driven current lockdown regime is enabled by symmetry, and suggests unexplored quantum control strategies. In addition, proceeding as in section 6, one can easily show [67, 108] that the twin kink branches in μ(λ, ) correspond (after the inverse Legendre transform) to two different regions in the (q, a)-plane where the G(q, a) is affine or nonconvex, signaling the dynamical coexistence of the two different transport channels in these current-activity regions. These affine or nonconvex zones are clearly visible in the main panel of figure 9, and comprise two well-defined bands (−|⟨q_S ⟩|, 0) ∪ [a_c, +∞) and (0, |⟨q_S ⟩|) ∪ [a_c, +∞) in (q, a)-space. Note that between these two zones there is a q = 0 line ∀a corresponding to the antisymmetric manifold. Note also that negative currents are more probable than positive ones, see color legend in left panel of figure 9, and hence the average current is negative. For clarity, we have represented by a thick orange solid curve the isoline λ = 0 in the (q, a)-plane, that marks the average current ⟨q_a ⟩ for a given activity a, while the dashed orange line represents the = 0 isoline capturing the average activity ⟨a_q ⟩ for a given current q. The λ = 0 isoline exists only in the negative current half-plane, as it is there where the typical behavior occurs in the absence of bias on the current, while the = 0 isoline propagates through both the positive and negative currents half-planes since a given typical activity can be associated with both positive or negative currents.

Zoom In Zoom Out Reset image size

Using the measured joint LDF G(q, a), see figure 9, and the univariate LDFs F(q) and I(a) obtained in section 6, see figure 3, it is now possible to study the conditional LDFs G_Q (q|a) and G_A (a|q) defined in equation (24). This is shown in the right panels of figure 9. In particular, the top-right panel shows the current LDF conditioned to a fixed value of the activity, G_Q (q|a) = G(q, a) − I(a), while the bottom-right panel displays the activity LDF conditioned to a fixed value of the current, G_A (a|q) = G(q, a) − F(q). The LDF G_Q (q|a) exhibits for a > a_c two symmetrical current intervals around q = 0 where it is affine or nonconvex, while for a < a_c it is only defined for q = 0, as expected from the joint activity-current fluctuation behavior observed in G(q, a), see left panel in figure 9. Interestingly, biasing the activity to high values beyond its average behavior leads to an increase in the probability of negative current fluctuations. On the other hand, the probability of high positive current fluctuations is very small and almost independent of the activity a, see the tails of G_Q (q|a) in the top-right panel of figure 9. The statistics of the activity conditioned on a given current is shown in the bottom-right panel of figure 9. The gray area in this plot represents the current intervals around q = 0 where the bivariate G(q, a) is affine or nonconvex. As q increases, the mean activity conditioned on this value of the current grows, as does the probability of high activity fluctuations. Note also that the activity constraint on the current, −a ⩽ q ⩽ a, implies that G_A (a < q|q) → −∞ so G_A (a|q) jumps discontinuously from a finite value to −∞ at a = |q|. Finally, due to the time-symmetric character of the activity, the statistics of activity conditioned on a current q is the same when conditioned on a current −q, and hence G_A (a|q) = G_A (a| − q). This is another instance of the Gallavotti–Cohen theorem based on microscopic time-reversibility.

Even if these results are specific for a three-qubit system, we emphasize that similar features will arise for bigger systems. Interestingly, it has been shown that systems with more than three qubits can exhibit multiple strong symmetries and associated DPTs [108, 109, 129], as well as dynamical symmetries [130]. In this more complex scenario the system of interest would present several invariant subspaces, and the relation between current and activity might be more complicated than in the three-qubit case, though most of the phenomenology presented here may still hold.

To end this section, we now characterize the distinct dynamical patterns in the different symmetry phases of the dynamics, and how they coexist dynamically, with a distinctive intermittent pattern, in the presence of a weak dephasing noise channel. For that, we perform quantum Monte Carlo simulations of individual quantum jump trajectories, as in section 3, and measure an appropriate order parameter capable of distinguishing between both types of dynamics. As discussed in previous sections, in the antisymmetric state the system is reduced to a single qubit (as qubits 1 and 2 fall into a dark state and dynamically decouple) and the current is exactly zero while there is a net activity in the system. This is only possible because there is always one excitation entering and one coming out, a sort of excitonic blinking pattern. This type of locked blinking dynamics does not happen in general in the symmetric sector. To capture this essential dynamical signature we now define two different observables, C^±, such that C⁺ = +1 whenever two consecutive quantum jumps introduce excitations in the system (C⁺ = 0 otherwise), while C⁻ = −1 whenever two consecutive quantum jumps remove excitations from the system (C⁻ = 0 otherwise). Figure 10 shows the time evolution of these two observables for different situations. In particular, the top panels show the dynamics of the three-qubits system in the absence of dephasing channel (i.e. when γ = 0, see equation (3) and section 2), starting with a symmetric initial state (left) or an antisymmetric one (right). Clearly, dynamics in the symmetric case is characterized by many consecutive double jumps both up and down, but with net prevalence of double exciton removal jumps that leads to a negative average current in the symmetric steady state. On the other hand, the locked blinking dynamics in the antisymmetric sector implies that both C^± remain exactly 0 along the whole time evolution, see top-right panel in figure 10, and the current is always zero. This clear difference in the dynamics of C^± confirms the validity of these observables as dynamical order parameters for the different symmetry sectors. The presence of dephasing noise, on the other hand, allows the mixing between both symmetry subspaces. Bottom panel in figure 10 shows the time evolution of our dynamical order parameters C^± for a weak dephasing amplitude γ = 0.01. Interestingly, the system evolution in this case exhibits intermittent behavior, with periods of time where the system is trapped in the symmetric manifold, interrupted by jumps to the antisymmetric manifold allowed by the weak mixing introduced by the dephasing channel. This intermittent evolution is typical of a dynamical coexistence between phases of distinct activity, and is a direct consequence and a dynamical signature of the underlying symmetry of the three-qubits system. Such dynamical signatures can be harnessed to infer molecular symmetries from nonequilibrium transport experiments [110].

Zoom In Zoom Out Reset image size

In this paper we have investigated how strong symmetries affect both the transport properties and the activity patterns of a particular class of Markovian open quantum system, a three-qubits model under the action of a magnetic field and in contact with a thermal bath. Strong symmetries in open quantum systems lead to broken ergodicity and the emergence of multiple degenerate steady states [102]. A first observation in this work is that, interestingly, for initial states overlapping with several symmetry subspaces, individual quantum jump trajectories select randomly one of the symmetry sectors, collapsing in a finite time to the corresponding subspace and remaining there from that time on. This a particular instance of the dissipative freezing phenomenon recently observed in [124], which implies a breakdown of a conservation law (associated to the underlying symmetry) at the individual trajectory level [124].

From a quantum jump perspective, the appearance of multiple steady states mentioned above is related to underlying DPTs at the fluctuating level, that lead to a dynamical coexistence of different transport/activity channels classified by symmetry. We have studied these DPTs in the univariate LDFs associated to the exciton current (a key time-antisymmetric observable characterizing transport out of equilibrium) and the dynamical activity (a time-symmetric magnitude of direct experimental relevance which may constraint the range of current fluctuations). In particular, we find a pair of twin DPTs in exciton current statistics, induced by the strong symmetry and related by time reversibility, where a zero-current antisymmetric phase (under the exchange of qubits 1 and 2) coexists with a symmetric phase of negative exciton current. On the other hand, the activity statistics exhibits a single DPT (since the activity is a time-symmetric observable) where the symmetric and antisymmetric phases of different but nonzero activities dynamically coexists. Interestingly, the maximum current and maximum activity subspaces do not coincide for the three-qubits model studied here, as the maximum current subspace corresponds to the antisymmetric one while the maximum activity subspace corresponds to the symmetric sector.

In addition, this work also discusses how symmetries are reflected in the joint large deviation statistics of the activity and the current, a central observable in order to fully characterize the complex, coupled quantum jump dynamics both in the time-symmetric and time-antisymmetric sectors. The presence of a strong symmetry under nonequilibrium conditions implies non-analyticities in the dynamical free energy in the dual activity-current plane (or equivalently in the joint activity-current LDF). Remarkably, the DPT predicted around the steady state and its Gallavotti–Cohen twin dual are extended into lines of first order DPTs in the current-activity plane, with a nontrivial structure which depends on the transport and activity properties of each of the symmetry phases (in particular on the average current and activity of each of these phases). We further find that activity constraints the range of current fluctuations, leading in particular to an activity-driven current lockdown phase for activities below some critical threshold, a new route to control quantum transport enabled by symmetry. The presence of a noisy dephasing channel acting on all qubits breaks the exchange symmetry of the three-qubits system, restoring global ergodicity and leading to an unique steady state. This dephasing noise also washes out the symmetry-induced DPTs, although the underlying topological symmetry leaves a dynamical fingerprint in the form of an intermittent, bursty on/off dynamics between the different symmetry sectors when the dephasing amplitude is weak, a phenomenon observed in quantum Monte Carlo simulations of individual quantum jump trajectories, and well-captured by some blinking order parameters C^±.

We acknowledge the Spanish Ministry and Agencia Estatal de Investigación (AEI) through Grant FIS2017-84256-P (European Regional Development Fund), as well as the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund, Ref. A-FQM-175-UGR18 and SOMM17/6105/UGR for financial support. We are also grateful for the computational resources and assistance provided by PROTEUS, the supercomputing center of Institute Carlos I for Theoretical and Computational Physics at the University of Granada, Spain.

The data that support the findings of this study are available upon reasonable request from the authors.