Extracting dynamical maps of non-Markovian open quantum systems

Strachan, David J.; Purkayastha, Archak; Clark, Stephen R.

doi:10.1063/5.0228428

The most general description of quantum evolution up to a time τ is a completely positive tracing preserving map known as a dynamical map $\hat{Λ} (τ)$ ⁠. Here, we consider $\hat{Λ} (τ)$ arising from suddenly coupling a system to one or more thermal baths with a strength that is neither weak nor strong. Given no clear separation of characteristic system/bath time scales, $\hat{Λ} (τ)$ is generically expected to be non-Markovian; however, we do assume the ensuing dynamics has a unique steady state, implying the baths possess a finite memory time τ_m. By combining several techniques within a tensor network framework, we directly and accurately extract $\hat{Λ} (τ)$ for a small number of interacting fermionic modes coupled to infinite non-interacting Fermi baths. First, we use an orthogonal polynomial mapping and thermofield doubling to arrive at a purified chain representation of the baths whose length directly equates to a time over which the dynamics of the infinite baths is faithfully captured. Second, we employ the Choi–Jamiolkowski isomorphism so that $\hat{Λ} (τ)$ can be fully reconstructed from a single pure state calculation of the unitary dynamics of the system, bath and their replica auxiliary modes up to time τ. From $\hat{Λ} (τ)$ ⁠, we also compute the time local propagator $\hat{L} (τ)$ ⁠. By examining the convergence with τ of the instantaneous fixed points of these objects, we establish their respective memory times $τ_{m}^{Λ}$ and $τ_{m}^{L}$ ⁠. Beyond these times, the propagator $\hat{L} (τ)$ and dynamical map $\hat{Λ} (τ)$ accurately describe all the subsequent long-time relaxation dynamics up to stationarity. These timescales form a hierarchy $τ_{m}^{L} \leq τ_{m}^{Λ} \leq τ_{re}$ ⁠, where τ_re is a characteristic relaxation time of the dynamics. Our numerical examples of interacting spinless Fermi chains and the single impurity Anderson model demonstrate regimes where τ_re ≫ τ_m, where our approach can offer a significant speedup in determining the stationary state compared to directly simulating the long-time limit. Our results also show that having access to $\hat{Λ} (τ)$ affords a number of insightful analyses of the open system thus far not commonly exploited.

I. INTRODUCTION

The problem of theoretically obtaining transient relaxation dynamics and steady states of open many-body quantum systems, both in and out of equilibrium, is a notoriously challenging problem. Yet, a better understanding of dissipation and decoherence is crucial for understanding real-world quantum systems. The recent developments in nanoscale devices have motivated interest in technologies that harness quantum effects to greatly outperform their classical analogs,^1–4 while the existence of non-trivial quantum effects in biological systems is being increasingly investigated.^5–7 Existing approaches to modeling open quantum systems often use perturbative arguments, requiring a clear separation of energy or time scales. A popular approach is the master equation formalism, based on the assumptions of weak system–bath coupling and the Markov approximation. However, their application is limited and can lead to drastically incorrect predictions for composite systems, where local descriptions of dissipation violate thermodynamic laws^8,9 while global descriptions¹⁰ fail to capture the coherences necessary to describe non-equilibrium steady states, even in the limit of weak coupling.¹¹

Alternatively, non-equilibrium Green functions¹² can be used to model energy transport under strong system–reservoir coupling, but interactions within the systems must be treated perturbatively.^12–14 To proceed with less severe approximations, one must consider the environment more explicitly, a task that appears formidable due to the many, possibly infinite environmental degrees of freedom. However, the true state of a combined system–bath setup is often characterized by an amount of information much smaller than the maximum capacity of the corresponding Hilbert space. This is exploited in the time-evolving matrix product operator (TEMPO) formalism where the effect of the environment on the system is described by Feynman–Vernon influence functionals as matrix-product states in the temporal domain^15–18 or more generally a process tensor.^19–22 Alternatively, one can use the time derivative hierarchical equations of motion (HEOM) approach.²³ In this work, we follow the approach of Kohn and Santoro²⁴ and describe both the system and environment as a single matrix product state, which is then time-evolved directly. This is motivated by the fact that many open-system Hamiltonians can be mapped onto 1D structures containing only nearest-neighbor interactions.²⁵ As such 1D tensor networks such as matrix product states (MPSs) can often be used to accurately reproduce the dynamics of the system-environment setup by compressing the time-evolved state up to a finite bond-dimension χ, defining the dominant numerical error. A crucial result of these 1D chain mappings is the fact that even in the absence of compression, the size of the bath required scales almost linearly with the evolution time required.^26–28 Therefore, the finite time evolution of a system coupled to an infinite bath can be modeled exactly with a finite bath.

This observation is exploited in the periodically refreshed baths (PReB) formalism,²⁹ which allows full reconstruction of the open-system dynamics up to time τ, provided there exists a unique steady state. It works by recursively evolving up to a much smaller time $τ_{m}^{Λ}$ associated with the memory of the bath correlations, before tracing out the bath and replacing it with one in its initial product state. This bares much similarity to a collisional model of an open system,^30–33 where the bath is represented by a reservoir of identically prepared ancilla each of which interacts with the system once for a short time before being traced out. Typically the ancilla are single sites, whereas in the PReB approach, finite-sized chains are used, which allows for significant non-Markovian effects and more accurate modeling of generic bath spectral functions. Performing a direct single evolution of the system + bath evolution using tensor networks often results in an intractably large bond dimension in the long-time limit owing to the generic dynamical growth of entanglement. This makes the stationary properties hard to evaluate. The main advantage of PReB is that the entanglement growth within the baths is removed with each refreshment, limiting the growth of χ.

Motivated by PReB, in this paper, we investigate the combination of thermofield doubling³⁴ and the Choi–Jamiolkowski isomorphism^35,36 to extract the effect of one PReB cycle on a fermionic system as a dynamical map. Assuming the PReB cycle has evolved well beyond $τ_{m}^{Λ}$ ⁠, the steady state and all transient dynamics from any initial state can then be extracted from a single simulation of an enlarged system described by a pure state. In addition, we calculate the time local propagator $\frac{\partial \hat{ρ}}{\partial t} = \hat{L} (t) [\hat{ρ} (t)]$ and also extract the steady state when it first becomes invariant. This is directly related to the memory kernel in the Nakajima–Zwanzig equation, $\frac{\partial \hat{ρ}}{\partial t} = - i \int_{t_{0}}^{t} d s \hat{K} (t, s) \hat{ρ} (s)$ ⁠, via a fixed point relation.³⁷ A key point here is that $\hat{K} (t, s)$ is often highly nontrivial and intractable to calculate.^38–40 There have been two notable branches of work on calculating memory kernels, the first of which involves numerically solving a Volterra equation of the second kind^41–44 and another is based on the transfer matrix approach.^45–47 Here, we provide a method to numerically extract $\hat{L} (t)$ ⁠, which gives an equivalent description to $\hat{K} (t, s)$ ⁠, and $\hat{Λ} (t)$ ⁠, allowing us direct access to these rich objects.

The structure of this paper is as follows: in Sec. II, we outline the general description of open quantum systems via dynamical maps $\hat{Λ} (t)$ and the time local propagator $\hat{L} (t)$ ⁠, culminating in the two approaches for computing the transient relaxation dynamics and stationary properties. This is followed in Sec. III, where we introduce the thermofield purification and the orthogonal polynomial chain mapping of an infinite non-interacting fermionic bath. In Sec. IV, we introduce the Choi–Jamiolkowski isomorphism, outline the modifications needed to apply this to a fermionic system, and discuss solving the dynamics for non-interacting systems and interacting ones using matrix product states. Building on this, Sec. V reports the main results of this work demonstrating the extraction of $\hat{Λ} (t)$ and $\hat{L} (t)$ for a small spinless Fermi chain both with and without interactions and the single impurity Anderson model. In Sec. VI, we conclude and provide an outlook for future work.

II. NON-MARKOVIAN DYNAMICS

The generic setup that we analyze is a fermionic system with Hilbert space

H_{S}

connected to two non-interacting electronic reservoirs α = {L, R}, defined by equilibrium temperatures T_α = 1/β_α and chemical potentials μ_α (taking k_B = 1 = ℏ throughout) with Hilbert space

H_{α}

⁠. The total Hamiltonian governing the dynamics of this system is given by

\hat{H} = {\hat{H}}_{S} + \sum_{α = L, R} ({\hat{H}}_{S α} + {\hat{H}}_{α}),

(1)

where

{\hat{H}}_{S}

is the system Hamiltonian,

{\hat{H}}_{α}

is the Hamiltonian of bath α, and

{\hat{H}}_{S α}

describes the coupling between them. We restrict our analysis to Hamiltonians

\hat{H}

that conserve fermion number

\hat{N} + \sum_{α} {\hat{N}}_{α}

⁠, where

\hat{N}

and

{\hat{N}}_{α}

are the total particle number operators for the system and bath α, respectively. We study quench dynamics where the initial state of the full setup is taken to be of the product form

{\hat{ρ}}_{tot} (t_{0}) = \hat{ρ} (t_{0}) {\hat{ρ}}_{B}

⁠, where

{\hat{ρ}}_{B}

is the composite initial thermal state of both baths and t₀ is the initial time. The state of the system at a later time t is then found by unitarily time evolving the system + bath and then tracing out the system bath as

\begin{align} \hat{ρ} (t) & = \hat{Λ} (t - t_{0}) [\hat{ρ} (t_{0})], \\ = {Tr}_{B} (\hat{U} (t, t_{0}) \hat{ρ} (t_{0}) {\hat{ρ}}_{B} {\hat{U}}^{†} (t, t_{0})), \end{align}

(2)

where

U (t, t_{0}) = e^{- i \hat{H} (t - t_{0})}

⁠, as shown in Fig. 1. From this construction,

\hat{Λ} (t - t_{0})

is manifestly a completely positive trace preserving (CPTP) map. A related exact description of open quantum system dynamics is given by the time convolutionless master equation,^49–51

\frac{\partial \hat{ρ} (τ)}{\partial τ} = \hat{L} (τ) [\hat{ρ} (τ)],

(3)

with time interval τ = t − t₀ and

\hat{L} (τ) = \frac{d}{d τ} [\hat{Λ} (τ)] {\hat{Λ}}^{- 1} (τ),

(4)

from which

\hat{Λ} (τ) = T e^{\int_{0}^{τ} \hat{L} (τ^{'}) d τ^{'}}

⁠, where

T

is the time-ordering operator. It should be noted that

\hat{L} (τ)

becomes undefined in the long-time limit as

\hat{Λ} (τ)

converges to a projector onto the steady state and thus has no inverse.

FIG. 1.

Schematic of the dynamical map Λ̂(τ). The system S and bath B start in a product initial state ρ̂(t0)ρ̂B. The system and bath then undergo joint unitary evolution for a time duration τ after which the bath is traced out. This construction manifestly generates a CPTP map on the system’s initial state.48

View large Download slide

Schematic of the dynamical map $\hat{Λ} (τ)$ ⁠. The system S and bath B start in a product initial state $\hat{ρ} (t_{0}) {\hat{ρ}}_{B}$ ⁠. The system and bath then undergo joint unitary evolution for a time duration τ after which the bath is traced out. This construction manifestly generates a CPTP map on the system’s initial state.⁴⁸

We consider the situation where the system approaches a unique steady state in the long-time limit, i.e.,

\lim_{τ \to \infty} \hat{Λ} (τ) [\hat{ρ} (0)] = \hat{ρ} (\infty)

for any initial state

\hat{ρ} (0)

⁠. A difference between the above-mentioned non-Markovian description and Markovian quantum dynamics is that both

\hat{Λ} (τ)

and

\hat{L} (τ)

are time-dependent, despite the global Hamiltonian

\hat{H}

being time-independent. The instantaneous steady state, or the time-dependent fixed point, can be defined either as the eigenoperator of

\hat{Λ} (τ)

with eigenvalue 1 or as the eigenoperator of

\hat{L} (τ)

with eigenvalue 0. Neither time-dependent fixed point may correspond to

\hat{ρ} (\infty)

at short times, but both approach

\hat{ρ} (\infty)

at long times. This leads to the existence of three dynamical timescales,

τ_{m}^{Λ}

⁠,

τ_{m}^{L}

⁠, and τ_re, given by

‖ \hat{Λ} (τ) [\hat{ρ} (0)] - \hat{ρ} (\infty) ‖ < ϵ, \forall τ \geq τ_{re},

(5)

‖ \hat{Λ} (τ) [\hat{ρ} (\infty)] - \hat{ρ} (\infty) ‖ < ϵ, \forall τ \geq τ_{m}^{Λ},

(6)

‖ \hat{L} (τ) [\hat{ρ} (\infty)] ‖ < ϵ \forall τ \geq τ_{m}^{L},

(7)

where ϵ is some arbitrarily small tolerance and

‖ \hat{A} ‖

is the norm of

\hat{A}

⁠. Here, τ_re defines the timescale of relaxation for an arbitrary initial state

\hat{ρ} (0)

⁠. In contrast,

τ_{m}^{L}

defines the timescale over which a steady state can be well defined and

τ_{m}^{Λ}

defines the timescale for which

\hat{ρ} (\infty)

is dynamically invariant. These timescales are ordered as

τ_{m}^{L} \leq τ_{m}^{Λ} \leq τ_{re}

⁠, such that given access to

\hat{Λ} (τ)

⁠, and thus

\hat{L} (τ)

⁠, for

τ > τ_{m}^{L}

⁠, we can extract

\hat{ρ} (\infty)

from dynamics over a memory timescale, which can be orders of magnitude smaller than τ_re. Often, it is numerically unfeasible to time-evolve up to τ_re due to the dynamical growth of entanglement, but by exploiting

\hat{Λ} (τ)

directly, we can bypass this and accelerate the extraction of the steady state with significantly shorter time evolution.

Crucially, we can also extract full transient dynamics with minimal computational cost in this dynamical map approach given the following assumption. If we assume the full propagator converges up to an error O(ϵ) on the timescale

τ_{m}^{L}

⁠,

\hat{L} (τ) = {\hat{L}}_{m} + O (ϵ)

⁠, and

\forall τ \geq τ_{m}^{L}

⁠, then we have the following factorized form for the dynamical map:

\hat{Λ} (τ) \approx e^{(τ - τ_{m}^{L}) {\hat{L}}_{m}} \hat{Λ} (τ_{m}^{L}) \forall τ \geq τ_{m}^{L} .

(8)

This decomposition represents a phenomenon called initial slippage, which has been widely studied^52–56 and is closely associated with the assumption of a finite memory time for the environment correlations.^52,57 It also leads to an interesting anomalous thermal relaxation called the non-Markovian quantum Mpemba effect.⁵⁸ If we assume the dynamics is analytic and the system is finite-sized, then we can use the key result in Ref. 29,

{\hat{ρ}}_{n τ + t_{1}} = \underset{n t i m e s}{\underset{︸}{\hat{Λ} (τ) [\dots [\hat{Λ} (τ)}} [\hat{Λ} (t_{1} - t_{0}) \hat{ρ} (t_{0})]]] \dots],

‖ \hat{ρ} (n τ + t_{1}) - {\hat{ρ}}_{n τ + t_{1}} ‖ = ϵ (τ), ϵ (τ) decays with τ,

(9)

where t₀ ≤ t₁ < τ and

τ_{m}^{Λ} ≪ τ \leq τ_{r e}

⁠. This can be thought of as an extension of repeated interaction models,²⁹ where each interaction of the system with the bath takes place over an extended timescale to capture non-Markovian effects rather than the standard instantaneous interaction repetition.

In our approach, we can directly extract the dynamical map

\hat{Λ} (τ)

via the Choi–Jamiolkowski isomorphism and through Eq. , the propagator

\hat{L} (τ)

⁠. Our focus in Sec. V is on small systems S comprising only a few fermionic modes, so given

\hat{Λ} (τ)

and

\hat{L} (τ)

are superoperators acting only on system operators both can be constructed fully as a matrices. With full access to

\hat{Λ} (τ)

and

\hat{L} (τ)

⁠, it is instructive to compute their spectral decomposition. Using the notation

(A | ◦ ≔ {Tr}_{S} A^{†} ◦

and

| B) ≔ B

to denote the Hilbert–Schmidt scalar product (A|B) = Tr_S A^†B, their spectral decomposition follow as

\hat{Λ} (τ) = \sum_{i = 1}^{d^{2}} Λ_{i} (τ) | g_{i} (τ)) ({\bar{g}}_{i} (τ) |,

(10)

\hat{L} (τ) = \sum_{i = 1}^{d^{2}} L_{i} (τ) | h_{i} (τ)) ({\bar{h}}_{i} (τ) |,

(11)

where

({\bar{g}}_{i} | g_{j}) = δ_{ij} = ({\bar{h}}_{i} | h_{j})

and d is the dimension of the system S Hilbert space. Here, (⁠

{\bar{g}}_{i} |

and |g_i) are distinct left and right eigenvectors of

\hat{Λ} (τ)

with the same eigenvalue Λ_i and likewise is true of (⁠

{\bar{h}}_{i} |

and |h_i) for

\hat{L} (τ)

⁠. As Λ(τ) is a CPTP map, we have |Λ_i| ≤ 1 with a leading eigenvalue Λ₁ = 1 with

| g_{1} (τ)) = {\hat{ρ}}^{Λ}_{FP} (τ)

corresponding to its time-dependent fixed point,

\hat{Λ} (τ) [{\hat{ρ}}_{FP}^{Λ} (τ)] = {\hat{ρ}}_{FP}^{Λ} (τ) .

(12)

We also have

R e (L_{i} (τ)) \leq 0

⁠, with a leading eigenvalue

L_{1} (τ) = 0

corresponding to

| h_{i} (τ)) = {\hat{ρ}}^{L}_{FP} (τ)

⁠, obeying

\hat{L} (τ) [{\hat{ρ}}_{FP}^{L} (τ)] = 0 .

(13)

Here, Λ_i(τ) can be interpreted as how much each mode

| g_{i} (τ))

has been damped by time τ, with

L_{i} (τ)

giving the corresponding instantaneous damping rate. It should be noted that the modes themselves are also evolving, a behavior not present in Markovian dynamics. In some cases where the dynamics is strongly non-Markovian, the decay of |Λ_i(τ)| is non monotonic and Λ_i(τ) = 0 can oscillate about zero before fully decaying away. This corresponds to a mode being suppressed and re-excited caused by information backflow from the bath to the system. At these times,

\hat{L} (τ)

becomes ill-defined, but this poses no issues for times other than these isolated points,^50,59 and in fact,

\hat{L} (τ) \hat{Λ} (τ)

is finite for all times.⁶⁰ In these cases, the memory times must be chosen with care. More details regarding this subtle issue can be found in Refs. 92 and 93.

While the stationary state can also be formally computed as the τ → ∞ limit of either of these fixed points, a key observation based on the PReB approach is that practically we only require that $τ > τ_{m}^{Λ}$ or $τ > τ_{m}^{L}$ ⁠, respectively. Equation (8) or (9) can then be used to calculate the transient dynamics for $τ > τ_{m}^{L}, τ ≫ τ_{m}^{Λ}$ to stationarity. We now proceed to outline the details of this approach.

III. MODELLING AN INFINITE BATH

We consider a setup as shown in Fig. 2 with a system coupled to two baths each being a continuum of non-interacting spinless fermionic modes⁶¹ and together governed by

\begin{align} {\hat{H}}_{B} + {\hat{H}}_{S B} & = \sum_{α = L, R} (\int_{- D}^{D} ω {\hat{b}}_{α}^{†} (ω) {\hat{b}}_{α} (ω) d ω \\ + \int_{- D}^{D} \sqrt{J_{α} (ω)} ({\hat{b}}_{α}^{†} (ω) {\hat{Q}}_{α} + {\hat{Q}}_{α}^{†} {\hat{b}}_{α} (ω)) d ω) . \end{align}

(14)

Here, for bath α, we have its spectral density

J_{α} (ω)

⁠, the system operator coupling to it

{\hat{Q}}_{α}

⁠, and its canonical fermionic creation and annihilation operators

{\hat{b}}_{α}^{†} (ω), {\hat{b}}_{α} (ω)

obeying

{{\hat{b}}_{α}^{†} (ω), {\hat{b}}_{α} (ω^{'})} = δ (ω - ω^{'})

⁠. It should be noted that

{\hat{Q}}_{α}

and

{\hat{b}}_{α} (ω)

also anticommute,

{{\hat{Q}}_{α}, {\hat{b}}_{α^{'}} (ω} = 0

⁠. We parameterize

J_{α} (ω)

via its total coupling strength,

Γ_{α} = \frac{1}{2 D} \int_{- D}^{D} 2 π J_{α} (ω) d ω,

(15)

where D is its bandwidth of the baths, assumed to be identical. The initial state of the baths is a product state

{\hat{ρ}}_{B} = \prod_{α} {\hat{ρ}}_{α}

⁠, where

{\hat{ρ}}_{α} = e^{- β_{α} ({\hat{H}}_{α} - μ_{α} {\hat{N}}_{α})} / Z_{α}

(16)

is a thermal state and Z_α is the partition function for bath α.

FIG. 2.

View large Download slide

We consider a system S comprising a small number of fermionic modes coupled to two infinite thermal Fermi baths α = {L, R} with inverse temperatures β_α, spectral density $J_{α}$ ⁠, and coupling strength Γ_α.

In order to represent the equilibrium state of the bath as a pure state, we use thermofield purification in which the finite temperature of any given bath is encoded in two different baths at zero temperature.³⁴ For simplicity, consider a single bath comprising a set of discrete modes governed by a non-interacting Hamiltonian

{\hat{H}}_{B} = \sum_{k} ϵ_{k} {\hat{b}}_{k}^{†} {\hat{b}}_{k}

at inverse temperature β and chemical potential μ, with each mode coupled to some system operator

\hat{Q}

via a hybridization Hamiltonian

{\hat{H}}_{S B} = \sum_{k} (v_{k} {\hat{Q}}^{†} {\hat{b}}_{k} + v_{k}^{*} {\hat{b}}_{k}^{†} \hat{Q})

⁠. The thermal state follows as

\begin{align} {\hat{ρ}}_{β} & = \frac{e^{- β ({\hat{H}}_{B} - μ {\hat{N}}_{B})}}{Tr (e^{- β ({\hat{H}}_{B} - μ {\hat{N}}_{B})})} \\ = \prod_{k} ((1 - f_{k}) (1 - {\hat{b}}_{k}^{†} {\hat{b}}_{k}) + f_{k} {\hat{b}}_{k}^{†} {\hat{b}}_{k}), \end{align}

(17)

where

f_{k} = {(1 + e^{β (ϵ_{k} - μ)})}^{- 1}

is the Fermi factor for the kth mode. By introducing an auxiliary mode

{\hat{A}}_{k}

for each bath mode

{\hat{b}}_{k}

⁠, we can define the thermofield double state |Ω_β⟩ such that the thermal expectation value of an operator

\hat{O}

acting on the system is given by

〈 Ω_{β} | \hat{O} | Ω_{β} 〉 = Tr ({\hat{ρ}}_{β} \hat{O})

⁠. This required state is then

| Ω_{β} 〉 = \prod_{k} (\sqrt{1 - f_{k}} {\hat{A}}_{k}^{†} + \sqrt{f_{k}} {\hat{b}}_{k}^{†}) | v a c 〉,

(18)

for a bath with a star-geometry, as shown in Fig. 3(a). Since purification is invariant to unitary transformations, |Ω_β⟩ is not unique. The form in Eq. has the advantage of being an eigenstate of the total particle number allowing its conservation to be exploited in numerical calculations.

FIG. 3.

View large Download slide

Illustration for a single bath, depicted as discrete modes, of the transformation into two chains. (a) Thermofield purification creates entangled states of each bath eigenmode and its corresponding auxiliary mode. Each bath eigenmode interacts with the system giving a star geometry, while the auxiliary modes are uncoupled. (b) Bath and auxiliary modes are mixed creating a set of empty “0” and filled “1” modes. All these modes couple to the system. (c) After tridiagonalizing the empty and filled modes separately, we arrive at a two chain geometry better suited for tensor network calculations.

By mixing the bath and auxiliary modes via a rotation,

(\begin{matrix} {\hat{C}}_{0, k} \\ {\hat{C}}_{1, k} \end{matrix}) = (\begin{matrix} \cos (θ_{k}) & - \sin (θ_{k}) \\ \sin (θ_{k}) & \cos (θ_{k}) \end{matrix}) (\begin{matrix} {\hat{b}}_{k} \\ {\hat{A}}_{k} \end{matrix}),

(19)

where

\sin (θ_{k}) = \sqrt{f_{k}}

⁠, we map |Ω_β⟩ to a Fock state,

| Ω_{β} 〉 = \prod_{k} {\hat{C}}_{1, k}^{†} | v a c 〉,

(20)

in which the “0” modes

{\hat{C}}_{0, k}

are completely empty while the “1” modes

{\hat{C}}_{1, k}

are completely filled. As the physical and auxiliary modes are decoupled in the Hamiltonian, we are free to choose an arbitrary Hamiltonian for the auxiliary modes themselves. A convenient choice is to set

{\hat{H}}_{A} = \sum_{k} ϵ_{k} {\hat{A}}_{k}^{†} {\hat{A}}_{k}

⁠, which leads to a mixed mode Hamiltonian,

{\hat{H}}_{A} + {\hat{H}}_{B} = \sum_{k} ϵ_{k} ({\hat{C}}_{0, k}^{†} {\hat{C}}_{0, k} + {\hat{C}}_{1, k}^{†} {\hat{C}}_{1, k}) .

(21)

Due to the mixing, all modes now couple to the system as

\begin{align} {\hat{H}}_{S B} & = \sum_{k} \sqrt{f_{k}} (v_{k} {\hat{Q}}^{†} {\hat{C}}_{1, k} + v_{k}^{*} {\hat{C}}_{1, k}^{†} \hat{Q}) \\ - \sum_{k} \sqrt{1 - f_{k}} (v_{k} {\hat{Q}}^{†} {\hat{C}}_{0, k} + v_{k}^{*} {\hat{C}}_{0, k}^{†} \hat{Q}), \end{align}

(22)

giving a modified star geometry bath shown in Fig. 3(b) and in which the thermal properties of the initial state have now been encoded.

Moving back to the continuum with two baths α, we have

v_{k} \to \sqrt{J_{α} (ω)}

⁠, f_k → f_α(ω), and

{\hat{C}}_{j, k} \to {\hat{C}}_{α j} (ω)

⁠, with the equivalent thermofield mixed Hamiltonian,

\begin{align} \hat{H} & = {\hat{H}}_{S} + \sum_{α j} \int_{- D}^{D} d ω ω {\hat{C}}_{α j}^{†} (ω) {\hat{C}}_{α j} (ω) \\ + \sqrt{J_{α j} (ω)} [{\hat{C}}_{α j}^{†} (ω) {\hat{Q}}_{α} + {\hat{Q}}_{α}^{†} {\hat{C}}_{α j} (ω)], \end{align}

(23)

where the empty and filled spectral densities are

J_{α 0} (ω) = (1 - f_{α} (ω)) J_{α} (ω),

(24)

J_{α 1} (ω) = f_{α} (ω) J_{α} (ω),

(25)

and

| Ω_{β} 〉 = \prod_{α} \prod_{- D}^{D} {\hat{C}}_{α 1}^{†} (ω) | v a c 〉

is a continuous Fock state. To render this problem practicable for numerical calculations, the continuous baths are discretized and truncated into a finite number of modes. In some applications, this may be accomplished by assuming some linear or logarithmic binning in the energy axis of the baths.⁶² Here, our focus is on capturing the time-evolution as accurately as possible up to some prescribed time. This is accomplished by performing a continuous mode tridiagonalization of the two zero temperature star geometry baths into two one-dimensional tight-binding chains, as shown in Fig. 3(c). In doing this, we define a new set of fermionic operators

{\hat{c}}_{α j, n}

as⁶³^,

{\hat{C}}_{α j} (ω) = \sum_{n} U_{α j, n} (ω) {\hat{c}}_{α j, n},

(26)

where

U_{α j, n} (ω) = \frac{\sqrt{J_{α j} (ω)} π_{α j, n} (ω)}{ρ_{α j, n}},

(27)

for j = 0, 1 and π_αj,n(ω) are monic orthogonal polynomials that obey

\int_{- D}^{D} d ω J_{α j} (ω) π_{α j, n} (ω) π_{α j, m} (ω) = ρ_{α j, n}^{2} δ_{n m},

(28)

with

ρ_{α j, n}^{2} = \int_{- D}^{D} d ω J_{α j} (ω) π_{α j, n}^{2} (ω) .

(29)

Using a finite cutoff of M modes for the mapping, we arrive at

\begin{align} \hat{H} & = {\hat{H}}_{S} + \sum_{α j} [ρ_{α j, 0} ({\hat{Q}}_{α}^{†} {\hat{c}}_{α j, 0} + {\hat{c}}_{α j, 0}^{†} {\hat{Q}}_{α}) \\ + \sum_{n = 0}^{M} γ_{α j, n} {\hat{c}}_{α j, n}^{†} {\hat{c}}_{α j, n} + (\sqrt{β_{α j, n + 1}} {\hat{c}}_{α j, n + 1}^{†} {\hat{c}}_{α j, n} + H . c .)], \end{align}

(30)

where the coefficients γ_αj,n and β_αj,n are defined through the recurrence relation,

π_{α j, n + 1} (k) = (k - γ_{α j, n}) π_{α j, n} (k) - β_{α j, n} π_{α j, n - 1} (k),

with π_αj,−1(k) = 0. These chain parameters are generated using the ORTHPOL package.⁶⁴ Generically, they are found to quickly converge to constants γ_αj,n → γ_αj and β_αj,n → β_αj. For a weighting function

J_{α j} (ω)

with support [a, b], it can be shown⁶⁵ that these coefficients converge as γ_αj,n → (a + b)/2 and β_αj,n → (b − a)²/16 for n → ∞. In Fig. 4, the chain coefficients β_j,n, γ_j,n are shown for a single bath with a semi-elliptical spectral function defined as Eq. and μ = 0 for various temperatures. At zero temperature,

J_{j} (ω)

has support [0, D], [−D, 0] for j = 0, 1, respectively, and for T > 0 both have support [−D, D]. For the case shown, we see that the chain coefficients converge after only a few sites to the asymptotic values expected from theory. By employing the Lieb–Robinson bound²⁸ sites further than ∼τβ_αj in the αj chain have a negligible effect on the system dynamics up to time τ, giving a well-defined measure of the length of bath chains needed to accurately capture the dynamics. In this sense, the discretization of the baths generated by orthogonal polynomials up to some length M is an exact description of the infinite bath up to this finite time.

FIG. 4.

View large Download slide

Convergence of chain coefficients for the orthogonal polynomial mapping for a single bath with μ = 0, Γ = 0.05D and a semi-elliptical spectral function defined as Eq. (43). (a) Energies and (b) couplings along the empty “0” (solid) and fully occupied “1” (dashed) chains. Due to the symmetric spectral density $J (ω) = J (- ω)$ ⁠, we have β_0,n = β_1,n and γ_0,n = −γ_1,n.

IV. EXTRACTING THE DYNAMICAL MAP

Our key methodological proposal here is to compute and exploit directly the dynamical map $\hat{Λ} (τ)$ for open systems. To accomplish this, we use the Choi–Jamiolkowski isomorphism, which provides equivalence between quantum states and superoperators.^35,36

A. Choi–Jamiolkowski isomorphism

Mirroring the thermofield construction for the bath, the isomorphism starts by adding auxilliary replica degrees of freedom A_S for the system S as well. Consider a d-level system each with basis states |s⟩ for which the application of any CPTP map

\hat{Λ}

to an initial state

{\hat{ρ}}_{0}

can be drawn, as shown in Fig. 5(a). The trick of the isomorphism is to initialize the system and its auxiliary replica in a perfectly correlated maximally entangled state,

| Φ^{+} 〉 = \frac{1}{\sqrt{d}} \sum_{s = 1}^{d} {| s 〉}_{S} \otimes {| s 〉}_{A_{S}},

(31)

whose diagrammatic form is shown in Fig. 5(b). Since

\hat{Λ}

is completely positive, applying it to the system alone is guaranteed to generate a valid density operator in the space of system + replica operators⁶⁶ given by

{\hat{ρ}}^{Λ} = (\hat{Λ} \otimes 1) {| Φ^{+} 〉 〈 Φ^{+} |},

(32)

as shown in Fig. 5(c). Armed with

{\hat{ρ}}^{Λ}

⁠, the action of the map

\hat{Λ}

then follows

\hat{Λ} [\hat{ρ_{0}}] = d {Tr}_{A_{S}} (1 \otimes {\hat{ρ}}_{0}^{T} {\hat{ρ}}^{Λ}),

(33)

as shown in Fig. 5(d). The isomorphism straightforwardly extends to maps applied to the N system sites by using the state,

| Φ_{N}^{+} 〉 = \frac{1}{\sqrt{d^{N}}} \otimes_{j = 1}^{N} (\sum_{s_{j} = 1}^{d} {| s_{j} 〉}_{S} \otimes {| s_{j} 〉}_{A_{S}}),

(34)

where each site is maximally entangled with its own auxiliary replica. The application of the isomorphism to many-body fermionic systems also follows analogously once some subtle details are taken care of as we now describe.

FIG. 5.

Tensor network type diagrams of the Choi–Jamiolkowski isomorphism. (a) The application of a CPTP map Λ̂[ρ̂0] can be viewed as a cap shaped tensor encompassing ρ̂. Internally, the cap tensor could be decomposed into a conjugation by Kraus operators K̂l as shown. (b) Owing to its perfect correlations, the maximally entangled state |Φ+⟩ from Eq. (31) is a bent line. (c) The Choi–Jamiolkowski isomorphism follows from inputting |Φ+⟩ into Λ̂, which diagrammatically corresponds to bending the input legs of Λ̂ outward to form ρ̂Λ as in Eq. (32). (d) The application of the map Λ̂[ρ̂0] is extracted from ρ̂Λ by contracting the appropriate legs consistent with panel (a) and inserting a factor of the dimension d to get the graphical version of Eq. (33).

View large Download slide

Tensor network type diagrams of the Choi–Jamiolkowski isomorphism. (a) The application of a CPTP map $\hat{Λ} [{\hat{ρ}}_{0}]$ can be viewed as a cap shaped tensor encompassing $\hat{ρ}$ ⁠. Internally, the cap tensor could be decomposed into a conjugation by Kraus operators ${\hat{K}}_{l}$ as shown. (b) Owing to its perfect correlations, the maximally entangled state |Φ⁺⟩ from Eq. (31) is a bent line. (c) The Choi–Jamiolkowski isomorphism follows from inputting |Φ⁺⟩ into $\hat{Λ}$ ⁠, which diagrammatically corresponds to bending the input legs of $\hat{Λ}$ outward to form ${\hat{ρ}}^{Λ}$ as in Eq. (32). (d) The application of the map $\hat{Λ} [{\hat{ρ}}_{0}]$ is extracted from ${\hat{ρ}}^{Λ}$ by contracting the appropriate legs consistent with panel (a) and inserting a factor of the dimension d to get the graphical version of Eq. (33).

B. Application to fermionic systems

Take the system S as being described by L spinless fermionic modes

{\hat{s}}_{j}

⁠. We then additionally introduce corresponding auxiliary modes

{\hat{a}}_{j}

for each of them and anticipate that these pairs of modes will be initialized in an entangled state. For a single bath, the setup has the form shown in Fig. 6(a). To construct fermionic Fock states for the setup, a specific one-dimensional mode ordering is needed. We will consider two such orderings. The first is separated ordering,

\begin{array}{l} {{\hat{d}}_{i}}_{sep} & = \{{\hat{s}}_{1}, \dots, {\hat{s}}_{L}, {\hat{a}}_{1}, \dots, {\hat{a}}_{L}, {\hat{c}}_{0, 1},, \dots, {\hat{c}}_{0, M}, {\hat{c}}_{1, 1}, \dots, {\hat{c}}_{1, M}\}, \end{array}

where system, auxiliary replica, and each set of thermofield bath modes are grouped together. The, Fock states of the setup follow as

\begin{array}{l} | n^{S}, n^{A_{S}}, n^{0}, n^{1} 〉 & = {({\hat{s}}_{1}^{†})}^{n_{1}^{S}} \dots {({\hat{s}}_{L}^{†})}^{n_{L}^{S}} {({\hat{a}}_{1}^{†})}^{n_{1}^{A_{S}}} \dots {({\hat{a}}_{L}^{†})}^{n_{L}^{A_{S}}} \\ \times {({\hat{c}}_{0, 1}^{†})}^{n_{1}^{0}} \dots {({\hat{c}}_{0, M}^{†})}^{n_{M}^{0}} {({\hat{c}}_{1, 1}^{†})}^{n_{1}^{1}} \dots {({\hat{c}}_{1, M}^{†})}^{n_{M}^{1}} | v a c 〉, \end{array}

with n^S = {0, 1}^L being the occupation numbers of the system modes, while

n^{A_{S}}

and n⁰, n¹ are likewise for the replica modes and the “0”/“1” thermofield chains.

FIG. 6.

View large Download slide

(a) For the case of a single bath, our setup comprises the system S, its auxiliary replica A_S and two thermofield chains “0” and “1,” which are empty and filled, respectively. (b) The assumed fermionic mode ordering used in this work. The system and its replicas are interleaved, but separated from the thermofield chains, which are also interleaved together. In this ordering, local interactions between modes now become 3-local (next-nearest-neighbor) and 4-local interactions along the one-dimensional ordering of modes.

One issue for the application of the Choi–Jamiolkowski to fermions is that it formally requires a partial trace over fermionic modes. Since fermionic Fock states are constructed from mode operators with a nonlocal canonical anticommutator algebra, they lack an underlying tensor product structure associated with the conventional partial trace.^67,68 However, it can be shown that if the modes to be retained are a contiguous subset of the assumed mode ordering, and the overall state is number parity symmetric, then the conventional partial trace coincides with the fermionic mode-reduced states, as shown explicitly in Ref. 67. The separated mode ordering accomplishes this for the thermofield bath modes, which are traced out in the definition of $\hat{Λ} (τ)$ in Eq. (2), and for the replica modes of the system which are traced out in Eq. (33) when reconstructing $\hat{Λ} (τ)$ from ${\hat{ρ}}^{Λ} (τ)$ ⁠.

With this ordering in place, the Choi–Jamiolkowski isomorphism proceeds by initializing the setup in the state,

| Ψ_{CJ} 〉 = \sum_{n \in {0, 1}^{L}} \frac{1}{\sqrt{2^{L}}} | n, n, 0, \dots, 0, 1, \dots, 1 〉,

(35)

where the occupation of the system and replica modes are perfectly correlated while the thermofield chains are empty and filled. To extract

\hat{Λ} (τ)

⁠, this initial state is time-evolved as

| Ψ_{CJ} (τ) 〉 = \exp (- i \hat{H} τ) | Ψ_{CJ} 〉

from which we obtain

{\hat{ρ}}^{Λ} (τ) = {Tr}_{01} (| Ψ_{CJ} (τ) 〉 〈 Ψ_{CJ} (τ) |)

⁠.

While the state |Ψ_CJ⟩ is number-parity symmetric, it does not possess a fixed total particle number. When

\hat{H}

itself is number conserving, as it is here, it is numerically advantageous to fully exploit this symmetry. To do this, we use a modified initial state with the system and replica modes having an anti-correlated occupation,

| Ψ_{AC} 〉 = \sum_{n \in {0, 1}^{L}} \frac{1}{\sqrt{2^{L}}} | n, \bar{n}, 0, \dots, 0, 1, \dots, 1 〉,

(36)

where

{\bar{n}}_{j} = n_{j} \oplus 1

denotes the occupation flipping of a mode. This state over the system, replica, and thermofield modes is precisely half-filled. We are permitted to use |Ψ_AC⟩ in place of |Ψ_CJ⟩ and compute

{\hat{ρ}}^{AC} (τ) = {Tr}_{01} (| Ψ_{AC} (τ) 〉 〈 Ψ_{AC} (τ) |)

since there exists a unitary transformation

\hat{P} | Ψ_{AC} 〉 = | Ψ_{CJ} 〉

⁠, which acts exclusively on the replica modes such that

[\hat{P}, \hat{H}] = 0

⁠, as outlined in Appendix A. Consequently,

\begin{array}{l} {\hat{ρ}}^{Λ} (τ) & = {Tr}_{01} (e^{- i \hat{H} τ} | Ψ_{CJ} 〉 〈 Ψ_{CJ} | e^{i \hat{H} τ}), \\ = {Tr}_{01} (e^{- i \hat{H} τ} \hat{P} | Ψ_{AC} 〉 〈 Ψ_{AC} | {\hat{P}}^{†} e^{i \hat{H} τ}), \\ = \hat{P} {\hat{ρ}}^{AC} (τ) {\hat{P}}^{†}, \end{array}

and hence

\hat{P}

can be applied to the reduced state of the system and its replica after the time-evolution with the bath has been performed.

The overall setup generalizes straightforwardly to two baths α = {L, R} by introducing the two thermofield chains for the left bath at the start as

\begin{align} {{\hat{d}}_{i}}_{sep} & = \{{\hat{c}}_{L 0, 1},, \dots, {\hat{c}}_{L 0, M}, {\hat{c}}_{L 1, 1}, \dots, {\hat{c}}_{L 1, M}, {\hat{s}}_{1}, \dots, {\hat{s}}_{L}, {\hat{a}}_{1}, \dots, {\hat{a}}_{L}, \\ {\hat{c}}_{R 0, 1}, \dots, {\hat{c}}_{R 0, M}, {\hat{c}}_{R 1, 1}, \dots, {\hat{c}}_{R 1, M}\}, \end{align}

(37)

Otherwise, the calculation proceeds identically.

C. Non-interacting fermions

If

{\hat{H}}_{S}

is quadratic in mode operators and particle number conserving, then the exact dynamics can be obtained via unitary evolution of the single-particle correlation matrix

C_{ij} (τ) = Tr (\hat{ρ} (τ) {\hat{d}}_{j}^{†} {\hat{d}}_{i})

over all the system, bath, and auxiliary modes

{{\hat{d}}_{i}}

⁠. Given

\hat{H} = \sum_{ij} h_{ij} {\hat{d}}_{i}^{†} {\hat{d}}_{j}

⁠, the time-evolved correlation matrix follows from

C (τ) = e^{i h τ} C (t_{0}) e^{- i h τ},

(38)

as shown in detail in Appendix C. The reduced density matrix of the system and its auxiliary modes is directly obtained from C(τ) via⁶⁹

\begin{align} {\hat{ρ}}^{AC} (τ) & = \det (1 - C^{T} (τ)) \\ \times \exp \{\sum_{ij \in S A_{S}} {[\log (C^{T} (τ)) {(1 - C^{T} (τ))}^{- 1}]}_{ij} {\hat{d}}_{i}^{†} {\hat{d}}_{j}\}, \end{align}

(39)

where SA_S denotes the 2L modes spanning the system and its replica modes. This equation holds provided

{\hat{ρ}}^{AC} (τ)

is block diagonal in the number basis, as is the case here for the particle number symmetric initial state |Ψ^AC⟩. After extracting

{\hat{ρ}}^{Λ} (τ)

⁠, a highly non-trivial verification of its correctness is achieved by comparing to the prediction of Landauer–Büttiker theory, as presented in Sec. V.

D. Matrix product state time evolution

To handle interacting fermions, we employ matrix product state (MPS) techniques. Given a many-body system of N qubits each described by Pauli-Z eigenstates |σ⟩ with σ = ±1, a pure state in its tensor product Hilbert space

{(C^{2})}^{\otimes N}

is represented by an MPS as^70,71

| ψ 〉 = \sum_{σ_{1}, \dots, σ_{N}} A_{1}^{σ_{1}} A_{2}^{σ_{2}} \dots A_{N - 1}^{σ_{N - 1}} A_{N}^{σ_{N}} | σ_{1}, \dots, σ_{N} 〉,

(40)

where

A_{j}^{σ_{j}}

is a tensor formed from the χ_j−1 × χ_j matrix (with χ₀ = χ_N = 1 fixed) indexed by the state σ_j of the jth qubit. A crucial feature of MPS is that the bond dimension χ_j can be directly linked to the entanglement between the two halves of the system it connects.⁷² The one-dimensional connectivity of MPS makes them naturally suited to describing systems with Hamiltonians sharing this geometry, where an area law of entanglement is expected to severely limit the entanglement in its ground state and low-lying excitations.⁷³

The chain-like geometry of the open systems that we arrived at in Eq. and Fig. 6(a) thus appears to be well suited to MPS. However, the separated mode ordering introduced in Sec. IV B is suboptimal for MPS, so we instead adopt interleaved ordering. For a two-bath setup, this is given by

\begin{align} {{\hat{d}}_{i}}_{int} & = \{{\hat{c}}_{L 1, M}, {\hat{c}}_{L 0, M}, \dots, {\hat{c}}_{L 1, 1}, {\hat{c}}_{L 0, 1}, \\ {\hat{s}}_{1}, {\hat{a}}_{1}, \dots, {\hat{s}}_{L}, {\hat{a}}_{L}, {\hat{c}}_{R 0, 1}, {\hat{c}}_{R 1, 1}, \dots, {\hat{c}}_{R 0, M}, {\hat{c}}_{R 1, M}\} . \end{align}

(41)

In this ordering, the thermofield chains for each baths are interleaved, as are the system and replica modes. This mode ordering has two key advantages. First, it minimizes the length of interactions by avoiding terms coupling modes spanning across separated chains. This results in the nearest-neighbor interactions in the chains becoming next-nearest-neighbor interactions in the MPS, as shown in Fig. 6(b). Locality is especially important since spinless fermionic modes are described in systems of qubits using the Jordan–Wigner transformation.^74,75 Consequently long-ranged fermionic couplings become long-ranged multi-body operators in the equivalent qubit Hamiltonian that are potentially deleterious due to dynamical entanglement growth. Second, interleaved ordering minimizes the MPS bond dimension required to capture the initial system + replica mode entangled state in Eq. by locating the entangled mode pairs adjacent to one another. In contrast, the separated ordering results in “rainbow”-like state with a bond dimension χ_j ∼ 2^L.

The MPS computation of the time-evolution of |Ψ_AC⟩ in interleaved mode ordering proceeds using the two-site time-dependent variational principle (2TDVP) algorithm^76–78 exploiting particle number conservation. While the effective qubit Hamiltonian is short-ranged, it does have terms beyond nearest-neighbor so it is essential to use a global subspace expansion to reduce the projection error.⁷⁶ Our calculations were all performed using the ITensor package.⁷⁷ Since the interleaved ordering still separates the bath modes from the system and its replica modes, the computation of ${\hat{ρ}}^{AC} (τ)$ is unchanged. The final extraction of $\hat{Λ} (τ)$ although requires that the system and replica modes are swapped back into separated ordering, as detailed in Appendix A. This step is potentially inefficient for large L and avoiding it is the subject of ongoing work.⁷⁹ For the small systems considered here, this presents no issues.

V. NUMERICAL RESULTS

We now demonstrate the effectiveness of the dynamical map approach for paradigmatic fermionic systems. For the tensor network calculations, we use a time step of δτ = 0.1/D, a 2TDVP truncation cutoff of 10⁻¹⁰–10⁻¹² to time-evolve baths comprising of M ∼ 10² modes, which generates bond dimensions on the order χ ∼ 10³.

Spinless Fermi chain: consider a chain of spinless fermions with a Hamiltonian given by

{\hat{H}}_{S} = \sum_{i = 1}^{L - 1} t_{c} ({\hat{s}}_{i + 1}^{†} {\hat{s}}_{i} + {\hat{s}}_{i}^{†} {\hat{s}}_{i + 1}) + U {\hat{s}}_{i}^{†} {\hat{s}}_{i} {\hat{s}}_{i + 1}^{†} {\hat{s}}_{i + 1},

(42)

where each end of the chain is coupled to a bath,

{\hat{Q}}_{L} = {\hat{s}}_{1}, {\hat{Q}}_{R} = {\hat{s}}_{L}

⁠. We study the dynamics of a small L = 3 site tight-binding chain with the system–bath couplings switched on at τ = 0. For simplicity, we assume a semi-elliptical spectral function for both baths,

J_{α} (ω) = \frac{2 Γ_{α}}{π^{2}} \sqrt{1 - {(ω / D)}^{2}},

(43)

with Γ_L = Γ_R and β_L = β_R, but with μ_L = −μ_R giving a non-equilibrium setup.

We begin by considering the spectral properties of the map and propagator. In Fig. 7, the time-dependent eigenvalues |Λ_i(τ)| for $\hat{Λ} (τ)$ and $R e (L_{i} (τ))$ for $\hat{L} (τ)$ are plotted for both non-interacting and interacting cases, respectively. They satisfy the properties of complete positivity outlined earlier in Sec. II. For |Λ_i(τ)|, we see the fixed unit eigenvalue, corresponding to the time-dependent fixed point ${\hat{ρ}}_{FP}^{Λ} (τ)$ ⁠, with the other eigenvalues decaying to zero with τ as $\hat{Λ} (τ)$ converges to a projector onto the stationary state. For the decay rates $R e (L_{i} (τ))$ ⁠, we also see a fixed zero eigenvalue corresponding to its time-dependent fixed point ${\hat{ρ}}_{FP}^{L} (τ)$ ⁠. The initial slippage at early times is clear in the transient behavior of $R e (L_{i} (τ))$ ⁠, after which they rapidly converge to constant values. The remaining oscillations observed correspond precisely to the bandwidth D and hence occur on a much faster timescale than all other timescales of the problem. These are an artifact of the bath spectral function having a sharp band edge and disappear if this edge is smoothed over a larger bandwidth. Aside from degeneracies among the decay modes being lifted, we observe the same generic behavior shown in Figs. 7(c) and 7(d) for the interacting case. Given the non-interacting results are computed exactly, this serves as a first demonstration of our MPS methodology.

FIG. 7.

Eigenvalues |Λi(τ)| of the dynamical map Λ̂(τ) and Re(Li(τ)) of the propagator L̂(τ) shown against time τ. Panels (a) and (b) show the non-interacting case and panels (c) and (d) show the interacting case U = 0.05D. Parameters: L = 3, tc = 0.02D, ΓL = ΓR = 0.05D, βL = βR = 1/D, μL = −μR = 0.2D.

View large Download slide

Eigenvalues |Λ_i(τ)| of the dynamical map $\hat{Λ} (τ)$ and $R e (L_{i} (τ))$ of the propagator $\hat{L} (τ)$ shown against time τ. Panels (a) and (b) show the non-interacting case and panels (c) and (d) show the interacting case U = 0.05D. Parameters: L = 3, t_c = 0.02D, Γ_L = Γ_R = 0.05D, β_L = β_R = 1/D, μ_L = −μ_R = 0.2D.

Next, we consider some simple observables of the system, namely, the density and the current, defined by

{\hat{N}}_{i} = {\hat{s}}_{i}^{†} {\hat{s}}_{i}

and

{\hat{J}}_{i} = i ({\hat{s}}_{i}^{†} {\hat{s}}_{i + 1} - {\hat{s}}_{i + 1}^{†} {\hat{s}}_{i})

⁠. The steady state current for non-interacting systems can be computed from Landauer–Büttiker theory as

⟨ {\hat{J}}_{\infty} ⟩ = \frac{1}{2 π} \int_{- D}^{D} d ω τ (ω) [f_{L} (ω) - f_{R} (ω)],

(44)

where f_α(ω) denotes the Fermi–Dirac distribution for bath α and τ(ω) is the transmission function. For details on the computation of τ(ω), see Appendix B. Figure 8(a) shows the evolution of the fixed point currents between site 2 and 3 of the system

{⟨ \hat{J} ⟩}_{L} \equiv Tr ({\hat{J}}_{2} {\hat{ρ}}_{FP}^{L} (τ))

and

{⟨ \hat{J} ⟩}_{Λ} \equiv Tr ({\hat{J}}_{2} {\hat{ρ}}_{FP}^{Λ} (τ))

⁠, which show fast convergence to

{⟨ \hat{J} ⟩}_{\infty}

within τ ≈ 20/D. This is in contrast to the much longer thermalization of a totally mixed initial state

{⟨ \hat{J} ⟩}_{ex} \equiv Tr ({\hat{J}}_{2} \hat{ρ} (τ))

shown in Fig. 8(b), which is still far from equilibrium by the end of the simulation at τ = 60/D. This demonstrates how our approach can offer a significant speedup in determining the stationary state compared to direct evolution.

FIG. 8.

View large Download slide

Evolution of a non-interacting Fermi chain in a non-equilibrium setup. (a) Convergence of the two time-dependent fixed point currents along with the exact steady state solution ${⟨ \hat{J} ⟩}_{\infty}$ ⁠, calculated using Landauer–Büttiker theory. (b) The evolution of a totally mixed initial state (blue) is compared to the solution obtained from Eqs. (8) and (9) using memory times $τ_{m}^{L} D = 20$ and $τ_{m}^{Λ} D = 40$ ⁠. Panels (c) and (d) show the same analysis for the density on the first site. Parameters: L = 3, U = 0, t_c = 0.02D, Γ_L = Γ_R = 0.05D, β_L = β_R = 1/D, μ_L = −μ_R = 0.2D.

In addition to efficiently extracting steady states, we can extract all transient dynamics accurately using Eq. (8) once $\hat{L} (τ)$ has converged. As shown in Fig. 8(b), this approximation works remarkably well using $τ_{m}^{L} D = 20$ recovering the exact dynamics up to τ = 60/D and extrapolating beyond this demonstrating the steady state limit is attained. Applying the PreB approach of Eq. (9) with a larger memory time, $τ_{m}^{Λ} D = 40$ recovers the transient dynamics reasonably well, but with a visible error expected due to the requirement that PReB should only be applied for times $τ ≫ τ_{m}^{Λ}$ ⁠. Figures 8(c) and 8(d) show the same analysis for the density, where the hierarchy of timescales is again clear. The error in the PreB approach shown in Fig. 8(b) is several orders of magnitude smaller than the scale of Fig. 8(d) and so is not visible, showing both Eqs. (8) and (9) recovering the correct transient dynamics and steady state limit.

The non-interacting results reported were computed using the exact numerical solution outlined in Sec. IV C. We have confirmed that these are fully reproduced by tensor network calculations (not shown). We now move to an identical setting with U > 0 interacting chain, where tensor network methods are a necessity and proceed with the same analysis. Figures 9(a) and 9(b) show the current with interactions identically to Figs. 8(a) and 8(b) for U = 0. The values of the currents are suppressed due to interactions, but otherwise the behavior is qualitatively identical. In particular, there is a fast convergence for ${⟨ \hat{J} ⟩}_{L}, {⟨ \hat{J} ⟩}_{Λ}$ and a contrasting slow relaxation of ${⟨ \hat{J} ⟩}_{ex} \equiv Tr ({\hat{J}}_{2} \hat{ρ} (τ))$ ⁠. Again, the extrapolations using Eqs. (8) and (9) recover the transient dynamics accurately, with Eq. (9) showing a small but visible slippage error. Together, these results demonstrate that our dynamical map approach is viable for interacting Fermi systems. Our next example considers a seminal interacting problem with a well-known emergent timescale.

FIG. 9.

View large Download slide

Evolution of an interacting Fermi chain in a non-equilibrium setup. (a) Convergence of the two time-dependent fixed point currents along with the steady state solution ${⟨ \hat{J} ⟩}_{\infty}$ inferred from the largest τ calculated. (b) The evolution of a totally mixed initial state (blue) is compared to the predicted solution obtained from Eqs. (8) and (9) using memory times $τ_{m}^{L} D = 20$ and $τ_{m}^{Λ} D = 40$ ⁠. Parameters: L = 3, U = 0.05D, t_c = 0.02D, Γ_L = Γ_R = 0.05D, β_L = β_R = 1/D, μ_L = −μ_R = 0.2D.

Anderson impurity model: we now consider the particle-hole symmetric single impurity Anderson model (SIAM), with Hamiltonian,

{\hat{H}}_{S} = U {\hat{s}}_{↑}^{†} {\hat{s}}_{↑} {\hat{s}}_{↓}^{†} {\hat{s}}_{↓} - \frac{U}{2} {\hat{s}}_{↑}^{†} {\hat{s}}_{↑} - \frac{U}{2} {\hat{s}}_{↓}^{†} {\hat{s}}_{↓} .

(45)

The Coulomb interaction U ≠ 0 gives rise to strong correlations both in and out of equilibrium causing complex behavior, such as the Kondo effect, where the impurity entangles with the reservoir electrons screening its local moment.^80–82 The crossover between local moments and singlet formation is given by the Kondo temperature T_K, which for a single bath with spectral density

J (ω)

is given by⁸²

T_{K} = \sqrt{2 U J (0)} \exp (- \frac{U}{8 J (0)}) .

(46)

Correspondingly, the Kondo cloud bound state forms over a non-perturbatively long timescale τ_K ∼ 1/T_K.

To simulate the SIAM, we treat it as a system of two spinless fermion modes, one for each spin projections σ = {↑, ↓}, which are coupled only through the Coulomb interaction.²⁴ Each spin projection σ is then coupled via

{\hat{Q}}_{σ} = {\hat{s}}_{σ}

to its own bath. We focus on the same equilibrium setup as in Ref. 83, so the σ baths have identical spectral functions

J (ω)

⁠, taken to be a flat spectral density with smoothed edges,

J (ω) = \frac{Γ}{2 π (1 + e^{ν (ω - D)}) (1 + e^{- ν (ω + D)})},

(47)

characterized by ν, which we fix to be ν = 100/D, and identical coupling strength Γ, inverse temperature β, and chemical potential μ. We initialize the impurity in the state

\hat{ρ} (t_{0}) = | ↑ 〉 〈 ↑ |

and fix the interaction to be strong as U = 4Γ.

The first scenario that we consider is a moderately coupled bath⁸⁴ where the associated Kondo temperature is T_K ≈ 0.02D and with a bath temperature well below this. Being in the Kondo regime, we expect the timescale τ_K = 50/D to manifest as a long memory time. In Fig. 10(a), the time-dependent fixed points ${\hat{ρ}}_{FP}^{L}, {\hat{ρ}}_{FP}^{Λ}$ are shown to converge very rapidly by τ ≈ 10/D to the thermal state suggesting a short memory time. In Fig. 10(b), we recover and match the accuracy of the results from the influence tensor approach used in Ref. 83 and similarly in Ref. 85 for the relaxation of all diagonal components of the impurities reduced density matrix. Even by τ = 100/D, the evolution of the spin polarized initial state is still far from equilibrium, reflecting the slow dynamical formation of a spin singlet in the strongly interacting Kondo regime expected for β > 1/T_K. None of this physics is visible in the behavior of the time-dependent fixed points but is instead revealed by the time dependence of the decay modes. In Fig. 10(c), the decay modes fail to converge by τ = 100/D, showing a long-lived time dependence of the propagator. It should be noted that only the slowest degenerate decay modes are shown as all other modes rapidly decay away. This is reflected by the accuracy of the dynamical predictions given by Eqs. (8) and (9). In Fig. 10(b), we show that when a small memory time $τ_{m}^{L}$ is used, such that $\hat{L}$ is far from converged, neither recover the transient dynamics. By increasing the memory time $τ_{m}^{L}$ ⁠, we can achieve accurate predictions, as seen by the solution of Eq. (8) in Fig. 10(d), while Eq. (9) still suffers visible slippage error. This points to the decay modes of $\hat{L} (τ), \hat{Λ} (τ)$ having a strong dependence on system–bath correlations associated with the formation of the Kondo singlet screening cloud, which is controlled by the much longer timescale τ_K. Similar findings are reported in Refs. 43 and 86 when analyzing the memory kernel in the Kondo regime.

FIG. 10.

Evolution of the diagonal elements of the impurity density matrix ρ̂kk, where k = {∅, ↑, ↓, ⇕} for the single impurity Anderson model. For strong interactions, U = 4Γ at (i) moderate coupling Γ = 0.2D and low temperature β = 500/D Kondo physics is expected, while at (ii) weak coupling Γ = 0.02D and high-temperature β = 0.1/D, the dynamics is expected to be more Markovian. (a) Convergence of the two time-dependent fixed points ρ̂FPΛ,ρ̂FPL is shown for case (i). It should be noted that due to particle-hole symmetry, ρ̂↑↑=ρ̂↓↓,ρ̂∅∅=ρ̂⇕⇕ for both fixed points. The totally mixed high-temperature fixed point for (ii) is also shown. (b) The direct evolution of the initial state |↑⟩⟨↑| for case (i) compared to the solution obtained from Eq. (8) (dashed) and Eq. (9) (dotted–dashed) using memory times τmLD=15 and τmΛD=15, respectively. The result of repeated applications of Λ̂(τmΛ) on the initial state is plotted stroboscopically (circles). (c) The evolution of the decay modes Re(Li(τ)) for case (i), showing a long lived time dependence. (d) The same as panel (b) but a longer memory time τmLD=75=τmΛD is used, where Eq. (8) now accurately extrapolates the correct time evolution. (e) The evolution of the decay modes for case (ii), showing convergence to fixed points up to bandwidth oscillation. (f) For case (ii), the direct solution of the initial state is again compared to the solution obtained from Eq. (8) (dashed) obtained from using τm LD=20. The solution using Eq. (9) is only plotted stroboscopically (circles) as it is indistinguishable from Eq. (8) for this case.

View large Download slide

Evolution of the diagonal elements of the impurity density matrix ${\hat{ρ}}_{k k}$ ⁠, where k = {∅, ↑, ↓, ⇕} for the single impurity Anderson model. For strong interactions, U = 4Γ at (i) moderate coupling Γ = 0.2D and low temperature β = 500/D Kondo physics is expected, while at (ii) weak coupling Γ = 0.02D and high-temperature β = 0.1/D, the dynamics is expected to be more Markovian. (a) Convergence of the two time-dependent fixed points ${\hat{ρ}}_{FP}^{Λ}, {\hat{ρ}}_{FP}^{L}$ is shown for case (i). It should be noted that due to particle-hole symmetry, ${\hat{ρ}}_{↑ ↑} = {\hat{ρ}}_{↓ ↓}, {\hat{ρ}}_{\emptyset \emptyset} = {\hat{ρ}}_{⇕ ⇕}$ for both fixed points. The totally mixed high-temperature fixed point for (ii) is also shown. (b) The direct evolution of the initial state |↑⟩⟨↑| for case (i) compared to the solution obtained from Eq. (8) (dashed) and Eq. (9) (dotted–dashed) using memory times $τ_{m}^{L} D = 15$ and $τ_{m}^{Λ} D = 15$ ⁠, respectively. The result of repeated applications of $\hat{Λ} (τ_{m}^{Λ})$ on the initial state is plotted stroboscopically (circles). (c) The evolution of the decay modes $R e (L_{i} (τ))$ for case (i), showing a long lived time dependence. (d) The same as panel (b) but a longer memory time $τ_{m}^{L} D = 75 = τ_{m}^{Λ} D$ is used, where Eq. (8) now accurately extrapolates the correct time evolution. (e) The evolution of the decay modes for case (ii), showing convergence to fixed points up to bandwidth oscillation. (f) For case (ii), the direct solution of the initial state is again compared to the solution obtained from Eq. (8) (dashed) obtained from using $τ_{m}^{L} D = 20$ ⁠. The solution using Eq. (9) is only plotted stroboscopically (circles) as it is indistinguishable from Eq. (8) for this case.

In the second scenario, we consider a weakly coupled bath where the associated Kondo temperature is T_K ≈ 0.001D and the bath temperature is far above this. In this high-temperature regime of a high-symmetry point, the eventual thermal state is extremely close to the totally mixed state ${\hat{ρ}}_{\infty} = (| \emptyset 〉 〈 \emptyset | + | ↑ 〉 〈 ↑ | + | ↓ 〉 〈 ↓ | + | ⇕ 〉 〈 ⇕ |) / 4$ ⁠. Correspondingly, we find (not shown) that the time-dependent fixed points ${\hat{ρ}}_{FP}^{Λ} (τ), {\hat{ρ}}_{FP}^{L} (τ)$ essentially immediately converge. While this is suggestive of a vanishingly short memory time, again, like with the Kondo regime, the predictive abilities and spectral properties of $\hat{Λ} (τ)$ and $\hat{L} (τ)$ give a sharper characterization. Figure 10(e) shows the decay rates $R e (L_{i} (τ))$ ⁠, showing an initial slippage before converging to fixed values up to bandwidth oscillations. The lack of a longer timescale is indicative that Kondo physics is not present. From this, we take $τ_{m}^{L} D = 20$ and see that Eq. (8) accurately reproduces the true dynamics with minimal error, as shown in Fig. 10(f). There is also strong agreement with Eq. (9), whose stroboscopic predictions are displayed confirming that the physics of this regime can be fully captured within a very short time. We, therefore, see that the success or failure of the slippage approximations Eqs. (8) and (9) in capturing the transient dynamics serves as a very effective indicator of non-Markovian behavior, allowing the Kondo regime to be revealed.

VI. CONCLUSIONS

Dynamical maps $\hat{Λ} (τ)$ and time-local propagators $\hat{L} (τ)$ are central but often formal objects in the theory of open systems. Here, we have shown that by combining the thermofield doubling and orthogonal polynomial chain mappings with the Choi–Jamiolkowski isomorphism, both $\hat{Λ} (τ)$ and $\hat{L} (τ)$ are readily accessible objects for numerical calculations of open Fermi systems. We work with pure state unitary evolution for which the numerical errors are well understood and positivity is manifestly preserved, in contrast to influence functional tensor methods.⁸⁵ For systems that are mildly non-Markovian with a short memory time, this dynamical map approach can accelerate the computation of stationary properties compared to directly simulating relaxation in the long-time limit. Moreover, converged $\hat{Λ} (τ)$ and $\hat{L} (τ)$ allow for predictions of the transient dynamics for all longer times up to stationarity fully characterizing the open system. We have demonstrated this for a spinless Fermi chain driven out of equilibrium by two baths, both with and without interactions, confirming the viability and effectiveness of this approach. We also investigated the equilibration of the SIAM hosting more complex Kondo physics. Here, the thermal stationary state, even at strong coupling could be computed rapidly, but attempts to capture the transient relaxation dynamics were found to be challenging in the Kondo regime. Consequently, this provides a robust means of quantifying the memory time and revealing non-trivial emergent effects known to distinguish different regimes of the SIAM. In the examples, we reach the steady state rapidly, but there could be additional benefits to working directly in the steady state limit as done in Refs. 87 and 88. The Choi–Jamiolkowski isomorphism is more general than the MPS methodology presented here and could also be combined with other methods, such as the HEOM approach and process tensor approaches. Here, we have focused exclusively on small systems comprising only L = 2–3 spinless fermionic modes, where both $\hat{Λ} (τ)$ and $\hat{L} (τ)$ can be fully constructed as matrices. A future work⁷⁹ is looking into utilizing tensor network methods for the entire calculations, where both these objects are computed with a compressed description. This should enable much larger open systems to be addressed with a dynamical map approach.

ACKNOWLEDGMENTS

S.R.C. acknowledges the financial support from UK’s Engineering and Physical Sciences Research Council (EPSRC) under Grant No. EP/T028424/1. A.P. acknowledges funding from Seed Grant from IIT Hyderabad, Project No. SG/IITH/F331/2023-24/SG-169. A.P. also acknowledges funding from Japan International Cooperation Agency (JICA) Friendship 2.0 Research Grant and from Finnish Indian Consortia for Research and Education (FICORE).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

David J. Strachan: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Archak Purkayastha: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal). Stephen R. Clark: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

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

APPENDIX A: FINDING THE FERMIONIC CORRECTION $\hat{P}$

In terms of qubit operators, the unitary

\hat{P}

is given by

\hat{P} = \prod_{i = 1}^{N_{s}} {\hat{σ}}_{a_{i}}^{x}

⁠, where

{\hat{σ}}_{a_{i}}^{x}

is the bit flip operator on the ith replica mode. In terms of fermionic operators, this is given by

\hat{P} = \prod_{i = 1}^{L} \sum_{k = 1}^{2 i - 1} e^{i π {\hat{N}}_{k}} ({\hat{a}}_{i}^{†} + {\hat{a}}_{i})

⁠, where k is summed over both system and replica modes. As each system mode is perfectly anti-correlated with its replica mode initially, we can express any system number operator in terms of the corresponding replica’s number operator,

{\hat{N}}_{k} = 1 - {\hat{N}}_{k + 1}

⁠, where k is even, such that

\hat{P}

can be fully represented through replica modes. As the Hamiltonian only contains quadratic or quartic terms of system operators, we have

[\hat{P}, \hat{H}] = 0

⁠. If the interleaved mode ordering

{{\hat{d}}_{i}}_{int}

is used, the initial state is

| Ψ_{CJ}^{int} 〉 = \prod_{i = 1}^{L} ({\hat{s}}_{i}^{†} + {\hat{a}}_{i}^{†}) | 0, 0, 0, \dots 0, 1 \dots, 1 〉 .

(A1)

We map

| Ψ_{CJ}^{int} 〉

to a state |Ψ_CJ⟩ using an additional unitary operation

| Ψ_{CJ} 〉 = {\hat{P}}_{2} | Ψ_{CJ}^{int} 〉

⁠, which is given by a series of swap gates that are not fermionic. In practice, we do not swap the bath modes back into separated ordering for the isomorphism as they are all traced out as a contiguous block in either case. In terms of fermionic operators, we have

{\hat{P}}_{2} = \prod_{i = 1}^{L} {\hat{T}}_{i}^{q} = \prod_{i = 1}^{L} e^{- i π {\hat{N}}_{i} {\hat{N}}_{i + 1}} {\hat{T}}_{i}^{f}

⁠, where

T_{i}^{q}, {\hat{T}}_{i}^{f}

are qubit and fermionic swap gates respectively, acting on modes i, i + 1, where

{\hat{N}}_{i}

is the number operator for mode i. These swaps will be between a replica and a system mode, or two replica modes. In either case, the phase correction can be expressed in terms of only replica modes due to the perfect correlation between

{\hat{s}}_{i}

and

{\hat{a}}_{i}

such that

[{\hat{P}}_{2}, \hat{H}] = 0

⁠.

APPENDIX B: LANDAUER–BÜTTIKER THEORY

Here, we briefly outline how to compute the currents in Landauer–Büttiker theory, which act as our point of comparison for non-interacting system steady states. In the continuum limit for macroscopic baths, the particle and energy currents in the absence of system interactions are given by^89,90

J_{L B}^{P} = \frac{1}{2 π} \int_{- D}^{D} d ω τ (ω) [f_{L} (ω) - f_{R} (ω)],

(B1)

and

J_{L B}^{E} = \frac{1}{2 π} \int_{- D}^{D} d ω ω τ (ω) [f_{L} (ω) - f_{R} (ω)],

(B2)

where f_α(ω) denotes the Fermi–Dirac distribution for bath α and τ(ω) is the transmission function of the system. This can be calculated in terms of the non-equilibrium Green’s function.^90,91 In our case, this is given by

G (ω) = M^{- 1} (ω),

(B3)

with

M (ω) = ω 1 - h_{S} - Σ^{(1)} (ω) - Σ^{(L)} (ω),

(B4)

where the only nonzero elements of the self-energy matrices of the leads Σ^(j) (ω) are

{[Σ^{(j)}]}_{j j} (ω) = P \int_{- D}^{D} d ω^{'} \frac{J_{j} (ω^{'})}{ω^{'} - ω} - i π J_{j} (ω),

(B5)

using

J_{1} (ω) = J_{L} (ω)

⁠,

J_{L} (ω) = J_{R} (ω)

⁠, and h_S is defined via

{\hat{H}}_{S} = \sum_{i, j = 1}^{L} {(h_{S})}_{ij} {\hat{s}}_{i}^{†} {\hat{s}}_{j}

⁠. If the system Hamiltonian is of the form

h_{S} = \sum_{j = 1}^{L} ϵ_{j} {\hat{s}}_{j}^{†} {\hat{s}}_{j} + \sum_{j = 1}^{L - 1} t_{i} ({\hat{s}}_{j + 1}^{†} {\hat{s}}_{j} + h.c.),

(B6)

the transmission function is given by

τ (ω) = 4 π^{2} J_{L} (ω) J_{R} (ω) | {[G (ω)]}_{1 D} |^{2} = \frac{J_{L} (ω) J_{R} (ω)}{| \det [M] |^{2}} \prod_{i = 1}^{L - 1} | t_{i} |^{2} .

(B7)

APPENDIX C: CORRELATION MATRIX PROPAGATION

Here, we prove Eq. . Defining

U (τ) = e^{- i \hat{H} τ}

⁠, we have

C_{ij} (τ) = 〈 ψ (τ) | {\hat{d}}_{j}^{†} {\hat{d}}_{i} | ψ (τ) 〉 = 〈 ψ (0) | U^{†} (τ) {\hat{d}}_{j}^{†} {\hat{d}}_{i} U (τ) | ψ (0) 〉 .

To evaluate this, first consider how

\hat{H}

acts on

{\hat{d}}_{j}^{†}

⁠,

\hat{H} {\hat{d}}_{j}^{†} = \sum_{k l} H_{k l} {\hat{d}}_{k}^{†} {\hat{d}}_{l} {\hat{d}}_{j}^{†} = \sum_{k l} H_{k l} {\hat{d}}_{k}^{†} (δ_{l j} - {\hat{d}}_{j}^{†} {\hat{d}}_{l}) = \sum_{k} H_{kj} {\hat{d}}_{k}^{†} .

(C1)

We can now evaluate

U^{†} {\hat{d}}_{j}^{†}

as

U^{†} {\hat{d}}_{j}^{†} = e^{iHt} {\hat{d}}_{j}^{†} = e^{i \sum_{k} H_{kj} t} {\hat{d}}_{k}^{†} = \sum_{k} U_{kj}^{†} {\hat{d}}_{k}^{†} .

(C2)

Substituting this into Eq. gives the result,

\begin{align} C_{ij} (t) & = 〈 ψ (0) | \sum_{k l} U_{kj}^{†} {\hat{d}}_{k}^{†} {(U_{l i}^{†} {\hat{d}}_{l}^{†})}^{†} | ψ (0) 〉 = \sum_{k l} U_{kj}^{†} U_{i l} 〈 ψ (0) | {\hat{d}}_{k}^{†} {\hat{d}}_{l} | ψ (0) 〉 \\ = \sum_{k l} U_{kj}^{†} U_{i l} C_{l k} (0) . \end{align}

(C3)

Therefore, we have

C (t) = U C (0) U^{†} .

(C4)

REFERENCES

1.

M. A.

Nielsen

and

I. L.

Chuang

,

Quantum Computation and Quantum Information

(Cambridge University Press,

2001

).

Google Scholar

2.

M.

Josefsson

,

A.

Svilans

,

A. M.

Burke

,

E. A.

Hoffmann

,

S.

Fahlvik

,

C.

Thelander

,

M.

Leijnse

, and

H.

Linke

, “

A quantum-dot heat engine operating close to the thermodynamic efficiency limits

,”

Nat. Nanotechnol.

13

,

920

–

924

(

2018

).

https://doi.org/10.1038/s41565-018-0200-5

Google Scholar

Crossref

PubMed

3.

A.

Ronzani

,

B.

Karimi

,

J.

Senior

,

Y.-C.

Chang

,

J. T.

Peltonen

,

C.

Chen

, and

J. P.

Pekola

, “

Tunable photonic heat transport in a quantum heat valve

,”

Nat. Phys.

14

,

991

–

995

(

2018

).

https://doi.org/10.1038/s41567-018-0199-4

Google Scholar

Crossref

4.

N.

Mosso

,

H.

Sadeghi

,

A.

Gemma

,

S.

Sangtarash

,

U.

Drechsler

,

C.

Lambert

, and

B.

Gotsmann

, “

Thermal transport through single-molecule junctions

,”

Nano Lett.

19

,

7614

–

7622

(

2019

).

https://doi.org/10.1021/acs.nanolett.9b02089

Google Scholar

Crossref

PubMed

5.

N.

Lambert

,

Y.-N.

Chen

,

Y.-C.

Cheng

,

C.-M.

Li

,

G.-Y.

Chen

, and

F.

Nori

, “

Quantum biology

,”

Nat. Phys.

9

,

10

–

18

(

2013

).

https://doi.org/10.1038/nphys2474

Google Scholar

Crossref

6.

H.

Van Amerongen

,

R.

Van Grondelle

et al,

Photosynthetic Excitons

(

World Scientific

,

2000

).

Google Scholar

Crossref

7.

V.

May

and

O.

Kühn

,

Charge and Energy Transfer Dynamics in Molecular Systems

(

John Wiley & Sons

,

2023

).

Google Scholar

8.

A.

Levy

and

R.

Kosloff

, “

The local approach to quantum transport may violate the second law of thermodynamics

,”

Europhys. Lett.

107

,

20004

(

2014

).

https://doi.org/10.1209/0295-5075/107/20004

Google Scholar

Crossref

9.

J. T.

Stockburger

and

T.

Motz

, “

Thermodynamic deficiencies of some simple lindblad operators

,”

Fortschr. Phys.

65

,

1600067

(

2016

).

https://doi.org/10.1002/prop.201600067

Google Scholar

Crossref

10.

H.

Wichterich

,

M. J.

Henrich

,

H.-P.

Breuer

,

J.

Gemmer

, and

M.

Michel

, “

Modeling heat transport through completely positive maps

,”

Phys. Rev. E

76

,

031115

(

2007

).

https://doi.org/10.1103/physreve.76.031115

Google Scholar

Crossref

11.

M. T.

Mitchison

and

M. B.

Plenio

, “

Non-additive dissipation in open quantum networks out of equilibrium

,”

New J. Phys.

20

,

033005

(

2018

).

https://doi.org/10.1088/1367-2630/aa9f70

Google Scholar

Crossref

12.

G.

Stefanucci

and

R.

van Leeuwen

,

Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction

(

Cambridge University Press

,

2013

).

Google Scholar

Crossref

13.

J.-S.

Wang

,

B. K.

Agarwalla

,

H.

Li

, and

J.

Thingna

, “

Nonequilibrium Green’s function method for quantum thermal transport

,”

Front. Phys.

9

,

673

–

697

(

2014

).

https://doi.org/10.1007/s11467-013-0340-x

Google Scholar

Crossref

14.

N.

Talarico

,

S.

Maniscalco

, and

N. L.

Gullo

, “

Study of the energy variation in many-body open quantum systems: Role of interactions in the weak and strong coupling regimes

,”

Phys. Rev. B

101

,

045103

(

2020

).

https://doi.org/10.1103/physrevb.101.045103

Google Scholar

Crossref

15.

A.

Strathearn

,

P.

Kirton

,

D.

Kilda

,

J.

Keeling

, and

B. W.

Lovett

, “

Efficient non-Markovian quantum dynamics using time-evolving matrix product operators

,”

Nat. Commun.

9

,

3322

(

2018

).

https://doi.org/10.1038/s41467-018-05617-3

Google Scholar

Crossref

PubMed

16.

J.

Thoenniss

,

A.

Lerose

, and

D. A.

Abanin

, “

Nonequilibrium quantum impurity problems via matrix-product states in the temporal domain

,”

Phys. Rev. B

107

,

195101

(

2023

).

https://doi.org/10.1103/physrevb.107.195101

Google Scholar

Crossref

17.

R.

Feynman

and

F.

Vernon

, “

The theory of a general quantum system interacting with a linear dissipative system

,”

Ann. Phys.

24

,

118

–

173

(

1963

).

https://doi.org/10.1016/0003-4916(63)90068-x

Google Scholar

Crossref

18.

A.

Lerose

,

M.

Sonner

, and

D. A.

Abanin

, “

Influence matrix approach to many-body Floquet dynamics

,”

Phys. Rev. X

11

,

021040

(

2021

).

https://doi.org/10.1103/physrevx.11.021040

Google Scholar

Crossref

19.

M.

Cygorek

,

M.

Cosacchi

,

A.

Vagov

,

V. M.

Axt

,

B. W.

Lovett

,

J.

Keeling

, and

E. M.

Gauger

, “

Simulation of open quantum systems by automated compression of arbitrary environments

,”

Nat. Phys.

18

,

662

–

668

(

2022

).

https://doi.org/10.1038/s41567-022-01544-9

Google Scholar

Crossref

20.

F. A.

Pollock

,

C.

Rodríguez-Rosario

,

T.

Frauenheim

,

M.

Paternostro

, and

K.

Modi

, “

Non-Markovian quantum processes: Complete framework and efficient characterization

,”

Phys. Rev. A

97

,

012127

(

2018

).

https://doi.org/10.1103/physreva.97.012127

Google Scholar

Crossref

21.

G. E.

Fux

,

D.

Kilda

,

B. W.

Lovett

, and

J.

Keeling

, “

Tensor network simulation of chains of non-Markovian open quantum systems

,”

Phys. Rev. Res.

5

,

033078

(

2023

).

https://doi.org/10.1103/physrevresearch.5.033078

Google Scholar

Crossref

22.

M.

Cygorek

and

E. M.

Gauger

, “

ACE: A general-purpose non-Markovian open quantum systems simulation toolkit based on process tensors

,”

J. Chem. Phys.

161

,

074111

(

2024

).

https://doi.org/10.1063/5.0221182

Google Scholar

Crossref

PubMed

23.

Y.

Tanimura

, “

Numerically ‘exact’ approach to open quantum dynamics: The hierarchical equations of motion (HEOM)

,”

J. Chem. Phys.

153

,

020901

(

2020

).

https://doi.org/10.1063/5.0011599

Google Scholar

Crossref

24.

L.

Kohn

and

G. E.

Santoro

, “

Quench dynamics of the Anderson impurity model at finite temperature using matrix product states: Entanglement and bath dynamics

,”

J. Stat. Mech.

2022

,

063102

(

2022

).

https://doi.org/10.1088/1742-5468/ac729b

Google Scholar

Crossref

25.

J.

Eisert

,

M.

Cramer

, and

M. B.

Plenio

, “

Colloquium: Area laws for the entanglement entropy

,”

Rev. Mod. Phys.

82

,

277

(

2010

).

https://doi.org/10.1103/revmodphys.82.277

Google Scholar

Crossref

26.

I.

de Vega

,

U.

Schollwöck

, and

F. A.

Wolf

, “

How to discretize a quantum bath for real-time evolution

,”

Phys. Rev. B

92

,

155126

(

2015

).

https://doi.org/10.1103/physrevb.92.155126

Google Scholar

Crossref

27.

M. P.

Woods

,

M.

Cramer

, and

M. B.

Plenio

, “

Simulating bosonic baths with error bars

,”

Phys. Rev. Lett.

115

,

130401

(

2015

).

https://doi.org/10.1103/physrevlett.115.130401

Google Scholar

Crossref

PubMed

28.

M. P.

Woods

and

M. B.

Plenio

, “

Dynamical error bounds for continuum discretisation via Gauss quadrature rules—A Lieb-Robinson bound approach

,”

J. Math. Phys.

57

,

022105

(

2016

).

https://doi.org/10.1063/1.4940436

Google Scholar

Crossref

29.

A.

Purkayastha

,

G.

Guarnieri

,

S.

Campbell

,

J.

Prior

, and

J.

Goold

, “

Periodically refreshed baths to simulate open quantum many-body dynamics

,”

Phys. Rev. B

104

,

045417

(

2021

).

https://doi.org/10.1103/physrevb.104.045417

Google Scholar

Crossref

30.

J.

Rau

, “

Relaxation phenomena in spin and harmonic oscillator systems

,”

Phys. Rev.

129

,

1880

–

1888

(

1963

).

https://doi.org/10.1103/physrev.129.1880

Google Scholar

Crossref

31.

V.

Scarani

,

M.

Ziman

,

P.

Štelmachovič

,

N.

Gisin

, and

V.

Bužek

, “

Thermalizing quantum machines: Dissipation and entanglement

,”

Phys. Rev. Lett.

88

,

097905

(

2002

).

https://doi.org/10.1103/physrevlett.88.097905

Google Scholar

Crossref

PubMed

32.

M.

Ziman

,

P.

Štelmachovič

,

V.

Bužek

,

M.

Hillery

,

V.

Scarani

, and

N.

Gisin

, “

Diluting quantum information: An analysis of information transfer in system-reservoir interactions

,”

Phys. Rev. A

65

,

042105

(

2002

).

https://doi.org/10.1103/physreva.65.042105

Google Scholar

Crossref

33.

S.

Campbell

and

B.

Vacchini

, “

Collision models in open system dynamics: A versatile tool for deeper insights?

,”

Europhys. Lett.

133

,

60001

(

2021

).

https://doi.org/10.1209/0295-5075/133/60001

Google Scholar

Crossref

34.

Y.

Takahashi

and

H.

Umezawa

, “

Thermo field dynamics

,”

Int. J. Mod. Phys. B

10

,

1755

–

1805

(

1996

).

https://doi.org/10.1142/s0217979296000817

Google Scholar

Crossref

35.

M.-D.

Choi

, “

Completely positive linear maps on complex matrices

,”

Linear Algebra Appl.

10

,

285

–

290

(

1975

).

https://doi.org/10.1016/0024-3795(75)90075-0

Google Scholar

Crossref

36.

A.

Jamiołkowski

, “

Linear transformations which preserve trace and positive semidefiniteness of operators

,”

Rep. Math. Phys.

3

,

275

–

278

(

1972

).

https://doi.org/10.1016/0034-4877(72)90011-0

Google Scholar

Crossref

37.

K.

Nestmann

and

M. R.

Wegewijs

, “

General connection between time-local and time-nonlocal perturbation expansions

,”

Phys. Rev. B

104

,

155407

(

2021

).

https://doi.org/10.1103/physrevb.104.155407

Google Scholar

Crossref

38.

D.

Chruściński

, “

Dynamical maps beyond Markovian regime

,”

Phys. Rep.

992

,

1

–

85

(

1992

).

https://doi.org/10.1016/j.physrep.2022.09.003

Google Scholar

Crossref

39.

R.

Zwanzig

, “

Ensemble method in the theory of irreversibility

,”

J. Chem. Phys.

33

,

1338

–

1341

(

1960

).

https://doi.org/10.1063/1.1731409

Google Scholar

Crossref

40.

S.

Nakajima

, “

On quantum theory of transport phenomena: Steady diffusion

,”

Prog. Theor. Phys.

20

,

948

–

959

(

1958

).

https://doi.org/10.1143/ptp.20.948

Google Scholar

Crossref

41.

Q.

Shi

and

E.

Geva

, “

A new approach to calculating the memory kernel of the generalized quantum master equation for an arbitrary system–bath coupling

,”

J. Chem. Phys.

119

,

12063

–

12076

(

2003

).

https://doi.org/10.1063/1.1624830

Google Scholar

Crossref

42.

M.-L.

Zhang

,

B. J.

Ka

, and

E.

Geva

, “

Nonequilibrium quantum dynamics in the condensed phase via the generalized quantum master equation

,”

J. Chem. Phys.

125

,

044106

(

2006

).

https://doi.org/10.1063/1.2218342

Google Scholar

Crossref

43.

G.

Cohen

and

E.

Rabani

, “

Memory effects in nonequilibrium quantum impurity models

,”

Phys. Rev. B

84

,

075150

(

2011

).

https://doi.org/10.1103/physrevb.84.075150

Google Scholar

Crossref

44.

G.

Cohen

,

E. Y.

Wilner

, and

E.

Rabani

, “

Generalized projected dynamics for non-system observables of non-equilibrium quantum impurity models

,”

New J. Phys.

15

,

073018

(

2013

).

https://doi.org/10.1088/1367-2630/15/7/073018

Google Scholar

Crossref

45.

J.

Cerrillo

and

J.

Cao

, “

Non-Markovian dynamical maps: Numerical processing of open quantum trajectories

,”

Phys. Rev. Lett.

112

,

110401

(

2014

).

https://doi.org/10.1103/physrevlett.112.110401

Google Scholar

Crossref

PubMed

46.

F. A.

Pollock

and

K.

Modi

, “

Tomographically reconstructed master equations for any open quantum dynamics

,”

Quantum

2

,

76

(

2018

).

https://doi.org/10.22331/q-2018-07-11-76

Google Scholar

Crossref

47.

F. A.

Pollock

,

E.

Gull

,

K.

Modi

, and

G.

Cohen

, “

Reduced dynamics of full counting statistics

,”

SciPost Phys.

13

,

027

(

2022

).

https://doi.org/10.21468/scipostphys.13.2.027

Google Scholar

Crossref

48.

W. F.

Stinespring

, “

Positive functions on C*-algebras

,”

Proc. Am. Math. Soc.

6

,

211

–

216

(

1955

).

https://doi.org/10.1090/s0002-9939-1955-0069403-4

Google Scholar

Crossref

49.

K.

Nestmann

and

C.

Timm

, “

Time-convolutionless master equation: Perturbative expansions to arbitrary order and application to quantum dots

,” arXiv:1903.05132 [cond-mat.mes-hall] (

2019

).

Google Scholar

50.

D.

Chruściński

and

A.

Kossakowski

, “

Non-Markovian quantum dynamics: Local versus nonlocal

,”

Phys. Rev. Lett.

104

,

070406

(

2010

).

https://doi.org/10.1103/physrevlett.104.070406

Google Scholar

Crossref

PubMed

51.

D.

Chruściński

, “

Dynamical maps beyond Markovian regime

,”

Phys. Rep.

992

,

1

–

85

(

2022

).

https://doi.org/10.1016/j.physrep.2022.09.003

Google Scholar

Crossref

52.

V.

Bruch

,

K.

Nestmann

,

J.

Schulenborg

, and

M.

Wegewijs

, “

Fermionic duality: General symmetry of open systems with strong dissipation and memory

,”

SciPost Phys.

11

,

053

(

2021

).

https://doi.org/10.21468/scipostphys.11.3.053

Google Scholar

Crossref

53.

U.

Geigenmüller

,

U.

Titulaer

, and

B.

Felderhof

, “

Systematic elimination of fast variables in linear systems

,”

Physica A

119

,

41

–

52

(

1983

).

https://doi.org/10.1016/0378-4371(83)90144-9

Google Scholar

Crossref

54.

F.

Haake

and

M.

Lewenstein

, “

Adiabatic drag and initial slip in random processes

,”

Phys. Rev. A

28

,

3606

–

3612

(

1983

).

https://doi.org/10.1103/physreva.28.3606

Google Scholar

Crossref

55.

F.

Haake

and

R.

Reibold

, “

Strong damping and low-temperature anomalies for the harmonic oscillator

,”

Phys. Rev. A

32

,

2462

–

2475

(

1985

).

https://doi.org/10.1103/physreva.32.2462

Google Scholar

Crossref

56.

P.

Gaspard

and

M.

Nagaoka

, “

Slippage of initial conditions for the Redfield master equation

,”

J. Chem. Phys.

111

,

5668

–

5675

(

1999

).

https://doi.org/10.1063/1.479867

Google Scholar

Crossref

57.

K.

Nestmann

,

V.

Bruch

, and

M. R.

Wegewijs

, “

How quantum evolution with memory is generated in a time-local way

,”

Phys. Rev. X

11

,

021041

(

2021

).

https://doi.org/10.1103/physrevx.11.021041

Google Scholar

Crossref

58.

D. J.

Strachan

,

A.

Purkayastha

, and

S. R.

Clark

, “

Non-Markovian quantum Mpemba effect

,” arXiv:2402.05756 [quant-ph] (

2024

).

Google Scholar

59.

A. S.

Hegde

,

K. P.

Athulya

,

V.

Pathak

,

J.

Piilo

, and

A.

Shaji

, “

Open quantum dynamics with singularities: Master equations and degree of non-Markovianity

,”

Phys. Rev. A

104

,

062403

(

2021

).

https://doi.org/10.1103/physreva.104.062403

Google Scholar

Crossref

60.

K.

Nestmann

, “

Open quantum systems: Time (non)-locality, fixed points, and renormalization groups

,” Ph.D. dissertation (

RWTH Aachen University

,

Aachen

,

2022

).

Google Scholar

61.

The generalization to spin-1/2 fermions is straightforward, either by carrying an additional spin quantum number or doubling the number of spinless modes. The latter approach will be applied later in Sec. V.

62.

R.

Bulla

,

T. A.

Costi

, and

T.

Pruschke

, “

Numerical renormalization group method for quantum impurity systems

,”

Rev. Mod. Phys.

80

,

395

–

450

(

2008

).

https://doi.org/10.1103/revmodphys.80.395

Google Scholar

Crossref

63.

A. W.

Chin

,

A.

Rivas

,

S. F.

Huelga

, and

M. B.

Plenio

, “

Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials

,”

J. Math. Phys.

51

,

092109

(

2010

).

https://doi.org/10.1063/1.3490188

Google Scholar

Crossref

64.

W.

Gautschi

, “

Orthogonal polynomials (in Matlab)

,”

J. Comput. Appl. Math.

178

,

215

–

234

(

2005

).

https://doi.org/10.1016/j.cam.2004.03.029

Google Scholar

Crossref

65.

R.

Rosenbach

,

J.

Cerrillo

,

S. F.

Huelga

,

J.

Cao

, and

M. B.

Plenio

, “

Efficient simulation of non-Markovian system-environment interaction

,”

New J. Phys.

18

,

023035

(

2016

).

https://doi.org/10.1088/1367-2630/18/2/023035

Google Scholar

Crossref

66.

Conversely, for any nonnegative operator we can also to a corresponding quantum map.

67.

G. G.

Amosov

and

S. N.

Filippov

, “

Spectral properties of reduced fermionic density operators and parity superselection rule

,”

Quantum Inf. Process.

16

,

2

(

2016

).

https://doi.org/10.1007/s11128-016-1467-9

Google Scholar

Crossref

68.

M.

Cirio

,

P.-C.

Kuo

,

Y.-N.

Chen

,

F.

Nori

, and

N.

Lambert

, “

Canonical derivation of the fermionic influence superoperator

,”

Phys. Rev. B

105

,

035121

(

2022

).

https://doi.org/10.1103/physrevb.105.035121

Google Scholar

Crossref

69.

S.-A.

Cheong

and

C. L.

Henley

, “

Many-body density matrices for free fermions

,”

Phys. Rev. B

69

,

075111

(

2004

).

https://doi.org/10.1103/physrevb.69.075111

Google Scholar

Crossref

70.

G.

Vidal

, “

Efficient classical simulation of slightly entangled quantum computations

,”

Phys. Rev. Lett.

91

,

147902

(

2003

).

https://doi.org/10.1103/physrevlett.91.147902

Google Scholar

Crossref

PubMed

71.

D.

Perez-Garcia

,

F.

Verstraete

,

M. M.

Wolf

, and

J. I.

Cirac

, “

Matrix product state representations

,”

Quantum Inf. Comput.

7

,

401

–

430

(

2007

).

https://doi.org/10.26421/qic7.5-6-1

Google Scholar

Crossref

72.

R.

Orús

, “

A practical introduction to tensor networks: Matrix product states and projected entangled pair states

,”

Ann. Phys.

349

,

117

–

158

(

2014

).

https://doi.org/10.1016/j.aop.2014.06.013

Google Scholar

Crossref

73.

M. B.

Hastings

, “

An area law for one-dimensional quantum systems

,”

J. Stat. Mech.: Theory Exp.

2007

,

P08024

.

https://doi.org/10.1088/1742-5468/2007/08/p08024

Crossref

74.

R. H.

Romer

, “

The spin-statistics theorem

,”

Am. J. Phys.

70

,

791

(

2002

).

https://doi.org/10.1119/1.1482064

Google Scholar

Crossref

75.

M. A.

Nielsen

, “

The fermionic canonical commutation relations and the Jordan-Wigner transform

,” https://futureofmatter.com/assets/fermions_and_jordan_wigner.pdf (

2005

).

Google Scholar

76.

M.

Yang

and

S. R.

White

, “

Time-dependent variational principle with ancillary Krylov subspace

,”

Phys. Rev. B

102

,

094315

(

2020

).

https://doi.org/10.1103/physrevb.102.094315

Google Scholar

Crossref

77.

M.

Fishman

,

S.

White

, and

E.

Stoudenmire

, “

The ITensor software library for tensor network calculations

,”

SciPost Phys. Codebases

4

(

2022

).

https://doi.org/10.21468/SciPostPhysCodeb.4

Google Scholar

78.

J.

Haegeman

,

C.

Lubich

,

I.

Oseledets

,

B.

Vandereycken

, and

F.

Verstraete

, “

Unifying time evolution and optimization with matrix product states

,”

Phys. Rev. B

94

,

165116

(

2016

).

https://doi.org/10.1103/physrevb.94.165116

Google Scholar

Crossref

79.

D. J.

Strachan

,

A.

Purkayastha

, and

S. R.

Clark

, “

Tensor network methods for dynamical map extraction of open quantum systems

” (unpublished).

80.

M.

Nuss

,

M.

Ganahl

,

E.

Arrigoni

,

W.

von der Linden

, and

H. G.

Evertz

, “

Nonequilibrium spatiotemporal formation of the Kondo screening cloud on a lattice

,”

Phys. Rev. B

91

,

085127

(

2015

).

https://doi.org/10.1103/physrevb.91.085127

Google Scholar

Crossref

81.

M.

Medvedyeva

,

A.

Hoffmann

, and

S.

Kehrein

, “

Spatiotemporal buildup of the Kondo screening cloud

,”

Phys. Rev. B

88

,

094306

(

2013

).

https://doi.org/10.1103/physrevb.88.094306

Google Scholar

Crossref

82.

P.

Coleman

,

Introduction to Many-Body Physics

(

Cambridge University Press

,

2015

).

Google Scholar

Crossref

83.

J.

Thoenniss

,

M.

Sonner

,

A.

Lerose

, and

D. A.

Abanin

, “

Efficient method for quantum impurity problems out of equilibrium

,”

Phys. Rev. B

107

,

L201115

(

2023

).

https://doi.org/10.1103/physrevb.107.l201115

Google Scholar

Crossref

84.

In Ref. 83, they use two identical baths and a coupling strength Γ_L = Γ_R = D/10 which has an equivalent total coupling strength as a single bath with Γ = D/5.

85.

N.

Ng

,

G.

Park

,

A. J.

Millis

,

G. K.-L.

Chan

, and

D. R.

Reichman

, “

Real-time evolution of Anderson impurity models via tensor network influence functionals

,”

Phys. Rev. B

107

,

125103

(

2023

).

https://doi.org/10.1103/physrevb.107.125103

Google Scholar

Crossref

86.

G.

Cohen

,

E.

Gull

,

D. R.

Reichman

,

A. J.

Millis

, and

E.

Rabani

, “

Numerically exact long-time magnetization dynamics at the nonequilibrium Kondo crossover of the Anderson impurity model

,”

Phys. Rev. B

87

,

195108

(

2013

).

https://doi.org/10.1103/physrevb.87.195108

Google Scholar

Crossref

87.

F.

Schwarz

,

I.

Weymann

,

J.

von Delft

, and

A.

Weichselbaum

, “

Nonequilibrium steady-state transport in quantum impurity models: A thermofield and quantum quench approach using matrix product states

,”

Phys. Rev. Lett.

121

,

137702

(

2018

).

https://doi.org/10.1103/physrevlett.121.137702

Google Scholar

Crossref

PubMed

88.

A.

Erpenbeck

,

E.

Gull

, and

G.

Cohen

, “

Quantum Monte Carlo method in the steady state

,”

Phys. Rev. Lett.

130

,

186301

(

2023

).

https://doi.org/10.1103/physrevlett.130.186301

Google Scholar

Crossref

PubMed

89.

D. A.

Ryndyk

,

Theory of Quantum Transport at Nanoscale: An Introduction

(

Springer

,

2016

), p.

246

, ISBN: 978-3-319-24086-2.

Google Scholar

Crossref

90.

A.

Purkayastha

, “

Classifying transport behavior via current fluctuations in open quantum systems

,”

J. Stat. Mech.: Theory Exp.

2019

,

043101

.

https://doi.org/10.1088/1742-5468/ab02f4

Crossref

91.

A.

Dorda

,

M.

Nuss

,

W.

von der Linden

, and

E.

Arrigoni

, “

Auxiliary master equation approach to nonequilibrium correlated impurities

,”

Phys. Rev. B

89

,

165105

(

2014

).

https://doi.org/10.1103/physrevb.89.165105

Google Scholar

Crossref

92.

V.

Jagadish

,

R.

Srikanth

, and

F.

Petruccione

, “

Noninvertibility and non-Markovianity of quantum dynamical maps

,”

Phys. Rev. A

108

(

4

),

042202

(

2023

).

https://doi.org/10.1103/PhysRevA.108.042202

Google Scholar

Crossref

93.

F.

vom Ende

, “

Almost all quantum channels are diagonalizable

,”

Open Syst. Inf. Dynam.

31

,

2450012

(

2024

); arXiv:2403.19643.

https://doi.org/10.1142/s1230161224500124

Google Scholar

Crossref

Published open access through an agreement with JISC Collections

Extracting dynamical maps of non-Markovian open quantum systems

I. INTRODUCTION

II. NON-MARKOVIAN DYNAMICS

III. MODELLING AN INFINITE BATH

IV. EXTRACTING THE DYNAMICAL MAP

A. Choi–Jamiolkowski isomorphism

B. Application to fermionic systems

C. Non-interacting fermions

D. Matrix product state time evolution

V. NUMERICAL RESULTS

VI. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: FINDING THE FERMIONIC CORRECTION $\hat{P}$

APPENDIX B: LANDAUER–BÜTTIKER THEORY

APPENDIX C: CORRELATION MATRIX PROPAGATION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Extracting dynamical maps of non-Markovian open quantum systems

I. INTRODUCTION

II. NON-MARKOVIAN DYNAMICS

III. MODELLING AN INFINITE BATH

IV. EXTRACTING THE DYNAMICAL MAP

A. Choi–Jamiolkowski isomorphism

B. Application to fermionic systems

C. Non-interacting fermions

D. Matrix product state time evolution

V. NUMERICAL RESULTS

VI. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: FINDING THE FERMIONIC CORRECTION P̂

APPENDIX B: LANDAUER–BÜTTIKER THEORY

APPENDIX C: CORRELATION MATRIX PROPAGATION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

APPENDIX A: FINDING THE FERMIONIC CORRECTION $\hat{P}$