Quantum Impurity Models coupled to Markovian and Non Markovian Baths

TL;DR

This paper introduces a new method for analyzing quantum impurity models with both Markovian and non-Markovian baths, enabling the calculation of the impurity's reduced density matrix evolution in complex environments.

Contribution

The authors develop an exact real-time hybridization expansion and a Dyson equation approach with Non-Crossing-Approximation for impurity models coupled to diverse baths, including Markovian dissipation.

Findings

Successfully applied to a fermionic impurity with Markovian pump, losses, and dephasing.

Provides a framework for stochastic sampling via diagrammatic Monte Carlo.

Extends existing methods to include non-linear couplings and Markovian environments.

Abstract

We develop a method to study quantum impurity models, small interacting quantum systems linearly coupled to an environment, in presence of an additional Markovian quantum bath, with a generic non-linear coupling to the impurity. We aim at computing the evolution operator of the reduced density matrix of the impurity, obtained after tracing out all the environmental degrees of freedom. First, we derive an exact real-time hybridization expansion for this quantity, which generalizes the result obtained in absence of the additional Markovian dissipation, and which could be amenable to stochastic sampling through diagrammatic Monte Carlo. Then, we obtain a Dyson equation for this quantity and we evaluate its self-energy with a resummation technique known as the Non-Crossing-Approximation. We apply this novel approach to a simple fermionic impurity coupled to a zero temperature fermionic bath and in presence of Markovian pump, losses and dephasing.

Equations

$$H_M$$

$$H=H_I+H_M+H_{IM}+H_{\bar{M}}+H_{I\bar{M}}$$

$$H_{I\bar{M}}=\sum _{ka}{V}_{ka}\left(d_a^{†}{c}_k+\text{hc}\right)$$

$$H_{IM}=\sum _{\alpha }{X}_{\alpha }{B}_{\alpha }$$

$${\rho }_I(t)={\text{tr}}_{M\bar{M}}\left[U(t,0)\rho (0)U^{†}(t,0)\right]$$

$${\rho }_I(t)=\mathcal{V}(t,0){\rho }_I(0)$$

$$\begin{array}{rl}\rho (t)=& e^{-i\left(H_0+H_{\bar{M}}\right)t}{T}_t{e}^{-i{\int }_0^{t}d{t}^{\prime }{H}_{I\bar{M}}(t^{\prime })}\rho (0)\times \\ \times & {\check{T}}_t{e}^{i{\int }_0^{t}d{t}^{\prime }{H}_{I\bar{M}}(t^{\prime })}{e}^{i\left(H_0+H_{\bar{M}}\right)t}\end{array}$$

$$t_{\gamma }>t_{{\gamma }^{\prime }}^{\prime }\text{if}\{\begin{array}{ll}t>t^{\prime }& \gamma ={\gamma }^{\prime }=+\\ t<t^{\prime }& \gamma ={\gamma }^{\prime }=-\\ \gamma =-& {\gamma }^{\prime }=+\end{array}$$

$$\xi$$

$$T_C{X}_1(t_{\gamma })X_2(t_{{\gamma }^{\prime }})=\{\begin{array}{ll}{X}_1(t_{\gamma })X_2(t_{{\gamma }^{\prime }})& \text{for tγ>tγ′′}\\ \xi X_2(t_{{\gamma }^{\prime }})X_1(t_{\gamma })& \text{for tγ<tγ′′}\end{array}$$

$$\rho (t)=T_C{e}^{-i\left(H_0+H_{\bar{M}}\right)(t_{+}-t_{-})}{e}^{-i{\int }_{C(0,t)}d{t}^{\prime }{H}_{I\bar{M}}(t^{\prime })}\rho (0)$$

$$\begin{array}{r}\mathcal{V}(t,0){\rho }_I(0)={\text{tr}}_{M\bar{M}}\left[T_C{e}^{-i{H}_0(t_{+}-t_{-})}{e}^{-i{\int }_{C(0,t)}d{t}^{\prime }{H}_{I\bar{M}}(t^{\prime })}\rho (0)\right]\end{array}$$

$${\Delta }_{a^{\prime }a}^{{\gamma }^{\prime }\gamma }(t^{\prime },t)\equiv \sum _{b^{\prime }}{V}_{a^{\prime }{b}^{\prime }}{V}_{a{b}^{\prime }}^{\ast }{G}_{b^{\prime }}(t_{{\gamma }^{\prime }}^{\prime },t_{\gamma })$$

$$\begin{array}{rl}\mathcal{V}(t,0){\rho }_I(0)& =\sum _{k=0}\frac{(-i{)}^k}{k{!}^2}\sum _{{\gamma }_1\dots {\gamma }_k^{\prime }}\prod _i{\gamma }_i{\gamma }_i^{\prime }\sum _{a_1\dots a_k^{\prime }}{\int }_0^{t}d{t}_1...{\int }_0^{t}d{t}_k^{\prime }{\text{tr}}_M\left[T_C{e}^{-i{H}_0(t_{+}-t_{-})}{d}_{a_k^{\prime }}^{†}(t_k^{\prime },{\gamma }_k^{\prime })d_{a_k}(t_k,{\gamma }_k)...d_{a_1}(t_1,{\gamma }_1){\rho }_{IM}(0)\right]\times \\ & \times \sum _{\sigma \in P}{\xi }^{\text{sign}(\sigma )}{\Delta }_{a_1^{\prime }{a}_{\sigma (1)}}^{{\gamma }_1^{\prime }{\gamma }_{\sigma (1)}}(t_1^{\prime },t_{\sigma (1)})\dots {\Delta }_{a_k^{\prime }{a}_{\sigma (k)}}^{{\gamma }_k^{\prime }{\gamma }_{\sigma (k)}}(t_k^{\prime },t_{\sigma (k)})\end{array}$$

$$\rho (t)=U(t,0)\rho (0)U(0,t)\equiv \mathcal{U}(t,0)\rho (0)$$

$$d_{\gamma }^{†}\left[\bullet \right]=\{\begin{array}{ll}{d}^{†}\bullet & \text{if }\gamma =+\\ \bullet d^{†}& \text{if }\gamma =-\end{array}$$

$$\begin{array}{rl}& T_C{X}_{1(t,\gamma )}{\mathcal{U}}_0(t,t^{\prime })X_{2(t^{\prime },{\gamma }^{\prime })}=\\ & =\{\begin{array}{ll}{X}_{1(t,\gamma )}{\mathcal{U}}_0(t,t^{\prime })X_{2(t^{\prime },{\gamma }^{\prime })}& \text{for (t,γ)>(t′,γ′)}\\ \xi X_{2(t^{\prime },{\gamma }^{\prime })}{\mathcal{U}}_0(t^{\prime },t)X_{1(t,\gamma )}& \text{for (t,γ)<(t′,γ′)}\end{array}\end{array}$$

$$\begin{array}{rl}& T_F{X}_{1(t,\gamma )}{\mathcal{U}}_0(t,t^{\prime })X_{2(t^{\prime },{\gamma }^{\prime })}=\\ & =\{\begin{array}{ll}{X}_{1(t,\gamma )}{\mathcal{U}}_0(t,t^{\prime })X_{2(t^{\prime },{\gamma }^{\prime })}& \text{for t>t′}\\ X_{2(t^{\prime },{\gamma }^{\prime })}{\mathcal{U}}_0(t^{\prime },t)X_{1(t,\gamma )}& \text{for t<t′}\end{array}\end{array}$$

$$\begin{array}{rl}& T_C{e}^{-i{H}_0(t_{+}-t_{-})}{d}^{†}(t_k^{\prime })\dots d(t_1){\rho }_{IM}(0)=\\ =& T_F{T}_C{\mathcal{U}}_0(t,t_k^{\prime })d_{t_k^{\prime }{\gamma }_k^{\prime }}^{†}{\mathcal{U}}_0(t_k^{\prime },t_{k-1})\dots d_{t_1{\gamma }_1}{\mathcal{U}}_0(t_1,0){\rho }_{IM}(0)\end{array}$$

$${\partial }_t{\rho }_I^0=-i\left[H_I,{\rho }_I^0\right]+\sum _{\alpha }{\gamma }_{\alpha }\left(L_{\alpha }{\rho }_I^0{L}_{\alpha }^{†}-\frac{1}{2}\left\{L_{\alpha }^{†}{L}_{\alpha },{\rho }_I^0\right\}\right)$$

$${\rho }_I^0(t)={\mathcal{V}}_0(t-t^{\prime }){\rho }_I(0)$$

$${\text{tr}}_M{\rho }_{IM}^0(t)={\text{tr}}_M\left[{\mathcal{U}}_0(t,t^{\prime }){\rho }_{IM}^0(t^{\prime })\right]={\mathcal{V}}_0(t-t^{\prime }){\text{tr}}_M{\rho }_{IM}^0(t^{\prime })$$

$$\begin{array}{rl}& {\text{tr}}_M\left[{\mathcal{U}}_0(t,t_k^{\prime })d_{t_k^{\prime }{\gamma }_k^{\prime }}^{†}\dots d_{t_1{\gamma }_1}{\mathcal{U}}_0(t_1,0){\rho }_{IM}(0)\right]=\\ & {\text{tr}}_M{r}_1(t)={\text{tr}}_M\left[{\mathcal{U}}_0(t,t_k^{\prime })r_1(t_k^{\prime })\right]={\mathcal{V}}_0(t-t^{\prime }){\text{tr}}_M{r}_1(t_k^{\prime })\end{array}$$

$$\begin{array}{rl}\mathcal{V}(t,0)& =\sum _{k=0}\frac{(-i{)}^k}{k{!}^2}\sum _{{\gamma }_1\dots {\gamma }_k^{\prime }}\prod _i{\gamma }_i{\gamma }_i^{\prime }\sum _{a_1\dots a_k^{\prime }}{\int }_0^{t}d{t}_1\cdots {\int }_0^{t}d{t}_k^{\prime }{T}_F{T}_C{\mathcal{V}}_0(t,t_k^{\prime })d_{a_k^{\prime }(t_k^{\prime }{\gamma }_k^{\prime })}^{†}{\mathcal{V}}_0(t_k^{\prime },t_{k-1})\dots d_{a_1(t_1{\gamma }_1)}{\mathcal{V}}_0(t_1,0)\times \\ & \times \sum _{\sigma \in P}{\xi }^{\text{sign}(\sigma )}{\Delta }_{a_1^{\prime }{a}_{\sigma (1)}}^{{\gamma }_1^{\prime }{\gamma }_{\sigma (1)}}(t_1^{\prime },t_{\sigma (1)})\dots {\Delta }_{a_k^{\prime }{a}_{\sigma (k)}}^{{\gamma }_k^{\prime }{\gamma }_{\sigma (k)}}(t_k^{\prime },t_{\sigma (k)})\end{array}$$

$$\begin{array}{rl}& {\mathcal{V}}_0(t,t_k^{\prime })d_{a_k^{\prime }(t_k^{\prime }{\gamma }_k^{\prime })}^{†}\dots d_{a_1(t_1{\gamma }_1)}{\mathcal{V}}_0(t_1,0){\rho }_I(0)\\ \to & {\overline{\overline{\mathcal{V}}}}_0(t,t_k^{\prime }){\overline{\overline{d}}}_{a_k^{\prime }(t_k^{\prime }{\gamma }_k^{\prime })}^{†}\dots {\overline{\overline{d}}}_{a_1(t_1{\gamma }_1)}{\overline{\overline{\mathcal{V}}}}_0(t_1,0)|{\rho }_I(0)⟩\end{array}$$

$$\rho =\sum _{n,m}{\rho }_{nm}|n⟩⟨m|\to \sum _{n,m}{\rho }_{nm}|n⟩\otimes |m⟩\equiv |\rho ⟩$$

$${\overline{\overline{d}}}_{+}$$

$${\overline{\overline{d}}}_{-}$$

$$i{\int }_0^{t}d{t}_1{\int }_0^{t_1}d{t}_1^{\prime }{\mathcal{V}}_0(t,t_1^{\prime })d_{a_1^{\prime }-}^{†}{\mathcal{V}}_0(t_1^{\prime },t_1)d_{a_1+}{\mathcal{V}}_0(t_1,0){\Delta }_{a_1^{\prime },a_1}^{-+}(t_1^{\prime },t_1)$$

$$\mathcal{V}={\mathcal{V}}_0+{\mathcal{V}}_0\circ \Sigma \circ {\mathcal{V}}_0+{\mathcal{V}}_0\circ \Sigma \circ {\mathcal{V}}_0\circ \Sigma \circ {\mathcal{V}}_0+\dots$$

Full text

I Introduction

Small interacting quantum systems coupled to external environments represent basic paradigms of transport, dissipation and non equilibrium phenomena. Understanding the dynamical behavior of these open quantum systems is therefore crucial in many different physical contexts where the idealization of an isolated quantum system obeying perfectly unitary quantum dynamics is either to restrictive or unable to capture the fundamental physics.

In condensed matter physics the motivation comes from studying models of quantum dissipation and macroscopic quantum tunneling in the early days of Caldeira-Leggett and spin-boson models Caldeira and Leggett (1981); Leggett et al. (1987) which keep attracting a lot of interest Hur et al. (2018) or from diluted magnetic impurities in metals Hewson (1993) and transport through quantum dots and single molecules attached to leads Goldhaber-Gordon et al. (1998); Grobis et al. (2008); Roch et al. (2008); Iftikhar et al. (2018) leading to fermionic realizations of so called quantum impurity models. These consist of a small quantum systems with few interacting degrees of freedom, the impurity, coupled via hybridization to a gapless reservoir of fermionic or bosonic excitations. The dynamical correlations of such reservoirs, which decay in time as a power law at zero temperature and feature strong memory effects Cohen and Rabani (2011), together with a local many body interaction, make quantum impurity physics highly non-trivial. Nevertheless, methods to solve the dynamics of quantum impurity models have flourished in recent years, mainly driven by the developments of Diagrammatic Monte Carlo Gull et al. (2011); Mühlbacher and Rabani (2008); Schiró and Fabrizio (2009); Werner, Oka, and Millis (2009); Cohen et al. (2013); Profumo et al. (2015); Cohen et al. (2015); Chen, Cohen, and Reichman (2017a, b).

On a different front, recent advances in quantum optics, quantum electronics and quantum information science have brought forth novel classes of driven open quantum systems in which excitations are characterized by finite lifetime due to unavoidable losses, dephasing and decoherence processes originating from their coupling to an external electromagnetic environment. Examples include atomic and optical systems such as ultracold gases in optical lattices Bloch, Dalibard, and Nascimbène (2012) or trapped ions Blatt and Roos (2012), as well as solid state systems such as arrays of nonlinear superconducting microwave cavities Houck, Tureci, and Koch (2012); Hur et al. (2016). In these settings, the dissipative processes associated to the external environment can be very well described in terms of a Lindblad master equation for the system density matrix Breuer and Petruccione (2002). A major effort here is to conceive situations in which coupling to a quantum environment can involve non-linear combinations of system operators thus mediating effective interactions which act as a resource for quantum state preparation and to engineer a desired steady state. Such a dissipation engineering is actively investigated in quantum optics Diehl et al. (2008); Verstraete, Wolf, and Ignacio Cirac (2009); Murch et al. (2012). This has stimulated a new wave of interest around open Markovian quantum systems at the interface between quantum optics and condensed matter physics.

The examples above represent two well studied, yet substantially separated, paradigms of open quantum systems. At the same time much less is known about the interface between those two settings, namely the interplay between Markovian dissipation and the coupling to a fully structured, frequency dependent, Non-Markovian quantum bath, especially for what concerns the emergent many body physics. Interest around this new kind of quantum impurity problems has recently grown Nakagawa, Kawakami, and Ueda (2018); Tonielli et al. (2019); Fröml et al. (2019). Such a question is potentially relevant in a number of contexts. From one side, mesoscopic quantum devices have been successfully coupled to electromagnetic resonators hosting dissipative photon fields  Delbecq et al. (2011); Schiró and Le Hur (2014); Bruhat et al. (2016); Cottet et al. (2017) offering the possibility to investigate the fate exotic many body phases such as the Kondo effect in presence of Markovian dissipation.

On the other hand, in the context of quantum optics and quantum information the role of Non-Markovian bath correlations has been recently attracting enormous interest Liu et al. (2011); Ángel Rivas, Huelga, and Plenio (2014); Breuer et al. (2016) and there is urgent need to develop novel theoretical approaches to address this question.

With these motivations, in this work we focus on a model for a quantum impurity coupled to two kind of external environments, as depicted in figure 1: a quantum bath described by a set of non interacting modes bilinearly coupled to our impurity degree of freedom, as in conventional quantum impurity models, resulting in a frequency dependent Non-Markovian evolution for the impurity and a second environment, where the non interacting bath modes can be also non-linearly coupled to the impurity, but for which the reduced impurity dynamics would be well described by a Markovian master equation. The resulting impurity problem encodes therefore the interplay between those two kind of dissipative mechanisms and the local interaction on the impurity.

Taking inspiration from recent developments in quantum impurity physics Gull et al. (2011); Mühlbacher and Rabani (2008); Schiró and Fabrizio (2009); Werner, Oka, and Millis (2009) we develop here an hybridization expansion for the real-time evolution operator of the impurity, obtained after tracing out both the non-Markovian environment, as usually done in the literature, as well as the Markovian bath. The final result naturally generalizes the well known real-time hybridization expansions to this mixed Markovian/Non-Markovian context. In addition to its own interest and potential for the development of diagrammatic Monte Carlo sampling, this expansion allows us to formulate a self-consistent resummation technique for the real-time impurity evolution operator based on the Non-Crossing Approximation used in the context of quantum impurity models Bickers (1987); Nordlander et al. (1999); Eckstein and Werner (2010); Rüegg et al. (2013); Strand, Eckstein, and Werner (2015); Peronaci, Schiró, and Parcollet (2018). We derive and discuss in details this approach and test it on a simple fermionic model coupled to a zero temperature bath and in presence of Markovian dissipation.

The paper is organized as follows. In Section II we introduce the model and the general formulation of the problem. Section III is devoted to derive the hybridzation expansion for the real-time impurity propagator, after tracing out the non-Markovian bath (section III.A) and the Markovian environment (Section III.B). In Section IV we develop a self-consistent resummation based on the Non-Crossing Approximation, while in Section V we apply this method to a simple model.

II Model and General Formulation of the Problem

We consider a model of a quantum impurity, a small quantum system with a finite number of bosonic/fermionic degrees of freedom $[d_{a},d^{\dagger}_{b}]_{\pm}=\delta_{ab}$ and with Hamiltonian $H_{I}[d_{a},d^{\dagger}_{a}]$, coupled to two different quantum baths (see figure 1). We will denote the Hamiltonian of the baths with $H_{M}$ and $H_{\bar{M}}$, where the subscripts refer to the fact that ${\bar{M}}$ is a non-Markovian bath and $M$ is a Markovian one. We describe the two environments as a collection of non-interacting bosonic/fermionic modes, (respectively if the impurity is bosonic/fermionic)

[TABLE]

The total Hamiltonian therefore reads

[TABLE]

where we have introduced the two coupling terms between the impurity and the $M$ and $\bar{M}$ baths. We will consider the impurity to be bilinearly coupled to the $\bar{M}$ bath, i.e. through a coupling Hamiltonian of the form

[TABLE]

while the coupling between the impurity and the $M$ bath is taken of the most general form for which one can derive a Lindblad master-equation Breuer and Petruccione (2007):

[TABLE]

with $X_{\alpha}=X_{\alpha}^{\dagger}$, $B_{\alpha}=B_{\alpha}^{\dagger}$ generic operators respectively of the impurity and of the Markovian bath.

Defining the time evolution operator of the entire system as $U(t,0)=e^{-iHt}$ and given an initial condition for the system density matrix $\rho(0)$ we can formally write down the reduced density matrix of the impurity at time $t$, tracing out the degrees of freedom of the two environments

[TABLE]

from which the dynamics of simple impurity observables can be readily obtained as $O_{I}(t)=\text{tr}\left[\rho_{I}(t)O_{I}\right]$. With the assumption that the initial density operator of the environment and the impurity factorizes Breuer and Petruccione (2007), we can define the evolution operator of the reduced dynamics

[TABLE]

This reduced density operator and its evolution operator are the key quantities over which we will focus our attention throughout the manuscript. Performing the trace over the environment degrees of freedom is a highly non-trivial problem. In the following we will obtain two main results. The first one is a formal series from $\mathcal{V}$ to all orders in the coupling with the non-Markovian bath, called hybridization expansion, and the second one is a closed equation for $\mathcal{V}$, based on a self-consistent resummation of the series.

We stress that the order in which the trace over the two environments is taken in Eq. (4) is not crucial. While in this paper we will proceed by first taking the trace over the non-Markovian environment, resulting in an hybridization expansion, and then over the Markovian one, we could have equally reversed this choice and still obtain the same final result.

III Hybridization Expansion

In this section we derive a formal hybridization expansion (21) for the reduced density matrix of the impurity. Such an expansion is usually derived in the context of quantum impurity models coupled to a single bath, as a starting point to develop exact Monte-Carlo sampling Mühlbacher and Rabani (2008); Schiró and Fabrizio (2009) or approximated resummation techniques Aoki et al. (2014); Rüegg et al. (2013) to solve the problem. There the formulation is typically done at the level of the partition function, i.e. tracing out also the impurity degrees of freedom, while we are interested in the reduced density matrix and the evolution operator, see Eq. (5), therefore we will not perform such a trace, a fact that will result in some formal difference in the approach. More importantly, here the quantum impurity is also coupled to a second Markovian environment that we will need to trace out as well and this can be done exactly under the assumption that the $IM$ subsystem obeys a Markovian Lindblad master equation Breuer and Petruccione (2007); Carmichael (1999).

III.1 Tracing Over the Non-Markovian Bath

We begin by performing the trace over the non-Markovian environment which is quadratic in terms of bath operators and bilinearly coupled to the impurity. Such a trace could be performed exactly within a path integral formulation leading to an effective Keldysh action which is non-local in time. Here we proceed instead at the operator level by noticing that the trace could be taken exactly order by order in an expansion in the the coupling between the non-Markovian environment and the impurity.

In order to derive this expansion, we write down the full Hamiltonian of the system as $H=H_{0}+H_{\bar{M}}+H_{I\bar{M}}$, describing respectively the impurity embedded in the Markovian bath ($H_{0}=H_{I}+H_{M}+H_{IM}$), the Non-Markovian environment and its coupling to the impurity. We then move to the interaction picture with respect to the Hamiltonian $H_{0}+H_{\bar{M}}$. Introducing the standard time-ordering and anti-time-ordering operators $T_{t}$ and $\check{T}_{t}$, the density operator becomes

[TABLE]

We will perform a simultaneous expansion in powers of $H_{I{\bar{M}}}(t^{\prime})$ both on the left and on the right of the initial density operator $\rho(0)$. A formal way to manage a single series expansion and to write all the operators on the left side of the density operator, is to use the formalism of the Schwinger/Keldysh double contour $C(t,0)$ Schwinger (1961); Keldysh (1964) (fig. 2). Operators on the left (right) side of $\rho(0)$ are assigned a + or - label, so that the couple $(t,\gamma)\equiv t_{\gamma}$ with $\gamma\in\{+,-\}$, allows to locate one operator on this double time-axis.

We introduce the standard Keldysh time-ordering on the contour as follows:

[TABLE]

This ordering allows to define a time-ordering operator $T_{C}$, such that two operators, $X_{1}$ and $X_{2}$, being $X$ a creation or annihilation fermionic (bosonic) operator, anticommute (commute) under time-ordering: $T_{C}X_{1}{(t_{\gamma})}X_{2}{(t_{\gamma^{\prime}})}=\xi T_{C}X_{2}{(t_{\gamma^{\prime}})}X_{1}{(t_{\gamma})}$, with

[TABLE]

The time ordering operator is defined as follows

[TABLE]

which naturally extends to the case of $n$ operators. Once time ordered, the operators belonging to the - branch of the contour have to be brought on the right side of the density matrix, as if there we were exploiting the cyclic property of a trace.

By defining contour integrals as $\int_{C(0,t)}dt\,\equiv\int_{0}^{t_{+}}dt_{+}-\int_{0}^{t_{-}}dt_{-}$, one can show that the density operator evolution (6) can be written in the compact form

[TABLE]

and, accordingly, the evolution operator defined in (5) can be written as

[TABLE]

In order to perform the partial trace on the non-Markovian environment, we assume that at time $t=0$ there’s no entanglement between the non-Markovian bath and the rest of the system, such that the density operator factorizes $\rho(0)=\rho_{IM}(0)\otimes\rho_{\bar{M}}(0)$, with $\rho_{\bar{M}}(0)$ quadratic in bosonic/fermionic operators. Initial thermal states could be taken into account considering a third, imaginary time branch of the non-equilibrium contour Danielewicz (1984a, b); Wagner (1991), but this is beyond the interest of this work.

We then Taylor-expand the time-ordered exponential in power of the impurity-bath hybridization $H_{I\bar{M}}$ and perform the trace over the bath degrees of freedom, which immediately reduces the expansion to only even order terms. Then using Wick’s theorem and performing the sums over $\{b_{i},b_{i}^{\prime}\}$, we can write the final result in terms of the bath hybridization function

[TABLE]

where $G_{b^{\prime}}(t^{\prime}_{\gamma^{\prime}},t_{\gamma})=-i\langle T_{C}c_{b^{\prime}}(t^{\prime}_{\gamma^{\prime}})c_{b^{\prime}}^{\dagger}(t_{\gamma})\rangle$ is the contour ordered bath Green’s function. Finally, we obtain the hybridization expansion Aoki et al. (2014); Schiró and Fabrizio (2009); Schiró (2010):

[TABLE]

$\gamma_{i},\gamma^{\prime}_{i}$ are contour indices $\gamma\in\{+,-\}$. We notice that the hybridization function $\Delta^{\gamma_{i}^{\prime},\gamma_{\sigma(i)}}_{a_{i}^{\prime}a_{\sigma(i)}}(t^{\prime}_{i},t_{\sigma(i)})$ connects the $d_{a_{\sigma(i)}}({t_{\sigma(i)}},{\gamma_{\sigma(i)}})$ operator with the $d_{a_{i}^{\prime}}^{\dagger}({t^{\prime}_{i}},{\gamma^{\prime}_{i}})$ one. We can interpret this construction as follows. The $d$ operator creates a ”hole” in the impurity, which is propagated through the system and then annihilated by a $d$ operator. To this hole it corresponds (from the definition of $\Delta$) a particle of the environment which is created, propagated and annihilated. Thus, the series eventually describes processes in which particles hop from the impurity to the environment and back to the impurity.

III.2 Tracing over the Markovian bath

III.2.1 Super-operators formalism

It is useful to describe time-evolution using super-operators, as these are natural objects to describe the dynamics of open systems and since they provide a useful framework to work out the trace on the Markovian environment in eq. (11). We call super-operator an operator that acts on an operator, rather then on quantum state. The focus is shifted from the standard evolution operator $U(t,0)=e^{-iHt}$, which evolves a pure state (a ket) in time, to the super-operator $\mathcal{U}(t,0)$ which time-evolves a density operator and is defined by

[TABLE]

We can write a generic time-ordered string of operators, like it appears in eq. (11), in the Schrödinger’s picture and in a compact form, using the super-operators notation. This comes at the price of introducing some notation.

We promote $d$, $d^{\dagger}$ operators in the Schrödinger’s picture to super-operators $d_{\gamma}$, $d^{\dagger}_{\gamma}$, with an contour index $\gamma$ reminiscent of the branch the original operators belonged to:

[TABLE]

We trivially generalize the contour time-ordering operator $T_{C}$ to the super-operators notation

[TABLE]

The $X_{t,\gamma}$ super-operators are objects in Schrödinger’s picture and their time label $t$ is just meant to know how to order them.

We also need to introduce a further ”forward” time-ordering operator $T_{F}$, that orders two super-operators according to their time labels $t,t^{\prime}$, regardless of their contour index:

[TABLE]

This definition is the same for both fermions and bosons, with no extra minus signs for fermions.

Using these definitions, we can write the following identity

[TABLE]

The second line is a chain of subsequent time-evolutions operated by $\mathcal{U}_{0}$, going overall from time [math] to time $t$, alternated with the application of $d_{\gamma},d^{\dagger}_{\gamma}$ super-operators. We remark that the two time-order operators $T_{C}$ and $T_{F}$ do not commute. In order to evaluate the second line of eq. (16), one has first to order the super-operators according to $T_{C}$; this first ordering is necessary in order to compute the non-trivial sign factor obtained by swapping fermionic operators. Then the super-operators must be re-ordered according to the ”forward” time-ordering operator $T_{F}$. This insures that, in order to evaluate eq. (16), one has to apply only forward in time evolution super-operators.

III.2.2 Performing the partial trace

We now aim at performing the partial trace on the Markovian environment which is left in eq. (11). This trace is taken on a contour time-ordered string of impurity operators. The latter are nevertheless evolved by the joint dynamics of the impurity plus the remaining bath, making the partial trace non-trivial to evaluate. Assuming that the impurity-bath dynamics is governed by a Lindblad master equation, then the partial trace becomes trivial Gardiner and Collett (1985); Carmichael (1999). We report a proof here as this is a crucial step to obtain the hybridization expansion (21). We recall the Lindblad master equation Breuer and Petruccione (2007); Carmichael (1999):

[TABLE]

where $L_{\alpha}$ are the jump operators, microscopically determined by the environment-impurity coupling (3). $\rho_{I}^{0}(t)$ must not be confused with $\rho_{I}(t)=\text{tr}_{M}\rho_{IM}(t)$, as the former is the density operator obtained by evolving $\rho(0)$ in presence of the Markovian environment alone. $\rho_{I}(t)$ instead is obtained by evolving $\rho(0)$ with a dynamics that includes both the Markovian and non-Markovian environments. Defining the Markovian evolution super-operator, $\mathcal{V}_{0}(t-t^{\prime})=e^{\mathcal{L}\left(t-t^{\prime}\right)}$, with $t>t^{\prime}$, then

[TABLE]

We remark that $\mathcal{V}_{0}$ depends only on time differences as it satisfies (17). This is equivalent to

[TABLE]

In order to show how to perform the trace of the string of super-operators in (16), let’s assume time ordering is already enforced so that we don’t have to care about it. Defining $r_{1}(t)=\mathcal{U}_{0}(t,t_{k}^{\prime})r_{1}(t_{k}^{\prime})$ and $r_{1}(t_{k}^{\prime})=d^{\dagger}_{t_{k}^{\prime}\gamma_{k}^{\prime}}\mathcal{U}_{0}(t_{k}^{\prime},t_{k-1})\dots d_{t_{1}\gamma_{1}}\mathcal{U}_{0}(t_{1},0)\rho_{IM}(0)$, we can break down the tracing operation as follows:

[TABLE]

The last equality is analogous to eq. (19) and holds under the same assumptions leading to Lindblad master equation. One can iterate this procedure, as now $\text{tr}_{M}r_{1}(t_{k}^{\prime})=d^{\dagger}_{t_{k}^{\prime}\gamma_{k}^{\prime}}\text{tr}_{M}\left[\mathcal{U}_{0}(t_{k}^{\prime},t_{k-1})r_{2}(t_{k-1})\right]$, to turn all the $\mathcal{U}_{0}$ super-operators in eq. (11) in $\mathcal{V}_{0}$ ones.

III.2.3 Generalized hybridization expansion

We then get to the final form of the hybridization expansion in presence of both a non-Markovian and Markovian environment, that is one of the main results of this work:

[TABLE]

This series can be sampled using stochastic sampling techniques Mühlbacher and Rabani (2008); Schiró and Fabrizio (2009); Werner et al. (2006) or approximately resummed Aoki et al. (2014); Rüegg et al. (2013). For both purposes, it is useful to define the Feynman rules for the series (21).

III.2.4 Matrix representation

Each term of the hybridization expansion (21) must be understood as a composition of applications of super-operators, from the right-most to the left-most one, on the initial density operator. We remark that in the usual representation in which states of the Hilbert space are vectors and operators are matrices, super-operators are rank-4 tensors. Instead, if we write operators as vectors, then super-operators become matrices. This is convenient to evaluate terms of the hybridization expansion (21) as matrix products. Namely

[TABLE]

using double bars to indicate matrices. Operators are represented as vectors in a space which is the tensor product (indicated with $\otimes$) of two copies of the original Hilbert space. The vectorization procedure, taking the example of the density operator, reads:

[TABLE]

The matrices $\overline{\overline{{\mathcal{V}}}},\overline{\overline{{d}}},\overline{\overline{{d}}}^{\dagger}$ corresponding to the super-operators $\mathcal{V},d,d^{\dagger}$ are defined by the following simple procedure. Let’s consider the super-operator $\mathcal{S}\bullet=A\bullet B$. Representing $\rho$ as a vector $|\rho\rangle$, then also $\mathcal{S}\rho$ will be a represented as a vector according to eq. (23). It’s simple to show that $|\mathcal{S}\rho\rangle=\overline{\overline{{\mathcal{S}}}}|\rho\rangle$ defining $\overline{\overline{{\mathcal{S}}}}=\overline{\overline{{A}}}\otimes\overline{\overline{{B}}}^{T}$. The matrix form of $d_{\pm},d_{\pm}^{\dagger}$ in the doubled Hilbert space is

[TABLE]

with $\overline{\overline{{d}}}_{\pm}^{\dagger}$ the hermitian conjugate of $\overline{\overline{{d}}}_{\pm}$. These rules allow to obtain the matrix representation of the Liouvillian and then of $\mathcal{V}_{0}$.

III.3 Feynman Rules

the Feynman rules to draw the hybridization expansion (21) are represented in figure 3. We will use these rules to draw a term with $2k$ annihilation and creation super-operators, with a particular ordering for the times $\{t_{i},t_{i}^{\prime}\dots t_{1},t_{1}^{\prime}\}$ and a choice of a permutation $\{\sigma(1),\sigma(2)\dots\sigma(k)\}$. To do that, we draw a couple of parallel axes representing the double contour from time $t=0$ to time $t$. $d_{\gamma}$ ($d^{\dagger}_{\gamma}$) super-operators are represented as a dashed half-line with outwards (inwards) arrows, stemming from the contour branch $\gamma$. The dashed half-lines corresponding to the super-operators $d_{a_{i}^{\prime}(t_{i}^{\prime}\gamma_{i}^{\prime})}^{\dagger}$ and $d_{a_{\sigma(i)}(t_{\sigma(i)}\gamma_{\sigma(i)})}$ are joined together to form a hybridization line, representing the hybridization function $\Delta^{\gamma_{i}^{\prime}\gamma_{\sigma(i)}}_{a^{\prime}_{i}a_{\sigma}(i)}(t_{i}^{\prime},t_{\sigma}(i))$, which has an arrow going from $d$ to $d^{\dagger}$. Then, each part of the double contour between two integration times, drawn as two parallel solid segments, represents a time-propagation super-operator $\mathcal{V}_{0}$. The dressed evolution operator $\mathcal{V}$ is drawn by replacing the contour solid lines by double lines. As an example, the diagram corresponding to

[TABLE]

is shown in figure 4.

All the diagrams with $2k$ annihilation and creation super-operators are generated by connecting $d^{\dagger}$ super-operators to $d$ ones in all possible choices of permutations $\sigma$ and considering all possible time orderings of integration times.

IV Self-Consistent Diagrammatic Resummations Techniques

In this section we start from the hybridization expansion derived in section III, which involves bare diagrams to all orders, and use the Feynman rules to introduce diagrammatic resummation techniques.

To proceed further it is useful to draw more compact diagrams where the double contour is collapsed on a single time-axis and thus a time propagation $\mathcal{V}_{0}$ is represented by a single line, as we show in figure 5. These compact diagrams represent an ensemble of diagrams drawn with the rules we introduced in 3.

The advantage of this notation is that all the diagrams represented by a single compact diagram have the same topology in terms of being 1-particle irreducible or non-crossing. Figure 6 shows the hybridization expansion drawn using these compact diagrams.

IV.1 Dyson Equation

As a first step it is useful to distinguish diagrams which are one-particle irreducibles, i.e. compact diagrams which cannot be separated, by cutting a solid line, in two parts that are not connected by any hybridization line, as indicated in figure 6. Then, we introduce the self-energy $\Sigma$ as the sum of one particle irreducible (1PI) diagrams. All the non-1PI diagrams can be obtained by joining some 1PI diagrams with solid lines, thus the whole series can be written as

[TABLE]

We remark that the objects composing this series, $\mathcal{V}_{0}$ and $\Sigma$, are super-operators and the series must be understood as a composition of applications of super-operators, from the right-most to the left-most one, on a target operator. Self-energies and propagators are joined by the circle operation, $\circ$, standing for a super-operator application and a partial time convolution. Using brackets to stress that we refer to a super-operator application and the symbol $\bullet$ to indicate a target operator, we have

[TABLE]

The series above sums up to the Dyson equation $\mathcal{V}=\mathcal{V}_{0}+\mathcal{V}_{0}\circ\Sigma\circ\mathcal{V}=\mathcal{V}_{0}+\mathcal{V}\circ\Sigma\circ\mathcal{V}_{0}$, or equivalently, in integro-differential form

[TABLE]

When the self-energy of the non-Markovian environment $\Sigma$ is set to zero, this equation yields $\mathcal{V}(t)=e^{\mathcal{L}t}$, which is the Lindblad evolution.

One of the main effects of dissipative dynamics is that the system may forget about initial conditions and reach the same stationary state for any initial condition. Assuming a stationary state exists for a non-Markovian map $\mathcal{V}$ defined by the Dyson equation (25), then it satisfies

[TABLE]

Setting the non-Markovian self-energy to zero, this equation reduces to the Lindblad condition for the stationary state. The derivation of this equation is in the supplementary material.

IV.2 The Non-Crossing Approximation

The non-crossing approximation (NCA) corresponds to approximating the series for $\mathcal{V}$, and thus also for $\Sigma$, by considering only the compact diagrams in which the hybridization lines do not cross Bickers (1987); Nordlander et al. (1999); Eckstein and Werner (2010); Rüegg et al. (2013). The NCA diagrams composing the self-energy are shown in figure 7.

In order to prove the second equality in figure 7, we remark that the first and last times of a self-energy diagram must be connected together by an hybridization line. If it’s not the case, in fact, the resulting diagram is either non-1PI or it’s crossed. Then all the diagrams of $\Sigma$ (in NCA) are obtained connecting the intermediate times to form all the possible non-crossing diagrams, (not only the 1PI ones this time). But the latter diagrams in turn define the NCA series for $\mathcal{V}$. This proves the second equality in figure 7. We remark that then, the NCA self-energy coincides with its contributions (with $k=1$), where the bare propagator $\mathcal{V}_{0}$ is replaced with the dressed one $\mathcal{V}$.

To obtain an analytic expression of the self-energy, we have to cast the $k=1$ term of the hybridization expansion (21) in a form in which the innermost integration time is lower than the outermost, that is with integrals of the form $\int_{t^{\prime}}^{t}dt_{1}\int_{t^{\prime}}^{t_{1}}dt_{2}$. In doing so, one must deal with the signs coming out of the time ordering. We report the calculation in the supplementary material.

The expression for the self-energy eventually reads

[TABLE]

where $\alpha,\beta\in\{+,-\}$ are contour indices, $a,b$ are the fermionic generic indices.

We can interpret the two terms in eq. (27) as follows. The first term propagates a hole in the impurity (applies $d^{\dagger}$ first and then $d$) and a particle in the bath, the latter being described by a hybridization function with the same time arguments of $\mathcal{V}$; The second term propagates a particle in the impurity and a hole in the bath with a hybridization function with opposite time arguments than $\mathcal{V}$.

Few comments are in order here, concerning the above result. First, in absence of the Markovian environment, that is by replacing $\mathcal{V}(t,t^{\prime})$ and $\mathcal{V}_{0}(t,t^{\prime})$ with $\mathcal{U}(t,t^{\prime})$ and $\mathcal{U}_{0}(t,t^{\prime})$, our results are equivalent to non-equilibium NCA schemes for unitary dynamics Eckstein and Werner (2010); Strand, Eckstein, and Werner (2015); Peronaci, Schiró, and Parcollet (2018). There is a formal difference consisting in our super-operators formulation of the hybridization expansion and of the Dyson equation, that is necessary to consider the additional Markovian environment without further approximations, but this difference is only formal and does not affect the results. In fact, if $N$ is the dimensionality of the impurity Hilbert space, the usual non-equilibrium NCA propagator Strand, Eckstein, and Werner (2015); Peronaci, Schiró, and Parcollet (2018) has different Keldysh components, each of them being an $N\times N$ matrix, while our $\mathcal{V}$ propagator is a $N^{2}\times N^{2}$ matrix with no Keldysh components.

Furthermore, the result obtained for the NCA self-energy in Eq. (27) makes clear that for a bath hybridization which is delta-correlated in time the resulting self-energy contribution to the Dyson equation takes the form of an additional Lindblad dissipator. In other words, one can recover the Lindblad master equation from our diagrammatic NCA approach in the Markovian limit, and possibly discuss corrections to the master equation from higher order terms, as recently done Müller and Stace (2017).

IV.3 Properties of the NCA propagator

The propagator $\mathcal{V}(t,t^{\prime})$ obtained in NCA is the time-evolution super-operator of the reduced density operator of the impurity. Assuming to switch on the interaction with the baths at time $t=0$, then the density operator of the impurity at time $t$ is given by $\rho_{I}(t)=\mathcal{V}(t,0)\rho_{I}(0)$. A time evolution super-operator must be a convex-linear, completely positive and trace-preserving map Breuer and Petruccione (2007). It is natural to ask which of these properties are preserved by the NCA approximation. $\mathcal{V}$ is a obviously a linear map, implying it is also convex linear. We proved that it is trace-preserving and that it preserves hermiticity (see supplementary material), while proving or disproving whether the map is completely positive is a tough task Reimer and Wegewijs (2018) that will be addressed in the future. We stress that $\mathcal{V}(t,t^{\prime})$ describes a non-Markovian evolution, so it does not form a semi-group, that is $\mathcal{V}(t,t^{\prime})\neq\mathcal{V}(t,t_{1})\mathcal{V}(t_{1},t^{\prime})$ with $t^{\prime}<t_{1}<t$. Time-evolution super-operators have also interesting spectral properties following from trace preservation. We refer to the supplementary material for the proof. We call $\lambda_{i}(t,t^{\prime}),v^{R}_{i}(t,t^{\prime})$ the eigenvalues and right eigenvectors of $\mathcal{V}(t,t^{\prime})$, depending on time. As it preserves the trace, $\mathcal{V}(t,t^{\prime})$ must have at least one eigenvalue equal to one, say $\lambda_{0}\equiv 1$. If we assume this eigenvalue is non-degenerate, then all the others eigenvectors with $i\neq 0$, are traceless. As a consequence of these properties, if one evolves an initial state $\rho(0)$ and expands $\rho(t)$ on the instantaneous eigenvectors $v^{R}_{i}(t,0)$, then those eigenvectors with $i\neq 0$ will represent decay modes of the dynamics as they will be suppressed by their corresponding vanishing eigenvalues for long times, while $v^{R}_{i=0}(t,0)$ will evolve in time undumped until reaching a stationary value, representing the stationary state of the non-Markovian evolution. We will numerically check these properties in figure 8 of the next section, where we will apply our NCA algorithm to a specific example.

V Case study: spin-less fermionic impurity

As a non-trivial application of the NCA approach for open system described so far, we consider here a model of a single-mode, spin-less fermionic impurity with Markovian losses, pump and dephasing and further coupled to a non-Markovian fermionic environment, as described in equation (1). We notice that the model in absence of dephasing, also known as Resonant Level Model, is quadratic in all the fermionic degrees of freedom and therefore easily solvable, with analytical expressions known for the wide band limit. At finite dephasing this is no longer the case and the model cannot, to the best of our knowledge, be solved by simple means. This could be understood naturally in the Keldysh approach, where the dephasing would result in density-density type of coupling between different Keldysh branches.

The Markovian dynamics is described by a Lindblad master equation

[TABLE]

where $\epsilon_{0}$ is the energy of the fermionic level.

The effect of the non-Markovian environment on the impurity is completely determined by its hybridization function (10). Here we choose a zero temperature, particle hole symmetric, fermionic bath with constant density of states of bandwidth $2w$ and with coupling strength to the impurity $\eta$. In this case the hybridization function depends only on time-differences. As a consequence, one can show that also $\mathcal{V}$ and $\Sigma$ depend only on time differences and we will set $t^{\prime}=0$. With these definitions we get for the hybridization functions

[TABLE]

To solve the Dyson equation (25) numerically we use the simple discretization scheme $\partial f(t)=\left[f(t+\Delta t)-f(t)\right]/\Delta t$, $\int_{0}^{t}dt_{1}f(t_{1})=\Delta t/2\sum_{l=0}^{t/\Delta t-1}\left[f((l+1)\Delta t)+f(l\Delta t)\right]$, with time-step $\Delta t$. More refined integration methods are explained in detail in Aoki et al. (2014). The Hilbert space of the impurity has size $N=2$, so that the super-operator $\mathcal{V}$ has size $N^{2}\times N^{2}=4\times 4$ .

V.1 Results

We start analyzing the spectral properties of the propagator $\mathcal{V}(t)$, which have been discussed generically in the previous section, and are reported in figure 8. In the main panel we plot the time dependence of the absolute value of the eigenvalues of $\mathcal{V}(t)$, both in the purely Markovian case (dashed lines) as well in presence of both kind of dissipations. In both cases there is an eigenvalue which remains equal to one, while the others decay to zero at long times, as pointed out in section IV.3. However the nature of this decay is rather different in the two cases, showing a faster dynamics and long time oscillations in the non-Markovian case as opposed to a pure exponential decay in the Markovian one. The inset of figure 8 shows instead that all the right eigenstates of $\mathcal{V}(t)$ with different-from-one eigenvalues are traceless, while the right eigenstate with eigenvalue one has a finite trace (that we could normalize to one at every time).

We then consider the dynamics of a simple observable, such as the density of fermions in the impurity level, as a function of time and for different parameters (see panel figure 9). In the left figures we plot the dynamics for different values of the coupling $\eta$ (top) and bandwidth $w$ (bottom) of the non-Markovian environment, in presence of fixed Markovian losses, pump and dephasing. We see that with respect to the purely Markovian dynamics, characterized by a simple exponential relaxation, the NCA approach captures aspects related to the non-Markovian nature of the environment. In particular the dynamics becomes characterized by oscillations whose amplitude and frequency increase with the coupling $\eta$. Similarly, increasing the bandwidth of the non-Markovian environment reduces the oscillations in the population dynamics, which disappear in the large bandwidth limit, as it is the case for the unitary dynamics of the Resonant Level Model. This is not surprising, since oscillations at short-times $t\sim 1/w$ come from high-energy modes of the bath.

Overall the Non-Markovian environment makes the dynamics substantially faster. In the top-right plot we discuss the role of dephasing, that is actually very interesting as it shows an effect of the combined Markovian and non-Markovian environments. The dashed line shows that, for Markovian dissipation only, the dephasing does not affect the population dynamics; this is well understood as the dephasing dissipator commutes with the number operator. It’s interesting to see that, instead, combined with a non-Markovian environment the dephasing has an impact on the dynamics of populations; this is a smaller effect as it involves both the Markovian and the non-Markovian environments. This effect can be understood as follows: let’s consider a process in the time-evolution in which the non-Markovian environment applies a $d^{\dagger}_{+}$ operator on the density matrix and then a $d_{+}$ operator after some time. The application of the creation operator converts the populations of the density matrix into coherences. Then the Markovian dephasing dumps those coherences, which are then converted back to populations when the non-Markovian environment applies the annihilation operator. As a net effect, the dephasing has produced a change into populations, as it is shown in the top-right plot. We also note that not only the dynamics, but also the stationary values of the occupation change with the dephasing. A more direct effect of the dephasing appears in the coherences, i.e. in the off-diagonal elements of the density matrix, which decay to zero faster as $\gamma_{d}$ is increased, as we show in figure 10.

For what concerns the stationary state, we notice from the bottom-right plot that the average density would be independent of the energy of the fermionic level in the purely Markovian case, which leads to a infinite-temperature fully-mixed stationary state for the chosen dissipation rates. This makes sense as a Markovian bath has no energy structure, thus the level effectively sees always the same bath even if it’s shifted in energy. On the other hand the coupling to the non-Markovian bath makes the population depend strongly on the position of the energy level and gives a result which is in good agreement with exact analytical calculations (dashed line); to justify the quantitative discrepancy with this analytical result, we stress that for the non-interacting model we consider here the NCA approximation, which is based on a strong coupling expansion, is not expected to be exact.

VI Conclusions

In this work we have focused on a model for a quantum impurity coupled simultaneously to a Markovian and a non-Markovian environment. We derived a formal hybridization expansion for the evolution super-operator of the impurity, obtained after tracing out all the bath degrees of freedom. This result generalizes to non-unitary, Markovian case the hybridization expansion obtained for unitary quantum impurity models. As such it provides the natural starting point for the development of stochastic sampling techniques of the dissipative real-time dynamics of the impurity based on Diagrammatic Monte Carlo, that we leave for future studies.

Starting from this expansion we define real-time diagrammatic rules and write down a Dyson Equation for the impurity propagator that we evaluate retaining only non-crossing diagrams, an approximation which is known to capture some aspects of the impurity physics at strong coupling. The resulting approach leads to a trace and hermiticity preserving Non-Markovian dynamical map, with consequences on the spectral properties of the evolution super-operator, while proving its complete positivity in full generality remains an open question.

It is interesting to comment on the relation between our approach and related methods to deal with impurity models coupled to multiple baths. While in principle both the hybridization expansion Golež, Eckstein, and Werner (2015) as well as the strong coupling diagrammatic resummation Chen et al. (2016) can be generalized in presence of multiple environments, taking the Markovian limit from the start has some practical and conceptual advantage. In particular, we can take direct advantage of the local nature of the Markovian evolution and perform an expansion around an atomic-limit which now contains not only interaction but also drive and dissipation. This limit can be solved exactly by direct diagonalization of a Lindbladian, as opposed to treating exactly the degrees of freedom of the environment which are not expected to introduce new physics in presence of memory-less Markovian correlations. The key idea is therefore to treat on equal footing all the energy scales related to fast processes, while resorting to perturbation theory when dealing with processes leading to slowing decay correlations such as the coupling to gapless reservoirs.

As an application, we solved numerically the Dyson equation for the simple model of a fermionic, single-mode impurity, with Markovian losses, pump and dephasing and a non-Markovian, zero temperature environment. This model is non-trivial for the presence of dephasing, which is a quartic term in fermionic operators. This simple implementation allowed to check the spectral properties of the evolution super-operator and to study how Markovian dynamics gets modified by coupling to a non-Markovian environment. In particular our method allowed to show a physical consequence of coupling simultaneously to Markovian and non-Markovian environments: Markovian dephasing combined with non-Markovian processes leads to a change in impurity occupations. Future directions include the exploration of more complex impurity models involving internal degrees of freedom such as the Anderson Impurity model as well as bosonic extensions and to use of this method as an impurity solver within a dynamical mean field theory approach to driven-dissipative systems.

VII Supplementary Material

The Supplementary Material includes (i) a derivation of the equation for the steady state density matrix starting from the Dyson Equation for the non-Markovian map (ii) the detailed derivation of the NCA self-energy and (iii) the proofs of trace and hermiticity preservation of the NCA dynamics.

VIII Stationary state

Assuming a stationary state exists for a non-Markovian map $\mathcal{V}$ defined by the Dyson equation

[TABLE]

then, setting $t^{\prime}=0$, it satisfies

[TABLE]

In order to derive this equation, we focus on the propagator $V(\infty,0)$, that projects any initial state on the stationary state: $\rho_{ss}=V(\infty,0)\rho(0)$. This propagator obeys the Dyson equation

[TABLE]

At this point we need to make some assumptions based on physical arguments. $\mathcal{V}(t_{1},0)$ is expected to have a transient dynamics in a finite time interval of duration $t_{\text{tr}}$, starting at time $t_{1}=0$, and then to become stationary, i.e. $\lim_{t\rightarrow\infty}\partial_{t}\mathcal{V}(t,0)=0$. In addition, the system is supposed to lose memory of initials conditions, thus the convolution in the above Dyson equation must be cutoffed by the self-energy $\Sigma(\infty,t_{1})$ for $\infty-t_{1}>t_{\text{mem}}$. Then, in the region where $\Sigma(\infty,t_{1})$ is non-zero, $\mathcal{V}(\infty,t_{1})$ is stationary and we can replace it with $\mathcal{V}(\infty,0)$. With these arguments we get $\left(\mathcal{L}+\int_{0}^{\infty}dt_{1}\Sigma(\infty,t_{1})\right)\mathcal{V}(\infty,0)=0$ and applying it to any initial state we find (S2).

IX Derivation of NCA self-energy

To obtain an analytic expression of the self-energy, we have to cast the $k=1$ term of the hybridization expansion obtained in the main text in a form in which the innermost integration time is lower than the outermost, that is with integrals of the form $\int_{t^{\prime}}^{t}dt_{1}\int_{t^{\prime}}^{t_{1}}dt_{2}$. In doing so, sign of each contribution will be determined by the time-ordering $T_{C}$. The $k=1$ term of the hybridization expansion, omitting the fermionic indices to simplify the expressions, reads

[TABLE]

where we have summed on contour indices $\gamma_{1}$ and $\gamma_{2}$. Now we brake the integrals in two pieces, for $t_{2}>t_{1}$ and $t_{2}<t_{1}$, which allows to put the operators in a time-ordered fashion according to $T_{C}$.

[TABLE]

where we omitted the time arguments of the $\mathcal{V}_{0}$ operators. The $T_{F}$ time ordering will now sort the operators in ascending order in time from right to left:

[TABLE]

Using the fact that $\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}=\int_{0}^{t}dt_{2}\int_{0}^{t_{2}}dt_{1}$ and exchanging the integration times in the second column we get

[TABLE]

[TABLE]

which defines the NCA self-energy

[TABLE]

X Trace preservation

The evolution super-operator $\mathcal{V}$ obtained in the non-crossing approximation preserves the trace of the density operator, that is

[TABLE]

Writing $\mathcal{V}(t,t^{\prime})\bullet$ means applying the super-operator $\mathcal{V}$ on a generic operator $\bullet$; we recall that in the usual representation where the operator $\bullet$ is a matrix, this is not a matrix product. To prove this, we take the trace of the NCA Dyson equation and we use the fact that $\mathcal{V}_{0}$ does preserve the trace

[TABLE]

Then, using the expression for the self-energy eq. (S7), we can prove that the integrand vanishes, which completes the proof. In order to show that, we remark that for the cyclic property of the trace it holds that $\text{tr}\left[X_{+}\bullet\right]=\text{tr}\left[X_{-}\bullet\right]$, $X$ being a generic super-operator. Then we can fix the $d^{\dagger}_{\beta},d_{\beta}$ super-operators to be $d^{\dagger}_{+},d_{+}$ under the trace, getting

[TABLE]

By summing over $\beta$, one gets the two terms $\Delta_{ba}^{+\alpha}(t_{1},t_{2})-\Delta_{ba}^{-\alpha}(t_{1},t_{2})$ and $\Delta_{ab}^{\alpha+}(t_{2},t_{1})-\Delta_{ab}^{\alpha-}(t_{2},t_{1})$, which vanish because of the following identities, holding for $t_{1}>t_{2}$:

[TABLE]

These identities hold because of the definition of $\Delta$, given in the main text, in terms of contour time-ordered Green’s functions, i.e. $\Delta^{\alpha,\beta}({t_{1}},{t_{2}})\sim-i\langle{T_{C}c(t_{1},\alpha)c^{\dagger}(t_{2},\beta)}\rangle$, and remembering that times on the $-$ contour branch come after times on the $+$ one.

X.1 Spectral properties the propagator $\mathcal{V}$

We call $\lambda_{i}(t,t^{\prime}),v^{R}_{i}(t,t^{\prime}),v^{L}_{i}(t,t^{\prime})$ the eigenvalues and right and left eigenvectors of $\mathcal{V}(t,t^{\prime})$, which depend on time. As it preserves the trace, $\mathcal{V}(t,t^{\prime})$ must have at least one eigenvalue equal to one, say $\lambda_{0}\equiv 1$. If we assume this eigenvalue is non-degenerate, then all the others eigenvectors with $i\neq 0$, are traceless. The proof of these properties goes as follows. $\text{tr}{\left[\mathcal{V}(t,t^{\prime})\rho\right]}=\text{tr}\left[\rho\right]$ in the matrix notation reads $\langle\mathbb{1}|\bar{\bar{\mathcal{V}}}(t,t^{\prime})|\rho\rangle=\langle{\mathbb{1}|\rho}\rangle$ which holds for every $|\rho\rangle$ as the trace is preserved; then $\langle\mathbb{1}|$ must be a left eigenvector of $\mathcal{V}(t,t^{\prime})$, $\langle v^{L}_{0}|\equiv\langle\mathbb{1}|$, with eigenvalue $\lambda_{0}=1$. If we assume that there is only one eigenvector with eigenvalue 1, then all the others right-eigenvectors of $\mathcal{V}(t,t^{\prime})$ must be orthogonal to $\langle\mathbb{1}|$, that is they must have zero trace: $\langle{\mathbb{1}|v_{i}^{R}(t,t^{\prime})}\rangle=\text{tr}{\left[v_{i}^{R}(t,t^{\prime})\right]}=0$, for $i\neq 0$.

XI Hermiticity preservation

A quantum dynamical map evolving the density operator should preserve its Hermiticity: will show that this property is not spoiled by the NCA approximation. The proof is inductive and goes as showing that, if $\mathcal{V}(t,t^{\prime})$ preserves Hermiticity, then $\mathcal{V}(t+dt,t^{\prime})$ does; given the initial condition $\mathcal{V}(t^{\prime},t^{\prime})=\mathbb{1}$, that obviously preserves Hermiticity, then it follows that $\mathcal{V}(t,t^{\prime})$ is Hermiticity preserving $\forall t$. Assuming $\mathcal{V}(t,t^{\prime})$ is analytic in $t$, then its increment is given by its Taylor series

[TABLE]

From the Dyson equation (S1), we can show that if $\mathcal{V}(t,t^{\prime})$ preserves Hermiticity, then all its derivatives do, which in turn implies, from the Taylor expansion, that $\mathcal{V}(t+dt,t^{\prime})$ does. This ultimately comes for the causal structure of the Dyson equation.

We will restrict to show that $\partial_{t}\mathcal{V}(t,t^{\prime})$ is Hermiticity preserving, assuming $\mathcal{V}(t,t^{\prime})$ is. This result can be generalized to higher order derivatives, obtained by taking derivatives of the Dyson equation (S1), with two observations: the n-th derivative of $\mathcal{V}$ depends only on its lower order derivatives; the structure of the equation for the n-th derivative is such that, if the lower order derivatives are Hermiticity preserving, then also the n-th derivative is.

We now assume $\mathcal{V}(t,t^{\prime})$ Hermiticity preserving and show that this implies $\partial_{t}\mathcal{V}(t,t^{\prime})$ also is. With $\bullet$ an Hermitian operator, Hermiticity preservation of $\mathcal{V}$ means $\left(\mathcal{V}\bullet\right)^{\dagger}=\mathcal{V}\bullet$. $\mathcal{L}$ preserves Hermiticity as it is a Lindblad generator. Then, taking the hermitian conjugate of the Dyson equation (S1) we get

[TABLE]

We need to determine the hermitian conjugate of $\Sigma(t_{1},t_{2})\,\bullet$. $\Sigma(t_{1},t_{2})$ in the NCA approximation is given by (S7) and it depends on the hybridization function $\Delta$. Defining the Keldysh indices $\bar{\alpha}=-\alpha$, $\bar{\beta}=-\beta$ and with ∗ meaning complex conjugation, the following property holds

[TABLE]

that can be proven from the definition of $\Delta^{\alpha,\beta}({t_{1}},{t_{2}})\propto-i\langle{T_{C}c(t_{1},\alpha)c^{\dagger}(t_{2},\beta)}\rangle$ and writing down explicitly its Keldysh components. It also holds that

[TABLE]

as $\left(X_{+}\bullet\right)^{\dagger}=\left(X\bullet\right)^{\dagger}=\bullet^{\dagger}X^{\dagger}=X_{-}^{\dagger}\bullet^{\dagger}$ and $\left(X_{-}\bullet\right)^{\dagger}=\left(\bullet X\right)^{\dagger}=X^{\dagger}\bullet^{\dagger}=X_{+}^{\dagger}\bullet^{\dagger}$. For a nested application of super-operators as it appears in the self-energy, this property gives $\left(X_{1\alpha}\mathcal{V}X_{2\beta}\bullet\right)^{\dagger}=X^{\dagger}_{1\bar{\alpha}}\mathcal{V}X^{\dagger}_{2\bar{\beta}}\,\bullet$, where we have used that the our ansaz for $\mathcal{V}$ preserves hermiticity and that $\bullet$ is hermitian. Using these two results it follows that $\left(\Sigma(t_{1},t_{2})\bullet\right)^{\dagger}=\Sigma(t_{1},t_{2})\bullet$ as

[TABLE]

In the one but last equality, $\xi\left(-1\right)^{{(1+\xi)/2}}=-1$, for both bosons and fermions $\xi=\pm 1$. This completes the proof as

[TABLE]