# Steady state in ultrastrong coupling regime: perturbative expansion and
first orders

Camille L Latune  Quantum Research Group, School of Chemistry and Physics,
University of KwaZulu-Natal, Durban, KwaZulu-Natal, 4001, South Africa
National Institute for Theoretical Physics (NITheP), KwaZulu-Natal, 4001,
South Africa

###### Abstract

Understanding better the dynamics and steady states of systems strongly
coupled to thermal baths is a great theoretical challenge with promising
applications in several fields of quantum technologies. Among several
interesting strategies to gain access to the steady state, one consists in
obtaining approximate expressions of the mean force Gibbs state, the reduced
state of the global system-bath steady state. Here, we present explicit
analytical expressions of corrective terms to the ultrastrong coupling limit
recently derived. Considering the traditional spin-boson model as
illustration, we compare our results with predictions obtained through a
refined singular coupling master equation established in a very recent work,
and find excellent agreement. The comparison of the general expressions is
also briefly commented.

The dynamics of quantum systems strongly coupled to thermal baths have
recently received a lot of interest fuelled by hopes of understanding or even
discovering new phenomena in quantum transport Chin_2012 ; Ribeiro_2015 ;
Strasberg_2018 ; Correa_2019 ; Moreira_2020 ; Harush_2020 ; Dwiputra_2021 ;
Anto_2021 , quantum thermodynamics Gelbwaser_2015 ; Strasberg_2016 ;
Newman_2017 ; Perarnau_2018 ; Wertinek_2018 ; Newman_2020 ; Wiedmann_2020 ,
quantum sensing Correa_2017 ; Mehboudi_2019 ; Mejia_2021 , as well as
understanding better the underlying physics of some essential biological
functions Kolli_2012 ; Lambert_2013 ; Scholes_2017 ; Lambert_2020 . In most of
these applications, to know the steady state of the system strongly coupled to
the bath is often essential and sufficient. Such steady states greatly depart
from usual equilibrium steady states Smith_2014 ; Purkayastha_2020 ;
Cresser_2021 ; CLL_2021 . While the properties of the most general open
quantum evolutions have been known for long time thanks to the seminal paper
Sudarshan_1961 , their precise time evolutions and steady states are still a
current challenge of the theory of open quantum system FrancescoBook .

To obtain some information about strongly coupled steady states, several
strategies have been developed, including embedding techniques like reaction
coordinate Garg_1985 ; Smith_2014 ; Smith_2016 ; Strasberg_2016  and pseudo-
mode Garraway_1997 ; Pleasance_2017 ; Teretenkov_2019 ; Pleasance_2020 , or
numerical techniques (Hierarchical Equation of Motion) Tanimura_2019 ;
Lambert_2019 ; Lambert_2020 . One alternative strategy focusing directly on
the steady state without going through the description of the whole dynamics
is to consider the global steady state of the system and bath, namely the
global thermal state at the bath temperature Bach_2000 ; Frohlich_2004 ;
Merkli_2007 ; Mori_2008 ; Konenberg_2016 ; Merkli_2020  (see more details in
the following). The steady state of the system is then given by tracing out
the bath, which is sometimes referred to as the mean force Gibbs state
Cresser_2021 ; Trushechkin_2021 . Such partial trace is usually very
challenging, but can be done at least approximately assuming for instance a
weak coupling Mori_2008 ; Subasi_2012 ; Cresser_2021 ; Purkayastha_2020 . One
other interesting regime, and potentially containing more novelty, is the
ultrastrong coupling regime, when the strength of the coupling is larger than
the system’s energy scale. However, only few papers considered such
situations. In Cresser_2021 , Cresser and Anders provide the explicit
expression of the mean force Gibbs state in the limit of infinite coupling. In
particular, they show that it is of the same form as the steady state
Orman_2019 ; Goyal_2020  obtained using arguments from eisenselection
Zurek_2003 , although there are also some slight differences. Despite being an
interesting result, it would be welcome to have also information on how the
transition from the weak coupling limit to the ultrastrong coupling limit
happens, as well as on the steady state in intermediate regimes which are
experimentally more accessible.

In this perspective, introducing a technique inspired from the displaced
oscillator picture Irish_2005 ; Li_2021  for diagonalization of the quantum
Rabi model, we recover the infinite coupling limit of Cresser_2021  and go
beyond by providing first order corrections. In the regime of high bath
temperature or low bath frequency, we derive a very simple approximate
expression. We compare our results with the very recent derivation of
Trushechkin Trushechkin_2021 , in which he derives a master equation in the
so-called strong-decoherence regime, a refined version of the traditional
singular-coupling regime Gorini_1978 ; FrancescoBook , and find very good
agreement for the spin-boson model. Additionally, we show that the provided
corrective terms allow to sketch the transition from the ultrastrong coupling
limit to the weak coupling limit.

For arbitrary system, we briefly show in section IV of the Supplementary
Material that our general expression shares some strong similarities with the
one in Trushechkin_2021 , although there are also some unexplained
differences. .

##  I Mean force Gibbs state

We consider the following total Hamiltonian

${\cal H}_{SB}=H_{S}+H_{B}+\lambda H_{I},$ (1)

where $H_{B}:=\sum_{k}\omega_{k}a_{k}^{{\dagger}}a_{k}$ is the bath
Hamiltonian composed of the bosonic creation and annihilation operators
$a_{k}^{{\dagger}}$ and $a_{k}$, $\lambda$ is a dimensionless parameter
representing the strength of the system-bath coupling, $H_{I}:=AB$ is the
system-bath coupling term, and $B:=\sum_{k}g_{k}(a_{k}+a_{k}^{{\dagger}})$ is
the usual coupling bath operator. As in the quantum Brownian motion
FrancescoBook , there is an extra term $\lambda^{2}QA^{2}$ appearing during
the derivation of the solution which corresponds to the renormalization of the
system’s energies due to the interaction with the bath. This is also taken
into account explicitly in Cresser_2021 . The quantity $Q$, the “re-
organization energy” May_2011 ; Cresser_2021 , is defined as
$Q:=\int_{0}^{\infty}d\omega J(\omega)/\omega$, where $J(\omega)$ is the bath
spectral density, $J(\omega)=\sum_{k}g_{k}^{2}\delta(\omega-\omega_{k})$. Such
extra term is usually added initially, “by hand”, so that it cancels out later
on during the calculation when the renormalization takes place. Following this
procedure, we add the extra term $\lambda^{2}QA^{2}$ so that the total
Hamiltonian is now

$H_{SB}:=H_{S}+\lambda^{2}QA^{2}+\lambda AB+H_{B}.$ (2)

Note the difference of notation between the “natural” (1) and renormalized (2)
Hamiltonians. As a matter of completeness, we also mention in section III of
the Supplementary Material the derivation starting from the “natural”
Hamiltonian (1), which leads to renormalized final energies (or pseudo-
energies). This shows in particular that both derivations are equivalent, as
expected, but that the renormalization has to take place at some stage, either
initially or finally.

As mentioned in the introduction, the global steady state of the system and
bath is the global thermal state (assming $[H_{S},A]\neq 0$) Bach_2000 ;
Frohlich_2004 ; Merkli_2007 ; Mori_2008 ; Konenberg_2016 ; Merkli_2020

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

where $Z_{SB}:={\rm Tr}_{SB}\left(e^{-\beta H_{SB}}\right)$ is the partition
function. Then, the reduced steady state of $S$ is obtained by

$\rho_{S}^{\rm ss}:=Z_{SB}^{-1}{\rm Tr}_{B}[e^{-\beta H_{SB}}],$ (4)

which is sometimes called the mean force Gibbs state Cresser_2021 ;
Trushechkin_2021 . We are going to approximate this state in the regime of
ultrastrong coupling, when $\lambda$ becomes larger than the energy scale of
$H_{S}$.

##  II General derivation

We first split $H_{SB}$ in two, $H_{S}$ on one hand, and
$\lambda^{2}QA^{2}+\lambda AB+H_{B}$ on the other hand, and take it out of the
exponential using the usual identities Feynman_1951 . It leads to

$\displaystyle e^{-\beta H_{SB}}=e^{-\beta(H_{B}+\lambda
AB+\lambda^{2}QA^{2})}e^{-{\cal T}\int_{0}^{\beta}du\tilde{H}_{S}(u)},$ (5)

with ${\cal T}$ representing the “$\beta$-ordering operator”, acting on
inverse temperatures in the same way as the usual time ordering operator acts
on exponential time integrals, and

$\tilde{H}_{S}(u):=e^{u(H_{B}+\lambda
AB+\lambda^{2}QA^{2})}H_{S}e^{-u(H_{B}+\lambda AB+\lambda^{2}QA^{2})}.$ (6)

Now, by noticing that $H_{B}+\lambda
AB+\lambda^{2}QA^{2}=\sum_{k}\omega_{k}{\cal D}_{k}b_{k}^{{\dagger}}b_{k}{\cal
D}_{k}^{{\dagger}}$, where ${\cal D}_{k}:=e^{-\frac{\lambda
g_{k}}{\omega_{k}}(b_{k}^{{\dagger}}-b_{k})A}$ is a displacement operator by
an “amount” $\lambda g_{k}A/\omega_{k}$ acting on the mode $k$, we can re-
write $H_{B}+\lambda AB+\lambda^{2}QA^{2}$ as a “mixture of displaced baths”.
To see that, we denote by $a_{l}$ and $|a_{l}\rangle$ the eigenvalues and
corresponding eigenvectors of the observable $A$, and by ${\cal
D}_{k,l}:=e^{-\frac{\lambda g_{k}}{\omega_{k}}(b_{k}^{{\dagger}}-b_{k})a_{l}}$
the displacement operator by the quantity $\lambda g_{k}a_{l}/\omega_{k}$
acting only on the mode $k$. Then, we have

$\displaystyle H_{B}+\lambda AB+\lambda^{2}QA^{2}$ $\displaystyle=$
$\displaystyle\sum_{k}\omega_{k}{\cal D}_{k}b_{k}^{{\dagger}}b_{k}{\cal
D}_{k}^{{\dagger}}$ (7) $\displaystyle=$
$\displaystyle\sum_{k}\omega_{k}\sum_{l}|a_{l}\rangle\langle a_{l}|{\cal
D}_{k}b_{k}^{{\dagger}}b_{k}{\cal D}_{k}^{{\dagger}}$ $\displaystyle=$
$\displaystyle\sum_{k}\omega_{k}\sum_{l}|a_{l}\rangle\langle a_{l}|{\cal
D}_{k,l}b_{k}^{{\dagger}}b_{k}{\cal D}_{k,l}^{{\dagger}}$ $\displaystyle=$
$\displaystyle\sum_{l}|a_{l}\rangle\langle a_{l}|H_{B,l},$

where we defined the “displaced bath” $H_{B,l}:=\sum_{k}\omega_{k}{\cal
D}_{k,l}b_{k}^{{\dagger}}b_{k}{\cal D}_{k,l}^{{\dagger}}=H_{B}+\lambda
a_{l}B+\lambda^{2}a_{l}^{2}Q$, reminiscent of the displaced oscillator picture
Irish_2005 ; Li_2021 . From there we obtain,

$\displaystyle e^{u(H_{B}+\lambda AB+\lambda^{2}QA^{2})}$ $\displaystyle=$
$\displaystyle e^{u\sum_{l}|a_{l}\rangle\langle a_{l}|H_{B,l}}$ (8)
$\displaystyle=$ $\displaystyle\sum_{l}|a_{l}\rangle\langle
a_{l}|e^{uH_{B,l}},$

so that

$\displaystyle\tilde{H}_{S}(u)$ $\displaystyle=$
$\displaystyle\sum_{l,l^{\prime}}|a_{l}\rangle\langle
a_{l}|e^{uH_{B,l}}H_{S}|a_{l^{\prime}}\rangle\langle
a_{l^{\prime}}|e^{-uH_{B,l^{\prime}}}$ (9) $\displaystyle=$ $\displaystyle
H_{S}^{\rm pop}+\tilde{H}_{S}^{\rm coh}(u),$

where $H_{S}^{\rm pop}:=\sum_{l}h_{l}|a_{l}\rangle\langle a_{l}|$ and
$\tilde{H}_{S}^{\rm coh}(u):=\sum_{l\neq
l^{\prime}}h_{l,l^{\prime}}|a_{l}\rangle\langle
a_{l^{\prime}}|e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}$, defining $h_{l}:=\langle
a_{l}|H_{S}|a_{l}\rangle$, and $h_{l,l^{\prime}}:=\langle
a_{l}|H_{S}|a_{l^{\prime}}\rangle$.

Taking out $H_{S}^{\rm pop}$ from the exponential $e^{-{\cal
T}\int_{0}^{\beta}du\tilde{H}_{S}(u)}$, we obtain,

$\displaystyle e^{-{\cal T}\int_{0}^{\beta}du\tilde{H}_{S}(u)}$
$\displaystyle=$ $\displaystyle e^{-\beta H_{S}^{\rm pop}}e^{-{\cal
T}\int_{0}^{\beta}du\accentset{\approx}{H}_{S}^{\rm coh}(u)},$ (10)

with

$\displaystyle\accentset{\approx}{H}_{S}^{\rm coh}(u)$ $\displaystyle:=$
$\displaystyle e^{uH_{S}^{\rm pop}}\tilde{H}_{S}^{\rm coh}(u)e^{-uH_{S}^{\rm
pop}}$ (11) $\displaystyle=$ $\displaystyle\sum_{l\neq
l^{\prime}}h_{l,l^{\prime}}e^{u\omega_{l,l^{\prime}}}|a_{l}\rangle\langle
a_{l^{\prime}}|e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}},$

defining $\omega_{l,l^{\prime}}:=h_{l}-h_{l^{\prime}}$. Combining the
identities (5), (8) and (10), we arrive at

$\displaystyle\rho_{SB}^{\rm th}$ $\displaystyle=$ $\displaystyle
Z_{SB}^{-1}\sum_{l}|a_{l}\rangle\langle a_{l}|e^{-\beta H_{B,l}}e^{-\beta
H_{S}^{\rm pop}}e^{-{\cal T}\int_{0}^{\beta}du\accentset{\approx}{H}_{S}^{\rm
coh}(u)}$ $\displaystyle=$ $\displaystyle Z_{SB}^{-1}\sum_{l}e^{-\beta
h_{l}}|a_{l}\rangle\langle a_{l}|e^{-\beta H_{B,l}}e^{-{\cal
T}\int_{0}^{\beta}du\accentset{\approx}{H}_{S}^{\rm coh}(u)}.$

This leads to the following expression for the reduced steady state of $S$,

$\displaystyle\rho_{S}^{\rm ss}$ $\displaystyle=$ $\displaystyle
Z_{SB}^{-1}\sum_{l}e^{-\beta h_{l}}|a_{l}\rangle\langle a_{l}|{\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}e^{-{\cal
T}\int_{0}^{\beta}du\accentset{\approx}{H}_{S}^{\rm coh}(u)}\right].$

The toughest part is, as expected,

${\rm Tr}_{B}\left[e^{-\beta H_{B,l}}e^{-{\cal
T}\int_{0}^{\beta}du\accentset{\approx}{H}_{S}^{\rm
coh}(u)}\right]=\sum_{n=0}^{\infty}(-1)^{n}\int_{0}^{\beta}du_{1}\int_{0}^{u_{1}}du_{2}...\int_{0}^{u_{n-1}}du_{n}{\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}\accentset{\approx}{H}_{S}^{\rm
coh}(u_{1})\accentset{\approx}{H}_{S}^{\rm
coh}(u_{2})...\accentset{\approx}{H}_{S}^{\rm coh}(u_{n})\right].$ (14)

So far, we made no approximation, making (II) an exact expression. Since this
problem is not exactly solvable, we have to make some approximations in order
to reach an explicit form. Before that, let us see how one can recover the
inifinite coupling limit.

###  II.1 Recovering the inifinite coupling limit

From expression (II), it might not be obvious how one recovers the inifinite
coupling limit. It actually comes from the terms
$\accentset{\approx}{H}_{S}^{\rm coh}(u)$ which contain overlaps between
displaced baths, $e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}$. The strength of the
coupling between $S$ and each mode of the bath is represented by $g_{k}$.
Thus, when the strength of the coupling goes to infinity, the displaced baths
$H_{B,l}$ and $H_{B,l^{\prime}}$ tend to be displaced infinitely far apart
from each other (for$l\neq l^{\prime}$). Consequently, the overlap and the
expectation value ${\rm Tr}_{B}[e^{-\beta
H_{B,l}}e^{uH_{B,l^{\prime}}}e^{-uH_{B,l^{\prime\prime}}}]$ tends to zero for
increasing coupling strength.

Extending this reasoning to higher order terms, we see that all terms ${\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}\accentset{\approx}{H}_{S}^{\rm
coh}(u_{1})\accentset{\approx}{H}_{S}^{\rm
coh}(u_{2})...\accentset{\approx}{H}_{S}^{\rm coh}(u_{n})\right]$ contain
multiple overlaps of different displaced baths, and therefore tends to zero as
the coupling strength increases. Thus, in the infinite coupling limit, only
the first term of the sum in (II) is different from zero. It leads to

$\displaystyle\rho_{S}^{\rm ss,\infty}$ $\displaystyle=$ $\displaystyle
Z_{SB}^{-1}\sum_{l}e^{-\beta h_{l}}|a_{l}\rangle\langle a_{l}|{\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}\right]$ (15) $\displaystyle=$
$\displaystyle(Z_{S}^{\rm ss})^{-1}\sum_{l}e^{-\beta
h_{l}}|a_{l}\rangle\langle a_{l}|,$

where $Z_{S}^{\rm ss}:=\frac{Z_{SB}}{Z_{B}}$ and ${\rm Tr}_{B}\left[e^{-\beta
H_{B,l}}\right]={\rm Tr}_{B}\left[e^{-\beta H_{B}}\right]:=Z_{B}$ is the
partition function of the uncoupled bath. Note that since $\rho_{S}^{\rm
ss,\infty}$ is a normalized state, we also have the identity $Z_{S}^{\rm
ss}=\sum_{l}e^{-\beta h_{l}}$. The above expression (15) is exactly equal to
the one derived in Cresser_2021 , and recently re-obtained in Trushechkin_2021
.

###  II.2 First approximation

As seen in the previous section, the term $\accentset{\approx}{H}_{S}^{\rm
coh}(u)$ contains overlaps of displaced baths, so that terms ${\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}\accentset{\approx}{H}_{S}^{\rm
coh}(u_{1})\accentset{\approx}{H}_{S}^{\rm
coh}(u_{2})...\accentset{\approx}{H}_{S}^{\rm coh}(u_{n})\right]$ of
increasing order contain overlaps of increasing order and are therefore
significantly smaller than terms of lower orders (for large coupling
strength). With this assumption, we are going to retain only the first and
second term,

$\displaystyle\rho_{S}^{\rm ss}$ $\displaystyle\simeq$ $\displaystyle
Z_{SB}^{-1}\sum_{l}e^{-\beta h_{l}}|a_{l}\rangle\langle a_{l}|\Bigg{\\{}{\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}\right]$ (16) $\displaystyle-{\rm
Tr}_{B}\left[e^{-\beta
H_{B,l}}\int_{0}^{\beta}du\accentset{\approx}{H}_{S}^{\rm
coh}(u)\right]\Bigg{\\}}.$

The first term gives ${\rm Tr}_{B}\left[e^{-\beta H_{B,l}}\right]=Z_{B}$, as
already seen in the previous section. The second term gives,

$\displaystyle{\rm Tr}_{B}\left[e^{-\beta
H_{B,l}}\int_{0}^{\beta}du\accentset{\approx}{H}_{S}^{\rm
coh}(u)\right]=\sum_{l^{\prime}\neq
l^{\prime\prime}}h_{l^{\prime},l^{\prime\prime}}|a_{l^{\prime}}\rangle\langle
a_{l^{\prime\prime}}|$
$\displaystyle\times\int_{0}^{\beta}due^{u\omega_{l^{\prime},l^{\prime\prime}}}{\rm
Tr}_{B}\left[e^{-\beta
H_{B,l}}e^{uH_{B,l^{\prime}}}e^{-uH_{B,l^{\prime\prime}}}\right].$ (17)

Since, when injected in (16), this expression will be multiplied by
$|a_{l}\rangle\langle a_{l}|$ on the left-hand side, it is enough to compute

$\displaystyle\sum_{l^{\prime}\neq l}h_{l,l^{\prime}}|a_{l}\rangle\langle
a_{l^{\prime}}|\int_{0}^{\beta}due^{u\omega_{l,l^{\prime}}}{\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}\right],$

where the index $l$ and $l^{\prime}$ were made equal. After some
manipulations, we can obtain the following expression for ${\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}\right]$
(details provided in section I of the Supplementary Material),

$\displaystyle{\rm Tr}_{B}\left[e^{-\beta
H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}\right]$
$\displaystyle=Z_{B}e^{-\lambda^{2}a_{l^{\prime},l}^{2}\int_{0}^{\infty}\frac{J(\omega)}{\omega^{2}}(e^{u\omega}-1)\left(1-\frac{1-e^{-u\omega}}{1-e^{-\omega\beta}}\right)},$
(19)

where $a_{l^{\prime},l}:=a_{l^{\prime}}-a_{l}$. We finally obtain, up to
second order,

$\displaystyle\rho_{S}^{\rm ss}$ $\displaystyle=$
$\displaystyle\sum_{l}p_{l}^{\rm ss}|a_{l}\rangle\langle
a_{l}|-\sum_{l,l^{\prime};l\neq l^{\prime}}p_{l}^{\rm
ss}h_{l,l^{\prime}}f_{l,l^{\prime}}(\beta)|a_{l}\rangle\langle
a_{l^{\prime}}|,$

with $p_{l}^{\rm ss}:=e^{-\beta h_{l}}/Z_{S}^{\rm ss}$ and

$\displaystyle f_{l,l^{\prime}}(\beta)$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}due^{u\omega_{l,l^{\prime}}}e^{-\lambda^{2}a_{l^{\prime},l}^{2}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}(e^{u\omega}-1)\left(1-\frac{1-e^{-\omega
u}}{1-e^{-\omega\beta}}\right)}.$

Note that we have the following identity (see section II of the Supplementary
Material) $p_{l}^{\rm ss}f_{l,l^{\prime}}(\beta)=p_{l^{\prime}}^{\rm
ss}f_{l^{\prime},l}(\beta)$, implying $\langle a_{l}|\rho_{S}^{\rm
ss}|a_{l^{\prime}}\rangle^{*}=\langle a_{l^{\prime}}|\rho_{S}^{\rm
ss}|a_{l}\rangle$, as it should be.

Additionally, $f_{+,-}(\beta)$ tends to zero when $\lambda$ goes to infinity,
so that we recover the infinite coupling limit (15). One can also see that the
first order corrections affect only the coherences (in the eigenbasis of $A$).
These observations coincide with the ones in Trushechkin_2021 . Finally, the
expression obtained starting from the natural Hamiltonian (1) instead of the
renormalized one (2) are the same as (II.2) and (II.2) but substituting the
pseudo-energies $h_{l}$ by the renormalized ones $h_{l}-\lambda^{2}a_{l}^{2}Q$
(see section III of the Supplementary Material).

###  II.3 Approximate expression of $f_{l,l^{\prime}}(\beta)$

Depending on the bath spectral density, it might not be possible to obtain an
exact analytical expression of $f_{l,l^{\prime}}(\beta)$, so that some
approximations would have to be made. However, assuming the bath temperature
is high, or equivalently that the dominant frequencies in the bath are low, we
can obtain an approximate expression of $f_{l,l^{\prime}}(\beta)$ witout even
specifying the form of the spectral density. More precisely, we assume that
$J(\omega)$ vanishes for $\omega\geq\omega_{c}$, where
$\omega_{c}\leq\beta^{-1}$. Then, for $\omega\in[0;\beta^{-1}]$, the factor
$(e^{u\omega}-1)\left(1-\frac{1-e^{-u\omega}}{1-e^{-\omega\beta}}\right)$ can
be approximated by

$\frac{e^{u\omega}-1}{\omega^{2}}\left(1-\frac{1-e^{-u\omega}}{1-e^{-\omega\beta}}\right)=\frac{u}{\omega}\left(1-\frac{u}{\beta}\right)+{\cal
O}(\omega\beta),$ (22)

which is actually a very good approximation as soon as $\omega\beta\leq 1$.
With that we obtain

$\displaystyle
e^{-\lambda^{2}a_{l^{\prime},l}^{2}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}(e^{u\omega}-1)\left(1-\frac{1-e^{-\omega
u}}{1-e^{-\omega\beta}}\right)}\simeq
e^{-\lambda^{2}a_{l^{\prime},l}^{2}u\left(1-\frac{u}{\beta}\right)Q},$

and

$\displaystyle f_{l,l^{\prime}}(\beta)$ $\displaystyle\simeq$
$\displaystyle\int_{0}^{\beta}due^{u\omega_{l,l^{\prime}}}e^{-\lambda^{2}a_{l^{\prime},l}^{2}u\left(1-\frac{u}{\beta}\right)Q}$
$\displaystyle=\frac{1}{\lambda|a_{l^{\prime},l}|}\sqrt{\frac{\beta}{Q}}\Bigg{\\{}{\rm
DF}\left[\frac{1}{2\lambda|a_{l^{\prime},l}|}\sqrt{\frac{\beta}{Q}}(\lambda^{2}a_{l^{\prime},l}^{2}Q-\omega_{l,l^{\prime}})\right]$
$\displaystyle+e^{\beta\omega_{l,l^{\prime}}}{\rm
DF}\left[\frac{1}{2\lambda|a_{l^{\prime},l}|}\sqrt{\frac{\beta}{Q}}(\lambda^{2}a_{l^{\prime},l}^{2}Q+\omega_{l,l^{\prime}})\right]\Bigg{\\}},$

where ${\rm DF}(x):=e^{-x^{2}}\int_{0}^{x}{\rm d}u~{}e^{u^{2}}$ is sometimes
referred to as the Dawson function.

The above expression can be further simplified as follows. The strong coupling
regime can be characterized by a re-organization energy $\lambda^{2}Q$
comparable to, or larger than, the energy scale of the system Cresser_2021 ;
CLL_2021 . Thus, in the strong coupling regime one can expect to have
$\lambda^{2}Q\gg{\rm max}_{l}~{}|h_{l}|$, implying

$\displaystyle\frac{1}{2\lambda|a_{l^{\prime},l}|}\sqrt{\frac{\beta}{Q}}(\lambda^{2}a_{l^{\prime},l}^{2}Q\pm\omega_{l,l^{\prime}})$
$\displaystyle\sim$
$\displaystyle\frac{1}{2\lambda|a_{l^{\prime},l}|}\sqrt{\frac{\beta}{Q}}\lambda^{2}a_{l^{\prime},l}^{2}Q$
(25) $\displaystyle\sim$ $\displaystyle\sqrt{\beta\lambda^{2}Q}.$

Finally, since $DF(x)=\frac{1}{2x}+{\cal O}(x^{-3})$ for $x\gg 1$ (the
approximation is actaully very good for $x\geq 3$), we can find the following
simple approximate expression for $f_{l,l^{\prime}}(\beta)$ asuming
$\lambda^{2}Q\gg\beta^{-1}$,

$f_{l,l^{\prime}}(\beta)=\frac{1+e^{\omega_{l,l^{\prime}}\beta}}{\lambda^{2}a_{l^{\prime},l}^{2}Q}+\frac{\omega_{l,l^{\prime}}(1-e^{\omega_{l,l^{\prime}}\beta})}{\lambda^{4}a_{l^{\prime},l}^{4}Q^{2}}+{\cal
O}[(\lambda^{2}Q\beta)^{-3}].$ (26)

##  III Example: spin-boson model

As illustration of our results, we consider the versatile and famous spin-
boson model Leggett_1987  characterised by the following total Hamiltonian

${\cal
H}_{SB}=\frac{\epsilon}{2}\sigma_{z}+\frac{\Delta}{2}\sigma_{x}+\lambda\sigma_{z}B+H_{B},$
(27)

where $\sigma_{x}$ and $\sigma_{z}$ are the Pauli matrices. The spin
Hamiltonian $H_{S}$ can be re-written as

$H_{S}=\frac{\epsilon}{2}\sigma_{z}+\frac{\Delta}{2}\sigma_{x}=\frac{\omega_{S}}{2}\big{(}|e\rangle\langle
e|-|g\rangle\langle g|\big{)},$ (28)

with $\omega_{S}:=\sqrt{\epsilon^{2}+\Delta^{2}}$,

$\displaystyle|e\rangle:=\frac{(\omega_{s}+\epsilon)|+\rangle+\Delta|-\rangle}{\sqrt{2\omega_{S}(\omega_{S}+\epsilon)}},$
$\displaystyle|g\rangle:=\frac{-\Delta|+\rangle+(\omega_{S}+\epsilon)|-\rangle}{\sqrt{2\omega_{S}(\omega_{S}+\epsilon)}},$
(29)

and $|\pm\rangle$ denotes the eigenstates of $\sigma_{z}$. Since coupling
observable $A$ is equal to $\sigma_{z}$, we have $a_{l=\pm}=\pm 1$,
$a_{l,l^{\prime}}=a_{+,-}=-a_{-,+}=2$, $h_{l=\pm}=\pm\epsilon/2$,
$h_{l,l^{\prime}}=h_{+,-}=h_{-,+}=\Delta/2$, and
$\omega_{l,l^{\prime}}=h_{l}-h_{l^{\prime}}=\omega_{+-}=-\omega_{-+}=\epsilon$.
Additionally, since $\sigma_{z}^{2}={\mathbb{I}}$, the renormalized
Hamiltonian is equal to

$H_{SB}=\frac{\epsilon}{2}\sigma_{z}+\frac{\Delta}{2}\sigma_{x}+\lambda\sigma_{z}B+H_{B}+\lambda^{2}Q,$
(30)

which simply corresponds to redefining the origin of the spin energy. This can
also be verified from the renormalized pseudo-energies defined in section III
of the Supplementary Material, which are given by ${\tt
h}_{l=\pm}:=h_{l=\pm}-a_{l}^{2}\lambda^{2}Q=\pm\epsilon/2-\lambda^{2}Q$. Thus,
for the spin-boson model, the energy renormalization induced by the
interaction with the bath has no impact on the reduced steady state.

Applying expression (II.2), we obtain

$\displaystyle\rho_{S}^{\rm ss}$ $\displaystyle=$
$\displaystyle\frac{1}{e^{-\beta\epsilon/2}+e^{\beta\epsilon/2}}\Big{[}e^{-\beta\epsilon/2}|+\rangle\langle+|+e^{\beta\epsilon/2}|-\rangle\langle-|$
$\displaystyle-\frac{\Delta}{2}e^{-\beta\epsilon/2}f_{+,-}(\beta)|+\rangle\langle-|-\frac{\Delta}{2}e^{\beta\epsilon/2}f_{-,+}(\beta)|-\rangle\langle+|\Big{]},$

with

$\displaystyle f_{+,-}(\beta)$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}due^{\epsilon
u}e^{-4\lambda^{2}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}(e^{u\omega}-1)\left(1-\frac{1-e^{-\omega
u}}{1-e^{\omega\beta}}\right)},$ $\displaystyle f_{-,+}(\beta)$
$\displaystyle=$ $\displaystyle\int_{0}^{\beta}due^{-\epsilon
u}e^{-4\lambda^{2}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}(e^{u\omega}-1)\left(1-\frac{1-e^{-\omega
u}}{1-e^{\omega\beta}}\right)}.$

###  III.1 High temperature approximation

Without specifying explicitly the bath spectral density, if we consider that
$J(\omega)$ vanishes for frequencies smaller than $\beta^{-1}$, then the
approximation (II.3) of section II.3 applies, and we have,

$\displaystyle f_{+,-}(\beta)$ $\displaystyle=$
$\displaystyle\frac{1}{2\lambda}\sqrt{\frac{\beta}{Q}}\Bigg{\\{}{\rm
DF}\left[\frac{1}{4\lambda}\sqrt{\frac{\beta}{Q}}(4\lambda^{2}Q-\epsilon)\right]$
(33) $\displaystyle+e^{\epsilon\beta}{\rm
DF}\left[\frac{1}{4\lambda}\sqrt{\frac{\beta}{Q}}(4\lambda^{2}Q+\epsilon)\right]\Bigg{\\}}.$

Assuming furthermore $\lambda^{2}Q\gg\epsilon$ and $\lambda^{2}Q\beta\gg 1$,
characteristics of the ultrastrong coupling regime, we can use expression
(26), leading to

$f_{+,-}(\beta)=\frac{1+e^{\epsilon\beta}}{4\lambda^{2}Q}+\frac{\epsilon(1-e^{\epsilon\beta})}{16\lambda^{4}Q^{2}}+{\cal
O}[(\lambda^{2}Q\beta)^{-3}].$ (34)

(a)
(b)

Figure 1: Plots of the steady state coherences (a) $c^{\rm
ss}:=\langle+|\rho_{S}^{\rm ss}|-\rangle$ in the eigenbasis of $A$, and (b)
$c^{\rm ss}_{e,g}:=\langle e|\rho_{S}^{\rm ss}|g\rangle$ in the eigenbasis of
$H_{S}$, both in function of $\lambda^{2}Q$, which represents the coupling
strength, in unit of $\epsilon$, for $\Delta=0.7\epsilon$,
$\omega_{c}=0.25\epsilon$ and $\beta=1/\epsilon$. The orange solid curve
corresponds to the coherence of the steady state (III) in the basis
$\\{|+\rangle,|-\rangle\\}$ and $\\{|e\rangle,|g\rangle\\}$ for panel (a) and
(b), respectively, and using the high temperature approximation (33) of
$f_{+,-}(\beta)$. The blue dashed line corresponds to the same coherences
using the further simplified expression (34) of $f_{+,-}(\beta)$. The purple
dotted line is the same coherences of the steady state given by the high
temperature expression derived in Trushechkin_2021 . The vertical red dashed
line denotes the approximate limit of strong coupling, when
$\lambda^{2}Q\sim\epsilon$, and so the approximate limit of validity of the
curves. Additionally, although the simplified expression (34) of
$f_{+,-}(\beta)$ (blue dashed curve) is supposed to be valid only for
$\lambda^{2}Q\gg\epsilon$ and $\lambda^{2}Q\gg\beta^{-1}$, delimited by the
grey vertical dotted line, one can see that it approximates well (33) (orange
solid line) even for $\lambda^{2}Q\sim\epsilon$ and
$\lambda^{2}Q\sim\beta^{-1}$. One can see in (a) that for increasing coupling
strength, the steady state tends to a diagonal state in the eigenbasis of $A$,
as expected from the ultrastrong coupling limit (15) and Orman_2019 ;
Goyal_2020 ; Cresser_2021 . Conversely, for decreasing coupling strength, one
can see in (b) the beginning of a progressive transition to a diagonal state
in the eigenbasis of $H_{S}$, as expected from the weak coupling limit.

###  III.2 Comparison with previous result

In a recent paper Trushechkin_2021 , Trushechkin obtains corrections to the
ultrastrong coupling limit by deriving an solving a refined singular coupling
master equation. In this section, we compare our result with Trushechkin’s
one, who also considered the spin-boson model and the high temperature regime
as illustrative example. Although our approximate expression (33) of
$f_{+,-}(\beta)$ is valid for any form of bath spectral density, Trushechkin
considers an over-damped spectral density, also sometimes called Lorentz-Drude
spectral density,

$J(\omega)=\frac{2Q}{\pi}\frac{\omega_{c}\omega}{\omega_{c}^{2}+\omega^{2}},$
(35)

where $Q$ is precisely the re-organization energy associated with the over-
damped spectral density, and $\omega_{c}$ represents the cutoff frequency.

Note that the dimensionless factor $\lambda$ is not explicitly present in
Trushechkin_2021 , but one can make it appearing by multiplying the spectral
density by $\lambda^{2}$. Then, the reorganisation energy becomes
$\lambda^{2}Q$. Additionally, since the steady state populations do not
contain correction with respect to the infinite coupling limit (15), we only
need to compare the steady state coherence.

In Fig. 1 (a), we have plotted the steady state coherence in the eigenbasis of
$A$, $c^{\rm ss}:=\langle+|\rho_{S}^{\rm ss}|-\rangle$, in function of
$\lambda^{2}Q$, for $\delta=0.7\epsilon$, $\omega_{c}=0.25\epsilon$ and
$\beta=1/\epsilon$. The orange solid line is the plot using the high
temperature approximation of $f_{+,-}(\beta)$ (33), while the blue dashed line
is using the further simplified expression of $f_{+,-}(\beta)$ (34). The first
important observation is that the coherence indeed vanishes as the coupling
strength increases, but after a fast decay around $\lambda^{2}Q/\epsilon\sim
1$, a slow decay takes place. Secondly, one can see that (34) (blue dashed
line) is a good approximation of (33) (orange solid line) even for coupling
strength such that $\lambda^{2}Q\sim\epsilon$ and
$\lambda^{2}Q\sim\beta^{-1}$, while it is supposed to happen only for
$\lambda^{2}Q\gg\epsilon$ and $\lambda^{2}Q\gg\beta^{-1}$, denoted by the grey
vertical dotted line.

The purple dotted line represents the steady state coherence as given by the
high temperature expression of Trushechkin_2021 . One can see a perfect
agreement with our prediction, even beyond the strong coupling regime,
delimited by the red vertical dashed line. However, one should not conclude
from this that the expressions’ validities go beyond the strong coupling
regime.

In Fig.1 (b), the plots represent the steady state coherences $c^{\rm
ss}_{e,g}:=\langle e|\rho_{S}^{\rm ss}|g\rangle$ in the eigenbasis of $H_{S}$,
still in function of $\lambda^{2}Q$. The value of the other parameters and the
colour conventions are the same as in Fig. 1 (a). From the weak coupling
theory, we expect the steady state to be diagonal in the eigenbasis of $H_{S}$
for vanishing coupling strength. Fig. 1 (b) shows precisely this progressive
transition to a diagonal state in the eigenbasis of $H_{S}$ as the coupling
strength decreases. However, the regime of validity of the derived expression
is limited from below by the red vertical dashed line. As expected, below this
limit, the behaviour of the curves start deviating from the expected one,
namely convergence to zero coherence.

Figure 2: Plots of the steady state coherences $c^{\rm
ss}=\langle+|\rho_{S}^{\rm ss}|-\rangle$ in function of the inverse
temperature $\beta$, in unit of $\epsilon^{-1}$, for $\Delta=0.7\epsilon$,
$\omega_{c}=0.5\epsilon$ and $\lambda^{2}Q=5\epsilon$. The colour conventions
of the plots are the same as in Fig. 1. The vertical red dashed line denotes
approximately the upper bound validity of the high temperature approximation,
when $\omega_{c}\beta\sim 0.5$, and therefore the validity of the three
expressions. The vertical grey dotted line represents the approximate limit
for $\lambda^{2}Q\beta\gg 1$, and the validity of the simplified expression
(34) of $f_{+,-}(\beta)$.

In the next graph, Fig. 2, the steady state coherence $c^{\rm
ss}=\langle+|\rho_{S}^{\rm ss}|-\rangle$ is plotted in function of the inverse
temperature $\beta$ for $\delta=0.7\epsilon$, $\omega_{c}=0.1\epsilon$ and
$\lambda^{2}Q=5\epsilon$. The colour conventions are the same as in Fig. 1.
One can see again that both our high temperature expression and the one
provided in Trushechkin_2021  coincide exactly for the whole range of the
plot, even though their validity is only below the red vertical dashed line,
for which $\omega_{c}\beta\leq 0.5$. Note that Lorentz-Drude spectral
densities have a very long high frequency “tail”, implying that we cannot
rigorously consider that $J(\omega)$ vanishes for $\omega\geq 2\omega_{c}$ but
more for $\omega\geq 10\omega_{c}$. This would imply that the high temperature
approximation is actually valid only for $\omega_{c}\beta\leq 1/10$. However,
as shown in Fig. 3, the criterion $\omega_{c}\beta\leq 0.5$ seems acceptable.
Fig. 2 shows that at high temperature (low $\beta$), the coherence $c^{\rm
ss}=\langle+|\rho_{S}^{\rm ss}|-\rangle$ vanishes, reproducing a behaviour of
the coherence $c^{\rm ss}=\langle e|\rho_{S}^{\rm ss}|g\rangle$ in weak
coupling regime Cresser_2021 ; CLL_2021 . Finally, the simplified expression
(34) of $f_{+,-}(\beta)$ is a good approximation of the high temperature
expression (33) only in the expected validity range, beyond the grey vertical
dotted line, for which $\lambda^{2}Q\gg\beta^{-1}$.

Figure 3: Plots of the steady state coherences $c^{\rm
ss}=\langle+|\rho_{S}^{\rm ss}|-\rangle$ in function of the cutoff frequency
$\omega_{c}$, in unit of $\epsilon^{-1}$, for $\Delta=0.7\epsilon$,
$\beta=1/\epsilon$ and $\lambda^{2}Q=5\epsilon$. The colour conventions of the
plots are the same as in the previous figures. The vertical red dashed line
denotes approximately the validity of the high temperature approximation, when
$\omega_{c}\beta\sim 0.5$ (see further explanation in the text).

In the last graph, Fig. 3, the steady state coherence $c^{\rm
ss}=\langle+|\rho_{S}^{\rm ss}|-\rangle$ is plotted in function of
$\omega_{c}$, for $\Delta=0.7\epsilon$, $\beta=1/\epsilon$ and
$\lambda^{2}Q=5\epsilon$. For these parameters, the conditions for the
validity of the simplified expression (34) of $f_{+,-}(\beta)$ are not fully
satisfied, reason why the blue dotted curve is significantly distant from the
orange solid curve. The main point of this graph is the validity condition of
the high temperature (or low bath frequency) approximation. As commented
above, due to the long high frequency “tail” of Lorentz-Drude spectral
densities, one would expect the validity condition to be something like
$\omega_{c}\beta\leq 1/10$. However, we can see in Fig. 3 that both curves
become constant and coincide below the red vertical dashed line, which could
suggest a larger range of validity than expected.

##  IV Conclusion

Introducing an alternative technique inspired from the displaced oscillator
picture used in Irish_2005 ; Li_2021 , we obtained a perturbative expansion
(II) and (14) of the mean force Gibbs state, the reduced steady state of $S$.
In the ultrastrong coupling regime, we obtain analytical expressions for the
contributions of the first orders, (II.2) and (II.2). The zero order term
confirms the result of Cresser_2021 ; Trushechkin_2021 . As briefly commented
section IV of the Supplementary Material, we left for future research the full
resolution of the exact connections between (II.2) and the corresponding
expression in Trushechkin_2021 , and in particular the finite temperature
integral versus the infinite temporal integral.

In the high temperature regime, a simple and explicit general expression is
derived (II.3) for the first order corrections. This expression is benchmarked
for the spin-boson model against the results of Trushechkin_2021  obtained
through a strong-decoherence-limit master equation. The agreement is perfect
within the regime of validity. An even simpler expression (26) is derived for
very large coupling strength and also benchmarked. Thanks to these first order
corrections, we are able to establish the beginning of the transition from the
ultrastrong coupling regime to the weak coupling regime, Fig. 1. Still, higher
order corrections must be derived to be able to reconstruct the full
transition and provide a good approximation of the steady state in all
regimes.

This might be possible by directly considering the next term in the expansion
(14), or by combining the present technique with some approximate
diagonalization of the full Hamiltonian. Additionally, benchmarking with other
techniques like reaction coordinate of pseudo-modes would be interesting and
instructive.

## Acknowledgements

I would like to thank Anton Trushechkin for interesting ongoing discussions on
this problem. I acknowledge the support of the National Institute for
Theoretical Physics (NITheP) of the Republic of South Africa.

## References

  * (1) A. W. Chin, S. F. Huelga ans 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. doi:10.1098/rsta.2011.0224
  * (2) P. Ribeiro and V. R. Vieira, Non-Markovian effects in electronic and spin transport, Phys. Rev. B 92, 100302(R) (2015). doi:10.1103/PhysRevB.92.100302
  * (3) P. Strasberg, G. Schaller, T. L. Schmidt, and M. Esposito, Fermionic reaction coordinates and their application to an autonomous Maxwell demon in the strong-coupling regime, Phys. Rev. B 97, 205405 (2018). doi:10.1103/PhysRevB.97.205405
  * (4) L. A. Correa, B. Xu, B. Morris, and G. Adesso, Pushing the limits of the reaction-coordinate mapping, J. Chem. Phys. 151, 094107 (2019). doi:10.1063/1.5114690
  * (5) S. V. Moreira, B. Marques, R. R. Paiva, L. S. Cruz, D. O. Soares-Pinto, and F. L. Semião, Enhancing quantum transport efficiency by tuning non-Markovian dephasing, Phys. Rev. A 101, 012123 (2020). doi:10.1103/PhysRevA.101.012123
  * (6) E. Zerah-Harush and Y. Dubi, Effects of disorder and interactions in environment assisted quantum transport, Phys. Rev. Research 2, 023294 (2020). doi:10.1103/PhysRevResearch.2.023294
  * (7) D. Dwiputra and F. P. Zen, Environment-assisted quantum transport and mobility edges, Phys. Rev. A 104, 022205 (2021). doi:10.1103/PhysRevA.104.022205
  * (8) Nicholas Anto-Sztrikacs and Dvira Segal, Strong coupling effects in quantum thermal transport with the reaction coordinate method, arXiv:2103.05670.
  * (9) D. Gelbwaser-Klimovsky and A. Aspuru-Guzik, Strongly Coupled Quantum Heat Machines, J. Phys. Chem. Lett. 6, 3477-3482 (2015). doi:10.1021/acs.jpclett.5b01404
  * (10) 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). doi:10.1088/1367-2630/18/7/073007
  * (11) D. Newman, F. Mintert, and A. Nazir, Performance of a quantum heat engine at strong reservoir coupling, Phys. Rev. E 95, 032139 (2017). doi:10.1103/PhysRevE.95.032139
  * (12) M. Perarnau-Llobet, H. Wilming, A. Riera, R. Gallego, and J. Eisert, Strong Coupling Corrections in Quantum Thermodynamics, Phys. Rev. Lett. 120, 120602 (2018). doi:10.1103/PhysRevLett.120.120602
  * (13) 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). doi:10.1063/1.5040898
  * (14) D. Newman, F. Mintert, and A. Nazir, Quantum limit to nonequilibrium heat-engine performance imposed by strong system-reservoir coupling, Phys. Rev. E 101, 052129 (2020). doi:10.1103/PhysRevE.101.052129
  * (15) 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). doi:10.1088/1367-2630/ab725a
  * (16) L. A. Correa, M. Perarnau-Llobet, K. V. Hovhannisyan, S. Hernández-Santana, M. Mehboudi, and A. Sanpera, Enhancement of low-temperature thermometry by strong coupling, Phys. Rev. A 96, 062103 (2017). doi:10.1103/PhysRevA.96.062103
  * (17) M. Mehboudi, A. Lampo, C. Charalambous, L. A. Correa, M. A. García-March and M. Lewenstein, Using Polarons for sub-nK Quantum Nondemolition Thermometry in a Bose-Einstein Condensate, Phys. Rev. Lett. 122, 030403 (2019). doi:10.1103/PhysRevLett.122.030403
  * (18) 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). doi:10.1088/2058-9565/abdca5
  * (19) 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). doi:10.1063/1.4764100
  * (20) N. Lambert, Y.-N. Chen, Y.-C. Cheng, C.-M. Li, G.-Y. Chen and Franco Nori, Quantum biology, Nat. Physics 9, 10-18 (2013). doi:10.1038/nphys2474
  * (21) G. Scholes, G. Fleming, L. Chen, et al., Using coherence to enhance function in chemical and biophysical systems, Nature 543, 647–656 (2017). doi:10.1038/nature21425
  * (22) 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.
  * (23) J. Iles-Smith, N. Lambert, and A. Nazir, Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems, Phys. Rev. A 90, 032114 (2014). doi:10.1103/PhysRevA.90.032114
  * (24) A. Purkayastha, G. Guarnieri, M. T. Mitchison, R. Filip, and J. Goold, Tunable phonon-induced steady-state coherence in a double-quantum-dot charge qubit, npj Quant. Inf. 6, 27 (2020). doi:10.1038/s41534-020-0256-6
  * (25) J. Cresser and J. Anders, Weak and ultrastrong coupling limits of the quantum mean force Gibbs state, arXiv:2104.12606.
  * (26) C. L. Latune, Steady state in strong bath coupling: perturbative expansion versus reaction coordinate, arXiv:2110.
  * (27) E. C. G. Sudarshan, P. M. Mathews, and Jayaseetha Rau, Stochastic Dynamics of Quantum-Mechanical Systems, Phys. Rev. 121, 920 (1961). doi:10.1103/PhysRev.121.920
  * (28) H. Breuer and F. Petruccione, Theory of Open Quantum Systems, (Oxford, Oxford, 2002).
  * (29) A. Garg, J. N. Onuchic, and V. Ambegaokar, Effect of friction on electron transfer in biomolecules, J. Chem. Phys. 83, 4491 (1985). doi:10.1063/1.449017
  * (30) 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). doi:10.1063/1.4940218
  * (31) B. M. Garraway, Nonperturbative decay of an atomic system in a cavity, Phys. Rev. A. 55, 2290 (1997). doi:10.1103/PhysRevA.55.2290
  * (32) G. Pleasance and B. M. Garraway, Application of quantum Darwinism to a structured environment, Phys. Rev. A 96, 062105 (2017). doi:10.1103/PhysRevA.96.062105
  * (33) A.E. Teretenkov, Pseudomode Approach and Vibronic Non-Markovian Phenomena in Light-Harvesting Complexes, Proc. Steklov Inst. Math. 306, 242–256 (2019). doi:10.1134/S0081543819050201
  * (34) 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). doi:10.1103/PhysRevResearch.2.043058
  * (35) Yoshitaka Tanimura, Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM), J. Chem. Phys. 153, 020901 (2020). doi:10.1063/5.0011599
  * (36) N. Lambert, S. Ahmed, M. Cirio, F. Nori, Modelling the ultra-strongly coupled spin-boson model with unphysical modes, Nat Commun 10, 3721 (2019). doi:10.1038/s41467-019-11656-1
  * (37) V. Bach, J. Fröhlich, and I. M. Sigal, Return to Equilibrium, J. Math. Phys. 41, 3985 (2000). doi:10.1063/1.533334
  * (38) J. Frohlich, Marco Merkli, Another Return of “Return to Equilibrium”, Commun. Math. Phys. 251, 235–262 (2004). doi:10.1007/s00220-004-1176-6
  * (39) M. Merkli, I. M. Sigal, and G. P. Berman, Decoherence and Thermalization, Phys. Rev. Lett. 98, 130401 (2007). doi:10.1103/PhysRevLett.98.130401
  * (40) 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). doi:10.1143/JPSJ.77.124005
  * (41) M. Könenberg and M. Merkli, On the irreversible dynamics emerging from quantum resonances, J. Math. Phys. 57, 033302 (2016). doi:10.1063/1.4944614
  * (42) Marco Merkli, Quantum Markovian master equations: Resonance theory shows validity for all time scales, Ann. Phys. 412, 16799 (2020). doi:10.1016/j.aop.2019.167996
  * (43) A. Trushechkin, Quantum master equations and steady states for the ultrastrong-coupling limit and the strong-decoherence limit, arxiv:2109.01888
  * (44) Y. Subaşı, C. H. Fleming, J. M. Taylor, and B. L. Hu, Equilibrium states of open quantum systems in the strong coupling regime, Phys. Rev. E 86, 061132 (2012). doi:10.1103/PhysRevE.86.061132
  * (45) P. L. Orman and R. Kawai, A qubit strongly interacting with a bosonic environment: Geometry of thermal states, arXiv:2010.09201v1 (2020).
  * (46) K. Goyal and R. Kawai, Steady state thermodynamics of two qubits strongly coupled to bosonic environments, Phys. Rev. Research 1, 033018 (2019). doi:10.1103/PhysRevResearch.1.033018
  * (47) W. H. Zurek, Decoherence, einselection, and the quantum origins of the classical, Rev. Mod. Phys. 75, 715 (2003). doi:10.1103/RevModPhys.75.715
  * (48) 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). doi:10.1103/PhysRevB.72.195410
  * (49) Zi-Min Li and Murray T. Batchelor, Generalized adiabatic approximation to the quantum Rabi model, Phys. Rev. A 104, 033712 (2021). doi:10.1103/PhysRevA.104.033712
  * (50) V. Gorini, A. Frigerio, M.Verri, A. Kossakowski, and E. C. G. Sudarshan, Properties of quantum Markovian master equations, Rep. Math. Phys. 13, 149-173 (1978). doi:10.1016/0034-4877(78)90050-2
  * (51) V. May and O Kühn, Charge and Energy Transfer Dynamics in Molecular Systems (Wiley-VCH, Weinheim, 2011).
  * (52) Richard P. Feynman, An Operator Calculus Having Applications in Quantum Electrodynamics, Phys. Rev. 84, 108-128 (1951). doi:10.1103/PhysRev.84.108
  * (53) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987). doi:10.1103/RevModPhys.59.1
  * (54) W. H. Louisell, Radiation and Noise in Quantum Electronics (New York: McGraw-Hill), p. 116, (1964).

##  Appendix A Details of the computation of ${\rm Tr}_{B}\left[e^{-\beta
H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}\right]$

In order to compute ${\rm Tr}_{B}\left[e^{-\beta
H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}\right]$, one can first rearranged
it in the following way,

$\displaystyle{\rm Tr}_{B}\left[e^{-\beta
H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}\right]={\rm
Tr}_{B}\left[e^{-(\beta-u)H_{B,l}}e^{-uH_{B,l^{\prime}}}\right]$
$\displaystyle\hskip 14.22636pt={\rm
Tr}_{B}\left[\underset{k}{\otimes}\left({\cal
D}_{k,l}e^{-(\beta-u)b_{k}^{{\dagger}}b_{k}}{\cal
D}^{{\dagger}}_{k,l}\right)e^{-uH_{B,l^{\prime}}}\right]$ $\displaystyle\hskip
14.22636pt={\rm Tr}_{B}\left[e^{-(\beta-u)H_{B}}e^{-u\sum_{k}\omega_{k}{\cal
D}_{k,l}^{{\dagger}}{\cal D}_{k,l^{\prime}}b_{k}^{{\dagger}}b_{k}{\cal
D}_{k,l^{\prime}}^{{\dagger}}{\cal D}_{k,l}}\right]$ $\displaystyle\hskip
14.22636pt={\rm Tr}_{B}\left[e^{-(\beta-u)H_{B}}e^{-u(H_{B}+\lambda
a_{l^{\prime},l}B+\lambda^{2}a_{l^{\prime},l}^{2}Q)}\right],$ (36)

where $a_{l^{\prime},l}:=a_{l^{\prime}}-a_{l}$. “Taking out” of the second
exponential the Hamiltonian $H_{B}$, we obtain,

$\displaystyle{\rm Tr}_{B}\left[e^{-(\beta-u)H_{B}}e^{-u(H_{B}+\lambda
a_{l^{\prime},l}B+\lambda^{2}a_{l^{\prime},l}^{2}Q)}\right]$
$\displaystyle=e^{-u(\lambda^{2}a_{l^{\prime},l}^{2}Q)}{\rm
Tr}_{B}\left[e^{-\beta H_{B}}e^{-\lambda a_{l^{\prime},l}{\cal
T}\int_{0}^{u}ds\tilde{B}(s)}\right],$ (37)

with
$\tilde{B}(s):=e^{sH_{B}}Be^{-sH_{B}}=\sum_{k}g_{k}(e^{s\omega_{k}}b_{k}^{{\dagger}}+e^{-s\omega_{k}}b_{k})$.
After splitting $e^{-\lambda a_{l^{\prime},l}{\cal
T}\int_{0}^{u}ds\tilde{B}(s)}$ in two, we obtain

$e^{-\lambda a_{l^{\prime},l}{\cal
T}\int_{0}^{u}ds\tilde{B}(s)}=e^{-\lambda^{2}a_{l^{\prime},l}^{2}f(u)}e^{\lambda
a_{l^{\prime},l}{\cal B}(-u)}e^{-\lambda a_{l^{\prime},l}{\cal
B}^{{\dagger}}(u)},$ (38)

with $f(u):=-uQ+\sum_{k}\frac{g_{k}^{2}}{\omega_{k}^{2}}(e^{u\omega_{k}}-1)$,
${\cal B}(-u):=\sum_{k}\frac{g_{k}}{\omega_{k}}(e^{-u\omega_{k}}-1)b_{k}$, and
${\cal
B}^{{\dagger}}(u):=\sum_{k}\frac{g_{k}}{\omega_{k}}(e^{u\omega_{k}}-1)b_{k}^{{\dagger}}$.
A simple way to conclude the computation of the trace is using the normal
order representation of the operators. Using the following formula
Louisell_1964

$e^{\alpha
a^{{\dagger}}a}=\sum_{n=0}^{\infty}\frac{(e^{\alpha}-1)^{n}}{n!}(a^{{\dagger}})^{n}a^{n},$
(39)

where $\alpha$ is a c-number and $a^{{\dagger}}$, $a$ are bosonic operators,
we obtain,

$\displaystyle{\rm Tr}_{B}\left[e^{-(\beta-u)H_{B}}e^{-u(H_{B}+\lambda
a_{l^{\prime},l}B+\lambda^{2}a_{l^{\prime},l}^{2}Q)}\right]$
$\displaystyle\hskip
56.9055pt=e^{-u(\lambda^{2}a_{l^{\prime},l}^{2}Q)}e^{-\lambda^{2}a_{l^{\prime},l}^{2}f(u)}{\rm
Tr}_{B}\left[e^{-\lambda a_{l^{\prime},l}{\cal
B}^{{\dagger}}(u)}\left(\otimes_{k}\sum_{n=0}^{\infty}\frac{(e^{-\beta\omega_{k}}-1)^{n}}{n!}(a_{k}^{{\dagger}})^{n}a_{k}^{n}\right)e^{\lambda
a_{l^{\prime},l}{\cal B}(-u)}\right]$ $\displaystyle\hskip
56.9055pt=e^{-\lambda^{2}a_{l^{\prime},l}^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}^{2}}(e^{u\omega_{k}}-1)}\prod_{k}\left(\int_{\mathbb{R}^{2}}\frac{{\rm
d}\alpha_{k}^{2}}{\pi}e^{-\lambda
a_{l^{\prime},l}\frac{g_{k}}{\omega_{k}}(e^{u\omega_{k}}-1)\alpha_{k}^{*}}\sum_{n=0}^{\infty}\frac{(e^{-\beta\omega_{k}}-1)^{n}}{n!}|\alpha_{k}|^{2n}e^{\lambda
a_{l^{\prime},l}\frac{g_{k}}{\omega_{k}}(e^{-u\omega_{k}}-1)\alpha_{k}}\right)$
(40) $\displaystyle\hskip
56.9055pt=e^{-\lambda^{2}a_{l^{\prime},l}^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}^{2}}(e^{u\omega_{k}}-1)}\prod_{k}\left(\int_{\mathbb{R}^{2}}\frac{{\rm
d}\alpha_{k}^{2}}{\pi}e^{-\lambda
a_{l^{\prime},l}\frac{g_{k}}{\omega_{k}}(e^{u\omega_{k}}-1)\alpha_{k}^{*}}e^{(e^{-\beta\omega_{k}}-1)|\alpha_{k}|^{2}}e^{\lambda
a_{l^{\prime},l}\frac{g_{k}}{\omega_{k}}(e^{-u\omega_{k}}-1)\alpha_{k}}\right)$
$\displaystyle\hskip
56.9055pt=e^{-\lambda^{2}a_{l^{\prime},l}^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}^{2}}(e^{u\omega_{k}}-1)}\prod_{k}\left(\int_{\mathbb{R}^{2}}\frac{{\rm
d}x_{k}{\rm d}y_{k}}{\pi}e^{-\lambda
a_{l^{\prime},l}\frac{g_{k}}{\omega_{k}}(e^{u\omega_{k}}-1)(x_{k}-iy_{k})+(e^{-\beta\omega_{k}}-1)(x_{k}^{2}+y_{k}^{2})+\lambda
a_{l^{\prime},l}\frac{g_{k}}{\omega_{k}}(e^{-u\omega_{k}}-1)(x_{k}+iy_{k})}\right)$
$\displaystyle\hskip
56.9055pt=e^{-\lambda^{2}a_{l^{\prime},l}^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}^{2}}(e^{u\omega_{k}}-1)}\prod_{k}\left(\int_{\mathbb{R}^{2}}\frac{{\rm
d}x_{k}{\rm
d}y_{k}}{\pi}e^{(e^{-\beta\omega_{k}}-1)(x_{k}^{2}+y_{k}^{2})+\lambda
a_{l^{\prime},l}\frac{g_{k}}{\omega_{k}}(e^{-u\omega_{k}}-e^{u\omega_{k}})x_{k}+i\lambda
a_{l^{\prime},l}\frac{g_{k}}{\omega_{k}}(e^{-u\omega_{k}}+e^{u\omega_{k}}-2)y_{k}}\right)$
$\displaystyle\hskip
56.9055pt=e^{-\lambda^{2}a_{l^{\prime},l}^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}^{2}}(e^{u\omega_{k}}-1)}\prod_{k}\left(\frac{1}{1-e^{-\omega_{k}\beta}}\right)e^{-\lambda^{2}a_{l^{\prime},l}^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}^{2}}\frac{2-e^{u\omega_{k}}-e^{-u\omega_{k}}}{1-e^{-\omega_{k}\beta}}}$
$\displaystyle\hskip
56.9055pt=Z_{B}e^{-\lambda^{2}a_{l^{\prime},l}^{2}\sum_{k}\frac{g_{k}^{2}}{\omega_{k}^{2}}(e^{u\omega_{k}}-1)\left(1-\frac{1-e^{-u\omega_{k}}}{1-e^{-\omega_{k}\beta}}\right)}$
$\displaystyle\hskip
56.9055pt=Z_{B}e^{-\lambda^{2}a_{l^{\prime},l}^{2}\int_{0}^{\infty}\frac{J(\omega)}{\omega^{2}}(e^{u\omega}-1)\left(1-\frac{1-e^{-u\omega}}{1-e^{-\omega\beta}}\right)},$
(41)

where $Z_{B}:=\Pi_{k}\left(\frac{1}{1-e^{-\omega_{k}\beta}}\right)={\rm
Tr}_{B}\left(e^{-\beta H_{B}}\right)$. In line (40), we used the property of
the coherent states $|\alpha_{k}\rangle$ which form an over-complete basis of
the Hilbert space of the $k$-mode, so that the trace of any operator ${\cal
O}_{k}$ acting in that Hilbert space can be computed as ${\rm Tr}_{k}({\cal
O}_{k})=\int_{\mathbb{R}^{2}}\frac{{\rm
d}\alpha_{k}^{2}}{\pi}\langle\alpha_{k}|{\cal O}_{k}|\alpha_{k}\rangle$.
Additionally, we expressed in the last lines the complex variable $\alpha_{k}$
explicitly in term of its real and imaginary part $\alpha_{k}=x_{k}+iy_{k}$.

##  Appendix B Important property of $f_{l,l^{\prime}}(\beta)$

By applying the change of variable $v=\beta-u$, one can re-write
$f_{l,l^{\prime}}(\beta)$ in the form

$\displaystyle f_{l,l^{\prime}}(\beta)$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}due^{u\omega_{l,l^{\prime}}}e^{-\lambda^{2}a_{l^{\prime},l}^{2}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}\frac{e^{\omega\beta/2}+e^{-\omega\beta/2}-e^{\omega(u-\beta/2)}-e^{-\omega(u-\beta/2)}}{e^{\omega\beta/2}-e^{-\omega\beta/2}}}$
(42) $\displaystyle=$ $\displaystyle
e^{\beta\omega_{l,l^{\prime}}}\int_{0}^{\beta}dve^{-v\omega_{l,l^{\prime}}}e^{-\lambda^{2}a_{l^{\prime},l}^{2}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}\frac{e^{\omega\beta/2}+e^{-\omega\beta/2}-e^{\omega(v-\beta/2)}-e^{-\omega(v-\beta/2)}}{e^{\omega\beta/2}-e^{-\omega\beta/2}}}.$

Then, one can easily verify that $p_{l}^{\rm
ss}f_{l,l^{\prime}}(\beta)=p_{l^{\prime}}^{\rm ss}f_{l^{\prime},l}(\beta)$.
This implies in particular that $\langle a_{l}|\rho_{S}^{\rm
ss}|a_{l^{\prime}}\rangle^{*}=\langle a_{l^{\prime}}|\rho_{S}^{\rm
ss}|a_{l}\rangle$, as required by the Hermicity of $\rho_{S}^{\rm ss}$.

##  Appendix C Without initial renormalization

As commented above, one can choose to perform the same derivation starting
from the “natural” Hamiltonian

${\cal H}_{SB}=H_{S}+H_{B}+\lambda AB,$ (43)

instead of the renormalized one $H_{SB}$ defined in the main text. Proceeding
as previously, we have

$\displaystyle e^{-\beta{\cal H}_{SB}}=e^{-\beta(H_{B}+\lambda AB)}e^{-{\cal
T}\int_{0}^{\beta}du\tilde{\cal H}_{S}(u)},$ (44)

with

$\tilde{\cal H}_{S}(u):=e^{u(H_{B}+\lambda AB)}H_{S}e^{-u(H_{B}+\lambda AB)}.$
(45)

We now have

$\displaystyle H_{B}+\lambda AB$ $\displaystyle=$
$\displaystyle\sum_{k}\omega_{k}{\cal D}_{k}b_{k}^{{\dagger}}b_{k}{\cal
D}_{k}^{{\dagger}}-\lambda^{2}QA^{2}$ (46) $\displaystyle=$
$\displaystyle\sum_{l}|a_{l}\rangle\langle a_{l}|{\cal H}_{B,l},$

with ${\cal H}_{B,l}:=H_{B,l}-\lambda^{2}a_{l}^{2}Q=H_{B}+\lambda a_{l}B$,
leading to

$\displaystyle e^{u(H_{B}+\lambda AB)}$ $\displaystyle=$ $\displaystyle
e^{u\sum_{l}|a_{l}\rangle\langle a_{l}|{\cal H}_{B,l}}$ (47) $\displaystyle=$
$\displaystyle\sum_{l}|a_{l}\rangle\langle a_{l}|e^{u{\cal H}_{B,l}}.$

Then,

$\displaystyle\tilde{H}_{S}(u)$ $\displaystyle=$
$\displaystyle\sum_{l,l^{\prime}}|a_{l}\rangle\langle
a_{l}|e^{uH_{B,l}}H_{S}|a_{l^{\prime}}\rangle\langle
a_{l^{\prime}}|e^{-uH_{B,l^{\prime}}}$ (48) $\displaystyle=$ $\displaystyle
H_{S}^{\rm pop}+\tilde{H}_{S}^{\rm coh}(u),$

where ${\cal H}_{S}^{\rm pop}=H_{S}^{\rm
pop}=\sum_{l}h_{l}|a_{l}\rangle\langle a_{l}|$ and $\tilde{\cal H}_{S}^{\rm
coh}(u):=\sum_{l\neq l^{\prime}}h_{l,l^{\prime}}|a_{l}\rangle\langle
a_{l^{\prime}}|e^{u(H_{B,l}-\lambda^{2}a_{l}^{2}Q)}e^{-u(H_{B,l^{\prime}}-\lambda^{2}a_{l^{\prime}}^{2}Q)}=\\\
\sum_{l\neq l^{\prime}}h_{l,l^{\prime}}|a_{l}\rangle\langle
a_{l^{\prime}}|e^{u{\cal H}_{B,l}}e^{-u{\cal H}_{B,l^{\prime}}}$. Again,
similarly as previously, we obtain,

$\displaystyle e^{-{\cal T}\int_{0}^{\beta}du\tilde{\cal H}_{S}(u)}$
$\displaystyle=$ $\displaystyle e^{-\beta H_{S}^{\rm pop}}e^{-{\cal
T}\int_{0}^{\beta}du\accentset{\approx}{\cal H}_{S}^{\rm coh}(u)},$ (49)

with

$\displaystyle\accentset{\approx}{\cal H}_{S}^{\rm coh}(u)$ $\displaystyle:=$
$\displaystyle e^{uH_{S}^{\rm pop}}\tilde{\cal H}_{S}^{\rm
coh}(u)e^{-uH_{S}^{\rm pop}}$ (50) $\displaystyle=$ $\displaystyle\sum_{l\neq
l^{\prime}}h_{l,l^{\prime}}e^{u\omega_{l,l^{\prime}}}|a_{l}\rangle\langle
a_{l^{\prime}}|e^{u{\cal H}_{B,l}}e^{-u{\cal H}_{B,l^{\prime}}},$

arriving at

$\displaystyle\rho_{SB}^{\rm th}$ $\displaystyle=$ $\displaystyle
Z_{SB}^{-1}\sum_{l}e^{-\beta h_{l}}|a_{l}\rangle\langle a_{l}|e^{-\beta{\cal
H}_{B,l}}e^{-{\cal T}\int_{0}^{\beta}du\accentset{\approx}{\cal H}_{S}^{\rm
coh}(u)}.$

Finally, the main change is that we are led to compute ${\rm
Tr}_{B}\left[e^{-\beta{\cal H}_{B,l}}\right]$ for the zeroth order, and ${\rm
Tr}_{B}\left[e^{-\beta{\cal H}_{B,l}}e^{u{\cal H}_{B,l}}e^{-u{\cal
H}_{B,l^{\prime}}}\right]$ for the second order, instead of ${\rm
Tr}_{B}\left[e^{-\beta H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l^{\prime}}}\right]$. We
have,

$\displaystyle{\rm Tr}_{B}\left[e^{-\beta{\cal H}_{B,l}}e^{u{\cal
H}_{B,l}}e^{-u{\cal H}_{B,l^{\prime}}}\right]$ (52) $\displaystyle=$
$\displaystyle{\rm
Tr}_{B}\left[e^{-\beta(H_{B,l}-\lambda^{2}a_{l}^{2}Q)}e^{u(H_{B,l}-\lambda^{2}a_{l}^{2}Q)}e^{-u(H_{B,l^{\prime}}-\lambda^{2}a_{l^{\prime}}^{2}Q)}\right]$
$\displaystyle=$ $\displaystyle{\rm
Tr}_{B}\left[e^{-(\beta-u)H_{B}}e^{-u(H_{B}+\lambda
a_{l^{\prime},l}B+\lambda^{2}a_{l^{\prime},l}^{2}Q)}\right]$
$\displaystyle\times e^{\lambda^{2}Q[\beta
a_{l}^{2}+u(a_{l^{\prime}}^{2}-a_{l}^{2})]}.$

With that, we obtain

$\displaystyle\rho_{S}^{\rm ss}$ $\displaystyle=$ $\displaystyle\sum_{l}{\tt
p}_{l}^{\rm ss}|a_{l}\rangle\langle a_{l}|-\sum_{l,l^{\prime};l\neq
l^{\prime}}{\tt p}_{l}^{\rm ss}h_{l,l^{\prime}}{\tt
f}_{l,l^{\prime}}(\beta)|a_{l}\rangle\langle a_{l^{\prime}}|,$

where ${\tt p}_{l}^{\rm ss}:=e^{-\beta(h_{l}-a_{l}^{2}\lambda^{2}Q)}/{\cal
Z}_{S}^{\rm ss}=e^{-\beta{\tt h}_{l}}/{\cal Z}_{S}^{\rm ss}$ is the
renormalized population, with ${\tt h}_{l}:=h_{l}-a_{l}^{2}\lambda^{2}Q$ the
renormalized “pseudo-energies” (diagonal elements of $H_{S}$ in the eigenbasis
of $A$), ${\cal Z}^{\rm ss}:=\sum_{l}e^{-\beta{\tt h}_{l}}$, and

$\displaystyle{\tt f}_{l,l^{\prime}}(\beta):=\int_{0}^{\beta}du$
$\displaystyle\times
e^{u\bar{\omega}_{l,l^{\prime}}}e^{-\lambda^{2}a_{l^{\prime},l}^{2}\int_{0}^{\infty}d\omega\frac{J(\omega)}{\omega^{2}}\frac{e^{\omega\beta/2}+e^{-\omega\beta/2}-e^{\omega(u-\beta/2)}-e^{-\omega(u-\beta/2)}}{e^{\omega\beta/2}-e^{-\omega\beta/2}}},$

with $\bar{\omega}_{l,l^{\prime}}:={\tt h}_{l}-{\tt h}_{l^{\prime}}$.
Additionally, one can also verify the identity ${\tt p}_{l}^{\rm ss}{\tt
f}_{l,l^{\prime}}(\beta)={\tt p}_{l^{\prime}}^{\rm ss}{\tt
f}_{l^{\prime},l}(\beta)$.

Conclusion.– If we do not renormalize the Hamiltonian initially, the
expressions are un-changed up to the substitution of $h_{l}$ by ${\tt h}_{l}$.
The renormalization has to happen, either initially, either finally. Note
however that $h_{l,l^{\prime}}$ is not changed. One must be aware of these
differences of choice especially when defining the strong coupling regime.

###  C.1 Approximation of ${\tt f}_{l,l^{\prime}}(\beta)$

Similarly as for ${f}_{l,l^{\prime}}(\beta)$, when the bath spectral density
is such that $J(\omega)$ vanishes for $\omega\geq\omega_{c}$, where
$\omega_{c}\leq\beta^{-1}$, we have

$\displaystyle{\tt f}_{l,l^{\prime}}(\beta)$ $\displaystyle\simeq$
$\displaystyle\int_{0}^{\beta}due^{u\bar{\omega}_{l,l^{\prime}}}e^{-\lambda^{2}a_{l^{\prime},l}^{2}u\left(1-\frac{u}{\beta}\right)Q}$
$\displaystyle=\frac{1}{\lambda|a_{l^{\prime},l}|}\sqrt{\frac{\beta}{Q}}\Bigg{\\{}{\rm
DF}\left[\frac{1}{2\lambda|a_{l^{\prime},l}|}\sqrt{\frac{\beta}{Q}}(\lambda^{2}a_{l^{\prime},l}^{2}Q-\bar{\omega}_{l,l^{\prime}})\right]$
$\displaystyle+e^{\beta\bar{\omega}_{l,l^{\prime}}}{\rm
DF}\left[\frac{1}{2\lambda|a_{l^{\prime},l}|}\sqrt{\frac{\beta}{Q}}(\lambda^{2}a_{l^{\prime},l}^{2}Q+\bar{\omega}_{l,l^{\prime}})\right]\Bigg{\\}}.$

Note that
$\lambda^{2}a_{l^{\prime},l}^{2}Q-\bar{\omega}_{l,l^{\prime}}=2\lambda^{2}Qa_{l}(a_{l}-a_{l^{\prime}})-\omega_{l,l^{\prime}}$
and
$\lambda^{2}a_{l^{\prime},l}^{2}Q+\bar{\omega}_{l,l^{\prime}}=2\lambda^{2}Qa_{l^{\prime}}(a_{l^{\prime}}-a_{l})+\omega_{l,l^{\prime}}$.
Thus, in the strong coupling regime when $\lambda^{2}Q\gg{\rm
max}_{l}~{}|h_{l}|$ and $\lambda^{2}Q\beta\gg 1$, it is still legitimate to
approximate the function $DF(x)$ by $1/2x$, which gives,

$\displaystyle{\tt f}_{l,l^{\prime}}(\beta)$ $\displaystyle=$
$\displaystyle\frac{1}{2\lambda^{2}Qa_{l,l^{\prime}}}\left(\frac{1}{a_{l}}-\frac{e^{\bar{\omega}_{l,l^{\prime}}\beta}}{a_{l^{\prime}}}\right)$
$\displaystyle+\frac{\omega_{l,l^{\prime}}}{4\lambda^{4}Q^{2}a_{l,l^{\prime}}^{2}}\left(\frac{1}{a_{l}^{2}}-\frac{e^{\bar{\omega}_{l,l^{\prime}}\beta}}{a_{l^{\prime}}^{2}}\right)+{\cal
O}[(\lambda^{2}Q/\omega_{l,l^{\prime}})^{-3}].$

##  Appendix D Comparison of the general expressions

In this section, we compare our general result, Eq.(20) and (21) of the main
text, with the general expression for the steady state coherences obtained in
Trushechkin_2021 . To make the comparison easier, we re-write in the following
the results of Trushechkin_2021  using the notations we have been using here.
Additionally, in order to simplify the presentation, we consider a simple
coupling of the form $V_{I}=\lambda AB$ instead of the general one
$V_{I}=\sum_{\alpha}A_{\alpha}B_{\alpha}$ considered in Trushechkin_2021 .
With that, in the eigenbasis $\\{|a_{l}\rangle\\}_{l}$ of $A$, the result of
Trushechkin_2021  takes the form

$\displaystyle\rho_{l,l^{\prime}}^{\rm ss,ME}$ $\displaystyle:=$
$\displaystyle\langle a_{l}|\rho_{S}^{\rm ss,ME}|a_{l^{\prime}}\rangle$
$\displaystyle=$ $\displaystyle h_{l,l^{\prime}}\Big{[}i{\tt p}_{l}^{\rm
ss}\int_{0}^{\infty}d\tau
e^{-\lambda^{2}a_{l^{\prime},l}^{2}[G^{*}(\tau)-i\tau
Q]}e^{-i\bar{\omega}_{l,l^{\prime}}\tau}$ $\displaystyle-i{\tt
p}_{l^{\prime}}^{\rm ss}\int_{0}^{\infty}d\tau
e^{-\lambda^{2}a_{l,l^{\prime}}^{2}[G(\tau)+i\tau
Q]}e^{-i\bar{\omega}_{l,l^{\prime}}\tau}\Big{]}$

with ${\tt p}_{l}^{\rm ss}:=e^{-\beta(h_{l}-a_{l}^{2}\lambda^{2}Q)}/{\cal
Z}_{S}^{\rm ss}=e^{-\beta{\tt h}_{l}}/{\cal Z}_{S}^{\rm ss}$ is the
renormalized population introduced in Appendix C, with ${\tt
h}_{l}:=h_{l}-a_{l}^{2}\lambda^{2}Q$ the renormalized “pseudo-energies”, and
${\cal Z}^{\rm ss}:=\sum_{l}e^{-\beta{\tt h}_{l}}$. Additionally, we defined
$G(\tau):=\int_{0}^{\tau}ds_{1}\int_{0}^{s_{1}}ds_{2}c_{B}(s_{2})$, and

$\displaystyle c_{B}(s)={\rm Tr}_{B}[e^{iH_{B}s}Be^{-iH_{B}s}B\rho_{B}^{\rm
th}].$ (58)

The superscript “ME” refers to the master equation nature of the derivation
Trushechkin_2021 .

Our result can be expressed as (using the results of the derivation with no
initial renormalization of the Hamiltonian, detailed in section C of this
Supplementary Material, since it is the choice made in Trushechkin_2021 ),

$\displaystyle\rho_{l,l^{\prime}}^{\rm ss}$ $\displaystyle:=$
$\displaystyle\langle a_{l}|\rho_{S}^{\rm ss}|a_{l^{\prime}}\rangle$ (59)
$\displaystyle=$ $\displaystyle-h_{l,l^{\prime}}{\tt p}_{l}^{\rm
ss}\int_{0}^{\beta}due^{-\lambda^{2}a_{l,l^{\prime}}^{2}[G(-iu)+uQ]}e^{u\bar{\omega}_{l,l^{\prime}}}$
$\displaystyle=$ $\displaystyle-\frac{1}{2}h_{l,l^{\prime}}{\tt p}_{l}^{\rm
ss}\int_{0}^{\beta}due^{-\lambda^{2}a_{l,l^{\prime}}^{2}[G(-iu)+uQ]}e^{u\bar{\omega}_{l,l^{\prime}}}$
$\displaystyle-\frac{1}{2}h_{l,l^{\prime}}{\tt p}_{l^{\prime}}^{\rm
ss}\int_{0}^{\beta}due^{-\lambda^{2}a_{l^{\prime},l}^{2}[G(-iu)+uQ]}e^{u\bar{\omega}_{l^{\prime},l}}$

where we used in the second line the identity ${\tt p}_{l}^{\rm ss}{\tt
f}_{l,l^{\prime}}(\beta)={\tt p}_{l^{\prime}}^{\rm ss}{\tt
f}_{l^{\prime},l}(\beta)$. Note also that
$G(-iu)=\int_{0}^{-iu}ds_{1}\int_{0}^{s_{1}}ds_{2}c_{B}(s_{2})=-\int_{0}^{u}ds_{1}\int_{0}^{s_{1}}ds_{2}c_{B}(-is_{2})$,
and

$\displaystyle c_{B}(-is)={\rm Tr}_{B}[e^{H_{B}s}Be^{-H_{B}s}B\rho_{B}^{\rm
th}].$ (60)

Substituting the variable of integration $u$ by $i\tau$, we obtain,

$\displaystyle\rho_{l,l^{\prime}}^{\rm ss}$ $\displaystyle=$
$\displaystyle\frac{h_{l,l^{\prime}}}{2}\Big{[}-i{\tt p}_{l}^{\rm
ss}\int_{0\rightarrow-i\beta}d\tau
e^{-\lambda^{2}a_{l,l^{\prime}}^{2}[G(\tau)+i\tau
Q]}e^{i\bar{\omega}_{l,l^{\prime}}\tau}$ (61) $\displaystyle-i{\tt
p}_{l^{\prime}}^{\rm ss}\int_{0\rightarrow-i\beta}d\tau
e^{-\lambda^{2}a_{l^{\prime},l}^{2}[G(\tau)+i\tau
Q]}e^{i\bar{\omega}_{l^{\prime},l}\tau}\Big{]},$

where $0\rightarrow-i\beta$ denotes the integration along the imaginary axis,
from 0 to $-i\beta$. There are three differences with respect to
$\rho_{l,l^{\prime}}^{\rm ss,ME}$. Firstly, we have a factor $1/2$, then, the
first term in the square brackets is the complex conjugate of the first one in
$\rho_{l,l^{\prime}}^{\rm ss,ME}$, and finally, the integration are from $0$
to $-i\beta$ along the imaginary axis while in $\rho_{l,l^{\prime}}^{\rm
ss,ME}$ integrations run from $0$ to $+\infty$.

As shown in section III B of the main text, the high temperature prediction of
both expressions are in excellent agreement. But beyond this regime, are these
two expressions still equivalent or equal? We left the question open for now,
as deeper analysis and alternative methods might be required. Note that the
discrepancy appearing beyond the regime of validity in Fig. 3 might have been
generated during the application of the high temperature approximation, so
that such discrepancy is not necessarily an evidence that (59) and (61) are
not equal.