# Non-Markovian Quantum Mpemba effect

David J. Strachan<EMAIL_ADDRESS>H. H. Wills Physics
Laboratory, University of Bristol, Bristol BS8 1TL, United Kingdom Archak
Purkayastha<EMAIL_ADDRESS>Department of Physics, Indian Institute
of Technology, Hyderabad 502284, India  Stephen R. Clark
<EMAIL_ADDRESS>H. H. Wills Physics Laboratory, University of
Bristol, Bristol BS8 1TL, United Kingdom

###### Abstract

Since it’s rediscovery in the twentieth century, the Mpemba effect, where a
far-from-equilibrium state may relax faster than a state closer to
equilibrium, has been extensively studied in classical systems and has
recently received significant attention in quantum systems. Many theories
explaining this counter-intuitive behavior in classical systems rely on memory
effects. However, in quantum systems, the relation between the Mpemba effect
and memory has remained unexplored. In this work, we consider a general non-
Markovian open quantum setting and reveal new classes of quantum Mpemba
effects, with no analog in Markovian quantum dynamics. Generically, open
quantum dynamics possess a finite memory time and a unique steady state. Due
to non-Markovian dynamics, even if the system is initialized in the steady
state it can take a long time to relax back. We find other initial states that
reach the steady state much faster. Most notably, we demonstrate that there
can be an initial state in which the system reaches the steady state within
the finite memory time itself, therefore giving the fastest possible
relaxation to stationarity. We verify the effect for quantum dot systems
coupled to electronic reservoirs in equilibrium and non-equilibrium setups at
weak, intermediate and strong coupling, and both with and without
interactions. Our work provides new insights into the rich physics underlying
accelerated relaxation in quantum systems.

In 1963, high school student Erasto B. Mpemba [1] rediscovered an intriguing
phenomenon, first proposed by Aristotle [2] and later discussed by Descartes
[3], where a hot liquid mixture freezes faster than an identical cold mixture.
The term Mpemba effect (MpE) has since been coined to describe the phenomenon
that a far-from-equilibrium state can relax to equilibrium faster than a state
closer to equilibrium. It has been studied extensively in many classical
systems [4, 5, 6, 7, 8, 9]. Many different theories exist in various contexts
explaining this behavior, suggesting that it is likely not one effect but
rather a broad umbrella for many different mechanisms of anomalous relaxation.
In classical systems, often this behavior is attributed to memory effects in
dynamics, although it has been shown to exist in Markovian classical systems
[10] also, where memory effects are absent. Recently, numerous theoretical and
experimental works have generalized MpEs to quantum systems [11, 12, 13, 14,
15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], considering either Markovian open
quantum dynamics or isolated systems at zero temperature. For both quantum and
classical systems, MpE has been characterised as strong when the state further
from the steady state relaxes exponentially faster and weak if it relaxes
faster but not exponentially. However, the interplay of memory and MpE in
quantum systems has remained unexplored. In this work, we investigate the
possibility of fast relaxation of specific initial states in general non-
Markovian open quantum systems, which have a finite, non-negligible memory
time and a unique steady state. This reveals new classes of quantum MpE, with
no analog in Markovian quantum dynamics.

We find that if such non-Markovian dynamics is initialized with the steady
state itself, it can take a long time to relax back, whereas there are other
special initial states which relax to the steady state much faster. We call
this the non-Markovian quantum Mpemba effect (NMQMpE). Most interestingly, in
addition to analogs of the strong and weak MpE, we also reveal the possibility
of an extreme case, where a specific initial state reaches the steady state
within the finite memory time itself, giving the fastest possible relaxation.
Our analysis is completely model-independent, with few physically motivated
assumptions. We verify our findings, focusing on the extreme NMQMpE, via
numerically exact calculations for single and double quantum dots coupled to
electronic reservoirs, as depicted in Fig. 1(a), for equilibrium and non-
equilibrium settings, at weak, intermediate and strong couplings to baths, and
both with and without many-body interactions in the system.

Figure 1: Setup and decomposition. (a) A system $S$ comprising a double
quantum dot setup with energies $\epsilon$, coupling $g$, interactions $U$ and
with each dot coupled to its own electronic thermal bath. (b) The
decomposition of the time evolution of the system density operator
$\hat{\rho}(\tau)$ for $\tau>\tau^{\mathcal{L}}_{\rm m}$ into the slippage
$\hat{\mathcal{S}}$ and time-independent propagator $\hat{\mathcal{L}}_{\rm
m}$.

General non-Markovian open quantum systems— A quantum system coupled to a bath
can in general be described by a Hamiltonian of the form
$\hat{H}=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{SB}$, where $\hat{H}_{S}$ is the
system Hamiltonian, $\hat{H}_{B}$ is the Hamiltonian of the environment and
$\hat{H}_{SB}$ describes the coupling between them. We consider quench
dynamics, assuming a factorised initial state of the form
$\hat{\rho}_{SB}(0)=\hat{\rho}(0)\hat{\rho}_{B}$, where $\hat{\rho}(0)$ is the
initial state of the system and $\hat{\rho}_{B}$ is the initial state of the
environment. After setting $k_{B}=1=\hbar$ and taking the Hamiltonian
$\hat{H}$ to be rescaled by some appropriate energy scale the time $\tau$
becomes dimensionless. The state of the system $\hat{\rho}(\tau)$ at later
times $\tau$ is obtained by the completely positive trace preserving (CPTP)
map, $\hat{\rho}(\tau)=\hat{\Lambda}(\tau)[\hat{\rho}(0)]$ where,
$\hat{\Lambda}(\tau)[\bullet]=\textrm{Tr}_{\textrm{B}}\big{(}e^{-i\hat{H}\tau}\bullet\hat{\rho}_{B}e^{i\hat{H}\tau}\big{)}.$
A related exact description of open quantum system dynamics is given by the
time convolutionless master equation [26, 27, 28]
$\frac{\partial\hat{\rho}(\tau)}{\partial\tau}=\hat{\mathcal{L}}(\tau)[\hat{\rho}(\tau)]$
with
$\hat{\mathcal{L}}(\tau)=\frac{d}{d\tau}[\hat{\Lambda}(\tau)]\hat{\Lambda}^{-1}(\tau)$
from which
$\hat{\Lambda}(\tau)=\mathcal{T}e^{\int_{0}^{\tau}\hat{\mathcal{L}}(\tau^{\prime})d\tau^{\prime}}$,
where $\mathcal{T}$ is the time-ordering operator. The effective propagator
$\hat{\mathcal{L}}(\tau)$ may be ill-defined because $\hat{\Lambda}(\tau)$ may
not be invertible. However, generically, at finite $\tau$, this happens only
at some isolated time points [29, 26, 30, 27]. Avoiding such times is enough
for the following analysis to hold. We consider the generic situation where
the system approaches a unique steady state in the long-time limit, i.e,
$\lim_{\tau\to\infty}\hat{\Lambda}(\tau)[\hat{\rho}(0)]=\hat{\rho}_{\textrm{SS}}$
for any initial state $\hat{\rho}(0)$. Without further approximations, we can
now show the possibility of NMQMpE.

Two different memory times— A crucial difference between the above general
non-Markovian description and Markovian quantum dynamics is that, both
$\hat{\Lambda}(\tau)$ and $\hat{\mathcal{L}}(\tau)$ are time-dependent,
despite the global Hamiltonian being time-independent. The instantaneous
steady state, or the time-dependent fixed point, can be defined either as the
eigenoperator of $\hat{\Lambda}(\tau)$ with eigenvalue $1$, or as the
eigenoperator of $\hat{\mathcal{L}}(\tau)$ with eigenvalue $0$. Neither time-
dependent fixed point may correspond to $\hat{\rho}_{\textrm{SS}}$ at short
times, but both converge to $\hat{\rho}_{\textrm{SS}}$ at long times. This
guarantees existence of two memory time scales $\tau^{\Lambda}_{\rm m}$,
$\tau^{\mathcal{L}}_{\rm m}$, such that

$\displaystyle||\hat{\Lambda}(\tau)[\hat{\rho}_{\textrm{SS}}]-\hat{\rho}_{\textrm{SS}}||<\epsilon,~{}~{}\forall~{}~{}\tau\geq\tau^{\Lambda}_{\rm
m},$ (1)
$\displaystyle||\hat{\mathcal{L}}(\tau)[\hat{\rho}_{\textrm{SS}}]||<\epsilon~{}~{}\forall~{}~{}\tau\geq\tau^{\mathcal{L}}_{\rm
m},$ (2)

where $\epsilon$ is some arbitrarily small tolerance, and $||\hat{A}||$ gives
the norm of $\hat{A}$. Equation (1) shows that, due to non-Markovianity, even
when the dynamics is initialized with $\hat{\rho}_{\rm SS}$, it takes a time
$\tau^{\Lambda}_{\rm m}$ to relax back. Equation (2) defines the time scale
$\tau^{\mathcal{L}}_{\rm m}$ in which the time-dependent fixed point of the
propagator effectively converges to $\hat{\rho}_{\rm SS}$. These two time
scales can be different in general, which then leads to NMQMpE as we discuss
below.

Unveiling the NMQMpE— Although not strictly required 111The NMQME doesn’t rely
on the convergence of $\hat{\mathcal{L}}(\tau)$, it only requires the
convergence of its fixed point
$\hat{\rho}_{\textrm{TDFP}}(\tau)\to\hat{\rho}_{\textrm{SS}}$ where
$\hat{\mathcal{L}}(\tau)[\hat{\rho}_{\textrm{TDFP}}(\tau)]=0$., it is possible
to argue on general grounds that $\hat{\mathcal{L}}(\tau)$ itself converges up
to an error $O(\epsilon)$ on the timescale $\tau^{\mathcal{L}}_{\rm m}$, i.e,
$\hat{\mathcal{L}}(\tau)=\hat{\mathcal{L}}_{m}+O(\epsilon),~{}~{}\forall~{}~{}\tau\geq\tau^{\mathcal{L}}_{\rm
m}$, where $\hat{\mathcal{L}}_{m}$ is the converged propagator. In this case
the following decomposition of the map applies

$\hat{\Lambda}(\tau)\approx e^{(\tau-\tau^{\mathcal{L}}_{\rm
m})\hat{\mathcal{L}}_{\rm
m}}\hat{\mathcal{S}}\quad\forall~{}~{}\tau\geq\tau^{\mathcal{L}}_{\rm m},$ (3)

where $\hat{\mathcal{S}}=\hat{\Lambda}(\tau^{\mathcal{L}}_{\rm m})$, as
depicted in Fig. 1(b). This decomposition represents a phenomenon called
initial slippage, which has been investigated in a wide class of systems,
usually with weak system-bath coupling [32, 33, 34, 35, 36]. It is closely
associated with the assumption of a finite memory time for the environment
correlations [37, 32]. From this decomposition, its clear
$\tau^{\mathcal{L}}_{\rm m}\leq\tau^{\Lambda}_{\rm m}$. Once the memory time
$\tau^{\mathcal{L}}_{\rm m}$ is reached the induced dynamics can be understood
in terms of the spectral decomposition of $\hat{\mathcal{L}}_{\rm m}$,
$\hat{\mathcal{L}}_{\rm
m}[\hat{\rho}]~{}=~{}\sum_{\mu}\lambda_{\mu}\hat{F}_{\mu}\textrm{Tr}(\hat{G}^{{\dagger}}_{\mu}\hat{\rho}),$
where $\hat{G}_{\mu}$, $\hat{F}_{\mu}$ define the damping basis for
$\hat{\mathcal{L}}_{\rm m}$ which satisfy $\hat{\mathcal{L}}_{\rm
m}\hat{F}_{\mu}=\lambda_{\mu}\hat{F}_{\mu}$,
$\hat{\mathcal{L}}^{{\dagger}}_{\rm
m}\hat{G}_{\mu}=\lambda^{*}_{\mu}\hat{G}_{\mu}$ and the normalisation
condition $\textrm{Tr}(\hat{G}_{\mu}\hat{F}_{\nu})=\delta_{\mu\nu}$. Due to
the complete positivity of $\hat{\mathcal{L}}_{\rm m}$,
$\textrm{Re}({\lambda_{\mu}})\leq 0$ with $\lambda_{1}=0$,
$\textrm{Tr}(\hat{F}_{\mu})=0$ for $\lambda_{\mu}\neq 0$ [38, 39, 40], where
the only eigenoperator corresponding to a physical state is the steady state
$\hat{F}_{1}=\hat{\rho}_{\textrm{SS}}$. We then have

$\displaystyle\hat{\rho}(\tau)\approx\hat{\rho}_{\textrm{SS}}+\sum_{\mu=2}^{d^{2}}e^{\lambda_{\mu}(\tau-\tau^{\mathcal{L}}_{\rm
m})}\textrm{Tr}(\hat{G}_{\mu}\hat{\mathcal{S}}[\hat{\rho}(0)]])\hat{F}_{\mu},$
(4)

where $d$ is the dimension of the system Hilbert space. The ordering of
eigenvalues is given by
$|\textrm{Re}(\lambda_{2})|\leq|\textrm{Re}(\lambda_{3})|\leq...\leq|\textrm{Re}(\lambda_{N})|$,
such that the timescale for the slowest relaxation is given by
$\tau_{\textrm{re}}=1/|\textrm{Re}(\lambda_{2})|$. The relaxation process is
then determined by the decay mode components of the state after the slippage,
$\hat{\alpha}_{\mu}(\rho_{0})=\textrm{Tr}(\hat{G}_{\mu}\hat{\mathcal{S}}[\hat{\rho}(0)])$.

In general, $[\hat{\mathcal{S}},\hat{\mathcal{L}}_{\rm m}]\neq 0$, meaning
$\hat{\mathcal{S}}$ will perturb $\hat{\rho}_{\textrm{SS}}$ itself, a property
unique to non-Markovian systems. This again guarantees
$\tau^{\mathcal{L}}_{\rm m}<\tau^{\Lambda}_{\rm m}$ and allows us to define
the NMQMpE as the existence of a state $\hat{\rho}_{\rm f}$ which relaxes
faster to the steady-state than $\hat{\rho}_{\textrm{SS}}$ itself. Since
NMQMpE is defined by a steady state it applies both in and out of equilibrium.

Weak, strong and extreme NMQMpE— There are three possible types of NMQMpE
which we call weak, strong and extreme. Let $\alpha_{\nu_{0}}(\hat{\rho}_{\rm
f})$ and $\alpha_{\kappa_{0}}(\hat{\rho}_{\textrm{SS}})$ be the components of
the slowest non-zero modes excited by $\hat{\mathcal{S}}$ when starting with
$\hat{\rho}_{\rm f}$ and $\hat{\rho}_{\textrm{SS}}$, respectively. The weak
NMQMpE arises when $\nu_{0}=\kappa_{0}$ and $\alpha_{\nu_{0}}(\hat{\rho}_{\rm
f})<\alpha_{\kappa_{0}}(\hat{\rho}_{\textrm{SS}})$. This means the slowest
non-zero modes excited by $\hat{\mathcal{S}}$ are same for both initial
conditions, but $\hat{\rho}_{\rm f}$ has a smaller amplitude in that mode.
Then, $\hat{\rho}_{\rm f}$ will relax faster, but not exponentially faster,
which is the hallmark of weak MpE.

If $\nu_{0}>\kappa_{0}$, there exists a strong NMQMpE where $\hat{\rho_{\rm
f}}$ relaxes exponentially faster than $\hat{\rho}_{\textrm{SS}}$ by a rate
given by the spectral gap
$\textrm{Re}|\lambda_{\nu_{0}}-\lambda_{\kappa_{0}}|$. The weak and strong
NMQMpE are analogs of weak and strong MpE observed in Markovian quantum
dynamics [14, 13, 15].

However, for NMQMpE, we can have an additional extreme case if
$\hat{\alpha}_{\nu}(\rho_{\rm f})=0$ for all $\nu$ in which
$\hat{\mathcal{S}}$ removes all decay components from $\rho_{\rm f}$ to give
the steady state $\hat{\rho}(\tau>\tau^{\mathcal{L}}_{\rm
m})=\hat{\mathcal{S}}[\hat{\rho}_{\rm f}]=\hat{\rho}_{\textrm{SS}}$. Thus, the
initial state converges to the steady state within the memory time
$\tau^{\mathcal{L}}_{\rm m}$. This is the fastest possible relaxation to the
steady state. Formally, this extremely fast relaxing state is given by

$\hat{\rho}_{\rm
f}=\frac{\hat{\mathcal{S}}^{-1}[\hat{\rho}_{\textrm{SS}}]}{\textrm{Tr}(\hat{\mathcal{S}}^{-1}[\hat{\rho}_{\textrm{SS}}])}.$
(5)

The above state may not be physically valid in general because (i)
$\mathcal{S}^{-1}$ may not exist and (ii) the resulting state may not be
positive semidefinite. However, in practice these conditions don’t pose an
issue. First, as mentioned before, at finite $\tau$, $\hat{\Lambda}(\tau)$ is
non-invertible only at isolated times. So we are free to shift
$\tau^{\mathcal{L}}_{\rm m}$ slightly to give a well defined
$\mathcal{S}^{-1}$ [41]. Second, if the system-bath coupling is not too
strong, $\mathcal{S}^{-1}$ can be found perturbatively [42, 43] and the
hermiticity of $\hat{\rho}_{\rm f}$ is guaranteed, as shown in the
supplementary material. In the following, we show the possibility of the
extreme NMQMpE in quantum dot setups and demonstrate that $\hat{\rho}_{\rm f}$
satisfying Eq. (5) is positive semidefinite over a wide range of parameters.

Baths for numerical examples— To explore the existence of $\hat{\rho}_{\rm
f}$, we consider a system $S$ coupled to two electronic thermal reservoirs,
allowing for both equilibrium and non-equilibrium dynamics, as shown in Fig.
1(a). We take the baths to be a continuum of modes, so that the total bath and
system-bath coupling Hamiltonians become
$\hat{H}_{B}+\hat{H}_{SB}=\sum_{\alpha=L,R}\big{(}\int_{-D}^{D}\omega\hat{c}_{\alpha}^{{\dagger}}(\omega)\hat{c}_{\alpha}(\omega)d\omega+\int_{-D}^{D}\sqrt{\mathcal{J}_{\alpha}(\omega)}(\hat{c}_{\alpha}^{{\dagger}}(\omega)\hat{Q}_{\alpha}+\hat{Q}^{{\dagger}}_{\alpha}\hat{c}_{\alpha}(\omega))d\omega\big{)},$
where $\mathcal{J}_{\alpha}(\omega)$ is the spectral density of the bath,
$\hat{Q}_{\alpha}$ is the system operator coupling to the bath and
$\hat{c}_{\alpha}^{{\dagger}}(\omega),\hat{c}_{\alpha}(\omega)$ are canonical
fermionic creation and annihilation operators for the bath modes obeying
$\\{\hat{c}_{\alpha}^{{\dagger}}(\omega),\hat{c}_{\alpha}(\omega^{\prime})\\}=\delta(\omega-\omega^{\prime})$,
with $\alpha=L,R$ corresponding to left and right baths. We parameterise
$\mathcal{J}_{\alpha}(\omega)$ via the total coupling strength
$\Gamma_{\alpha}$ defined as
$\Gamma_{\alpha}/D=\frac{1}{2}\int_{-D}^{D}2\pi\mathcal{J}_{\alpha}(\omega)d\omega$.
For simplicity, we assume an elliptical spectral function for both baths
$\mathcal{J}_{\alpha}(\omega)=(2\Gamma_{\alpha}/\pi^{2})\sqrt{1-(\omega/D)^{2}}$
222Note that the NMQME is not dependent on this choice, we have confirmed it
for a flat and Lorentzian band also.. The initial state of each bath is given
by a Gibbs state $\hat{\rho}_{B}=\prod_{\alpha}\hat{\rho}_{\alpha}$, where
$\hat{\rho}_{\alpha}=e^{-\beta_{\alpha}(\hat{H}_{\alpha}-\mu_{\alpha}\hat{N}_{\alpha})}/Z_{\alpha}$
is a thermal state and $Z_{\alpha}$ is the partition function for bath
$\alpha$.

Numerical methods— For our calculations the continuum of bath modes is
discretized using a chain mapping [45, 46] combined with the thermofield
transformation [47, 48] to encode the finite temperature effects exactly over
a finite time. In the absence of interactions, the Hamiltonian is quadratic so
can be solved exactly using unitary evolution of the single-particle
correlation matrix and compared to the analytical solution of the continuous
system via the Landauer-Büttiker formalism [49, 50, 51]. For interacting
systems, where $\hat{H}$ is no longer quadratic, we instead use the time
dependent variational principle [52, 53, 54] applied to matrix product states
[53] to efficiently simulate the real-time dynamics of the system and baths.
We extract $\hat{\mathcal{S}}$ from the evolved system state via the Choi-
Jamiolkowski isomorphism, as detailed in the supplementary material.

Quantum dot example— We first demonstrate our result for a quantum dot
described by a single fermionic mode $\hat{s}$ with an energy $\epsilon$,
giving $\hat{H}_{S}=\epsilon\hat{s}^{{\dagger}}\hat{s}$, and coupled via
$\hat{Q}_{\alpha}=\hat{s}$ to a single bath, i.e. one half of Fig. 1(a) with
$g=U=0$. In this case $\hat{\mathcal{S}}^{-1}$ is guaranteed to exist except
at the symmetric point where
$\hat{\rho}_{\textrm{SS}}=\frac{1}{2}(\ket{0}\bra{0}+\ket{1}\bra{1})$, where
$\ket{1}=\hat{s}^{{\dagger}}\ket{0}$.

Figure 2: NMQMpE for a single quantum dot. (a) The relaxation of population of
the single-dot towards steady state when starting from (i) the fastest
relaxing state, $\textrm{p}_{\rm f}(\tau)$, (ii) the steady state itself,
$\textrm{p}_{\rm SS}(\tau)$. (b) The population of the time-dependent fixed
point of $\hat{\rho}_{\rm TDFP}(\tau)$ satisfying
$\hat{\mathcal{L}}(\tau)[\hat{\rho}_{\rm TDFP}(\tau)]=0$, along with the time
dependent relaxation rate $\lambda_{2}(\tau)$. The convergence of both
quantities is consistent with $\tau^{\mathcal{L}}_{\rm m}\approx 20/D$. (c)
Phase diagram of $\delta\textrm{p}\equiv\textrm{p}_{\rm
f}-\textrm{p}_{\textrm{SS}}$, with contours of constant
$\textrm{p}_{\textrm{SS}}$ overlaid. The white represents regions with no
extreme NMQMpE. Parameters: $\epsilon=0,\beta_{\rm L}=\beta_{\rm
R}=10/D,\mu_{\rm L}=\mu_{\rm R}=0.1D,\Gamma_{\rm L}=\Gamma_{\rm R}=0.01D.$

For a single quantum dot, the only allowed states for the system are spanned
by a single parameter
$\hat{\rho}(\textrm{p})=(1-\textrm{p})\ket{0}\bra{0}+\textrm{p}\ket{1}\bra{1}$
with $0\leq\textrm{p}\leq 1$. The optimal state for a given
$\hat{\rho}_{\textrm{SS}}(\tau^{\mathcal{L}}_{\rm
m})=\hat{\rho}(\textrm{p}_{\textrm{SS}})$ extracted from
$\hat{\mathcal{L}}_{\rm m}$ is given by $\hat{\rho}_{\rm
f}=\hat{\rho}(\textrm{p}_{\rm f})$ where

$\textrm{p}_{\rm f}=\frac{\nu-(1-\textrm{p}_{\textrm{SS}})}{\nu-\sigma},$ (6)

with $\nu=\bra{0}\hat{\mathcal{S}}[\ket{0}\bra{0}]\ket{0}$,
$\sigma=\bra{0}\hat{\mathcal{S}}[\ket{1}\bra{1}]\ket{0}$, $0\leq\nu,\sigma\leq
1$. Eq. (6) is completely general for all bath and system parameters. But, it
corresponds to a physical solution for the fastest relaxing state only for the
parameter range where $0\leq\textrm{p}_{\rm f}\leq 1$. We also define the
population of the time-dependent fixed point as $\textrm{p}_{\rm TDFP}(\tau)$,
corresponding to the density operator $\hat{\rho}_{\rm TDFP}(\tau)$ satisfying
$\hat{\mathcal{L}}(\tau)[\hat{\rho}_{\rm TDFP}(\tau)]=0$, and define the
relaxation rate $\lambda_{2}(\tau)$ as the instantaneous non-zero eigenvalue
of $\hat{\mathcal{L}}(\tau)$.

As shown in Fig. 2(a), the fast initial state $\hat{\rho}_{\rm f}$ converges
within $\tau\approx 20/D$, which commensurate with the convergence of
$\hat{\mathcal{L}}(\tau)$, as seen in the plots of $\textrm{p}_{\rm TDFP}$ and
$\lambda_{2}(\tau)$ shown in Fig. 2(b). On the other hand in Fig. 2(a), the
steady state is seen to be significantly perturbed by the slippage and
subsequently decays back to itself on a much longer timescale. Fig. 2(c) shows
the deviation in density of $\hat{\rho}_{\rm f}$ from
$\hat{\rho}_{\textrm{SS}}$, serving as a measure of the strength of the NMQMpE
across the $(\mu,\beta)$-phase diagram of this system. Note that the anti-
symmetry about $\mu=0$ represents the particle-hole symmetry of the setup. The
effect disappears at the boundaries $\beta=0,\beta\to\infty$ due to the
following. At $\beta=0$, $[\hat{\mathcal{L}}_{\rm m},\hat{\mathcal{S}}]=0$
such that $\hat{\rho}_{\textrm{SS}}$ is invariant and there is no effect. For
$\beta\to\infty$, $p_{\textrm{SS}}=0,1$ for $\mu<0$, $\mu>1$ giving
$\textrm{p}_{\rm f}>1$, $\textrm{p}_{\rm f}<0$, respectively. A similar
reasoning can be used to explain the regions of high $|\mu|$, $\beta$ which
give no effect as $\textrm{p}_{\textrm{SS}}\approx 1$ for $\mu>0$ and
$\textrm{p}_{\textrm{SS}}\approx 0$ for $\mu<0$, again giving $\textrm{p}_{\rm
f}$ values outside the physically allowed range. The increase in the deviation
$\delta{\rm p}=\textrm{p}_{\rm f}-\textrm{p}_{\textrm{SS}}$ for increasing
$\beta$ can be seen as a measure of non-Markovianity of the dynamics. As
expected, lower temperatures create more non-Markovian behaviour leading to a
stronger NMQMpE. For a single quantum dot the effect is always an effective
equilibrium one, since any $\hat{\rho}$ can be expressed as a thermal state,
naturally connecting to the original effect.

Double quantum dot example— To explore a non-equilibrium setup we next
consider two quantum dots described by modes $\hat{s}_{1}$ and $\hat{s}_{2}$,
respectively. The system Hamiltonian is given by
$\hat{H}_{S}=\epsilon(\hat{s}^{{\dagger}}_{1}\hat{s}_{1}+\hat{s}^{{\dagger}}_{2}\hat{s}_{2})+g(\hat{s}^{{\dagger}}_{2}\hat{s}_{1}+\hat{s}^{{\dagger}}_{1}\hat{s}_{2})+U\hat{n}_{1}\hat{n}_{2},$
where $\hat{n}_{i}=\hat{s}^{{\dagger}}_{i}\hat{s}_{i}$. Each mode is connected
to its own bath via $\hat{Q}_{\rm L}=\hat{s}_{1}$ and $\hat{Q}_{\rm
R}=\hat{s}_{2}$. We look at a non-equilibrium setup parameterised as
$\delta\mu=(\mu_{\rm L}-\mu_{\rm R})/2$, $\bar{\mu}=(\mu_{L}+\mu_{R})/2$ and
$\beta_{L}=\beta_{R}=\beta$. To simplify the analysis we set
$\bar{\mu}=\epsilon=0$, and consider the dynamics of the particle current
between the two system modes, quantified by $\langle J(\tau)\rangle=i{\rm
Tr}[\hat{\rho}(\tau)(\hat{s}^{{\dagger}}_{1}\hat{s}_{2}-\hat{s}^{{\dagger}}_{2}\hat{s}_{1})]$.

Figure 3: NMQMpE for two quantum dots with $\Gamma_{\rm L}=\Gamma_{\rm
R}=\Gamma$ and $\epsilon=0,g=0.1D$. Panels (a), (b), (c) shows relaxation of
current between the two dots towards the non-equilibrium steady state value
when starting from (i) the fastest relaxing state, $\langle\hat{J}_{\rm
f}(\tau)\rangle$ and (ii) the steady state itself $\langle\hat{J}_{\rm
SS}(\tau)\rangle$. The steady state current is denoted $\langle\hat{J}_{\rm
SS}(\tau)\rangle$ and shown by the dashed line. In all cases
$\tau^{\mathcal{L}}_{\rm m}\approx 20/D$. Parameters in (a) $U=0$,
$\Gamma=0.02D$, (b) $U=0.3D$, $\Gamma=0.02D$, (c) $U=0$, $\Gamma=0.2D$, while
the thermal parameters for (a)-(c) are
$\beta_{L}=\beta_{R}=10D,\bar{\mu}=0,\delta\mu=0.1D$. (d) Phase diagram of
NMQMpE in terms of the chemical potential difference $\delta\mu$ against
inverse temperature $\beta_{L}=\beta_{R}=\beta$, $\Gamma=0.02D$, $U=0$. Here
the metric used to measure the deviation of $\hat{\rho}_{\rm f}$ from
$\hat{\rho}_{\textrm{SS}}$ is the trace distance $T(\hat{\rho}_{\rm
f},\hat{\rho}_{\textrm{SS}})$.

Figure 3(a) shows the NMQMpE for weak coupling and $U=0$. Here the fast state
$\hat{\rho}_{\rm f}$ possesses an initial current $\langle\hat{J}_{\rm
f}(0)\rangle$ very similar to the steady state current $\langle\hat{J}_{\rm
SS}\rangle$ of $\hat{\rho}_{\textrm{SS}}$. Yet the subsequent evolution shows
that $\langle\hat{J}_{\rm f}(\tau)\rangle$ increases slightly before
converging by $\tau\approx 20/D$ to $\langle\hat{J}_{\rm SS}\rangle$. In stark
contrast, the time-dependent current $\langle\hat{J}_{\rm SS}(\tau)\rangle$
obtained from using $\hat{\rho}_{\textrm{SS}}$ as the initial state displays a
damped large amplitude oscillation. Its angular frequency is $\approx 2g$,
corresponding to particles sloshing between the dots, and its relaxation back
to $\langle\hat{J}_{\rm SS}\rangle$ occurs with a decay time $\tau\approx
40/D$ consistent with $\tau_{\rm re}=1/|\rm Re(\lambda_{2})|$ extracted from
$\hat{\mathcal{L}}_{\rm m}$ for this problem. In Fig. 3(b) we consider the
same setup with $U=0.3D$. We see that all currents are reduced due to
interactions and $\langle\hat{J}_{\rm f}(0)\rangle$ now differs more
appreciably from $\langle\hat{J}_{\rm SS}\rangle$. Otherwise the NMQMpE is
quantitatively unchanged in the presence of interactions.

In Fig. 3(c) we consider a stronger coupling. In this case the slippage causes
a large perturbation resulting in the initial current of $\hat{\rho}_{\rm f}$
being reversed from that of $\hat{\rho}_{\textrm{SS}}$. However, the NMQMpE is
less pronounced in the subsequent time evolution because stronger coupling
also results in faster decay rates of the converged generator
$\hat{\mathcal{L}}_{\rm m}$ substantially reducing the difference in
timescales $\tau^{\mathcal{L}}_{\rm m}-\tau^{\Lambda}_{\rm m}$.

Concentrating on the weak coupling regime Fig. 3(d) shows the trace distance
$T(\hat{\rho}_{\rm
f},\hat{\rho}_{\textrm{SS}})\equiv\frac{1}{2}\textrm{Tr}[\sqrt{(\hat{\rho}_{\rm
f}-\hat{\rho}_{\textrm{SS}})^{{\dagger}}(\hat{\rho}_{\rm
f}-\hat{\rho}_{\textrm{SS}})}]$ for a range of $\delta\mu$ and non-zero
$\beta$ with $\bar{\mu}=0$. In contrast to Fig. 2(c) for a single dot, here
$\hat{\rho}_{\rm f}$ has a physical solution for the entire parameter regime
of the symmetric non-equilibrium double-dot setup shown. The reason for this
is that $\hat{\rho}_{\rm SS}$ in this region always lies near the centre of
the simplex of physical density operators [55]. So intuitively given
$\hat{\rho}_{\rm f}\propto\hat{\mathcal{S}}^{-1}[\hat{\rho}_{\rm SS}]$ is an
$O(\Gamma\tau^{\mathcal{L}}_{\rm m})$ perturbation of $\hat{\rho}_{\rm SS}$
then $\hat{\rho}_{\rm f}$ will be physical if it is more than
$O(\Gamma\tau^{\mathcal{L}}_{\rm m})$ away from a simplex boundary. Figure
3(d) displays increasing trace distance deviation $T(\hat{\rho}_{\rm
f},\hat{\rho}_{\textrm{SS}})$ with lower temperature consistent with the
dynamics being more non-Markovian. Much like the $\mu$ dependence for a single
dot in Fig. 2(c), the maximum deviation for the double-dot at a given $\beta$
occurs at an at optimal value of $\delta\mu$ which decreases with increasing
$\beta$. However, the increased state space of the double-dot yields very
different behaviour at $\delta\mu=0$ compared to a single dot at $\mu=0$. At
this point a single dot has $\hat{\rho}_{\rm f}=\hat{\rho}_{\textrm{SS}}$ and
no NMQMpE. For the double dot $\langle\hat{J}(\tau)\rangle=0$ and
$\langle\hat{s}_{1}^{{\dagger}}\hat{s}_{1}\rangle=\langle\hat{s}_{2}^{{\dagger}}\hat{s}_{2}\rangle=0.5$
for the time dependent fixed point $\hat{\mathcal{L}}(\tau)[\hat{\rho_{\rm
TDFP}}(\tau)]=0$ for all $\tau$. Yet despite $\hat{\rho}_{\rm f}$ and
$\hat{\rho}_{\textrm{SS}}$ having identical densities and currents at all
times the states can still differ in the real parts of their coherences
$\textrm{Re}(\langle\hat{s}^{{\dagger}}_{1}\hat{s}_{2}\rangle)$. A surprising
feature of Fig. 3(d) is that the greatest deviation occurs at equilibrium
$\delta\mu=0$ as $\bar{\beta}\to\infty$ where the relaxation dynamics is
entirely determined by the quantum coherence.

Possible experimental verification— The single dot setup lends itself to
experimental verification of the NMQMpE by measuring the dynamics of the dot’s
occupation using a quantum point contact charge sensor. Assuming a measurement
time resolution of $\approx 10~{}\mu{\rm s}$ [56] requires a memory time
$\tau^{\mathcal{L}}_{\rm m}\geq 20~{}\mu{\rm s}$ for the slippage to be
observable. Further assuming an occupation resolution of $1\%$ for the
measurement then requires $\delta\rm p\geq 0.01$ so that $\hat{\rho}_{\rm f}$
can be reliably distinguished from $\hat{\rho}_{\rm SS}$. This is satisfied
for the case shown in Fig. 2(a) for a narrow bandwidth bath $D=1~{}\mu{\rm
s}^{-1}$ implying a low temperature $T\approx 5~{}\mu{\rm K}$. In an realistic
setup the control of dot and bath couplings will not be a sudden quench. This
poses no fundamental issue for NMQMpE as finite time quenches can simply be
incorporated into the slippage $\hat{\mathcal{S}}$ and thus accounted for in
$\hat{\rho}_{\rm f}\propto\hat{\mathcal{S}}^{-1}[\hat{\rho}_{\textrm{SS}}]$,
as shown in the supplementary material.

Conclusion— We have shown that for a generic open system a novel Mpemba effect
can emerge within the slippage dynamics describing the memory time of the
bath. Specifically, we show that a fast initial state $\hat{\rho}_{\rm f}$ can
exist which relaxes to the steady state $\hat{\rho}_{\rm SS}$ within the bath
memory time. In contrast $\hat{\rho}_{\rm SS}$ initial state itself is
perturbed by the slippage and has yet to begin its slow relaxation back to
stationarity on the same timescale. We coin this an extreme Mpemba effect and
show its manifestation is widespread within quantum impurity setups both in
and out of equilibrium, both with and without interactions and for weak,
intermediate and strong coupling. The generality of our theoretical framework
lends itself to many different experimental platforms from quantum dots to
cold-atoms. We have shown that while the Mpemba effect isn’t exclusively
created by memory effects, it’s an inherent feature of them.

D.S and S.R.C. acknowledge insightful discussions with Francesco Turci, John
Goold and Tony Short. S.R.C. also gratefully acknowledges financial support
from UK’s Engineering and Physical Sciences Research Council (EPSRC) under
grant EP/T028424/1.

## References

  * Mpemba and Osborne [1969] E. B. Mpemba and D. G. Osborne, Cool?, Physics Education 4, 172 (1969).
  * W.D.Ross [1931] W.D.Ross, The Works of Aristotle, translated into English under the editorship of W.D. Ross (Oxford Clarendon Press, 1931).
  * Descartes [1950] R. Descartes,  _Discourse on Method_ (Harmondsworth, Penguin, Harmondsworth,, 1950).
  * Brownridge [2011] J. Brownridge, When does hot water freeze faster then cold water? a search for the mpemba effect, American Journal of Physics 79, 78 (2011).
  * Baity-Jesi  _et al._ [2019] M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvión, A. Gordillo-Guerrero, D. Iñiguez, A. Lasanta, A. Maiorano, E. Marinari, _et al._ , The mpemba effect in spin glasses is a persistent memory effect, Proceedings of the National Academy of Sciences 116, 15350 (2019).
  * Biswas  _et al._ [2020] A. Biswas, V. V. Prasad, O. Raz, and R. Rajesh, Mpemba effect in driven granular maxwell gases, Phys. Rev. E 102, 012906 (2020).
  * Lasanta  _et al._ [2017] A. Lasanta, F. Vega Reyes, A. Prados, and A. Santos, When the hotter cools more quickly: Mpemba effect in granular fluids, Phys. Rev. Lett. 119, 148001 (2017).
  * Greaney  _et al._ [2011] P. Greaney, G. Lani, G. Cicero, and J. Grossman, Mpemba-like behavior in carbon nanotube resonators, Metallurgical and Materials Transactions A 42, 3907 (2011).
  * Biswas  _et al._ [2023] A. Biswas, R. Rajesh, and A. K. Pal, Mpemba effect in a langevin system: Population statistics, metastability, and other exact results., The Journal of chemical physics 159 4 (2023).
  * Lu and Raz [2017] Z. Lu and O. Raz, Nonequilibrium thermodynamics of the markovian mpemba effect and its inverse, Proceedings of the National Academy of Sciences 114, 201701264 (2017).
  * Manikandan [2021] S. K. Manikandan, Equidistant quenches in few-level quantum systems, Phys. Rev. Res. 3, 043108 (2021).
  * Nava and Fabrizio [2019] A. Nava and M. Fabrizio, Lindblad dissipative dynamics in the presence of phase coexistence, Physical Review B 100, 10.1103/physrevb.100.125102 (2019).
  * Chatterjee  _et al._ [2023a] A. K. Chatterjee, S. Takada, and H. Hayakawa, Quantum mpemba effect in a quantum dot with reservoirs, Phys. Rev. Lett. 131, 080402 (2023a).
  * Carollo  _et al._ [2021] F. Carollo, A. Lasanta, and I. Lesanovsky, Exponentially accelerated approach to stationarity in markovian open quantum systems through the mpemba effect, Phys. Rev. Lett. 127, 060401 (2021).
  * Ivander  _et al._ [2023] F. Ivander, N. Anto-Sztrikacs, and D. Segal, Hyperacceleration of quantum thermalization dynamics by bypassing long-lived coherences: An analytical treatment, Phys. Rev. E 108, 014130 (2023).
  * Chatterjee  _et al._ [2023b] A. K. Chatterjee, S. Takada, and H. Hayakawa, Multiple quantum mpemba effect: exceptional points and oscillations (2023b), arXiv:2311.01347 [quant-ph] .
  * Zhang  _et al._ [2024] J. Zhang, G. Xia, C.-W. Wu, T. Chen, Q. Zhang, Y. Xie, W.-B. Su, W. Wu, C.-W. Qiu, P. xing Chen, W. Li, H. Jing, and Y.-L. Zhou, Observation of quantum strong mpemba effect (2024), arXiv:2401.15951 [quant-ph] .
  * Wang and Wang [2024] X. Wang and J. Wang, Mpemba effects in nonequilibrium open quantum systems (2024), arXiv:2401.14259 [quant-ph] .
  * Shapira  _et al._ [2024] S. A. Shapira, Y. Shapira, J. Markov, G. Teza, N. Akerman, O. Raz, and R. Ozeri, The mpemba effect demonstrated on a single trapped ion qubit (2024), arXiv:2401.05830 [quant-ph] .
  * Joshi  _et al._ [2024] L. K. Joshi, J. Franke, A. Rath, F. Ares, S. Murciano, F. Kranzl, R. Blatt, P. Zoller, B. Vermersch, P. Calabrese, C. F. Roos, and M. K. Joshi, Observing the quantum mpemba effect in quantum simulations (2024), arXiv:2401.04270 [quant-ph] .
  * Rylands  _et al._ [2023] C. Rylands, K. Klobas, F. Ares, P. Calabrese, S. Murciano, and B. Bertini, Microscopic origin of the quantum mpemba effect in integrable systems (2023), arXiv:2310.04419 [cond-mat.stat-mech] .
  * Bao and Hou [2022] R. Bao and Z. Hou, Accelerating relaxation in markovian open quantum systems through quantum reset processes (2022), arXiv:2212.11170 [quant-ph] .
  * Ares  _et al._ [2023] F. Ares, S. Murciano, and P. Calabrese, Entanglement asymmetry as a probe of symmetry breaking, Nature Communications 14, 10.1038/s41467-023-37747-8 (2023).
  * Kochsiek  _et al._ [2022] S. Kochsiek, F. Carollo, and I. Lesanovsky, Accelerating the approach of dissipative quantum spin systems towards stationarity through global spin rotations, Physical Review A 106, 10.1103/physreva.106.012207 (2022).
  * Caceffo  _et al._ [2024] F. Caceffo, S. Murciano, and V. Alba, Entangled multiplets, asymmetry, and quantum mpemba effect in dissipative systems (2024), arXiv:2402.02918 [cond-mat.stat-mech] .
  * Nestmann and Timm [2019] K. Nestmann and C. Timm, Time-convolutionless master equation: Perturbative expansions to arbitrary order and application to quantum dots (2019), arXiv:1903.05132 [cond-mat.mes-hall] .
  * Chruściński and Kossakowski [2010] D. Chruściński and A. Kossakowski, Non-markovian quantum dynamics: local versus nonlocal, Physical review letters 104, 070406 (2010).
  * Chruściński [2022] D. Chruściński, Dynamical maps beyond markovian regime, Physics Reports 992, 1 (2022).
  * Reimer  _et al._ [2019] V. Reimer, M. R. Wegewijs, K. Nestmann, and M. Pletyukhov, Five approaches to exact open-system dynamics: Complete positivity, divisibility, and time-dependent observables, J. Chem. Phys. 151, 044101 (2019).
  * Li  _et al._ [2012] Y. Li, B. Luo, J.-W. Zhang, and H. Guo, Time-convolutionless non-markovian master equation in strong-coupling regime (2012), arXiv:1202.3288 [quant-ph] .
  * Note [1] The NMQME doesn’t rely on the convergence of $\hat{\mathcal{L}}(\tau)$, it only requires the convergence of its fixed point $\hat{\rho}_{\textrm{TDFP}}(\tau)\to\hat{\rho}_{\textrm{SS}}$ where $\hat{\mathcal{L}}(\tau)[\hat{\rho}_{\textrm{TDFP}}(\tau)]=0$.
  * Bruch  _et al._ [2021] V. Bruch, K. Nestmann, J. Schulenborg, and M. Wegewijs, Fermionic duality: General symmetry of open systems with strong dissipation and memory, SciPost Physics 11, 10.21468/scipostphys.11.3.053 (2021).
  * Geigenmüller  _et al._ [1983] U. Geigenmüller, U. Titulaer, and B. Felderhof, Systematic elimination of fast variables in linear systems, Physica A: Statistical Mechanics and its Applications 119, 41 (1983).
  * Haake and Lewenstein [1983] F. Haake and M. Lewenstein, Adiabatic drag and initial slip in random processes, Phys. Rev. A 28, 3606 (1983).
  * Haake and Reibold [1985] F. Haake and R. Reibold, Strong damping and low-temperature anomalies for the harmonic oscillator, Phys. Rev. A 32, 2462 (1985).
  * Gaspard and Nagaoka [1999a] P. Gaspard and M. Nagaoka, Slippage of initial conditions for the redfield master equation, Journal of Chemical Physics 111, 5668 (1999a).
  * Nestmann  _et al._ [2021] 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).
  * Minganti  _et al._ [2019] F. Minganti, A. Miranowicz, R. W. Chhajlany, and F. Nori, Quantum exceptional points of non-hermitian hamiltonians and liouvillians: The effects of quantum jumps, Physical Review A 100, 10.1103/physreva.100.062131 (2019).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Rivas and Huelga [2012] ´. Rivas and S. F. Huelga,  _Open Quantum Systems_ (Springer Berlin Heidelberg, 2012).
  * Hegde  _et al._ [2021] 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).
  * Gaspard and Nagaoka [1999b] P. Gaspard and M. Nagaoka, Slippage of initial conditions for the redfield master equation, The Journal of chemical physics 111, 5668 (1999b).
  * Purkayastha  _et al._ [2016] A. Purkayastha, A. Dhar, and M. Kulkarni, Out-of-equilibrium open quantum systems: A comparison of approximate quantum master equation approaches with exact results, Physical Review A 93, 10.1103/physreva.93.062114 (2016).
  * Note [2] Note that the NMQME is not dependent on this choice, we have confirmed it for a flat and Lorentzian band also.
  * Chin  _et al._ [2011] A. W. Chin, S. F. Huelga, and M. B. Plenio, Chain representations of open quantum systems and their numerical simulation with time-adaptive density matrix renormalisation group methods, in  _Semiconductors and Semimetals_ , Vol. 85 (Elsevier, 2011) pp. 115–143.
  * Mühlpfordt  _et al._ [2020] T. Mühlpfordt, F. Zahn, V. Hagenmeyer, and T. Faulwasser, PolyChaos.jl – A Julia Package for Polynomial Chaos in Systems and Control, arXiv e-prints , arXiv:2004.03970 (2020), arXiv:2004.03970 [eess.SY] .
  * Takahashi and Umezawa [1996] Y. Takahashi and H. Umezawa, Thermo field dynamics, International journal of modern Physics B 10, 1755 (1996).
  * Borrelli and Gelin [2021] R. Borrelli and M. F. Gelin, Finite temperature quantum dynamics of complex systems: Integrating thermo-field theories and tensor-train methods, Wiley Interdisciplinary Reviews: Computational Molecular Science 11, e1539 (2021).
  * Ryndyk [2017] D. A. Ryndyk, Theory of quantum transport at nanoscale: An introduction dmitry a. ryndyk: Springer, 2016 246 pages, $129.00(e-book$99.00) isbn 978-3-319-24086-2, MRS Bulletin 42, 74–75 (2017).
  * Purkayastha [2019] A. Purkayastha, Classifying transport behavior via current fluctuations in open quantum systems, Journal of Statistical Mechanics: Theory and Experiment 2019, 043101 (2019).
  * Cornean  _et al._ [2005] H. D. Cornean, A. Jensen, and V. Moldoveanu, A rigorous proof of the landauer–büttiker formula, Journal of mathematical physics 46 (2005).
  * Yang and White [2020] M. Yang and S. R. White, Time-dependent variational principle with ancillary krylov subspace, Physical Review B 102, 10.1103/physrevb.102.094315 (2020).
  * Fishman  _et al._ [2022] M. Fishman, S. White, and E. Stoudenmire, The ITensor software library for tensor network calculations, SciPost Physics Codebases 10.21468/scipostphyscodeb.4 (2022).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Physical Review B 94, 10.1103/physrevb.94.165116 (2016).
  * Thirring  _et al._ [2011] W. Thirring, R. A. Bertlmann, P. Köhler, and H. Narnhofer, Entanglement or separability: the choice of how to factorize the algebra of a density matrix, The European Physical Journal D 64, 181–196 (2011).
  * Elzerman  _et al._ [2004] J. M. Elzerman, R. Hanson, L. H. Willems van Beveren, B. Witkamp, L. M. K. Vandersypen, and L. P. Kouwenhoven, Single-shot read-out of an individual electron spin in a quantum dot, Nature 430, 431–435 (2004).
  * Choi [1975] M.-D. Choi, Completely positive linear maps on complex matrices, Linear Algebra and its Applications 10, 285 (1975).
  * Jamiołkowski [1972] A. Jamiołkowski, Linear transformations which preserve trace and positive semidefiniteness of operators, Reports on Mathematical Physics 3, 275 (1972).
  * Amosov and Filippov [2016] G. G. Amosov and S. N. Filippov, Spectral properties of reduced fermionic density operators and parity superselection rule, Quantum Information Processing 16, 10.1007/s11128-016-1467-9 (2016).
  * Cheong and Henley [2004] S.-A. Cheong and C. L. Henley, Many-body density matrices for free fermions, Physical Review B 69, 075111 (2004).
  * Sánchez-Barquilla and Feist [2021] M. Sánchez-Barquilla and J. Feist, Accurate truncations of chain mapping models for open quantum systems, Nanomaterials 11, 10.3390/nano11082104 (2021).
  * Anto-Sztrikacs  _et al._ [2022] N. Anto-Sztrikacs, F. Ivander, and D. Segal, Quantum thermal transport beyond second order with the reaction coordinate mapping, The Journal of Chemical Physics 156, 10.1063/5.0091133 (2022).
  * Dorda  _et al._ [2014] 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).
  * Kohn and Santoro [2022] L. Kohn and G. E. Santoro, Quench dynamics of the anderson impurity model at finite temperature using matrix product states: entanglement and bath dynamics, Journal of Statistical Mechanics: Theory and Experiment 2022, 063102 (2022).
  * Chin  _et al._ [2010] A. W. Chin, Á . Rivas, S. F. Huelga, and M. B. Plenio, Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials, Journal of Mathematical Physics 51, 10.1063/1.3490188 (2010).

Supplementary Material for: Non-Markovian Quantum Mpemba effect
David J. Strachan,1 Archak Purkayastha,2 and Stephen R. Clark1

1H. H. Wills Physics Laboratory, University of Bristol, Bristol BS8 1TL,
United Kingdom
2Department of Physics, Indian Institute of Technology, Hyderabad 502284,
India
(Dated: )

##  I Redfield Equation

Here we show the slippage $\hat{\mathcal{S}}$ is always invertible in the weak
coupling regime. Consider a Hamiltonian
$\hat{H}=\hat{H}_{S}+\hat{H}_{B}+\gamma\hat{H}_{SB}$ with
$\hat{H}_{SB}=\sum_{\alpha}\hat{Q}_{\alpha}\hat{B}_{\alpha}^{{\dagger}}+{\rm
h.c.}$, $\hat{H}_{B}=\sum_{\alpha}\hat{H}_{B_{\alpha}}$ where the sum over
$\alpha$ includes all baths in contact with the system and
$||\hat{H}_{SB}||=O(1)$. The Redfield equation is then the second order
approximation to the Nakajima-Zwanzig master equation [43, 39],

$\displaystyle\frac{d\hat{\rho}}{d\tau}$
$\displaystyle=i[\hat{\rho},\hat{H}]-\gamma^{2}\sum_{\alpha}([\hat{Q}^{{\dagger}}_{\alpha},\hat{Q}^{1}_{\alpha}(\tau)\,\hat{\rho}]+[\hat{\rho}\,\hat{Q}^{2}_{\alpha}(\tau),\hat{Q}^{{\dagger}}_{\alpha}]+{\rm
h.c.}),$ $\displaystyle=\hat{\mathcal{L}}_{\rm RE}(\tau)\hat{\rho}(\tau),$
(S1)

with

$\displaystyle\hat{Q}^{1}_{\alpha}(\tau)=\int_{0}^{\tau}d\tau^{\prime}\hat{Q}^{I}_{\alpha}(-\tau^{\prime})\langle\hat{B}_{\alpha}(0)\hat{B}^{I^{{\dagger}}}_{\alpha}(-\tau^{\prime})\rangle_{B},\quad$
$\displaystyle\hat{Q}^{2}_{\alpha}(\tau)=\int_{0}^{\tau}d\tau^{\prime}\hat{Q}^{I}_{\alpha}(-\tau^{\prime})\langle\hat{B}^{I}_{\alpha}(-\tau^{\prime})\hat{B}^{{\dagger}}_{\alpha}(0)\rangle_{B},$
(S2)

where
$\hat{Q}^{I}_{\alpha}(\tau)=e^{i\hat{H}_{S}\tau}\hat{Q}_{\alpha}e^{-i\hat{H}_{S}\tau}$,
$\hat{B}^{I}_{\alpha}(\tau)=e^{i\hat{H}_{B_{\alpha}}\tau}\hat{B}_{\alpha}e^{-i\hat{H}_{B_{\alpha}}\tau}$
and the average is taken over bath $\alpha$. This description works well when
$\gamma$ is small compared to the other energy scales present. The memory
timescale is then given by the time such that

$\displaystyle|\langle\hat{B}_{\alpha}(0)\hat{B}^{I^{{\dagger}}}_{\alpha}(-\tau^{\prime})\rangle_{B}|<O(\epsilon)\quad\forall\tau>\tau^{\mathcal{L}}_{\rm
m},$
$\displaystyle|\langle\hat{B}^{I}_{\alpha}(-\tau^{\prime})\hat{B}^{{\dagger}}_{\alpha}(0)\rangle_{B}|<O(\epsilon)\quad\forall\tau>\tau^{\mathcal{L}}_{\rm
m},$ (S3)

where $\epsilon$ is an arbitrarily small error. This then gives a time
independent generator for $\tau>\tau_{\rm m}^{\mathcal{L}}$

$\displaystyle\frac{d\hat{\rho}}{d\tau}$
$\displaystyle=i[\hat{\rho},\hat{H}]-\gamma^{2}\sum_{\alpha}([\hat{Q}^{{\dagger}}_{\alpha},\hat{Q}^{1}_{\alpha}\,\hat{\rho}]+[\hat{\rho}\,\hat{Q}^{2}_{\alpha},\hat{Q}^{{\dagger}}_{\alpha}]+{\rm
h.c.}),$ (S4)

with

$\displaystyle\hat{Q}^{1}_{\alpha}=\int_{0}^{\infty}d\tau^{\prime}\hat{Q}^{I}_{\alpha}(-\tau^{\prime})\langle\hat{B}_{\alpha}(0)\hat{B}^{I^{{\dagger}}}_{\alpha}(-\tau^{\prime})\rangle_{B},\quad$
$\displaystyle\hat{Q}^{2}_{\alpha}=\int_{0}^{\infty}d\tau^{\prime}\hat{Q}^{I}_{\alpha}(-\tau^{\prime})\langle\hat{B}^{I}_{\alpha}(-\tau^{\prime})\hat{B}^{{\dagger}}_{\alpha}(0)\rangle_{B}.$
(S5)

We now calculate the fast state $\hat{\rho}_{\rm
f}=\hat{\mathcal{S}}^{-1}[\hat{\rho}_{\rm SS}$]. Solving Eq.(S4) up to
$O(\gamma^{2})$ gives

$\hat{\rho}(\tau)=e^{-i\hat{H}_{S}(\tau)}\hat{\rho}(0)e^{i\hat{H}_{S}(\tau)}-\gamma^{2}\sum_{\alpha}\int_{0}^{\tau}d\tau^{\prime}e^{-i\hat{H}_{S}(\tau-\tau^{\prime})}\bigg{(}[\hat{Q}^{{\dagger}}_{\alpha},\hat{Q}^{1}_{\alpha}(\tau^{\prime})\,\hat{\rho}(\tau^{\prime})]+[\hat{\rho}(\tau^{\prime})\,\hat{Q}^{2}_{\alpha}(\tau^{\prime}),\hat{Q}^{{\dagger}}_{\alpha}]+{\rm
h.c.}\bigg{)}e^{i\hat{H}_{S}(\tau-\tau^{\prime})},$ (S6)

and hence

$\hat{\rho}(0)=e^{i\hat{H}_{S}\tau}\hat{\rho}(\tau)e^{-i\hat{H}_{S}\tau}-\gamma^{2}\sum_{\alpha}\int_{0}^{\tau}d\tau^{\prime}e^{i\hat{H}_{S}\tau^{\prime}}\bigg{(}[\hat{Q}^{{\dagger}}_{\alpha},\hat{Q}^{1}_{\alpha}(\tau^{\prime})\,\hat{\rho}]+[\hat{\rho}\,\hat{Q}^{2}_{\alpha}(\tau^{\prime}),\hat{Q}^{{\dagger}}_{\alpha}]+{\rm
h.c.}\bigg{)}e^{-i\hat{H}_{S}\tau^{\prime}},$ (S7)

with
$\hat{\rho}(\tau^{\prime})=e^{-i\hat{H}_{S}(\tau^{\prime}-\tau)}\hat{\rho}(\tau)e^{i\hat{H}_{S}(\tau^{\prime}-\tau)}+O(\gamma^{2})$.
Setting $\tau=\tau^{\mathcal{L}}_{\rm m}$ and
$\hat{\rho}(\tau^{\mathcal{L}}_{\rm m})=\hat{\rho}_{\rm SS}$, we find the
$\hat{\rho}_{\rm f}$ up to $O(\gamma^{2})$

$\displaystyle\hat{\rho}_{\rm f}$
$\displaystyle=\hat{\rho}(0)=e^{i\hat{H}_{S}\tau}\hat{\rho}(\tau)e^{-i\hat{H}_{S}\tau}$
$\displaystyle-\gamma^{2}\sum_{\alpha}\int_{0}^{\tau}d\tau^{\prime}e^{i\hat{H}_{S}\tau^{\prime}}\bigg{(}[\hat{Q}^{{\dagger}}_{\alpha},\hat{Q}^{1}_{\alpha}(\tau^{\prime})e^{-i\hat{H}_{S}(\tau^{\prime}-\tau)}\hat{\rho}(\tau)e^{i\hat{H}_{S}(\tau^{\prime}-\tau)}]+[e^{-i\hat{H}_{S}(\tau^{\prime}-\tau)}\hat{\rho}(\tau)e^{i\hat{H}_{S}(\tau^{\prime}-\tau)}\hat{Q}^{2}_{\alpha}(\tau^{\prime}),\hat{Q}^{{\dagger}}_{\alpha}]+{\rm
h.c.}\bigg{)}e^{-i\hat{H}_{S}\tau^{\prime}}.$ (S8)

This formulation shows $\hat{\rho}_{\rm f}$ is always defined in the weak
coupling regime. This state may not be physical however, as the Redfield
description preserves hermicity and trace but not positivity. If this solution
gives non-negative eigenvalues for $\hat{\rho}_{\rm f}$, there exists a
NMQMpE.

##  II Calculation of $\hat{\Lambda}(\tau)$

To extract $\hat{\Lambda}(\tau)$, we use the Choi-Jamiolkowski isomorphism
[57, 58] which describes the correspondence between quantum maps and quantum
states. Consider maximally entangling the system with an auxillary system $A$
with the same Hilbert space $\mathcal{H}_{S}$ as $S$,
$\ket{\Phi^{+}}=\frac{1}{\sqrt{d}}\sum_{i=0}^{d-1}\ket{i}_{S}\otimes\ket{i}_{A}$.
Since $\hat{\Lambda}(\tau)$ is a completely positive trace preserving map,
$\hat{\rho}^{\Lambda}(\tau)=(\hat{\Lambda}(\tau)\otimes\mathds{1})\\{\ket{\Psi^{+}}\bra{\Psi^{+}}\\}$
is a nonnegative operator. Conversely, for any nonnegative operator on
$L(\mathcal{H}_{S})\otimes L(\mathcal{H}_{S})$, we can associate a quantum map
from operators on $L(\mathcal{H}_{S})$ to $L(\mathcal{H}_{S})$. This
isomorphism assumes a tensor product Hilbert space structure and can be
applied to 1D fermionic systems once they have been Jordan Wigner transformed
into effective spins. The fermionic map can then be extracted once the anti-
commuting behaviour of fermions is accounted for [59].

##  III Time evolution

For numerical calculations the continuum of modes for each bath $\alpha$ was
approximated using a finite number of modes $N_{\alpha}$. We use the same
number of modes for each bath $N_{\alpha}=N_{b}$. To do this we employ a
finite chain mapping using orthogonal polynomials and a thermofield
purification scheme to describe the finite temperature bath initial states.
This results in the following Hamiltonian,

$\displaystyle\hat{H}$
$\displaystyle=\sum_{j=1}^{N_{S}}\epsilon_{j}\hat{s}^{{\dagger}}_{j}\hat{s}_{j}+\sum_{j=1}^{N_{S}-1}g_{i}(\hat{s}^{{\dagger}}_{j+1}\hat{s}_{j}+\hat{s}^{{\dagger}}_{j}\hat{s}_{j+1})$
$\displaystyle+\sum_{\tilde{\alpha}}\bigg{[}\kappa_{\tilde{\alpha},0}^{-1}(\hat{s}^{{\dagger}}\hat{b}_{\tilde{\alpha},0}+\hat{b}^{{\dagger}}_{\tilde{\alpha},0}\hat{s})+\sum_{n=0}^{N_{b}}\bigg{(}\gamma_{\tilde{\alpha},n}\hat{b}^{{\dagger}}_{\tilde{\alpha},n}\hat{b}_{\tilde{\alpha},n}+\sqrt{\beta_{\tilde{\alpha},n+1}}\hat{b}^{{\dagger}}_{\tilde{\alpha},n}\hat{b}_{\tilde{\alpha},n+1}+\sqrt{\beta_{\tilde{\alpha},n+1}}\hat{b}^{{\dagger}}_{\tilde{\alpha},n+1}\hat{b}_{\tilde{\alpha},n}\bigg{)}\bigg{]}.$
(S9)

Here $\tilde{\alpha}=L_{f},L_{e},R_{f},R_{e}$ now denotes the four baths, left
and right correspondingly filled and empty, with all their parameters fully
defined by the bath spectral functions, temperatures and chemical potentials.
The 1D mode ordering used is given by

$\displaystyle\\{\hat{o}_{i}\\}=\\{\hat{b}_{L_{f}1},\hat{b}_{L_{e}1},...,\hat{b}_{L_{f}N_{b}},\hat{b}_{L_{e}N_{b}},\hat{s}_{1},...\hat{s}_{N_{s}},\hat{a}_{1},...,\hat{a}_{N_{s}},\hat{b}_{R_{f}1},\hat{b}_{R_{e}1},...,\hat{b}_{R_{f}N_{b}},\hat{b}_{R_{e}N_{b}}\\},$
(S10)

where $\hat{a}_{i}$ are the system ancilla modes. The details of this scheme
is outlined in Sec. VIII.

Given Eq. (S9) is non-interacting the exact dynamics can then be obtained via
unitary evolution of the single-particle correlation matrix
$\bm{C}_{ij}(t)=\textrm{Tr}(\hat{\rho}(t)\hat{o}^{{\dagger}}_{j}\hat{o}_{i})$
using the quadratic Hamiltonian defined via
$\hat{H}=\sum_{ij}\bm{h}_{ij}\hat{o}^{{\dagger}}_{i}\hat{o}_{j}$ and

$\bm{C}(t)=e^{i\bm{h}t}\bm{C}(t_{0})e^{-i\bm{h}t},$ (S11)

as shown in detail in Sec. IX. The reduced density matrix of the combined
system + ancilla setup is related to $C(t)$ via [60]

$\hat{\rho}^{\Lambda}(t)=\mathrm{det}(\mathbb{1}-\bm{C}^{\rm
T}(t))\mathrm{exp}\bigg{\\{}\sum_{ij\in SA}[\mathrm{log}(\bm{C}^{\rm
T}(t))(\mathbb{1}-\bm{C}^{\rm
T}(t))^{-1}]_{ij}\hat{o}^{{\dagger}}_{i}\hat{o}_{j}\bigg{\\}},$ (S12)

where $SA$ denotes the modes spanning the system + ancilla. This equation
holds provided $\hat{\rho}^{\Lambda}(t)$ is block diagonal in the number
basis. For this reason, we perform a particle-hole transformation on the
entangled states used in the Choi-Jamiolkowski isomorphism.

##  IV NMQMpE across coupling strength

Figure S1: NMQMpE for weak, intermediate and strong couplings for two non-
interacting modes in a non-equilibrium setup,
$\delta\mu=0.1D,\bar{\mu}=0,\beta_{L}=\beta_{R},\Gamma_{L}=\Gamma_{R}=\Gamma,U=0,\epsilon=0,g=0.1D$.

Here we show the NMQMpE for a number of parameter regimes, displaying the
effect of coupling and temperature for a non-interacting two mode setup. For
each coupling strength, the colder setup $\beta D=10$ always gives a larger
difference $T(\hat{\rho}_{\rm f},\hat{\rho}_{\textrm{SS}})$ than the hotter
case $\beta D=1$, reflecting the general intuition that colder systems display
more non-Markovian behaviour. For $\beta D=1$, as $\Gamma$ increases the
difference in convergence $\tau^{\Lambda}_{\rm m}-\tau^{\mathcal{L}}_{\rm m}$
decreases reflecting the increased decay rates of the converged generator
$\hat{\mathcal{L}}(\tau)$. In the case of strong coupling,
$\tau^{\Lambda}_{\rm m}-\tau^{\mathcal{L}}_{\rm m}\to 0$ as can be seen in
Figs. S1(d) and S1(h). For $\beta D=1$, the deviation $T(\hat{\rho}_{\rm
f},\hat{\rho}_{\textrm{SS}})$ increases for increasing $\Gamma$. This is
naively expected as $\Gamma$ controls the extent to which $\hat{\mathcal{S}}$
can deviate from the identity. However, this intuition breaks down at low
temperatures and strong couplings as seen in Fig. S1(d), which shows a smaller
deviation than Fig. S1(a)-(c).

##  V Finite Quench

Throughout this paper, we make the conventional assumption that both the bath
and system are disconnected until a sudden quench at $t=t_{0}$, i.e.
$\Gamma_{\alpha}(\tau)=\Gamma_{\alpha}\Theta(\tau)$ where $\Theta(\tau)$ is
the Heaviside function and $\alpha=L,R$. Realistically this will happen over a
finite timescale $\tau_{Q}$ such that
$\Gamma_{\alpha}(\tau)=\Gamma_{\alpha}f(\tau)$ for $\tau<\tau_{Q}$, where
$f(\tau)$ continuously ramps up from $0$ at $\tau=0$ to $\Gamma_{\alpha}$ at
$\tau=\tau_{Q}$. This poses no fundamental issues for NMQMpE as the existence
of a state $\hat{\rho}_{\rm f}\propto\hat{\mathcal{S}}^{-1}[\hat{\rho}_{\rm
SS}]$ doesn’t rely on either a sudden quench or time independent bath
couplings $\Gamma_{\alpha}$, and the finite quench can simply be incorporated
into the slippage $\hat{\mathcal{S}}$. The robustness of NMQME to a finite
$\tau_{Q}$ is shown in Fig. S2(b) with $f(\tau)=\Theta(\tau)\tau/\tau_{Q}$,
where its clear the effect is independent of $\tau_{Q}$, if taken into account
through $\hat{\Lambda}(\tau^{\mathcal{L}}_{\rm m})$. The existence of
$\hat{\rho}_{\rm f}$ may be invariant to $\tau_{Q}$ but its form will not,
i.e. $\hat{\rho}_{\rm f}$ = $\hat{\rho}_{\rm f}(\tau_{Q})$. To show this, we
look at the convergence of $\hat{\rho}_{f}(\tau_{Q}=0)$ for systems with
various ramp up times $\tau_{Q}$, as reported in Fig. S2(c). For small
$\tau_{Q}$, $\hat{\rho}_{\rm f}(\tau_{Q}=0)$ remains a fast converging state
but as $\tau_{Q}$ increases it deviates from the true $\hat{\rho}_{\rm f}$ and
thus doesn’t converge in the memory time $\tau^{\mathcal{L}}_{\rm m}$. Figure
S2(a) shows the convergence of the corresponding steady state $\hat{\rho}_{\rm
SS}$.

Figure S2: Simulating a finite quench for two modes using
$\tau^{\mathcal{L}}_{\rm m}=20D+\tau_{Q}$. (a) Convergence of
$\hat{\rho}_{\textrm{SS}}$. (b) Convergence of $\hat{\rho}_{\rm f}(\tau_{Q})$.
(c) Convergence of $\hat{\rho}_{\rm f}(\tau_{Q}=0)$. Parameters:
$\epsilon=0,g=0.1D,\beta_{L}=\beta_{R}=10/D,\mu_{L}=\mu_{R}=0.1D,\Gamma_{L}=\Gamma_{R}=0.03D,U=0$.

To do this, we only rescale the first bath modes of each of the four chains
$L_{f},L_{e},R_{f},R_{e}$ (see Sec. VIII) rather than re-calculating the chain
mapping at each time step in $0\leq\tau\leq\tau_{Q}$, as the rest of modes are
left invariant by a rescaling of $\Gamma$. The proof is as follows. We first
note that the chain mapping using orthogonal polynomials is equivalent to the
reaction coordinate (RC) chain mapping. Using the reaction coordinate
description for the bath chains in Eq.(S9) gives [61, 62]

$\displaystyle\beta_{\tilde{\alpha},n}=\int_{-D}^{D}J_{\tilde{\alpha},n}(\omega)d\omega,\quad$
$\displaystyle\gamma_{\tilde{\alpha},n}=\frac{1}{\beta_{\tilde{\alpha},n-1}}\int_{-D}^{D}J_{\tilde{\alpha},n-1}(\omega)d\omega,$
(S13)

where $\beta_{\tilde{\alpha},n}=\kappa_{\tilde{\alpha},0}^{-1}$ and
$J_{\tilde{\alpha}}(\omega)$ is the spectral function for $\tilde{\alpha}$.
The recursive spectral functions are given by

$J_{\tilde{\alpha},n+1}(\omega)=\frac{\beta_{\tilde{\alpha},n}J_{\tilde{\alpha},n}(\omega)}{\left|\pi
J_{\tilde{\alpha},n}(\omega)+\mathcal{P}\int_{-D}^{D}\frac{J_{\tilde{\alpha},n}(\omega^{\prime})d\omega^{\prime}}{\omega^{\prime}-\omega}\right|^{2}},$
(S14)

where $\mathcal{P}$ denotes Cauchy’s principle value. Changing the coupling
strength $\Gamma\to\Gamma^{\prime}$ is thus equivalent to
$J_{\tilde{\alpha},0}(\omega)\to(\Gamma^{\prime}/\Gamma)J_{\tilde{\alpha},0}(\omega)$.
This leaves $J_{\tilde{\alpha},n}(\omega)$ invariant for $n\geq 1$ as can be
seen in Eq. (S14) such that only the first site in the the chain is modified.

##  VI 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 [49, 50]

$J_{LB}^{P}=\frac{1}{2\pi}\int_{-D}^{D}d\omega\tau(\omega)[f_{L}(\omega)-f_{R}(\omega)],$
(S15)

and

$J_{LB}^{E}=\frac{1}{2\pi}\int_{-D}^{D}d\omega\omega\tau(\omega)[f_{L}(\omega)-f_{R}(\omega)],$
(S16)

where $f_{\alpha}(\omega)$ denotes the Fermi-Dirac distribution for bath
$\alpha$ and $\tau(\omega)$ is the transmission function of the system. This
can be calculated in terms of the non-equilibrium Green’s function [63, 50].
In our case this is given by

$G(\omega)=\bm{M}^{-1}(\omega),$ (S17)

with

$\bm{M}(\omega)=\omega\mathds{1}-\bm{h}_{S}-\bm{\Sigma}^{(1)}(\omega)-\bm{\Sigma}^{(N_{S})}(\omega),$
(S18)

where the only nonzero elements of the self-energy matrices of the leads
$\bm{\Sigma}^{(j)}(\omega)$ are

$[\bm{\Sigma}^{(j)}]_{jj}(\omega)=\mathcal{P}\int_{-D}^{D}d\omega^{\prime}\frac{\mathcal{J}_{j}(\omega^{\prime})}{\omega^{\prime}-\omega}-i\pi\mathcal{J}_{j}(\omega),$
(S19)

using $\mathcal{J}_{1}(\omega)=\mathcal{J}_{L}(\omega)$,
$\mathcal{J}_{N_{S}}(\omega)=\mathcal{J}_{R}(\omega)$ and $\bm{h}_{S}$ is
defined via
$\hat{H}_{S}=\sum_{ij}(\bm{h}_{S})_{ij}\hat{o}^{{\dagger}}_{i}\hat{o}_{j}$. If
the system Hamiltonian is of the form

$\bm{h}_{S}=\sum_{j=1}^{N_{S}}\epsilon_{j}\hat{s}^{{\dagger}}_{j}\hat{s}_{j}+\sum_{j=1}^{N_{S}-1}t_{i}(\hat{s}^{{\dagger}}_{j+1}\hat{s}_{j}+\text{h.c.}),$
(S20)

the transmission function is given by

$\tau(\omega)=4\pi^{2}\mathcal{J}_{L}(\omega)\mathcal{J}_{R}(\omega)|[G(\omega)]_{1D}|^{2}=\frac{\mathcal{J}_{L}(\omega)\mathcal{J}_{R}(\omega)}{|\text{det}[\bm{M}]|^{2}}\prod_{i=1}^{N_{S}-1}|t_{i}|^{2}.$
(S21)

##  VII Single quantum dot case

A single quantum dot is described by a spinless fermion mode as

$\hat{H}_{S}=\epsilon\hat{s}^{{\dagger}}\hat{s},$ (S22)

such that its density operator has to be block diagonal in the Fock basis
$\ket{0},\ket{1}$ as

$\hat{\rho}(\textrm{p})=(1-\textrm{p})\ket{0}\bra{0}+\textrm{p}\ket{1}\bra{1},$
(S23)

with an occupation $0\leq{\rm p}\leq 1$. Using $\textrm{Tr}(\hat{F}_{2})=0$
and $\textrm{Tr}(\hat{G}_{\mu}\hat{F}_{\nu})=\delta_{\mu\nu}$ for a single
mode we have

$\displaystyle\hat{F}_{2}=\ket{0}\bra{0}-\ket{1}\bra{1},\quad$
$\displaystyle\hat{F}_{1}=(1-\textrm{p}_{\rm
SS})\ket{0}\bra{0}+\textrm{p}_{\rm SS}\ket{1}\bra{1},$
$\displaystyle\hat{G}_{1}=\ket{0}\bra{0}+\ket{1}\bra{1},\quad$
$\displaystyle\hat{G}_{2}=\textrm{p}_{SS}\ket{0}\bra{0}+(\textrm{p}_{\rm
SS}-1)\ket{1}\bra{1},$

where $\hat{\rho}_{\rm SS}=\hat{\rho}(\textrm{p}_{\rm SS})$. Now consider the
slippage $\hat{\mathcal{S}}=\hat{\Lambda}(\tau^{\mathcal{L}}_{\rm m})$. After
vectorizing $\hat{\rho}$ we can express $\hat{\mathcal{S}}$ as a $2\times 2$
matrix acting in the physical subspace $\\{\ket{0}\bra{0},\ket{1}\bra{1}\\}$.
For this map to be CPTP, it must be of the form

$\hat{\mathcal{S}}=\begin{pmatrix}\nu&\sigma\\\ 1-\nu&1-\sigma\end{pmatrix}$
(S24)

with $\nu=\bra{0}\hat{\mathcal{S}}\\{\ket{0}\bra{0}\\}\ket{0}$,
$\sigma=\bra{0}\hat{\mathcal{S}}\\{\ket{1}\bra{1}\\}\ket{0}$, with
$0\leq\nu,\sigma\leq 1$. Now $\hat{\mathcal{S}}^{-1}$ is manifestly trace
preserving for the single mode,
$\textrm{Tr}(\hat{\mathcal{S}}^{-1}[\hat{\rho}_{\rm SS}])=1$. The fast state
is then given by $\hat{\rho}_{\rm f}=\hat{\rho}_{\rm f}(\textrm{p}_{\rm
f})=\hat{\mathcal{S}}^{-1}\hat{\rho}_{\rm SS}$ where

$\textrm{p}_{\rm f}=\frac{\nu-(1-\textrm{p}_{\rm SS})}{\nu-\sigma}.$ (S25)

If $\textrm{p}_{\rm f}\neq\textrm{p}_{\textrm{\rm SS}}$ and $0\leq{\rm p}_{\rm
f}\leq 1$ then we have a NMQMpE. The evolution of $\hat{\rho}(\tau)$ is given
by
$\frac{\partial\hat{\rho}(\tau)}{\partial\tau}=\mathcal{L}(\tau)\hat{\rho}(\tau)$
where $\mathcal{L}(\tau)=\sum_{\mu}\hat{F}_{\mu}\hat{G}^{{\dagger}}_{\mu}$ in
vectorised form, giving

$\frac{\partial\textrm{p}(\tau)}{\partial\tau}=\lambda_{2}(\tau)\big{(}\textrm{p}(\tau)-\textrm{p}_{\text{TDSS}}(\tau)\big{)}.$
(S26)

where $\textrm{p}_{\text{TDSS}}(\tau)$ is the occupation of the time-dependent
steady state of $\hat{\mathcal{L}}(\tau)$. This shows that the evolution is
controlled by the instantaneous steady state of $\hat{\mathcal{L}}(\tau)$,
acting as an attractor.

##  VIII Thermofield Transformation and chain mapping

The thermofield formalism is a commonly used technique in quantum field theory
[47, 48, 64] and statistical mechanics to relate the properties of a quantum
system at finite temperature (impure state) to a doubled system at zero
temperature (pure state). For a single fermionic bath mode $\hat{c}$ with
Hamiltonian $\hat{H}_{B}=\epsilon\hat{c}^{{\dagger}}\hat{c}$ it’s thermal
state at inverse temperature $\beta$ and chemical potential $\mu$ is

$\rho_{\beta}=\frac{1}{1+e^{-\beta(\epsilon-\mu)}}\big{(}\ket{0}\bra{0}+e^{-\beta(\epsilon-\mu)}\ket{1}\bra{1}\big{)},$
(S27)

where $\ket{1}=\hat{c}^{{\dagger}}\ket{0}$. Introducing an ancilla mode
$\hat{a}$ with Hamiltonian $\hat{H}_{A}=\epsilon\hat{a}^{{\dagger}}\hat{a}$,
we can express the bath’s thermal state as a partial trace of the thermofield
double state in the enlarged system as

$\ket{\Omega_{\beta}}=\sqrt{1-f}\ket{0}_{B}\otimes\ket{0}_{A}+\sqrt{f}\ket{1}_{B}\otimes\ket{1}_{A},$
(S28)

where $f=(1+e^{\beta(\epsilon-\mu)})^{-1}$ is the fermi factor. Performing a
particle hole transformation on the ancilla mode, we have

$\displaystyle\ket{\Omega_{\beta}}$
$\displaystyle=\sqrt{1-f}\ket{0}_{B}\otimes\ket{1}_{A}+\sqrt{f}\ket{1}_{B}\otimes\ket{0}_{A},$
$\displaystyle=\big{(}\sqrt{1-f}\hat{a}^{{\dagger}}+\sqrt{f}\hat{c}^{{\dagger}}\big{)}\ket{\text{vac}}.$

Defining two new fermionic mode operators,

$\displaystyle\hat{f}^{{\dagger}}=\sqrt{1-f}\hat{a}^{{\dagger}}+\sqrt{f}\hat{c}^{{\dagger}},\quad$
$\displaystyle\hat{e}^{{\dagger}}=\sqrt{f}\hat{a}^{{\dagger}}-\sqrt{1-f}\hat{c}^{{\dagger}},$

we see that the thermal state can then be expressed as a single particle
product state $\ket{\Omega_{\beta}}=\hat{f}^{{\dagger}}\ket{\text{vac}}$. We
now transform our Hamiltonian into this basis. The self energy terms take a
simple form,

$\hat{{H}}_{A}+\hat{H}_{B}=\epsilon\big{(}\hat{a}^{{\dagger}}\hat{a}+\hat{c}^{{\dagger}}\hat{c}\big{)}=\epsilon\big{(}\hat{f}^{{\dagger}}\hat{f}+\hat{e}^{{\dagger}}\hat{e}\big{)},$
(S29)

and if we assume the mode couples to the system via a hybridisation term
$\hat{H}_{SB}$ with one system mode $\hat{s}$ we then have

$\displaystyle\hat{H}_{SB}=v\hat{b}^{{\dagger}}\hat{s}+v^{*}\hat{s}^{{\dagger}}\hat{b}=v\sqrt{f}\hat{f}^{{\dagger}}\hat{s}-v\sqrt{1-f}\hat{e}^{{\dagger}}\hat{s}+v^{*}\sqrt{f}\hat{s}^{{\dagger}}\hat{f}-v^{*}\sqrt{1-f}\hat{s}^{{\dagger}}\hat{e}.$
(S30)

Moving back to the continuum, $v\to\sqrt{\mathcal{J}_{\alpha}(\omega)}$,$f\to
f_{\alpha}(\omega)$, $\hat{c}\to\hat{c}_{\alpha}(\omega)$, we have

$\hat{H}=\hat{H}_{S}+\sum_{\alpha=L,R}\int_{-D}^{D}\omega\big{(}\hat{f}_{\alpha}^{{\dagger}}(\omega)\hat{f}_{\alpha}(\omega)+\hat{e}_{\alpha}^{{\dagger}}(\omega)\hat{e}_{\alpha}(\omega)\big{)}+\sqrt{\mathcal{J}_{\alpha
f}(\omega)}\big{[}\hat{f}_{\alpha}^{{\dagger}}(\omega)\hat{s}+\hat{s}^{{\dagger}}\hat{f}_{\alpha}(\omega)\big{]}-\sqrt{\mathcal{J}_{\alpha
e}(\omega)}\big{[}\hat{e}_{\alpha}^{{\dagger}}(\omega)\hat{s}+\hat{s}^{{\dagger}}\hat{e}_{\alpha}(\omega)\big{]}\mathrm{d}\omega,$
(S31)

with

$\displaystyle\mathcal{J}_{\alpha
f}(\omega)=f_{\alpha}(\omega)\mathcal{J}_{\alpha}(\omega),\quad$
$\displaystyle\mathcal{J}_{\alpha
e}(\omega)=\big{(}1-f_{\alpha}(\omega)\big{)}\mathcal{J}_{\alpha}(\omega).$

This mapping has moved the dependence on temperature from the state into the
Hamiltonian, where two separate baths are coupled to the system, one filled
and one empty. This separation between the filled and empty modes gives us
much greater freedom in optimising our Hamiltonian for matrix product state
calculations. The thermal state $\hat{\rho}_{B}$ has become a purified
thermofield state
$\ket{\Omega_{\beta}}=\prod_{k}\ket{\Omega_{\beta,k}}=\prod_{k}\hat{f}^{{\dagger}}_{k}\ket{\textrm{vac}}$.

Once the thermofield transformation is applied, we then map each of the
continuous baths to finite chains using orthogonal polynomials [65, 61]. This
is done numerically using the orthopol package [46] which implements the
following protocol. For simplicity we only consider one bath, but the analysis
generalises straightforwardly. We have the following bath Hamiltonian terms

$\hat{H}_{SB}=\sum_{c=1,2}\int_{-D}^{D}d\omega\sqrt{J_{c}(\omega)}(\hat{s}^{{\dagger}}\hat{f}_{c}(\omega)+{\rm
h.c.}),$ (S32)
$\hat{H}_{A}+\hat{H}_{B}=\sum_{c=1,2}\int_{-D}^{D}d\omega\omega\hat{f}_{c}^{{\dagger}}(\omega)\hat{f}_{c}(\omega),$
(S33)

where we have denoted empty modes as $\hat{e}(\omega)=\hat{f}_{1}(\omega)$ and
the filled modes as $\hat{f}(\omega)=\hat{f}_{2}(\omega)$. To carry out the
chain mapping, we define new fermionic operators

$\hat{b}_{c,n}=\int_{-D}^{D}d\omega\sqrt{J_{c}(\omega)}\kappa_{c,n}\pi_{c,n}(\omega)\hat{f}_{c}(\omega),$
(S34)

with an inverse transformation

$\hat{f}_{c}(\omega)=\sum_{n=0}^{\infty}\sqrt{J_{c}(\omega)}\kappa_{c,n}\pi_{c,n}(\omega)\hat{b}_{c,n},$
(S35)

where $\pi_{c,n}(x)$ is an $n$th monic polynomial with a corresponding
normalisation constant $\kappa_{c,n}$ (defined below). We then have
$\pi_{c,n}(x)=\sum_{j=0}^{n}c_{nj}x^{j}$, where the monic condition means
$c_{nn}=1$. We define $\pi_{c,n}(x)$ such that they obey the following
orthogonality condition

$\int_{-D}^{D}d\omega
J_{c}(\omega)\pi_{c,m}(\omega)\pi_{c,n}(\omega)=\kappa_{c,n}^{-2}\delta_{n,m},$
(S36)

which also defines the normalisation constants $\kappa_{c,n}$. Note that this
transformation will leave the state invariant as the new creation
(annihilation) operators are linear combinations of creation (annihilation)
operators only, so a filled (empty) bath state is mapped to a filled (empty)
chain. This is major advantages of the thermofield method for tackling finite
temperature. This gives

$\displaystyle\hat{H}_{SB}$
$\displaystyle=\sum_{c}\sum_{n=0}^{\infty}\kappa_{c,n}(\hat{s}^{{\dagger}}\hat{b}_{c,n}+{\rm
h.c.})\int_{-D}^{D}d\omega J_{c}(\omega)\pi_{c,n}(\omega)$
$\displaystyle=\sum_{n=0}^{\infty}\kappa_{c,n}(\hat{s}^{{\dagger}}\hat{b}_{c,n}+{\rm
h.c.})\int_{-D}^{D}d\omega J_{c}(\omega)\pi_{c,n}(\omega)\pi_{c,0}(\omega)$
$\displaystyle=\sum_{c}\kappa_{c,0}^{-1}\hat{s}^{{\dagger}}\hat{b}_{c,0}+{\rm
h.c.}.$ (S37)

Now consider the bath Hamiltonian

$\displaystyle\hat{H}_{A}+\hat{H}_{B}=\sum_{c=1,2}\sum_{n,m=0}^{\infty}\hat{b}^{{\dagger}}_{c,n}\hat{b}_{c,m}\int_{-D}^{D}d\omega\omega\pi_{c,n}(\omega)\pi_{c,m}(\omega).$
(S38)

To progress, we make use of the following recurrence relation for the monic
polynomials.

$\pi_{c,n+1}(\omega)=(x-\gamma_{c,n})\pi_{c,n}(\omega)-\beta_{c,n}\pi_{c,n-1}(\omega),$
(S39)

where $\gamma_{c,n}$ and $\beta_{c,n}$ are uniquely determined by the weight
function as

$\gamma_{c,n}=\kappa_{c,n}^{2}\int_{-D}^{D}d\omega\omega
J_{c}(\omega)\pi^{2}_{c,n}(\omega),$ (S40)

and

$\beta_{c,n}=\kappa_{c,n}\kappa_{c,n+1}\int_{-D}^{D}d\omega
J_{c}(\omega)\omega\pi_{c,n}(\omega)\pi_{c,n-1}(\omega).$ (S41)

Using this, we have

$\displaystyle\hat{H}_{A}+\hat{H}_{B}$
$\displaystyle=\sum_{c=1,2}\sum_{n,m=0}^{\infty}\kappa_{c,n}^{2}\hat{b}^{{\dagger}}_{c,n}\hat{b}_{c,m}\int_{-D}^{D}d\omega
J_{c}(\omega)\pi_{c,m}(\omega)\big{[}\pi_{c,n+1}(\omega)+\beta_{c,n}\pi_{c,n-1}(\omega)+\gamma_{c,n}\pi_{c,n}(\omega)\big{]},$
$\displaystyle=\sum_{c=1,2}\sum_{n=0}^{\infty}\bigg{(}\gamma_{c,n}\hat{b}^{{\dagger}}_{c,n}\hat{b}_{c,n}+\frac{\kappa_{c,n+1}\beta_{c,n+1}}{\kappa_{c,n}}\hat{b}^{{\dagger}}_{c,n}\hat{b}_{c,n+1}+\frac{\kappa_{c,n}}{\kappa_{c,n+1}}\hat{b}^{{\dagger}}_{c,n+1}\hat{b}_{c,n}\bigg{)},$
$\displaystyle=\sum_{c=1,2}\sum_{n=0}^{\infty}\bigg{(}\gamma_{c,n}\hat{b}^{{\dagger}}_{c,n}\hat{b}_{c,n}+\sqrt{\beta_{c,n+1}}\hat{b}^{{\dagger}}_{c,n}\hat{b}_{c,n+1}+\sqrt{\beta_{c,n+1}}\hat{b}^{{\dagger}}_{c,n+1}\hat{b}_{c,n}\bigg{)}.$
(S42)

Combining this with the rest of the Hamiltonian and choosing a finite cutoff
$N_{b}$ for both baths, we obtain

$\displaystyle\hat{H}$
$\displaystyle=\hat{H}_{S}+\sum_{\tilde{\alpha}}\kappa_{\tilde{\alpha},0}^{-1}(\hat{s}^{{\dagger}}\hat{b}_{\tilde{\alpha},0}+\hat{b}^{{\dagger}}_{\tilde{\alpha},0}\hat{s})$
$\displaystyle\qquad\qquad+\sum_{\tilde{\alpha}}\sum_{n=0}^{N_{n}}\bigg{(}\gamma_{\tilde{\alpha},n}\hat{b}^{{\dagger}}_{\tilde{\alpha},n}\hat{b}_{\tilde{\alpha},n}+\sqrt{\beta_{\tilde{\alpha},n+1}}\hat{b}^{{\dagger}}_{\tilde{\alpha},n}\hat{b}_{\tilde{\alpha},n+1}+\sqrt{\beta_{\tilde{\alpha},n+1}}\hat{b}^{{\dagger}}_{\tilde{\alpha},n+1}\hat{b}_{\tilde{\alpha},n}\bigg{)},$
(S43)

where $\tilde{\alpha}$ runs over all combinations of $\alpha,c$.

##  IX Correlation matrix propagation

Here we prove Eq. (S11). Defining $\mathcal{U}(\tau)=e^{-i\hat{H}\tau}$, we
have

$\displaystyle\bm{C}_{ij}(\tau)$
$\displaystyle=\bra{\psi(\tau)}\hat{c}^{{\dagger}}_{j}\hat{c}_{i}\ket{\psi(\tau)}=\bra{\psi(0)}\mathcal{U}^{{\dagger}}(\tau)\hat{c}^{{\dagger}}_{j}\hat{c}_{i}\mathcal{U}(\tau)\ket{\psi(0)}.$

To evaluate this, first consider how $\hat{H}$ acts on
$\hat{c}^{{\dagger}}_{j}$,

$\displaystyle\hat{H}\hat{c}^{{\dagger}}_{j}$
$\displaystyle=\sum_{kl}\bm{H}_{kl}\hat{c}^{{\dagger}}_{k}\hat{c}_{l}\hat{c}^{{\dagger}}_{j}=\sum_{kl}\bm{H}_{kl}\hat{c}^{{\dagger}}_{k}(\delta_{lj}-\hat{c}^{{\dagger}}_{j}\hat{c}_{l})=\sum_{k}\bm{H}_{kj}\hat{c}^{{\dagger}}_{k}.$
(S44)

We can now evaluate $\mathcal{U}^{{\dagger}}\hat{c}^{{\dagger}}_{j}$ as

$\displaystyle\mathcal{U}^{{\dagger}}\hat{c}^{{\dagger}}_{j}$
$\displaystyle=e^{iHt}\hat{c}^{{\dagger}}_{j}=e^{i\sum_{k}\bm{H}_{kj}t}\hat{c}^{{\dagger}}_{k}=\sum_{k}\bm{U}^{{\dagger}}_{kj}\hat{c}^{{\dagger}}_{k},$
(S45)

Substituting this into Eq. (IX) gives the result

$\displaystyle\bm{C}_{ij}(t)$
$\displaystyle=\bra{\psi(0)}\sum_{kl}\bm{U}^{\dagger}_{kj}\hat{c}^{\dagger}_{k}(\bm{U}^{\dagger}_{li}\hat{c}^{\dagger}_{l})^{\dagger}\ket{\psi(0)}=\sum_{kl}\bm{U}^{\dagger}_{kj}\bm{U}_{il}\bra{\psi(0)}\hat{c}^{\dagger}_{k}\hat{c}_{l}\ket{\psi(0)}=\sum_{kl}\bm{U}^{\dagger}_{kj}\bm{U}_{il}\bm{C}_{lk}(0).$
(S46)

Thus we have

$\bm{C}(t)=\bm{U}\bm{C}(0)\bm{U}^{{\dagger}}.$ (S47)