Dynamical signatures of molecular symmetries in nonequilibrium quantum transport

Symmetries play a crucial role in ubiquitous systems found in Nature. In this work, we propose an elegant approach to detect symmetries by measuring quantum currents. Our detection scheme relies on initiating the system in an anti-symmetric initial condition, with respect to the symmetric sites, and using a probe that acts like a local noise. Depending on the position of the probe the currents exhibit unique signatures such as a quasi-stationary plateau indicating the presence of metastability and multi-exponential decays in case of multiple symmetries. The signatures are sensitive to the characteristics of the probe and vanish completely when the timescale of the coherent system dynamics is much longer than the timescale of the probe. These results are demonstrated using a 4-site model and an archetypal example of the para-benzene ring and are shown to be robust under a weak disorder.

Symmetries are a vital ingredient to many molecular systems and play an important role in their transport properties [1][2][3][4].They lead to intriguing non-trivial consequences like the existence of multiple steady states [5,6] and could help design smart devices [7,8].Despite these advantages, the effect of symmetries in the transport properties of molecular systems have not been fully explored.Several studies on symmetric systems have explored the effects of dephasing [2,10] and interference [11,12].
A thorough knowledge of symmetries could lead to a better understanding of molecular junctions that have become promising candidates to build transport devices [13,14].They have been in the limelight due to the sizable quantum effects that have led to exciting effects in electronic [15][16][17], heat [18][19][20][21], and excitonic [22,23] transport.Despite the plethora of studies, transport signatures that arise solely due to the inherent molecular symmetry have not been put forward.These signatures could help identify the symmetric sites within the molecule thus providing a clear atomic picture of the system of interest.Using these signatures, we could explore molecular symmetry to design nanodevices with tunable transport properties.
In this study we propose an efficient approach to detect symmetries in molecular junctions by measuring the excitonic currents.The main idea relies on using a probe that could act locally on each site [24] thus breaking the symmetry of the system.Therefore, by varying the position of the probe and initiating the system in an appropriately chosen initial condition we find the excitonic currents vary from zero (for symmetric system) to a finite non-zero value (for systems with broken symmetry).The non-zero currents exhibit a long quasi-stationary plateau akin to * jianshu@mit.edua classical metastable state [25].If there are multiple symmetries present within a system, like in the case of para-benzene, then the multiple symmetries can also be distinguished due to the presence of a multi-exponential decay in the current dynamics.The signatures we obtain are exclusive to symmetric systems and are found to be sensitive to the timescales of the system and the probe.Thus, we provide a comprehensive detection scheme of molecular symmetries and elucidate its dynamical effects.

Results
Detecting symmetries.Symmetry detection forms the first crucial step to manipulate molecular junctions in order to build novel nanodevices, such as thermal switches [8].One of the key properties of symmetric systems is the presence of multiple steady states.These occur when the reduced dynamics is confined to an invariant subspace that is not disturbed by the dissipative leads (see methods section for a general proof).Our proof shows that when the interaction of the reduced system doesnt perturb the system out of the invariant subspace, the symmetries are preserved and reflected in the reduced dynamics in the form of multiple steady states.In other words irrespective of the coupling strength or the form of the master equation that governs the reduced dynamics a symmetric system will possess multiple steady states.
Detecting such symmetries is highly non-trivial, especially in experimental set-ups [26], and in this work we propose a scheme that would help achieve this objective and provide robust signatures in the transient currents.Our scheme consists of two steps: i) Initiating the molecular system in a dark state and ii) Using a probe that acts locally on the molecular sites to detect which sites are symmetric.The probe when placed on different sites could break or preserve the symmetries of the system.The symmetry breaking would produce a non-zero current whereas when symmetries are preserved the current would be zero (due to the dark state initial condition that belongs to an invariant subspace).
For perfectly symmetric systems the switching phenomenon can be observed even in the steady state.But perfect symmetry is difficult to obtain in experimental set-ups and most systems are subject to conformational disorder.Even for a 'close to' symmetric system with a weak disorder (either static or dephasing) there would be a single unique steady state.Hence the steady-state currents will not switch depending on the probe position.In this case the dynamical regime could provide essential signatures.If the disorder is weak it would affect the dynamical currents on the longest timescale and thus during the intermediate times we expect that a weakly disordered system would mimic a perfectly symmetric system.
The breaking of symmetries could also affect the relaxation times and help unravel the multiplicity of the steady-state degeneracy.Thus, our goal of detecting symmetries could be achieved using the two step approach that could even help design a control over the timescale of the relaxation process.Inspired by these ideas and the timely application of detecting symmetries we explore the points presented in this section using two concrete examples of a 4-site model and the archetypal benzene molecule.4-site model.The minimal nonequilibrium model that exhibits molecular symmetries is a 4-site model as sketched in Figure 1.The Hamiltonian of the system is given by, where the 4-site Hamiltonian H ′ S reads In the above the sum over i, j indicates nearest neighbor interactions with cyclic boundary conditions.This Hamiltonian within the Hückle theory of atomic molecular orbitals can be interpreted as having 4 atomic sites with nearest neighbor interaction [27,28].The diagonal elements ε represent the combination of kinetic energy and Coulomb integrals whereas the off-diagonal coefficient h represents the nearest-neighbor overlap between atomic sites.The system is connected to left and right leads that govern the transitions of excitons between the ground state e g and the excited states e 1 and e 3 .The presence of two leads at different temperatures creates a nonequilibrium setup that generates a flow of excitons within the molecular system.The effect of the leads is modeled via a Lindblad term in the master equation (see methods section) whose Lindblad operators read The nonequilibrium 4-site model exhibits two nonequilibrium steady states [1,5,7]  Our goal in this work is to devise a scheme to detect molecular symmetries using the excitonic currents.In order to achieve this objective we introduce a probe that acts on the system as a site specific dynamical noise without injecting or extracting excitons.When the probe acts on sites 1 (S = |e 1 e 1 |) or 3 (S = |e 3 e 3 |) the Hamiltonian of the probe-system interaction commutes with the original symmetry operator (see methods section).On the other hand the action of the probe on site 2 (S = |e 2 e 2 |) or 4 (S = |e 4 e 4 |) breaks the symmetry of the system leading to a single unique nonequilibrium steady state.We incorporate the effects of the probe via a Redfield tensor as shown in the methods section.Similar probes, e.g. the Büttiker probe [30], have been used to investigate effects of coherence [30], interference [12] and also to measure temperatures at a molecular level [31].
The behavior of excitonic currents for the 4-site model starting from different initial conditions and different probe positions is shown in Fig. 1.The probe is acting either on site 1 (unbroken symmetry case) or on site 2 (broken symmetry case).The three initial conditions correspond to the symmetric exchange of sites 2 and 4, ρ s = 0. Figure 1a shows the steady-state currents as a function of the mixing angle (θ) between the symmetric and antisymmetric initial conditions.We set the initial condition ρ(0) = cos 2 (θ)ρ s + sin 2 (θ)ρ a . ( For the perfectly symmetric system when the probe is placed at site 1 (unbroken symmetry case, black solid line) the excitonic currents can be controlled by tuning the mixing angle θ.On the other hand when the probe is placed at site 2 (broken symmetry case, green dashed-dotted line) we obtain a unique steady state and hence there is no modulation of the steady state currents due to the mixing angle.Thus, in case of perfectly symmetric systems the steady state currents exhibit a clear signature of the underlying symmetry in terms of the switching behavior depending on the probe position.Unfortunately, perfectly symmetric systems are rarely feasible and in order to describe a realistic system we introduce a static random disorder ∆ε i at each site e i (∀i = g) chosen from a uniform distribution with width . The excitonic currents are averaged over 1000 samples (enough to achieve convergence) of disorder to obtain the closed circles in Fig. 1a.Clearly, in this case the probe position (site 1 ≡ black closed circles and site 2 ≡ green closed circles) plays no role since we have only one unique steady state.Therefore, in a realistic system with conformational disorder it is impossible to observe any steady-state signature of the underlying symmetry.Hence, this observation serves as our main motivation to study the dynamics of the excitonic currents, due to the separation of time scale of the disorder and the symmetries, in order to observe the dynamical signatures related to symmetries.
In order to achieve our objective we first study the transient currents in perfectly symmetric systems.We express the time evolution of the reduced density matrix in terms of an eigenfunction expansion [32] of the dissipative Liouvillian L (see Eq. ( 19) from methods section) as, where we are working in the Liouville-Fock space [33], mapping the density matrix to a vector and the Liouvillian L to a matrix L that follows Here Φ l | and |Φ r are the left and right eigenoperators of the dissipative non-Hermitian Liouvillian that form a dual basis, and Λ k are their corresponding eigenvalues.The dynamics then clearly depends on two important factors, namely, the eigenvalues Λ k and the weight determined by the overlap of the corresponding left eigenoperator with the initial density matrix Φ l k |ρ(0) .The eigenspectrum for the dissipative Liouvillian for our 4-site model is shown in Fig. 2. The probe, when placed at site 2, produces a subtle effect by splitting the degenerated zero eigenvalues as indicated in Fig. 2b.The eigenvalue that separates from the zero (indicated by Λ 1 ) corresponds to the unstable steady state and has an antisymmetric structure.Thus, its left eigenoperator overlap with the symmetric initial condition is zero and the dynamics of the excitonic currents is governed by the remaining eigenvalues and eigenoperators that are not drastically affected due to the probe position.Therefore, we find that for the symmetric initial condition there are no dynamical signatures of symmetries in the excitonic currents as shown in Fig. 1b.
The antisymmetric initial condition (also known as the dark state) [34] is one of the multiple steady states (unstable manifold).In this state the excitonic currents are exactly zero as shown by the black solid line in Fig. 1c.As soon as we break the symmetry we observe a slow relaxation for the excitonic currents due to the transition from the unstable to the unique steady state.The timescale of this transition is now governed by the dissipative Liouvillian eigenvalue corresponding to the unstable steady state τ R ∝ Λ −1 1 which is several magnitudes larger than the relaxation time for the symmetric initial condition (note the x-axis in Figs.1b and c).We also find that there is no change in the time-dependent excitonic currents if the probe is placed at site 2 or 4. Thus, the slow quasi-stationary relaxation, when initiating from a dark state, provides a clear signature of symmetries in the dynamical currents that inherently depend on the probe position.
The canonical initial condition has overlap with both the stable and unstable steady states.The unstable steady state does not contribute to the dynamics because its relative overlap with the canonical initial condition is much smaller than the overlap with the stable steady state.Hence, the mechanism leading to dynamics that is independent of the probe position is very similar to the symmetric initial condition case.Thus, the canonical initial condition, that is extensively used to study transport using nonequilibrium Greens function [35], hierarchy equation of motion [12], and polaron transformation [20], makes it impossible to detect molecular symmetries.
Our results described above are robust even in presence of a weak disorder as shown by closed circles in Fig. 1c.If the disorder is weak (much weaker than the effect of the probe, i.e., δ 2 ≪ γ), then the eigenspectrum of the dissipative Liouvillian (Fig. 2 asterisk), except the eigenvalue corresponding to the steady state, is not affected by the disorder.Since the weak disorder acts on the longest timescale τ dis ∝ δ −2 of the problem, we ob-serve the dynamics at an intermediate timescale by ignoring the longer times that are affected by the disorder.In this intermediate timescale regime the disorder is insignificant and the dynamics is similar to the perfectly symmetric scenario (black solid and green dashed-dotted line in Fig. 1c).Hence, a priori if it is known that a weak disorder affects the system of interest, its effects could be neglected by observing the current dynamics at an intermediate timescale.
Next we investigate the effects of probe properties on the dynamics when initiated from the antisymmetric state and the probe is placed at site 2. Due to the broken symmetry the unstable steady state forms a Hermitian decay mode, i.e., right eigenoperator of the dissipative Liouvillian with non-zero real eigenvalue [36], that remain closest to the zero eigenvalue (unique steady state).Since this mode is formed due to the presence of the probe its eigenvalue is sensitive to the probe properties thus influencing the relaxation time [37].
Although this argument seems straightforward there is a caveat in terms of the dominant timescales of the problem.The largest frequency modes of the probe corresponding to the cut-off frequency ω D , Eq. ( 20), dictate a timescale τ P ∝ ω −1 D [38], herein termed as the probe timescale.On the other hand the coherent system dynamics timescale is given by τ S ∝ ∆ where ∆ is the smallest finite energy difference in the system eigenspectrum.These timescales play an influential role in determining the relaxation time of the system τ R .Figure 3 shows the effect of interplay between these two timescales.When the system dynamics timescale dominates over the probe τ S ≫ τ P , the probe has very little influence on the relaxation time τ M .The probe merely creates a unique steady state and a change in probe parameters, e.g.temperature, does not strongly influence the eigenvalue of the unstable steady state Λ 1 .Thus, as a result the system relaxes quickly to the unique steady state and a change in probe temperature shows no effect as seen in the inset of Fig. 3a.On the other hand when the probe timescale dominates the system τ P ≫ τ S the probe can strongly influence the relaxation time τ R .In this case a change in the parameters, e.g.probe temperature, causes a proportional shift of the unstable steady state along the real axis thus influencing the relaxation time (inset of Fig. 3a).Since the shift is proportional to the temperature of the probe it is expected that the eigenvalue Λ 1 would be very close to zero for small values of temperature.In other words when τ P ≫ τ S the reorganization energy, γ reorg = ∞ 0 J(ω)/(πω) ≡ γω D with J(ω) being the spectral density defined in Eq. ( 20), is extremely small and the probe acts like a weak perturbation to the system, whereas when τ S ≫ τ P the probe acts as a strong perturbation wiping out the symmetry information.Therefore, initiating from the antisymmetric state at low temperatures leads to a slow quasi-stationary relaxation towards the unique nonequilibrium steady state due to the strong interaction of the initial condition with the unstable steady state.The solid curve in the inset is for the case when τP ≫ τS, whereas the dotted line is for τP ≪ τS.In both panels the system is initiated in an antisymmetric state.All other parameters are the same as in Fig. 1.
Thus, overall by changing the probe position and carefully engineering the initial state to a dark state [7,8,36] we detect the signatures of molecular symmetries.These signatures correspond to observing exactly zero excitonic currents for symmetric systems and a finite timedependent current with a long quasi-stationary plateau for systems with broken symmetry via the probe.Interestingly, the signatures completely vanish for the symmetric and canonical initial conditions making these an unacceptable choice for symmetry detection.

Benzene molecule. The Hückle theory Hamiltonian H ′
S for benzene is given by wherein we choose the parameters ε = −11.2eVand h 1 = −0.7eV[39].The benzene molecule is connected to Lindblad leads similar to Eqs.The steady states that depend only on the symmetries can be obtained analytically and they are given by wherein |ψ i are eigenvectors of the benzene Hamiltonian.The nonequilibrium steady state ρ 0 depends on all the lead properties and is nontrivial to obtain.In the presence of a probe ρ 1 and ρ 2 become unstable.In order to understand how the probe affects the multiple steady states we project the full dissipative Liouvillian, Eq. ( 19), into the subspace of the multiple steady states.This is achieved by transforming the steady states into vectors and the dissipative Liouvillian into a matrix as done in Eq. ( 6).Thus, the 3 × 3 perturbation matrix reads where the elements θ n (n = 0, 1, 2) contain information about the nonequilibrium steady state and the probes and σ k ∝ j =1,2 R jk .The matrix elements R jk are the elements of the Redfield tensor R[ρ] = ∞ 0 dt[S, ρS(t)]C(t) + h.c.(see methods section Eq. ( 19) and also Refs.[40,41] for more details) and are given by Above we have taken into consideration that several independent probes could act on the system whose actions are taken into account via the summed index p.
The operator elements are Here E k is the eigenenergy of the benzene Hamiltonian corresponding to the state |ψ k .The perturbation matrix δL has one zero eigenvalue corresponding to the unique nonequilibrium steady state and two non-zero eigenvalues that dictate the relaxation time.The structure of the matrix δL clearly indicates that the elements R jk govern the effect of the probe on the symmetries of the system.Since the probe obeys detailed balance the backward transition rates Re[W p 12 ] are extremely small.Thus, we can neglect the matrix element R 12 ≈ 0 and simplify the perturbation matrix.This simplification occurs only due to the probe properties and does not depend on how the probe connects to the system.
The matrix element R 21 controls the structure of the eigenvalues of the perturbation matrix and depends only on how the probe connects to the system (via operator S).In particular, we find that when R 21 = 0 the nonzero eigenvalues of δL form a complex conjugate pair.On the other hand when R 21 = 0 we break all the symmetries and obtain three distinct eigenvalues of the perturbation matrix.In order to validate our claim we take the simplest example of a single probe p = 1 that breaks the symmetry of the system and connects via a general operator S = 6 l=g,1 C l |e l e l | with C l being constants.Thus, in order to achieve R 21 = 0 we require If C 2 = C 3 = C 5 = C 6 then the probe is affecting all the symmetric sites of benzene equally without breaking any symmetry.Thus, even though Eq. ( 12) is satisfied we obtain 3 multiple steady-states.Hence it is of utmost importance that the probe breaks at least one of the symmetries of benzene.If we choose a local probe acting on site 2, i.e., C 3 = C 5 = C 6 = 0, it breaks some but not all the symmetries [since Eq. ( 12) is not satisfied].This case leads to a unique steady state even if not all the symmetries are broken.The eigenvalues of the perturbation matrix (black plus) and the full dissipative Liouvillian (green triangles) are plotted in Fig. 4a for this case.Due to the fact that we are not breaking all the symmetries, we obtain a couple of non-zero eigenvalues that are conjugated and there is only one relaxation time (given by the real part of the non-zero eigenvalue).The dynamical currents are unable to show any exclusive signatures of the three steady states.The current dynamics are qualitatively similar to the 4-site model (see supplementary material).In order to break all the symmetries we take the simplest case C 2 = C 3 and C 5 = C 6 = 0.In this case the probe acts non locally on sites 2 and 3 [Eq.( 12) is satisfied].This leads to three distinct real eigenvalues of the perturbation matrix as shown in Fig. 4b.The left-eigenvectors corresponding to the distinct non-zero eigenvalues inherit the symmetry properties of the unstable steady states.Particularly the left-eigenvector corresponding to the eigenvalue closest to zero retains the symmetry of ρ 1 whereas the one further away transforms into a linear combination of ρ 1 and ρ 2 .Thus, we initiate the system in a state that is antisymmetric in sites 2 and 5 (equivalent to ρ 1 ) and find that at long times both the non-zero eigenvalues closest to the steady state (Λ 1 and Λ 2 ) will contribute to the relaxation time.Hence, due to the remnants of the unstable steady states in the Hermitian decay modes we can uncover the signatures of the multiple steady state as seen in the current dynamics in Fig. 5.At long times the current dynamics exhibits two distinct exponential relaxation times τ 1 ∝ Λ −1 1 and τ 2 ∝ Λ −1 2 .The corresponding coefficients of the exponential decays are are opposite in sign and hence the excitonic current shows a peak indicating the presence of two unstable steady states.Here once again similar to the 4-site model a weak static disorder in the system, as shown by the closed green circles in Fig. 5, plays no role at the relevant timescale.Thus, by redesigning the probe to act non locally it is possible to detect the number of multiple steady states by observing the multi-exponential decay at long times.

Discussion and Conclusions
Several molecular systems possess symmetries that could play a crucial role in the nonequilibrium properties of the quantum devices.Naturally, one wonders if it is possible to detect symmetries without any prior knowledge of the molecular system.In the present study, we provide a toolbox that detects molecular symmetries via the measurement of excitonic currents.
Our goal is achieved by introducing a probe that acts on various sites of the molecular junction and depending on the probe position the currents differ suggesting the presence of symmetries.In particular, we find signatures of molecular symmetries if we start from an antisymmetric (dark) state as our initial condition.If the probe is then positioned on a site that preserves the molecular symmetry we find zero excitonic currents because the dark state belongs to an invariant subspace of the nonequilibrium set-up.On the other hand if the probe is placed on one of the symmetric sites it breaks the inherent symmetry of the system leading to a quasi-stationary plateau in the excitonic currents.When the probe is placed on any of the symmetric sites we see identical time-dependence of the excitonic current indicating that these sites are interchangeable, i.e., symmetric.These signatures vanish if we initiate the system in the symmetric or canonical initial condition making the antisymmetric state a primary choice for symmetry detection.
Our detection scheme is based on dynamical signatures in transient currents and is robust against the weak disorder.We specifically show that in the presence of disorder the steady-state currents do not show any symmetry related signatures, whereas the transients exhibit dynamical signatures related to the the perfectly symmetric systems.The strength of the disorder must be weaker than that of the probe, such that the effect of disorder does not appear except for the longest timescale.Thus, by ignoring the longest timescale τ dis ∝ δ −2 the dynamics at intermediate timescales becomes identical to that without disorder.In other words, even for a weakly disordered symmetric system we could detect the underlying symmetries using the probe, since the dynamics at intermediate timescale mimics the long-time dynamics of the perfectly symmetric system.
Since the antisymmetric state is our prime choice for symmetry detection we also study the effect of probe properties on the dynamics, if we initiate the system from this state.We find that the interplay between the timescales of the probe and the system dynamics play a significant role in determining the dynamics and the relaxation time of the system.In particular, when the molecular symmetry is broken, the system possess a long relaxation time that scales with temperature if the timescale of the probe dynamics is much larger than the system (τ P ≫ τ S ).In the opposite case temperature shows no effect on the relaxation time indicating that the probe does not affect the unstable steady state leading to quick relaxation.The long relaxation time indicates the presence of a long-lived quantum quasi-stationary distribution [42] that survives due to its interaction with the unstable steady state.
In case of the symmetric para-benzene molecule we found three nonequilibrium steady states.In order to break all the symmetries we analyzed the system by treating the probe as a perturbation and showed that the multiple symmetries can be broken via a redesigned probe.Our new probe acted on two sites of the parabenzene ring simultaneously and led to a bi-exponential decay in the current dynamics.The bi-exponential decay exhibited long relaxation times and the exponents corresponded to the eigenvalues of the unstable manifold.This trend of multi-exponent decays corresponding to the number of symmetries present was then recognized as the second key signature to detect symmetries.
Overall, our work provides a first glimpse into symmetry detection in nonequilibrium transport set-ups and could be extended to systems such as conjugated dendrimers [43], light-harvesting complexes [44], C60 bucky balls [45], etc. to detect naturally occurring symmetries.Also since quantum coherence plays a key role in the performance of light-harvesting systems and quantum heat engines [46,47], it naturally leads to the speculation that structural symmetry studied here will lead to interesting quantum effects in these systems Methods Invariant subspaces.It has been previously proved that a Lindblad master equation has invariant subspaces if there is a unitary operator that commutes with all the master equation elements [1,5].In this section we extend this proof to a general open system.We begin with a microscopic description and decompose the total Hamiltonian as H = H S + H B + H BS .The Hilbert space of the total system is decomposed in the direct product of our system of interest (S), and the bath (B),

and
The bath B in our description could comprise of several leads and probes, e.g., H B = H L + H R + H P , where the sub-scripts L, R, and P correspond to the left-lead, right-lead, and the probe.
We assume that there is a unitary operator Π ∈ O H S s.t.
[Π, H S ] = [Π, S] = 0.The total Hilbert space can be decomposed in mutually orthogonal eigenspaces of Π, i.e., H = ⊕ ns i=1 H i with n s being the number of eigenvalues of Π.In an equivalent way the Hamiltonian of the system can be decomposed As the operator Π is compatible with all the components of the Hamiltonian they share an eigenbasis and each part of the Hamiltonian can be blockdiagonalised.This implies that if |ψ α ∈ H k the following statements hold To deal with density matrices we define the adjoint operator Ṽ ∈ O (H) as Ṽ (̺) = V ̺V † (̺ is the density matrix for the total Hamiltonian).If the eigenvalues of V are ν k = e iΩ k , (k = 1, 2, ..., n s ), the adjoint operator eigenvalues are just the product of these eigenvalues ν k ν l = e i(Ω k −Ω l ) .The adjoint space O (H) can be decomposed as O (H) = ⊕ ns k=1 ⊕ ns l=1 V kl .If a density matrix has the form ̺ kl = |ψ k ψ l |, it means that ̺ kl ∈ V kl .From the spectral decomposition, together with Eq. ( 13), we can derive that if We work in the interaction picture, where the dynamics of the system is given by the interaction Hamiltonian, ̺(t) = −i[H BS (t), ̺(t)], and the time dependence of the interaction Hamiltonian is given by H BS (t) = e i(HS +HP )t H BS e −i(HS +HB )t .
Combining Eqs. ( 14) and (15) it can be proved that It is clear from Eq. ( 16) and the von Neumann equation that if ̺(0) ∈ V kl then ̺(t) ∈ V kl .Consequently, since tracing over the probe degrees of freedom doesn't cause a change in subspaces =⇒ ρ(t) ∈ V kl .This completes our general proof of the invariant subspaces in open quantum systems.
Redfield-Lindblad quantum master equation.The nonequilibrium transport set-up described within this work can be modeled using the total Hamiltonian, wherein H S describes the system and H α describes the left-lead, right-lead and the probe with α = L, R and P respectively.The interaction of the system with the leads and the probe is given via H αS .The leads and the probe Hamiltonian are given by an infinite set of independent harmonic oscillators that read We assume the system-lead/probe interaction to be weak (Born-approximation) and the leads/probe to be memory-less (Markov-approximation).Additionally, only for the leads, we assume that the dissipative effects of the leads should be relevant on a timescale much longer than all finite times of the problem (secularapproximation).Under these assumptions and initiating the total system in a product state of the system, leads, and probes, with no cross-correlations between the leads and the probe, the master equation for the reduced density matrix reads [48,49], Throughout this work we set = 1 for notational simplicity.The secular approximation causes the effect of the leads on the system to be described via the Gorini-Kossakowski-Sudarshan-Lindblad master equation [50,51].The Lindblad operators A αk corresponding to the α-th lead allow local injection (k = 1) or extraction (k = 2) of one exciton to the molecular system H ′ S (see Eq. ( 1)).The rates Γ αi dictate the temperature of each lead, namely, Γ α1 = Γn α and Γ α2 = Γ(n α + 1) with Γ being the dissipation strength of each lead and n α = [exp(β α ω 0 ) − 1] −1 the Bose-Einstein distribution function.The ratio of the rates Γ α1 /Γ α2 = exp(−β α ω 0 ) obeys detailed balance where β α = 1/k B T α represents the inverse temperature of the α-th lead and ω 0 is the characteristic phonon frequency of the lead.
The effect of the probe on the system is described via the Redfield tensor [52] that allows us to study the interplay of various additional timescales (due to the absence of secular approximation), e.g. the system and the probe dynamics timescales.The operator S in the Redfield tensor originates from the system-probe interaction Hamiltonian that we choose to be of a general form H P S = S ⊗ Y .The operator Y is then encapsulated in the correlator C(t) = Tr P Ỹ (t)Y (0) exp (βH P ) where β is the inverse temperature of the probe and Ỹ (t) is the free-evolution of the probe operator Y with respect to the probe Hamiltonian.We choose the probe operator to be the collective position operator Y = − k c k x k with c k being the coupling strength of each harmonic mode of the probe to the system.Throughout this work, the explicit super-script P will be suppressed for convenience.All parameters of the probe are described via a spectral density is chosen to be of the ohmic form with a Lorentz-Drude cut-off frequency ω D and dissipation strength γ.The corresponding correlator C(t) is given by, Full counting statistic for excitonic transport.In order to quantify the excitation flux through the system we include a counting field in the quantum master Eq. ( 19) [7,53,54].We first introduce the reduced density matrix ρ q (t) that is the projection of the density matrix into the subspace of q excitations interchanged between the system and the left-lead at time t.The probability of observing an exciton current J q is thus given by P (q) = Tr [ρ q (t)].The evaluation of this probability is easier after a change of ensemble, made by introducing the Laplace transform ρ λ (t) = q ρ q (t) exp(−λq) with λ being a counting field.The evolution of the density matrix in the Laplace transformed ensemble ρ λ (t) is given by ρλ the above density matrix as a function of the counting field λ the moment generation function can be calculated as Z λ (t) ≡ Tr [ρ λ (t)].Thus, the average exciton current flowing into the system from the left-lead J L can be obtained from the generation function as J L = (1/t) [∂ λ log Z λ (t)] λ=0 .Considering that ρ λ (t) = exp (+tL λ ) ρ(0) we obtain the left-lead expected current as,

4-site model revisited
The results we show in the main text correspond to a system where diagonal and off-diagonal elements of the Hamiltonian H ′ S are chosen to be the same.This requirement can be relaxed.Symmetries could exist even if these parameters are mismatched.In this section we explore a more general case with the help of the 4-site model and show that the signatures of molecular symmetries are robust.Thus, if the system possess a molecular symmetry in any form our probe based approach could help detect these symmetries.We begin by choosing a general 4-site Hamiltonian where the parameters ε i and h ij = h ji are different for each site.In order to maintain the mirror symmetry of the system we choose ε 2 = ε 4 , h 12 = h 14 , and h 23 = h 34 .
Figure 6 shows the behavior of the excitonic currents when the probe is placed either at site 1 (black solid lines) or at site 2 (green dashed-dotted lines) for symmetric (panel a), antisymmetric (panel b) or canonical (panel c) initial conditions.The conclusions remain the same as in the main text, i.e., the signatures of molecular symmetries is observed only in the case of a dark state initial condition.One main difference we observe when all parameters are not the same is the change in the dissipative Liouvillian spectrum (Fig. 6d).Contrary to the fully identical case (see Fig. 2 of main text) the eigenspectrum changes drastically.Although the splitting of the degenerate eigenvalue still persists and the unstable manifold still remains closest to the zero eigenvalue thus influencing the long relaxation times of the system.Effect of probe position in para-Benzene Similar to the 4-site model shown in the main text the effect of a local probe on benzene is shown in Fig. 7.The qualitative features remain the same as in the 4-site system with the anti-symmetric state showing the clear signatures of multiple steady states.The Liouvillian spectrum also shows exotic modes like oscillating coherences [1] that correspond to eigenvalues of the Liouvillian with zero real part but finite imaginary part.These could be the modes responsible for the observation of circular currents observed in benzene [2].

FIG. 1 .
FIG. 1.Effect of probe position on symmetry detection.Panel a depicts the modulation of the steady-state excitonic currents in the 4-site model as a function of the mixing angle, Eq. (5), that determines the initial condition.Panels b, c, and d show the time evolution of the excitonic current with symmetric (panel b), antisymmetric (panel c), and canonical (panel d) initial conditions.The probe is positioned at sites 1 (black solid line) and 2 (green dashed-dotted line) for all panels.The black crosses and open green circles correspond to the nonequilibrium steady state values of the excitonic currents.The steady state values for panel b and d are small but non-zero.Note the different time unit in panel (c) indicating a very slow relaxation from the antisymmetric initial condition.The closed circles in panel (a) and (b) represent results for a weakly disordered system with disorder strength δ = 0.13meV.The currents are averaged over 1000 samples of uniformly distributed random disorder [standard error of mean (SEM) is contained within each closed circle].The black closed circles correspond to the probe being placed at site 1, whereas the green closed circle correspond to the probe at site 2. The system parameters are: ε = −142.2meVand h = −9.35meV.The lead and probe parameters are chosen as: TL = 330K, TR = 270K, T = 300K, Γ = 196GHz, γ = 19.6GHz,ω0 = 78.55THz,and ωD = 1.96THz.
due to the existence of a unitary operator Π = exp (|e 2 e 4 | + |e 4 e 2 |).The operator Π, also known as the symmetry operator, obeys [Π, H S ] = [Π, A αk ] = 0 (α = L, R; k = 1, 2) due to which the time evolution of a general density matrix ρ can be decomposed into invariant subspaces (see methods section for a general proof).It is worth noting here that the Lindblad description provides a minimalist setup and does not affect the overall results of this work.

FIG. 2 .
FIG. 2. Dissipative Liouvillian spectrum for the 4-site model.(a) Eigenspectrum of the dissipative Liouvillian for the 4-site model.The asterisk correspond to the eigenvalues Λ when the probe is placed at site 1, whereas the open green triangles correspond to probe positioned at site 2. (b) shows a magnification around the zero eigenvalue.There are 2 degenerated eigenvalues at zero when the probe is positioned at site 1.These eigenvalues split when the probe is at site 2 and the resultant gap dictates the relaxation time τR.The dashed black line marks the Im[Λ] = 0 axis.All other parameters are the same as in Fig. 1.

FIG. 3 .
FIG. 3. Probe temperature T dependence on the dynamics of excitonic currents.(a) shows the current as a function of time for different temperatures [T = 200K (blue curve), 300K (black curve), and 400K (red curve)] when the cutoff frequency of the probe is very large ωD = 58.80PHz, i.e. τP ≪ τS, for any temperature the behavior is the same and all the lines collapse, whereas (b) depicts the small cutoff frequency case ωD = 1.96THz, i.e., τP ≫ τS.Inset in panel (a) depicts the normalized relaxation time τnorm = τ R|T /τ R|T 0 where τ R|T and τ R|T 0 are the relaxation times evaluated at temperature T and T0 = 100K.The solid curve in the inset is for the case when τP ≫ τS, whereas the dotted line is for τP ≪ τS.In both panels the system is initiated in an antisymmetric state.All other parameters are the same as in Fig.1.

( 3 )
and (4).The symmetry operator for the para-benzene ring Π = exp[(|2 6| + |6 2|) ⊗ (|3 5| + |5 3|)] due to which the system has 3 multiple steady states.Only one of these multiple steady states contains nonequilibrium information (e.g.temperature T L , T R or coupling strength Γ) about the leads.The other two are pure states are decoupled from the leads and depend only on the symmetries.

FIG. 4 .
FIG. 4. Eigenvalues of the dissipative Liouvillian for the para-benzene model.(a) Eigenvalues closest to zero when the para-benzene symmetry is broken via a local probe acting on site 2 (as depicted in the illustration above).(b) the probe is acting non-locally on sites 2 and 3 (illustrated above panel b) to break all the symmetries of the system.The black plus corresponds to the eigenvalues of the perturbation matrix δL [Eq.(10)], whereas the green open triangles correspond to the smallest (in magnitude) eigenvalues of the full dissipative Liouvillian.The dashed black line marks the Im[Λ] = 0 axis.The system parameters are: ε = −11.2eVand h = −0.7eV.The lead and probe parameters are chosen as: TL = 330K, TR = 270K, TP = 300K, Γ = 151.9THz,γ = 15.19THz,ωD = 151.92THzand ω0 = 78.55THz.

FIG. 5 .
FIG. 5. Bi-exponential decay of currents as a signature of 3 steady states.Long-time dynamics of the para-benzene model exhibiting a bi-exponential decay, due to which a peak is formed, when the system is initiated in an antisymmetric state w.r.t.sites 2 and 5.The probe acts non-locally on sites 2 and 3 as illustrated above panel (b) in Fig. 4. The open green circle corresponds to the steady state value for the current.The closed green circles are for a weakly disordered system with disorder strength δ = 1meV.The currents are averaged over 1000 samples of uniformly distributed random disorder [standard error of mean (SEM) is contained within each closed circle].All other parameters are same as Fig. 4.

1 FIG. 6 . 4 FIG. 7 .
FIG. 6.Effect of probe position on symmetry detection.Time evolution of the excitonic current in the 4-site model with symmetric (panel a), antisymmetric (panel b), and canonical (panel c) initial conditions.(d) shows the eigenspectrum of the dissipative Liouvillian and (e) depicts the magnification around the zero eigenvalue.The dashed black line in panels (d) and (e) marks the Im[Λ] = 0 axis.The probe is positioned at sites 1 (black solid line in panels a, b and c; black asterisk in panels d and e) and 2 (green dashed-dotted line in panels a, b, and c; green triangles in panels d and e) for all panels.The crosses and circles correspond to the nonequilibrium steady-state values of the excitonic currents.The system parameters are: ε1 = −64.6meV,ε3 = −193.9meV, ε2 = ε4 = −129.3meV, h12 = h14 = −14.52meV, and h23 = h34 = −1.59meV.The lead and probe parameters are chosen as: TL = 330K, TR = 270K, T = 300K, Γ = 196GHz, γ = 19.6GHz,ω0 = 78.55THz,and ωD = 1.96THz.