Metastability and quantum coherence-assisted sensing in interacting parallel quantum dots

We study the transient dynamics subject to quantum coherence effects of two interacting parallel quantum dots weakly coupled to macroscopic leads. The stationary particle current of this quantum system is sensitive to perturbations much smaller than any other energy scale, specifically compared to the system-lead coupling and the temperature. We show that this is due to the presence of a parity-like symmetry in the dynamics, as a consequence of which, two distinct stationary states arise. In the presence of small perturbations breaking this symmetry, the system exhibits metastability with two metastable phases that can be approximated by a combination of states corresponding to stationary states in the unperturbed limit. Furthermore, the long-time dynamics can be described as classical dynamics between those phases, leading to a unique stationary state. In particular, the competition of those two metastable phases explains the sensitive behavior of the stationary current towards small perturbations. We show that this behavior bears the potential of utilizing the parallel dots as a charge sensor which makes use of quantum coherence effects to achieve a signal to noise ratio that is not limited by the temperature. As a consequence, the parallel dots outperform an analogous single-dot charge sensor for a wide range of temperatures.


I. INTRODUCTION
The coherent control of electronic quantum devices is a challenging necessity for many novel device concepts [1]. Some device concepts are based on quantum coherent non-equilibrium charge transport, while in other cases measurements of charge transport can be used to gain information about the quantum system. In many cases, the stationary properties are of interest, but with increasing control over quantum systems understanding the full transient behavior is of importance. For any such application the dynamics beyond the stationary state becomes relevant and attracted much attention in recent years [2][3][4][5][6][7]. Generally, the relaxation of a quantum system can be complicated and it can experience a quasistationary state before relaxing into its true stationary state. This phenomenon of metastability [8][9][10][11] occurs when the timescales dictating the system's dynamics are well separated. Importantly, metastability always occurs for systems in the proximity to multi-stable points where the dynamics features multiple stationary states, which may arise, e.g., due to a symmetry of the dynamical equations. It can also be observed in constrained systems such as quantum spin glasses [12,13], quantum gases [14,15] and superconducting nanojunctions [16].
Another ingredient to understand a quantum system's behavior lies with quantum interference and coherence effects and how they affect the system's properties and dynamics. Different types of quantum dot systems constitute particularly well-controlled and versatile platforms to study and use various aspects of interference effects [17][18][19][20]. It has been suggested that they can, e.g., reduce power fluctuations and rectify heat transport [21,22], or affect thermoelectric properties [23] and electronic trans-port properties [24] in molecules. One can also exploit the quantum interference to construct a transistor [25]. Particularly, interference effects in the transport properties of parallel double quantum dots have been a longstanding topic of interest. On the one hand, in the strong coupling regime, strongly-correlated physics dominates the transport exhibiting the Kondo effect and leading to defined signatures in the conductance as well as population switching [26][27][28][29][30].
On the other hand, the coherence effects present in parallel double quantum dots, understood as a superposition state of a single fermion on either the upper or lower dot, have been shown to play a role for quantum thermodynamics [31], and impact the thermoelectric current [32], thermal conductance [33] and electric transport [34,35] in the weak coupling limit.
In this paper we use a quantum master equation approach to study the non-equilibrium transport properties of such a parallel dots system, sketched in Fig. 1(a). The stationary state of this quantum system has been shown, in certain parameter regimes, to be sensitive to small perturbations in the system-lead coupling and the detuning of the dot energies [34,35]. Remarkably, due to quantum coherence effects, changes that are much smaller than all other energy scales in the system can lead to large changes in the stationary particle current [35].
We show that the sensitive response results from the presence of two stationary states when both the dot energies and their tunnel couplings to the leads are identical. Perturbations breaking the corresponding symmetry introduce a large timescale which is well separated from any other timescale in the system's dynamics. For intermediate times, metastability occurs, and the state of the system is well approximated by a probabilistic mix-ture of two metastable phases [9][10][11]. In the long time limit, those probabilities evolve according to classical dynamics dominated by the emergent timescale due to the broken symmetry. As the long-time dynamics depend not only on the size but also on the structure of perturbations breaking the symmetry, the stationary state does as well, and the stationary current displays large changes resulting from small parameter changes. In particular, we focus on the regime of large Coulomb interaction where one of the metastable phases features suppressed particle current through the system, while the other phase supports a much larger current, so that the stationary current varies from suppressed to larger values.
Furthermore, we investigate how the sensitive behavior of the current can be used to enable the system to act as a charge sensor. To quantify the accuracy of the sensor we also need the current noise, which we calculate based on counting statistics [36,37]. One problem for sensing applications is that metastability typically results in large current noise. Another problem is the long relaxation times associated with metastability. Nonetheless, we show that there is a large parameter regime where the parallel dots by far outperform a single dot used as a charge sensor.
The paper is organized as follows. After introducing the model in Section II, we discuss the Lindblad dynamics of the parallel dots in Section III. In Section IV we discuss the transient dynamics. This includes a discussion of the unperturbed and perturbed dynamics in the context of symmetry breaking and the resulting metastability and the long time dynamics towards a unique stationary state. Finally, in Section V we investigate how the parallel dots could be used as a charge sensor.

II. MODEL
The setup under consideration consists of two interacting parallel quantum dots weakly coupled to macroscopic leads, see Fig. 1(a). To simplify the analytic treatment we consider spinless electrons in the following, but we have verified by direct comparison that the qualitative physics and results remain the same for spin-degenerate dot orbitals (our spinless model can be realized in a system where the Zeeman energy is larger than the applied bias voltage).
The Hamiltonian of the entire setup splits into three parts, for the parallel dots, the leads and the interaction between the subsystems. The dot Hamiltonian is given by with the fermionic creation and annihilation operators, d † j and d j , where j = 1, 2 is the dot label, the energy levels of the dots are j and U denotes the Coulomb interaction; we set = 1 throughout the paper. The leads are described by non-interacting fermions, The operators c † ks , c sk create or annihilate an electron in the left (s = L) or right (s = R) lead at momentum k, with their corresponding energy dispersion given by ω sk . Finally, the tunneling between the parallel dots and the leads is governed by with the tunneling amplitude t jsk . We note that H T assumes both dots to couple to the same lead channel. This is a crucial assumption for the results of this work, which will hold for one-dimensional leads, but is also justified when the dots are close on a length scale set by the Fermi wavelength in the leads. We furthermore take t jsk to be real and positive (we can always make this choice if relative phases between t 1sk and t 2sk are independent of k and s). If a magnetic field is present, or the leads are ferromagnetic or superconducting, one might need to consider complex t jsk , which possibly leads to additional phenomena, such as phase lapses in the conductance [27,28,38,39]. The electrons in the leads are assumed to be described by the grand canonical ensemble with the chemical potentials µ s and temperatures T s ; we set Boltzmann's constant k B = 1 and the elementary charge e = 1. The chemical potentials can be controlled by the bias voltage V B = 2µ L = −2µ R , while the gate allows control of the dot energy levels, V G = −( 1 + 2 )/2.

III. LINDBLAD DYNAMICS OF PARALLEL DOTS
A. Master equation In this work, we focus on the limit of weak tunnel coupling. In this case, the dynamics of the reduced density matrix ρ PD (t) of the two quantum dots can be well approximated by a quantum master equation of Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) form [40,41].
We assume that tunneling amplitudes are independent of the momentum t jsk = t js , and that the density of states is constant, ν sk = ν. The tunneling amplitudes t js define the tunneling rates as Γ js = 2πν|t js | 2 . The weak coupling limit where our master equation is valid is defined by Γ js T s . Following [42][43][44], we consider terms quadratic in tunneling amplitudes to arrive at the master equation beyond the secular approximation, where [·, ·] and {·, ·} stand for the commutator and anticommutator, respectively. Here, the effective Hamiltonian H eff = H PD +H LS includes a Lamb shift H LS , which renormalizes the energies of the parallel dots and can be identified as an effective tunnel splitting [26]. We exclude a direct interdot tunnel coupling of the form (2), but such a term can easily be added to H PD and the qualitative physics described in the following does not change for Ω 12 T s . The jump operators J αs describe the exchange of an electron between the parallel dots and lead s with α = + representing an electron entering the dots and α = − an electron leaving. The closed form expressions for the Lamb shift and the jump operators are derived in Appendix A.

B. Spectral decomposition and metastability
The equation of motion in Eq. (5) can be recast as where L is the Liouvillian. Therefore, the full time evolution, which is a completely positive, trace-preserving map, can be formally solved as It follows that the evolution can be decomposed in terms of the Liouville operator spectrum, that is, its eigenvalues λ i and left and right eigenmatrices, L i and R i , as with the coefficients c i = Tr[L i ρ PD (0)] and the left and right eigenmatrices normalized so that Tr(L i R j ) = δ ij .
Here, the eigenvalues are ordered with a decreasing real part, so that λ 1 = 0 corresponds to a stationary state R 1 = ρ ss PD , while L 1 = 1, where 1 is the identity operator on the dots. When the stationary state is unique, the sum runs over the decay modes with Re(λ i ) < 0, so that lim t→∞ ρ PD (t) = ρ ss PD for any initial state. If there exists a large difference in the real parts of the second and third eigenvalues, −λ 2 −Re(λ 3 ), metastability arises followed by long-time dynamics dominated by the single low-lying eigenmode [9,10]. Here, λ 2 is necessarily real as L preserves the Hermiticity of ρ PD (t) and therefore complex eigenvalues need to appear as complex conjugate pairs. Indeed, for times t such that −Re(λ 3 )t 1, the state of the system can be approximated as For times t such that −λ 2 t 1, the system is metastable, with its state approximated by a linear combination of the stationary state and the low-lying eigenmode, ρ PD (t) ≈ ρ ss PD + c 2 R 2 , where c 2 carries the information about the initial condition. For longer times, the decay of the low-lying mode in Eq. (9) can no longer be neglected. In particular, when −λ 2 t 1, the system state approaches its asymptotic limit and is well approximated by the stationary state, ρ ss PD , independently of the initial condition. Thus, −1/Re(λ 3 ) and −1/λ 2 can be considered as the timescales of the initial and final relaxation, respectively.
In this work, we show that for the dynamics in Eq. (5), such metastability emerges as a consequence of breaking a parity-like symmetry originating from the Hamiltonian in Eq. (1). For the case of a single low-lying eigenvalue, metastable states correspond to probabilistic mixtures of two metastable phases and can be investigated numerically [9,10], see also Appendix B. Here, we uncover the metastable phases, together with the unique stationary state analytically by means of non-Hermitian perturbation theory [45]. C. Dynamics of particle currents The average particle current leaving lead s at time t is given by I s (t) = −i [H, N s ] , with the electron number operator in the lead N s = k c † sk c sk . Within the Lindblad dynamics, the current I s (t) is given in terms of the jump operator by [37,46] Asymptotically, the currents equilibrate, I L + I R = 0, with I s = lim t→∞ I s (t) denoting the stationary current leaving lead s. In the remainder of the paper, we will therefore study the current I L (t) leaving the left lead and we will omit the lead index. Beyond its average the current dynamics can be investigated in terms of full counting statistics [36,47,48] used in Section V, see also Appendix C. Figure 1(b) shows the differential conductance dI/dV B as a function of V G and V B . It has a much richer structure than the results for a double dot system with a density matrix assumed to be diagonal in the eigenbasis of H PD and thus evolving according to a Pauli rate equation, rather than by Eq. (5); see Appendix D. This is due to coherences between singly occupied states playing a nonnegligible role in the properties of both the dynamics and the stationary state, especially when model parameters are chosen in the proximity to those for which strong symmetries are present. But we also clarify in which limits the Pauli rate equation reproduces the dynamics.
The stationary limit of the transport properties of the parallel dots system has been studied before [34,35] and it has been reported that the stationary current may be highly sensitive to perturbations in the tunneling rates and to detuning of the dot energies [35], see also Fig. 4(a). In this work, we explain this phenomenon in relation to symmetry breaking and demonstrate that the sensitivity of the stationary current can in fact be arbitrarily large. We then verify its usefulness for sensing applications.

IV. SYMMETRIES AND DYNAMICS
We now discuss symmetries of the Hamiltonian dynamics for the total setup and the resulting properties of the Lindblad dynamics of the dots. In particular, we show how two distinct stationary states occur as a result of a swap symmetry between the dots present in the Hamiltonian. We then analyze the metastability arising by perturbatively breaking this symmetry, the long-time dynamics that follows, and the resulting unique stationary state and the associated current. Crucially, the structure of perturbations affects the stationary state already in the leading order.

A. Weak and strong symmetries
The Hamiltonian H in Eq. (1) conserves the total number of electrons, i.e., [H, N PD + s=L,R N s ] = 0, with N PD and N s the number operator for the parallel dots and for the lead s, respectively. Since the leads feature no coherences between states with different numbers of electrons N s , any ρ PD that is diagonal in charge at the intial time will remains so at all times.
In the approximation of the Lindblad dynamics this symmetry is inherited as a weak symmetry of the Liouvillian with respect to N PD , that is, [49,50], see also Appendix A. Corresponding density matrices feature at most six non-zero entries in the basis of , which is the eigenbasis of H PD in Eq. (2) and will be referred to as the local basis. In that case, at most six modes contribute in Eq. (8).
We now discuss symmetries of the Hamiltonian and Lindblad dynamics originating from degenerate energies of the dots and their identical couplings to the two leads. Let us consider tunneling amplitudes such that and consider degenerate dot energies, In the local basis, the Hamiltonian in Eq. (1), which we denote by H (0) to indicate that the above conditions are fulfilled, remains the same when swapping the dot labels, up to the change of the sign for the doubly occupied state for the Coulomb interaction term. Thus, it is left invariant by the swap operator S, [H (0) , S] = 0, which exchanges the fermionic excitations between the dots, We have that S 2 = 1, so S is a parity operator. Indeed, the basis {|00 , . We refer to |+ and |− as the bonding and anti-bonding states, even though they remain degenerate here because of the absence of hybridization between the dots, and to the basis of |00 , |+ , |− , and |11 as bonding/anti-bonding basis.
For t jsk = t js assumed in the derivation of the Lindblad dynamics in Eq. (5), the condition in Eq. (11) can be expressed as the tunneling rates being equal [35] Γ js = Γ.
The Liouvillian inherits the symmetry of the Hamiltonian as a strong symmetry [49,50] with respect to S, i.e., the symmetries of the effective Hamiltonian, [S, H eff ] = 0, and, in contrast to a weak symmetry, additionally the jump operators, [S, J (0) αs ] = 0. Here, we used the index (0) to indicate that the conditions in Eq. (12) and (14) are fulfilled. Indeed, we obtain the effective Hamiltonian Here,B( ) = s=L,R B s ( )/2 is the average of the function B s ( ) that arises in the Lamb shift. This contribution lifts the degeneracy of the dot Hamiltonian caused by Eq. (12). Furthermore, the jump operators are given by

B. Strong symmetry implications for dynamics
In the bonding/anti-bonding basis for the dots, for parameters chosen as in Eqs. (11) and (12), the Hamiltonian separates into two distinct sectors which are not coupled by the tunneling processes (see [29,34,35]). Any initial state of the eigenspace corresponding to the eigenvalue 1 (or the eigenvalue −1) of the swap operator, that is, in the subspace of |00 and |+ (or the subspace of |− and |11 ), remains supported there at all times.
The dynamics preserves the eigenspaces of S due to the block-diagonal structures of the effective Hamiltonian in Eq. (15) and the jump operators in Eq. (16). Due to the strong symmetry, the parallel dots split into two independent two-dimensional systems in the subspaces of |00 and |+ , and of |− and |11 .
Indeed, J +s and J −s facilitate classical transitions between |00 00| and |+ +| at the respective rates 2Γf s ( ) and 2Γ[1 − f s ( )], and between |− −| and |11 11| at the respective rates 2Γf s ( + U ) and 2Γ[1 − f s ( + U )], while H eff does not contribute. In fact, these dynamics can be obtained as the Pauli rate dynamics in the bonding-antibonding basis. In contrast, a Pauli approach in the local basis (Appendix D) fails [35], in particular not capturing the degeneracy of stationary states.
The two stationary states in the eigenspaces of S are wheref ( ) = s=L,R f s ( )/2. While we have limited ourselves to considering initial states symmetric with respect to N PD , there are no other stationary states for U = 0. The two stationary states resemble the two degenerate ground states in the corresponding equilibrium system, which exhibits a quantum critical point in the zero-temperature limit [29].
The stationary currents from the left lead corresponding to the stationary states in Eq. (17) are Identical currents are found only if the Coulomb interaction vanishes (U = 0), or the Fermi distributions for the two leads are the same (T L = T R , and V B = 0). In the former case, the connected configurations feature the same energy difference, such that the system behaves as two identical independent single dots, see also Appendix B 1. In the latter case, the particle currents must asymptotically vanish as there is no directionality induced in the dot system with the leads at equilibrium. In the limit of infinite Coulomb interaction (U → ∞), f s ( + U ) → 0, and therefore I 2 → 0. The dynamics between |− −| and |11 11| corresponds then to a decay towards the lower energy state |− −|, so it becomes stationary while |11 11| is prohibited energetically [Eq. (17)]. But the asymptotic current I 1 remains unchanged as it is independent from U . The two stationary states in Eq. (17) fix the choice of the right eigenmatrices for two zero eigenvalues of the unperturbed Liouvillian. The corresponding left eigenmatrices are given by the projections on their supports, These determine the asymptotic state for a general initial state ρ PD (0) as p 1 ρ ss 1 + p 2 ρ ss 2 , where p j = Tr[P j ρ PD (0)], j = 1, 2, and thus the stationary current to p 1 I 1 + p 2 I 2 .
Next to the stationary states, there are two decay modes corresponding to the classical dynamics with the degenerate eigenvalues The other two modes describe the decay of coherences in the bonding/anti-bonding basis with the eigenvalues where the oscillation frequency arises from the Lamb shift. For the eigenmatrices, see Appendix B 1.

C. Breaking of strong symmetry and metastability
We now consider perturbations in the dynamical parameters that break the strong swap symmetry.
As a consequence, the two-fold degeneracy of the zeroeigenvalue of the Liouville operator is lifted, and a unique stationary state arises together with a new timescale in the dynamics for the system's final relaxation. Using non-Hermitian perturbation theory [45], we investigate those aspects of the dynamics with a focus on how the current is affected.
We examine perturbations that break the degeneracy of the dot energy levels as and for definiteness choose to alter the tunneling rates as In this work, we focus on small perturbations in the sense that δ , δΓ Γ, which also implies that δ T . This allows us to exploit non-Hermitian perturbation theory to characterize the eigenvalues and eigenmatrices of the perturbed Liouvilian, Above, we have expanded L in the perturbation parameters, where the superscript indicates the order of the perturbation. We consider corrections to L within the stationary state manifold of L (0) , which consists of probabilistic mixtures of the stationary states in Eq. (17). The focus of the following discussion is on the physical aspects; technical details can be found in Appendix B 2, see also Supplemental Materials in Refs. [9,11]. The first-order correction is necessarily zero, due to the effective classicality of the manifold of the stationary states of L (0) . The second-order correction corresponds to the classical dynamics of the probabilistic mixtures of the unperturbed system's stationary states in Eq. (17), where p 1,2 (0) = Tr[P 1,2 ρ PD (0)], so that p 1 (t) + p 2 (t) = 1 as P 1 + P 2 = 1 and the dynamics in Eq. (25) conserves the total probability. The decay rates γ 1,2 can be found by expressing the second-order corrections to the reduced dynamics in the basis of Eqs. (17) and (19), and they are quadratic functions of δ and δΓ. We now exploit this result in two ways.
The higher order corrections, indicated by ..., are of the first order. They can be understood to arise as the corrections to the structure of the two states corresponding to the stationary states in the unperturbed limit, which now constitute the two metastable phases. It is important to note that because γ 1 and γ 2 depend on two rather than a single perturbation, the stationary state depends on the perturbations already in lowest order, but only via their ratios, which can be seen as a competition between the metastable phases. This leads to distinct stationary current values when perturbations are varied, which will be discussed in more detail later, see also Section V. Indeed, the stationary current is where the corrections are at least of the second-order. Finally, since γ 1 = γ 2 in general, the stationary state will feature coherences in the local basis. Therefore, the approximation of the dynamics by a Pauli rate equation in that basis will fail, see also Appendix D.
Second, these results can be used to understand the long-time dynamics. For times longer than the initial relaxation, t −1/Re(λ 3,4 ), which in the leading order is determined by the eigenvalues λ (0) 3,4 of the unperturbed dynamics, the dynamics can be understood as follows.
The two-fold degeneracy of the zero-eigenvalue is now lifted with the second eigenvalue so that a new timescale emerges in the dynamics. This low-lying eigenvalue is accompanied by the right and left eigenmatrices where P 1 and P 2 are defined in Eq. (19) and corrections are at least of the first order. For times t such that t|δλ 2 | 1, with δλ 2 denoting the corrections in Eq. (28), which are at least fourth order, the state of the dots (see Eq. (9)) can be approximated as a probabilistic mixture of the states corresponding to stationary states in the unperturbed system, with the corrections of the first order. Access to later times can be gained by including higher order corrections in Eq. (25) (or modifying the generator even further [51]). The numerical methods introduced in Ref. [9,10] allow for the study of the long-time dynamics to all orders, simply by considering the low-lying part of the spectrum of the Liouvillian L, as obtained by its diagonalization; for a short summary, see Appendix B 3.
For times t such that −tλ 2 1, the effective dynamics in Eq. (25) can be neglected and the system is approximately stationary, i.e., metastable, with the leading contribution given by the asymptotic state of the unperturbed dynamics. In this metastable regime, the unperturbed system's stationary states take the role of metastable phases, and the system can be in any probabilistic mixture of these depending on the initial state [9,10]. As a consequence, a whole range of approximately constant average currents can be supported in this regime, In contrast to Eq. (27), the current values are in the leading order determined by the initial system state and thus independent from the perturbations.
We demonstrate our analytical findings of the parallel dots' dynamics for concrete parameter choices in Fig. 2. In Fig. 2(a), we characterize perturbation strengths for which the perturbative approach is applicable, with the metastability criterion [9,10] Re(λ 3 )/λ 2 1. In this regime, the stationary state is well approximated by the zeroth order terms in Eq. (26). To show this, the stationary probabilities of Eq. (25) are plotted as a function of the detuning δ in Fig. 2(b). They are compared with the decomposition into two metastable phases constructed numerically to all orders, see also Appendix B 3. Changing the dot energy perturbation while keeping a fixed tunneling rate perturbation, the distribution varies non-monotonically. In turn the stationary current is impacted, and can be suppressed when the metastable phase with the vanishing current (due to a large Coulomb interaction) predominantly contributes.
The transient dynamics in the regime where the system is metastable is qualitatively different to where the system does not exhibit metastability. In Fig. 2(c) the transient current calculated from the full dynamics [Eqs. (8) and (10)] is plotted for different choices of the detuning δ . For small perturbations, metastability occurs and the transient current remains approximately constant over a long time after the initial dynamics, see Eq. (32), before it finally evolves into its true stationary value, see also Eq. (27). That long-time evolution is well approximated by the classical effective dynamics of Eq. (25). For larger perturbations, the current evolves towards its stationary value continuously as the perturbative approach breaks down and to capture the dynamics correctly, the full expression for the evolution of Eq. (8) is needed.

D. Dynamics beyond small perturbations
Further insight into the system dynamics is provided by the eigenvalue spectrum of the Liouville operator in Fig. 3.
For increasing perturbation strengths, the difference between λ 2 and Re(λ 3 ) decreases, so there is no clear separation between the corresponding decay rates and the metastability is absent. Due to the Lamb shift, the eigenvalues λ 3,4 acquire an imaginary part which scales with the difference of the (renormalized) energy levels of the system, while λ 2 remains real. Disregarding the Lamb shift leads to a fundamental difference for δΓ = 0 and varying δ [ Fig. 3(a)]. In this case, the spectrum remains real below a threshold in δ , above which the eigenvalues λ 2 and λ 3 merge and their corresponding eigenvectors are identical, so that the spectrum exhibits an exceptional point [52]. For larger δ , the merged branches acquire an imaginary contribution forming a complex conjugate pair. The differences between the evolutions with and without the Lamb shift are less pronounced along constant δ but varying δΓ.
For large detuning δ /Γ 1, we observe the loss of coherence in the system as the real parts of the eigenvalues approach the values 0, −Γ, −3Γ, −4Γ. These correspond to the eigenvalues of the Pauli rate equation that describes the dynamics of a density matrix assumed diagonal in the local basis [37]; see also Appendix D.

V. QUANTUM COHERENCE ASSISTED SENSING
The coherences in the parallel dots lead to a stationary current which changes significantly as the perturbations δ and δΓ vary. From Sec. IV C, in the perturbative regime we understand this as the result of the competition of two metastable phases, but even for larger perturbations the current remains sensitive to changes in the perturbations [35]. We now investigate the possibility of using the parallel dots as a sensor when a parameter quench is detected through its effect on the particle current. In particular, we consider sensing a change in the nearby charge distribution which is assumed to lead to a shift in the perturbation of the dot energies δ . The change in charge distribution could, for example, be due to a single electron being added to or removed from another nearby quantum dot [53][54][55].
We consider a continuous measurement of the current in Eq. (10) during a time τ in order to detect the parameter quench that has occurred at t = 0. The statistics of such a measurement can be accessed using full counting statistics [36,47,48], see Appendix C.
The experimentally feasible measurement time τ sets a lower bound on the perturbation strength we consider. Indeed, the typical measurement time for charge detection in quantum dots is τ ∼ µs [56,57], while a typical value for tunneling rate sets Γ −1 ∼ ns. We assume that not only the initial but also the final relaxation takes place within the measurement time, and we therefore calculate all quantities in Fig. 4 in the stationary state. This means that the system cannot be too far into the regime where relaxation becomes extremely slow, which puts some lower bound on the perturbation strength, see Fig. 4.
In Fig. 4(a), the stationary current as a function of the two perturbations δ , δΓ clearly depends only on the perturbation ratio, but not visibly on the perturbation strength. For the regime where the system exhibits metastability this is expected in terms of the dependence of the current on p ss 1 , see Eq. (27). The signal of the (time-)integrated current is asymptotically linear in time with the rate equal to the response of the stationary current to a change in δ , which is shown in Fig. 4(b). In contrast to the current it depends on the perturbation strength and diverges as that is reduced. In the perturbative regime, this response is dominated by the change in the stationary probabilities, where we used ∂ δ p ss 1 = −∂ δ p ss 2 as p ss 1 + p ss 2 = 1. Indeed, ∂ δ γ j is of the first order while γ j is of the second order for j = 1, 2, so that ∂ δ I diverges with the inverse of the perturbations. This is true except for when ∂ δ p ss 1 = 0, which occurs in two cases. First, for the perturbation at δΓ = 0 since p ss 1 is then independent from δ . Second, when the perturbation ratio corresponds to the minimal current in Fig. 4(a). In those cases, the signal rate appears to actually vanish [as the corrections to the stationary current in Eq. (27) are of the second-order, the signal rate is then of the first-order].
The variance of the integrated current is asymptotically linear in time with the rate S(0) equal to the zero- frequency noise. The divergence of the fluctuation rate is evident in Fig. 4(c) and agrees with the behaviour of the Fano factor F = S(0)/I observed for δ , δΓ → 0 in Ref. [34]. Indeed, for smaller perturbations, metastability arises leading to long-lived correlations in the current so that its fluctuation rate becomes large [34,46] For fixed γ 1 + γ 2 , the smallest multiplicative factor corresponds to the stationary probability p ss 1 being minimal or maximal. Outside the parameter regime where metatability occurs, S(0) saturates to a constant value.
Estimation errors are determined via the standard error propagation formula as the measurement variance rescaled by the square of its signal. As the measurement time τ is assumed much longer than the relaxation time, the variance is dominated by τ S(0), and the signal by τ ∂ δ I. Therefore, the error is given by where the second line holds for the perturbative regime where metastability occurs. We see that the divergence of the signal rate in Eq. (33) and the fluctuations rate in Eq. (34) cancel out, and the difference in currents between the metastable phases simplifies as well.
For the best sensing setup, the errors in Eq. (35) should be minimised by the choice of the perturbation values before the quench. This corresponds to a trade-off between the fast divergence of the signal and the slow divergence of the fluctuations, which is non-trivial as, for fixed γ 1 + γ 2 , the slowest divergence of the fluctuations occurs exactly when the current is maximal or minimal and the divergence of the signal is actually absent.
We benchmark the performance of the parallel dots sensor with a single quantum dot setup for charge sensing [1,54,58], where a change in the charge distribution is assumed to affect the dot energy . The parameters for the single quantum dot are chosen such that it operates at a conductance peak where sensitivity is maximal, see Appendix D. The width of the conductance peak depends on the temperature where ∂ I ∝ 1/T , while S(0) is independent of T . Consequently the errors σ 2 show a quadratic temperature dependence. The advantage of the parallel dots system operating as a sensor is that it is not limited by temperature in contrast to the single dot setup. In Fig. 4(d) the temperature dependence of the error in Eq. (35) is shown for both setups. For the double quantum dot, the parameters are chosen close to the minimal value of the current (marker 1 in Fig. 4), where the signal rate is numerically observed to diverge fast, and small δ which results in the corresponding errors σ 2 δ remaining small over a wide range of temperatures. For example, at T /Γ ∼ 60 the ratio of error for single and parallel dots σ 2 /σ 2 δ ∼ 55. In the limit T /Γ → 0, the error of parallels dots approaches a constant value.
The parameter region in which the parallel dots can be operated as a sensor is not limited to the chosen set of parameters indicated in the stability diagram in Fig. 1(b). In fact, the trade-off between the noise and the signal anywhere between high conductance lines for positive bias voltages shows similar behaviour, see Appendix E.

VI. CONCLUSIONS
We have analyzed the transient dynamics and nonequilibrium transport properties of two interacting parallel quantum dots coupled to macroscopic leads.
A swap symmetry between the dots is present when the dots energies are degenerate and tunneling rates to both dots are identical. This parity-like symmetry leads to the existence of two stationary states distinguishable by their currents values. In equilibrium, this parity-like symme-try translates to a SU (2) symmetry of a pseudo-spin [30], and its diverging susceptibility indicates a quantum critical point [29]. Perturbations of dot energies and tunneling rates introduce metastability into the dots dynamics, and long-time dynamics towards a unique stationary state arises with a rate that is quadratic in the perturbations. Crucially, the stationary state depends already in the leading order on the ratio of the perturbations in dot energies and tunneling rates. This leads to the diverging signal (change in the stationary current in response to a small change in the perturbations). Since the current fluctuations also diverge in this limit, we found that the signal to noise ratio remains finite, but dependent on the perturbation ratios. In the context of charge sensing, a comparison with a single dot showed that the parallel dots may perform significantly better.
While the dynamics of the parallel dots was considered in this work as GKLS dynamics beyond the secular approximation, the discussed aspects were directly connected to those of the Hamiltonian dynamics of the dots and the leads, especially in the context of symmetry breaking. Therefore, our results should be qualitatively valid for other approximations for the dot dynamics [35].
The physics observed in this work crucially depends on coherent dynamics. It would therefore be interesting to investigate how the results change when including additional decoherence mechanisms, for example due to charge fluctuations in the environment. Additionally, how strong-correlation signatures in the conductance [27,28] are altered in the non-equilibrium setup, and how the corresponding current noise influences the sensitivity in a potential sensing application, remain open questions for future studies.

ACKNOWLEDGMENTS
We would like to thank R. Seoane Souto and V. Kashcheyevs for fruitful discussions. We acknowledge financial support from the Swedish Research Council (VR), the Knut and Alice Wallenberg Foundation (project 2016.0089), the Wallenberg Center for Quantum Technologies (WACQT), and from NanoLund. K.M. acknowledges the support from the Henslow Research Fellowship and, as a visiting researcher, from the Department of Physics, University of Cambridge.

Effective Hamiltonian
The dot Hamiltonian of Eq. (2), in its eigenbasis ordered as |00 , |10 , |01 and |11 , takes the following form, To find the effective Hamiltonian, we need to determine the Lamb shift H LS describing the renormalization of the system's eigenenergies due to the coupling to the leads. The latter is beyond the secular approximation given by [44] The operators X (a) , X (b) , a, b = 1, . . . , 4, correspond to the physical processes generated in the dots by the tunneling Hamiltonian in Eq. (4), i.e., Furthermore, l, m, n label the eigenbasis of H PD , while ω mn = E m − E n denote the corresponding energy differences. Finally, S αβ (ω) is determined by the odd Fourier transform, where D is the bandwidth. The function C ab (ω) = dτ C ab (τ )e iωτ , is in turn the Fourier transform of the lead correlation function Here, ρ L is the state of electrons in the leads, assumed to be given by the grand canonical ensemble and the evolution in the interaction picture [i.e., with the lead H L in Eq. In the continuum limit, with the assumption of tunneling amplitudes being independent from momentum k, the Fourier transforms of non-zero correlation functions are given by Here, ν is the density of states (assumed constant) and the Fermi distribution with the temperature T s and the chemical potential µ s for the lead s.
In the limit of a large bandwidth D → ∞, for S αβ (ω) we make use of the principal value integral [37,44,59] with the inverse temperature β s and the digamma function Ψ, to arrive at the Lamb shift Hamiltonian H LS with the following non-zero entries The Lindblad dynamics feature a weak symmetry with respect to the number of electrons N PD on the dots. From Eqs. (A1) and (A10), the effective Hamiltonian conserves the number of electrons, [N PD , H eff ] = 0.

Jump operators
To construct the jump operators for Eq. (5), we follow the approach presented in [42]. For the model of two parallel quantum dots we proceed as follows.

Each jump operator is identified with a physical
jump process between the quantum dot system and the leads, see Eq. (A3). For the parallel dots there are eight jump processes in total.
2. The processes for each lead are combined [42], which yields 3. Finally, each jump operator is reweighted in the eigenbasis of H PD of Eq. (2) according to the energy differences, by the corresponding distribution of the resonant energy in the lead before the electron exchange [cf. Eq. (A8)], and the corresponding density of states ν (assumed constant), where f s ( ) is the Fermi distribution for energy in lead s = L, R, so that in the basis of |00 , |10 , |01 and |11 .
The jump operators of Eq. (A13) increase or decrease the number of electrons only by 1, where α = +, − and s = L, R. In particular, Eq. (A14) leads to [N PD , J † αs J αs ] = 0, so that the average particle current in Eq. (10) is determined only by the components ρ PD (t) diagonal in charge. Therefore, the left and right eigenmatrices of the Liouville operator can be chosen as eigenmatrices of N PD , see also Appendix B 1. Here, we discuss further the dynamics in the presence of the strong swap symmetry, cf. Eqs. (15) and (16).
Next to the stationary states in Eq. (17) and the projections in Eq. (19), there are two decay modes corresponding to the classical dynamics, and with the degenerate pair of eigenvalues given by Eq. (20) There is also a decay of the quantum coherences in the bonding/anti-bonding basis, with the conjugate pair of eigenvalues in Eq. (21).

Perturbation theory for strong symmetry breaking
We now investigate dynamics of the Liouville operator L using the non-Hermitian perturbation theory with respect to perturbations away from dynamics featuring the strong swap symmetry, see Eq. (24). Those arise due to perturbations in the effective Hamiltonian and the jump operators, cf. Eqs. (15) and (16), when dynamical parameters are changed according to Eqs. (22) and (23).

a. Perturbations of Liouvillian
The resulting perturbations to the Liouville operator caused by perturbations of the dynamical parameters are of all orders. This is due to the fact that the effective Hamiltonian and the jump operators are non-linear functions of the dot energies and the tunneling rates, see Appendix A. In particular, we have that the first-order perturbation of the Liouvillian [cf. Eqs. (5) and (24)] is given by which stems from the first-order perturbations to the effective Hamiltonian in Eq. (15) and the jump operators in Eq. (16). Similarly, the second-order perturbation of L is given by αs , (B6) where both the first-and second-order perturbations to the effective Hamiltonian and the jump operators contribute. We will denote the contribution from the first-order perturbations only as Below, we only give first-order perturbations to the effective Hamiltonian and the jump operators, as they fully determine the leading second-order corrections to the long-time dynamics, via L (1) and L (2) , which is argued in Appendix B 2 b.
For the dot Hamiltonian of Eq. (2) the perturbation is linear in δ We note that the effective Hamiltonian and the jump operators feature symmetry-breaking perturbations only in the first order. In fact, it can be shown that symmetrybreaking perturbations of those operators appear in odd orders, while symmetry-preserving perturbations appear in even orders. This is a consequence of the fact that choosing perturbations with the opposite signs in Eqs. (22) and (23), directly corresponds to the dynamics with the dots swapped. Under this transformation, |− is replaced by −|− and |11 by −|11 in the bonding/antibonding basis. Since the simultaneous change of all perturbation signs changes the sign of odd order corrections, those must correspond to the symmetry-breaking contributions, while even orders must be accompanied only by symmetry-preserving contributions.
b. Perturbative corrections to dynamics As the dynamics, no matter the size of perturbation, preserve the weak symmetry with respect to N PD , the only unperturbed modes that contribute to the perturbed dynamics ofρ PD , as considered in the main text, are the modes diagonal in charge [see Section IV A and Appendix A].
We consider the perturbation theory for the reduced dynamics of the first two eigenmodes of the dynamics, LP, where P(ρ PD ) = ρ ss + Tr(L 2 ρ PD )R 2 [cf. Eqs. (8) and (9)]. We have that P = P (0) +P (1) +... where P (0) is the projection on the zero-eigenspace of the unperturbed dynamics L (0) , The first-order corrections to the reduced dynamics are always within that subspace and are formally given by P (0) L (1) P (0) [45]. We now show that these corrections are zero for the dynamics considered in this work. This can be seen as the consequence of the classicality of the zero-eigenspace of the unperturbed dynamics (cf. Ref. [11]).
In the first order, the perturbations of the effective Hamiltonian and the jump operators are symmetry-breaking and give rise to the first-order corrections to the Liouvillian L (1) , given in Eq. (B5), which break the strong symmetry. That is, L (1) (ρ ss i ) is a linear combination of coherences |+ −| and |− +|. Since the coherences decay to 0 under the unperturbed dynamics, Tr(P i |+ −|) = 0 = Tr(P i |− +|), the first-order corrections vanish, P (0) L (1) P (0) = 0.
For P (0) L (1) P (0) = 0, the second-order corrections to the reduced dynamics, are found within the zerosubspace of L (0) and formally given by P (0) L (2) P (0) − P (0) L (1) R (0) L (1) P (0) [45], where R (0) is the reduced resolvent of L (0) at 0, so that R (0) L (0) = L (0) R (0) = I −P (0) , with the identity map I(ρ P D ) = ρ P D . In terms of the eigenmatrices of L (0) we can write [Eqs. (20), (21), (B1) and (B3)] We now show that the contribution to the second-order dynamics stem only from the first-order corrections to the effective Hamiltonian and the jump operators. Indeed, let us note that L (2) − L (2) is of an analogous form to L (1) but with the first-order perturbations H We conclude that the second-order corrections to the reduced dynamics are given by P (0) L (2) P (0) − P (0) L (1) R (0) L (1) P (0) . We now use this results to calculate the decay rates in Eq. (25). In the operator basis of |+ −|, |− +|, |00 00|, |+ +|, |− −|, |11 11| (which we index by I, II, III, IV, V, VI, respectively), the firstorder perturbation L (1) of the Liouvillian has the following structurê Here, the complex conjugation relations between the first and second columns, and between the first and second rows, follow from the fact that L (1) is Hermiticity preserving. We also haveL V,II ,L IV,II ,L V,I , and IV,I , as those contributions arise only due to the effective Hamiltonian. The structure of the second-order perturbation L (2) due to the first-order perturbations of the jump operators, can be understood as corresponding to the strong symmetry with |00 00| − |+ +| + |− −| − |11 11|, with which the first-order perturbations of the jump operators commute. We then use the trace-preservation of L (2) to connect the diagonal terms to the off-diagonal ones, and its Hermiticity-preservation to note that [L (2) ] II,II = [L (2) ] * I,I . Therefore, the decay rates in Eq. (25) are whereL (1) using Eqs. (B8) to (B10).
Third-order corrections to the reduced dynamics projected on the zero subspace of L (0) in general contribute to the first-order perturbations of the stationary state in Eq. (26) and the eigenmatrices in Eq. (28) corresponding to the second eigenvalue. But here the third-order corrections (see Supplemental Material of Ref. [9] and [45]) are projected to 0 as, when acting on ρ ss 1 and ρ ss 2 , they give rise to linear combinations of |+ −| and |− +|. In particular, we have, λ 2 = 0. These results can be seen to hold for all odd-order corrections, as the opposite sign of the perturbations in Eqs. (22) and (23) corresponds to the dynamics with the dots swapped, which however leaves ρ ss 1 and ρ ss 2 unchanged.

c. Perturbative corrections to metastable and stationary states
Beyond the zero-subspace spanned by the stationary states of L (0) in Eq. (17), the first-order corrections to the projection P on the stationary state and the second eigenmode are given by P (1) = −R (0) L (1) P (0) [45]. These corrections determine the first-order corrections to the metastable phases as −R (0) L (1) (ρ ss 1 ) and −R (0) L (1) (ρ ss 2 ); see Supplemental Material of Ref. [9]. We now assess the first-order corrections to the unique stationary state for L, cf. Eq. (24). Again, due to the weak symmetry with respect to N PD , the only unperturbed modes that contribute are the symmetric ones, cf. Appendix B 1. The first-order corrections to the stationary state in Eq. (26) are given by ρ ss(1) Here, the projected third-order corrections to the reduced dynamics should also contribute (cf. Supplemental Material of Ref. [9] and [45]), but they vanish for the considered perturbations, as we explained above.

d. Perturbative corrections to initial dynamics
Generally, the dynamics taking place before the metastable regime can be analyzed in terms of the perturbative corrections to the remaining fast eigenmodes. The leading corrections for the eigenvalues are of second order, but for the eigematrices of first order.
The eigenvalues corresponding to the coherence decay in Eq. (21) will acquire second-order corrections as in the first order λ The degeneracy of the classical decay of eigenvalues in Eq. (20) will only be lifted in the second order, with the first-order corrections being zero (cf. Fig. 3). The first order corrections to the corresponding matrices in Eqs. (B1) and (B2) will be present.

Metastable phases and long-time dynamics -all orders
Formally, the long-time dynamics in Eq. (9) is a projection onto the subspace of the first two eigenmatrices of the Liouvillian. Here, we review the construction first introduced in [9] that allows for considering it in a physical basis, and thus considering Eq. (25) not only up the second, but to all orders, as in Fig. 2(b).

a. Metastable phases
To explain how to construct the metastable states in terms of the first two eigenmodes of the dynamics, as used in Fig. 2(b), we follow [9,10].
The metastable manifold is spanned by two extreme metastable statesρ 1 = ρ ss PD + c max 2 R 2 , ρ 2 = ρ ss PD + c min 2 R 2 . (B19) The coefficients c min 2 and c max 2 are the smallest and largest eigenvalues of the left eigenmatrix L 2 . Using the results of Appendix B 2, up to the first-order corrections, we obtainρ 1 = ρ ss 1 + ... andρ 2 = ρ ss 2 + .... The approximation in Eq. (9) for any state during the metastable regime corresponds to the projection P on the stationary state and the second eigenmode, and can be equivalently expressed as a linear combination withp i (t) defined via the observables P 1 = L 2 − c min 2 1 /∆c 2 , P 2 = (−L 2 + c max 2 1) /∆c 2 , with ∆c 2 = c max 2 − c min 2 asp i (t) = Tr[P i ρ PD (t)] for i = 1, 2. The metastable phasesρ i defined above feature trace 1, but are in general not positive and thus are not described by density matrices. In contrast,p i (t) always correspond to probabilities [9,10]. Up to the first-order corrections, we haveP i = P i + ..., so that p i (0) = p i (0) + .... p ss 1 p ss 2 (I s1 − I s2 ) 2 + ..., so that the fluctuation rate diverges inversely with the square of perturbation strength [11], see also Eqs. (26) and (29).

Appendix D: Pauli rate equation for parallel dots
Here, we compare the Pauli rate equation from Eq. (5) and consider the resulting stationary distributions. A Pauli rate equation for the diagonal entries of the density matrix in any basis can be obtained by neglecting the contribution from coherences.
In the basis |00 , |10 , |01 , and |11 , the Pauli rate equation [37] features a single stationary probability distribution. Even for the parameters chosen as in Eqs. (12) and (14), which lead to stationary state degeneracy in the Lindblad dynamics of Eq. (5), the distribution remains unique.
In Fig. 5(a) and (b), the stationary differential conductance and stationary current for the rate equation are plotted [37]. The stability diagrams significantly differ from those for the stationary state of Eq. (5), cf. Fig. 1(b) and Fig. 6(a). Here, the differential conductance recovers the typical Coulomb diamond structure for a single spinful dot. In fact, when 1 = 2 = and Γ 1s = Γ 2s = Γ s , additionally with U , the dynamics corresponds to a single spin-degenerate dot, where coherences are suppressed because spin is a good quantum number in both leads and dot. Such systems are used as charge sensors [1], where the parameters are chosen along the high conductance lines, where small changes in the gate voltage V G results in large response in the current.
In the lowest order, classical dynamics arises between |00 00|, |01 10|, |10 10|, and |11 11| that are left invariant by the Hamiltonian, and is given by the Pauli rate equation in the local basis. The real parts of the eigenvalues for the rapidly oscillating coherences |10 01| and |01 10| are both given by − j=1,2 s=L,R Γ js [1 − f s ( j ) + f s ( j + U )]/2. The corresponding eigenvalue spectrum in Fig. 5(c) indicate such behavior.  Fig. 1(b) as functions of VB and VG. The structure of the error in (d) within the Coulomb diamond is due to the current, noise and sensitivity not being exactly zero, but exponentially suppressed in this region. The star indicates the position in the stability diagram used for Fig. 4.