Long-lived circulating currents in strongly correlated nanorings

TL;DR

This paper investigates the dynamics of currents in strongly correlated nanorings, revealing a regime where circulating currents decay slowly and persist far from equilibrium despite stationary transmitted current.

Contribution

It uncovers the existence of long-lived circulating currents in nanorings and analyzes their decay rates as a function of interaction strength, highlighting a novel non-equilibrium regime.

Findings

Circulating currents decay at a rate much slower than lead-induced broadening.

Long-lived ring transients can exist even when transmitted current is stationary.

Decay rate of circulating currents depends strongly on interaction strength.

Abstract

We study the time evolving currents flowing in an interacting, ring-shaped nanostructure after a bias voltage has been switched on. The source-to-drain current exhibits the expected relaxation towards its quasi-static equilibrium value at a rate $\Gamma_0$ reflecting the lead-induced broadening of the ring states. In contrast, the current circulating within the ring decays with a different rate $\Gamma$, which is a rapidly decaying function of the interaction strength and thus can take values orders of magnitude below $\Gamma_0$. This implies the existence of a regime in which the nanostructure is far from equilibrium even though the transmitted current is already stationary. We discuss experimental setups to observe the long-lived ring transients.

Equations

$$H_{\text{r}}=$$

$$+U\sum _{⟨i,j⟩}\left(n_i{n}_j-\frac{n_i+n_j}{2}\right)+{\epsilon }_T{n}_2,$$

$$H_{\text{l}}=-J\sum_{\alpha=\text{L,R}}\sum_{n\geq 0}\bigl{(}c^{\dagger}_{\alpha,n+1}c_{\alpha,n}+c^{\dagger}_{\alpha,n}c_{\alpha,n+1}\bigr{)},$$

$$H_{\text{c}}=-J_{\text{c}}\left(d_1^{†}{c}_{\text{L},0}+c_{\text{L},0}^{†}{d}_1+d_4^{†}{c}_{\text{R},0}+c_{\text{R},0}^{†}{d}_4\right),$$

$$\begin{split}H_{\text{SW}}=&\sum_{k,\sigma}\epsilon_{k}c_{\text{res},k\sigma}^{\dagger}c_{\text{res},k\sigma}+hS^{z}+\tilde{h}S_{\text{res}}^{z}\\ &+J_{\perp}\bigl{(}S^{x}S_{\text{res}}^{x}+S^{y}S_{\text{res}}^{y}\bigr{)}+J_{z}S^{z}S_{\text{res}}^{z},\end{split}$$

$$\Gamma=\frac{\pi\rho_{0}^{2}J_{\perp}^{2}}{16}\bigl{(}|\varepsilon_{T}+V|+|\varepsilon_{T}-V|+2|\varepsilon_{T}|\bigr{)}+\frac{\pi\rho_{0}^{2}J_{z}^{2}}{8}V.$$

$$|{\psi }_0⟩=\sqrt{N}\left(\sqrt{{\rho }_{gg}}|2,g⟩\pm \sqrt{{\rho }_{ee}}|2,e⟩\right),$$

$$|\psi (t)⟩=\exp (-i{H}_{\text{r}}t)|{\psi }_0⟩,$$

$$H(t=0)=H+\frac{V}{2}\left(\sum _i{n}_{\text{L},i}-\sum _i{n}_{\text{R},i}\right),$$

$$I_{\text{t}}$$

$$I_{\text{u}}$$

$$I_{\text{l}}$$

$$f(t)=a\exp (-\Gamma t)\cos (\epsilon t+b)+c,$$

$$H_{\text{r},8}=-J\sum _{⟨i,j⟩}\left(d_i^{†}{d}_j+\text{h.c.}\right)+U\sum _{⟨i,j⟩}\left(n_i{n}_j-\frac{n_i+n_j}{2}\right)+{\epsilon }_T{n}_4,$$

$$H_{\text{c}}=J_{\text{c}}\left(d_1^{†}{c}_{\text{L},0}+d_8^{†}{c}_{\text{R},0}+\text{h.c.}\right),$$

$$LA\equiv {\left[H,A\right]}_{-}.$$

$$\dot{\rho }(t)=-i{\left[H,\rho (t)\right]}_{-}=-iL\rho (t),$$

$$\rho (t)=\exp \left[-iL(t-t_0)\right]\rho (t_0).$$

$${\tilde{\rho }}_{\text{ns}}(E)=$$

$$=$$

$${\tilde{\rho }}_{\text{ns}}(E)=i\sum {\text{tr}}_{\text{l}}\frac{1}{E-L_{\text{l}}-L_{\text{ns}}}{L}_V\dots L_V\frac{1}{E-L_{\text{l}}-L_{\text{ns}}}{\rho }_{\text{l}}{\rho }_{\text{ns}}(t_0).$$

$$L_{\text{eff}}(E)=L_{\text{ns}}+\Sigma (E).$$

$$\frac{i}{E-L_{\text{eff}}(E)}.$$

$$z-L_{\text{eff}}(z).$$

$$L_V=G_1^{p_1}:J_1^{p_1}:,$$

$$\displaystyle J^{p}_{1}A=\left\{\begin{array}[]{ll}c_{1}A&p=+\\ Ac_{1}&p=-\end{array}\right.\,,$$

$$\displaystyle c_{1}=\left\{\begin{array}[]{ll}c^{\dagger}_{\nu,\omega}&\eta=+\\ c_{\nu,\omega}&\eta=-\end{array}\right.\,.$$

$$\displaystyle G^{p}_{1}A=\left\{\begin{array}[]{ll}d_{1}A&p=+\\ -\sigma^{p}Ad_{1}&p=-\end{array}\right.\,.$$

$$|l)=G_1^{-}|m)=|i⟩⟨j|,$$

$$\text{mod}\left[\left(\sum _i{d}_i^{†}{d}_i|i⟩-\sum _i{d}_i^{†}{d}_i|j⟩\right),2\right]=1.$$

$${\gamma }_{11^{\prime }}^{p{p}^{\prime }}=⟨J_1^p{J}_{1^{\prime }}^{p^{\prime }}{⟩}_{\text{eq.}},$$

$${\gamma }_{11^{\prime }}^{p{p}^{\prime }}={\delta }_{1{\bar{1}}^{\prime }}{p}^{\prime }{\gamma }_1^s+{\delta }_{1{\bar{1}}^{\prime }}{\gamma }_1^a,$$

Full text

Introduction.—Isolated quantum systems, such as small molecules, feature a discrete set of energy levels. When brought to contact with two electrodes, a nano-junction can form and a current begins to flow. At weak coupling, the associated level broadening, $\Gamma_{0}$, is still small as compared to the typical energy spacing, $\Delta E$, of the isolated system. One might perhaps suspect that these energies by themselves set the only relevant time scales. But in fact a prominent exception is known, the Kondo phenomenon Hewson93 , which occurs in a situation where $\Delta E$ is dominated by strong on-site repulsion between the charge carriers. This suppresses charge fluctuations but allows for spin fluctuations, leading to an emergent energy scale, the Kondo temperature $T_{\text{K}}$, which is parametrically small compared to the native scales $\Gamma_{0}$ and $\Delta E$.

In this work, we report another example of an emergent energy scale, $\Gamma$; it manifests in the relaxation of circulating currents in mesoscopic nanostructures. Like the Kondo temperature, the new scale is a many-body phenomenon, originating from interactions between particles on the nanostructure. However, the manifestation of the new relaxation rate $\Gamma$ requires the nanostructure to be brought out of equilibrium.

A sketch of a minimal model system that exhibits the novel scale $\Gamma$ is displayed in Fig. 1. Originally, similar ring-shaped devices served as a toy-models to study the interplay of interaction and interference BohrSchmitteckert12 ; Schmitteckert13 and to explain quantum-interference effects in transport through functionalised graphene ribbons Walz-14 . The ring geometry supports stationary circulating (“orbital”) currents that can exceed the source-drain (“transport”) current by orders of magnitude at Fermi-energies situated close to a Fano-resonance.

Strong circulating currents in ring-shaped devices, Fig. 1, generically arise as transients after a voltage quench. They then carry an oscillating amplitude with a frequency resembling the lowest lying excitation gap of the nanostructure. We here report results from time-dependent density matrix renomalisation group (tdDMRG) Vidal04 ; WhiteFeiguin04 ; Daley-04 ; Schmitteckert04 simulations showing that in situations where the interaction $U$ is the dominating native scale of the nanoring, these oscillations can be very pronounced and very long-lived. They exhibit a lifetime $\Gamma^{-1}$ that exceeds the transients in transport currents, $\Gamma_{0}^{-1}$, by orders of magnitude if the repulsive interaction $U$ becomes strong. The computational finding is complemented with perturbative arguments that explain this effect and clarify the relevant physical processes. In essence, the strong suppression of $\Gamma$ originates in a large energy gap between the two low-lying states and the rest of the spectrum on the nanoring (see Fig. 2), with ring current connecting the low-lying states. Possible experimental signatures of the effect proposed here are discussed. We note that in contrast to previously discussed MedenSchollwoeck03 ; Molina-03 ; RejekRamsak03 persistent ring currents driven by magnetic fields, the oscillating ring current we observe is a pure non-equilibrium effect.

Nanostructure.—The model associated with Fig. 1 is represented by the Hamiltonian $H=H_{\text{r}}+H_{\text{l}}+H_{\text{c}}$ describing the ring, the leads and their mutual coupling, respectively. The ring Hamiltonian is given by

[TABLE]

with operators $d^{\dagger}_{j}$ and $d_{j}$ creating/annihilating spinless fermions at site $j$ and $n_{j}=d^{\dagger}_{j}d_{j}$ denoting the corresponding density. The first term describes hopping of the fermions between nearest neighbours, while the second represents the repulsive nearest-neighbour interaction. The last term is an external potential at the top site which breaks the symmetry between the upper and lower path through the ring. The lead Hamiltonian reads

[TABLE]

where $c^{\dagger}_{n,\alpha}$ and $c_{n,\alpha}$ create and annihilate a spinless fermion at site $n$ in the lead $\alpha{=}\text{L,R}$. For simplicity we assume the hopping parameter $J$ in the ring and lead to be equal. Finally, the coupling between both subsystems is facilitated by

[TABLE]

coupling the outer sites on the ring to the leads.

In the following we analyse the non-equilibrium currents in the nanostructure by three different methods: (i) tdDMRG simulations, (ii) a reduced density-operator transport theory (RDTT) Schoeller09 ; SaptsovWegewijs12 , and (iii) mapping to an effective two-state nanostructure Bravyi-11 .

tdDMRG simulations.—First, we study the time evolution after a voltage quench using the tdDMRG algorithm BohrSchmitteckert07 ; Boulat-08 ; Kirino-08 ; DiasdaSilva-08 ; Heidrich-Meisner-09 ; Branschadel-10 ; Schwarz-18 . Specically we use the time evolution scheme outlined in Refs. Schmitteckert04 ; Branschadel-10 ; supplement performing the evaluation of the time evolution via matrix exponentials within the framework of Krylov spaces. At times $t<0$ the system is prepared in the ground state of the model with an additional charge excess induced by a stationary gating with $V/2\left(\sum_{i}n_{\text{L},i}-\sum_{i}n_{\text{R},i}\right)$. At $t=0$ the gate is switched off, so the electrodes begin to discharge and currents start to flow through the system. We simulate the time evolution with finite leads which are long enough to be able to study the transient regime all the way into the quasi-stationary, non-equilibrium limit. Finite-size effects will interfere only at times exceeding the recurrence time $\tau_{\text{rec}}{=}L/(2v_{\text{F}})$, at which the electrons reach the boundary of the leads. (For details of the quenching protocol see Ref. Branschadel-10 .) Here $L$ denotes the total number of sites, ie, the length of the leads is given by $(L-4)/2\approx L/2$, and $v_{\text{F}}=2J$ is the Fermi velocity of the lead electrons.

During the time evolution we determine the expectation values of the local currents $I_{\text{t}}\propto\text{Im}(c^{\dagger}_{i}c_{i-1})$ and $I_{\text{r}}\propto\text{Im}(d^{\dagger}_{k}d_{l})$ flowing in the leads and the impurity, respectively, where $l$ and $k$ are neighbouring sites. The local current densities after quenching are displayed in Fig. 2. The transport (“transmitted”) current $I_{\text{t}}$ initially fluctuates in response to the quench for times $\Gamma_{0}t\leq 3$, where we use $\Gamma_{0}=2\pi\rho_{0}J_{\text{c}}^{2}$ with the density of states in the leads $\rho_{0}=1/(2\pi J)$ as our time unit. After this transient the transmitted current appears to have reached a largely time-independent steady state in line with predictions from non-equilibrium Green function formalism Jauho-94 ; Tuovinen-13 .

In contrast, for the local currents in the ring $I_{\text{r}}$ we observe a drastically different behaviour. Although some transient features decay quickly, the ring currents oscillate with a distinct frequency $\omega$ for long times. In fact, for sufficiently strong Coulomb repulsions $U$ we do not observe a significant reduction of the oscillation amplitude within the observation times accessible to our simulations. Qualitatively similar results were obtained for a ring structure with eight sites supplement .

The frequency of the oscillations can be understood based on the spectrum supplement of the uncoupled ring $H_{\text{r}}$ shown in Fig. 2. We find that the frequency $\omega$ extracted from the tdDMRG simulations matches the energy gap between the two lowest-lying states on the ring. These two states can be identified as charge-density wave (CDW) states with $N=2$ particles on the ring, one being the ground state $|2,g\rangle$ and the other the first excited state $|2,e\rangle$. Thus we confirm that the ring current originates from the mixing of these two states by the time evolution, which is driven by the coupling of the ring to the leads as exemplified by the proportionality of the ring current to the coupling $\Gamma_{0}$ shown in the inset of Fig. 5.

The decay rate $\Gamma$ of the ring currents is rapidly decreasing with the interaction strength $U$, see Fig. 5, exhibiting a wide regime with $\Gamma\ll\Gamma_{0}$. To understand the origin of this regime, we proceed with the RDTT analysis.

RDTT analysis.—The RDTT Schoeller09 ; SaptsovWegewijs12 method aims at determining the time evolution of the reduced density matrix of the nanostructure, $\rho_{\text{ns}}(t)=\text{tr}_{\text{l}}\rho(t)$, where the trace is taken over the lead degrees of freedom in the density matrix $\rho(t)$ of the full system. The time evolution of $\rho_{\text{ns}}(t)$ can be cast in the form $\dot{\rho}_{\text{ns}}(t)=-\text{i}L_{\text{ns}}\rho_{\text{ns}}(t)$, with the effective Liouvillian $L_{\text{ns}}$ governing the relaxation of the nanostructure. Since the ring current originates from the mixing of the two CDW states $|2,g\rangle$ and $|2,e\rangle$, its decay is related to the decay of the off-diagonal elements $\rho_{ge}$ and $\rho_{eg}$ of $\rho_{\text{ns}}$. We have determined the corresponding decay rate from the Liouvillian $L_{\text{ns}}$ calculated supplement to first order in the bare coupling rate $\Gamma_{0}$, with the perturbative regime set by $\Gamma_{0}\ll T$ with the temperature $T$.

The obtained results for the decay rate $\Gamma$ of the ring current are shown in Fig. 9(a). The results are qualitatively similar to the ones obtained via tdDMRG shown in Fig. 5 in the sense that the rate is strongly suppressed at large $U$. The quantitative differences between the RDTT and tdDMRG results reflect the fact that both methods operate in different parameter regimes.

Furthermore, the RDTT allows us to identify supplement the relaxation processes contributing to the decay rate, which are visualised in Fig. 9(b). The dominant processes are shown in sketches (1) and (2), which involve the tunneling of a particle off or onto the ring, while the sub-leading processes are shown in sketches (3) and (4). All processes are constraint by energetics: (1) and (2) only contribute in the regions (i) and (ii) in Fig. 9(a), (3) only in regions (i) and (iii), and (4) is relevant in the regions (i)–(iv). We stress that in region (v) no relaxation processes in order $\Gamma_{0}$ exist. Thus at sufficiently large interaction strengths $U$ the rate $\Gamma$ essentially drops to zero (to order $\Gamma_{0}^{2}$), explaining the very slow decay of the ring current.

Schrieffer–Wolff transformation.—Finally we focus on the regime of strong interactions, $U/\text{max}(\varepsilon_{T},J)\rightarrow\infty$, where we can derive the analytic dependence $\Gamma\sim U^{-6}$ consistent with our computational results, Fig. 5. As can be seen from the spectrum of the bare ring (Fig. 2), in this limit the two CDW states $|2,g\rangle$ and $|2,e\rangle$ will be well separated by an energy splitting $\Delta(U)\sim U$ from the higher excited states. It is thus instructive to construct an effective two-level system containing only these states, where the couplings to the higher excited states are treated using a Schrieffer–Wolff transformation Bravyi-11 in fourth order in the couplings $J$ and $J_{\text{c}}$. Going to this order in the expansion is necessary since all off-diagonal matrix elements exactly cancel in second order due to the mirror symmetry of the isolated ring structure supplement .

The resulting two-level system can be written in the form of an electronic Kondo model, with the localised spin identified with the CDW states as $|\!\downarrow\rangle=|2,g\rangle$ and ${|\!\uparrow\rangle}=|2,e\rangle$ and the corresponding spin operator denoted by $\vec{S}$. An effective reservoir electronic degree of freedom can be formed via $c_{\text{res},\uparrow\downarrow}=(c_{\text{L}}\pm c_{\text{R}})/\sqrt{2}$ from the leads (2) of the original model; the effective spin operator formed from the first sites ($n=0$) is denoted by $\vec{S}_{\text{res}}$. With this notation the effective model reads supplement

[TABLE]

where the first term is the energy of the electronic reservoir, the second and third are effective magnetic fields $h\approx\varepsilon_{T}$ and $\tilde{h}=\mathcal{O}(U^{-4})\ll h$ acting on the two-level system and spin of the electron reservoir, and the fourth and fifth term represent a Kondo coupling between the two, with the coupling being strongly anisotropic with $J_{z}\simeq 10J^{2}J_{\text{c}}^{2}/U^{3}$ and $J_{\perp}=\mathcal{O}(U^{-5})$, and thus $|J_{\perp}|\ll|J_{z}|\ll J_{\text{c}},J$.

Due to the formation of the effective reservoir electron spin from the leads L,R the bias voltage $V$ enters the effective Kondo model in the form of a transverse field in the reservoir, ie, as $V/2\sum_{k\sigma\sigma^{\prime}}c_{\text{res},k\sigma}\tau^{x}_{\sigma\sigma^{\prime}}c_{\text{res},k\sigma^{\prime}}$ with $\tau^{x}$ being the x-component of the Pauli matrices. Finally, the ring current corresponds to oscillations between the two CDW states and thus is related to the localised spin via $I_{\text{r}}\sim S^{y}$. Performing a suitable spin rotation in the electronic reservoir we calculated supplement the corresponding relaxation rate using standard perturbation theory in the Kondo system Schoeller09 ; Rosch-03prl ; SchoellerReininghaus09 with the result

[TABLE]

We stress that in the considered regime of strong interactions this rate is vanishingly small, $\Gamma\sim\rho_{0}^{2}J_{z}^{2}V\sim\rho_{0}^{2}J^{4}J_{\text{c}}^{4}V/U^{6}$, in accordance with our finding of long-lived oscillations in the ring current. In particular, the predicted behaviour $\Gamma\sim U^{-6}$ is consistent with our tdDMRG simulations shown in Fig. 5. We note that the result (5) is applicable deep in region (v) of Fig. 9(a), where we found that processes of order $\Gamma_{0}$ vanish. Furthermore, the effective model (4) will show the Kondo effect, however, the relevant energy scale $T_{\text{K}}$ will be much smaller than the energy scales we consider here, in particular $T_{\text{K}}\ll\varepsilon_{T}$. Thus the equilibrium Kondo effect is not observable in our setup.

Finally we note that nanostructures with two energetically well separated low-lying states can generically be approximated by an effective Kondo model using a Schrieffer–Wolff transformation. In the absence of the above mentioned mirror symmetry the exchange couplings will be of the order $J_{z},J_{\perp}\sim J_{\text{c}}^{2}/U\ll J_{\text{c}}$, resulting in a relaxation rate $\Gamma\sim\rho_{0}^{2}J_{z}^{2}V\sim U^{-2}\ll\Gamma_{0}$. Thus ring currents that couple to these low-lying states are still expected to decay very slowly.

Experimental verification.—We see a possible experimental realisation of the ring-shaped model system, Fig. 1, in molecules such as porphyrines or phthalocyanines. Single molecule conductance measurements have indeed been performed at these systems Sedghi-11 ; Schmaus-11 ; Bagrets-12 so the possibility for bias-ramping has also been demonstrated already. As an observable indicating the slow decay of the ring currents we propose to measure the photons that are emitted when these currents decay via coupling to the radiation field. In this context we note that single-molecule electroluminscence measurements have been performed Marquardt-10 ; Reecht-14 already and thus are indeed experimentally feasible. An alternative realisation of our ring-shaped model may be provided by quantum dot arrays Mukhopadhyay-18 , which in particular offer a high level of control of the couplings and allow to enter the regime of strong interactions essential for the long-lived ring currents.

Conclusion.—We have studied the relaxation of transport processes in an interacting ring-shaped nanostructure. Owing to a mirror symmetry of the Hamiltonian, the system supports oscillating ring currents long after the transmitted current has died out, with the ratio $\Gamma/\Gamma_{0}$ of the respective relaxation rates being strongly suppressed by the interactions. Our work provides a striking example for an untypical situation in thermodynamic relaxation processes: Two observable currents approach their equilibrium values on timescales that are parametrically separated with rates differing by orders of magnitude. In addition, our system provides new insight into the field of quantum devices as we show that internal oscillations can be longer-lived than observed in currents through the system. While we have focused on a ring-shaped nanostructure, the appearance of the suppressed relaxation rate $\Gamma$ is generally expected in systems that can be effectively described by a two-level model with the ring current connecting the low-lying states.

We thank Theo Costi, Mikhail Pletyukhov and Peter Wölfle for useful discussions. This work is part of the D-ITP consortium, a program of the Netherlands Organisation for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW). BMS and PS thank the HPC project QWHISTLE at the Steinbuch Centre of Computing at Karlsruhe Institute of Technology (KIT). PS was supported by ERC-StG-Thomale-TOPOLECTRICS-336012. FE thanks the DFG for support under grant EV30/08-1 and SFB 1277 project A03. BMS and DS were supported by the Netherlands Organisation for Scientific Research (NWO) under FOM 14PR3168.

I Exact diagonalization of the decoupled ring impurity

Spectrum and particle densities

We have performed an exact numerical diagonalization of the Hamiltonian matrix $H_{\text{r}}(U,\varepsilon_{T},J)$ of the ring impurity in the absence of the leads. In figure 2 we plot the relative spectrum $(E-E_{0})$ for the for eigenstates with the lowest energy. The energy of these states is shown as a function of the interaction strength $U$ and a gate potential $\varepsilon_{T}=J/2$. The ground state features half-filling of the ring ($n=2$) and is indicated by the red line. The other eigenstate in the spectrum with half-filling is shown as the blue line. The state marked by the orange line features only a single electron in the ring while the state indicated by the magenta line has three electrons in the ring. For interaction $U/\varepsilon_{T}\geq 1$ we observe an increasing energy separation between the two eigenstates at half-filling and the rest of the spectrum. When comparing the frequency of the observed oscillations of the local currents in the ring with the relative spectrum of the ring, we find an excellent agreement of the frequencies with the energy gap between the ground state $|2,g\rangle$ and the second eigenstate at half-filling $|2,e\rangle$. The frequencies that we have obtained from the fit of a cosine function to the data of the ring current are displayed as black dots in figure 2. We show the local electron density on the ring sites for the four low energy eigenstates in figure 2. We find that the two eigenstates at half-filling exhibit characteristics of charge density waves. The ground state has a significantly increased electron density on site $1$ and $4$ of the ring, while the excited state features an increased density on sites $2$ and $3$. The other two states have a more evenly distributed electron density. We will therefore refor to the states $|2,g\rangle$ and $|2,e\rangle$ as charge density wave (CDW) states from now on.

Time evolution of an initial superposition

We have performed DMRG calculations of the time-dependent reduced density matrix of the ring impurity. We find finite occupation probabilities for both CDW states at time $t=0$. With increasing bias voltage, the occupation probability of the excited CDW state tends to grow as well. We have used these occupation probabilities from the DMRG to construct an initial pure state

[TABLE]

where $\sqrt{N}$ is a normalization factor, $\rho_{gg}$ refers to the ground state occupation probability and $\rho_{ee}$ to the occupation probability of the excited CDW state. Using exact diagonalization we then perform the time evolution of this initial state in the bare ring impurity as

[TABLE]

and calculate the expectation values $\langle I_{\text{u}}\rangle$ and $\langle I_{\text{l}}\rangle$ of the local currents in the ring. The results of this calculation are in good agreement with our DMRG results in both amplitude and frequency.

II DMRG

DMRG implementation

For our numerical calculation of the time evolution of the complete system including ring impurity and leads we have employed a typical finite lattice Density Matrix Renormalization Group (DMRG) algorithm. We keep a maximum of $N_{\text{cut}}=2800$ states per block and set the maximum amount of discarded entanglement entropy to $\delta S_{\text{max}}=10^{-7}$ in each DMRG step. We use a Krylov subspace method to calculate the matrix exponential, allowing us to chose larger time steps up to $\Delta t$ of order one. Each state that is reached through application of the matrix exponential onto the initial state $|\psi_{0}\rangle$ is included into the density matrix from which we determine the subspace of the Hilbert that we project onto in each DMRG step. At each time step we measure the observables of interest as $\langle\psi(t)|\mathcal{O}|\psi(t)\rangle$ where the operator $\mathcal{O}$ has also been projected onto the retained subspace of the Hilbert space.

Quench protocol

At time $t=0$ we prepare the system in the ground state of

[TABLE]

and perform the time evolution using $H(t>0)=H$. We simulate time evolution up $t\leq L/2v_{F}$, where $L$ is the length of the chain (usually $L\geq 72$) and $v_{F}=2J$ is the Fermi velocity of the fermions in the leads. During the time evolution we measure the expectation value of the local currents in the leads as

[TABLE]

and on specific bonds $1\to 2$ and $1\to 3$ in the ring (see Fig. 1 in the main paper) as

[TABLE]

For the majority of our calculations we have used a set of default parameters, namely $L=72$, $\varepsilon_{T}=J/2$, $J_{\text{c}}=J/2$ and $V=0.4\,J/e$.

Detailed discussion of the DMRG time evolution results

In figure 4 we plot the time-dependent expectation values of the operators $I_{\text{t}}$, $I_{\text{u}}$ and $I_{\text{l}}$ using our default parameters and interaction strengths $U/J\in\left\{0.1,0.5,1.0,2.0\right\}$. We begin by discussing the results for weak interaction $U/J=0.1$ shown in figure 4 (a). For the transmitted current $\langle I_{\text{t}}\rangle(t)$ we observe significant initial oscillations inside the typical transient regime $\Gamma_{0}t\leq 1$ that appear to have decayed for $\Gamma_{0}t>1$ while a weak periodic oscillation remains even for large times. This periodic oscillation is not physical but a known finite size effect with a frequency $\omega\equiv V$. For the local currents in the ring we first verify that $I_{\text{u}}+I_{\text{l}}=I_{\text{t}}$ as a consistency check of our results. For times $\Gamma_{0}t\leq 1$ we find the oscillations of $\langle I_{\text{u}}\rangle(t)$ and $\langle I_{\text{l}}\rangle(t)$ small when compared to the oscillations of $\langle I_{\text{t}}\rangle(t)$. The finite size effect with $\omega=V$ for the the local currents in the ring on the other hand is large when compared to the transmitted current. We also indicate $(\langle I_{\text{u}}\rangle-\langle I_{\text{l}}\rangle)/2$ as a dashed black line in fig. 4. This observable corresponds to a ring current in clockwise direction. For interaction strength $U/J=0.5$, shown in fig. 4 (b), we solely observe quantitative differences for $\langle I_{\text{t}}\rangle(t)$. While the initial transient features remain largely unchanged, the steady state current for $\Gamma_{0}t>1$ is reduced. For $(\langle I_{\text{u}}\rangle-\langle I_{\text{l}}\rangle)/2$ we observe what seems to be an initial oscillatory feature that is not due to finite size effect for $\Gamma_{0}t<4$. Due to the small window $1\leq\Gamma_{0}t\leq 4$ a fit does not yield reliable results for frequency and decay rate. For $U/J=1$ the steady state value of the transmitted current experiences yet another significant reduction, whereas the transient features remain of similar size as for $U/J=0.1$. We stil observe that the transient features of the transmitted current have largely decayed by $\Gamma_{0}t=1$. For the ring currents we find a qualitatively different behavior. The ring current exhibits periodic oscillations with a distinct frequency and a visible decay rate $\Gamma$ which is an order of magnitude smaller than $\Gamma_{0}$. For the directional ring current $(\langle I_{\text{u}}\rangle-\langle I_{\text{l}}\rangle)/2$ there is even a window in which the direction of the current has changed. By increasing the interaction strength to $U/J=2$ we find yet another decrease of the steady state trnamsitted current. In the transient regime $\Gamma_{0}t\leq 1$ we now also observe an additional sign change of the transmitted current. We also no longer see the oscillations due to the finite system size. The oscillations of the local currents in the ring $I_{\text{u}}$ and $I_{\text{l}}$ become even more pronounced and feature a periodic change of direction. Through a fit we find that the decay rate of these oscillations is an order of magnitude smaller than in the case $U/J=1$ and now amounts to $\Gamma/\Gamma_{0}\approx 1/100$. There is a clear separation of scales between the typical decay rate $\Gamma_{0}$ which holds for the transmitted current and the decay rate $\Gamma$ of the local currents in the ring impurity. Calculations for stronger interaction $U/J>2$ show a continuation of this trend.

Fitting procedure for the computation of

$\Gamma$

To determine the oscillation frequency $\varepsilon$ and decay rate $\Gamma$ we fit a function

[TABLE]

to our DMRG data for the local currents where $\Gamma$, $\varepsilon$, $a$, $b$ and $c$ are fitting parameters. The fit is performed for $\Gamma_{0}<t<L/2v_{F}$. This fitting procedure only yields reliable results for $0.5<U/J\leq 5$. For weak interaction $U/J\leq 0.5$ the decay time is too short to observe the amount of sine waves necessary to reliably determine the decay rate. For very strong interaction the decay rate becomes so small that it does not lead to a visible reduction in oscillation amplitude for $t<L/2v_{F}$. As a result, the fitting error associated with decay rate becomes larger than the decay rate itself. These limitation of the fitting procedure limit the application of our DMRG method as a tool to determine the decay rate $\Gamma$ to a parameter range $0.5<U/J\leq 5$ as indicated in figure 5.

II.1 DMRG calculations for the decay rate $\Gamma$

We have performed a set of DMRG calculations to study the behavior of the decay rate $\Gamma$ as a function of $U/J$ for a range of specifically chosen parameters $U$, $\varepsilon_{T}$ and $V$. The results of these calculations are shown in figures 5 (a) and (b). Due to the aforementioned limitations of our fitting procedure it is not possible to quantify $\Gamma$ for $0.5<U/J\leq 5$. In the vicinity of $U/\varepsilon_{T}\simeq 1$ a comparison of the log-linear and log-log plots indicates a small region of exponential suppression. For stronger interactions $U/\varepsilon_{T}>1$ we observe a power law behavior of the decay rates as a function of $U/J$. The fit of a power law to the data indicates a smaller exponent for smaller values of $\varepsilon_{T}$. In the case of $\epsilon_{T}=0.5$ we are safely in the regime $U/\epsilon_{T}\gg 1$ for $U/J\geq 4$. In this regime one could consider the data comparable to results obtained in the limit $U/\varepsilon_{T}\rightarrow\infty$. The fit of a power law finds an exponent $\alpha=6.0\pm 0.4$ in this case.

II.2 DMRG calculations for an eight-site ring

As a test of generality of the ring current oscillations, we have performed additional calculations for an asymmetric ring consisting of eight lattice sites. The corresponding Hamiltonian reads

[TABLE]

and

[TABLE]

where $\langle i,j\rangle$ again denotes neighboring sites, and the gate potential $\varepsilon_{T}$ is now applied to the site with index $j=4$. A sketch of the ring is shown in figure 6 (a). We plot the low-energy spectrum of the uncoupled ring, which was obtained by means of exact diagonalization, as a function of the interaction strength $U/J$ in figure. 6 (b). Once again, we find a large separation in energy between the two lowest eigenstates ($|4,g\rangle$ and $|4,e\rangle$) and the remainder of the spectrum for $U/\varepsilon_{T}\gg 1$. The particular eigenstates again correspond to CDWs at half-filling, namely $N=4$. We have performed several td-DMRG calculations in the same fashion as for the four site ring. We have chosen parameters for interaction strength $U$ and coupling $J_{\text{c}}$ as well as lead sizes $\ell_{\text{L},\text{R}}=(L-8)/2$, for which we have previously observed slowly decaying ring current oscillations in the four site ring.

We show the results for the time-dependent currents that have been measured in the td-DMRG calculations in figure 7. In each calculation we observe an overall increase in transient features, both for the transmitted current and the ring current $I_{\text{r}}=(I_{u}-I_{l})/2$, where $I_{u}=2\text{e}\,\text{Im}(d^{\dagger}_{4}d_{2})$ and $I_{l}=2\text{e}\,\text{Im}(d^{\dagger}_{5}d_{3})$. In case of the ring current, these transient features are modulated on top of a single dominant oscillation. These modulations increase the difficulty of fitting an exponentially decaying cosine function to the data, such that the obtained values are less reliable than in the case of the four site ring. In figure 7 (a) we find a noticable decay $\Gamma/J\simeq 1/100$ of the ring current oscillation within the simulation time for $U/J=3/2$. For interaction strength $U/J=2$ (see figure 7 (b)), the decay rate diminishes further to $\Gamma/J\simeq 1/200$. In the case $U/J=3$, shown in figure 7 (c), other transient features of higher frequency have become significantly more pronounced. One can nevertheless still observe the underlying oscillation of frequency $\omega\simeq E_{4,e}-E_{4,g}$. A fit of an exponential decay to the data is however no longer feasible due to the other transient features. We show the fitted values for the oscillation frequencies in figure 6 (b). For $U/J=2$ we find very good agreement. The deviation for the other values of the interaction strength, can be explained with the deficiencies of the fitting procedure.

III Perturbation theory in the limit of small hybridization

For our computation of the decay rate $\Gamma$, which is associated with the oscillation of the local currents in the ring, we make use of the reduced density-operator transport theory (RDTT). We mainly follow Schoeller, Eur. Phys. J. Special Topics 168, 179 (2009). The RDTT approach is exact in the Hilbert space $\mathcal{H}_{\text{r}}$ of the impurity and perturbative in the hybridization between impurity and reservoirs. It may be applied if the associated perturbative scale $\Gamma_{0}=2\pi\rho_{0}J_{\text{c}}^{2}$ satisfies $\Gamma_{0}\ll T$, where $T$ denotes the system temperature. The RDTT determines the time-dependent reduced density-matrix $\rho_{\text{ns}}(t)$ of an impurity by calculating corrections to the Liouvillian $L$ of the impurity caused by tunneling processes between impurity and leads. The Liouvillian can be understood as a superoperator that corresponds to the action of the commutator between the Hamiltonian $H$ and a second operator $A\in\mathcal{H}$,

[TABLE]

The von Neumann equation, which governs the time-evolution of the density matrix $\rho$, can be written in terms of the Liouvillian as

[TABLE]

and is in turn solved by

[TABLE]

A Laplace transform and a subsequent trace over the reservoir degrees of freedom of solution (17) yields the expression

[TABLE]

where $L_{\text{ns}}$ denotes the original Liouvillian of the impurity and $L_{\text{l}}$ the Liouvillian of the reservoirs, i.e., the total Liouvillian is decomposed as $L=L_{\text{l}}+L_{\text{ns}}+L_{V}$. Expression (18) can conveniently be expanded in powers of $L_{V}$, the contribution to the Liouvillian containing the coupling between the impurity and the reservoirs. The resulting series expansion for $\tilde{\rho}_{\text{ns}}$ reads

[TABLE]

In the limit $J_{\text{c}}^{2}\ll T$ we can set up a perturbation theory in $L_{V}$. The effective Liouvillian $L_{\text{eff}}$ of the impurity then obtains perturbative corrections $\Sigma(E)$ that are functions of the Laplace variable $E$. It reads

[TABLE]

The transient features of the reduced density matrix $\rho_{\text{ns}}(t)$ are encoded in non-zero poles of

[TABLE]

To obtain these poles we solve for the complex roots of

[TABLE]

The Laplace variables $z^{*}_{\pm}$, that are roots of equation (22), have a real part that corresponds to an oscillation frequency $\varepsilon$ of the associated transient feature and an imaginary part denoting its decay rate $\Gamma$. Our aim is to compute the particular $\Gamma$ of the transient features whose frequency coincide with the frequency $\varepsilon$ of the oscillation of the currents in the ring shaped impurity.

Properties of the Liouville space

In order to represent the Liouvillians $L_{\text{ns}}$ and $L_{V}$ as well as other superoperators $G$ as matrices we introduce a new vector space $\mathcal{L}$ that we refer to as Liouville space. Objects that act as matrices in the Hilbert space $\mathcal{H}_{\text{r}}$ of the impurity can be thought of as vectors in this Liouville space $\mathcal{L}$. The most relevant example of such an object is the reduced density matrix $\rho_{\text{ns}}$ of the impurity. Each matrix element $(\rho_{\text{ns}})_{i,j}=|i\rangle\langle j|$ of $\rho_{\text{ns}}$ corresponds to a basis vector $|m)$ of the Liouville space $\mathcal{L}$. We will subsequently denote vectors in $\mathcal{H}_{\text{r}}$ as $|i\rangle$ and vectors in $\mathcal{L}$ as $|j)$. To represent each element of an operator $O\in\mathcal{H}_{\text{r}}$ as a basis vector of $\mathcal{L}$, the size of the vector space $\mathcal{L}$ has to be chosen such that $\text{dim}(\mathcal{L})=\text{dim}(\mathcal{H}_{\text{r}})^{2}$.

Definition of the superoperators

The coupling Liovilliain $L_{V}$, which can be interpreted as the interaction vertex of the perturbation theory, induces charge fluctuations on the impurity. It has the form

[TABLE]

where $G_{1}^{p_{1}}$ denotes the superoperator acting on the impurity and $:J_{1}^{p_{1}}:$ the normal ordered field superoperator acting on the reservoirs. The reservoir field superoperator is defined by its action on operators $A$ acting in the reservoir Hilbert space and reads

[TABLE]

where $1\equiv\eta,\nu,\omega$ is a collection of indices classifying the field operator $c_{1}$ such that

[TABLE]

Similarly we define $\bar{1}\equiv-\eta,\nu,\omega$. The action of the impurity vertex superoperator on this specific eigenvector is given by

[TABLE]

The index $p$, that appears in the definition of both superoperators, determines whether the respective field operator acts on the second operator $A$ from the left ($p=+$) or from the right ($p=-$). It can be interpreted as indicating the position of the field operator on the Keldysh contour and is thus sometimes referred to as Keldysh index. The operator $\sigma^{p}$ accounts for fermionic sign factors. It returns a negative sign if

[TABLE]

such that

[TABLE]

Reservoir contractions

We perform the trace $\text{tr}_{\text{l}}$ over the lead degrees of freedom by contracting pairs of reservoir field superoperators in our series expansion of $\tilde{\rho}_{\text{ns}}(E)$. We denote these contractions

[TABLE]

where $\langle\dots\rangle_{\text{eq.}}$ indicates that we assume the semi-infinite reservoirs to be in thermal equilibrium. The contractions are thus proportional to the equilibrium distribution function $f(\omega)$ at temperature $T$. We can simplify the subsequent calculations by separating the distribution function $f(\omega)$ into a symmetric and an antisymmetric contribution. The reservoir contraction then reads

[TABLE]

with the symmetric contribution

[TABLE]

and the antisymmetric contribution

[TABLE]

where $\rho_{0}$ is the density of states in the reservoir. It is possible to absorb the Keldysh index appearing in the contraction (36) by introducing the vertices

[TABLE]

Definition of the perturbative corrections

The leading order correction $\Sigma^{(1)}(E)$ to the effective impurity Liouvillian $L_{\text{eff}}$, which derives from charge fluctuations, has the form

[TABLE]

As with the reservoir contractions we can separate $\Sigma^{(1)}(E)$ in a symmetric and an antisymmetric term,

[TABLE]

where $\Sigma_{s}$ does not depend on the Laplace variable $E$. When using the redefined vertices $\bar{G}_{1}$ and $\tilde{G}_{1}$ we can write $\Sigma_{s}$ as

[TABLE]

where we have integrated over all reservoir frequencies $\omega_{1}$ ranging from the lower to the upper reservoir band edge $D$. The symmetric contribution $\Sigma_{s}$ turns out to be entirely imaginary. It thus adds only to the decay rate of transient features but not to their oscillation frequency. The antisymmetric contribution $\Sigma_{a}(E)$ is a function of the Laplace variable. It reads

[TABLE]

where $\psi(x)=\partial_{x}\log(\Gamma(x))$ is the Digamma function and $|v_{j})$ are the eigenvectors of the initial impurity Liouvillian $L_{\text{ns}}$ associated with the eigenvalues $\lambda_{j}$ of $L_{\text{ns}}$. The imaginary part of $\Sigma_{a}(E)$, which is the part contributing to the decay rate, takes a more simple, intuitive form. It reads

[TABLE]

In the basis spanned by the eigenvectors $|l\rangle$ of the impurity Hamiltonian $H_{\text{r}}$ the initial impurity Liouvillian $L_{\text{ns}}$ is diagonal as well and one can easily establish a one-to-one correspondence between an eigenvector $|v_{j})$ of $L_{\text{ns}}$ and a matrix element of $\rho_{\text{ns}}$ in this eigenbasis through

[TABLE]

with the associated eigenvalue

[TABLE]

where $\lambda_{j}$ is the energy difference between the two eigenstates $|l\rangle$ and $|m\rangle$ of the Hamiltonian $H_{\text{r}}$. There are two eigenvalues $\lambda_{\varepsilon,\pm}$ of the impurity Liovilliain $L_{\text{ns}}$ that correspond to the energy difference between the two charge density wave eigenstates $|g\rangle$ and $|e\rangle$. We denote the eigenvector that corresponds to the positive eigenvalue $\lambda_{\varepsilon,+}$ as

[TABLE]

The action of the impurity vertex superoperators on this eigenvector is given by

[TABLE]

where $d_{1}$ creates or annihilates a particle on lattice sites $x=1$ or $x=4$ of the impurity.

Perturbative diagonalization of $L_{\text{eff}}(E)$

While $L_{\text{ns}}$ is diagonal in the eigenbasis of $H_{\text{r}}$, the corrections $\Sigma_{s}$ and $\Sigma_{a}(E)$ are not. Due to the large size of the Liouville space, $\text{dim}(\mathcal{L})=256$, an analytical diagonalization of the effective Liouvillian $L_{\text{eff}}(E)=L_{\text{ns}}+\Sigma_{s}+\Sigma_{a}(E)$ is not feasible. To determine the eigenvalues of $L_{\text{eff}}(E)$ we therefore treat $\Sigma_{s}+\Sigma_{a}(E)$ as perturbations to the initial Liouvillian $L_{\text{ns}}$ and calculate the leading order corrections to its eigenvalues $\lambda_{j}$. This approximation is reasonable because $\|\Sigma_{s}+\Sigma_{a}(E)\|\leq J_{\text{c}}^{2}\ll\varepsilon_{T}\approx\lambda_{j}$. The eigenvalue corresponding to $|v_{\varepsilon})$ is then given by

[TABLE]

The particle number symmetry of the impurity Hamiltonian ($[H_{\text{r}},\sum_{j}n_{j}]=0$) guarantees that $\langle 2,e|d_{1}|2,e\rangle\equiv 0\equiv\langle 2,g|d_{1}|2,g\rangle$. We therefore find

[TABLE]

Using (52) and (53) the perturbation theory corrections from the symmetric contribution $\Sigma_{s}$ reduce to

[TABLE]

We see that the symmetric contribution from the leading order tunneling processes between reservoirs and impurity causes a decay rate $\Gamma$ equal to the typical decay rate $\Gamma_{0}$ of transient features. However, this contribution does not yet factor in the fermion distribution function $f(\omega)$ in the reservoirs, meaning that each tunneling process is treated equally. The information about the distribution function is encoded in the antisymmetric correction $\Sigma_{a}(E)$.

Matrix elements of the antisymmetric contribution

$\Sigma_{a}(E)$

The evaluation of the antisymmetric corrections is more involved as the contribution from each intermediate eigenstate $|v_{j})$ of the Liouvillian $L_{\text{ns}}$ is individually weighted by $\psi[1/2-i(E\pm V-\lambda_{j})/(2\pi T)]$. First, we identify the intermediate states $|v_{j})$ that feature in the finite contributions

[TABLE]

The impurity vertex superoperator $G^{\pm}_{1}$ either creates or annihilates a fermion on the impurity. Finite contributions thus only involve eigenstates $|v_{j})$ which satisfy $|v_{j})=|m\rangle\langle 2,g|$ or $|v_{j})=|2,e\rangle\langle m|$ such that $\sum_{i=1}^{4}n_{i}|m\rangle\in\{1,3\}$. One finds in total $N=8$ eigenstates in the Hilbert space $\mathcal{H}_{\text{r}}$ with particle number $n=1$ or $n=3$, implying $N_{L}\leq 16$ finite matrix elements. A quantitative study of the matrix elements $(G^{\pm}_{1})_{j\varepsilon}$ reveals large contributions $|(v_{j}|G^{\pm}_{1}|v_{\varepsilon})|^{2}\simeq J_{\text{c}}^{2}$ for two eigenstates $|v_{j})\in\mathcal{L}$. The two particular eigenstates are

[TABLE]

where $|1,g\rangle$ and $|3,g\rangle$ are the two low-energy eigenstates of the impurity Hamiltonian that do not exhibit a charge density wave character and for which the energy, $E_{1}=\langle 1,g|H_{\text{r}}|1,g\rangle<\langle 3,g|H_{\text{r}}|3,g\rangle=E_{3}$, is plotted in figure 2. The matrix elements read

[TABLE]

where we note that the difference in Keldysh index $p=\pm$ of the vertex superoperators $G^{p}_{1}$ causes the opposite sign of the matrix elements. An analysis of the remaining matrix elements reveals a third matrix element that gives a sizable contribution to the self energy. Here, the intermediate state is

[TABLE]

and the matrix element evaluates to

[TABLE]

The contribution from this matrix element becomes particularly relevant in the vicinity of $U=\varepsilon_{T}$ due to the analytic structure of its associated weight function.

Weight function

The weight function $f(\omega)$ contains the information about the fermionic distribution function $n(\omega)$ in the leads, which details the probability for an eigenstate of the lead Hmailtonian $H_{\text{l}}$ with energy $\omega$ to be occupied by a fermion. Our aim is to determine the decay rate $\Gamma$ which directly corresponds to the imaginary part of the the eigenvalue $\lambda_{\varepsilon}$ of the effective Liouvillian that satisfies

[TABLE]

The imaginary part of $\lambda_{\varepsilon}$ originates entirely from the imaginary part of the self-energy $\Sigma(\omega)$ correction, which for the asymmetric correction $\Sigma_{a}(\omega)$ stems from the weight function $f_{j}^{\pm}(\omega)$. The imaginary part of the weight function evaluated at the eigenvalue $\varepsilon^{*}$ has the simple form

[TABLE]

We know that $\varepsilon^{*}\simeq E_{e}-E_{g}$. We can thus also express the weigth function as

[TABLE]

The eigenvalues $\lambda_{j}$ that correspond to the three largest matrix elements read

[TABLE]

It is easy to see that, depending on the eigenvalue $\lambda_{j}$, either the dependence on $E_{g}$ or $E_{e}$ is removed from the argument of the weight function. To simplify the expression we introduce $v\equiv|V/2|$. The weight function then reads

[TABLE]

In the limit $T\ll\{J,U,\varepsilon,D\}$ we can approximate the weight function as

[TABLE]

When evaluating this approximation for the three relevant eigenvalues $\lambda_{j}^{\pm}$ one finds

[TABLE]

where we note that $\varepsilon,\Delta\geq 0$.

Decay channels

The three matrix elements $(v_{\varepsilon}|G^{-}_{1}|v^{-}_{1})(v^{-}_{1}|G^{-}_{\bar{1}}|v_{\varepsilon})$, $(v_{\varepsilon}|G^{+}_{1}|v_{3}^{+})(v_{3}^{+}|G^{+}_{\bar{1}}|v_{\varepsilon})$ and $(v_{\varepsilon}|G^{+}_{1}|v^{+}_{1})(v^{+}_{1}|G^{+}_{\bar{1}}|v_{\varepsilon})$ correspond to four different decay channels that cause the decoherence of a state of the form $|\psi\rangle=\alpha|2,g\rangle+\beta|2,e\rangle$. A schematic of these decay channels is shown in figure 10. We now turn to the discussion of the decay channels and why they become suppressed for specific sets of parameters $U$, $\varepsilon_{T}$ and $V$.

$(v_{\varepsilon}|G^{-}_{1}|v^{-}_{1})(v^{-}_{1}|G^{-}_{\bar{1}}|v_{\varepsilon})$:

An electron is ejected from the ring impurity, which has initially been in the ground state $|2,g\rangle$. Due to the particle hole symmetry of the repulsive nearest neighbor interaction $U$, this requires the energy $\Delta(U)$. Said energy needs to be supplied by the increase of chemical potential energy $\mu$, which the electron gains by entering the metallic lead. The process is thus only possible if $\mu=-V/2<-\Delta(U)$. Here, we have assumed that the electron can only enter the lead to which a negative chemical potential $-V/2\sum_{k}n_{k}$ was applied. The condition is reflected by the weigth function $f^{-}_{1}(\varepsilon^{*})$, which evaluates to zero if $v$ surpasses $\Delta$. Then the imaginary part of the asymmetric correction $\Sigma_{a}(\omega)$ does not compensate the constant imaginary part of the symmetric correction $\Sigma_{s}$ for this decay channel. A sketch of the decay process is shown figure 10 (1), where red in indicates the initial and green the final configuration of the decay process.

$(v_{\varepsilon}|G^{+}_{1}|v^{+}_{3})(v^{+}_{3}|G^{+}_{\bar{1}}|v_{\varepsilon})$:

An electron tunnels onto the ring impurity, which has initially been in the excited charge density wave state $|2,e\rangle$. The additional electron increases the interaction energy on the ring impurity by $\Delta(U)$. This energy has to be supplied by the additional electron. The process is thus only possible if the chemical potential in the lead that the electron originates from satisfies $\mu=V/2>\Delta(U)$. As with the previous matrix element, $f^{+}_{3}(\varepsilon^{*})$ vanishes once $v>\Delta$ such that the constant negative imaginary part of the symmetric correction $\Sigma_{s}$ is not compensated. We display a sketch of this decay channel in figure 10 (2).

The decay channels $(v_{\varepsilon}|G^{-}_{1}|v^{-}_{1})(v^{-}_{1}|G^{+}_{\bar{1}}|v_{\varepsilon})$ and $(v_{\varepsilon}|G^{+}_{1}|v^{+}_{3})(v^{+}_{3}|G^{+}_{\bar{1}}|v_{\varepsilon})$ are closely related - one involves the ground state while the other one involves the excited state - and are thus respectively allowed or suppressed in the same parameter regimes.

$(v_{\varepsilon}|G^{+}_{1}|v^{+}_{1})(v^{+}_{1}|G^{+}|v_{\varepsilon})$:

An electron tunnels out of the ring impurity, which has initially been in the excited charge density wave state $|2,e\rangle$. Depending on the interaction strength $U$, this tunneling process is energetically favorable or unfavorable. For $U<\varepsilon_{T}$ the one particle state $|1,g\rangle$ is lower in energy than $|2,e\rangle$. Since the electron can carry this excess energy it can tunnel into either lead as long as $\mu=V/2<\varepsilon-\Delta$. Having two effective decay channels, one for each lead, increases the decay rate as can be seen in sector (i) and (iii) of figure 9. The weight function reflects this as $f^{+}_{1}(\varepsilon^{*})=-i\pi\rho_{0}$, which adds to the imaginary part of the symmetric correction $\Sigma_{s}\propto-i\pi\rho_{0}$ instead of compensating for it. A schematic of this process is shown in figure 10 (3) and (4). For $U>\varepsilon_{T}$ the state $|1,g\rangle$ becomes higher in energy than $|2,e\rangle$. For an electron to tunnel out of the ring additional energy is now required. This energy needs to be provided by the increase in chemical potential energy $\mu$ that the electron gains by entering the lead. The tunneling process is thus only possible if $\mu=-V/2<-(\Delta-\varepsilon)$. We sketch this process in figure 10 (4).

For $U>\varepsilon_{T}$ and $v<\Delta-\varepsilon$ each decay channel becomes suppressed and we find

[TABLE]

The first order corrections to the imaginary part of the transient feature with oscillation frequency $\varepsilon^{*}$ therefore vanish entirely.

Discussion of the phase diagram

In figure 9 we plot the decay rate $\Gamma/\Gamma_{0}$ of the eigenvalue $\lambda_{\varepsilon}(\varepsilon^{*})$ as a function of the ratios $U/\varepsilon_{T}$ and $V/(2\varepsilon_{T})$. We identify five different sectors of these ratios in which the decay rate $\Gamma$ take different values due to the presence or absence of the previously outlined decay channels respectively. In sector (i) we find the presence of the decay channels (1), (2), (3) and (4). The decay channel (3) does not exist for many of the typical quantum dot systems. Its presence leads to a decay rate $\Gamma$ that exceeds the typical level broadening $\Gamma_{0}$. By increasing the bias voltage $V$ one crosses from sector (i) into sector (ii) where the decay channel (3) becomes suppressed as there is no remaining unoccupied state with energy $\omega=\varepsilon_{F}+(\varepsilon-\Delta)$ available in the left lead . In sector (ii) we find $\Gamma=\Gamma_{0}$. By increasing the interaction strength $U/\varepsilon_{T}$ sufficiently one passes from sector (ii) into sector (iv). The increase in interaction strength causes an increased energy gap $\Delta(U)$. As soon as $\Delta(U)>V/2$ both the decay channels (1) and (2) simultaneously become suppressed. This leads to a significant reduction of the decay rate by almost two orders of magnitude such that $\Gamma\leq\Gamma_{0}/80$. For $U/\varepsilon_{T}<1$ the reduction of the bias voltage facilitates a crossover from sector (iv) into sector (iii). In this sector, the decay channel (3) is no longer suppressed leading to small increase of $\Delta\Gamma\simeq\Gamma_{0}/80$. For $U/\varepsilon_{T}>1$ and $V/(2\varepsilon_{T})<(\Delta-\varepsilon)$ every decay channel is suppressed as is shown in sector (v). The corresponding decay rate becomes $\Gamma/\Gamma_{0}\rightarrow 0$. For a finite decay rate, higher order perturbation theory corrections would be required. However, these corrections can induce no more than a decay rate $\Gamma\propto\Gamma^{2}_{0}\ll T$.

Eigenvalue spectrum of the effective Liouvillian

The disappearance of the decay rate $\Gamma$ for an eigenvalue $\lambda(z)$ of the effective Liouvillian with finite real part is unique to the eigenvalues $\pm\varepsilon^{*}-i\Gamma$. In figure 11 we display the real and imaginary part of each root of

[TABLE]

where $\lambda_{j}(z)$ are the eigenvalues of the effective Liouvillian.

We find that the imaginary part of all but four roots is of order $\mathcal{O}(\Gamma_{0})$. This means that almost all transient features decay on the expected time scale. In figure 11 (d) we zoom in on the four extraordinary roots. We see that their imaginary part is $\text{Im}(z)\simeq\Gamma_{0}/100000$. This indicates a clear separation of scales between the decay rate of these four roots and each other root. Such a separation of scales in the solutions to (74) is sometimes referred to as dissipative phase transition. The four roots that feature a small imaginary part belong to $|2,g\rangle\langle 2,g|$ with root $z=0-i0^{+}$, $|2,e\rangle\langle 2,e|$ with root $z=0-i\Gamma_{1}$, $|2,e\rangle\langle 2,g|$ with root $z=\varepsilon-i\Gamma$, and $|2,g\rangle\langle 2,e|$ with root $z=-\varepsilon-i\Gamma$.

III.1 Inverse transformation to realtime

The effective reduced density matrix matrix of the impurity in Laplace space $\tilde{\rho}_{\text{ns}}(E)$ is given by

[TABLE]

where each pole of the resolvent

[TABLE]

has an imaginary part $\Gamma_{j}\leq 0$. We can therefore replace the inverse Laplace transform by a Fourier transform and close the integration contour in the lower half-plane such that

[TABLE]

Each pole of (75) corresponds to a transient feature with frequency $\lambda_{j}$ and decay rate $\Gamma_{j}$.

III.2 Coupling of the current operators to the pole

$\lambda_{\varepsilon}$

Ring current

The operator $I_{\text{r}}$ measuring the local current in the ring reads

[TABLE]

with

[TABLE]

After a transformation of the current operator and the reduced density matrix $\rho_{\text{ns}}$ to the basis of the eigenstates of Hamiltonian, this becomes

[TABLE]

where $\rho_{\text{ns},n}(t)$ denotes the time-dependent reduced density matrix expressed in the basis given by the eigenstates of the Hamiltonian. We are mainly interested in the matrix element of $\tilde{I}_{r}$ that couples to the matrix element $|v_{\varepsilon})=|2,e\rangle\langle 2,g|$ of the reduced density matrix

[TABLE]

In figure 12 (a) we plot the absolute value of the coupling of the current operator to the off-diagonal matrix element $|e\rangle\langle g|$ of the reduced density matrix. We find that after an initial increase with interaction strength, the coupling decreases with interaction strength. In the entire range of values for the interaction strength that we have studied, the matrix element $(\tilde{I}_{r})_{g,e}$ exceed every other matrix element of the current operator $\tilde{I}_{r}$.

Transmitted current

We determine the extent to which the operator $I_{\text{t}}$, measuring the transmitted current, couples to the matrix elements $|2,e\rangle\langle 2,g|$ and $|2,g\rangle\langle 2,e|$ of the reduced density matrix $\tilde{\rho}_{S}(E)$ directly from the perturbation theory. The expectation value of the transmitted current in Laplace space is given by

[TABLE]

where

[TABLE]

The modified vertex superoperator is defined as

[TABLE]

where the operator $I_{\text{t}}$ acts in the Hilbert space $\mathcal{H}_{\text{r}}$ of the ring as

[TABLE]

The part of the transmitted current that acquires the small decay rate $\Gamma$ is given by

[TABLE]

where the vectors $|n\rangle$ form a basis of the Hilbert space of the ring. In figure 12 (b) we plot the real and imaginary part of $(I_{\text{t}})_{e,g}$ as a function of $U/\varepsilon_{T}$. We find that the real part of $(I_{\text{t}})_{e,g}$ is small but features a reasonance at $U=\varepsilon_{T}$. We thus find that only a small part of the transmitted current $I_{\text{t}}$ decays with the decay rate $\Gamma$ while the majority relaxes with the decay rate $\Gamma_{0}$.

IV Schrieffer-Wolff transformation and perturbation theory

IV.1 Schrieffer-Wolff transformation of the impurity

system

In the limit $U/\varepsilon_{T}\rightarrow\infty$ the low-energy sector of the spectrum of the uncoupled ring impurity features only the two charge density wave eigenstates $|2,g\rangle$ and $|2,e\rangle$. From figure 2 it becomes obvious that for $U/J\gg 1$ the energy gap $\varepsilon$ between the two CDW eigenstates becomes small compared to the energy separation between the CDW states and the remainder of the spectrum. It is then intuitive to construct an effective low-energy Hamiltonian in the subspace of the Hilbert space $\mathcal{H}_{\text{r}}$, which is spanned by the two CDW eigenstates. In the limit $U/\varepsilon_{T}\rightarrow\infty$ and $U/J\rightarrow\infty$ the CDW eigenstates take the form of simple product states $|2,g\rangle=|1\rangle|0\rangle|0\rangle|1\rangle\equiv|1001\rangle$ and $|2,e\rangle=|0\rangle|1\rangle|1\rangle|0\rangle\equiv|0110\rangle$. We define the operator that projects onto this low-energy subspace as

[TABLE]

The Hamiltonian of the full system can be separated into a contribution that is diagonal in this new basis

[TABLE]

and a contribution that connects the subspace spanned by the CDW eigenstates with the rest of the Hilbert space, which reads

[TABLE]

The second contribution can be regarded as a small perturbation. We then perform a Schrieffer-Wolff transformation to project onto the subspace spanned by the CDW states and to virtually include transitions to states orthogonal to the CDW eigenstates. We follow the work by Bravyi, DiVincenzo and Loss [Ann. Phys. 326, 2793 (2011)] to expand the general expression for a transformation of the Hamiltonian

[TABLE]

where $\exp(S)$ is a unitary operator, into a power series up to fourth order in the perturbation $\hat{V}=V_{\text{od}}+V_{\text{d}}$ reading

[TABLE]

where $V_{\text{od}}$ denotes the part of the perturbation $\hat{V}$ that facilitates transitions between the low-energy Hilbert and the complementary Hilbert space and $V_{\text{d}}$ denotes the part of the perturbation that only connects states exclusively inside either subspace. We use the shorthand

[TABLE]

to denote the inverse energy difference between two states $|i\rangle$ and $|j\rangle$ that are connected through the operator $X$. After evaluation of the commutators we arrive at the expression for the effective Hamiltonian $H_{\text{eff}}^{(4)}$. It reads

[TABLE]

where the first line includes all contributions up to second order in the perturbation $\hat{V}$ and lines two through five contain the contributions up to fourth order.

Second order correction

Evaluating the diagonal contribution $H_{0}$ for the two CDW eigenstates yields the effective Hamiltonian in zeroth order as

[TABLE]

where from now on we treat the two CDW states like pseudo-spins defined as

[TABLE]

The leading order corrections to the effective Hamiltonian are of second order in the perturbation $\hat{V}$. The corrections encompass two consecutive tunneling processes, either tunneling within the ring impurity leading to corrections $\propto J^{2}$, or tunneling between the ring and the leads and back yielding corrections $\propto J^{2}_{\text{c}}$. In the follwing we show the calculation of each correction term featured in

[TABLE]

The first leading order correction term gives

[TABLE]

where $n_{\text{L},0}$ and $n_{\text{R},0}$ is the electron density on the site of the lead closest to the impurity for the left and the right lead respectively. The second term of the leading order correction yields

[TABLE]

Assuming that the mean electron density in both leads combined is $n_{\text{L},x}+n_{\text{R},x}=1$, the two correction terms are identical. The effective Hamiltonian in leading order then reads

[TABLE]

We see that the effective energy gap $\varepsilon^{(2)}$ between the ground state $|\downarrow\rangle$ and the excited state $|\uparrow\rangle$ is reduced as compared to the bare energy gap $\varepsilon_{T}$ by the perturbative corrections. In leading order we furthermore find no off-diagonal terms and as it turns out not in any higher order $(J^{2})^{n}$ of perturbations $\hat{V}$ which feature only in-ring hopping terms $\sim Jd_{i}^{\dagger}d_{j}$. Since the hopping between ring and leads alone cannot facilitate a pseudo-spin flip, one finds that they are not possible in leading order.

Mirror symmetry

We attribute the lack of the off-diagonal, pseudo-spin flip terms to a symmetry of the nanostructure associated with the mirror symmetry in the axis through lattice sites $2$ and $3$ or equivalently an exchange of lattice sites $1\leftrightarrow 4$. The operator $M$ corresponding to this symmetry reads

[TABLE]

The symmetry operator satisfies

[TABLE]

as well as

[TABLE]

and

[TABLE]

The CDW eigenstates of the Hamiltonian for $U/J\rightarrow\infty$ are also eigenstates of the symmetry operator $M$ with eigenvalues $m=\pm 1$. Since $M$ commutes with the Hamiltonian $H_{\text{r}}$ of the uncoupled impurity for arbitrary $U/J$, the in-ring hopping terms of $H_{\text{r}}$ cannot couple the different eigenstates of $M$. For the Hamiltonian $H_{\text{c}}$, which couples the ring to the leads, we instead find $\left[H_{\text{c}},M\right]\neq 0$. As a consequence we have $\left[(H_{\text{c}}+H_{\text{r}})^{2},M\right]\neq 0$, indicating that off-diagonal, pseudo-spin flip terms can occur in higher orders of the perturbation. We illustrate the connection between the mirror symmetry and the absence of off-diagonal terms in leading order in figure 13. Each process connecting $|\downarrow\rangle$ and $|\uparrow\rangle$ has a mirror image with opposite sign leading to pairwise cancellation of all terms.

Fourth order correction

In order to obtain finite off-diagonal terms in the effective Hamiltonian one needs to include the fourth order corrections. Here we show an example calculation of one such correction term. All other fourth order correction terms follow accordingly.

[TABLE]

The contributions that feature the underlined term in equation (IV.1) combine into

[TABLE]

We find that the pseudo-spins on the impurity couple to a second spin-like degree of freedom in the leads which can be associated with the symmetric and antisymmetric modes in the leads. We define the annihilation operator for an electron with pseudo-spin $\sigma$ in the leads as

[TABLE]

We use this pseudo-spin notation for the lead degrees of freedom and collect the different correction terms up to fourth order in the perturbation. The corrections read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

First, we identify the off-diagonal terms which couple $|1001\rangle=|\downarrow\rangle$ and $|0110\rangle=|\uparrow\rangle$. Now, unlike the leading order corrections, the off-diagonal terms are finite and cause a simultaneous pseudo-spin-flip on the impurity and in the leads. We can thus express the off-diagonal terms as

[TABLE]

We find that the off-diagonal terms describe a spin-spin interaction in the x and y direction with an amplitude $J_{\perp}(U,\varepsilon_{T})$. The diagonal correction terms are

[TABLE]

The diagonal correction terms contain three different couplings. A spin-spin interaction in the z-direction with amplitude $J_{z}(U,\varepsilon_{T})$, a correction to the effective magnetic field on the impurity and a small effective magnetic field $h^{*}(U,\varepsilon_{t})$ on the sites of the leads next to the impurity. The corrections expressed in the pseudo-spin notation are

[TABLE]

In the last step we collect the correction terms arising from lines three through five of (IV.1). These read

[TABLE]

[TABLE]

Adding the correction terms up to fourth order in the perturbation $\hat{V}$ to the Hamiltonian describing the lead degrees of freedom we arrive at the effective low-energy Hamiltonian $H^{(4)}_{\text{eff}}$ which reads

[TABLE]

where $h(U,\varepsilon_{T})$ constitutes an effective magnetic field on the impurity and incorporates all terms coupling to $S^{z}\otimes\mathbf{1}_{r}$, $\tilde{h}(U,\varepsilon_{T})$ denotes the terms proprotional to $\mathbf{1}_{d}\otimes S^{z}_{\text{res}}$, and $J_{z}(U,\varepsilon_{T})$ and $J_{\perp}(U,\varepsilon_{T})$ feature all terms coupling the impurity spin and the lead spins in the z-direction or x-y-direction respectively. This effective model is reminiscent of the anisotropic single-channel Kondo model with anisotropic coupling between lead spins and impurity spin in the z-direction and the x-y-plane as well as a magnetic field $h(\varepsilon_{T})\simeq\mathcal{O}(\varepsilon_{T})$ on the impurity. In our effective model the spin degrees of freedom do not correspond to physical spins. Instead, we identify pseudo-spins $|\downarrow\rangle\equiv|g\rangle$ and $|\uparrow\rangle\equiv|e\rangle$ on the impurity. For vanishing bias voltage, $V=0$, we can identify $|\downarrow\rangle$ with the antisymmetric and $|\uparrow\rangle$ with symmetric modes in the leads. The operators creating these modes read

[TABLE]

For finite bias voltage, $V\neq 0$, the term

[TABLE]

leads to a hybridization of these two modes, which means they are no longer eigenstates of the lead Hamiltonian. The conserved quantum number becomes the lead index $\alpha=L,R$ instead. The linear dependance between lead index $\alpha$ and pseudo-spin index $\sigma$ in the leads proves to be responsible for differing properties, e.g. decay rates, of our effective model as compared to the anisotropic Kondo model.

IV.2 Schrieffer-Wolff transformation of the current operators

To determine how the operator representing the ring current $I_{\text{r}}$ couples to to the matrix elements of the reduced impurity density matrix $\rho_{\text{ns}}$ we perform a second Schrieffer-Wolff transformation up to leading order in $J^{2}/U$. From this we obtain an effective ring current operator $I_{\text{eff}}$ acting in the subspace $P_{0}$. The operators measuring the local currents in the ring were previously defined as

[TABLE]

The Schrieffer-Wolff of the current operators has been performed in the same fashion as the transformation for the effective Hamiltonian. The resulting effective ring current operator in leading order reads

[TABLE]

One consequently finds that $(I_{\text{r,eff}})_{\uparrow,\downarrow}=(I_{\text{eff}})^{\dagger}_{\downarrow,\uparrow}\propto J^{2}/U$ and $(I_{\text{r,eff}})_{\downarrow,\downarrow}=(I_{\text{eff}})_{\uparrow,\uparrow}\equiv 0$. The effective current operator couples exclusively to the off-diagonal matrix elements $\rho_{\uparrow,\downarrow}$ and $\rho_{\downarrow,\uparrow}$ of the reduced impurity density matrix $\rho_{\text{ns}}$ in the context the effective low-energy decription. The transient decay rate of the ring current is therefore determined by the decay rate of these two particular matrix elements.

In contrast, the current, which is transmitted through the ring, is equivalent to

[TABLE]

as the sum of the two local currents amounts to the total current flowing from one lead through the ring to the second lead. Using the effective operators for the local currents $I_{\text{u}}$ and $I_{\text{l}}$ after a Schrieffer-Wolff transformation, we find for the transmitted current in leading order the effective operator

[TABLE]

It is immediately obvious that the effective operator for the transmitted current $I_{\text{t,eff}}\rightarrow 0$ for $U/\varepsilon_{T}\rightarrow\infty$. As a consequence, the coupling of the transmitted current to the off-diagonal elements of the effective reduced density matrix, which are the ones exhibiting the small decay rate, is strongly suppressed for strong interaction. This is consistent with our numerical data from tdDMRG, where we also do not observe a slow decay of the transmitted current, but a decay on the time scale given by the hybridization $\Gamma_{0}$.

IV.3 Limitations on the viability of the effective low-energy model

As a consistency check of the effective model $H_{\text{eff}}^{(4)}$ up to order $U^{-4}$ we perform a series expansion of the perturbative corrections in $U^{-n}$ around $U/J\rightarrow\infty$. The results of this series expansion for the amplitudes of the spin-spin interaction terms read

[TABLE]

When expanding the expression for $J_{\perp}$ up to fourth order in the inverse interaction strength $U^{-1}$ we encounter an inconsistency of our Schrieffer-Wolff transformation as $J_{\perp}$ vanishes up to this order. We can therefore not assume with certainty that the spin-spin interaction $J_{\perp}$ in the effective model is finite.

IV.4 Perturbation theory for the effective model

in the limit $T\rightarrow 0$

In our effective model the hybridization between the impurity and the leads satisfies

[TABLE]

which is small in the limit $U\gg\text{max}(J,\varepsilon)$ even in the case of a chosen bare coupling $J_{\text{c}}=\mathcal{O}(J)$. In the limit $U/\varepsilon\gg 1$ we can thus perform a perturbation theory calculation for the effective low-energy model whilst employing the same values for the bare model parameters $J$, $J_{\text{c}}$ as in our initial DMRG calculations. This way we can compare the results from both methods for the decay rate of the off-diagonal matrix elements $\rho_{\uparrow,\downarrow}$ and $\rho_{\downarrow,\uparrow}$ and thus the decay rate of the transient ring current in the strong interaction limit. In contrast to the earlier perturbation theory calculation, we no longer study charge fluctuations on the impurity but pseudo-spin fluctuations instead. This requires a few modifications to the procedure outlined in the previous section on the perturbation theory. The coupling Liovilliain $L_{V}$ now features two field superoperators for the leads and the impurity, instead of just one. It reads

[TABLE]

where

[TABLE]

The field superoperators for the leads remain unchanged

[TABLE]

As a consequence, the first order corrections to the Liouvillian now contain two reservoir contractions $\gamma^{pp^{\prime}}_{11^{\prime}}$ which in turn requires integration over two reservoir frequencies $\omega_{1}$ and $\omega_{2}$. The perturbative correction reads

[TABLE]

It is again possible to separate the corrections into a symmetric and an antisymmetric part. In the zero temperature limit the two contributions read

[TABLE]

First we discuss the integrals over the reservoir frequencies $\omega_{1}$ and $\omega_{2}$. For this we introduce the density of states

[TABLE]

The integral in the symmetric part of the self-energy correction evaluates to

[TABLE]

Similarly, the integral in the antisymmetric self-energy correction gives

[TABLE]

Next, we discuss the superoperators acting in the Liouville space of the impurity. We follow the notation introduced in Schoeller and Reininghaus [Phys. Rev. B 80, 045117 (2009)]. First we define the Liouville superoperators that act as the spin operators $\underline{S}=(S^{x},S^{y},S^{z})$ on the impurity. These Liouville superoperators are

[TABLE]

where the sign $p=\pm$ indicates the order of the operators as

[TABLE]

A matrix representation of these superoperators in the basis $|\uparrow\uparrow)$, $|\downarrow\downarrow)$, $|\uparrow\downarrow)$, $|\downarrow\uparrow)$ reads

[TABLE]

[TABLE]

From these superoperators we can construct a basis of superoperators sufficient to describe the spin-spin interaction processes between impurity and lead pseudo-spins. We further introduce the ’scalar’ superoperators

[TABLE]

as well as the vector superoperators

[TABLE]

Due to the anisotropy of the interactions we need to introduce a third set of superoperators, which reads

[TABLE]

In terms of these basis superoperators the bare Liouvillian is given as

[TABLE]

where $h(U,\varepsilon_{T})$ is the effective magnetic field on the impurity and $L^{h}$ represents the action of $[S^{z},\bullet]$ on the impurity. The effective Liouvillian for the impurity in first order perturbation theory reads

[TABLE]

The first order, energy-dependent self-energy corrections are

[TABLE]

where $\mu_{12}=\eta_{1}\mu_{1}+\eta_{2}\mu_{2}$ and

[TABLE]

as well as

[TABLE]

It is then straightforward to calculate the self energy corrections that are proportional to simple products $\bar{G}_{12}\bar{G}_{\bar{2}\bar{1}}$ and $\bar{G}_{12}\tilde{G}_{\bar{2}\bar{1}}$. The two different products of superoperators evaluate to

[TABLE]

and

[TABLE]

To calculate terms involving $\bar{G}f(\mu_{12})\bar{G}$ and $\bar{G}f(\mu_{12})\tilde{G}$ a rotation to a different basis is necessary.

IV.5 Modifications to the perturbation theory for spin

fluctuations due to the linear dependence between pseudo-spin and lead index

From eq. (IV.1) we see that the magnetic field $\tilde{h}(U,\varepsilon_{T})$, experienced by the pseudo-spins on the lead sites closest to the impurity, is small up to order $\mathcal{O}(U^{-4})$, ie $h^{*}\ll h,V$ and we can thus neglect it. It is then more useful to express the Hamiltonian in the basis of the lead index $\alpha\in\{L,R\}$ eigenstates, which corresponds to a rotation $\tau^{x}_{\sigma_{1}\sigma_{2}}\rightarrow\tau^{z}_{\alpha_{1}\alpha_{2}}$, $\tau^{y}_{\sigma_{1}\sigma_{2}}\rightarrow-\tau^{y}_{\alpha_{1}\alpha_{2}}$, $\tau^{z}_{\sigma_{1}\sigma_{2}}\rightarrow\tau^{x}_{\alpha_{1}\alpha_{2}}$ in the leads. In terms of the rotated operators, $c_{\alpha,k}=1/\sqrt{2}(c_{\downarrow,k}\pm c_{\uparrow,k})$, the effective Hamiltonian reads

[TABLE]

In this rotated basis the part of the Hamiltonian acting exclusively on the leads is diagonal so the reservoir contractions reduce to simple fermionic distribution functions. After rotation the vertex superoperators read

[TABLE]

and

[TABLE]

where we have dropped the factor $1/2$ resulting from the substitution $S_{\sigma_{1}\sigma_{2}}\rightarrow\tau_{\alpha_{1}\alpha_{2}}$ for convenience. We reintroduce the factor in the final result. The first set of self-energy corrections that are affected by the linear dependance between pseudo-spin index and lead index involve terms proportional to $\bar{G}_{12}\mu_{12}\bar{G}_{\bar{2}\bar{1}}$ and $\bar{G}_{12}\mu_{12}\tilde{G}_{\bar{2}\bar{1}}$. The first term reads

[TABLE]

In the following we evaluate the cases $l=k$ and $l\neq k$ separately. For each example calculation we set $\eta_{1}=-\eta_{2}=+$ without loss of generality. For $l=k$ we find

[TABLE]

and for $l\neq k$ we obtain

[TABLE]

This term, proportional to the bias voltage $V$, does not appear in the perturbation theory of the regular anisotropic Kondo model. The correction term still satisfies $\text{Tr}_{S}(L^{2x})=0$ such that $\text{Tr}_{S}(L_{\text{eff}})=0$, a necessary requirement for the validity of the perturbation theory. Similarly for $\bar{G}_{12}\mu_{12}\tilde{G}_{\bar{2}\bar{1}}$ we find

[TABLE]

which satisfies $\text{Tr}_{S}(L^{3x})=0$ as well. Next we discuss the correction terms proportional to $\bar{G}_{12}|E+\mu_{12}-L_{\text{ns}}|\bar{G}_{\bar{2}\bar{1}}$. We know that $\Gamma_{0}\ll h$ which means that the perturbative corrections to the roots of the unperturbed Liovillian $L_{\text{ns}}$ are small. We can thus safely assume $\lambda_{\pm}^{*}=L_{\text{eff}}(\lambda_{\pm}^{*})\simeq\pm h$. As such we evaluate the correction terms proportional to $\bar{G}_{12}|E+\mu_{12}-L_{\text{ns}}|\bar{G}_{\bar{2}\bar{1}}$ for $E\simeq\pm h$ and obtain

[TABLE]

where $\delta h=h-h_{0}$ and we verify that $\text{Tr}_{S}(L^{4}_{\pm})=\text{Tr}_{S}(L^{5}_{\pm})=\text{Tr}_{S}(L^{6}_{\pm})=0$. Lastly we discuss the correction term that involves the logarithm of the Liouvillian $L_{\text{ns}}$. We abbreviate $z=E+\mu_{12}-L_{\text{ns}}$ and approximate $E=\pm h$. The correction term then evaluates to

[TABLE]

With all the self-energy terms evaluated we can determine the eigenvalues of the effective Liouvillian $L_{\text{eff}}(E)$. To obtain an analytical result for the roots $\pm h$ we perform the diagonalization of $L_{\text{eff}}$ perturbatively as well. In first order

[TABLE]

we find

[TABLE]

where $\delta h=h-h_{0}\simeq 0$ and $h_{0}$ denotes the root of the bare Liouvillian $L_{\text{ns}}$. The imaginary part of the root $h$, which corresponds to its transient decay rate, reads

[TABLE]

We see that for $V\rightarrow 0$ the imaginary part of $h$ is essentially given by the terms proportional to $J_{\perp}^{2}$. For $V=0$ we thus find a power law decrease of the decay rate with $U^{-\alpha}$ and $\alpha\geq 8$. For finite bias voltage and $U\rightarrow\infty$ the terms proportional to $J_{z}^{2}$ become dominant and we observe a power law decrease of the decay rate with $U^{-\beta}$ and $\beta=6$.

In figure 14 we plot our numerical results for the decay rates $\text{Im}(\lambda^{*}_{\pm})$ with $\lambda^{*}_{\pm}=\pm h$. We find that for $V=0$ the decay rates obey a power law, $\text{Im}(\lambda^{*}_{\pm})(U)\propto U^{-8}$, the same as the spin-flip interaction $J_{\perp}^{2}(U)$. For finite bias voltage we observe a different power law, $\text{Im}(\lambda^{*}_{\pm})(U)\propto U^{-6}$, a behavior shared by $J_{z}^{2}(U)$. Our numerical findings support our perturbative result for the decay rates (187). The perturbation theory for the effective model finds that $\Gamma\rightarrow 0$ for $U\rightarrow\infty$ and supports our findings from DMRG calculations and perturbation theory in the limit of small coupling which see very long life times $\tau\gg\Gamma_{0}^{-1}$ of the ring current oscillations.