# Steady state in strong system-bath coupling regime: reaction coordinate
versus perturbative expansion

Camille L. Latune1,2 1Univ Lyon, ENS de Lyon, CNRS, Laboratoire de
Physique,F-69342 Lyon, France
2Quantum Research Group, School of Chemistry and Physics, University of
KwaZulu-Natal, Durban, KwaZulu-Natal, 4001, South Africa

###### Abstract

Motivated by the growing importance of strong system-bath coupling in several
branches of quantum information and related technological applications, we
analyze and compare two strategies currently used to obtain (approximately)
steady states in strong coupling regime. The first strategy is based on
perturbative expansions while the second one uses reaction coordinate mapping.
Focusing on the widely used spin-boson model, we show that the predictions of
these two strategies coincide in many situations. This confirms and
strengthens the relevance of both techniques. Beyond that, it is also crucial
to know precisely their respective range of validity. In that perspective,
thanks to their different limitations, we use one to benchmark the other. We
introduce and successfully test some very simple validity criteria for both
strategies, bringing some answers to the question of the validity range.

##  I Introduction

It is notoriously challenging to describe the dynamics and the steady states
of quantum systems coupled to noisy environment[1]. This is particularly true
when the coupling strength is such that system-environment correlations cannot
be neglected, invalidating the traditional Born and Markov approximations [2],
characterizing the so-called strong coupling regime.

Still, strong coupling effects are playing increasing role in quantum
transport [3, 4, 5, 6, 7, 8, 9, 10, 11], quantum sensing [12, 13, 14], quantum
thermal engines [15, 16, 17, 18, 19, 20, 21], magnets properties for memory
hard-drive [22], and possibly in biological systems [23, 24, 25, 26]. For such
applications, since one is usually interested in stationary properties and
performances, the knowledge of the steady state of the strongly dissiaptive
dynamics is enough.

Several techniques have been established to gain access to the strong coupling
dynamics where usual Markovian Master equations cannot be used
straightforwardly. One of them is the polaron transformation [27, 28, 29, 30,
15, 31], which consists in analyzing the problem in a rotating frame with
respect to the coupling Hamiltonian. In the following, however, we will focus
on two other approaches which have been recently used to obtain estimates of
the steady states in the strong coupling regime. The first one uses embedding
techniques like reaction coordinate [32, 33, 16, 5, 6, 19] and pseudo-mode
[34, 35, 36, 37] to obtain the dynamics of the system of interest and then
considers time going to infinity. The second one relies on a result sustained
by several studies [38, 39, 40, 43, 41, 42] establishing, under some generic
conditions, that a system $S$ interacting strongly with a thermal bath $B$
tends together with $B$ to a global thermal state. Greater details on this
technique can be found in the recent review [44].

Both approaches have some strengths and weaknesses. The limitations of the
first approach relying on embedding techniques typically come from the bath
spectral density. On the other hand, the results are expected to have a broad
range of validity in terms of coupling strength. The second approach requires
to trace out the bath in the global thermal state, which actually amounts to
similar difficulties as computing the exact dynamics in the first place. Thus,
one is left with perturbative expansions, with limited range of validity.
However, within the range of validity of these expansions, the obtained
expressions are expected to provide very good description of the steady
states.

Even though one expects these two approaches to coincide, at least for some
regions of parameter, this has not been tested. This is the first aim of this
paper. Secondly, we will use the strength of each approach to benchmark and
define more precisely the range of validity of the other approach. This allows
us to introduce and successfully test some simple validity criteria for both
approaches, providing some answers to the crucial question of validity range
for each approach.

##  II Perturbative expansion approach

We consider a system $S$, of self Hamiltonian $H_{S}$, interacting with a
thermal bosonic bath $B$ of self Hamiltonian $H_{B}$ at inverse temperature
$\beta:=1/k_{B}T$ ($T$ being the usual temperature). The interaction is of the
form $V=AB$, where $A$ is an observable of $S$, and $B$ is the standard
bosonic operator $B=\sum_{k}g_{k}(b_{k}^{{\dagger}}+b_{k})$, where $g_{k}$ is
the coupling coefficient between $S$ and the $k^{\rm th}$ mode of the bath
(setting $\hbar=1$), with creation and annihilation operators
$b_{k}^{{\dagger}}$ and $b_{k}$, respectively. The starting point of the
perturbative expansion approach is the convergence of the system together with
the bath towards the global thermal state

$\rho_{SB}^{\rm th}=Z_{SB}^{-1}e^{-\beta H_{SB}},$ (1)

where $H_{SB}:=H_{S}+H_{B}+V$ is the total Hamiltonian generating the dynamics
of $SB$, and $Z_{SB}:={\rm Tr}_{SB}[e^{-\beta H_{SB}}]$ is the partition
function. This fundamental result has been widely used for classical and
quantum systems, and is supported by several important studies focusing on
quantum systems [38, 39, 40, 43, 41, 42, 44] under the important condition
$[H_{S},V]\neq 0$ (which can be seen as pure dephasing and can be simply
solved). From (1), the reduced steady state of $S$ is obtained by tracing out
$B$,

$\rho_{S}^{\rm ss}:={\rm Tr}_{B}[\rho_{SB}^{\rm th}].$ (2)

As usually, tracing out the bath is a very challenging task which can only be
done approximately. A common approach is perturbative expansion, whose main
steps are presented below (for more details see Appendix A and [43, 44, 45,
47, 46]). By “taking out” the local contributions in (2) and then expanding up
to second order, we obtain,

$\displaystyle\rho_{S}^{\rm ss}$ $\displaystyle=$ $\displaystyle
Z_{SB}^{-1}{\rm Tr}_{B}\left[e^{-\beta(H_{S}+H_{B})}e^{-{\cal
T}\int_{0}^{\beta}du\tilde{A}(u)\tilde{B}(u)}\right]$ (3)
$\displaystyle\underset{{\rm 2^{d}~{}order}}{\simeq}$ $\displaystyle
Z_{SB}^{-1}e^{-\beta H_{S}}\left[1-\int_{0}^{\beta}du\tilde{A}(u){\rm
Tr}_{B}[e^{-\beta
H_{B}}\tilde{B}(u)]+\int_{0}^{\beta}du_{1}\int_{0}^{u_{1}}du_{2}\tilde{A}(u_{1})\tilde{A}(u_{2}){\rm
Tr}_{B}[e^{-\beta H_{B}}\tilde{B}(u_{1})\tilde{B}(u_{2})]\right]$
$\displaystyle=$ $\displaystyle\frac{Z_{S}Z_{B}}{Z_{SB}}\rho_{S}^{\rm
th}\left[1+\int_{0}^{\beta}du_{1}\int_{0}^{u_{1}}du_{2}\tilde{A}(u_{1})\tilde{A}(u_{2})c_{B}(u_{1}-u_{2})\right],$

where we used in the first line the usual “splitting” formula [48], and
defined the operators $\tilde{X}(u):=e^{u(H_{S}+H_{B})}Xe^{-u(H_{S}+H_{B})}$,
the bath correlation function $c_{B}(u_{1}-u_{2}):={\rm Tr}_{B}[\rho_{B}^{\rm
th}\tilde{B}(u_{1})\tilde{B}(u_{2})]={\rm Tr}_{B}[\rho_{B}^{\rm
th}\tilde{B}(u_{1}-u_{2})\tilde{B}]$ taken in the thermal state $\rho_{B}^{\rm
th}:=Z_{B}^{-1}e^{-\beta H_{B}}$ with $Z_{B}:={\rm Tr}_{B}[e^{-\beta H_{B}}]$,
and the local thermal state of $S$, $\rho_{S}^{\rm th}=Z_{S}^{-1}e^{-\beta
H_{S}}$, with the local partition function $Z_{S}:={\rm Tr}_{S}[e^{-\beta
H_{S}}]$. We also use the property of stationary baths, namely, ${\rm
Tr}_{B}[\rho_{B}^{\rm th}\tilde{B}(u)]=0$. We obtain for the bath correlation
function,

$c_{B}(u)=\int_{0}^{\infty}d\omega J(\omega)\left[e^{-\omega
u}(n_{\omega}+1)+e^{\omega u}n_{\omega}\right]$ (4)

where $n_{\omega}^{\rm th}=(e^{\omega\beta}-1)^{-1}$ is the thermal occupation
at the frequency $\omega$, and the bath spectral density is defined as

$J(\omega):=\sum_{k}g_{k}^{2}\delta(\omega-\omega_{k}).$ (5)

Introducing the eigen-decomposition of the coupling observable
$A=\sum_{\nu}A(\nu)$ such that $[A(\nu),H_{S}]=\nu A(\nu)$,
$A^{{\dagger}}(\nu)=A(-\nu)$, and $\tilde{A}(u)=\sum_{\nu}e^{-\nu u}A(\nu)$,
we have, up to the second order,

$\displaystyle\rho_{S}^{\rm ss,PE}\underset{{\rm
2^{d}~{}order}}{=}\frac{Z_{S}Z_{B}}{Z_{SB}}\rho_{S}^{\rm
th}\left[1+\sum_{\nu,\nu^{\prime}}A(\nu)A^{{\dagger}}(\nu^{\prime})g(\nu,\nu^{\prime})\right],$

where

$g(\nu,\nu^{\prime}):=\int_{0}^{\beta}du_{1}\int_{0}^{u_{1}}du_{2}e^{-\nu
u_{1}+\nu^{\prime}u_{2}}c_{B}(u_{1}-u_{2}),$ (7)

and the superscript “PE” stands for “Perturbative Expansion”. Addtionally,

$\displaystyle Z_{SB}$ $\displaystyle=$ $\displaystyle{\rm
Tr}_{SB}[e^{-\beta(H_{S}+H_{B}+V)}]$ $\displaystyle\underset{{\rm
2^{d}~{}order}}{=}$ $\displaystyle Z_{S}Z_{B}\left[1+\sum_{\nu}{\rm
Tr}_{S}[\rho_{S}^{\rm th}A(\nu)A^{{\dagger}}(\nu)]g(\nu,\nu)\right].$

Expression (II) is Hermitian and is equivalent to the expression obtained in
[47] (up to the initial renormalisation term). An explicit, analytical
expression of the function $g(\nu,\nu^{\prime})$ in term of usual functions is
provided in Appendix A for under-damped (22) and over-damped (23) bath
spectral densities.

###  II.1 Conditions of validity

The above expression (II) is a second order expansion. It provides a good
approximation of the exact steady state $\rho_{S}^{\rm ss}$ Eq.(2) in the
limit of small corrections, when the second order contribution is much larger
than higher order contributions. Based on that, one can build validity
criteria by requiring that the corrections brought by (II) with respect to the
thermal state $\rho_{S}^{\rm th}=Z_{S}^{-1}e^{-\beta H_{S}}$ remain small.
Expressed in different mathematical ways, we present in the following several
potential validity criteria to be tested later with numerical simulations
(section IV).

  * •

A first criterion can be obtained by considering that small corrections should
imply that the global partition function $Z_{SB}$ is close the product of the
local partition functions $Z_{S}Z_{B}$,

${\rm cr}_{1}:=\Big{|}\frac{Z_{SB}}{Z_{B}Z_{S}}-1\Big{|}\ll 1,$ (9)

which turns out to be equivalent to the condition suggested in [47].

  * •

Alternatively, one could consider that the expansion is valid as long as the
corrections to the populations (in the eigenbasis of $H_{S}$) are small,
resulting in the following criterion

${\rm cr}_{2}:=\frac{\Big{|}p_{n}^{\rm ss}-p_{n}^{\rm th}\Big{|}}{p_{n}^{\rm
ss}}\ll 1,$ (10)

for all energy level $n$ of $S$, where $p_{n}^{\rm ss}$ stands for the
populations of the steady state (II) and $p_{n}^{\rm th}$ corresponds to the
populations of the thermal state $\rho_{S}^{\rm th}$, reached in the weak
coupling limit. Note that a coherence-based criterion would typically be
equivalent to the above population-based one.

  * •

Additionally, the quantity

$Q:=\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega},$ (11)

known as the “re-organisation energy” [49, 50, 51], gives a figure of merit of
the coupling energy. Therefore, one can expects the expansion to be valid for
$Q\ll\omega_{S}$, where $\omega_{S}$ stands for the typical energy difference
in $H_{S}$. Thus, we define a third candidate for the validity criterion as

${\rm cr}_{3}:=\frac{Q}{\omega_{S}}\ll 1.$ (12)

Anticipating sections III and IV, we can obtain explicit expressions of $Q$ in
term of the bath parameters for the bath spectral densities used there. For
the under-damped spectral density
$J_{UD}(\omega):=\omega\frac{2}{\pi}\frac{\gamma_{UD}\Omega^{2}\lambda^{2}}{(\Omega^{2}-\omega^{2})^{2}+(\gamma_{UD}\Omega\omega)^{2}}$
(see more details in section III), we obtained $Q_{UD}=\lambda^{2}/\Omega$,
while for the over-damped spectral density
$J_{OD}(\omega):=\alpha\omega\frac{\omega_{c}^{2}}{\omega_{c}^{2}+\omega^{2}}$,
we have $Q_{OD}=\frac{\pi}{2}\alpha\omega_{c}$.

  * •

Finally, considering a two-level system, focus of the comparison section IV,
we can come up with an additional criteria obtained from a special choice of
the system parameters for which the partition function can be easily computed.
More precisely, if we take $r_{x}=r_{y}=0$ (see next section II.2), we can
diagonalise the total Hamiltonian and compute exactly the partition function,
giving a simple but non-trivial expression, $Z_{SB}=e^{\beta Q}Z_{S}Z_{B}$.
Then, in the same spirit as in the first criterion, we can consider that small
corrections imply $|Z_{SB}/Z_{S}Z_{B}-1|\ll 1$, which leads to the simple
criterion

${\rm cr}_{4}:=\beta Q\ll 1.$ (13)

Note that at this stage we only use the setting $r_{x}=r_{y}=0$ as a
mathematical trick to compute $Z_{SB}$, while in the reminder of the paper we
consider $r_{x}\neq 0$, implying $[H_{S},V]\neq 0$, necessary condition for
the applicability of (1). Although it might not seem totally justified to
approximate the order of magnitude of $Z_{SB}$ for arbitrary $r_{x}$ and
$r_{y}$ by its value for $r_{x}=r_{y}=0$, we will see in the following
(section IV), that ${\rm cr}_{4}$ is always very close to ${\rm cr}_{1}$,
justifying afterwards this approximation. Additionally, the factor $e^{\beta
Q}$ is reminiscent of the renormalization factor $e^{\beta QA^{2}}$ due to the
bath interaction [47, 52], so one can conjecture that this criterion could be
extended to arbitrary systems in the form ${\rm cr}_{4}:=\beta Q|A^{2}|\ll 1$.
Finally, this last criterion seems promising because it involves the bath
inverse temperature $\beta$. Indeed, the expansion (II) becomes trivially
valid when the energy scale set by the bath temperature, $k_{B}T=\beta^{-1}$,
is much larger than the system-bath coupling [47] (infinite temperature
limit), suggesting that $\beta$ should play a role in the validity criterion.

We will test these criteria in section IV and see that two of them, ${\rm
cr}_{1}$ and ${\rm cr}_{4}$, seems to indicate particularly well the validity
range of expression (II).

###  II.2 Spin-boson model

In order to obtain explicit comparison with embedding techniques (reaction
coordinate), we choose a specific system, namely the spin-boson model, for
being a widely used system, experimentally as well as theoretically. The
Hamiltonian of the two-level system is of the form

$H_{S}=\frac{\omega_{S}}{2}(r_{x}\sigma_{x}+r_{y}\sigma_{y}+r_{z}\sigma_{z})=\frac{\omega_{s}}{2}\vec{r}.\vec{\sigma},$
(14)

where $\vec{r}$ is a real unit vector of component $r_{x}$, $r_{y}$, and
$r_{z}$ (such that $r_{x}^{2}+r_{y}^{2}+r_{z}^{2}=1$), and $\vec{\sigma}$ is
the Pauli vector of component the Pauli matrices $\sigma_{x}$, $\sigma_{y}$,
and $\sigma_{z}$. Importantly, we will use in the following the notation
$r:=r_{x}+ir_{y}$. We consider a typical coupling with the bath, namely
$A=\sigma_{z}$. The eigen-decomposition takes the form
$A(u)=A(\omega_{s})e^{-u\omega_{s}}+A(-\omega_{s})e^{u\omega_{s}}+A(0)$, with

$\displaystyle A(\omega_{s})$ $\displaystyle=$
$\displaystyle-r|g\rangle\langle e|,$ $\displaystyle A(-\omega_{s})$
$\displaystyle=$ $\displaystyle-r^{*}|e\rangle\langle g|,$ $\displaystyle
A(0)$ $\displaystyle=$ $\displaystyle r_{z}(|e\rangle\langle
e|-|g\rangle\langle g|):=r_{z}\Sigma_{z},$ (15)

where

$\displaystyle|e\rangle:=\frac{(1+r_{z})|+\rangle+r|-\rangle}{\sqrt{2(1+r_{z})}},$
$\displaystyle|g\rangle:=\frac{-r^{*}|+\rangle+(1+r_{z})|-\rangle}{\sqrt{2(1+r_{z})}},$
(16)

are the excited and ground eigenstate of $H_{S}$, respectively. In the above
expression, we used the notation $|+\rangle$, $|-\rangle$ to denote
respectively the excited and ground state of $\sigma_{z}$. Injecting these
expressions in (II) with the use of the explicit expression of the function
$g(\nu,\nu^{\prime})$ provided in Appendix A, we obtain for the reduced steady
state of $S$ in the basis $\\{|e\rangle,|g\rangle\\}$,

$\rho_{S}^{\rm ss,PE}=\begin{pmatrix}p^{\rm ss}_{e}&c_{ge}^{\rm ss*}\\\
c_{ge}^{\rm ss}&p_{g}^{\rm ss}\end{pmatrix},$ (17)

with

$\displaystyle c_{ge}^{\rm ss,PE}$ $\displaystyle=$
$\displaystyle\frac{-2rr_{z}(\beta/\omega_{s})[G(\omega_{s},\beta)-(1+e^{-\omega_{s}\beta})Q/\beta]}{(1+e^{-\omega_{s}\beta})[1+r_{z}^{2}\beta
Q]+|r|^{2}\beta^{2}G(\omega_{s},\beta)},$

and the population,

$\displaystyle p_{e}^{\rm ss,PE}$ $\displaystyle=$
$\displaystyle\frac{e^{-\omega_{s}\beta}(1+r_{z}^{2}\beta Q)-|r|^{2}\beta
G^{\prime}(\omega_{s},\beta)}{(1+e^{-\omega_{s}\beta})[1+r_{z}^{2}\beta
Q]+|r|^{2}\beta^{2}G(\omega_{s},\beta)},$

where $Q$ is the re-organisation energy defined above,

$G(\omega_{s},\beta):=\int_{0}^{1}due^{-\omega_{s}\beta u}c_{B}(u\beta),$ (20)

and
$G^{\prime}(\omega_{s},\beta):=\frac{\partial}{\partial\omega_{s}}G(\omega_{s},\beta)$.
The explicit expression of $G(\omega_{s},\beta)$ and
$G^{\prime}(\omega_{s},\beta)$ in term of usual functions is provided in
Appendix A for both under-damped and over-damped spectral densities
$J_{UD}(\omega)$ (22) and $J_{OD}(\omega)$ (23).

##  III Reaction coordinate

In the perspective of comparing the perturbative expansion approach with
embedding approaches, we briefly review some important features of the
reaction coordinate mapping. Introduced in [53] and further developed in [55,
32, 33, 16, 5, 54], the archetypal application of reaction coordinate is for
the spin-boson model, although it can applied to other systems [6, 19, 5].
Thus, considering the two-level system of the previous section, the spin-boson
model of Hamiltonian
$H_{SB}=\frac{\omega_{s}}{2}\vec{r}.\vec{\sigma}+\sigma_{z}B+H_{B}$ can be
mapped onto [32, 33]

$\displaystyle H_{SB}$ $\displaystyle=$ $\displaystyle
H_{SRCE}:=\frac{\omega_{s}}{2}\vec{r}.\vec{\sigma}+\lambda\sigma_{z}(a^{{\dagger}}+a)+\Omega
a^{{\dagger}}a$ (21)
$\displaystyle+(a^{{\dagger}}+a)B_{E}+H_{E}+(a^{{\dagger}}+a)^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}},$

where $a$ and $a^{{\dagger}}$ are the annihilation and creation operators of
the collective bosonic mode, called the reaction coordinate (RC), defined as
$\lambda A(a^{{\dagger}}+a)=\sum_{k}g_{k}(b_{k}^{{\dagger}}+b_{k})$, and the
system $E$ is a residual bath (the original bath “minus” the collective mode)
of self Hamiltonian $H_{E}$ and coupling to the reaction coordinate through
the operator $B_{E}$. The detailed expressions of the residual bath’s modes
and parameters are not useful in our problem so we refer interested readers to
[32, 33] for further details.

Importantly, the mapping is exact when the original bath $B$ has an under-
damped spectral density,

$J_{UD}(\omega):=\omega\frac{2}{\pi}\frac{\gamma_{UD}\Omega^{2}\lambda^{2}}{(\Omega^{2}-\omega^{2})^{2}+(\gamma_{UD}\Omega\omega)^{2}},$
(22)

where $\lambda$ (frequency), $\Omega$ (frequency) and $\gamma_{UD}$
(dimensionless) characterise respectively the strength of the coupling, the
peak of the spectral density, and its width. According to the reaction
coordinate mapping, the parameters of the collective mode are given directly
by the parameters of the under-damped spectral density [32, 33]: $\lambda$
corresponds to the strength of the coupling between $S$ and the collective
mode, and $\Omega$ is its frequency.

It is also possible to find an approximate mapping when the original bath
spectral density is over-damped, namely of the form,

$J_{OD}(\omega)=\alpha\omega\frac{\omega_{c}^{2}}{\omega_{c}^{2}+\omega^{2}}.$
(23)

where $\omega_{c}$ is sometimes referred to as the cutoff frequency and
$\alpha$ is a dimensionless parameter determining the coupling strength. The
$RC$ coupling and frequency can be expressed in terms of the parameters of
$J_{OD}(\omega)$ as [32, 33]

$\Omega=\gamma\omega_{c}~{}~{}{\rm
and}~{}~{}\lambda=\sqrt{\frac{\pi}{2}\alpha\omega_{c}\Omega}.$ (24)

Written directly in term of the reaction coordinate parameters, the over-
damped spectral density takes the form
$J_{OD}(\omega)=\omega\frac{2}{\pi}\frac{\lambda^{2}\gamma}{\Omega^{2}+\gamma^{2}\omega^{2}}$.
While the expression (23) is a function with two parameters $\alpha$ and
$\omega_{c}$, this later expression contains three parameters. The extra
parameter $\gamma$ is introduced during the reaction coordinate mapping and is
required to be much larger than 1. To understand better the emergence of this
free parameter $\gamma$, one should mention that for over-damped spectral
density, the mapping is actually obtained from an asymptotic limit of the
under-damped spectral density case, as follows [32, 33]. For $\gamma_{UD}\gg
1$, we have $J_{UD}(\omega)\simeq J_{OD}(\omega)$ when setting
$\alpha=\frac{2\gamma_{UD}\lambda^{2}}{\pi\Omega^{2}}$ and
$\omega_{c}=\frac{\Omega}{\gamma_{UD}}$. Then, with these settings, we can use
the reaction coordinate mapping used for the under-damped case, and the
parameter $\gamma$ appearing in (24) is actually $\gamma_{UD}$ which must be
much larger than 1 in order to have the approximate identification
$J_{UD}(\omega)\simeq J_{OD}(\omega)$ . As a conclusion, the reaction
coordinate mapping is not exact for over-damped spectral density, and only
holds under the condition $\Omega\gg\omega_{c}$ (or $\gamma_{UD}\gg 1$).

Now that we have introduced the reaction coordinate mapping for the under-
damped and over-damped spectral densities, we can focus on the steady state.
For weak coupling between the reaction coordinate $RC$ and the residual bath
$E$, which is precisely the situation where the reaction coordinate mapping is
useful, one expects from weak dissipation theory that the extended system
$SRC$ ($S$ and the reaction coordinate) tends to the thermal state at inverse
temperature $\beta$, the inverse temperature of the original and residual
bath,

$\rho_{SRC}^{\rm th}=Z_{SRC}^{-1}e^{-\beta H_{SRC}},$ (25)

where

$H_{SRC}:=H_{S}+\lambda A(a^{{\dagger}}+a)+\Omega a^{{\dagger}}a.$ (26)

This conjecture was indeed benchmarked by numerical techniques (hierarchical
equation of motions) in [32, 33] and Redfield master equation [46] and used in
[16, 17]. However, when the residual coupling between $RC$ and $E$ is not
weak, one expects $\rho_{SRC}^{\rm th}$ to depart from the exact steady state
of $SRC$. Thus, $\rho_{SRC}^{\rm th}$ becomes an approximation of the exact
steady state. How good is this approximation and when exactly does it start
breaking down are the questions which motivated this paper. In the following,
we will refer to $\rho_{SRC}^{\rm th}$ as the “reaction coordinate mapping of
the steady state”, or “mapping of the steady state” in short. From
$\rho_{SRC}^{\rm th}$, the reduced steady state of $S$ is given by

$\rho_{S}^{\rm ss,RC}:={\rm Tr}_{RC}[\rho_{SRC}^{\rm th}].$ (27)

Again, we stress that since in general $\rho_{SRC}^{\rm th}$ is an
approximation of the exact steady state of $SRC$, $\rho_{S}^{\rm ss,RC}$ is
also in general an approximation of $\rho_{S}^{\rm ss}$ (2), the exact steady
state of $S$.

The partial trace over the RC mode can be realised numerically or analytically
via approximate diagonalization of $H_{SRC}$ (see for instance [56, 57]). Note
that the plots presented below were indeed realised using numerical
diagonalisation using QuTiP (with adequate truncation of the RC mode). The
remainder of the paper is mainly dedicated to the comparison of the
predictions of the two approaches, namely comparing (27) with (17), (II.2),
and (II.2). Before that, we introduce a third approximation of the steady
state which will help us in the comparison and is detailed in the next section
III.1.

###  III.1 Perturbative expansion applied to reaction coordinate

Beyond our prime objective to confirm that the perturbative expansion approach
and the reaction coordinate-based approach coincide, at least for some range
of parameters, we also aim at studying the validity range of each approach. In
that perspective, when some discrepancies appear between $\rho_{S}^{\rm
ss,RC}$ and $\rho_{S}^{\rm ss,PE}$, how can we tell that it is because the
reaction coordinate mapping of the steady state, $\rho_{SRC}^{\rm th}$, fails
to faithfully approximate the exact steady state of $RC$, or that it is
because the perturbative expansion stops being valid, or both? How can we
separate the two effects?

We can obtain some insights on these questions by considering a third state,
obtained by applying the general perturbative expansion of section II to
$\rho_{SRC}^{\rm th}$ (25), where the reaction coordinate $RC$ plays the role
of the bath $B$. We denote the resulting state by $\rho_{S}^{\rm ss,PRC}$,
where the superscript “PRC” stands for Perturbative expansion of the RC
mapping. Then, from the point of view of the perturbative expansion,
$\rho_{S}^{\rm ss,RC}$ is the “exact” (containing all orders) version of
$\rho_{S}^{\rm ss,PRC}$. Consequently, the discrepancy between $\rho_{S}^{\rm
ss,RC}$ and $\rho_{S}^{\rm ss,PRC}$ provides information on the validity of
the perturbative expansion. Conversely, $\rho_{S}^{\rm ss,PE}$ and
$\rho_{S}^{\rm ss,PRC}$ are both expansions to the same order, of the original
problem and of the reaction coordinate mapping, respectively. Then, from the
point of view of the reaction coordinate mapping, $\rho_{S}^{\rm ss,PE}$ is
the exact version of $\rho_{S}^{\rm ss,PRC}$. Thus, by observing the
discrepancy between $\rho_{S}^{\rm ss,PE}$ and $\rho_{S}^{\rm ss,PRC}$ we can
obtain information on the performance of the reaction coordinate mapping. We
will use in the next section these discrepancies $\rho_{S}^{\rm ss,RC}$ versus
$\rho_{S}^{\rm ss,PRC}$ and $\rho_{S}^{\rm ss,PE}$ versus $\rho_{S}^{\rm
ss,PRC}$ to gain precious information on the range of validity of each
approach.

Applying the general perturbative expansion of section II to $\rho_{SRC}^{\rm
th}$ we obtain the same form as (II), namely,

$\displaystyle\rho_{S}^{\rm ss,PRC}\underset{{\rm
2^{d}~{}order}}{=}\frac{Z_{S}Z_{B}}{Z_{SB}}\rho_{S}^{\rm
th}\left[1+\sum_{\nu,\nu^{\prime}}A(\nu)A^{{\dagger}}(\nu^{\prime})g(\nu,\nu^{\prime})\right].$

The functions $g(\nu,\nu^{\prime})$, $c_{B}(u)$, $G(\nu,\beta)$, and
$G^{\prime}(\nu,\beta)$, have the same general expression as the one detailed
in Appendix A but using the following spectral density

$J(\omega)=\lambda^{2}\delta(\omega-\Omega).$ (29)

Thus, for the steady state populations and coherences, it leads to the same
expressions as (II.2) and (II.2), respectively, substituting $G(\nu,\beta)$
and $G^{\prime}(\nu,\beta)$ by $G^{\rm PRC}(\nu,\beta):=C^{\rm
PRC}(-\nu)+e^{-\nu\beta}C^{\rm PRC}(\nu)$, and $G^{\rm
PRC^{\prime}}(\nu,\beta):=-C^{\rm PRC^{\prime}}(-\nu)-\beta
e^{-\nu\beta}C^{\rm PRC}(\nu)+e^{-\nu\beta}C^{\rm PRC^{\prime}}(\nu)$, with

$C^{\rm PRC}(\nu):=\frac{\lambda^{2}}{\beta}\left(\frac{n_{\Omega}^{\rm
th}+1}{\Omega-\nu}-\frac{n_{\Omega}^{\rm th}}{\Omega+\nu}\right)$ (30)

and

$C^{\rm
PRC^{\prime}}(\nu):=\frac{\lambda^{2}}{\beta}\left(\frac{n_{\Omega}^{\rm
th}+1}{(\Omega-\nu)^{2}}+\frac{n_{\Omega}^{\rm th}}{(\Omega+\nu)^{2}}\right).$
(31)

(a)
(b)

Figure 1: Steady state coherences $c_{ge}^{\rm PE,OD}$ (39) (red thick solid
line), $c_{ge}^{\rm PE,UD}$ (37) (orange thick solid line), $c_{ge}^{\rm PRC}$
(35) (blue dashed line), and $c_{ge}^{\rm RC}$ (33) (sparse dotted line), in
function of the inverse temperature $\beta$ in unit of $\omega_{S}^{-1}$, for
(a) $\lambda/\omega_{S}=0.5$ and (b) $\lambda/\omega_{S}=1.5$. The other
parameters are given by $\Omega/\omega_{S}=10$, $\gamma_{UD}=0.1$,
$r_{z}=\sqrt{0.75}$, $r=r_{x}+ir_{y}=0.5$.

##  IV Comparison

In this section, we compare the steady state $\rho_{S}^{\rm ss,RC}$ given by
the reaction coordinate (27), the steady state $\rho_{S}^{\rm ss,PRC}$ (III.1)
given by the perturbative expansion of the reaction coordinate, and the steady
states $\rho_{S}^{\rm ss,PE,UD}$ and $\rho_{S}^{\rm ss,PE,OD}$, given
respectively by the perturbative expansion of the original problem (17) for
under-damped bath spectral density $J_{UD}(\omega)$ (22) and over-damped bath
spectral density $J_{OD}(\omega)$ (23). We denote the coherence and excited
population in the eigenbasis $\\{|e\rangle,|g\rangle\\}$ of $H_{S}$ as

$\displaystyle p_{e}^{\rm RC}:=\langle e|\rho_{S}^{\rm ss,RC}|e\rangle,$ (32)
$\displaystyle c_{ge}^{\rm RC}:=\langle g|\rho_{S}^{\rm ss,RC}|e\rangle,$ (33)
$\displaystyle p_{e}^{\rm PRC}:=\langle e|\rho_{S}^{\rm ss,PRC}|e\rangle,$
(34) $\displaystyle c_{ge}^{\rm PRC}:=\langle g|\rho_{S}^{\rm
ss,PRC}|e\rangle,$ (35) $\displaystyle p_{e}^{\rm PE,UD}:=\langle
e|\rho_{S}^{\rm ss,PE,UD}|e\rangle,$ (36) $\displaystyle c_{ge}^{\rm
PE,UD}:=\langle g|\rho_{S}^{\rm ss,PE,UD}|e\rangle,$ (37) $\displaystyle
p_{e}^{\rm PE,OD}:=\langle e|\rho_{S}^{\rm ss,PE,OD}|e\rangle,$ (38)
$\displaystyle c_{ge}^{\rm PE,OD}:=\langle g|\rho_{S}^{\rm
ss,PE,OD}|e\rangle.$ (39)

Fig. 1 presents the plots of the steady state coherences as given by
$c_{ge}^{\rm PE,OD}$ (red thick solid line), $c_{ge}^{\rm PE,UD}$ (orange
thick solid line), $c_{ge}^{\rm PRC}$ (blue dashed line), and $c_{ge}^{\rm
RC}$(sparsely dotted line), in function of the inverse temperature $\beta$ (in
unit of $\omega_{S}^{-1}$). The panel (a) corresponds to a coupling
$\lambda/\omega_{S}=0.5$ and the panel (b) to $\lambda/\omega_{S}=1.5$. The
other parameters are chosen as follows, $\Omega/\omega_{S}=10$,
$\gamma_{UD}=0.1$ (dimensionless), $r_{z}=\sqrt{0.75}$, $r=r_{x}+ir_{y}=0.5$,
meaning that $r_{y}=0$ (note than one would obtain a similar behavior up to a
$\pi/2$ phase for the coherence with the alternative choice $r_{x}=0$, and
$r_{y}=0.5$).

One can see a very good agreement between $c_{ge}^{\rm RC}$ (sparse dots) and
$c_{ge}^{\rm PRC}$ (blue dashed line) at high temperature (low inverse
temperature $\beta$), but this agreement slowly deteriorates beyond $\beta\sim
4\omega_{S}^{-1}$ in panel (a), and beyond $\beta\sim\omega_{S}^{-1}$ in panel
(b). By contrast, the agreement between $c_{ge}^{\rm PE,OD}$ (red solid line)
and $c_{ge}^{\rm PRC}$ (blue dashed line) is only good at very high
temperature and it deteriorates quickly as the temperature decreases. Finally,
$c_{ge}^{\rm PE,UD}$ (orange solid line) and $c_{ge}^{\rm PRC}$ (blue dashed
line) coincide perfectly at high temperature, and only a very small
discrepancy appears at small temperatures.

(a)
(b)
(c)
(d)

Figure 2: Steady state populations in function of the inverse temperature
$\beta\in[5;6]$ (in unit of $\omega_{S}^{-1}$) for panels (a) and (b), and
$\beta\in[0.5;0.51]$ for (c) and (d). The coupling strength is
$\lambda/\omega_{S}=0.5$ for panels (a) and (c), and $\lambda/\omega_{S}=1.5$
for panel (b) and (d). Following the same color convention as in Fig. 1, the
red thick solid line corresponds to $p_{e}^{\rm PE,OD}$ (38), the orange thick
solid line corresponds to $p_{e}^{\rm PE,OD}$ (36), the blue dashed line
corresponds to $p_{e}^{\rm PRC}$ (34), and the sparse dotted line corresponds
to $p_{e}^{\rm RC}$ (32). The dotted grey line represents the thermal excited
population $p_{e}^{\rm th}:=\langle e|\rho_{S}^{\rm th}|e\rangle$ at inverse
temperature $\beta$. The other parameters are as in Fig 1.

Fig. 2 is the counter-part of Fig. 1 for the excited population. Note that
instead of plotting directly the populations for $\omega_{S}\beta$ from 0 to
10, we zoom in and consider two sections, otherwise all five curves would be
indistinguishable: the panels (a) and (b) represents the steady state excited
population as a function of the inverse temperature $\beta$ in the interval
$[5;6]$ (in unit of $\omega_{S}^{-1}$), while panels (c) and (d) correspond to
$\beta\in[0.5;0.51]$. Additionally, panels (a) and (c) correspond to a
coupling $\lambda/\omega_{S}=0.5$, while panel (b) and (d) to
$\lambda/\omega_{S}=1.5$. The color convention is the same as in Fig. 1,
namely, $p_{e}^{\rm PE,OD}$ (38) (red solid line), $p_{e}^{\rm PE,UD}$ (36)
(orange solid line), $p_{e}^{\rm PRC}$ (34) (the blue dashed line), and
$p_{e}^{\rm RC}$ (sparsely dotted line) (32). The grey dotted line represents
the thermal excited population $p_{e}^{\rm th}:=\langle e|\rho_{S}^{\rm
th}|e\rangle$ at inverse temperature $\beta$. The other parameters are chosen
as in Fig. 1.

The conclusions are the same as in Fig. 1, namely, the agreement between
$p_{e}^{\rm RC}$ and $p_{e}^{\rm PRC}$ is very good at high temperature and
deteriorates more at low temperature when the system-bath coupling is larger.
However, the agreement between $p_{e}^{\rm PE,OD}$ and $p_{e}^{\rm PRC}$ is
relatively good only at high temperature, and large discrepancies appear at
low temperatures, while $p_{e}^{\rm PE,UD}$ and $p_{e}^{\rm PRC}$ coincide
perfectly at both ranges of temperatures, even in strong coupling.

(a)
(b)
(c)
(d)

Figure 3: Steady state populations (a,c) and coherences (b,d) in function of
the coupling strength $\lambda$ (in unit of $\omega_{S}$) at inverse
temperature (a,b) $\omega_{S}\beta=0.5$ and (c,d) $\omega_{S}\beta=5$. The
color conventions are the same as in Figs. 1 and 2. The remainder of the
parameters are chosen as in previous figures, namely $\Omega/\omega_{S}=10$,
$\gamma_{UD}=0.1$, $r_{z}=\sqrt{0.75}$, $r=r_{x}+ir_{y}=0.5$.

Finally, Fig. 3 presents the plots of the steady state coherences and
populations in function of the coupling strength $\lambda$ (in unit of
$\omega_{S}$) for (a)-(b) $\omega_{S}\beta=0.5$ and (c)-(d)
$\omega_{S}\beta=5$. The color conventions are the same as in Fig. 1 and Fig.
2, as well as the remaining parameters.

As a first observation, one can see that the plots for the coherences and
populations have almost identical shapes, which will be confirmed in Figs. 4
and 7. Secondly, the agreement between $c_{ge}^{\rm ss,RC}$ (sparse dotted
line) and $c_{ge}^{\rm ss,PRC}$ (blue dashed line) (as well as $p_{e}^{\rm
ss,RC}$ and $p_{e}^{\rm ss,PRC}$) is excellent below $\lambda\sim 2\omega_{S}$
at high temperature $\omega_{S}\beta=0.5$, and below $\lambda\sim
0.75\omega_{S}$ at low temperature $\omega_{S}\beta=5$, while it starts
deteriorating beyond these values of coupling strength. For $c_{ge}^{\rm
PE,OD}$ (red thick line) and $c_{ge}^{\rm PRC}$ (blue dashed line) (as well as
$p_{e}^{\rm PE,OD}$ and $p_{e}^{\rm PRC}$), the agreement is good until
$\lambda\sim 1$ for $\omega_{S}\beta=0.5$, but for $\omega_{S}\beta=5$ the
agreement drops quickly beyond $\lambda\sim 0.2\omega_{S}$. By contrast, the
agreement between $c_{ge}^{\rm PE,UD}$ (orange thick line) and $c_{ge}^{\rm
PRC}$ (blue dashed line) (as well as $p_{e}^{\rm PE,UD}$ and $p_{e}^{\rm
PRC}$) is almost perfect for both values of $\beta$ and for all $\lambda$.

Conclusion–. From Fig. 1 \- 3, we can conclude that the two approaches do
coincide on an extended region of parameters for under-damped spectral
densities, namely from high to low temperatures for limited coupling
$\lambda=0.5\omega_{S}$, and even for large coupling $\lambda=1.5\omega_{S}$
at high temperature $\omega_{S}\beta=0.5$. However, for over-damped spectral
densities, the two approaches coincide only at high temperatures.

(a)
(b)

Figure 4: Panel (a): the yellow and green dots correspond to $d_{\rm coh}^{\rm
exp}$ in function of $\lambda$ (in unit of $\omega_{S}$) for
$\omega_{S}\beta=0.5$ and $\omega_{S}\beta=5$, respectively, and the black and
purple large dots (in the background of the yellow and green dots) correspond
to $d_{\rm pop}^{\rm exp}$ in function of $\lambda$ for $\omega_{S}\beta=0.5$
and $\omega_{S}\beta=5$, respectively. Panel (b): the yellow and green dots
correspond to $d_{\rm coh}^{\rm exp}$ in function of $\beta$ (in unit of
$\omega_{S}^{-1}$) for $\lambda=0.5\omega_{S}$ and $\lambda=1.5\omega_{S}$,
respectively, and the black and purple large dots correspond to $d_{\rm
pop}^{\rm exp}$ in function of $\beta$ for $\lambda=0.5\omega_{S}$ and
$\lambda=1.5\omega_{S}$, respectively. For both panels, the remainder of the
parameters are chosen as in previous figures, namely $\Omega/\omega_{S}=10$,
$\gamma_{UD}=0.1$, $r_{z}=\sqrt{0.75}$, $r_{x}+ir_{y}=0.5$

.

###  IV.1 Discrepancies due to the perturbative expansion

In order to obtain a more precise and quantitative criterion of “good” and
“bad” agreement, we introduce the relative discrepancies

$\displaystyle d_{\rm coh}^{\rm exp}:=(c_{ge}^{\rm RC}-c_{ge}^{\rm
PRC})/c_{ge}^{\rm RC},$ $\displaystyle d_{\rm pop}^{\rm exp}:=(p_{e}^{\rm
RC}-p_{e}^{\rm PRC})/(p_{e}^{\rm RC}-p_{e}^{\rm th}),$

where one should note that the relative discrepancy related to the population
is defined with respect to the deviation from the thermal excited population
$p_{e}^{\rm th}=\langle e|\rho_{S}^{\rm th}|e\rangle$ (at inverse temperature
$\beta$). As already discussed in section III.1, these relative discrepancies
provide information on the validity of the perturbative expansion.

Fig. 4 (a) provides the relative discrepancies associated with Fig. 3, namely,
the yellow and green dots correspond to $d_{\rm coh}^{\rm exp}$ for
$\omega_{S}\beta=0.5$ and $\omega_{S}\beta=5$, respectively, while the black
and purple large dots (in the background of the yellow and green dots)
correspond to $d_{\rm pop}^{\rm exp}$ for $\omega_{S}\beta=0.5$ and
$\omega_{S}\beta=5$, respectively. The remainder of the parameters are chosen
as in previous figures, namely $\Omega/\omega_{S}=10$, $\gamma_{UD}=0.1$,
$r_{z}=\sqrt{0.75}$, $r=r_{x}+ir_{y}=0.5$. Firstly, one can see that the
relative discrepancies is exactly the same for coherences and for populations
at both temperatures, confirming observations from Fig. 3. More importantly,
adopting the standard 10$\%$ error criterion, these plots testify that the
perturbative expansion is valid up to a coupling strength $\lambda\sim
2\omega_{S}$ at $\omega_{S}\beta=0.5$, and up to $\lambda\sim 0.5\omega_{S}$
at $\omega_{S}\beta=5$, which also coincides with what we can see from Fig. 3.

Fig. 4 (b) provides the relative discrepancies of Figs. 1 and 2. Adopting a
similar color convention as panel (a), the yellow and green dots correspond to
$d_{\rm coh}^{\rm exp}$ for $\lambda=0.5\omega_{S}$ and
$\lambda=1.5\omega_{S}$, respectively, while the black and purple large dots
correspond to $d_{\rm pop}^{\rm exp}$ for $\lambda=0.5\omega_{S}$ and
$\lambda=1.5\omega_{S}$, respectively. Thus, as for panel (a), the relative
discrepancies are also exactly the same for coherences and populations.
Conclusion–. These plots confirm the observations from Figs. 1 and 2, namely
that the perturbative expansion is roughly valid up to $\omega_{S}\beta\sim 5$
for $\lambda=0.5\omega_{S}$, and up to $\omega_{S}\beta\sim 1$ for
$\lambda=1.5\omega_{S}$.

###  IV.2 Benchmarking the validity criteria for perturbative expansions

Using the observations from the previous figures, we can benchmark the
capacity of the validity criteria introduced in section II.1 to pinpoint the
actual range of validity of the perturbative expansion. In Fig. 5, we plot the
four criteria ${\rm cr}_{1}$ (9) (black line), ${\rm cr}_{2}$ (10) (blue
line), ${\rm cr}_{3}$ (12) (green line), ${\rm cr}_{4}$ (13) (red dashed
line), in function of $\lambda$ (in unit of $\omega_{S}$), for
$\omega_{S}\beta=0.5$ in panel (a), and $\omega_{S}\beta=5$ in panel (b).
Considering that ${\rm cr}_{i}\ll 1$ means ${\rm cr}_{i}\lesssim 0.1$, one can
see that only criteria ${\rm
cr}_{1}=\Big{|}\frac{Z_{SB}}{Z_{B}Z_{S}}-1\Big{|}$ and ${\rm cr}_{4}=\beta Q$,
which, interestingly, coincide exactly, indicate a range of validity in
agreement with our conclusions from Figs. 3 and 4. More precisely, at inverse
temperature $\omega_{S}\beta=0.5$ ($\omega_{S}\beta=5$), ${\rm cr}_{1}$ and
${\rm cr}_{4}$ indicate a validity of the perturbative expansion up to a
coupling strength $\lambda\sim 1.5\omega_{S}$ ($\lambda\sim 0.5\omega_{S}$),
in agreement with the value $\lambda\sim 2\omega_{S}$ ($\lambda\sim
0.5\omega_{S}$ to $0.75\omega_{S}$) from Figs. 3 and 4.

(a)
(b)

Figure 5: Plots of the validity criteria introduced in section II.1 in
function of the coupling strength $\lambda$ (in unit of $\omega_{S}$) for (a)
$\omega_{S}\beta=0.5$ and (b) $\omega_{S}\beta=5$. Criteria ${\rm cr}_{1}$ and
${\rm cr}_{4}$ coincide almost exactly and correspond to the black solid line
and red dashed line, respectively; ${\rm cr}_{2}$ corresponds to the blue
solid line, and ${\rm cr}_{3}$ corresponds to the green solid line. The
remainder of the parameters are chosen as in previous figures, namely
$\Omega/\omega_{S}=10$, $\gamma_{UD}=0.1$, $r_{z}=\sqrt{0.75}$,
$r=r_{x}+ir_{y}=0.5$.

Similarly, in Fig. 6, we plot the same validity criteria, but now as functions
of the inverse temperature, using the same color convention as in the previous
figure 5. One can see that the conclusions are the same: only ${\rm cr}_{1}$
and ${\rm cr}_{4}$, which coincide almost perfectly, indicate a range of
validity in agreement with our previous conclusions from Figs. 1, 2 and 4.

(a)
(b)

Figure 6: Plots of the validity criteria introduced in section II.1 in
function of $\beta$ (in unit of $\omega_{S}^{-1}$) for (a)
$\lambda=0.5\omega_{S}$ and (b) $\lambda=5\omega_{S}$. The color conventions
are the same as in Fig. 5, namely: ${\rm cr}_{1}$ (black solid line), ${\rm
cr}_{4}$ (red dashed line), ${\rm cr}_{2}$ (blue solid line), and ${\rm
cr}_{3}$ (green solid line). The remainder of the parameters are chosen as in
previous figures, namely $\Omega/\omega_{S}=10$, $\gamma_{UD}=0.1$,
$r_{z}=\sqrt{0.75}$, $r=r_{x}+ir_{y}=0.5$.

Thus, we obtain a very simple criterion for the validity of the perturbative
expansion incorporating both the coupling strength and the temperature
aspects,

$\beta Q\lesssim 0.1$ (41)

This takes the explicit form $\beta\frac{\lambda^{2}}{\Omega}\lesssim 0.1$ for
an under-damped spectral density of parametrization (22), and
$\frac{\pi}{2}\beta\alpha\omega_{c}\lesssim 0.1$ for an over-damped spectral
density of parametrization (23).

Importantly, we also benchmarked this result for $\Omega\sim\omega_{S}$ and
$\Omega\ll\omega_{S}$, both in function of $\beta$ and $\lambda$ (see
additional plots provided in Appendix B). In all situations, we confirm that
$\beta Q\lesssim 0.1$ is a surprisingly good validity criterion. In
particular, even in regimes where the criterion $Q\ll\omega_{S}$ is totally
misleading, $\beta Q\lesssim 0.1$ indicates accurately the region where the
expansion stops being valid.

One should note that there is something unexpected in this result. It is
comparing the energy scale of the system with the energy scale of the coupling
that one sometimes defines weak/strong/ultrastrong coupling: the system’s
transition energy scale (here $\omega_{S}$) is the reference. However, our
results point at a criterion which is independent of $\omega_{S}$, namely
$\beta Q\lesssim 0.1$, while the criterion based on $\omega_{S}$ is wrong most
of the time. This could suggest that either the breakdown of the validity of
the perturbative expansion is not strictly related to strong coupling, or the
definition of strong coupling might not always be related to the system’s
energy scale.

As a side comment, ${\rm cr}_{1}$ and ${\rm cr}_{4}$ are exactly equal for
$r_{x}=r_{y}=0$, by definition. However, their exact agreement, at least for
the considered range of parameters (see also plots in Appendix B), is quite
surprising and justifies afterwards the assumption made in section II.1 that
the order of magnitude of $Z_{SB}$ does not depend on the spin orientation.

###  IV.3 Benchmarking the reaction coordinate mapping

In this section we focus on the other aspect of the problem: how well does the
reaction coordinate mapping of the steady state (25) approximate the steady
state of the original problem, $\rho_{S}^{\rm ss}$ (2)? We already saw in
Figs. 1, 2, and 3 that it depends strongly on the bath spectral density as
well as on the bath temperature. As in the previous section IV.2, in order to
obtain more quantitative information on the performance of the mapping of the
steady state, we introduce the following relative discrepancies,

$\displaystyle d_{\rm coh}^{\rm map,UD}:=(c_{ge}^{\rm PE,UD}-c_{ge}^{\rm
PRC})/c_{ge}^{\rm PE,UD},$ $\displaystyle d_{\rm pop}^{\rm
map,UD}:=(p_{e}^{\rm PE,UD}-p_{e}^{\rm PRC})/(p_{e}^{\rm PE,UD}-p_{e}^{\rm
th}),$ $\displaystyle d_{\rm coh}^{\rm map,OD}:=(c_{ge}^{\rm
PE,OD}-c_{ge}^{\rm PRC})/c_{ge}^{\rm PE,OD},$ $\displaystyle d_{\rm pop}^{\rm
map,OD}:=(p_{e}^{\rm PE,OD}-p_{e}^{\rm PRC})/(p_{e}^{\rm PE,OD}-p_{e}^{\rm
th}).$ (42)

and plot them in function of $\lambda$ and $\beta$ in Fig. 7. In Fig. 7 (a),
the yellow and green thin solid line represent $d_{\rm coh}^{\rm map,UD}$ as a
function of $\lambda$ for $\omega_{S}\beta=0.5$ and $\omega_{S}\beta=5$,
respectively, while the black and purple large solid line represent $d_{\rm
pop}^{\rm map,UD}$ in function of $\lambda$ for $\omega_{S}\beta=0.5$ and
$\omega_{S}\beta=5$, respectively. In Fig. 7 (c), the same quantities are
plotted for the over-damped spectral densities. Interestingly, we can see that
the relative discrepancies are independent of the coupling strength. In some
sense it means that both perturbative expansions $\rho_{S}^{\rm ss,PRC}$ and
$\rho_{S}^{\rm ss,PE}$ drift away in parallel from their respective exact
states $\rho_{S}^{\rm ss,RC}$ and $\rho_{S}^{\rm ss}$. This is a good
indication that the strategy of comparing $\rho_{S}^{\rm ss,PRC}$ and
$\rho_{S}^{\rm ss,PE}$ does result in getting rid of discrepancies stemming
from the perturbative expansion and thus retains only discrepancies stemming
from the reaction coordinate mapping, as if we were measuring directly the
discrepancies between $\rho_{S}^{\rm ss,RC}$ and $\rho_{S}^{\rm ss}$. From
Fig. 7 (a) we also have the confirmation that the mapping of the steady state
performs well for narrow under-damped spectral densities, and fails for over-
damped spectral densities, panel (c).

In panel (b), the yellow line corresponds to $d_{\rm coh}^{\rm map,UD}$ in
function of $\beta$ (in unit of $\omega_{S}^{-1}$) and the black thick line
corresponds to $d_{\rm pop}^{\rm map,UD}$ also in function of $\beta$ (both
for arbitrary $\lambda$ since $d_{\rm pop/coh}^{\rm map,UD}$ is independent of
$\lambda$). The same quantities are plotted in panel (d) for over-damped
spectral density. The main message of these plots is that the reaction
coordinate mapping seems to always perform well at high temperatures.

In the following, we briefly explain the reasons behind the failure of the
reaction coordinate mapping of the steady state for over-damped spectral
densities at arbitrary temperatures, and we also explain its universal success
at high temperature.

(a)
(b)
(c)
(d)

Figure 7: Plots of the relative discrepancies $d_{\rm pop/coh}^{\rm UD/OD}$ in
function of the coupling strength $\lambda$ and $\beta$. Panel (a) (panel
(c)): the yellow and green thin solid line represent $d_{\rm coh}^{\rm
map,UD}$ ($d_{\rm coh}^{\rm map,OD}$) as a function of $\lambda$ (in unit of
$\omega_{S}$) for $\omega_{S}\beta=0.5$ and $\omega_{S}\beta=5$, respectively,
while the black and purple large solid line represent $d_{\rm pop}^{\rm
map,UD}$ ($d_{\rm pop}^{\rm map,OD}$) in function of $\lambda$ for
$\omega_{S}\beta=0.5$ and $\omega_{S}\beta=5$, respectively. Panel (b) (panel
(d)): the yellow line corresponds to $d_{\rm coh}^{\rm map,UD}$ ($d_{\rm
coh}^{\rm map,OD}$) in function of $\beta$ (in unit of $\omega_{S}^{-1}$) and
the black thick line corresponds to $d_{\rm pop}^{\rm map,UD}$ ($d_{\rm
pop}^{\rm map,OD}$) also in function of $\beta$, both for arbitrary $\lambda$.
The other parameters are chosen as in previous figures, namely
$\Omega/\omega_{S}=10$, $\gamma_{UD}/\omega_{S}=0.1$, $r_{z}=\sqrt{0.75}$,
$r_{x}+ir_{y}=0.5$.

####  IV.3.1 Reasons for discrepancies

Since the reaction coordinate mapping is exact for under-damped bath spectral
densities [32, 33, 16], one might not be surprised that we observed a good
performance for such spectral densities. On the other hand, for under-damped
spectral densities of increasing width (determined by $\gamma_{UD}$), the
reaction coordinate mapping is still exact, but one can verify that the under-
damped spectral density becomes indistinguishable from an over-damped spectral
density, and the steady state coherences and populations also become
indistinguishable from the ones of an over-damped spectral density. Thus, how
can we have an exact mapping giving a wrong stead state? As already commented
in section III, the reason for this apparent contradiction is that the
reaction coordinate mapping of under-damped spectral densities is exact for
the dynamics, and the steady state $\rho_{SRC}^{\rm th}$ (25) is only an
approximation of the actual steady state of the dynamics [32, 33, 16, 17, 46].
Thus, as already stressed above, the observed failure of the reaction
coordinate mapping of the steady state is not a failure of the reaction
coordinate mapping per se, but is a break down of the approximation consisting
in equating the exact steady state of $SRC$ by $\rho_{SRC}^{\rm th}$ (25).
This breakdown can be understood from three related point of view.

  * •

Although the mapping is exact for under-damped spectral densities of arbitrary
width, the strength of the coupling between the residual bath and the reaction
coordinate grows as $\gamma_{UD}$. Therefore, for increasing spectral widths,
one should expect $\rho_{SRC}^{\rm th}$ to depart from the exact steady state.
In particular, one expects a steady state of the form $\rho_{SRC}^{\rm
ss}={\rm Tr}_{E}[\rho_{SRCE}^{\rm th}]\neq\rho_{SRC}^{\rm th}$. More
precisely, one can show [16] that ${\rm Tr}_{RC}[\rho_{SRC}^{\rm th}]$ is
equal to ${\rm Tr}_{B}[\rho_{SB}^{\rm th}]$ at lowest order in the residual
coupling between $RC$ and $E$. Thus, for increasing $\gamma_{UD}$, and
therefore increasing residual coupling, discrepancies between $\rho_{SRC}^{\rm
ss}$ and $\rho_{SRC}^{\rm th}$ as well as $\rho_{S}^{\rm ss}$ and
$\rho_{S}^{\rm ss,RC}$ increase.

  * •

Alternatively, this can be seen directly from the expressions we obtained for
the general perturbative expansion (II). The steady state (II) depends on the
function $g(\nu,\nu^{\prime})$, which is entirely determined by the bath
correlation function $c_{B}(u)$ (4), which is itself ultimately determined by
the bath spectral density $J(\omega)$. Thus, when approximating the steady
state $\rho_{S}^{\rm ss}={\rm Tr}_{B}[\rho_{SB}^{\rm th}]$ (2) by ${\rm
Tr}_{RC}[\rho_{SRC}^{\rm th}]$ (27), we are ultimately approximating the
original bath spectral density by a single mode, represented by the reaction
coordinate. In other words, we are approximating the original bath spectral
density $J(\omega)$ by $\lambda^{2}\delta(\omega-\Omega)$. This approximation
is reasonable if $J(\omega)$ is a narrow spectral density centered in
$\Omega$, but is not justified for a broad spectral density. Thus, one expects
that ${\rm Tr}_{RC}[\rho_{SRC}^{\rm th}]$ becomes increasingly distant from
$\rho_{S}^{\rm ss}={\rm Tr}_{B}[\rho_{SB}^{\rm th}]$ (2) as the spectral width
increases.

  * •

Final viewpoint, strongly related to the first one. One can show that the
steady state $\rho_{SRC}^{\rm th}$ (25) is actually the steady state of the
master equation derived in the Supplementary Material of [32] when applying
the secular approximation. However, the secular approximation is valid when
${\rm max}|\omega_{SRC}-\omega_{SRC}^{\prime}|^{-1}\ll\tau_{D}$, where
$\tau_{D}$ denotes the dissipation timescale induced by the action of the bath
and $\omega_{SRC}$ denotes the Bohr frequencies of the extended system $SRC$.
A rough analysis show that $\tau_{D}\sim(\pi\gamma_{UD}\omega_{S})^{-1}$, so
that one expects the secular approximation to become unjustified for growing
$\gamma_{UD}$, and thus a steady state increasingly distinct from
$\rho_{SRC}^{\rm th}$. This question has been analyzed in great details in
[6], and one of the conclusion actually limits the strength of this last
argument: the authors show that the reaction coordinate mapping of the steady
state might actually be valid way beyond the supposed validity of the secular
approximation.

As a rule of thumb, from observations coming from the above plots and
additional plots (not shown), one can consider that the reaction coordinate
steady state performs well, meaning $d_{\rm pop/coh}^{\rm map,UD}\leq 0.1$, as
long as the spectral width $\gamma_{UD}$ is smaller than $\sim 3/\Omega\beta$.

####  IV.3.2 Universal faithful mapping of the steady state at high
temperature

Contrasting with the breakdown of the reaction coordinate mapping of the
steady state for broad bath spectral densities at arbitrary temperatures, the
mapping seems to be always faithful at high temperatures (see Fig. 7 (d)).
This can be seen as follows. Assuming that the bath spectral density
$J(\omega)$ vanishes for $\omega\geq 2/\beta$, the correlation function can be
approximated by

$\displaystyle c_{B}(u)$ $\displaystyle:=$ $\displaystyle{\rm
Tr}_{B}\rho_{B}^{\rm th}B(u)B$ (43) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega J(\omega)\left[e^{-\omega
u}(n_{\omega}+1)+e^{\omega u}n_{\omega}\right]$ $\displaystyle\sim$
$\displaystyle\int_{0}^{\infty}d\omega
J(\omega)\frac{2}{\omega\beta}=\frac{2Q}{\beta},$

where we use the approximation $e^{-\omega u}(n_{\omega}+1)+e^{\omega
u}n_{\omega}\sim\frac{2}{\omega\beta}$ valid for $\omega\beta\leq 2$
(reminding that the variable $u$ belongs to $[0;\beta]$). Then, applying this
result to the bath spectral densities we have been considering,
$J_{OD}(\omega)$ and $J_{UD}(\omega)$, both vanishing for $\omega\gg\Omega$,
one expects to have $c_{B}(u)\sim\frac{2Q}{\beta}$ for both spectral densities
as soon as $\Omega\ll 2/\beta$. Additionally, considering the effective
spectral density representing the reaction coordinate
$J_{RC}(\omega)=\delta(\omega-\Omega)$, we have $c_{B}(u)\sim\frac{2Q}{\beta}$
when $\Omega\leq 2/\beta$. Thus, for $\Omega\ll 2/\beta$,
$G(\omega_{s},\beta)$, and therefore $G^{\prime}(\omega_{s},\beta)$, become
independent of the form of the bath spectral density, retaining only a
dependence on $Q$. Then, for $\Omega\ll 2/\beta$, we should have the same
steady state for any bath spectral densities of same re-organization energy
$Q$. This is what we observed in Fig. 7 (d), and confirmed in Fig. 8, as soon
as

$\Omega\beta\leq 1.$ (44)

This hols for arbitrary spectral width $\gamma_{UD}$, and might also hold for
arbitrary coupling strength (again, according to our conclusions from Fig. 7).

Figure 8: Plots of $d_{\rm coh}^{\rm map,UD/OD}$ (yellow thin line) and
$d_{\rm pop}^{\rm map,UD/OD}$ (black thick line) in function of the inverse
temperature $\beta$ (in unit of $\omega_{S}$) for $\Omega/\omega_{S}=1/\beta$,
$\gamma_{UD}=20$ (any value larger than 20 gives the same plot),
$r_{z}=\sqrt{0.75}$, $r=r_{x}+ir_{y}=0.5$. The plots of $d_{\rm coh}^{\rm
map,OD}$ and $d_{\rm coh}^{\rm map,UD}$ ($d_{\rm pop}^{\rm map,OD}$ and
$d_{\rm pop}^{\rm map,UD}$) are indistinguishable. Additionally, we chose
$\lambda=1.5\omega_{S}$, but the plots are actually independent of $\lambda$
as shown in Fig. 7.

##  V Conclusion

We compare the perturbative expansion of the mean force Gibbs state with the
approximate steady state (25) from the reaction coordinate mapping. We reach
our first objective by showing the agreement of these two approaches, for some
ranges of parameters and focusing on the spin-boson model, see Figs. 1 \- 3.

In a second time, we focus on the crucial task of exploring and understanding
their respective range of validity. To achieve that, we use one approach to
benchmark the other. We establish and test successfully a validity criterion
(41) for the perturbative expansion depending only on the inverse bath
temperature $\beta$ and on the reorganization energy $Q$ (11).

Regarding the reaction coordinate mapping and its approximate steady state
(27), we quantify its performance and derived a validity criterion (44)
involving only the inverse bath temperature $\beta$ and the reaction
coordinate frequency $\Omega$, holding for arbitrary spectral width
$\gamma_{UD}$ and arbitrary coupling strength. This criterion relies on
analytical arguments which were confirmed numerically.

Thanks to these validity criteria, one has in hand practical tools to assess
the validity range of these two techniques. Although these validity criteria
were numerically tested for the spin-boson model, they can be extended to
arbitrary systems, and the curious fact that they do not involve the system’s
energy scale might suggest that they do work for other systems. It would be
interesting to test that.

Additionally, it would also be instructive and useful to extend this
comparative analysis to other techniques like pseudo-mode [34, 35, 36, 37], as
well as to the ultrastrong coupling regime [47, 58, 52].

###### Acknowledgements.

I am grateful for on going discussions with Ilya Sinayskiy, Graeme Pleasance,
and Francesco Petruccione. I also would like to thank Patrice Camati for a
crash course on QuTiP, as well as all QuTiP contributors for setting up and
developing such a useful tool. The author acknowledges the support of the
French Agence Nationale de la Recherche (ANR), under grant ANR-20-ERC9-0010
(project QSTEAM).

##  Appendix A Expression of the function $g(\nu,\nu^{\prime})$

The function $g(\nu,\nu^{\prime})$ introduced in the main text is defined by
$g(\nu,\nu^{\prime}):=\int_{0}^{\beta}du_{1}\int_{0}^{u_{1}}du_{2}e^{-\nu
u_{1}+\nu^{\prime}u_{2}}c_{B}(u_{1}-u_{2})$. Re-writing its expression and
introducing the variable $v_{2}=u_{1}-u_{2}$, we obtain

$\displaystyle g(\nu,\nu^{\prime})$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}du_{1}\int_{0}^{u_{1}}du_{2}e^{-\nu
u_{1}+\nu^{\prime}u_{2}}c_{B}(u_{1}-u_{2})$ (45) $\displaystyle=$
$\displaystyle\int_{0}^{\beta}du_{1}\int_{0}^{u_{1}}dv_{2}e^{-\nu
u_{1}}e^{\nu^{\prime}(u_{1}-v_{2})}c_{B}(v_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}du_{1}\int_{0}^{u_{1}}dv_{2}e^{(\nu^{\prime}-\nu)u_{1}}e^{-\nu^{\prime}v_{2}}c_{B}(v_{2})$
$\displaystyle=$
$\displaystyle\int_{0}^{\beta}dv_{2}\int_{v_{2}}^{\beta}du_{1}e^{(\nu^{\prime}-\nu)u_{1}}e^{-\nu^{\prime}v_{2}}c_{B}(v_{2})$
$\displaystyle=$
$\displaystyle\int_{0}^{\beta}dv\frac{e^{(\nu^{\prime}-\nu)\beta}-e^{(\nu^{\prime}-\nu)v}}{\nu^{\prime}-\nu}e^{-\nu^{\prime}v}c_{B}(v)$
$\displaystyle=$
$\displaystyle\beta\int_{0}^{1}dv\frac{e^{(\nu^{\prime}-\nu)\beta}-e^{(\nu^{\prime}-\nu)v\beta}}{\nu^{\prime}-\nu}e^{-\nu^{\prime}v\beta}c_{B}(v\beta)$
$\displaystyle=$
$\displaystyle\frac{\beta}{\nu^{\prime}-\nu}\int_{0}^{1}dv[e^{(\nu^{\prime}-\nu)\beta}e^{-\nu^{\prime}v\beta}-e^{-\nu
v\beta}]c_{B}(v\beta)$ $\displaystyle=$
$\displaystyle\frac{\beta}{\nu^{\prime}-\nu}\Bigg{[}e^{(\nu^{\prime}-\nu)\beta}\int_{0}^{1}due^{-\nu^{\prime}\beta
u}c_{B}(u\beta)$ $\displaystyle\hskip 39.83368pt-\int_{0}^{1}due^{-\nu\beta
u}c_{B}(u\beta)\Bigg{]}.$

We are then led to compute

$G(\nu,\beta):=\int_{0}^{1}due^{-\nu\beta u}c_{B}(u\beta).$ (46)

In order to obtain analytical expressions for the under-damped and over-damped
bath spectral density, it will be convenient to decompose $G(\nu,\beta)$ in
the following way,

$G(\nu,\beta)=C(-\nu)+e^{-\beta\nu}C(\nu),$ (47)

where

$\displaystyle C(\nu)$ $\displaystyle:=$
$\displaystyle\int_{0}^{\infty}d\omega
J(\omega)\left[\frac{n_{\omega}+1}{\beta(\omega-\nu)}-\frac{n_{\omega}}{\beta(\omega+\nu)}\right]$
(48) $\displaystyle=$ $\displaystyle\int_{0}^{\infty}d\omega
J(\omega)\frac{\nu{\rm
coth}(\omega\beta/2)+\omega}{\beta(\omega^{2}-\nu^{2})}.$ (49)

Then, from (45), we obtain

$\displaystyle
g(\nu,\nu^{\prime})=\frac{\beta}{\nu^{\prime}-\nu}\left[e^{(\nu^{\prime}-\nu)\beta}G(\nu^{\prime},\beta)-G(\nu,\beta)\right],$
(50)

for $\nu\neq\nu^{\prime}$, and for $\nu^{\prime}=\nu$,

$g(\nu):=g(\nu,\nu)=\beta\left[\beta
G(\nu,\beta)+G^{\prime}(\nu,\beta)\right]$ (51)

where
$G^{\prime}(\nu,\beta):=\frac{\partial}{\partial\nu}G(\nu,\beta)=-C^{\prime}(-\nu)-\beta
e^{-\nu\beta}C(\nu)+e^{-\nu\beta}C^{\prime}(\nu)$ and $C^{\prime}(\nu)$ is the
partial derivative with respect to $\nu$,

$\displaystyle C^{\prime}(\nu)$ $\displaystyle:=$ $\displaystyle\frac{\partial
C(\nu)}{\partial\nu}$ $\displaystyle=$ $\displaystyle\int_{0}^{\infty}d\omega
J(\omega)\left[\frac{n_{\omega}+1}{\beta(\omega-\nu)^{2}}+\frac{n_{\omega}}{\beta(\omega+\nu)^{2}}\right].$

Alternatively, in term of the function $C(\nu)$, we have
$g(\nu)=\beta\left[\beta
C(-\nu)-C^{\prime}(-\nu)+e^{-\beta\nu}C^{\prime}(\nu)\right]$.

###  A.1 Exact expression of $C(\omega_{S})$ and $C^{\prime}(\omega_{S})$

####  A.1.1  Over-damped (Lorentz-Drude) spectral density

For
$J_{OD}(\omega)=\alpha\frac{\omega\omega_{c}^{2}}{\omega^{2}+\omega_{c}^{2}}$,
we have

$\displaystyle C_{OD}(\nu)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega J(\omega)\frac{\nu{\rm
coth}(\omega\beta/2)+\omega}{\beta(\omega^{2}-\nu^{2})}$ $\displaystyle=$
$\displaystyle\frac{\alpha\omega_{c}^{2}}{\beta}\int_{0}^{\infty}d\omega\frac{\omega}{\omega^{2}+\omega_{c}^{2}}\frac{\nu\coth{\omega\beta/2}+\omega}{\omega^{2}-\nu^{2}}$
$\displaystyle=$
$\displaystyle\frac{\alpha\omega_{c}^{2}}{\beta}\Bigg{[}\int_{0}^{\infty}d\omega\frac{\omega}{\omega^{2}+\omega_{c}^{2}}\frac{\omega}{\omega^{2}-\nu^{2}}$
$\displaystyle+\frac{2\nu}{\beta}\sum_{n=-\infty}^{+\infty}\int_{0}^{\infty}d\omega\frac{\omega}{\omega_{c}^{2}+\omega^{2}}\frac{1}{\omega^{2}-\nu^{2}}\frac{\omega}{\omega^{2}+\nu_{n}^{2}}\Bigg{]}$

with $\nu_{n}=2\pi n/\beta$, called the Matsubara frequencies [1]. For the
first term we have

$\displaystyle\int_{0}^{\infty}d\omega\frac{\omega}{\omega^{2}+\omega_{c}^{2}}\frac{\omega}{\omega^{2}-\nu^{2}}$
$\displaystyle=\frac{1}{\omega_{c}^{2}+\nu^{2}}\int_{0}^{\infty}d\omega\left(\frac{\omega_{c}^{2}}{\omega^{2}+\omega_{c}^{2}}+\frac{\nu^{2}}{\omega^{2}-\nu^{2}}\right)$
$\displaystyle=\frac{\omega_{c}^{2}}{\omega_{c}^{2}+\nu^{2}}\frac{\pi}{2\omega_{c}}.$
(54)

Generalising that to situations where $\omega_{c}$ is a complex number (which
will be useful for under-damped spectral densities, see in the following), we
have

$\displaystyle\int_{0}^{\infty}d\omega\frac{\omega}{\omega^{2}+\omega_{c}^{2}}\frac{\omega}{\omega^{2}-\nu^{2}}=\begin{cases}\frac{\omega_{c}^{2}}{\omega_{c}^{2}+\nu^{2}}\frac{\pi}{2\omega_{c}}&\text{if
$\Re\omega_{c}>0$},\\\
-\frac{\omega_{c}^{2}}{\omega_{c}^{2}+\nu^{2}}\frac{\pi}{2\omega_{c}}&\text{if
$\Re\omega_{c}<0$}.\end{cases}$

For the second term, we have

$\displaystyle\sum_{n=-\infty}^{\infty}\int_{0}^{\infty}d\omega\frac{\omega}{\omega_{c}^{2}+\omega^{2}}\frac{1}{\omega^{2}-\nu^{2}}\frac{\omega}{\omega^{2}+\nu_{n}^{2}}$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{1}{\nu^{2}+\nu_{n}^{2}}\int_{0}^{\infty}d\omega\left[\frac{\nu_{n}^{2}}{\omega_{c}^{2}-\nu_{n}^{2}}\left(\frac{1}{\omega^{2}+\nu_{n}^{2}}-\frac{1}{\omega^{2}+\omega_{c}^{2}}\right)+\frac{\nu^{2}}{\nu^{2}+\omega_{c}^{2}}\left(\frac{1}{\omega^{2}-\nu^{2}}-\frac{1}{\omega^{2}+\omega_{c}^{2}}\right)\right]$
$\displaystyle=\sum_{n=-\infty}^{\infty}\frac{1}{\nu^{2}+\nu_{n}^{2}}\left[\frac{\nu_{n}^{2}}{\omega_{c}^{2}-\nu_{n}^{2}}\left(\frac{\pi}{2|\nu_{n}|}-\frac{\pi}{2(\pm\omega_{c})}\right)+\frac{\nu^{2}}{\nu^{2}+\omega_{c}^{2}}\left(0-\frac{\pi}{2(\pm\omega_{c})}\right)\right]$
$\displaystyle=\frac{\pi}{2}\sum_{n=-\infty}^{\infty}\frac{1}{\nu^{2}+\nu_{n}^{2}}\frac{1}{\pm\omega_{c}}\left[\frac{|\nu_{n}|}{\pm\omega_{c}+|\nu_{n}|}-\frac{\nu^{2}}{\nu^{2}+\omega_{c}^{2}}\right]$
$\displaystyle=\frac{\pi}{\pm
2\omega_{c}}\left[\sum_{n=-\infty}^{\infty}\frac{1}{\nu^{2}+\nu_{n}^{2}}\frac{|\nu_{n}|}{\pm\omega_{c}+|\nu_{n}|}-\frac{\nu^{2}}{\nu^{2}+\omega_{c}^{2}}\sum_{n=-\infty}^{\infty}\frac{1}{\nu^{2}+\nu_{n}^{2}}\right]$
$\displaystyle=\frac{\pi}{\pm
2\omega_{c}}\left[2\sum_{n=1}^{\infty}\frac{1}{\nu^{2}+\nu_{n}^{2}}\frac{\nu_{n}}{\pm\omega_{c}+\nu_{n}}-\frac{1}{\nu^{2}+\omega_{c}^{2}}-\frac{2\nu^{2}}{\nu^{2}+\omega_{c}^{2}}\sum_{n=1}^{\infty}\frac{1}{\nu^{2}+\nu_{n}^{2}}\right]$
$\displaystyle=-\frac{\pi}{\pm
2\omega_{c}}\frac{1}{\nu^{2}+\omega_{c}^{2}}+\frac{\pi}{\pm
2\omega_{c}}\left[2\frac{\beta^{2}}{4\pi^{2}}\sum_{n=1}^{\infty}\frac{1}{n^{2}+\frac{\nu^{2}\beta^{2}}{4\pi^{2}}}\frac{n}{n\pm\frac{\omega_{c}\beta}{2\pi}}-\frac{2\nu^{2}}{\nu^{2}+\omega_{c}^{2}}\frac{\beta^{2}}{4\pi^{2}}\sum_{n=1}^{\infty}\frac{1}{n^{2}+\frac{\nu^{2}\beta^{2}}{4\pi^{2}}}\right]$
$\displaystyle=-\frac{\pi}{\pm
2\omega_{c}}\frac{1}{\nu^{2}+\omega_{c}^{2}}+\frac{\pi}{\pm\omega_{c}}\frac{\beta^{2}}{4\pi^{2}}\left[\frac{2\pi^{2}}{\beta^{2}}\frac{1}{\nu^{2}+\omega_{c}^{2}}F(\nu\beta,\pm\omega_{c}\beta)-\frac{\nu^{2}}{\nu^{2}+\omega_{c}^{2}}\frac{\frac{\nu\beta}{2}\coth\frac{\nu\beta}{2}-1}{\frac{\nu^{2}\beta^{2}}{2\pi^{2}}}\right]$
$\displaystyle=-\frac{\pi}{\pm
2\omega_{c}}\frac{1}{\nu^{2}+\omega_{c}^{2}}+\frac{\pi}{\pm
2\omega_{c}}\frac{1}{\nu^{2}+\omega_{c}^{2}}\left[F(\nu\beta,\pm\omega_{c}\beta)-\frac{\nu\beta}{2}\coth\frac{\nu\beta}{2}+1\right]$
$\displaystyle=\frac{\pi}{\pm
2\omega_{c}}\frac{1}{\nu^{2}+\omega_{c}^{2}}\left[F(\nu\beta,\pm\omega_{c}\beta)-\frac{\nu\beta}{2}\coth\frac{\nu\beta}{2}\right],$

where $\pm\omega_{c}$ stands for the possibility of $\omega_{c}$ being
complex, where in such case one has to choose the sign corresponding to
$\Re(\pm\omega_{c})>0$, and

$F(\nu\beta,\pm\omega_{c}\beta)=\frac{1}{2\pi}\left[-(\pm\omega_{c}+i\nu)\beta\Psi\left(1+i\frac{\nu\beta}{2\pi}\right)-(\pm\omega_{c}-i\nu)\beta\Psi\left(1-i\frac{\nu\beta}{2\pi}\right)\pm
2\omega_{c}\beta\Psi\left(1+\frac{\pm\omega_{c}\beta}{2\pi}\right)\right],$
(57)

with $\Psi(x)$ being the Digamma function. All together we obtain (changing
the variable in the argument from $\nu$ to $\omega$),

$C_{OD}(\omega)=\frac{\pi}{\pm
2\omega_{c}}\frac{\alpha}{\beta}\frac{\omega_{c}^{2}}{\omega_{c}^{2}+\omega^{2}}\left[\omega_{c}^{2}-\omega^{2}{\rm
coth}(\omega\beta/2)+\frac{2\omega}{\beta}F(\omega\beta,\pm\omega_{c}\beta)\right].$
(58)

For
$C(0)=\int_{0}^{\infty}d\omega\frac{J(\omega)}{\beta\omega}=\frac{Q}{\beta}$,
we simply have

$C_{OD}(0)=\frac{Q_{OD}}{\beta}=\frac{\pi}{2}\frac{\alpha(\pm\omega_{c})}{\beta}.$
(59)

Then, $C_{OD}^{\prime}(\omega)$ is ”just” the derivative of $C_{OD}(\omega)$,
which gives

$C_{OD}^{\prime}(\omega)=\frac{\pi}{2}\frac{\alpha}{\beta}\frac{\omega_{c}^{2}}{\omega_{c}^{2}+\omega^{2}}\left[\frac{-2\omega\omega_{c}}{\omega_{c}^{2}+\omega^{2}}[1+{\rm
coth}(\omega\beta/2)]-\frac{\omega^{2}\beta}{2\omega_{c}}{\rm
coth^{\prime}}(\omega\beta/2)+\frac{2}{\omega_{c}\beta}\frac{\omega_{c}^{2}-\omega^{2}}{\omega_{c}^{2}+\omega^{2}}F(\omega\beta,\omega_{c}\beta)+\frac{2\omega}{\omega_{c}\beta}F^{\prime}(\omega\beta,\omega_{c}\beta)\right],$
(60)

where ${\rm coth^{\prime}}(x):=\frac{e^{x}}{{\rm sinh}x}(1-{\rm coth}x)$ is
simply the derivative of ${\rm coth}(x)$, and

$\displaystyle F^{\prime}(\omega\beta,\omega_{c}\beta)$ $\displaystyle:=$
$\displaystyle\frac{\partial}{\partial\omega}F(\omega\beta,\omega_{c}\beta)$
$\displaystyle=$
$\displaystyle\frac{\beta}{2\pi}\left[-i\Psi\left(1+i\frac{\omega\beta}{2\pi}\right)+i\Psi\left(1-i\frac{\omega\beta}{2\pi}\right)+(\omega-i\omega_{c})\frac{\beta}{2\pi}\Psi^{\prime}\left(1+i\frac{\omega\beta}{2\pi}\right)+(\omega+i\omega_{c})\frac{\beta}{2\pi}\Psi^{\prime}\left(1-i\frac{\omega\beta}{2\pi}\right)\right]$

with $\Psi^{\prime}$ is the derivative of the Digamma function.

Again, if $\omega_{c}$ is complex, we simply have

$C_{OD}^{\prime}(\omega)=\frac{\pi}{\pm
2\omega_{c}}\frac{\alpha}{\beta}\frac{\omega_{c}^{2}}{\omega_{c}^{2}+\omega^{2}}\left[\frac{-2\omega\omega_{c}^{2}}{\omega_{c}^{2}+\omega^{2}}[1+{\rm
coth}(\omega\beta/2)]-\frac{\omega^{2}\beta}{2}{\rm
coth^{\prime}}(\omega\beta/2)+\frac{2}{\beta}\frac{\omega_{c}^{2}-\omega^{2}}{\omega_{c}^{2}+\omega^{2}}F(\omega\beta,\pm\omega_{c}\beta)+\frac{2\omega}{\beta}F^{\prime}(\omega\beta,\pm\omega_{c}\beta)\right].$
(62)

####  A.1.2 Under-damped spectral density

We now consider an under-damped spectral density of the form $J_{UD}(\omega)$
(22),

$J_{UD}(\omega):=\omega\frac{2}{\pi}\frac{\gamma_{UD}\Omega^{2}\lambda^{2}}{(\Omega^{2}-\omega^{2})^{2}+(\gamma_{UD}\Omega\omega)^{2}}.$
(63)

Such spectral densities can lead to difficulties related to analytical
integration. This can be circumvented by expressing $J_{UD}$ as the difference
of two over-damped spectral densities,

$J_{UD}(\omega)=J_{OD}^{-}(\omega)-J_{OD}^{+}(\omega)$ (64)

where

$J_{OD}^{\pm}(\omega):=\frac{2}{\pi}\frac{\gamma_{UD}\Omega^{2}\lambda^{2}}{\omega^{2}_{+}-\omega^{2}_{-}}\frac{\omega}{\omega^{2}+\omega_{\pm}^{2}},$
(65)

with
$\omega^{2}_{\pm}:=\Omega^{2}\left(\frac{\gamma_{UD}^{2}}{2}-1\pm\gamma_{UD}\sqrt{\frac{\gamma_{UD}^{2}}{4}-1}\right)$,
always positive for $\gamma_{UD}\geq 2$, and complex for $\gamma_{UD}<2$.
Using this mapping, we straightforwardly obtain

$\displaystyle C_{UD}(\omega)=C_{OD}^{-}(\omega)-C_{OD}^{+}(\omega),$
$\displaystyle
C^{\prime}_{UD}(\omega)={C^{\prime}}_{OD}^{-}(\omega)-{C^{\prime}}_{OD}^{+}(\omega),$
(66)

where $C_{OD}^{\pm}(\omega)$ and ${C^{\prime}}_{OD}^{\pm}(\omega)$ are given
by the above expressions (58) and (60) substituting $\omega_{c}^{2}$ by
$\omega_{\pm}^{2}$, and $\alpha$ by
$\alpha_{\pm}:=\frac{2}{\pi}\frac{\gamma_{UD}\Omega^{2}\lambda^{2}/\omega_{\pm}^{2}}{\omega^{2}_{+}-\omega^{2}_{-}}=\frac{\lambda^{2}}{\pi\Omega^{2}}f_{\pm}(\gamma_{UD})$
with

$f_{\pm}(\gamma_{UD})=\frac{1}{\left(\frac{\gamma_{UD}^{2}}{2}-1\pm\gamma_{UD}\sqrt{\frac{\gamma_{UD}^{2}}{4}-1}\right)\sqrt{\frac{\gamma_{UD}^{2}}{4}-1}}.$
(67)

However, one has to be careful (see (A.1.1)) for $\gamma_{UD}<2$ since
$\omega_{\pm}^{2}$ becomes complex. One can verifies that for all
$\gamma_{UD}>0$, $\Re\sqrt{\omega_{\pm}^{2}}>0$, so that the above expressions
(A.1.2) still hold without change of signs, namely,

$\displaystyle C_{UD}(\omega)$ $\displaystyle=$
$\displaystyle\frac{\pi}{2\omega_{-}}\frac{\alpha_{-}}{\beta}\frac{\omega_{-}^{2}}{\omega_{-}^{2}+\omega^{2}}\Big{[}\omega_{-}^{2}-\omega^{2}{\rm
coth}(\omega\beta/2)$ (68) $\displaystyle\hskip
73.97733pt+\frac{2\omega}{\beta}F(\omega\beta,\omega_{-}\beta)\Big{]}$
$\displaystyle-\frac{\pi}{2\omega_{+}}\frac{\alpha_{+}}{\beta}\frac{\omega_{+}^{2}}{\omega_{+}^{2}+\omega^{2}}\Big{[}\omega_{+}^{2}-\omega^{2}{\rm
coth}(\omega\beta/2)$ $\displaystyle\hskip
73.97733pt+\frac{2\omega}{\beta}F(\omega\beta,\omega_{+}\beta)\Big{]},$

and similarly for $C^{\prime}_{UD}(\omega)$.

As a side note, we also show that the reorganisation energy is given by
$Q_{UD}=\frac{\lambda^{2}}{\Omega}$. From the above mapping into over-damped
spectral densities, we have

$\displaystyle Q_{UD}$ $\displaystyle=$ $\displaystyle Q_{OD,-}-Q_{OD,+}$
$\displaystyle=$
$\displaystyle\frac{\pi}{2}\alpha_{-}\sqrt{\omega_{-}^{2}}-\frac{\pi}{2}\alpha_{+}\sqrt{\omega_{+}^{2}}$
$\displaystyle=$
$\displaystyle\frac{\gamma_{UD}\Omega^{2}\lambda^{2}}{\omega_{+}^{2}-\omega_{-}^{2}}\left(\frac{1}{\omega_{-}}-\frac{1}{\omega_{+}}\right).$

Since,

$\displaystyle\frac{1}{\omega^{2}_{+}-\omega^{2}_{-}}\left(\frac{1}{\omega_{-}}-\frac{1}{\omega_{+}}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{\omega^{2}_{+}-\omega^{2}_{-}}\frac{\omega_{+}-\omega_{-}}{\omega_{+}\omega_{-}}$
(70) $\displaystyle=$
$\displaystyle\frac{1}{\omega_{+}+\omega_{-}}\frac{1}{\omega_{+}\omega_{-}}$
$\displaystyle=$
$\displaystyle\frac{1}{\omega_{+}+\omega_{-}}\frac{1}{\Omega^{2}}$
$\displaystyle=\frac{1}{\Omega^{3}}\Bigg{[}\left(\frac{\gamma_{UD}^{2}}{2}-1+\gamma_{UD}\sqrt{\frac{\gamma_{UD}^{2}}{4}-1}\right)^{1/2}$
$\displaystyle+\left(\frac{\gamma_{UD}^{2}}{2}-1-\gamma_{UD}\sqrt{\frac{\gamma_{UD}^{2}}{4}-1}\right)^{1/2}\Bigg{]}^{-1}$
$\displaystyle=\frac{1}{\Omega^{3}}\frac{1}{\gamma_{UD}}.$

Note that one can easily see the last line by taking the square of what is in
the square bracket. Then, we finally obtain

$Q_{UD}=\frac{\lambda^{2}}{\Omega},$ (71)

Finally, note some useful identities with $\omega_{\pm}^{2}$,

$\displaystyle\omega_{+}\omega_{-}=\Omega^{2}$
$\displaystyle\omega_{+}+\omega_{-}=\Omega\gamma_{UD}.$ (72)

##  Appendix B Some additional plots of $d_{\rm coh/pop}^{\rm ext}$ and
criteria ${\rm cr}_{i}$

In this section we provide some additional plots in Figs. 9 and 10, showing
unambiguously that the criterion ${\rm cr}_{4}=\beta Q$ (as well as ${\rm
cr}_{1}$) accurately predicts the validity range of the perturbative
expansion.

(a)
(b)
(c)
(d)

Figure 9: On the left-hand-side, plots of $d_{\rm pop}^{\rm exp}$ (large
purple dots) and $d_{\rm coh}^{\rm exp}$ (small green dots); on the right-
hand-side plots of the criteria ${\rm cr}_{i}$, all in function of $\lambda$
(unit of $\omega_{S}$). The color convention for the criteria is the same as
in the main text, namely, ${\rm cr}_{1}$ (black solid line), ${\rm cr}_{2}$
(blue solid line), ${\rm cr}_{3}$ (green solid line), and ${\rm cr}_{4}$ (red
dashed line). The other parameters are as follows:(a) $\Omega/\omega_{S}=1$,
$\omega_{S}\beta=0.5$; (b) $\Omega/\omega_{S}=0.1$, $\omega_{S}\beta=0.5$; (c)
$\Omega/\omega_{S}=1$, $\omega_{S}\beta=20$; (d)$\Omega/\omega_{S}=0.1$,
$\omega_{S}\beta=20$. The remainder of the parameters are chosen as in figures
of the main text, namely $\gamma_{UD}=0.1$, $r_{z}=\sqrt{0.75}$,
$r=r_{x}+ir_{y}=0.5$.

(a)
(b)
(c)
(d)
(e)

Figure 10: On the left-hand-side, plots of $d_{\rm pop}^{\rm exp}$ (large
purple dots) and $d_{\rm coh}^{\rm exp}$ (small green dots); on the right-
hand-side plots of the criteria ${\rm cr}_{i}$, all in function of $\beta$
(unit of $\omega_{S}^{-1}$). The color convention for the criteria is the same
as in the main text, namely, ${\rm cr}_{1}$ (black solid line), ${\rm cr}_{2}$
(blue solid line), ${\rm cr}_{3}$ (green solid line), and ${\rm cr}_{4}$ (red
dashed line). The other parameters are as follows:(a) $\Omega/\omega_{S}=10$,
$\lambda/\omega_{S}=0.1$; (b) $\Omega/\omega_{S}=1$, $\lambda/\omega_{S}=0.1$;
(c) $\Omega/\omega_{S}=0.1$, $\lambda/\omega_{S}=0.1$;
(d)$\Omega/\omega_{S}=10$, $\lambda/\omega_{S}=1.5$; (e)$\Omega/\omega_{S}=1$,
$\lambda/\omega_{S}=1.5$. The remainder of the parameters are chosen as in
figures of the main text, namely $\gamma_{UD}=0.1$, $r_{z}=\sqrt{0.75}$,
$r=r_{x}+ir_{y}=0.5$.

## References

  * [1] H. Breuer and F. Petruccione, Theory of Open Quantum Systems, (Oxford, Oxford, 2002).
  * [2] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, Processus d’interaction entre photons et atomes, (EDP Science/CNRS Éditions, Paris, 2001).
  * [3] A. W. Chin, S. F. Huelga and M. B. Plenio, Coherence and decoherence in biological systems: principles of noise-assisted transport and the origin of long-lived coherences, Phil. Trans. R. Soc. A (2012) 370, 3638–3657.
  * [4] P. Ribeiro and V. R. Vieira, Non-Markovian effects in electronic and spin transport, Phys. Rev. B 92, 100302(R) (2015).
  * [5] P. Strasberg, G. Schaller, T. L. Schmidt, and M. Esposito, Phys. Rev. B 97, 205405 (2018).
  * [6] L. A. Correa, B. Xu, B. Morris, and G. Adesso, Pushing the limits of the reaction-coordinate mapping, J. Chem. Phys. 151, 094107 (2019).
  * [7] S. V. Moreira, B. Marques, R. R. Paiva, L. S. Cruz, D. O. Soares-Pinto, and F. L. Semião, Phys. Rev. A 101, 012123 (2020).
  * [8] E. Zerah-Harush and Y. Dubi, Phys. Rev. Research 2, 023294 (2020).
  * [9] D. Dwiputra and F. P. Zen, Phys. Rev. A 104, 022205 (2021).
  * [10] Nicholas Anto-Sztrikacs and Dvira Segal, Strong coupling effects in quantum thermal transport with the reaction coordinate method, New J. Phys. 23, 063036 (2021).
  * [11] A. Trushechkin, Calculation of coherences in Forster and modified Redfield theories of excitation energy transfer, J. Chem. Phys. 151, 074101 (2019).
  * [12] L. A. Correa, M. Perarnau-Llobet, K. V. Hovhannisyan, S. Hernández-Santana, M. Mehboudi, and A. Sanpera, Phys. Rev. A 96, 062103 (2017).
  * [13] M. Mehboudi, A. Lampo, C. Charalambous, L. A. Correa, M. A. García-March and M. Lewenstein, Phys. Rev. Lett. 122, 030403 (2019).
  * [14] M Salado-Mejía, R. Román-Ancheyta, F. Soto-Eguibar, and H. M. Moya-Cessa, Spectroscopy and critical quantum thermometry in the ultrastrong coupling regime, Quantum Sci. Technol. 6, 025010 (2021).
  * [15] D. Gelbwaser-Klimovsky and A. Aspuru-Guzik, Strongly Coupled Quantum Heat Machines, J. Phys. Chem. Lett. 6, 3477 (2015).
  * [16] P. Strasberg, G. Schaller, N. Lambert, and T. Brandes, Nonequilibrium thermodynamics in the strong coupling and non-Markovian regime based on a reaction coordinate mapping, New J. Phys. 18, 073007 (2016).
  * [17] D. Newman, F. Mintert, and A. Nazir, Phys. Rev. E 95, 032139 (2017).
  * [18] M. Perarnau-Llobet, H. Wilming, A. Riera, R. Gallego, and J. Eisert, Phys. Rev. Lett. 120, 120602 (2018).
  * [19] M.Wertnik, A. Chin, F. Nori, and N. Lambert, Optimizing co-operative multi-environment dynamics in a dark-state-enhanced photosynthetic heat engine, J. Chem. Phys. 149, 084112 (2018).
  * [20] D. Newman, F. Mintert, and A. Nazir, Phys. Rev. E 101, 052129 (2020).
  * [21] M Wiedmann, J T Stockburger, and J Ankerhold, Non-Markovian dynamics of a quantum heat engine: out-of-equilibrium operation and thermal coupling control, New J. Phys. 22, 033007 (2020).
  * [22] J. Anders, C.R.J. Sait, and S.A.R. Horsley, Quantum Brownian Motion for Magnets, arXiv:2009.00600
  * [23] A. Kolli, E. J. O’Reilly, G. D. Scholes, and A. Olaya-Castro, The fundamental role of quantized vibrations in coherent light harvesting by cryptophyte algae, J. Chem. Phys. 137, 174109 (2012).
  * [24] N. Lambert, Y.-N. Chen, Y.-C. Cheng, C.-M. Li, G.-Y. Chen and Franco Nori, Nat. Physics 9, Quantum biology, 10-18 (2013).
  * [25] G. Scholes, G. Fleming, L. Chen, et al., Using coherence to enhance function in chemical and biophysical systems, Nature 543, 647–656 (2017).
  * [26] N. Lambert, T. Raheja, S. Ahmed, A. Pitchford, F. Nori, BoFiN-HEOM: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics, arXiv:2010.10806.
  * [27] D. P S McCutcheon and A. Nazir, Quantum dot Rabi rotations beyond the weak exciton–phonon coupling regime, New Journal of Physics 12, 113042 (2012).
  * [28] D. P S McCutcheon and A. Nazir, Coherent and incoherent dynamics in excitonic energy transfer: correlated fluctuations and off-resonance effects, Phys. Rev. B 83, 165101 (2011).
  * [29] C. Roy and S. Hughes, Phonon-Dressed Mollow Triplet in the Regime of Cavity Quantum Electrodynamics: Excitation-Induced Dephasing and Nonperturbative Cavity Feeding Effects, Phys. Rev. Lett. 106, 247403 (2011).
  * [30] C. Kong Lee, J. Moix, J. Cao, Accuracy of Second Order Perturbation theory in the Polaron and Variational Polaron Frames, J. Chem. Phys. 136, 204120 (2012)
  * [31] A. Nazir and D. P S McCutcheon, Modelling exciton–phonon interactions in optically driven quantum dots, J. Phys.: Condens. Matter 28, 103002 (2016).
  * [32] Jake Iles-Smith, Neill Lambert, and Ahsan Nazir, Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems, Phys. Rev. A 90, 032114 (2014).
  * [33] Jake Iles-Smith, Arend G. Dijkstra, Neill Lambert, and Ahsan Nazir, Energy transfer in structured and unstructured environments: Master equations beyond the Born-Markov approximations, J. Chem. Phys. 144, 044110 (2016).
  * [34] B. M. Garraway, Nonperturbative decay of an atomic system in a cavity, Phys. Rev. A. 55, 2290 (1997).
  * [35] G. Pleasance and B. M. Garraway, Application of quantum Darwinism to a structured environment, Phys. Rev. A 96, 062105 (2017).
  * [36] A.E. Teretenkov, Pseudomode Approach and Vibronic Non-Markovian Phenomena in Light-Harvesting Complexes, Proc. Steklov Inst. Math. 306, 242–256 (2019).
  * [37] Graeme Pleasance, Barry M. Garraway, and Francesco Petruccione, Generalized theory of pseudomodes for exact descriptions of non-Markovian quantum processes, Phys. Rev. R. 2, 043058 (2020).
  * [38] V. Bach, J. Fröhlich, and I. M. Sigal, Return to Equilibrium, Journal of Mathematical Physics 41, 3985 (2000).
  * [39] J. Frohlich, Marco Merkli, Another Return of “Return to Equilibrium”, Commun. Math. Phys. 251, 235–262 (2004).
  * [40] M. Merkli, I. M. Sigal, and G. P. Berman, Decoherence and Thermalization, Phys. Rev. Lett. 98, 130401 (2007).
  * [41] M. Könenberg and M. Merkli, On the irreversible dynamics emerging from quantum resonances, J. Math. Phys. 57, 033302 (2016).
  * [42] Marco Merkli, Quantum Markovian master equations: Resonance theory shows validity for all time scales, Ann. Phys. 412, 16799 (2020).
  * [43] T. Mori and S. Miyashita, Dynamics of the Density Matrix in Contact with a Thermal Bath and the Quantum Master Equation, Journal of the Physical Society of Japan 7, 124005 (2008).
  * [44] A. S. Trushechkin, M. Merkli, J. D. Cresser, J. Anders, Open quantum system dynamics and the mean force Gibbs state, arXiv:2110.01671
  * [45] Y. Subaşı, C. H. Fleming, J. M. Taylor, and B. L. Hu, Phys. Rev. E 86, 061132 (2012).
  * [46] Archak Purkayastha, Giacomo Guarnieri, Mark T. Mitchison, Radim Filip and John Goold, Tunable phonon-induced steady-state coherence in a double-quantum-dot charge qubit, npj Quantum Information (2020) 6:27.
  * [47] J. D. Cresser, and J. Anders, Weak and ultrastrong coupling limits of the quantum mean force Gibbs state, Phys. Rev. Lett. 127, 250601 (2021).
  * [48] Richard P. Feynman, Phys. Rev. 84, 108-128 (1951).
  * [49] J. Wu, F. Liu, Y. Shen, J. Cao, and R. J. Silbey, Efficient energy transfer in light-harvesting systems, I: optimal temperature, reorganization energy and spatial–temporal correlations, New J. Phys. 12, 105012 (2010).
  * [50] G. Ritschel, J. Roden, W. T. Strunz, and A. Eisfeld, An efficient method to calculate excitation energy transfer in light-harvesting systems: application to the Fenna–Matthews–Olson complex, New J. Phys. 13, 113034 (2011).
  * [51] V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems (Wiley-VCH, Weinheim, 2011).
  * [52] C. L. Latune, Steady states in ultrastrong coupling regime: perturbative expansion and first orders, arXiv:2110.02186
  * [53] A. Garg, J. N. Onuchic, and V. Ambegaokar, Effect of friction on electron transfer in biomolecules, J. Chem. Phys. 83, 4491 (1985).
  * [54] A. Nazir, G. Schaller, The Reaction Coordinate Mapping in Quantum Thermodynamics (2018). In: Binder F., Correa L., Gogolin C., Anders J., Adesso G. (eds) Thermodynamics in the Quantum Regime. Fundamental Theories of Physics, vol 195. Springer, Cham.
  * [55] J. Huh, S. Mostame, T. Fujita, M.-H. Yung and A. Aspuru-Guzik, Linear-algebraic bath transformation for simulating complex open quantum systems, New J. Phys. 16, 123008 (2014).
  * [56] E. K. Irish, J. Gea-Banacloche, I. Martin, and K. C. Schwab, Dynamics of a two-level system strongly coupled to a high-frequency quantum oscillator, Phys. Rev. B 72, 195410 (2005).
  * [57] Zi-Min Li, Murray T. Batchelor, Phys. Rev. A 104, 033712 (2021).
  * [58] A. Trushechkin, Quantum master equations and steady states for the ultrastrong-coupling limit and the strong-decoherence limit, arXiv:2109.01888