# Steady State in Ultrastrong Coupling Regime:
Expansion and First Orders

Camille Lombard Latune

<figure><img src="image_1.png" /></figure>

► To cite this version:

Camille Lombard Latune. Steady State in Ultrastrong Coupling Regime: Expansion and First Orders. Fields and Quanta, 2022, 11 (1), pp.53 - 71. 10.12743/quanta.v11i1.167. hal-04924334

HAL Id: hal-04924334

https://cnrs.hal.science/hal-04924334v1

Submitted on 31 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Steady State in Ultrastrong Coupling Regime: Expansion and First Orders

Camille Lombard Latune

Laboratoire de Physique, École Normale Supérieure de Lyon, France. E-mail: cami-address@hotmail.com
Quantum Research Group, University of KwaZulu-Natal, Durban, South Africa
National Institute for Theoretical Physics (NITheP), Durban, South Africa

Editors: Vinayak Jagadish, Anton Trushechkin, James Cresser & Danko D. Georgiev

Article history: Submitted on September 29, 2021; Accepted on October 28, 2022; Published on November 20, 2022.

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 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 thermal state, largely credited to be the steady state. Here, we present analytical expressions of corrective terms to the ultrastrong coupling limit of the mean force Gibbs state, which has been recently derived. We find that the first order term precisely coincides with the first order correction obtained from a dynamical approach—master equation in the strong-decoherence regime. This strengthens the identification of the reduced steady state with the mean force Gibbs state. Additionally, we also compare our expressions with another recent result obtained from a high temperature expansion of the mean force Gibbs state. We observe numerically a good agreement for ultra strong coupling as well as for high temperatures. This confirms the validity of all these results. In particular, we show that, in term of coherences, all three results allow one to sketch the transition from ultrastrong coupling to weak coupling.

Quanta 2022; 11: 53–71.

## 1 Introduction

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 [1–8], quantum thermodynamics [9–15], quantum sensing [16–18], as well as understanding better the underlying physics of some essential biological functions [19–22]. 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 [23–26]. While the properties of the most general open quantum evolutions have been known for long time thanks to the seminal paper by Sudarshan, Mathews and Rau [27], their precise time evolutions and steady states are still a challenge of the theory of open quantum system [28].

To obtain some information about strongly coupled steady states, several strategies have been developed, including embedding techniques like reaction coordinate [10, 23, 29, 30] and pseudo-mode [31–34], or numerical techniques (Hierarchical Equation of Motion) [22, 35, 36].

<figure><img src="image_3.png" /></figure>

This is an open access article distributed under the terms of the Creative Commons Attribution License CC-BY-3.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
One alternative strategy consists in focusing directly on the steady state without going through the description of the whole dynamics. In this perspective, the global steady state of the system and bath is expected to be the global system-bath thermal state at the bath temperature [37–43]. The steady state of the system is then given by tracing out the bath, which is often referred to as the mean force Gibbs state [25, 43, 44]. Such partial trace is usually very challenging, but can be done at least approximately assuming for instance a weak coupling [24, 25, 40, 45]. Another 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 [25], 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 close to the form of the steady state obtained in [46, 47] using arguments from eisenselection [48], 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 [49, 50] for diagonalization of the quantum Rabi model as well as reminiscent of the polaron transformation [51, 52], we recover the infinite coupling limit of [25] 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 to two very recent derivations. The first one [44] was obtained from a master equation in the so-called strong-decoherence regime, a generalization of the ultra-strong coupling regime, and actually coincides with the first order expansion derived here. This confirms that a system interacting with a thermal bath does converge, at least up to first order, to the mean force Gibbs state even in the ultra-strong coupling regime. The second recent derivation [53] consists in a high temperature expansion of the mean force Gibbs state, a generalization of a derivation introduced in [54]. The comparison with this result is essentially numeric, and we find an overall good agreement in the expected regime of validity.

Additionally, we show that the provided corrective terms allow to sketch the transition from the ultrastrong coupling limit to the weak coupling limit. Finally, we provide some higher order corrections in Section 7.3.

## 2 Mean force Gibbs state

We consider the following total Hamiltonian

$$ \mathcal{H}_{SB} = H_S + H_B + \lambda H_I, \quad (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 [28], 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 [25]. The quantity $Q$, the “re-organization energy” [25,55], 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 often 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. \quad (2) $$

Note the difference of notation between the “natural” Hamiltonian $\mathcal{H}_{SB}$ (1) and the renormalized Hamiltonian $H_{SB}$ (2). As a matter of completeness, we also mention in Section 7.5 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 expected to be the global thermal state (assuming $[H_S, A] \neq 0$) [37–43]

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

where $Z_{SB} := \text{Tr}_{SB}(e^{-\beta H_{SB}})$ is the partition function. Then, the corresponding reduced steady state is

$$ \rho_S^{ss} := Z_{SB}^{-1} \text{Tr}_B[e^{-\beta H_{SB}}], \quad (4) $$

which is called mean force Gibbs state [25, 43, 44]. We are going to approximate this state in the regime of ultrastrong coupling, when $\lambda$ becomes larger than the energy scale of $H_S$.

## 3 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 expo-
nential using the usual identities [56]. It leads to

$$e^{-\beta H_{SB}} = e^{-\beta(H_B + \lambda AB + \lambda^2 QA^2)} e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S(u)}, \quad (5)$$

with $\mathcal{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)}. \quad (6)$$

Now, by noticing that $H_B + \lambda AB + \lambda^2 QA^2 = \sum_k \omega_k \mathcal{D}_k b_k^\dagger b_k \mathcal{D}_k^\dagger$, where $\mathcal{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 (assumed to be non-degenerate for simplicity) and corresponding eigenvectors of the observable $A$, and by $\mathcal{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

$$\begin{aligned}
H_B + \lambda AB + \lambda^2 QA^2 &= \sum_k \omega_k \mathcal{D}_k b_k^\dagger b_k \mathcal{D}_k^\dagger \\
&= \sum_k \omega_k \sum_l |a_l\rangle\langle a_l| \mathcal{D}_k b_k^\dagger b_k \mathcal{D}_k^\dagger \\
&= \sum_k \omega_k \sum_l |a_l\rangle\langle a_l| \mathcal{D}_{k,l} b_{k,l}^\dagger b_k \mathcal{D}_{k,l}^\dagger \\
&= \sum_l |a_l\rangle\langle a_l| H_{B,l},
\end{aligned} \quad (7)$$

where we defined the “displaced bath” $H_{B,l} := \sum_k \omega_k \mathcal{D}_{k,l} b_k^\dagger b_k \mathcal{D}_{k,l}^\dagger = H_B + \lambda a_l B + \lambda^2 a_l^2 Q$, reminiscent of the displaced oscillator picture [49, 50]. From there we obtain,

$$\begin{aligned}
e^{u(H_B + \lambda AB + \lambda^2 QA^2)} &= e^{u \sum_l |a_l\rangle\langle a_l| H_{B,l}} \\
&= \sum_l |a_l\rangle\langle a_l| e^{u H_{B,l}},
\end{aligned} \quad (8)$$

so that

$$\begin{aligned}
\tilde{H}_S(u) &= \sum_{l,l'} |a_l\rangle\langle a_l| e^{u H_{B,l}} H_S |a_{l'}\rangle\langle a_{l'}| e^{-u H_{B,l'}} \\
&= H_S^{\text{pop}} + \tilde{H}_S^{\text{coh}}(u),
\end{aligned} \quad (9)$$

where $H_S^{\text{pop}} := \sum_l h_l |a_l\rangle\langle a_l|$ and $\tilde{H}_S^{\text{coh}}(u) := \sum_{l \neq l'} h_{l,l'} |a_{l'}\rangle\langle a_{l'}| e^{u H_{B,l'}} e^{-u H_{B,l'}}$, defining $h_l := \langle a_l | H_S | a_l \rangle$, and $h_{l,l'} := \langle a_l | H_S | a_{l'} \rangle$.

Taking out $H_S^{\text{pop}}$ from the exponential $e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S(u)$, we obtain,

$$e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S(u)} = e^{-\beta H_S^{\text{pop}}} e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S^{\text{coh}}(u)}, \quad (10)$$

with

$$\begin{aligned}
\tilde{H}_S^{\text{coh}}(u) &= e^{u H_S^{\text{pop}}} \tilde{H}_S^{\text{coh}}(u) e^{-u H_S^{\text{pop}}} \\
&= \sum_{l \neq l'} h_{l,l'} e^{u \omega_{l,l'}} |a_l\rangle\langle a_l| e^{u H_{B,l}} e^{-u H_{B,l'}},
\end{aligned} \quad (11)$$

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

$$\begin{aligned}
\rho_{SB}^{\text{th}} &= Z_{SB}^{-1} \sum_l |a_l\rangle\langle a_l| e^{-\beta H_{B,l}} e^{-\beta H_S^{\text{pop}}} e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S^{\text{coh}}(u)} \\
&= Z_{SB}^{-1} \sum_l e^{-\beta h_l} |a_l\rangle\langle a_l| e^{-\beta H_{B,l}} e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S^{\text{coh}}(u)}.
\end{aligned} \quad (12)$$

This leads to the following expression for the mean force Gibbs state,

$$\rho_S^{\text{ss}} = Z_{SB}^{-1} \sum_l e^{-\beta h_l} |a_l\rangle\langle a_l| \text{Tr}_B \left[ e^{-\beta H_{B,l}} e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S^{\text{coh}}(u)} \right]. \quad (13)$$

The toughest part is, as expected,

$$\text{Tr}_B \left[ e^{-\beta H_{B,l}} e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S^{\text{coh}}(u)} \right] = \sum_{n=0}^{\infty} (-1)^n \int_0^\beta du_1 \int_0^{u_1} du_2 \dots \int_0^{u_{n-1}} du_n \text{Tr}_B \left[ e^{-\beta H_{B,l}} \tilde{H}_S^{\text{coh}}(u_1) \tilde{H}_S^{\text{coh}}(u_2) \dots \tilde{H}_S^{\text{coh}}(u_n) \right]. \quad (14)$$

In the following, when we mention “first order term” or “higher order terms”, we refer to the terms appearing in the above expansion (13). So far, we made no approximation, making (13) 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 infinite coupling limit.

## 3.1 Recovering the infinite coupling limit

From expression (13), it might not be obvious how one recovers the infinite coupling limit. It actually comes from the terms $\tilde{H}_S^{\text{coh}}(u)$ which contain overlaps between displaced baths, $e^{uH_{B,l}} e^{-uH_{B,l'}}$. The coupling strength between $S$ and each mode of the bath is given by $\lambda g_k$. Thus,
when $\lambda$ goes to infinity, the displaced baths $H_{B,l}$ and $H_{B,l'}$ tend to be displaced infinitely far apart from each other (for $l \neq l'$). Consequently, the overlap and the expectation value $\text{Tr}_B[e^{-\beta H_{B,l}}e^{uH_{B,l'}}e^{-uH_{B,l''}}]$ tends to zero for increasing coupling strength.

Extending this reasoning to higher order terms, we see that all terms $\text{Tr}_B[e^{-\beta H_{B,l}} \tilde{H}_S^{\text{coh}}(u_1) \tilde{H}_S^{\text{coh}}(u_2) \dots \tilde{H}_S^{\text{coh}}(u_n)]$ contain multiple overlaps of different displaced baths, and therefore tends to zero as the coupling strength increases. In Section 7.3, we show numerically (for a two-level systems) that the higher order terms tend to zero as the coupling strength increases. Thus, in the infinite coupling limit, only the first term of the sum in (13) is different from zero. It leads to

$$ \rho_S^{ss,\infty} = Z_{SB}^{-1} \sum_l e^{-\beta h_l} |a_l\rangle\langle a_l| \text{Tr}_B[e^{-\beta H_{B,l}}] \\ = (Z_S^{ss})^{-1} \sum_l e^{-\beta h_l} |a_l\rangle\langle a_l|, \quad (15) $$

where $Z_S^{ss} := \frac{Z_{SB}}{Z_B}$ and $\text{Tr}_B[e^{-\beta H_{B,l}}] = \text{Tr}_B[e^{-\beta H_B}] := Z_B$ is the partition function of the uncoupled bath. Note that since $\rho_S^{ss,\infty}$ is a normalized state, we also have the identity $Z_S^{ss} = \sum_l e^{-\beta h_l}$. The above expression (15) is exactly equal to the one derived in [25], and recently shown to coincide with the steady state in the ultrastrong coupling limit [44].

## 3.2 First approximation

As seen in the previous section, the term $\tilde{H}_S^{\text{coh}}(u)$ contains overlaps of displaced baths, so that terms $\text{Tr}_B[e^{-\beta H_{B,l}} \tilde{H}_S^{\text{coh}}(u_1) \tilde{H}_S^{\text{coh}}(u_2) \dots \tilde{H}_S^{\text{coh}}(u_n)]$ 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,

$$ \rho_S^{ss} \approx Z_{SB}^{-1} \sum_l e^{-\beta h_l} |a_l\rangle\langle a_l| \left\{ \text{Tr}_B[e^{-\beta H_{B,l}}] - \text{Tr}_B \left[ e^{-\beta H_{B,l}} \int_0^\beta du \tilde{H}_S^{\text{coh}}(u) \right] \right\}. \quad (16) $$

The first term gives $\text{Tr}_B[e^{-\beta H_{B,l}}] = Z_B$, as already seen in the previous section. The second term gives,

$$ \text{Tr}_B \left[ e^{-\beta H_{B,l}} \int_0^\beta du \tilde{H}_S^{\text{coh}}(u) \right] = \sum_{l',l''} h_{l',l''} |a_{l'}\rangle\langle a_{l''}| \\ \times \int_0^\beta du e^{u\omega_{l',l''}} \text{Tr}_B \left[ e^{-\beta H_{B,l}} e^{uH_{B,l'}} e^{-uH_{B,l''}} \right]. \quad (17) $$

Since, when injected in (16), this expression will be multiplied by $|a_l\rangle\langle a_l|$ on the left-hand side, one only needs to

compute

$$ \sum_{l' \neq l} h_{l',l} |a_l\rangle\langle a_{l'}| \int_0^\beta du e^{u\omega_{l,l'}} \text{Tr}_B \left[ e^{-\beta H_{B,l}} e^{uH_{B,l}} e^{-uH_{B,l'}} \right], \quad (18) $$

where the index $l$ and $l'$ were made equal. After some manipulations, we can obtain the following expression for $\text{Tr}_B[e^{-\beta H_{B,l}} e^{uH_{B,l}} e^{-uH_{B,l'}}]$ (details provided in Section 7.1),

$$ \text{Tr}_B \left[ e^{-\beta H_{B,l}} e^{uH_{B,l}} e^{-uH_{B,l'}} \right] = Z_B e^{-\lambda^2 a_{l',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)}, \quad (19) $$

where $a_{l',l} := a_{l'} - a_l$. We finally obtain, up to second order,

$$ \rho_S^{ss} = \sum_l p_l^{ss} |a_l\rangle\langle a_l| - \sum_{l,l';l \neq l'} p_l^{ss} h_{l,l'} f_{l,l'}(\beta) |a_l\rangle\langle a_{l'}|, \quad (20) $$

with $p_l^{ss} := e^{-\beta h_l}/Z_S^{ss}$ and

$$ f_{l,l'}(\beta) = \int_0^\beta du e^{u\omega_{l,l'}} e^{-\lambda^2 a_{l',l}^2 \int_0^\infty d\omega \frac{J(\omega)}{\omega^2} (e^{u\omega}-1) \left(1-\frac{1-e^{-u\omega}}{1-e^{-\omega\beta}}\right)}. \quad (21) $$

Note that we have the following identity (see Section 7.2) $p_l^{ss} f_{l,l'}(\beta) = p_{l'}^{ss} f_{l',l}(\beta)$, implying $\langle a_l|\rho_S^{ss}|a_{l'}\rangle^* = \langle a_{l'}|\rho_S^{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 [44]. Finally, the expression obtained starting from the natural Hamiltonian (1) instead of the renormalized one (2) are the same as (20) and (21) but substituting the pseudo-energies $h_l$ by the renormalized ones $h_l - \lambda^2 a_l^2 Q$ (see Section 7.5).

## 3.3 Approximate expression of $f_{l,l'}(\beta)$

Depending on the bath spectral density, it might not be possible to obtain an exact analytical expression of $f_{l,l'}(\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'}(\beta)$ without even specifying the form of the spectral density. More precisely, we assume that $J(\omega)$ vanishes for $\omega \ge \omega_c$, where $\omega_c \le \beta^{-1}$. Then, for $\omega \in [0; \beta^{-1}]$, the factor $(e^{u\omega} - 1) (1 - \frac{1-e^{-u\omega}}{1-e^{-\omega\beta}})$ 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) + \mathcal{O}(\omega\beta), \quad (22) $$
which is actually a very good approximation as soon as $\omega\beta \le 1$. With that we obtain

$$e^{-\lambda^2 a_{l',l}^2} \int_0^\infty d\omega \frac{J(\omega)}{\omega^2} (e^{iu\omega} - 1) \left(1 - \frac{1-e^{-\omega u}}{1-e^{-\omega\beta}}\right) \approx e^{-\lambda^2 a_{l',l}^2} u(1-\frac{u}{\beta})Q, \quad (23)$$

and

$$f_{l,l'}(\beta) \approx \int_0^\beta du e^{u\omega_{l,l'}} e^{-\lambda^2 a_{l',l}^2 u(1-\frac{u}{\beta})Q} \\ = \frac{1}{\lambda|a_{l',l}|} \sqrt{\frac{\beta}{Q}} \left\{ \text{DF} \left[ \frac{1}{2\lambda|a_{l',l}|} \sqrt{\frac{\beta}{Q}} (\lambda^2 a_{l',l}^2 Q - \omega_{l,l'}) \right] \right. \\ \left. + e^{\beta\omega_{l,l'}} \text{DF} \left[ \frac{1}{2\lambda|a_{l',l}|} \sqrt{\frac{\beta}{Q}} (\lambda^2 a_{l',l}^2 Q + \omega_{l,l'}) \right] \right\}, \quad (24)$$

where $DF(x) := e^{-x^2} \int_0^x du 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 [25, 26]. Thus, in the strong coupling regime one can expect to have $\lambda^2 Q \gg \max_l |h_l|$, implying

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

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

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

## 4 Example: spin-boson model

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

$$\mathcal{H}_{SB} = \frac{\epsilon}{2}\sigma_z + \frac{\Delta}{2}\sigma_x + \lambda\sigma_z B + H_B, \quad (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}(|e\rangle\langle e| - |g\rangle\langle g|), \quad (28)$$

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

$$|e\rangle := \frac{(\omega_S + \epsilon)|+\rangle + \Delta|-\rangle}{\sqrt{2\omega_S(\omega_S + \epsilon)}}, \\ |g\rangle := \frac{-\Delta|+\rangle + (\omega_S + \epsilon)|-\rangle}{\sqrt{2\omega_S(\omega_S + \epsilon)}}, \quad (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'} = a_{+,+} = -a_{-,+} = 2$, $h_{l=\pm} = \pm\epsilon/2$, $h_{l,l'} = h_{+,+} = h_{-,+} = \Delta/2$, and $\omega_{l,l'} = h_l - h_{l'} = \omega_{+-} = -\omega_{+-} = \epsilon$. Additionally, since $\sigma_z^2 = 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, \quad (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 7.5, which are given by $\mathbf{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 (20), we obtain

$$\rho_S^{ss} = \frac{1}{e^{-\beta\epsilon/2} + e^{\beta\epsilon/2}} \left[ e^{-\beta\epsilon/2}|+\rangle\langle +| + e^{\beta\epsilon/2}|-\rangle\langle -| \right. \\ \left. - \frac{\Delta}{2}e^{-\beta\epsilon/2}f_{+-,-(\beta)}|+\rangle\langle -| - \frac{\Delta}{2}e^{\beta\epsilon/2}f_{-+,(\beta)}|-\rangle\langle +| \right], \quad (31)$$

with

$$f_{+-,-(\beta)} = \int_0^\beta du e^{\epsilon u} e^{-4\lambda^2 \int_0^\infty d\omega \frac{J(\omega)}{\omega^2} (e^{iu\omega} - 1) \left(1 - \frac{1-e^{-\omega u}}{1-e^{-\omega\beta}}\right)}, \\ f_{-+,(\beta)} = \int_0^\beta du e^{-\epsilon u} e^{-4\lambda^2 \int_0^\infty d\omega \frac{J(\omega)}{\omega^2} (e^{iu\omega} - 1) \left(1 - \frac{1-e^{-\omega u}}{1-e^{-\omega\beta}}\right)}. \quad (32)$$

## 4.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 (24) of Section 3.3 applies, and we have,

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

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} + O[(\lambda^2 Q\beta)^{-3}]. \quad (34)$$
## 5 Comparison with previous results

In a recent paper [44], Trushechkin obtains corrections to the ultrastrong coupling limit by deriving and solving a strong-decoherence regime master equation (generalization of the ultra-strong coupling regime). Additionally, a high temperature expansion of the mean force Gibbs state was recently derived in [53]. In this section, we compare our result with the aforementioned ones. We start by briefly introducing them.

Trushechkin. In [44], the author derives the first order correction to the steady state in the ultra strong coupling limit, which actually coincides with the general expressions (20) and (21) (see details in Section 7.6). As a technical note, the results in [44] are actually valid when the bath spectral density satisfies $\lim_{\omega \to +\infty} \frac{J(\omega)}{\omega^2} > 0$, and the equivalence of the first order correction in [44] with (20) and (21) is guaranteed only within this condition on the bath spectral density (satisfied by usual spectral densities).

As a sanity check for the numerical simulations below, we compare, for the spin-boson model, our approximate expression (33) of $f_{+,-}(\beta)$ with the high temperature approximation in [44] based on 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}, \quad (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 [44], but one can make it appear by multiplying the bath spectral density by $\lambda^2$. Then, the reorganization energy becomes $\lambda^2 Q$. Additionally, the steady state populations obtained in [44] correspond to the one obtained for infinite coupling limit [25], as we also derived in (15) and (20).

Timofeev & Trushechkin. The derivation in [53], a generalization to arbitrary system-bath coupling of [54], consists in expressing the mean force Gibbs state through approximating the Hamiltonian of mean force. In other words, $\rho_S^{ss} = Z_{SB}^{-1} \text{Tr}_B[e^{-\beta H_{SB}}]$ is expressed in the form $Z_{MF}^{-1} e^{-\beta H_{MF}}$, and an approximated expression of $H_{MF}$, the Hamiltonian of mean force, is provided to second order in $\beta$ in the following form [53]

$$H_{MF} = H_S^{\text{diag}} - \Lambda A^2 + \sum_{n \neq m} J_{nm} e^{-\beta \Lambda (a_n - a_m)^2 / 6} |n\rangle\langle m| + \mathcal{O}(\beta^4), \quad (36)$$

where $\Lambda$ is the re-organization energy, $A = \sum_n a_n |n\rangle\langle n|$ is the system operator coupling with the bath, as previously,

$H_S^{\text{diag}} = \sum_n \langle n | H_S | n \rangle |n\rangle\langle n|$, and $J_{nm} := \langle n | H_S | m \rangle$. Reexpressed using the notations used throughout our paper, we have

$$H_{MF} = \sum_{l} h_{l} |a_{l}\rangle\langle a_{l}| - \lambda^{2}QA^{2} \\ + \sum_{l \neq l'} h_{l,l'} e^{-\beta\lambda^{2}Q(a_{l}-a_{l'})^{2}/6}|a_{l}\rangle\langle a_{l'}|. \quad (37)$$

The term $-\lambda^2 QA^2$ corresponds to the renormalization of the system's energies, which has to be performed initially as we did above, or finally, as shown in Section 7. The partition function is given by $Z_{MF} = \text{Tr}_B[e^{-\beta H_{MF}}]$.

Applying the above expression to the spin-boson model, we obtain

$$H_{MF} = \frac{\epsilon}{2}\sigma_z - \lambda^2 Q + \frac{\Delta}{2}e^{-2\beta\lambda^2 Q/3}\sigma_x \\ = \frac{\epsilon}{2}\sigma_z + \frac{\Delta}{2}e^{-2\beta\lambda^2 Q/3}\sigma_x, \quad (38)$$

where we dropped the renormalization term in the second line since in this situation it only corresponds to redefining the origin of the energies. Thus, the expression of the mean force Gibbs state derived in [53] and applied to the present spin-boson model is

$$\rho_{MF} = Z_{MF}^{-1} [e^{-\omega'_{S}\beta/2}|e'\rangle\langle e'| + e^{\omega'_{S}\beta/2}|g'\rangle\langle g'|], \quad (39)$$

with $\omega'_{S} := \sqrt{\epsilon^2 + \Delta'^2}$, $\Delta' := e^{-2\beta\lambda^2 Q/3}\Delta$, $|e'\rangle := (\omega'_{S} + \epsilon)|+\rangle + |\Delta' - \rangle$, $|g'\rangle := \frac{-\Delta' |+\rangle + (\omega'_{S} + \epsilon)|-\rangle}{\sqrt{2\omega'_{S}(\omega'_{S} + \epsilon)}}$, and $Z_{MF} = e^{-\omega'_{S}\beta/2} + e^{\omega'_{S}\beta/2}$. This expression has the merit of tending to the right limit when $\lambda \to 0$, namely $\rho_{MF} \to Z^{-1}e^{-\beta H_S}$, the usual thermal equilibrium state. This is not the case of the expression (31), which is expected since it is meant to be valid in the opposite regime, when $\lambda^2 Q\beta \gg 1$. In the eigenbasis of $\sigma_z$, the expression of $\rho_{MF}$ becomes

$$\rho_{MF} = \frac{1}{2} \left( 1 - \frac{\epsilon}{\omega'_{S}} \tanh \frac{\omega'_{S} \beta}{2} \right) |+\rangle\langle +| \\ + \frac{1}{2} \left( 1 + \frac{\epsilon}{\omega'_{S}} \tanh \frac{\omega'_{S} \beta}{2} \right) |-\rangle\langle -| \\ - \frac{\Delta'}{2\omega'_{S}} \tanh \frac{\omega'_{S} \beta}{2} (|+\rangle\langle -| + |-\rangle\langle +|). \quad (40)$$

A brief analytical comparison between (40) and (31) is provided in Section 7.4.

In Fig. 1, we have plotted in semi-log scale the coherences of the mean force Gibbs state (a) $c_S^{ss} := \langle +|\rho_S^{ss}|-\rangle$ in the eigenbasis of A, and (b) $c_{e,g}^{ss} := \langle e|\rho_S^{ss}|g\rangle$ in the eigenbasis of $H_S$, in function of the coupling strength $\lambda^2 Q$ in unit of $\epsilon$, for $\Delta = 1.5\epsilon$, $\omega_c = 2\epsilon$ and $\epsilon\beta = 0.1$. The orange
<figure><img src="image_4.png" /></figure>

**Figure 1:** Plots in the semi-log scale of the mean force Gibbs state coherences (a) $c^{ss} := \langle +|\rho_S^{ss}|-\rangle$ in the eigenbasis of A and (b) $c_{e,g}^{ss} := \langle e|\rho_S^{ss}|g\rangle$ in the eigenbasis of $H_S$, in function of the coupling strength $\lambda^2 Q$ in unit of $\epsilon$, for $\Delta = 1.5\epsilon$, $\omega_c = 2\epsilon$ and $\epsilon\beta = 0.1$. The orange thick solid curve corresponds to (31) using the high temperature approximation (33) of $f_{+,-}(\beta)$, while the gray thin dashed line corresponds to the further simplified expression (34) of $f_{+,-}(\beta)$. The purple thick dashed line is the mean force Gibbs state coherences given by the high temperature expression derived in [44]. The black thin line is the mean force Gibbs state coherences given by (40).

thick solid curve corresponds to (31) using the high temperature approximation (33) of $f_{+,-}(\beta)$, while the gray thin dashed line corresponds to the further simplified expression (34) of $f_{+,-}(\beta)$. The purple thick dashed line is the mean force Gibbs state coherences given by the high temperature expression derived in [44]. The black thin line is the mean force Gibbs state coherences given by (40).

As expected, on can see an excellent agreement between (33) (orange curve) and the expression given in [44] (purple dashed line). The very slight discrepancy appearing at intermediate coupling strength comes from a slight difference in the derivation of high temperature approximation between (22)-(23) and eqs. (56)-(57) of [44] (high temperature approximation before time integration of the bath correlation function). However, there is a significant discrepancy with (40) (black line) at intermediary and strong coupling strength. Beyond that, it appears from (a) that for increasing coupling strength, the mean

<figure><img src="image_5.png" /></figure>

**Figure 2:** Plots in the semi-log scale of the mean force Gibbs state coherences (a) $c^{ss} := \langle +|\rho_S^{ss}|-\rangle$ in the eigenbasis of A, and (b) $c_{e,g}^{ss} := \langle e|\rho_S^{ss}|g\rangle$ in the eigenbasis of $H_S$, in function of the inverse temperature $\epsilon\beta$, for $\Delta = 1.5\epsilon$, $\omega_c = 1\epsilon$ and $\lambda^2 Q = 10\epsilon$. The orange thick solid curve corresponds to (31) using the high temperature approximation (33) of $f_{+,-}(\beta)$, while the gray thin dashed line (almost indistinguishable from the orange curve) corresponds to the further simplified expression (34) of $f_{+,-}(\beta)$. The purple thick dashed line is the mean force Gibbs state coherences given by the high temperature expression derived in [44]. The black thin line is the mean force Gibbs state coherences given by (40).

force Gibbs state tends to a diagonal state in the eigenbasis of A, as expected from the ultrastrong coupling limit (15) and [25, 46, 47]. Conversely, for decreasing coupling strength, one can see in (b) the emergence of a progressive transition to a diagonal state in the eigenbasis of $H_S$, as expected from the weak coupling limit.

In Fig. 2, we have plotted in the semi-log scale of the mean force Gibbs state coherences (a) $c^{ss} := \langle +|\rho_S^{ss}|-\rangle$ in the eigenbasis of A, and (b) $c_{e,g}^{ss} := \langle e|\rho_S^{ss}|g\rangle$ in the eigenbasis of $H_S$, in function of the inverse temperature $\epsilon\beta$, for $\Delta = 1.5\epsilon$, $\omega_c = 1\epsilon$ and $\lambda^2 Q = 10\epsilon$. The colour convention is the same as in the previous figure Fig. 1. Again, as expected, we observe a very good agreement between (33) (orange curve) and the expression given in [44] (purple dashed line). We also observe a significant discrepancy with (40) (black line) out side the high temperature regime.
<figure><img src="image_6.png" /><figcaption>**Figure 3:** Plots of the mean force Gibbs state excited population (a) $p_+^{ss} := \langle +|\rho_S^{ss}|+\rangle$ in the eigenbasis of A, and (b) $p_e^{ss} := \langle e|\rho_S^{ss}|e\rangle$ in the eigenbasis of $H_S$, in function of the coupling strength $\lambda^2 Q$ in unit of $\epsilon$, for $\Delta = 1.5\epsilon$ and $\epsilon\beta = 0.5$. The orange thick solid curve corresponds to the mean force Gibbs state excited population obtained in (33), [44], and [25, 46, 47] (all coinciding). The black thin curve corresponds to the excited population given by (40). The blue dotted curve corresponds to the thermal population in the vanishing coupling limit, $p_e^{th} = e^{-\omega_s\beta}/(1+e^{-\omega_s\beta})$.</figcaption></figure>

In Figs. 3 and 4, we have plotted the mean force Gibbs state excited population (a) $p_+^{ss} := \langle +|\rho_S^{ss}|+\rangle$ in the eigenbasis of A, and (b) $p_e^{ss} := \langle e|\rho_S^{ss}|e\rangle$ in the eigenbasis of $H_S$. As for the coherence, Fig. 3 is in function of the coupling strength $\lambda^2 Q$ in unit of $\epsilon$, for $\Delta = 1.5\epsilon$ and $\epsilon\beta = 0.5$, while in Fig. 4, the plots are in function of the inverse temperature $\epsilon\beta$ in a semi-log scale, for $\Delta = 1.5\epsilon$ and $\lambda^2 Q = 5\epsilon$. For both figures, the orange thick solid curve corresponds to the excited populations of the mean force Gibbs state obtained in (33), [44], and [25, 46, 47] (all coinciding), while the black thin curve corresponds to the excited population given by (40). Finally, the blue dotted curve corresponds to the thermal population in the vanishing coupling limit, $p_e^{th} = e^{-\omega_s\beta}/(1+e^{-\omega_s\beta})$.

We can draw conclusions similar to Figs. 1 and 2, namely that both predictions coincide very well at high temperature as well as for ultra strong coupling. However, some discrepancies emerge for intermediate and large coupling strength as well as for large $\beta$.

<figure><img src="image_7.png" /><figcaption>**Figure 4:** Plots of the mean force Gibbs state excited population (a) $p_+^{ss} := \langle +|\rho_S^{ss}|+\rangle$ in the eigenbasis of A, and (b) $p_e^{ss} := \langle e|\rho_S^{ss}|e\rangle$ in the eigenbasis of $H_S$, in function of the inverse temperature $\epsilon\beta$ in a semi-log scale, for $\Delta = 1.5\epsilon$ and $\lambda^2 Q = 10\epsilon$. The orange thick solid curve corresponds to the mean force Gibbs state excited population obtained in (33), [44], and [25, 46, 47] (all coinciding). The black thin curve corresponds to the excited population given by (40). The blue dotted curve corresponds to the thermal population in the vanishing coupling limit, $p_e^{th} = e^{-\omega_s\beta}/(1+e^{-\omega_s\beta})$.</figcaption></figure>

# 6 Conclusion

We have obtained a perturbative expansion (13) and (14) of the mean force Gibbs state in the ultra strong coupling regime, using a technique inspired from the displaced oscillator picture [49,50] and polaron transformation. The zero order term confirms the result of [25, 44]. Analytical expression of the first order term is derived, (20) and (21), and is found to coincide with the first order correction derived in [44] (within the mild validity condition of [44], see Section 7.6 for more detail). This increases the range of the known situations in which the steady state of a system interacting with a thermal bath does converge to the mean force Gibbs state. This convergence was recently shown for weak coupling up to the first order [24, 40, 45], for the ultra strong coupling regime up to the zeroth order [25], and is now confirmed in the ultra-strong coupling regime up to the first order, thanks to the complementary results in [44].
In the high temperature regime, a simple and explicit general expression is derived (24) for the first order corrections. This expression is numerically compared for the spin-boson model with the result of [53] obtained by high temperature expansion of the mean force Gibbs state. For sanity check, a comparison of (24) with the high temperature approximation of [44] is also made, and the very good agreement in all regimes of parameters is confirmed. Very good agreement is also observed at ultra strong coupling as well as high temperatures between our result and the one from [53]. This strengthens the validity of all three approaches. Thanks to them, we are able to draft the transition from the ultra-strong coupling regime to the weak coupling regime, Fig. 1b.

However, going further, significant discrepancies emerge between our result and [53] at moderate coupling strength and moderate temperatures. For those values of parameters, further analysis with higher orders are needed in order to confirm which result is the more accurate, as well as to reconstruct accurately the full transition weak-to-ultra-strong coupling and provide a good approximation of the steady state in all regimes. This does not seem out-of-reach according to the preliminary results on higher order terms obtained in Section 7.3. Additionally, benchmarking the results with other techniques like reaction coordinate or pseudo-modes would be interesting and instructive.

# 7 Appendix

## 7.1 Details of the computation of

$$ \mathrm{Tr}_B[e^{-\beta H_{B,l}} e^{uH_{B,l}} e^{-uH_{B,l'}}] $$

In order to compute $\mathrm{Tr}_B[e^{-\beta H_{B,l}} e^{uH_{B,l}} e^{-uH_{B,l'}}]$, one can first rearrange it in the following way,

$$ \begin{aligned} \mathrm{Tr}_B[e^{-\beta H_{B,l}} e^{uH_{B,l}} e^{-uH_{B,l'}}] &= \mathrm{Tr}_B[e^{-(\beta-u)H_{B,l}} e^{-uH_{B,l'}}] \\ &= \mathrm{Tr}_B\left[\otimes_k \left(\mathcal{D}_{k,l} e^{-(\beta-u)b_k^\dagger b_k} \mathcal{D}_{k,l}^\dagger\right) e^{-uH_{B,l'}}\right] \\ &= \mathrm{Tr}_B\left[e^{-(\beta-u)H_B} e^{-u\sum_k \omega_k \mathcal{D}_{k,l}^\dagger \mathcal{D}_{k,l'} b_k^\dagger b_k \mathcal{D}_{k,l'}^\dagger \mathcal{D}_{k,l}}\right] \\ &= \mathrm{Tr}_B\left[e^{-(\beta-u)H_B} e^{-u(H_B+\lambda a_{l',l} B+\lambda^2 a_{l',l}^2 Q)}\right], \end{aligned} \quad (41) $$

where $a_{l',l} := a_{l'} - a_l$. “Taking out” of the second exponential the Hamiltonian $H_B$, we obtain,

$$ \begin{aligned} \mathrm{Tr}_B \left[ e^{-(\beta-u)H_B} e^{-u(H_B+\lambda a_{l',l} B+\lambda^2 a_{l',l}^2 Q)} \right] \\ = e^{-u(\lambda^2 a_{l',l}^2 Q) \mathrm{Tr}_B \left[ e^{-\beta H_B} e^{-\lambda a_{l',l} T \int_0^u ds \tilde{B}(s)} \right]}, \end{aligned} \quad (42) $$

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)$. The time ordered operator $e^{-\lambda a_{l',l} T \int_0^u ds \tilde{B}(s)}$ can be split in two

as follows,

$$ e^{-\lambda a_{l',l} T \int_0^u ds \tilde{B}(s)} = e^{-\lambda^2 a_{l',l}^2 f(u)} e^{\lambda a_{l',l} B(-u)} e^{-\lambda a_{l',l} \tilde{B}^\dagger(u)}, \quad (43) $$

with $f(u) := -uQ + \sum_k \frac{g_k^2}{\omega_k^2}(e^{u\omega_k} - 1)$, $B(-u) := \sum_k \frac{g_k}{\omega_k}(e^{-u\omega_k} - 1)b_k$, and $\tilde{B}^\dagger(u) := \sum_k \frac{g_k}{\omega_k}(e^{u\omega_k} - 1)b_k^\dagger$. The above decomposition can be demonstrated as follows. We denote the left-hand side of (43) as $L(u) := e^{-\lambda a_{l',l} T \int_0^u ds \tilde{B}(s)}$ and the right-hand side by $R(u) := e^{-\lambda^2 a_{l',l}^2 f(u)} e^{\lambda a_{l',l} B(-u)} e^{-\lambda a_{l',l} \tilde{B}^\dagger(u)}$. By definition of the time ordering operator, or alternatively by taking the time derivative with respect to $u$, the left-hand side of (43) satisfies the differential equation

$$ \frac{\partial}{\partial u} L(u) = -\lambda a_{l',l} \tilde{B}(u) L(u). \quad (44) $$

Now, by taking the time derivative with respect to $u$ of the right-hand side of (43), one obtains the following differential equation

$$ \begin{aligned} \frac{\partial}{\partial u} R(u) = & \left[ -\lambda^2 a_{l',l}^2 \frac{\partial}{\partial u} f(u) + \lambda a_{l',l} \frac{\partial}{\partial u} B(-u) \right] R(u) \\ & - \lambda a_{l',l} R(u) \frac{\partial}{\partial u} \tilde{B}^\dagger(u). \end{aligned} \quad (45) $$

Using the expressions of $f(u)$, $B(u)$, and $\tilde{B}^\dagger(u)$, we obtain,

$$ \begin{aligned} \frac{\partial}{\partial u} f(u) &= -Q + \sum_k \frac{g_k^2}{\omega_k} e^{u\omega_k}, \\ \frac{\partial}{\partial u} B(-u) &= -\sum_k g_k e^{-u\omega_k} b_k, \\ \frac{\partial}{\partial u} \tilde{B}^\dagger(u) &= \sum_k g_k e^{u\omega_k} b_k^\dagger. \end{aligned} \quad (46) $$

Then, with the help of the identity

$$ R(u) \frac{\partial}{\partial u} \tilde{B}^{\dagger}(u) = \left[ \frac{\partial}{\partial u} \tilde{B}^{\dagger}(u) + \lambda a_{l',l} \sum_k \frac{g_k^2}{\omega_k} (1 - e^{u\omega_k}) \right] R(u), \quad (47) $$

and the discrete version of the re-organization energy $Q = \sum_k \frac{g_k^2}{\omega_k} \equiv \int_0^\infty d\omega \frac{J(\omega)}{\omega}$, we arrive at

$$ \frac{\partial}{\partial u} R(u) = -\lambda a_{l',l} \left[ \sum_k g_k \left( e^{u\omega_k} b_k^\dagger + e^{-u\omega_k} b_k \right) \right] R(u), \quad (48) $$

which is exactly the same differential equation as (44). Since the initial conditions are the same, namely $R(0) = L(0) = I$, we conclude that $R(u) = L(u)$ for all times $u$.
A simple way to conclude the computation of the trace (42) is using the normal order representation of the operators. Using the following formula [58, p. 116]

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

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

$$\mathrm{Tr}_B \left[ e^{-(\beta-u)H_B} e^{-u(H_B+\lambda a_{\prime l} B+\lambda^2 a_{\prime l}^2 Q)} \right] = e^{-u(\lambda^2 a_{\prime l}^2 Q)} e^{-\lambda^2 a_{\prime l}^2 f(u) \mathrm{Tr}_B \left[ e^{-\lambda a_{\prime l} B \dagger(u)} \left( \otimes_k \sum_{n=0}^{\infty} \frac{(e^{-\beta \omega_k} - 1)^n}{n!} (b_k^{\dagger})^n b_k^n \right) e^{\lambda a_{\prime l} B(-u)} \right]} \quad (50)$$

$$= e^{-\lambda^2 a_{\prime l}^2 \sum_k \frac{g_k^2}{\omega_k^2} (e^{u\omega_k}-1)} \prod_k \left( \int_C \frac{d^2\alpha_k}{\pi} e^{-\lambda a_{\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_{\prime l} \frac{g_k}{\omega_k} (e^{-u\omega_k}-1)\alpha_k} \right) \quad (51)$$

$$= e^{-\lambda^2 a_{\prime l}^2 \sum_k \frac{g_k^2}{\omega_k^2} (e^{u\omega_k}-1)} \prod_k \left( \int_C \frac{d^2\alpha_k}{\pi} e^{-\lambda a_{\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_{\prime l} \frac{g_k}{\omega_k} (e^{-u\omega_k}-1)\alpha_k} \right)$$

$$= e^{-\lambda^2 a_{\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{dx_k dy_k}{\pi} e^{-\lambda a_{\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_{\prime l} \frac{g_k}{\omega_k} (e^{-u\omega_k} - e^{u\omega_k})x_k + i\lambda a_{\prime l} \frac{g_k}{\omega_k} (e^{-u\omega_k} + e^{u\omega_k}-2)y_k} \right)$$

$$= e^{-\lambda^2 a_{\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{dx_k dy_k}{\pi} e^{(e^{-\beta\omega_k}-1)(x_k^2+y_k^2)+\lambda a_{\prime l} \frac{g_k}{\omega_k} (e^{-u\omega_k}-e^{u\omega_k})x_k+i\lambda a_{\prime l} \frac{g_k}{\omega_k} (e^{-u\omega_k}+e^{u\omega_k}-2)y_k} \right)$$

$$= e^{-\lambda^2 a_{\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_{\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}}}$$

$$= Z_B e^{-\lambda^2 a_{\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)}$$

$$= Z_B e^{-\lambda^2 a_{\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)}, \quad (52)$$

where $Z_B := \Pi_k(\frac{1}{1-e^{-\omega_k\beta}}) = \mathrm{Tr}_B(e^{-\beta H_B})$. In line (51), 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 $O_k$ acting in that Hilbert space can be computed as $\mathrm{Tr}_k(O_k) = \int_{\mathbb{R}^2} \frac{d\omega_k^2}{\pi} \langle\alpha_k|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$.

## 7.2 Important property of $f_{l,l'}(\beta)$

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

$$\begin{aligned} f_{l,l'}(\beta) &= \int_0^\beta du e^{u\omega_{l,l'}} e^{-\lambda^2 a_{\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}}} } \\ &= e^{\beta\omega_{l,l'}} \int_0^\beta du e^{-v\omega_{l,l'}} e^{-\lambda^2 a_{\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}}} } . \end{aligned} \quad (53)$$

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

## 7.3 Higher order terms

As explained in Section 3, we consider only the first two terms of the infinite sum in Eq. (14). In Section 3.1, we gave an intuitive argument to justify that the second term as well as all following terms in Eq. (14) converge to zero when the coupling strength increases. We now come back to this point and provide a more rigorous argument.
For simplicity, we focus here on the spin-boson model used in Section 4, while this could be easily extended to arbitrary systems. Combining Eqs. (13) and (14), we start re-writing the exact expression of $\rho_S^{ss}$ as

$$ \rho_S^{ss} = (Z_S^{ss})^{-1} \left[ e^{-\beta h_+} \sum_{n=0}^{\infty} T_{+,n} + e^{-\beta h_-} \sum_{n=0}^{\infty} T_{-,n} \right], \quad (54) $$

with $h_{\pm} = \langle \pm | H_S | \pm \rangle = \pm \epsilon/2$, and where

$$ T_{\pm,n} := (-1)^n \int_0^\beta du_1 \int_0^{u_1} du_2 \dots \int_0^{u_{n-1}} du_n |\pm\rangle\langle\pm| \mathrm{Tr}_B [Z_B^{-1} e^{-\beta H_{B,\pm}} \tilde{H}_S^{\mathrm{coh}}(u_1) \tilde{H}_S^{\mathrm{coh}}(u_2) \dots \tilde{H}_S^{\mathrm{coh}}(u_n)], \quad (55) $$

with $H_{B,\pm} := H_B \pm \lambda B + \lambda^2 Q$, $T_{\pm,0} = |\pm\rangle\langle\pm|$ and

$$ Z_S^{ss} := \frac{Z_{SB}}{Z_B} = e^{-\beta\epsilon/2} \sum_{n=0}^{\infty} \langle +|T_{+,n}|+\rangle + e^{\beta\epsilon/2} \sum_{n=0}^{\infty} \langle -|T_{-,n}|-\rangle. \quad (56) $$

We then work on the integrand

$$ t_{\pm,n} := |\pm\rangle\langle\pm| \mathrm{Tr}_B [e^{-\beta H_{B,\pm}} \tilde{H}_S^{\mathrm{coh}}(u_1) \tilde{H}_S^{\mathrm{coh}}(u_2) \dots \tilde{H}_S^{\mathrm{coh}}(u_n)]. \quad (57) $$

Starting with $t_{+,n}$, we have

$$
\begin{aligned}
t_{+,n} &= |+\rangle\langle+| \mathrm{Tr}_B [Z_B^{-1} e^{-\beta H_{B,+}} \tilde{H}_S^{\mathrm{coh}}(u_1) \tilde{H}_S^{\mathrm{coh}}(u_2) \dots \tilde{H}_S^{\mathrm{coh}}(u_n)] \\
&= Z_B^{-1} |+\rangle\langle+| \mathrm{Tr}_B \left[ e^{-\beta H_{B,+}} (h_{+,+} e^{u_1 \omega_{+,+}} |+\rangle\langle-| e^{u_1 H_{B,+}} e^{-u_1 H_{B,-}} + h_{-,+} e^{u_1 \omega_{-,+}} |-\rangle\langle+| e^{u_1 H_{B,-}} e^{-u_1 H_{B,+}}) \right. \\
&\qquad \left. \times \tilde{H}_S^{\mathrm{coh}}(u_2) \dots \tilde{H}_S^{\mathrm{coh}}(u_n) \right] \\
&= Z_B^{-1} h_{+,+} e^{u_1 \omega_{+,+}} |-\rangle\langle+| \mathrm{Tr}_B \left[ e^{-(\beta-u_1)H_{B,+}} e^{-u_1 H_{B,-}} \tilde{H}_S^{\mathrm{coh}}(u_2) \dots \tilde{H}_S^{\mathrm{coh}}(u_n) \right] \\
&= Z_B^{-1} h_{+,+} e^{u_1 \omega_{+,+}} |-\rangle\langle+| \mathrm{Tr}_B \left[ e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}} e^{-u_2 H_{B,+}} \tilde{H}_S^{\mathrm{coh}}(u_3) \dots \tilde{H}_S^{\mathrm{coh}}(u_n) \right] \\
&= Z_B^{-1} h_{+,+} h_{-,+} e^{u_1 \omega_{+,+}} e^{u_2 \omega_{-,+}} |-\rangle\langle+| \mathrm{Tr}_B \left[ e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}} e^{-u_2 H_{B,+}} \tilde{H}_S^{\mathrm{coh}}(u_3) \dots \tilde{H}_S^{\mathrm{coh}}(u_n) \right] \\
&= Z_B^{-1} \left(\frac{\Delta}{2}\right)^n e^{\epsilon(u_1-u_2+\dots+(-1)^n u_n)} \\
&\qquad \times \mathrm{Tr}_B \left[ e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}} \dots e^{-(u_{n-1}-u_n)H_{B,(n-1)n-1}} e^{-u_n H_{B,(n-1)n}} \right] |-\rangle\langle(-1)^n|, \quad (58)
\end{aligned}
$$

since $h_{+,+} = h_{-,+} = \Delta/2$ and $\omega_{+,+} = -\omega_{-,+} = \epsilon$. Reproducing the decomposition used in Eq.(42) and (43), we obtain

$$ e^{-\alpha_i H_{B,\pm}} = e^{-\lambda^2 \bar{f}(\alpha_i)} e^{-\alpha_i H_B} e^{\pm \lambda B(-\alpha_i)} e^{\mp \lambda B^\dagger(\alpha_i)} \quad (59) $$

with $\bar{f}(\alpha_i) := \sum_k \frac{g_k^2}{\omega_k^2} (e^{\alpha_i \omega_k} - 1)$ and $\alpha_i = u_i - u_{i+1}$. The aim is now to re-write

$$ e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}} \dots e^{-(u_{n-1}-u_n)H_{B,(n-1)n-1}} e^{-u_n H_{B,(n-1)n}} \quad (60) $$

in normal order so that we can compute the trace. To have a better understanding of how we will proceed, we first consider the simpler situation where $n=1$. We have, using (59) twice,

$$ e^{-(\beta-u_1)H_{B,+}} e^{-u_1 H_{B,-}} = e^{-\lambda^2 \bar{f}(\beta-u_1)} e^{-\lambda^2 \bar{f}(u_1)} e^{-(\beta-u_1)H_B} e^{\lambda B(-\beta+u_1)} e^{-\lambda B^\dagger(\beta-u_1)} e^{-u_1 H_B} e^{-\lambda B(-u_1)} e^{\lambda B^\dagger(u_1)}. \quad (61) $$
$$ 
\begin{aligned}
e^{\chi_k b_k} e^{-\alpha_i \omega_k b_k^\dagger b_k} &= e^{-\alpha_i \omega_k b_k^\dagger b_k} e^{\chi_k} e^{-\alpha_i \omega_k b_k} \\
e^{\chi_k b_k^\dagger} e^{-\alpha_i \omega_k b_k^\dagger b_k} &= e^{-\alpha_i \omega_k b_k^\dagger b_k} e^{\chi_k} e^{\alpha_i \omega_k b_k^\dagger}
\end{aligned} 
\qquad (62) 
$$

for arbitrary coefficient $\chi_k$, we have

$$ e^{-(\beta-u_1)H_{B,+}} e^{-u_1 H_{B,-}} = e^{-\lambda^2 \bar{f}(\beta-u_1)} e^{-\lambda^2 \bar{f}(u_1)} e^{-\beta H_B} e^{\lambda \sum_k \frac{g_k}{\omega_k} e^{-u_1 \omega_k} (e^{-(\beta-u_1)\omega_k}-1)b_k} e^{-\lambda \sum_k \frac{g_k}{\omega_k} e^{u_1 \omega_k} (e^{\beta-u_1}\omega_k-1)b_k^\dagger} e^{-\lambda \mathcal{B}(-u_1)} e^{\lambda \mathcal{B}^\dagger(u_1)}. \qquad (63) $$

We finally use the identity

$$ e^{\chi_k b_k^\dagger} e^{\xi_k b_k} = e^{-\chi_k \xi_k} e^{\xi_k b_k} e^{\chi_k b_k^\dagger}, \qquad (64) $$

to obtain,

$$ 
\begin{aligned}
e^{-(\beta-u_1)H_{B,+}} e^{-u_1 H_{B,-}} &= e^{-\lambda^2 \bar{f}(\beta-u_1)} e^{-\lambda^2 \bar{f}(u_1)} e^{-\lambda^2 \sum_k \frac{g_k^2}{\omega_k^2} e^{u_1 \omega_k} (e^{-(\beta-u_1)\omega_k}-1)(e^{-u_1 \omega_k}-1)} e^{-\beta H_B} \\
&\quad \times e^{\lambda \sum_k \frac{g_k}{\omega_k} [e^{-u_1 \omega_k} (e^{-(\beta-u_1)\omega_k}-1)-(e^{-u_1 \omega_k}-1)]b_k} e^{-\lambda \sum_k \frac{g_k}{\omega_k} [e^{u_1 \omega_k} (e^{\beta-u_1}\omega_k-1)-(e^{u_1 \omega_k}-1)]b_k^\dagger}.
\end{aligned} 
\qquad (65) $$

The above expression is not yet in normal order, but once injected in the trace, one can use the permutation property as in the previous section to obtain an expression in normal order, exactly of the same form as in (50), the only difference being the coefficients in front of the annihilation and creation operators $b_k$ and $b_k^\dagger$. We then proceed to the computation of the trace exactly in the same way as in the previous section. We obtain something of the form,

$$ \mathrm{Tr}_B \left[ e^{-(\beta-u_1)H_{B,+}} e^{-u_1 H_{B,-}} \right] = Z_B e^{-\lambda^2 \Psi_{+,-}(\beta, u_1)}, \qquad (66) $$

where $\Psi_{+,-}(\beta, u_1)$ is a real function.

If now we consider the term $n = 2$, we can use the above form for $e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}}$, and then split $e^{-u_2 H_{B,+}}$ using (59). We obtain an expression of the same form as (61), and thus one can repeat the above steps to reach an expression of the later form (66), namely,

$$ \mathrm{Tr}_B \left[ e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}} e^{-u_2 H_{B,+}} \right] = Z_B e^{-\lambda^2 \Psi_{+,-,+}(\beta, u_1, u_2)}. \qquad (67) $$

We can repeat the same procedure for any $n$, obtaining

$$ \mathrm{Tr}_B \left[ e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}} \dots e^{-(u_{n-1}-u_n)H_{B,(-1)^{n-1}}} e^{-u_n H_{B,(-1)^n}} \right] = Z_B e^{-\lambda^2 \Psi_{+,-,+,...,(-1)^n}(\beta, u_1, u_2, \dots, u_n)}, \qquad (68) $$

where $\Psi_{+,-,+,...,(-1)^n}(\beta, u_1, u_2, \dots, u_n)$ is a real function.

Having this result in mind, one can obtain a more direct form of $\Psi_{+,-,+,...,(-1)^n}(\beta, u_1, u_2, \dots, u_n)$ simply by computing

$$ \left. \frac{\partial^2}{\partial \lambda^2} \mathrm{Tr}_B \left[ e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}} \dots e^{-(u_{n-1}-u_n)H_{B,(-1)^{n-1}}} e^{-u_n H_{B,(-1)^n}} \right] \right|_{\lambda=0} . \qquad (69) $$

After computing this second derivative and evaluating it in $\lambda = 0$, followed by some simple algebraic manipulations, we arrived at

$$ 
\begin{aligned}
\Psi_{+,-,+,...,(-1)^n}(\beta, u_1, u_2, \dots, u_n) &= -\frac{1}{2} \frac{\partial^2}{\partial \lambda^2} \mathrm{Tr}_B \left[ Z_B^{-1} e^{-(\beta-u_1)H_{B,+}} e^{-(u_1-u_2)H_{B,-}} \dots e^{-(u_{n-1}-u_n)H_{B,(-1)^{n-1}}} e^{-u_n H_{B,(-1)^n}} \right] \Big|_{\lambda=0} \\
&= \sum_{0 \le i \le j \le n} (-1)^{i+j} M_{i,j},
\end{aligned} 
\qquad (70) $$

with

$$ M_{i,j} = \sum_k \frac{g_k^2}{\omega_k^2} (n_{\omega_k} + 1) (e^{-u_i \omega_k} - e^{-u_{i+1} \omega_k}) (e^{-(\beta-u_i) \omega_k} - e^{u_{i+1} \omega_k}), \qquad (71) $$
and

$$M_{i,j} = \sum_{i<j} \frac{g_k^2}{\omega_k^2} \left[ (e^{-u_i \omega_k} - e^{-u_{i+1} \omega_k}) (e^{u_j \omega_k} - e^{u_{j+1} \omega_k}) (n_{\omega_k} + 1) + (e^{u_i \omega_k} - e^{u_{i+1} \omega_k}) (e^{-u_j \omega_k} - e^{-u_{j+1} \omega_k}) n_{\omega_k} \right], \quad (72)$$

with the convention $u_0 := \beta$, and $u_{n+1} = 0$. One can verify that for $n = 1$ we recover the expression of the previous section. We finally obtain for $t_{+,n}$,

$$t_{+,n} = \left(\frac{\Delta}{2}\right)^n e^{\epsilon(u_1-u_2+\dots+(-1)^n u_n)} e^{-\lambda^2 \sum_{0 \le i \le j \le n} (-1)^{j+M_{i,j}} |+\rangle\langle-1}|^n}. \quad (73)$$

Repeating a similar derivation we obtain for $t_{-,n}$,

$$t_{-,n} = \left(\frac{\Delta}{2}\right)^n e^{-\epsilon(u_1-u_2+\dots+(-1)^n u_n)} e^{-\lambda^2 \sum_{0 \le i \le j \le n} (-1)^{j+M_{i,j}} |+\rangle\langle-1}|^n}, \quad (74)$$

(only $\epsilon$ is changed to $-\epsilon$).

Although we now have an explicit expression of $t_{l,n}$, the nested integrals over the $u_i$ appearing in (55) are still challenging in general. Thus, as in the main text, we proceed by looking at the high temperature regime, characterized by $\beta\omega_c \ll 1$, where $\omega_c$ stands for the cut-off frequency of the bath spectral density. Consequently, we can expand all exponential functions appearing in the coefficients $M_{i,j}$ (since all $u_i$ are smaller than $\beta$). We obtain the following simple expressions

$$M_{i,i} = \frac{Q}{\beta}(u_i - u_{i+1})(\beta - u_i + u_{i+1}) \\ M_{i,j} = \sum_{i<j} -2 \frac{Q}{\beta}(u_i - u_{i+1})(u_j - u_{j+1}) \quad (75)$$

Using the above expressions, one can see easily that for $n = 1$, we have $\sum_{0 \le i \le j \le 1} (-1)^{i+j} M_{i,j} = 4 \frac{Q}{\beta} u_1 (\beta - u_1)$, recovering a result from the main text. For $n = 2$, one obtains $\sum_{0 \le i \le j \le 2} (-1)^{i+j} M_{i,j} = 4 \frac{Q}{\beta} (u_1 - u_2)(\beta - u_1 + u_2)$. More generally, one can show by iteration that for arbitrary $n$,

$$\sum_{0 \le i \le j \le n} (-1)^{i+j} M_{i,j} = 4 \frac{Q}{\beta} (u_1 - u_2 + u_3 - \dots - u_n) \\ \times (\beta - u_1 + u_2 - u_3 + \dots - u_n). \quad (76)$$

Finally, we obtain

$$T_{\pm,n} = \frac{(-\Delta)^n}{2^n} \int_0^\beta \int_0^{u_1} \dots \int_0^{u_{n-1}} du_1 du_2 \dots du_n \\ e^{\pm \epsilon(u_1 - u_2 + \dots + u_n)} e^{-4\lambda^2 Q \beta \frac{u_1 - u_2 + \dots + u_n}{\beta} \left(1 - \frac{u_1 - u_2 + \dots + u_n}{\beta}\right)} \\ \times |\pm\rangle\langle\pm(-1)^n| \\ = \frac{(-\beta\Delta)^n}{2^n} \int_0^1 \int_0^{u_1} \dots \int_0^{u_{n-1}} du_1 du_2 \dots du_n \\ e^{\beta(\pm\epsilon - 4\lambda^2 Q)(u_1 - u_2 + \dots + u_n)} e^{4\lambda^2 Q \beta (u_1 - u_2 + \dots + u_n)^2} |\pm\rangle\langle\pm(-1)^n| \quad (77)$$

Additionally, we can show that for $n$ impair, $T_{+,n} = e^{\epsilon\beta} T_{-,n}$, which guarantees the Hermicity of $\rho_S^{ss}$. Note that this relation does not hold for $n$ pair, which does not affect the Hermicity of $\rho_S^{ss}$ since corrections of even order affect only the populations.

Beyond that, since $(u_1 - u_2 + \dots + u_n) \in [0; 1]$, we have the following simple upper bound,

$$|T_{\pm,n}| \le \left(\frac{\beta\Delta}{2}\right)^n \int_0^1 \int_0^{u_1} \dots \int_0^{u_{n-1}} du_1 du_2 \dots du_n \\ \times e^{\pm \epsilon\beta(u_1 - u_2 + \dots + u_n)} \quad (78)$$

implying

$$|T_{+,n}| \le \frac{1}{n!} \left(\frac{\beta\Delta}{2}\right)^n e^{\epsilon\beta} \quad (79)$$

and

$$|T_{-,n}| \le \frac{1}{n!} \left(\frac{\beta\Delta}{2}\right)^n \quad (80)$$

which ensures that the higher order corrections $|T_{\pm,n}|$ vanish quickly as $n$ increases. Additionally, for a fixed $n$, we can show numerically that $|T_{\pm,n}|$ goes to zero as $\lambda$ increases (see plots in Fig. 5), confirming that for large $\lambda$, only the first few orders are enough to obtain a good approximation of $\rho_S^{ss}$.

## 7.4 Further analytical comparison between Eq. (33) and the result from Timofeev & Trushechkin

Here we provide some brief analytical comparison between our results and the one obtained in [53]. As detailed in Section 5, the result from [53] applied to the
<figure><img src="image_8.png" /><figcaption>**Figure 5:** (a) Plots of $T_{+,n}/(-\beta\Delta/2)^n$ in function of $\lambda^2 Q$ in unit of $\epsilon$, for $\epsilon\beta = 2$, and $n = 1$ (yellow dots), $n = 2$ (orange dots), $n = 3$ (red dots), $n = 4$ (purple dots), $n = 5$ (black dots). (b) Zoom in of the plot of $T_{+,5}/(-\beta\Delta/2)^5$.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

spin-boson model gives (40) in the basis of $\sigma_z$, $\{|+\rangle, |-\rangle\}$. with In the regime where $\omega'_S\beta \ll 1$, we obtain

$$ \rho_{MF} = \frac{1}{2} \left(1 - \frac{\epsilon\beta}{2}\right)|+\rangle\langle+| \\ +\frac{1}{2} \left(1 + \frac{\epsilon\beta}{2}\right)|-\rangle\langle-| \\ -\frac{\Delta\beta}{4} e^{-2\beta\lambda^2 Q/3} (|+\rangle\langle-| + |-\rangle\langle+|). \quad (81) $$

For comparison, if we take the limit $\beta\lambda^2 Q \ll 1$ (and $\epsilon^2\beta/\lambda^2 Q \ll 1$) in the expression (33) of $f_{+,-}(\beta)$, we obtain, using $DF(x) = x + O(x^3)$ for $x \ll 1$,

$$ \rho_S^{ss} = \frac{1}{e^{-\epsilon\beta/2} + e^{\epsilon\beta/2}} \left( e^{-\epsilon\beta/2} |+\rangle\langle+| + e^{\epsilon\beta/2} |-\rangle\langle-| \right) \\ -\frac{\Delta\beta}{4} \left( 1 + \frac{\epsilon^2\beta}{8\lambda^2 Q} \right) \sigma_x. \quad (82) $$

Both expressions are equivalent when applying $\beta\lambda^2\beta \ll 1$ to (81). However, in the opposite limit, when $\beta\lambda^2 Q \gg 1$, the expression (33) becomes

$$ \rho_S^{ss} = \frac{1}{e^{-\epsilon\beta/2} + e^{\epsilon\beta/2}} \left( e^{-\epsilon\beta/2} |+\rangle\langle+| + e^{\epsilon\beta/2} |-\rangle\langle-| \right) \\ -\frac{\Delta}{8\lambda^2 Q} \sigma_x. \quad (83) $$

This is significantly different from (40). While in the above expression the coherences vanish as $1/\lambda^2 Q$, they vanish exponentially in (40) and (81).

## 7.5 Without initial renormalization

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

$$ \mathcal{H}_{SB} = H_S + H_B + \lambda AB, \quad (84) $$

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

$$ e^{-\beta H_{SB}} = e^{-\beta(H_B + \lambda AB)} e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S(u)}, \quad (85) $$

$$ \tilde{H}_S(u) := e^{u(H_B + \lambda AB)} H_S e^{-u(H_B + \lambda AB)}. \quad (86) $$

We now have

$$ \begin{aligned} H_B + \lambda AB &= \sum_k \omega_k \mathcal{D}_k b_k^\dagger b_k \mathcal{D}_k^\dagger - \lambda^2 QA^2 \\ &= \sum_l |a_l\rangle\langle a_l| \mathcal{H}_{B,l}, \end{aligned} \quad (87) $$

with $\mathcal{H}_{B,l} := H_{B,l} - \lambda^2 a_l^\dagger Q = H_B + \lambda a_l B$, leading to

$$ \begin{aligned} e^{u(H_B + \lambda AB)} &= e^{u \sum_l |a_l\rangle\langle a_l| \mathcal{H}_{B,l}} \\ &= \sum_l |a_l\rangle\langle a_l| e^{u \mathcal{H}_{B,l}}. \end{aligned} \quad (88) $$

$$ \begin{aligned} \tilde{H}_S(u) &= \sum_{l,l'} |a_l\rangle\langle a_l| e^{uH_{B,l}} H_S |a_{l'}\rangle\langle a_{l'}| e^{-uH_{B,l'}} \\ &= H_S^{\text{pop}} + \tilde{H}_S^{\text{coh}}(u), \end{aligned} \quad (89) $$

where $H_S^{\text{pop}} = H_S^{\text{pop}} = \sum_l h_l |a_l\rangle\langle a_l|$ and $\tilde{H}_S^{\text{coh}}(u) := \sum_{l \neq l'} h_{l,l'} |a_l\rangle\langle a_{l'}| e^{u(H_{B,l}-\lambda^2 a_l^\dagger Q)} e^{-u(H_{B,l'}-\lambda^2 a_{l'}^\dagger Q)} = \sum_{l \neq l'} h_{l,l'} |a_l\rangle\langle a_{l'}| e^{u\mathcal{H}_{B,l}} e^{-u\mathcal{H}_{B,l'}}$. Again, similarly as previously, we obtain,

$$ e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S(u)} = e^{-\beta H_S^{\text{pop}}} e^{-\mathcal{T} \int_0^\beta du \tilde{H}_S^{\text{coh}}(u)}, \quad (90) $$

with

$$ \begin{aligned} \tilde{\mathcal{H}}_S^{\text{coh}}(u) &= e^{uH_S^{\text{pop}}} \tilde{\mathcal{H}}_S^{\text{coh}}(u) e^{-uH_S^{\text{pop}}} \\ &= \sum_{l \neq l'} h_{l,l'} e^{u\omega_{l,l'}} |a_l\rangle\langle a_{l'}| e^{u\mathcal{H}_{B,l}} e^{-u\mathcal{H}_{B,l'}}, \end{aligned} \quad (91) $$

arriving at

$$ \rho_{SB}^{\text{th}} = Z_{SB}^{-1} \sum_{l} e^{-\beta h_{l}} |a_{l}\rangle\langle a_{l}| e^{-\beta H_{B,l}} e^{-\mathcal{T} \int_{0}^{\beta} du \tilde{H}_{S}^{\text{coh}}(u)}. \quad (92) $$
Finally, the main change is that we are led to compute $\mathrm{Tr}_B[e^{-\beta H_{B,l}}]$ for the zeroth order, and $\mathrm{Tr}_B[e^{-\beta H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l}'}]$ for the second order, instead of $\mathrm{Tr}_B[e^{-\beta H_{B,l}}e^{uH_{B,l}}e^{-uH_{B,l}'}]$. We have,

$$ 
\begin{aligned}
\mathrm{Tr}_B & [e^{-\beta H_{B,l}} e^{uH_{B,l}} e^{-uH_{B,l}'}] \\
&= \mathrm{Tr}_B [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}'-\lambda^2 a_{l'}^2 Q)}] \\
&= \mathrm{Tr}_B [e^{-(\beta-u)H_B} e^{-u(H_B+\lambda a_{l'} B+\lambda^2 a_{l'}^2 Q)}] \\
&\quad \times e^{\lambda^2 Q [\beta a_l^2 + u(a_{l'}^2 - a_l^2)]}.
\end{aligned} 
\quad (93) $$

With that, we obtain

$$ \rho_S^{ss} = \sum_l \mathbf{p}_l^{ss} |a_l\rangle\langle a_l| - \sum_{l,l';l \neq l'} \mathbf{p}_l^{ss} h_{l,l'} \mathbf{f}_{l,l'}(\beta) |a_l\rangle\langle a_{l'}|, \quad (94) $$

where $\mathbf{p}_l^{ss} := e^{-\beta(h_l - a_l^2 \lambda^2 Q)}/Z_S^{ss} = e^{-\beta \mathbf{h}_l}/Z_S^{ss}$ is the renormalized population, with $\mathbf{h}_l := h_l - a_l^2 \lambda^2 Q$ the renormalized “pseudo-energies” (diagonal elements of $H_S$ in the eigenbasis of $A$), $Z_S^{ss} := \sum_l e^{-\beta \mathbf{h}_l}$, and

$$ \mathbf{f}_{l,l'}(\beta) := \int_0^\beta du \times e^{u\bar{\omega}_{l,l'}} e^{-\lambda^2 a_{l'}^2} \int_0^\infty d\omega \frac{J_l(\omega) e^{\omega\beta/2 + e^{-\omega\beta/2} - e^{\omega(u-\beta/2)} - e^{\omega(u-\beta/2)}}}{\omega^2 e^{\omega\beta/2} - e^{-\omega\beta/2}}, \quad (95) $$

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

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

### 7.5.1 Approximation of $\mathbf{f}_{l,l'}(\beta)$

Similarly as for $\mathbf{f}_{l,l'}(\beta)$, when the bath spectral density is such that $J(\omega)$ vanishes for $\omega \ge \omega_c$, where $\omega_c \le \beta^{-1}$, we have

$$ 
\begin{aligned}
\mathbf{f}_{l,l'}(\beta) &\approx \int_0^\beta du e^{u\bar{\omega}_{l,l'}} e^{-\lambda^2 a_{l'}^2} u(1-\frac{u}{\beta})Q \\
&= \frac{1}{\lambda|a_{l',l}|} \sqrt{\frac{\beta}{Q}} \left\{ \mathrm{DF} \left[ \frac{1}{2\lambda|a_{l',l}|} \sqrt{\frac{\beta}{Q}} (\lambda^2 a_{l',l}^2 Q - \bar{\omega}_{l,l'}) \right] \right. \\
&\qquad \left. + e^{\beta\bar{\omega}_{l,l'}} \mathrm{DF} \left[ \frac{1}{2\lambda|a_{l',l}|} \sqrt{\frac{\beta}{Q}} (\lambda^2 a_{l',l}^2 Q + \bar{\omega}_{l,l'}) \right] \right\}. 
\end{aligned} 
\quad (96) $$

Note that $\lambda^2 a_{l',l}^2 Q - \bar{\omega}_{l,l'} = 2\lambda^2 Q a_l(a_l - a_{l'}) - \omega_{l,l'}$ and $\lambda^2 a_{l',l}^2 Q + \bar{\omega}_{l,l'} = 2\lambda^2 Q a_l(a_l' - a_l) + \omega_{l,l'}$. Thus, in the strong coupling regime when $\lambda^2 Q \gg \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,

$$ 
\begin{aligned}
\mathbf{f}_{l,l'}(\beta) &= \frac{1}{2\lambda^2 Q a_{l,l'}} \left( \frac{1}{a_l} - \frac{e^{\bar{\omega}_{l,l'}\beta}}{a_{l'}} \right) \\
&\quad + \frac{\omega_{l,l'}}{4\lambda^4 Q^2 a_{l,l'}^2} \left( \frac{1}{a_l^2} - \frac{e^{\bar{\omega}_{l,l'}\beta}}{a_{l'}^2} \right) + O[(\lambda^2 Q / \omega_{l,l'}')^{-3}].
\end{aligned} 
\quad (97) $$

## 7.6 Identity of first orders

Here, we compare our general first order result, Eqs. (20) and (21), with the general first order expression for the steady state coherences obtained in Eq. (56) of [44]. This expression was obtained through projection operator techniques [28, 59, 60], by choosing the operator $\mathcal{P}$ projecting the system's state onto the diagonal subspace (in the eigenbasis of the unperturbed Hamiltonian). Thus, the non-diagonal elements can be obtained by establishing and solving the dynamics associated with the complementary operator $1-\mathcal{P}$ (section IV.A. of [44]). The steady state coherences are then deduced by taking the time to infinity. One should keep in mind that the master equation derived in the strong-decoherence limit in [44] is valid when the bath spectral density is such that $\lim_{\omega \to +\infty} \frac{J(\omega)}{\omega}$ is non-zero (and possibly infinite). This condition, satisfied by usual spectral densities, guarantees that the unperturbed dynamics in [44] leads to full decoherence, which is the starting point of the projection technique used therein.

To make the comparison easier, we re-write in the following the results of [44] using the notations we have been using here. Additionally, in order to simplify the presentation, we consider a simple coupling of the form $V_l = \lambda AB$ instead of the general one $V_l = \sum_\alpha A_\alpha B_\alpha$ considered in [44]. With that, in the eigenbasis $\{|a_l\rangle\}_l$ of $A$, the result of [44] takes the form

$$ 
\begin{aligned}
\rho_{l,l'}^{\mathrm{ss,ME}} &:= \langle a_l | \rho_S^{\mathrm{ss,ME}} | a_{l'} \rangle \\
&= h_{l,l'} \left[ i \mathbf{p}_l^{\mathrm{ss}} \int_0^\infty d\tau e^{-\lambda^2 a_{l',l}^2 [G^*(\tau) - i\tau Q]} e^{-i\bar{\omega}_{l,l'} \tau} - i \mathbf{p}_{l'}^{\mathrm{ss}} \int_0^\infty d\tau e^{-\lambda^2 a_{l',l}^2 [G(\tau) + i\tau Q]} e^{-i\bar{\omega}_{l,l'} \tau} \right]
\end{aligned} 
\quad (98) $$

with $\mathbf{p}_l^{ss} := e^{-\beta(h_l - a_l^2 \lambda^2 Q)}/Z_S^{ss} = e^{-\beta \mathbf{h}_l}/Z_S^{ss}$ is the renormalized population introduced in Section 7.5, with $\mathbf{h}_l := h_l - a_l^2 \lambda^2 Q$ the renormalized “pseudo-energies”, and $Z_S^{ss} := \sum_l e^{-\beta \mathbf{h}_l}$. Additionally, we defined $G(\tau) := \int_0^\tau ds_1 \int_0^{s_1} ds_2 c_B(s_2)$, and

$$ c_B(s) = \mathrm{Tr}_B[e^{iH_B s} B e^{-iH_B s} B \rho_B^{\mathrm{th}}]. \quad (99) $$
The superscript “ME” refers to the master equation nature of the derivation [44].

Our result can be expressed as (using the results of the derivation with no initial renormalization of the Hamiltonian, detailed in Section 7.5, since it is the choice made in [44]),

$$ 
\begin{aligned}
\rho_{l,l'}^{ss} &:= \langle a_l | \rho_S^{ss} | a_{l'} \rangle \\
&= -h_{l,l'} \mathbf{p}_l^{ss} \int_0^\beta du e^{-\lambda^2 a_{l,l'}^2 [G(-iu)+uQ]} e^{u\bar{\omega}_{l,l'}} \\
&= -h_{l,l'} \mathbf{p}_{l'}^{ss} \int_0^\beta du e^{-\lambda^2 a_{l',l'}^2 [G(-iu)+uQ]} e^{u\bar{\omega}_{l',l}},
\end{aligned} 
\quad (100) $$

where we used in the second line the identity $\mathbf{p}_l^{ss} \mathbf{f}_{l,l'}(\beta) = \mathbf{p}_{l'}^{ss} \mathbf{f}_{l',l}(\beta)$, shown in Section 7.2. Now, let us consider the integral of the function $e^{-\lambda^2 a_{l',l}^2 [G(-iu)+uQ]} e^{u\bar{\omega}_{l',l}}$ extended to the complex plan along the contour C defined as $u: 0 \to \beta \to \beta+ix \to ix \to 0$, where x is a positive number that later will be taken to infinity. Since it is a closed contour, the integral is equal to zero (as long as $G(u)$ is an analytic function). Additionally, when the bath spectral density satisfies $\lim_{\omega \to +\infty} \frac{J(\omega)}{\omega^2} > 0$, the real part of the function $G(x)$ tends to $+\infty$ as $x$ increases. The condition $\lim_{\omega \to +\infty} \frac{J(\omega)}{\omega^2} > 0$ is precisely the condition of validity of the results in [44]. This is not a coincidence. The validity of the master equation in [44] relies on the full decoherence of the unperturbed dynamics, which is guaranteed as long as $G(x)$ tends to $+\infty$ as $x$ increases, which itself results in the condition $\lim_{\omega \to +\infty} \frac{J(\omega)}{\omega^2} > 0$ on the bath spectral density. Then, the integral of $e^{-\lambda^2 a_{l',l}^2 [G(-iu)+uQ]} e^{u\bar{\omega}_{l',l}}$ tends to zero on the segment $\beta + ix \to ix$ for $x \to +\infty$. As a result, we have

$$ 
\begin{aligned}
\rho_{l,l'}^{ss} &= h_{l,l'} \mathbf{p}_l^{ss} \lim_{x \to +\infty} \int_{\beta}^{\beta+ix} du e^{-\lambda^2 a_{l',l}^2 [G(-iu)+uQ]} e^{u\bar{\omega}_{l',l}} \\
&\quad + h_{l,l'} \mathbf{p}_{l'}^{ss} \lim_{x \to +\infty} \int_{+ix}^{0} du e^{-\lambda^2 a_{l',l}^2 [G(-iu)+uQ]} e^{u\bar{\omega}_{l',l}} \\
&= ih_{l,l'} \mathbf{p}_l^{ss} \int_0^{+\infty} d\tau e^{-\lambda^2 a_{l',l}^2 [G(\tau-i\beta)+(i\tau+\beta)Q]} e^{i(\tau+\beta)\bar{\omega}_{l',l}} \\
&\quad -ih_{l,l'} \mathbf{p}_{l'}^{ss} \int_0^{+\infty} d\tau e^{-\lambda^2 a_{l',l}^2 [G(\tau)+i\tau Q]} e^{i\tau\bar{\omega}_{l',l}} \\
&= ih_{l,l'} \mathbf{p}_l^{ss} \int_0^{+\infty} d\tau e^{-\lambda^2 a_{l',l}^2 [G^*(\tau)-i\tau Q]} e^{-i\tau\bar{\omega}_{l,l'}} \\
&\quad -ih_{l,l'} \mathbf{p}_{l'}^{ss} \int_0^{+\infty} d\tau e^{-\lambda^2 a_{l',l}^2 [G(\tau)+i\tau Q]} e^{-i\tau\bar{\omega}_{l,l'}}
\end{aligned} 
\quad (101) $$

which is precisely equal to the expression (98) of $\rho_{l,l'}^{ss,ME}$. Note that in the last line we used the identities $\mathbf{p}_l^{ss} e^{\beta\bar{\omega}_{l',l}} = \mathbf{p}_l^{ss}$ and $G(\tau - i\beta) + (i\tau + \beta)Q = G^*(\tau) - i\tau Q$ (as well as $\omega_{l',l} = -\omega_{l,l'}$). This concludes the proof that the

expression derived in this paper coincides with the expression obtained in [44], although they have been obtained from very different methods and starting point. The equivalence of the expressions is guaranteed as soon as the results in [44] are valid (namely, as soon as $\lim_{\omega \to +\infty} \frac{J(\omega)}{\omega^2} > 0$). As reminded in the main text, this is an important step forward since it proves that the mean force Gibbs state is indeed the actual steady state (at least up the first order) even when the system interacts strongly with the thermal bath.

## Acknowledgments

I would like to thank Anton Trusheckin for very interesting and helpful discussions on this topic. I am also grateful to the Editors for their constructive comments, participating in improving the quality of the paper.

## Funding

This work was supported by funding from the National Institute for Theoretical Physics (NITheP) of the Republic of South Africa.

## References

[1] A. W. Chin, S. F. Huelga, M. B. Plenio. Coherence and decoherence in biological systems: principles of noise-assisted transport and the origin of long-lived coherences. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 2012; **370**(1972):3638–3657. doi:10.1098/rsta.2011.0224.

[2] P. Ribeiro, V. R. Vieira. Non-Markovian effects in electronic and spin transport. *Physical Review B* 2015; **92**(10):100302. doi:10.1103/PhysRevB.92.100302.

[3] P. Strasberg, G. Schaller, T. L. Schmidt, M. Esposito. Fermionic reaction coordinates and their application to an autonomous Maxwell demon in the strong-coupling regime. *Physical Review B* 2018; **97**(20):205405. doi:10.1103/PhysRevB.97.205405.

[4] L. A. Correa, B. Xu, B. Morris, G. Adesso. Pushing the limits of the reaction-coordinate mapping. *Journal of Chemical Physics* 2019; **151**(9):094107. doi:10.1063/1.5114690.
[5] S. V. Moreira, B. Marques, R. R. Paiva, L. S. Cruz, D. O. Soares-Pinto, F. L. Semião. Enhancing quantum transport efficiency by tuning non-Markovian dephasing. *Physical Review A* 2020; **101**(1):012123. doi:10.1103/PhysRevA.101.012123.

[6] E. Zerah-Harush, Y. Dubi. Effects of disorder and interactions in environment assisted quantum transport. *Physical Review Research* 2020; **2**(2):023294. doi:10.1103/PhysRevResearch.2.023294.

[7] D. Dwiputra, F. P. Zen. Environment-assisted quantum transport and mobility edges. *Physical Review A* 2021; **104**(2):022205. doi:10.1103/PhysRevA.104.022205.

[8] N. Anto-Sztrikacs, D. Segal. Strong coupling effects in quantum thermal transport with the reaction coordinate method. *New Journal of Physics* 2021; **23**(6):063036. doi:10.1088/1367-2630/ac02df.

[9] D. Gelbwaser-Klimovsky, A. Aspuru-Guzik. Strongly coupled quantum heat machines. *Journal of Physical Chemistry Letters* 2015; **6**(17):3477-3482. doi:10.1021/acs.jpclett.5b01404.

[10] P. Strasberg, G. Schaller, N. Lambert, T. Brandes. Nonequilibrium thermodynamics in the strong coupling and non-Markovian regime based on a reaction coordinate mapping. *New Journal of Physics* 2016; **18**(7):073007. doi:10.1088/1367-2630/18/7/073007.

[11] D. Newman, F. Mintert, A. Nazir. Performance of a quantum heat engine at strong reservoir coupling. *Physical Review E* 2017; **95**(3):032139. doi:10.1103/PhysRevE.95.032139.

[12] M. Perarnau-Llobet, H. Wilming, A. Riera, R. Gallego, J. Eisert. Strong coupling corrections in quantum thermodynamics. *Physical Review Letters* 2018; **120**(12):120602. doi:10.1103/PhysRevLett.120.120602.

[13] M. Wertnik, A. Chin, F. Nori, N. Lambert. Optimizing co-operative multi-environment dynamics in a dark-state-enhanced photosynthetic heat engine. *Journal of Chemical Physics* 2018; **149**(8):084112. doi:10.1063/1.5040898.

[14] D. Newman, F. Mintert, A. Nazir. Quantum limit to nonequilibrium heat-engine performance imposed by strong system-reservoir coupling. *Physical Review E* 2020; **101**(5):052129. doi:10.1103/PhysRevE.101.052129.

[15] M. Wiedmann, J. T. Stockburger, J. Ankerhold. Non-Markovian dynamics of a quantum heat engine: out-of-equilibrium operation and thermal coupling control. *New Journal of Physics* 2020; **22**(3):033007. doi:10.1088/1367-2630/ab725a.

[16] L. A. Correa, M. Perarnau-Llobet, K. V. Hohannisyan, S. Hernández-Santana, M. Mehboudi, A. Sanpera. Enhancement of low-temperature thermometry by strong coupling. *Physical Review A* 2017; **96**(6):062103. doi:10.1103/PhysRevA.96.062103.

[17] M. Mehboudi, A. Lampo, C. Charalambous, L. A. Correa, M. Á. García-March, M. Lewenstein. Using polarons for sub-nK quantum nondemolition thermometry in a Bose-Einstein condensate. *Physical Review Letters* 2019; **122**(3):030403. doi:10.1103/PhysRevLett.122.030403.

[18] M. Salado-Mejía, R. Román-Ancheyta, F. Soto-Eguibar, H. M. Moya-Cessa. Spectroscopy and critical quantum thermometry in the ultrastrong coupling regime. *Quantum Science and Technology* 2021; **6**(2):025010. doi:10.1088/2058-9565/abdca5.

[19] A. Kolli, E. J. O'Reilly, G. D. Scholes, A. Olaya-Castro. The fundamental role of quantized vibrations in coherent light harvesting by cryptophyte algae. *Journal of Chemical Physics* 2012; **137**(17):174109. doi:10.1063/1.4764100.

[20] N. Lambert, Y.-N. Chen, Y.-C. Cheng, C.-M. Li, G.-Y. Chen, F. Nori. Quantum biology. *Nature Physics* 2013; **9**(1):10–18. doi:10.1038/nphys2474.

[21] G. D. Scholes, G. R. Fleming, L. X. Chen, A. Aspuru-Guzik, A. Buchleitner, D. F. Coker, G. S. Engel, R. van Grondelle, A. Ishizaki, D. M. Jonas, J. S. Lundeen, J. K. McCusker, S. Mukamel, J. P. Ogilvie, A. Olaya-Castro, M. A. Ratner, F. C. Spano, K. B. Whaley, X. Zhu. Using coherence to enhance function in chemical and biophysical systems. *Nature* 2017; **543**(7647):647–656. doi:10.1038/nature21425.

[22] N. Lambert, T. Raheja, S. Cross, P. Menczel, S. Ahmed, A. Pitchford, D. Burgarth, F. Nori. QuTiP-BoFiN: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics 2020; arXiv:2010.10806.
[23] J. Iles-Smith, N. Lambert, A. Nazir. Environmen-
tal dynamics, correlations, and the emergence of
noncanonical equilibrium states in open quantum
systems. *Physical Review A* 2014; **90**(3):032114.
doi:10.1103/PhysRevA.90.032114.

[24] A. Purkayastha, G. Guarnieri, M. T. Mitchison, R. Filip, J. Goold. Tunable phonon-induced steady-state coherence in a double-quantum-dot charge qubit. *npj Quantum Information* 2020; **6**(1):27. doi:10.1038/s41534-020-0256-6.

[25] J. D. Cresser, J. Anders. Weak and ultrastrong cou-
pling limits of the quantum mean force Gibbs state.
*Physical Review Letters* 2021; **127**(25):250601.
doi:10.1103/PhysRevLett.127.250601.

[26] C. L. Latune. Steady state in strong system-bath coupling regime: Reaction coordinate versus perturbative expansion. *Physical Review E* 2022; **105**(2):024126. doi:10.1103/PhysRevE.105.024126.

[27] E. C. G. Sudarshan, P. M. Mathews, J. Rau. Stochas-
tic dynamics of quantum-mechanical systems. *Phys-
ical Review* 1961; **121**(3):920–924. doi:10.1103/
PhysRev.121.920.

[28] H.-P. Breuer, F. Petruccione. The Theory of Open Quantum Systems. Oxford University Press, Oxford, 2007. doi:10.1093/acprof:oso/9780199213900.001.0001.

[29] A. Garg, J. N. Onuchic, V. Ambegaokar. Effect of friction on electron transfer in biomolecules. *Journal of Chemical Physics* 1985; **83**(9):4491–4503. doi:10.1063/1.449017.

[30] J. Iles-Smith, A. G. Dijkstra, N. Lambert, A. Nazir. Energy transfer in structured and unstructured environments: Master equations beyond the Born–Markov approximations. *Journal of Chemical Physics* 2016; **144**(4):044110. doi:10.1063/1.4940218.

[31] B. M. Garraway. Nonperturbative decay of an atomic system in a cavity. *Physical Review A* 1997; **55**(3):2290–2303. doi:10.1103/PhysRevA.55.2290.

[32] G. Pleasance, B. M. Garraway. Application of quantum Darwinism to a structured environment. *Physical Review A* 2017; **96**(6):062105. doi:10.1103/PhysRevA.96.062105.

[33] A. E. Teretenkov. Pseudomode approach and vibronic non-Markovian phenomena in light-harvesting complexes. *Proceedings of the Steklov Institute of Mathematics* 2019; **306**(1):242–256. doi:10.1134/S0081543819050201.

[34] G. Pleasance, B. M. Garraway, F. Petruccione. Generalized theory of pseudomodes for exact descriptions of non-Markovian quantum processes. *Physical Review Research* 2020; **2**(4):043058. doi:10.1103/PhysRevResearch.2.043058.

[35] Y. Tanimura. Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM). *Journal of Chemical Physics* 2020; **153**(2):020901. 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. *Nature Communications* 2019; **10**(1):3721. doi:10.1038/s41467-019-11656-1.

[37] V. Bach, J. Fröhlich, I. M. Sigal. Return to equi-
librium. *Journal of Mathematical Physics* 2000;
**41**(6):3985-4060. doi:10.1063/1.533334.

[38] J. Fröhlich, M. Merkli. Another return of “return to equilibrium”. *Communications in Mathematical Physics* 2004; **251**(2):235–262. doi:10.1007/s00220-004-1176-6.

[39] M. Merkli, I. M. Sigal, G. P. Berman. Decoherence and thermalization. *Physical Review Letters* 2007; **98**(13):130401. doi:10.1103/PhysRevLett.98.130401.

[40] T. Mori, S. Miyashita. Dynamics of the density ma-
trix in contact with a thermal bath and the quantum
master equation. *Journal of the Physical Society of
Japan* 2008; **77**(12):124005. doi:10.1143/JPSJ.
77.124005.

[41] M. Könenberg, M. Merkli. On the irreversible dynamics emerging from quantum resonances. *Journal of Mathematical Physics* 2016; **57**(3):033302. doi:10.1063/1.4944614.

[42] M. Merkli. Quantum Markovian master equations: Resonance theory shows validity for all time scales. *Annals of Physics* 2020; **412**:167996. doi:10.1016/j.aop.2019.167996.

[43] A. S. Trushechkin, M. Merkli, J. D. Cresser, J. Anders. Open quantum system dynamics and the mean force Gibbs state. *AVS Quantum Science* 2022; **4**(1):012301. doi:10.1116/5.0073853.
[44] A. Trushechkin. Quantum master equations and steady states for the ultrastrong-coupling limit and the strong-decoherence limit. *Physical Review A* 2022; **106**(4):042209. doi:10.1103/PhysRevA.106.042209.

[45] Y. Subaşı, C. H. Fleming, J. M. Taylor, B. L. Hu. Equilibrium states of open quantum systems in the strong coupling regime. *Physical Review E* 2012; **86**(6):061132. doi:10.1103/PhysRevE.86.061132.

[46] P. L. Orman, R. Kawai. A qubit strongly interacting with a bosonic environment: Geometry of thermal states 2021; arXiv:2010.09201.

[47] K. Goyal, R. Kawai. Steady state thermodynamics of two qubits strongly coupled to bosonic environments. *Physical Review Research* 2019; **1**(3):033018. doi:10.1103/PhysRevResearch.1.033018.

[48] W. H. Zurek. Decoherence, einselection, and the quantum origins of the classical. *Reviews of Modern Physics* 2003; **75**(3):715-775. doi:10.1103/RevModPhys.75.715.

[49] E. K. Irish, J. Gea-Banacloche, I. Martin, K. C. Schwab. Dynamics of a two-level system strongly coupled to a high-frequency quantum oscillator. *Physical Review B* 2005; **72**(19):195410. doi:10.1103/PhysRevB.72.195410.

[50] Z.-M. Li, M. T. Batchelor. Generalized adiabatic approximation to the quantum Rabi model. *Physical Review A* 2021; **104**(3):033712. doi:10.1103/PhysRevA.104.033712.

[51] C. K. Lee, J. Moix, J. Cao. Accuracy of second order perturbation theory in the polaron and variational polaron frames. *Journal of Chemical Physics* 2012; **136**(20):204120. doi:10.1063/1.4722336.

[52] A. Kolli, A. Nazir, A. Olaya-Castro. Electronic excitation dynamics in multichromophoric systems described via a polaron-representation master equation. *Journal of Chemical Physics* 2011; **135**(15):154112. doi:10.1063/1.3652227.

[53] G. M. Timofeev, A. S. Trushechkin. Hamiltonian of mean force in the weak-coupling and high-temperature approximations and refined quantum master equations. *International Journal of Modern Physics A* 2022; **37**(20-21):2243021. doi:10.1142/s0217751x22430217.

[54] A. Gelzinis, L. Valkunas. Analytical derivation of equilibrium state for open quantum system. *Journal of Chemical Physics* 2020; **152**(5):051103. doi:10.1063/1.5141519.

[55] V. May, O. Kühn. Charge and Energy Transfer Dynamics in Molecular Systems. 3rd Edition. Wiley-VCH, Weinheim, Germany, 2011. doi:10.1002/9783527633791.

[56] R. P. Feynman. An operator calculus having applications in quantum electrodynamics. *Physical Review* 1951; **84**(1):108-128. doi:10.1103/PhysRev.84.108.

[57] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, W. Zwerger. Dynamics of the dissipative two-state system. *Reviews of Modern Physics* 1987; **59**(1):1-85. doi:10.1103/RevModPhys.59.1.

[58] W. H. Louisell. Radiation and Noise in Quantum Electronics. McGraw Hill, New York, 1964.

[59] S. Nakajima. On quantum theory of transport phenomena: steady diffusion. *Progress of Theoretical Physics* 1958; **20**(6):948-959. doi:10.1143/ptp.20.948.

[60] R. Zwanzig. Ensemble method in the theory of irreversibility. *Journal of Chemical Physics* 1960; **33**(5):1338-1341. doi:10.1063/1.1731409.