The dynamical convergence of a system to the thermal distribution, or Gibbs state, is a standard assumption across all of the physical sciences. The Gibbs state is determined just by temperature and energies of the system. However, at decreasing system sizes, i.e., for nanoscale and quantum systems, the interaction with their environments is not negligible. The question then arises: Is the system's steady state still the Gibbs state? If not, how may the steady state depend on the interaction details? Here, we provide an overview of recent progress on answering these questions. We expand on the state of the art along two general avenues: First, we take the static point-of-view, which postulates the so-called mean force Gibbs state. This view is commonly adopted in the field of *strong coupling thermodynamics*, where modified laws of thermodynamics and nonequilibrium fluctuation relations are established on the basis of this modified state. Second, we take the dynamical point of view, originating from the field of *open quantum systems*, which examines the time-asymptotic steady state within two paradigms. We describe the mathematical paradigm, which proves *return to equilibrium*, i.e., convergence to the mean force Gibbs state, and then discuss a number of microscopic physical methods, particularly master equations. We conclude with a summary of established links between statics and equilibration dynamics and provide an extensive list of open problems. This comprehensive overview will be of interest to researchers in the wider fields of quantum thermodynamics, open quantum systems, mesoscopic physics, statistical physics, and quantum optics and will find applications whenever energy is exchanged on the nanoscale, from quantum chemistry and biology to magnetism and nanoscale heat management.

## I. INTRODUCTION

Our everyday experience tells us that a macroscopic system, which is brought in contact with a much larger thermal environment at temperature *T*, such as a cup of coffee in a room, itself reaches a steady state characterized by the environment's temperature. Statistical physics argues that such an equilibrium state is determined by the system energetics, given by the Hamiltonian $HS$, as well as the temperature. Classically, this equilibrium state is known as the thermal distribution, while for quantum systems, it is known as the Gibbs state:

where $kB$ is the Boltzmann constant and *Z* is the partition function that normalizes the density matrix *τ*.

Taking the Gibbs state as the equilibrium or thermodynamically “free” state is a central assumption in much recent research on nanoscale and quantum thermodynamics. For example, it forms the basis of thermodynamic resource theory^{1–4} and is assumed in “thermal operations,” which investigate the properties of CPTP maps^{5} that have the Gibbs state as their fixed point.^{6,7} However, there is a (potentially serious) inconsistency here—the Gibbs state assumption can be problematic exactly for these “small” systems, as we will discuss below.

When a nanoscale or quantum system interacts with its environment, such as a molecule with the surrounding solution^{8} or a quantum spin with the phononic modes within a material,^{9} the system-environment interaction energy can become comparable in size to the system's bare (or self) energy. This is due to the fact that the surface-to-volume ratio of smaller systems can be much higher than that of macroscopic systems, e.g., scaling as $R2/R3=1/R$ for a spherical system with radius *R*. For short range interactions, the surface size of the system determines the strength of the interaction with its environment, while the self-energy of the system usually scales with volume. This implies that, while negligible for macroscopic systems, the interaction energy is relevant for systems of decreasing size.^{8–11}

Conventional Gibbs state statistical physics, which makes a tacit assumption that this interaction is vanishingly weak, then no longer applies. This motivates the core questions addressed in this overview article. Let $\rho (t)$ be the system density matrix at time *t* and denote the system steady state as

We ask

(Q) Is the system steady state $\rho ss$ the Gibbs state τ?

If not, how does $\rho ss$ depend on the interaction details?

As we will see, relaxation without recurrences is an irreversible effect, which, mathematically speaking, can only happen if the dynamics has no “eternal oscillations.” Such convergence to a steady state is attributed to the environment. However, not all environments cause such irreversible effects. For example, an environment consisting of a single qubit cannot make another qubit relax to a steady state. We will call an environment that *can* induce the convergence of the system (S) to a steady state a “*bath*” (B). Baths must have certain properties (large size, infinitely many degrees of freedom, a continuum of energies, etc.), which will be detailed in Sec. III. An illustration of the dynamical approach to a stationary state, as in Eq. (2), is provided in Fig. 1 for a qubit. A discussion of the quality of agreement of several predictions discussed in Sec. II, with the numerically solved steady state, is provided in Sec. IV F.

### A. Dynamic and static points of view

There are two points of view to answering the questions (Q): In the *dynamic* point of view, one considers the dynamics of the system, which continuously interacts with an environment, beginning from an initial state $\rho SB(0)$. One then asks whether the reduced state of the system alone, $\rho (t):=trB[\rho SB(t)]$, stops changing (significantly) at late times *t*. If this is so, then one may define (one or more) system steady states $\rho ss$. It is often assumed that the initial $SB$ state is uncorrelated, $\rho SB(0)=\rho (0)\u2297\tau B$, where $\tau B:=e\u2212HB/kBT/trB[e\u2212HB/kBT]$ is the bath Gibbs state (bath Hamiltonian $HB$ and temperature *T*). Furthermore, it is common to make the so-called *Born approximation*, $\rho SB(t)\u2248\rho (t)\u2297\tau B$ for all times $t\u22650$, implying that any effect of the system on the reduced bath state can be neglected as well as any correlations that may built up between $S$ and $B$.

In the *static* point of view, one postulates that at the end of any equilibration process, the combined system+bath complex is in the global Gibbs state $\tau SB\u221de\u2212Htot/kBT$, where $Htot$ is the total (interacting) Hamiltonian and the temperature is the same as the bath's temperature at the beginning of the equilibration process. The system equilibrium state is then “simply” the reduced state of the global state, $\tau MF=trB[\tau SB]$, called the *mean force Gibbs state*.

Each avenue has its merits and limitations. For instance, many analytical results from the dynamics point of view are based on the analysis of master equations (MEs), which are approximations of the true system dynamics. MEs are a powerful and widely used tool for the assessment of dissipation of quantum systems. The applicability of most MEs requires weak—while not negligible (see Fig. 2)—system–bath interaction strengths and additional further approximations. Different approximations can lead to different steady states $\rho ss$. Furthermore, quite often, approximations are *chosen* such that the steady state of the ME becomes the standard Gibbs state *τ*, turning the question regarding the steady state somewhat on its head.

On the other hand, the static point of view has given rise to the subfield of *strong coupling thermodynamics*, see Fig. 2, concerned with building a thermodynamic framework that correctly includes the bath's fingerprint.^{8,11–39} Based on the definition of an effective system Hamiltonian in equilibrium, called the mean force Hamiltonian, a general theory has been established, which includes thermodynamic laws and stochastic fluctuation relations for out-of-equilibrium processes. While there is some debate about the non-uniqueness of the mean force Hamiltonian, the mean force Gibbs state $\tau MF$ is generally accepted as the (unique) formal equilibrium state. An explicit expression of this state in terms of system operators alone is, however, only known in a handful of cases. This leaves much of the difficulty of including the environment in the reduced system state unresolved. Moreover, answering whether and/or when the dynamical steady state(s) $\rho ss$ and the reduced global equilibrium state $\tau MF$ are identical has been addressed only relatively recently for a handful of settings.

In this article, we assemble results—originally reported in many individual papers—into a focused, comparative overview describing both the static and the dynamic points of view. We begin with detailing expressions for the static mean force Gibbs state (MFG state) in Sec. II, followed by discussing mathematical results on the dynamical return to equilibrium (RtE) in Sec. III. Section IV gives a brief summary of key results on the dynamical steady states of microscopic master equations and other dynamical methods, and how these compare to the MFG state. We conclude in Sec. V with a summary of the outlined state of the art on the link between dynamics and statics and end with a (rather long) list of open questions.

Before embarking on the above topics, we will first set out the general setting of open quantum systems and clarify the naming convention that we will use for various coupling regimes.

### B. General setting and coupling strength regimes

The starting point for describing an open systems is to view it as a subsystem of a bigger, *closed* bipartite system $SB$, consisting of the system and the remaining part, called the bath. Deciding which part of an interacting complex is $S$ and which is $B$ may seem somewhat arbitrary. Intuitively, the system is understood to consist of degrees of freedom (DoFs) that one can manipulate and/or measure, such as the position and momentum of a pendulum, while the bath consists of degrees of freedom that are uncontrolled, such as air molecules that dampen the motion of the pendulum. The two components $S$ and $B$ are not equal partners: the bath influences the thermodynamic and dynamical properties of the system significantly, while $S$ cannot move $B$ too far from its initial state. This is modeled by taking bath Hamiltonians $HB$, which have a continuum of energies, while system Hamiltonians $HS$ have only discrete energies. In models where the bath has a spatial structure (say, the bath consists of a gas of quantum particles allowed to move in a region $R\u2282\mathbb{R}3$ of position space), continuous bath energies arise in the limit of infinite volume ($R\u2192\mathbb{R}3$), see also Sec. III A.

The total Hamiltonian of a general $SB$ complex has the following form:

where $\lambda \u2208\mathbb{R}$ is a dimensionless coupling constant and $VSB$ is the system-bath interaction operator. The latter is generally of the form

where $X(j)$ and $B(j)$ are operators acting on the system and bath Hilbert spaces, respectively.

Note that throughout the text, we will subindex states for system+bath with $SB$, e.g., the global Gibbs state is *τ _{SB}*, and we will subindex states for the bath alone with $B$, but for the system states, we will drop the index $S$, e.g., the system Gibbs state is denoted as

*τ*. We will, however, keep the index

*S*for the system Hamiltonian $HS$. Furthermore, unless otherwise stated, we set $\u210f=1$.

One distinguishes several regimes related to the magnitude of the coupling strength *λ*. Usually, these are set by a comparison of typical system ($HS$) energy differences and energies associated with the interaction $VSB$. We note that the current naming convention used in the theory of open quantum systems' literature and the strong coupling thermodynamics literature are not uniform.

Here, we propose a unified naming of the regimes, see Fig. 2, for a visual illustration. Our definition of the regimes is as follows: The *weak coupling regime* describes the regime where *λ* is small, and perturbation theory in *λ*, usually to second order in *λ*, is justified. Weak coupling regime examples include many master equation derivations, see Sec. IV C. The *ultraweak coupling regime* is taken when all terms of order *λ* and higher are neglected. In the other extreme, the *ultrastrong coupling regime* is achieved when *λ* is very large, a perturbation expansion in $\lambda \u22121$ may be performed, and all orders $\lambda \u22121$ and higher can be neglected. In the *intermediate coupling regime*, *λ* cannot be considered either large or small. In this challenging regime, either non-perturbative methods or other kinds of approximations are required.

## II. STATICS

Here, we discuss the static point of view, which arises from equilibrium statistical mechanics.

The classical Boltzmann distribution and the quantum Gibbs state *τ* are justified by a number of arguments.^{40} Gibbs' original derivation^{41} (see also Refs. 40 and 42–44) considers sharing of energies between two systems in equilibrium and makes use of the equal probability postulate for microstates. Many modern approaches to deriving the canonical equilibrium state in the classical and quantum regime follow the maximum entropy principle.^{45,46} This approach is justified by the second law of thermodynamics, which introduces the concept of irreversible entropy production and with it a direction toward higher entropy states. *Gibbs state statistical physics*, see Fig. 2, emerges when the equilibrium state of a system with fixed energy $\u27e8U\u27e9$ at temperature *T* is taken to be the state that maximizes entropy under the fixed energy constraint. To evaluate this maximum, one needs to know the energy operator, i.e., a system Hamiltonian $HS$, and chose an expression for the entropy, usually taken to be the Shannon or von Neumann entropy.^{47} An implicit assumption is that neither the Hamiltonian nor the entropy functional depends on the properties of the equilibrium state, such as the temperature. Maximization introduces a Lagrange multiplier *β*, which, upon equating the average statistical energy with $\u27e8U\u27e9$, becomes $\beta =1/kBT$. The resulting state of maximum entropy takes the form of the *Gibbs state τ*.

Thermodynamically, the system is postulated to reach this equilibrium state when it has been in weak thermal contact with a bath at temperature *T* for a long enough time. Within Gibbs state statistical physics, no further explicit mention is made of any bath. The only impact that the bath is assumed to have on the system is that it determines its temperature *T* and that the energy of the system is subject to (statistical) fluctuations around the fixed mean value $\u27e8U\u27e9$.

### A. Mean force Gibbs state $\tau MF$

We now consider the system and bath complex, $SB$, to be in the global Gibbs state $\tau SB$ associated with the total Hamiltonian $Htot$ [Eq. (3)]. The emergence of this state can be justified—for now—by considering that the $SB$ compound has been in very weak thermal contact for a long time with a *super-bath*^{14,48} *R* at inverse temperature $\beta =1/kBT$. Gibbs state statistical physics for the compound $SB$ then tells us that the equilibrium state for $SB$ is the Gibbs state

where $ZSB=trSB[e\u2212\beta Htot]$ is the global partition function. The *mean force Gibbs state* (MFG state) is defined as the system state obtained by taking the partial trace over the bath degrees of freedom,

Generally, $\tau MF$ differs—sometimes substantially—from the Gibbs state $\tau \u221de\u2212\beta HS$, as we will see below. The naming arises from casting $\tau MF$ in the Gibbsian (exponential) form

for an effective Hamiltonian $HMF$ called the Hamiltonian of mean force (HMF) or the potential of mean force. The MFG state, as well as the HMF, has found widespread use in chemistry^{49–55} since the 1930s.

Before proceeding with the discussion of the state $\tau MF$, let us first comment on the HMF. Unlike the bare system Hamiltonian $HS$ used throughout Gibbs state statistical physics, the HMF is temperature dependent. (As a consequence, identities of statistical physics will not necessarily hold and require corrections.) It will also depend on the coupling strength *λ* and some of the details of the interaction $VSB$ in Eq. (3). Note that the $HMF$ is not uniquely defined since, e.g., a constant may be added to it without changing $\tau MF$ in Eq. (7). This is due to the fact that such a constant would also change the partition function $ZMF=trS[e\u2212\beta HMF]$. (The constant also cancels when calculating energetic differences). A common choice is to set $ZMF=ZSB/ZB$ with $ZB$ the bare bath partition function, and to include certain strong coupling corrections into energetic and entropic potentials.^{8,10,11,19,26,39} This leads to an extensive (additive) behavior of the effective system and bare bath potentials for classical and quantum systems, which mirrors that of standard thermodynamics. Other thermodynamically consistent choices are being discussed,^{8,39} and approaches to determining the physically meaningful HMF are being explored.^{39,56,57}

Meanwhile, based on the above definitions, much progress has been made in constructing a comprehensive framework of “strong coupling thermodynamics”^{11} that includes corrections arising from the system's interaction with the environment. Strong coupling thermodynamic potentials have been identified,^{8,19,21,23,28} detailed entropy fluctuation relations have been shown to hold,^{10,26} and the validity of the Jarzynski equality^{13,14,58} and the Clausius inequality^{15–17} has been proven. The strong coupling impact on a Maxwell demons' operation has been elucidated,^{29,30} an extension of Bohr's energy-temperature uncertainty relation to the strong coupling limit has been proven,^{31} and quantum measurements have been included in a stochastic description of strongly coupled quantum systems.^{33} In quantum thermometry, strong coupling has been found to improve measurement precision,^{27,34} while it can be detrimental for the efficiency of quantum engines.^{35}

We now return to the MFG state $\tau MF$, which is uniquely defined by the formal identity [Eq. (6)]. However, unfortunately, giving explicit expressions for $\tau MF$ in terms of system operators alone is very often intractable—because it requires carrying out the trace over the (large number of) bath DoFs. Exact results have been obtained for the quantum harmonic oscillator interacting with a bath of oscillators, see Sec. II D. For more general systems $S$, still interacting with a bosonic bath, perturbative results have been established in the weak coupling limit, see Sec. II E, as well as the ultrastrong coupling limit, see Sec. II F.

### B. Open systems with a discrete bosonic environment

The paradigmatic open quantum system model is a system with Hamiltonian $HS$ coupled to a field of quantum harmonic oscillators according to the Hamiltonian (3) with^{59,60}

where *ω _{k}* are the frequencies of the oscillator (modes)

*k*, the creation $ak\u2020$ and annihilation operators

*a*obey the commutation relations $[ak,a\u2113\u2020]=\delta k,\u2113$, and

_{k}*X*is an arbitrary system operator.

*g*are complex numbers that weigh the strength of the interaction between $S$ and the oscillator mode at frequency

_{k}*ω*. Even though $HB$, Eq. (8), has infinitely many energy levels, those levels do not fill a continuum of values. Thus, technically, according to our bath definition, $HB$ is

_{k}*not*the Hamiltonian of a “bath.” Nevertheless, the discrete mode model [Eq. (8)] often serves as a starting point, see also Sec. II C. An additional, so-called counter term $\lambda 2\u2211k|gk|2X2/\omega k$ is often included in the Hamiltonian, which physically arises whenever coupling is introduced via the difference of coordinates, for example, $(xk\u2212X)2$, instead of a product, e.g., $xkX$. When included, the total Hamiltonian is

where the bath oscillators are now *displaced* by the system operators. Note that we here neglect the zero-point energy of the bath oscillators, i.e., $\u2211k\omega k/2$, which cancels in the reduced system state. This energy diverges in the continuum mode limit, and dropping it amounts to a renormalization of the bath energy.

When the system is a two level system, then Eq. (8) is the Hamiltonian of the ubiquitous *spin-Boson model*, which is used to describe a wide range of physical systems,^{59,61–66} ranging from qubits in quantum computers^{7,67,68} to electron transfer complexes in quantum chemistry and biology.^{69–72}

For a particle moving in an arbitrary potential *v*(*x*) with *x* the particle position operator, Eq. (9) is the well-known *Caldeira–Leggett (CL) model* of quantum Brownian motion.^{73} Here, the inclusion of the counterterm^{74} guarantees that the particle dynamics given by the Heisenberg equation of motion for *x* is determined by the bare potential *v*(*x*) and not by a renormalized potential $v(x)\u2212\lambda 2\u2211k|gk|2X2/\omega k$. Mathematically, this is significant since, without the counterterm and at sufficiently strong system–reservoir interaction, the energy spectrum of the global system can become unbounded from below leading to a thermodynamically unstable scenario.^{75,76} Of particular interest is the quantum harmonic oscillator model for which $v(x)=m\omega 02x2/2$, giving the Hamiltonian^{77}

where *x _{k}* and

*p*are the bath oscillator position and momentum operators, respectively, $gk=\u2212ck\u2009\u210f/(m\u2009mk\omega k)$, and the coupling operator is $X=x\u2009m/2$. The Caldeira–Leggett model is a key model for open quantum systems theory

_{k}^{21,78,79}with widespread application to quantum tunnelling

^{80}and studies of decoherence and the quantum-classical transition.

^{78}

### C. Continuum limit for bosonic baths

For the open system to actually exhibit irreversibility, one must take a bath with a continuous spectrum, as we will discuss in detail in Sec. III D on return to equilibrium. In addition, there is a practical advantage of taking the continuum limit: it means replacing sums with integrals and with it converting some intractable summations into analytically solvable integrals.

For bosonic baths, one wants to replace the discrete set ${\omega k}k$ by a continuum of frequencies, for instance, $\omega \u2208[0,\u221e)$. There are two common ways of implementing this continuum limit. An *ad hoc* way is to perform the limit in expressions for specific physical quantities (such as time-dependent population probabilities and coherences), which are obtained from calculations using a discrete mode model, such as Eq. (8). In this approach, the state of the continuous mode model is actually never constructed. The procedure may be rather easily implementable, but it has some disadvantages. For instance, often one has to consider the limits of continuous modes (infinite volume), small/large coupling and large time “simultaneously,” and it is not possible to control those limits in this setup.

The second approach is to immediately construct the continuous mode model^{9,81,82,173,176,275} and then analyze the full $SB$ statics and dynamics. This allows one, in particular, to control the perturbation theory for *all times*, even $t\u2192\u221e$. This is done in the quantum resonance theory, which we explain in Sec. IV C 1. In quantum optics models, the index *k* labeling the oscillators in Eq. (8) represents a wave vector in physical space of dimension *d* (usually, *d *=* *3).^{83} The continuous mode limit then leads to $k\u2208\mathbb{R}d$, and the continuous mode Hamiltonian associated with Eq. (8), becomes^{9,81}

where $a\u2020(g)=\u222b\mathbb{R}ddk\u2009gk\u2009ak\u2020$ is the creation operator smoothed out with *g _{k}*, which in this context is sometimes called the

*form factor*, a square-integrable complex function of $k\u2208\mathbb{R}d$. Taking the continuum limit here amounts to taking the quantization volume of the problem to infinity, which also redefines the creation and annihilation operators. In the continuous mode limit, they obey the continuous canonical commutation relations, $[ak,a\u2113\u2020]=\delta (k\u2212\u2113)$, where

*δ*is the Dirac delta-function in

*d*dimensions. The index

*k*does not need to have a physical meaning though, generally, beyond simply being a continuous index labeling the modes.

^{9,81}Often it is directly chosen to be the energy

*ω*of a mode, which means that $HB\u221d\u222b0\u221ed\omega \u2009\omega \u2009a\omega \u2020a\omega $ and $VSB\u221dX\u2297\u222b0\u221ed\omega \u2009g\omega \u2009a\omega \u2020+h.c.$

For a discrete model, an alternative to specifying the coupling constants *g _{k}* is to specify the (real) bath spectral density,

^{59,84,85}

a choice that allows modeling diverse physical situations. For continuous mode environments, the sum is naturally replaced by an integral $\u222bdk$.

If $J(\omega )\u221d\omega s$ for small *ω*, then the spectral density is called “Ohmic” for *s *=* *1, and “super-Ohmic” for *s *>* *1. In what follows, we will assume that the spectral density is either Ohmic or super-Ohmic. In both cases, the limit $J(\omega )/\omega $ as $\omega \u21920$ is finite. This will be important for obtaining finite damping rates in Secs. IV C 1 and IV C 2, which turn out proportional to $J(\omega )coth(\beta \omega /2)\u221dJ(\omega )/\omega $ at low *ω*. The sub-Ohmic case *s *<* *1 is studied by different theoretical methods, see, e.g., Refs. 86–93, and non-Ohmic densities have been found to be relevant in, e.g., opto-mechanical resonator experiments.^{94}

### D. Exactly solvable $\tau MF$ for the Caldeira–Leggett model

The properties of the MFG state for the damped quantum harmonic oscillator given by the Caldeira–Leggett Hamiltonian [Eq. (10)], for arbitrary coupling strengths *λ*, were obtained in Refs. 74, 95, and 96. For simplicity, let us formally put here *λ* = 1 (one can consider *λ* to be included in the coupling coefficients *c _{k}*). The most explicit results can be obtained in the case of the Drude–Lorentz spectral density,

in the limit of continuous modes. Here, $\omega D$ is the Drude frequency (which determines the time scale of the bath relaxation) and *γ* is a damping frequency. Since the CL model is quadratic, the system MFG state will be of Gaussian form and completely determined by the first and second moments^{97} of the oscillator position and momentum operators: $\u27e8x\u27e9,\u2009\u27e8p\u27e9,\u2009\u27e8x2\u27e9,\u2009\u27e8p2\u27e9$, and $\u27e8px\u27e9$. The first moments trivially vanish, while the second moments are given in terms of the partition function $Z(\omega 0,\gamma )$ at inverse temperature *β*,

where $\Gamma (z)$ denotes the gamma function and $\nu =2\pi /\beta $ is the first *Matsubara frequency*. The functions $\mu j(\omega 0,\gamma )$, for *j *=* *1, 2, 3, denote the roots of the cubic polynomial

With these expressions, the second moments given in unit-free form [$x\u0303$ and $p\u0303$, see Eq. (10)], and including $\u210f$ explicitly, are^{95}

Reference 95 implies the MFG state as the Gibbs state of an effective harmonic oscillator Hamiltonian (the HMF) $HSeff=p2/(2m(\lambda ))+m(\lambda )\omega 02(\lambda )x2/2$. They formally establish the mass $m(\lambda )$ and frequency $\omega 0(\lambda )$ of the oscillator, which are rescaled from their bare counterparts due to the interaction with the bath. Alternatively, the harmonic oscillator's MFG state can also be given in Gaussian integral form,^{97} expanded in the Weyl-basis, as

It is an open question to show if this expression indeed reduces to the Gibbs state of the effective Hamiltonian given in Ref. 95.

Another method of obtaining the oscillator moments, cf. Eq. (16), is via the closed expressions for the oscillator correlation functions in the MFG state. These have been derived^{21,98} using Heisenberg equations of motion for an arbitrary spectral density $J(\omega )$, e.g.,

where the mass was set to *m *=* *1 and $\u210f=1$ again. Here, $G(\omega )=\u22121/(\omega 2\u2212\omega 02[1\u2212\chi (\omega )])$ is a Green's function, and $\chi (\omega )$ is a complex susceptibility given by

where the integral is understood as the principal part integral. Setting $t=t\u2032$ makes the correlation function [Eq. (18)] time-independent, and it should then correspond to the first line in Eq. (16). This route allows us to give analytic expressions for thermal energies^{21} of damped quantum and classical harmonic oscillators as a functional of $\chi (\omega )$ and inverse temperature *β*.

### E. Expansion of $\tau MF$ for small coupling constant *λ*

*λ*

We are now interested in expressions for the MFG state in the weak coupling limit, when *λ* in Eq. (3) is assumed to be small, so that an expansion of $\tau MF$ can be considered. Such expansions are based on the Kubo identity^{99} and are sometimes referred to as “canonical perturbation theory,”^{100} equivalent to a standard time dependent perturbation expansion with *t* replaced by $\u2212i\beta $. In the mathematical literature, the expansion in *λ* is known as the perturbation theory of Kubo–Martin–Schwinger (KMS) states.^{101} For a bosonic bath and a linear $SB$ coupling as in Eq. (11), one can consider the expansion

Only even terms in *λ* are present because thermal equilibrium averages of products with any odd number of bath creation and annihilation operators vanish. Such expansions have been provided in Refs. 66, 82, 100, 102, 103, and 104 to second order in *λ*. In Ref. 104, the correction term $\tau MF(2)$ was obtained for the Caldeira–Leggett model [with an arbitrary potential *v*(*x*)] and in Ref. 66, for the spin-boson model.

Reference 82 considers an arbitrary system $S$ with system Hamiltonian $HS$ with a discrete spectrum coupled to a continuous bosonic bath via an arbitrary Hermitian system operator *X*. The Hamiltonian, including the counter term [cf. Eqs. (9) and (10)] is given by

Note that here a spectral density $J(\omega )$ is used immediately without specifying any coupling constants *c _{k}*. Here, $[q\omega ,p\omega \u2032]=i\u2009\delta (\omega \u2212\omega \u2032)$ are the commutation relations for the bath position and momentum operators. Using perturbative methods and evaluating the bath trace in Eq. (6) explicitly, the dominant correction $\tau MF(2)$ for an arbitrary system was found to be

Here, the decomposition of the Hermitian system operator *X* into a sum of energy eigenoperators *X _{m}* is used, where

*X*are defined by

_{m}with *ω _{m}* being the Bohr frequencies of the system (energy differences of $HS$). Furthermore, the function $D\beta (\omega m)$ is defined as

where the integral is understood as a principal part integral. Expression (22) evidences the appearance of coherences in the $HS$ basis in the system's equilibrium state $\tau MF$ (see Fig. 1 and Sec. IV F). Coherences are often considered a quantum “resource,”^{105} and beyond their significance in quantum thermodynamics,^{66,106–110} play an important role in some biological processes.^{111–115}

One can now also quantify what is “weak enough” for the weak coupling limit and expression [Eq. (22)] to be valid. Beyond the loose requirement that *λ* ought to be “small,” one finds (by comparing perturbative orders) that *λ* has to obey the inequality,^{82}

This condition gives a well-quantified limit for *λ* being in the weak coupling regime at a given *β*. Note that the range of *λ* for which the weak coupling regime and, hence, Eq. (22) is applicable changes as a function of temperature with larger temperature generally allowing larger *λ*.

### F. Ultrastrong coupling for general system and bosonic bath

One can also consider the opposite limit, when coupling is much stronger than other energy scales of the system, i.e., the “ultrastrong” coupling limit $\lambda \u2192\u221e$, see Fig. 2 and Refs. 116–120. Here, it is also possible^{82} to find an explicit expression for $\tau MF$ for a general system. This is done for the case that it couples to a bosonic bath as in Eq. (21) with a single system interaction operator *X* with a nondegenerate spectrum (extensions to degenerate situations should be straightforward),

where *x _{n}* are real numbers and

*P*are orthogonal projectors of rank one. By expanding the global Gibbs state $\tau SB$ to orders of $1/\lambda $ and explicitly integrating out the bath oscillators, the mean force Gibbs state simplifies to

_{n}^{82}

This is a surprisingly neat form for the open system equilibrium state at ultrastrong coupling. It implies that, in this limit, the equilibrium state of the system becomes diagonal in the basis of the system's coupling operator *X*. In the context of measurement and decoherence theory, it is referred to as the “pointer basis.”^{118,121–124} For $[HS,X]\u22600$, an immediate consequence is that $\tau MF$ will maintain coherences in the $HS$ energy basis $|em\u27e9$, i.e., $\u27e8em|\tau MF|em\u27e9\u22600$ for some *m* (see Fig. 1 and Sec. IV F). Corrections to Eq. (27) with respect to $\lambda \u22121$ have been obtained in Ref. 125.

It is worthwhile to build bridges between these results and discussions about localized and delocalized excitations in the theory of excitation energy transfer in biological photosynthetic complexes.^{69,126–128} Due to the dipole interaction between the chromophore molecules, the eigenstates of the system Hamiltonian are superpositions of *local excitations* and describe the so-called delocalized excitons. The *X* operator [or more generally $X(j)$ in Eq. (4)] is diagonal in the local excitation basis. Thus, the local excitation basis corresponds to the pointer basis.

In this context, “weak coupling” theory describes the relaxation in the basis of delocalized excitons. Non-negligible coupling to the phonon bath leads to the relaxation not in the exciton basis, but in a more localized basis.^{127} In the ultrastrong coupling limit, which corresponds to the Förster theory of excitation energy transfer,^{69} the relaxation occurs in the local excitation (pointer) basis. This is exactly the regime described by the mean force Gibbs state [Eq. (27)]. The dynamical aspects of the Förster regime of excitation energy transfer and the ultrastrong coupling regime for a general open quantum system will be discussed in Sec. IV D.

### G. Intermediate coupling: Polaron transformation

Under certain conditions, the ultrastrong and intermediate coupling regime, see Fig. 2(b), can be treated with the so-called *polaron transformation*. Originally, a polaron is a quantum quasiparticle in a solid material consisting of an electron and a field of elastic deformations of the crystal lattice (a phonon cloud).^{129} In the more general context of open quantum systems theory, a polaron is a state of the system “dressed” by the bath excitations. Mathematically, the polaron transformation is a certain unitary transformation acting on $SB$, which mixes the system and bath DoFs.^{130–134} The benefit of the polaron transformation is that one can apply weak coupling perturbation theory for the redefined system and bath; see also Sec. IV E 2.

As an illustration of this formalism, we consider the spin-boson model.^{67,135–137} The total Hamiltonian (21) here contains

and the system coupling operator is $X=\sigma z$, where $\sigma x,y,z$ are the usual Pauli matrices. Thus, the pointer basis is the *σ _{z}*-basis. Then, the polaron transformation is given by the unitary transformation

The polaron-transformed Hamiltonian (indicated with a tilde) is $H\u0303tot=UHtotU\u2020=H\u0303S+HB+V\u0303SB$, where $H\u0303S=\epsilon 2\sigma z+\kappa \u2009\Delta 2\sigma x$ with

The bath part remains unchanged, $HB=\u222b0\u221e\u2009\u2009d\omega \u2009\omega \u2009a\omega \u2020a\omega $, and the interaction becomes

where *λ* has now moved inside $V\u0303SB$. When the integral in the definition of *κ* converges, it turns out that $trB[\tau B\u2009e\xb1i2\lambda R\u0302]=\kappa $. However, convergence only happens for a subclass of super-Ohmic spectral densities. For example, the integral converges whenever for small *ω* and considering strictly positive temperatures, $J(\omega )$ is proportional to $\omega 3$, but it diverges whenever $J(\omega )$ is proportional to *ω ^{s}* for $s\u22642$. This is a restriction of the polaron transformation method.

The factor *κ* represents the above-mentioned “phonon cloud,” while the *B _{x}* and

*B*operators represent fluctuations around this cloud. One may hope that these fluctuations are not large and $V\u0303SB$ can be treated perturbatively. Thus, the benefit of the polaron transformation is that one can now apply the weak coupling perturbation theory for the rotated $SB$ complex. However, strictly speaking, the conjecture of applicability of the weak coupling theory to the polaron-transformed Hamiltonian is justified only in two opposite limits: (i) the weak system-bath limit [small $\lambda 2J(\omega )$], where the polaron transformation is trivially reduced to the identity transformation, and (ii) the limit of small Δ (weak tunneling limit).

_{y}^{137,138}Since $|\kappa \Delta |<|\Delta |$, the eigenvectors of $H\u0303S$ are more “localized” superpositions of the pointer basis vectors than the eigenvectors of $HS$. Such localization due to non-negligible system–bath interaction was mentioned at the end of Sec. II F.

Now, we can consider the total Gibbs state in the polaron frame $\tau \u0303SB=U\u2009\tau SB\u2009U\u2020\u221de\u2212\beta H\u0303tot$. One can show that the diagonal part (in the pointer basis) of the desired MFG state $\tau MF$ formally coincides with the diagonal part of the reduced system state $\eta \u0303=trB[\tau \u0303SB]$, i.e., $\tau MF,jj=\eta \u0303jj$ for *j *=* *1, 2. Approximate expressions for $\eta \u0303$ were obtained^{135–137} using second-order perturbation theory with respect to $V\u0303SB$,

where $\eta \u0303(0)=e\u2212\beta H\u0303S/Z\u0303$ with $Z\u0303=trS[e\u2212\beta H\u0303S]$ is the Gibbs state corresponding to the system Hamiltonian in the polaron frame. The next term in Eq. (32) is

where

and $Gmn(\beta )=trSB[Bm(\beta )Bn\tau B]$ are the imaginary-time bath correlation functions. The operators in imaginary time, such as $\sigma m(\beta )$ and $Bm(\beta )$, are defined as $O(\beta )=e\beta (H\u0303S+HB)Oe\u2212\beta (H\u0303S+HB)$. These expressions can be made more explicit using the methods of Ref. 82 presented in Sec. II E.

However, to determine the off-diagonal element $\tau MF,12$ (in pointer basis), $\tau \u0303MF$ does not suffice—the bath DoFs of the total polaron-transformed Gibbs state are also required. Up to the first order in $V\u0303SB$, it is found^{136,137} to be

where $\Lambda =\epsilon 2+(\kappa \Delta )2,\u2009Sm(\beta \u2032)=trS[\sigma m(\beta \u2032)\u2009\sigma \u2212\u2009\tau \u0303]$, and $Km(\beta \u2032)=trB[Vm(\beta \u2032)\u2009\u2009cos\u2009(2\lambda R\u0302)\u2009\tau B]$, where the functions *S _{m}* and

*K*can be explicitly evaluated and $\sigma \u2212=(\sigma x\u2212i\sigma y)/2$.

_{m}In Refs. 135–137, the following (super-Ohmic) spectral density is considered:

where *γ* (or, more precisely, $\lambda 2\gamma $) determines the system–bath coupling strength, while $\omega c$ is the cutoff frequency and determines the rate of relaxation of the bath correlation functions in time. The above expressions for the elements of $\tau MF$ are compared with numerically exact simulations. It turns out that the approximation works well for the cases of fast bath $\omega c>\Delta $ and ultrastrong coupling (large *λ*). It remains an open question to simplify expressions (32) for $\tau MF$ to a form similar to (22).

Finally, the so-called variational (partial) polaron transformation can be used to enlarge the range of applicability of this approach. In particular, the variational polaron transformation allows one to overcome the assumption of the super-Ohmic spectral density. Numerical calculations of the equilibrium state and equilibrium physical observables for the Ohmic spectral baths using the variational polaron approach are presented in Refs. 135–137 and 139.

### H. High-temperature expansion

An approximate expression for $\tau MF$ in the high temperature limit can be obtained Ref. 363 by expanding in powers of inverse temperature *β*.

In the notations introduced above, the model considered in Ref. 363 can be formulated as follows: a multi-state system with orthonormal basis ${|n\u27e9}$ interacts, via $Xn=|n\u27e9\u27e8n|$, with several baths $n=1,\u2026,N$ that all have the same temperature. Instead of (21), the Hamiltonian is

for independent baths,

The system Hamiltonian can be decomposed as $HS=H\epsilon +HJ$, where $H\epsilon $ is the part that is diagonal in the $|n\u27e9$ basis and *H _{J}* contains all off-diagonal contributions (cf. Ref. 63). Expanding Eq. (6) for the Hamiltonian Eq. (36) to second order in

*β*, evaluating the partial traces, and then re-summing into an exponential form, give a MFG state (Ref. 363) proportional to

Here, $\Lambda =\lambda 2\u2211n=1N\u222b0\u221ed\omega \u2009Jn(\omega )\omega \u2009Xn\u2261\u2211n=1N\u2113n\u2009Xn$ is an operator-valued reorganization term. The result [Eq. (37)] implies that the interstate coupling constants *J _{mn}* contained in

*H*, describing hopping between states $|m\u27e9$ and $|n\u27e9$, get rescaled by the bath interaction into effective coupling constants. The rescaling itself is temperature dependent and vanishes for $\beta \u21920$. For a dimer system with $\u21131=\u21132=\u2113$, expression (37) is found to be accurate (Ref. 363) for temperatures satisfying $\u2113\beta \u22722$. We emphasize that generally Eq. (37) is valid even at intermediate system–bath coupling strengths as long as the temperature is large enough (see Sec. IV F).

_{J}### I. Adequacy of the bosonic bath model

The bosonic bath model considered from Sec. II B onwards is very widely used in the theory of open quantum systems and quantum thermodynamics, though fermionic and spin bath models are also actively studied.^{140–147} In addition to it being bosonic, coupling to the system is assumed to be linear in creation and annihilation operators. In this case, the bath DoFs can, without loss of generality, be assumed to be noninteracting, as re-diagonalization can always bring it into such a normal mode form. Nevertheless, the resulting bath model is a very special case, and we here discuss its range of applicability.

One can distinguish three levels of validity of this model. The only case where this model is exact is for a bath consisting of photons (electromagnetic field).^{148–154} On the second level, the bosonic bath model is used as an approximation of the real physics. For example, in solid-state physics and chemistry, e.g., for modeling charge and energy transfer, the bath consists of phonons or vibrational modes that describe oscillatory degrees of freedom of nuclei. If the magnitudes of these oscillations are not too large, the harmonic approximation can be used^{153,155} implying linear coupling to the system (electronic degrees of freedom). This approximation is used for various system-bath (electron-phonon) coupling regimes.^{153} Even at strong coupling, it is often a reasonable assumption that the oscillations of the nuclei around the equilibria are small enough to warrant the harmonic approximation (though the equilibria themselves can be significantly shifted due to strong coupling).^{80}

The third level is the use of this model as a phenomenological model, which need not directly represent the real physics of the bath. Namely, let the bath be a complex system of many interacting particles that cannot be reduced to a set of harmonic oscillators. However, the details of the bath dynamics are not that important for the reduced dynamics of the system—only some aggregated properties of the bath dynamics, such as correlation functions, are required. In addition, Gibbs states of the bosonic bath are Gaussian, implying that all correlation functions can be expressed in terms of just the second-order correlation functions. Such a Gaussian property is likely to emerge also for rather general baths containing a large number of particles, as a consequence of the central limit theorem. Thus, there is hope that a bosonic bath with the same second-order correlation functions as that of a real bath may serve as a phenomenological model of the real bath, at least for qualitative analysis.

## III. RETURN TO EQUILIBRIUM

Return to equilibrium (RtE) is a basic and intuitive phenomenon, saying that initial states which do not deviate much from the equilibrium state will converge to the equilibrium state in the long time limit. As an analogy, a ripple created at some point on the surface of a still lake (equilibrium) will propagate away. Eventually, the lake's surface will return to be still. For this to happen, it seems clear that the total system under consideration has to be infinitely large to avoid recurrences for all times, and that the dynamics has to be dissipative in the sense that it propagates local disturbances away to infinity. Furthermore, the perturbation must be “small,” for instance localized in space (if initial ripples are created everywhere in space then at any fixed point, the surface will not remain still, even for large times, as ripples keep arriving from far away positions). In Sec. III A–III E, we formalize these intuitive notions in mathematical terms.

### A. Continuous spectrum and the emergence of irreversibility

In the static approach, see Sec. II A, we have assumed a super-bath to justify the system+bath Gibbs state $\tau SB$ as the equilibrium state. This immediately implied that the MFG state $\tau MF$ is the equilibrium state for the system alone. Now, we abandon the super-bath and adopt the point of view that $SB$ together forms a *closed* system complex. In the following, we will identify system and bath properties which lead to “self-thermalization,” that is, the convergence toward the global Gibbs state $\tau SB$ (return to equilibrium).

The Hamiltonian $Htot$ of the total, closed $SB$ complex determines the dynamics in time *t* from an initial $SB$ state $\rho SB(0)$ according to the Schrödinger equation,

$Htot$ also determines the global Gibbs state $\tau SB$ at inverse temperature *β*, see Eq. (5). We now discuss two mathematical intricacies relating to equilibration toward $\tau SB$, issues that have been kept under the carpet in Sec. II.

There are formal definitions of irreversibility^{156}—here, we mean by it, somewhat intuitively, that averages of suitable observables approach constant values in the limit of large times. The first point to recall is that a *closed* system exhibits truly irreversible dynamics *only* if its Hamiltonian *H* has a continuous energy spectrum.^{157} Indeed, if the energies are discrete, $E1,E2,\u2026\u2208\mathbb{R}$, with corresponding eigenstates $|\psi j\u27e9$, then the evolution is $e\u2212itH=\u2211je\u2212itEj|\psi j\u27e9\u27e8\psi j|$. This shows that the dynamics is simply oscillating for all times.

The connection between irreversible dynamics and continuity of modes (energy spectrum) can be illustrated on a bath consisting of noninteracting particles as follows. Consider a closed system of *n* particles in a region $V\u2282\mathbb{R}3$ with Hamiltonian (kinetic energy),

where $\Delta j=\u2202xj2$ is the Laplacian. For *finite V*, the energies of *H _{V}* are discrete and the dynamics generated by

*H*is quasi-periodic (a sum of oscillating terms). Particles are confined to

_{V}*V*, and boundary effects cause recurrence. However, for $V=\mathbb{R}3$, the spectrum of

*H*becomes continuous (equal to $[0,\u221e)$). The dynamics is now irreversible in the following sense: Given an arbitrary

_{V}*finite*region of observation $R\u2282\mathbb{R}3$, the probability of finding any of the particles inside

*R*converges to zero in the limit of large times (the particles travel to infinity). Of course, even in the infinite volume situation, the probability of finding the particles in all of space $\mathbb{R}3$ equals 100% at all times. This is simply a consequence of the global unitarity of $e\u2212itHV$. However, all observations made in any finite volume (such as a laboratory) show dynamical irreversibility.

^{158}

The second point to highlight is that in writing Eq. (5), one implicitly assumes that the matrix $e\u2212\beta Htot$ is trace-class (meaning that its trace is finite). However, for Hamiltonians $Htot$ which have continuous spectrum, we always have^{159} $trSB[e\u2212\beta Htot]=\u221e$, so the equilibrium state cannot be expressed by Eq. (5). In this situation, one must in fact use a limiting procedure to mathematically define the equilibrium state, as we illustrate in Subsection III C (see also Sec. 4 of Ref. 160 and Ref. 101).

From now, we consider the bath to be a very large and complex environment with a Hamiltonian $HB$ having a continuum of energies. In contrast, the system is “small and simple,” typically having a Hamiltonian $HS$ with discrete (i.e., not continuous) spectrum. The total interacting Hamiltonian $Htot$ for the $SB$ complex, see Eq. (3), then inherits the property of continuous energies. In a sense, if you add to a very complex physical system (the bath $B$) some more degrees of freedom (by coupling it to a small system $S$), then the overall characteristics of the energy spectrum does not change (the continuous spectrum remains a continuous spectrum under such coupling).

Despite our above convention of using the term “bath” in the case of continuous modes (infinite volume), we sometimes want to discuss the situation of large but finite “baths” in which we use the term “finite bath.”

Before discussing the continuum limit of $e\u2212\beta Htot$ in Sec. III C, we first comment on the decription of (apparent) irreversibility for noncontinuous systems.

### B. Finite baths and effective dimension

While any physical lab environment is large—but finite—taking the continuum limit is a meaningful approximation of most real situations, where a multitude of uncontrolled and spatially far extended bath modes may interact with a system of interest. The question in what way a very large, but finite environment can describe thermalization or, more generally, equilibration (convergence to an equilibrium state which is not necessarily thermal), is addressed in Refs. 161 and 162.

The authors consider a generic class of $Htot$ with nondegenerate eigenvalues and nondegenerate Bohr frequencies different from zero. (The Bohr frequencies are differences between the eigenvalues). They show that the average magnitude of fluctuations in time, around the equilibrium state (generally dependent on the initial state), is proportional to $1/deff$, where the effective dimension is defined by

in which $\rho \xafSB$ is the time-averaged state,

where $|Ek\u27e9$ are the eigenvectors of $Htot$, see Refs. 161 and 162. The quantity $trSB[\rho \xafSB2]$ is equal to the time average of the Loschmidt echo, or the survival probability^{163} of the initial state, i.e.,

The recurrence time grows exponentially with $deff$, see Ref. 164. There are estimates for interacting many-body systems,^{163} indicating that $deff$ is exponentially large in the joint system plus bath size. Indeed, strong evidence exists that $deff$ is exponentially large for almost any wavefunction.^{165,166} One may conjecture that $deff$ increases indefinitely with the number of modes; but so far, a rigorous proof of this for the bosonic bath model described above has not been achieved.

As the recurrence time increases with the bath size, taking the infinite volume, or continuous modes limit, corresponds to setting the recurrence time to infinity. This is physically not always realistic, and when it is not, then one has to study the system and bath dynamics for finite baths, which is an intricate task, as it necessitates simultaneously the analysis of the nontrivial (nonconstant) bath dynamics. It has been observed that for certain models, the dynamics converges numerically to a stationary regime on rather fast time scales, even if the bath consists of only a relatively small number of degrees of freedom.^{167–172}

### C. Constructing an infinite volume equilibrium state

The mathematical construction of the equilibrium state $\tau SB$ associated with an infinitely extended system proceeds via two steps: First, one takes the *“thermodynamic limit”* for the bath $B$, resulting in a continuous spectrum of $HB$, and a characterization of its own equilibrium state $\tau B$ at inverse temperature *β*. Second, the much “smaller” system is coupled to the bath, which results in a new system plus bath equilibrium state $\tau SB$.

We now explain the first step for a bath consisting of free bosons with Hamiltonian *H _{V}*, see Eq. (39). The procedure was first carried out in Ref. 173 and further explained in Refs. 160 and 174. Consider a finite volume $V\u2282\mathbb{R}3$ of position space, say a cube of side length

*L*, centered at the origin. The momenta

*k*and eigenstates $|\Psi j\u27e9$ of a single particle are explicitly known. [The single particle Hamiltonian is just the Laplacian $\u2212\Delta $; set $mj=1/2$ in Eq. (39)].

_{j}By applying a suitable selection of creation operators $a\u2020(\Psi j)$ [cf. definition after Eq. (11)] to the vacuum state $|0V\u27e9$ (*V* indicates finite volume), one builds

which is the state describing *n*_{1} particles of momentum *k*_{1} and *n*_{2} particles of momentum *k*_{2}, and so on, in the volume *V*. Now, one increases the volume *V*, keeping the density $\mu j=nj/|V|$ fixed, in such a way that the discrete distribution *μ _{j}* tends to a pre-selected function $\mu (k)$ of continuous momenta $k\u2208\mathbb{R}$. This $\mu (k)$ is called the (continuous) momentum density distribution because $\mu (k)dk$ is the number of particles per unit volume (in position space) having momenta in the volume $dk\u2282\mathbb{R}3$ around

*k*.

As it turns out, the limit cannot be taken directly on the states. Rather, it has to be taken on averages of local observables *A*, which are operators built from (e.g., integrals of) creation and annihilation operators $ax\u2020$, *a _{x}* with $x\u2208V$ for some finite (but arbitrary) $V\u2282\mathbb{R}3$ (the $ax\u2020,ax$ are the Fourier transforms of $ak\u2020,ak$). This limiting procedure defines the average $\u27e8A\u27e9\beta ,\u221e$ of the observable

*A*in the infinite volume state. The values of the expectation functional $A\u21a6\u27e8A\u27e9\beta ,\u221e$, for all observables

*A*, define the infinite volume state. Now, the question is how to represent this state as a vector (or density matrix) $\tau B$. One can find a suitable Hilbert space $H$ and a normalized state $|\Omega \u27e9\u2208H$, such that $\u27e8A\u27e9\beta ,\u221e=\u27e8\Omega |\pi (A)\u2009\Omega \u27e9H=trH[|\Omega \u27e9\u27e8\Omega |\pi (A)]$ (inner product and trace of $H$). Here,

*π*is a

*representation*of the observables, mapping each

*A*to an operator $\pi (A)$ acting on $H$. The vector $|\Omega \u27e9$, or equivalently, the density matrix $|\Omega \u27e9\u27e8\Omega |$, is often times called the

*purification*of the infinite volume state. The triple $(H,\pi ,\Omega )$ is called the

*Gelfand–Naimark–Segal*representation.

^{101,175}It is given explicitly as follows for any prescribed momentum density distribution $\mu (k)$.

^{160,173,176}The Hilbert space is $H=F\u2297F$, where $F$ is the usual Fock space for the bosonic gas in which a general

*N*particle state is given by

where $\Phi (k1,\u2026,kN)$ is the *N*-particle wave function in momentum representation. [Note that in Eq. (44), we integrate over all the possible continuous values $k1,\u2026,kN\u2208\mathbb{R}3$ of the momenta of the *N* particles; $k1,\u2026,kN$ are just integration variables, not to be confused with the fixed momenta chosen to build Eq. (43) in the finite volume situation.] The infinite volume bath state associated with the momentum density distribution $\mu (k)$ is represented as the vector $|\Omega \u27e9:=|0F\u27e9\u2297|0F\u27e9$, where $|0F\u27e9$ is the vacuum state of the Fock space $F$. The representation map is given by

The desired density $\mu (k)$ is correctly reproduced, as one can check easily $\u27e8\Omega |\pi (ak\u2020a\u2113)\u2009\Omega \u27e9H=\mu (k)\delta (k\u2212\u2113)$. The construction works for all momentum density distributions $\mu (k)$. Upon choosing Planck's black body distribution, $\mu (k)=1/(e\beta \omega (k)\u22121)$, the state $|\Omega \u27e9$ is the (purification of the) infinite volume equilibrium state $\tau B$ for the bosonic gas in $\mathbb{R}3$.

We now discuss the second step—introducing the (much smaller) system. For an uncoupled $SB$ complex, with the infinitely extended $B$, the global equilibrium state is $\rho \u2297|\Omega \u27e9\u27e8\Omega |$, where $\rho \u221de\u2212\beta HS$. The full (interacting) $SB$ equilibrium state, corresponding to an $SB$ interaction operator $VSB$, see Eq. (3), is given by^{101,176,177} $\tau SB\u221de\u2212\beta (L0+\lambda \pi (VSB))/2(\rho \u2297|\Omega \u27e9\u27e8\Omega |)e\u2212\beta (L0+\lambda \pi (VSB))/2$. Here, $L0=LS+\u2009LB$ is called the noninteracting Liouville operator with $LS\rho =[HS,\rho ]$ (defined on system density matrices *ρ*) and $LB=HB\u2297\u20091F\u22121F\u2297HB$ (acting on the purification Hilbert space $F\u2297F$), where $HB=\u222b\mathbb{R}3dk\omega (k)ak\u2020ak$.

The construction of the evolution of the infinitely extended bath and system complex follows the same procedure as the above infinite volume limit. Now, one takes the thermodynamic limit of the finite volume Heisenberg picture *evolution* of observables. The evolution is represented in the purified Hilbert space by $e\u2212itL$, where *L* is the Liouville operator. It plays the role of the Hamiltonian, but now this evolution acts in the infinite volume Hilbert space $H$, whose vectors represent states.

To conclude this somewhat technical section, we summarize: It is possible to explicitly construct the equilibrium state of a system–bath complex for a bath that is infinitely, spatially extended. The state is represented by a vector in a new Hilbert space, not simply Fock space $F$. In a way, when taking the volume of the bath to infinity and keeping the density of particles fixed and not zero, the usual Fock space is *not* suitable any longer to describe the equilibrium state. This is due to the fact that any density matrix acting on Fock state describes a state with only finitely many particles, which means a zero density at infinite volume! The explicit form of the infinitely extended $SB$ equilibrium state is an important ingredient in the rigorous analysis of the dynamics, such as for return to equilibrium, see Sec. III D. In the following, we will still use the notation $\tau SB$ given in Eq. (5), even if we mean that the limiting procedure has been performed. We further point out that this construction works for any value of the coupling parameter *λ*; no weak coupling regime is needed here.

### D. Long time $SB$ asymptotics and RtE

We say that the property of RtE holds for a class of $SB$ initial states $S$ and a class of $SB$ observables $O$ if

for all initial states $\rho SB(0)\u2208S$ and all observables $A\u2208O$. In dynamical systems' parlance, Eq. (46) means that the Gibbs state $\tau SB$ is dynamically attractive and has a basin of attraction containing the class of states $S$. The convergence is measured by limits of expectations of $SB$ observables $A\u2208O$. We cannot expect Eq. (46) to hold for all states and all observables. For instance, the initial state $\tau \u2032SB$, the equilibrium at a different temperature $\beta \u2032\u2260\beta $, is stationary, so it does not converge to $\tau SB$ (nor does any other initial state approaching $\tau \u2032SB$ in the long time limit). Also, for models in which the bath is a spatially extended physical system, bath observables which sample space locations arbitrarily far away will capture deviations from the equilibrium state at arbitrarily late moments in time, so $O$ should exclude such global observables.

In the pioneering papers,^{178–184} the property of RtE is shown to hold for an arbitrary *N*-level system coupled to a spatially infinitely extended bath of noninteracting bosons (as explained in Sec. III C). $SB$ coupling is $\lambda VSB$ [see Eq. (3)] with $VSB$ as in Eq. (11). The class of initial states $S$ consists of all states that can be obtained by a local modification of $\tau SB$, and the observable algebra $O$ contains all system observables and all spatially localized bath observables. It is shown that Eq. (46) holds provided the coupling constant *λ* in Eq. (8) is small enough, cf. Eq. (25), namely, $0<|\lambda |<\lambda 0$ for some (not very explicit) *λ*_{0}. As the temperature *T* becomes smaller, the upper bound *λ*_{0} on *λ* shrinks, and the method breaks down in the zero temperature case. Beyond the smallness condition on *λ*, there are two further assumptions: the bath correlation function decays in time (exponential decay is assumed in the above references, while in subsequent improvements,^{185–187} polynomial decay suffices), and the so-called *Fermi Golden Rule Condition* is assumed [$S$ and $B$ are well coupled so that relaxation effects are visible at $O(\lambda 2)$]. Two remarks are in order: (i) The result [Eq. (46)] holds for small enough *λ*, but the final state is *exactly* the coupled equilibrium state $\tau SB$, to *all orders in λ*. (ii) The result [Eq. (46)] is a statement about the full $SB$ dynamics, not merely the reduced system dynamics.

The class of $SB$ observables $O$ contains $OS$, the algebra of all observables acting on the system $S$ alone. For the system dynamics, the primary consequence of the RtE relation (46) is that for all system observables $AS$,

with $\tau MF$ being the MFG state defined in Eq. (6) for the total Hamiltonian (3). A special class of initial states for which (47) holds is that of uncorrelated states of the form

where *ρ* and $\tau B$ are an arbitrary system density matrix and the bath equilibrium state, respectively. In other words, under the conditions mentioned above, namely, that $0<|\lambda |<\lambda 0$, that the bath correlation function decays and that the Fermi Golden Rule holds, Refs. 178, 179, and 183–187 show the following: For any global initial state [Eq. (48)], regardless of the details of *ρ*, the system converges in the long time limit to the mean force Gibbs state $\tau MF$. The same holds for correlated initial states $\rho SB(0)$, which are not of the product form (48), as long as they remain in the class $S$ of states explained at the beginning of Sec. III D.

### E. Relation to non-integrable baths and eigenstate thermalization hypothesis (ETH)

In Subsection II I, we outlined how the bosonic bath can serve as a phenomenological model that can describe more complex, potentially nonintegrable, baths. This raises the issue of comparing the system–bath approaches with those actively studied in many-body physics.^{188–191} The question of convergence to some equilibrium state (and in particular, thermalization to a thermal state) and the identification of the correct form of this state are central questions also in this field.

Many-body physics usually considers systems of many identical particles. Typically, each particle interacts with its neighboring particles, either in physical space (for example, particles in a gas) or in a lattice (for example, spin chains). Such many-body models are often nonintegrable.^{192,193} In this context, the celebrated ETH^{194–197} answers the question about the steady state of a small subsystem as a part of a large isolated system. The ETH can be viewed as a quantum version of the ergodic hypothesis in classical mechanics. Though there is no rigorous proof of the ETH (the same holds for the ergodic hypothesis for most classical models), there is enormous numerical evidence that the ETH is satisfied for a large class of nonintegrable physical models.

A bath consisting of noninteracting particles is integrable and does not obey the ETH. It is not known whether the complex, obtained by coupling this bath to a system, satisfies the ETH. There are studies which show that a localized perturbation of an integrable system is often sufficient to obtain a nonintegrable and thermalizing system,^{198–201} while other results shows that this is not always the case.^{202} Note that the noninteracting bosonic bath coupled to a system, as detailed in Sec. II B, can be exactly mapped into a chain of interacting harmonic oscillators coupled to the system.^{203,204} This observation could serve as a starting point to compare and link the bosonic bath used in open systems' theory with many-body physics baths.

If we use the model of noninteracting bosons as a phenomenological model for a more complex nonintegrable bath satisfying ETH, then one may argue that our question (Q) about the steady state is answered directly by the ETH. However, just like the ergodic hypothesis in classical mechanics, the ETH does not say anything about the *rate* of equilibration. In contrast, with the bosonic bath, one does obtain this more detailed information, including decoherence rates. One may then ask whether the nonintegrability of the bath is responsible for thermalization, or whether thermalization can be explained from the viewpoint of a noninteracting bath model? Imagine a very weakly nonintegrable bath for which one may expect that the thermalization process caused by the ETH occurs very slowly. However, it might be that thermalization occurs much faster due to a mechanism independent of nonintegrability and the ETH, a mechanism which can be explained and understood in the framework of a simplified model of a noninteracting bosonic bath. It would be interesting to compare the predictions of both models and to establish further links between them in future works.

## IV. SYSTEM DYNAMICS AND STEADY STATE

The results presented in Sec. III D are concerned with the asymptotics of the full $SB$ dynamics at $t\u2192\u221e$ leading to Eq. (46), which also implies that the system converges to the stationary state $\tau MF$. Now, we consider the microscopic details of the dynamics of the system alone, starting from the combined system and bath complex evolving according to the unitary dynamics given by the total Hamiltonian $Htot=HS+HB+\lambda \u2009VSB$, as stated in Eq. (3).

### A. General dynamical setting

The dynamical point of view for an *open system* $S$ is generally concerned with describing its state evolution $\rho (t)$ when $S$ is brought into contact with a heat bath $B$ at temperature *T*. One tries to solve the dynamical equations of motion for the system or at least to determine the system's steady state $\rho ss=limt\u2192\u221e\rho (t)$, cf. Eq. (2). One may then address a more refined version of our initial questions (Q):

(Q′) Is the system's dynamical steady state $\rho ss$ equal to the mean force Gibbs state $\tau MF$ discussed in Sec. II?

Below, we summarize some key results on steady states of various dynamical systems and make a connection to $\tau MF$ where possible.

To proceed, consider the general interaction of the form (4). If the operators $X(j)$ commute with $HS$, then the interaction is called *energy conserving*. The system populations are constant in time but the bath still causes irreversible effects in the system, such as decoherence. Those models are suitable for situations in which decoherence happens much more quickly than thermalization and we are interested only in the decoherence time scale.^{68,205} If at least one $X(j)$ does not commute with $HS$, then energy exchange processes between the system and bath are enabled.

The initial $SB$ state, $\rho SB(0)$, is usually assumed to be of product form (48), although the evolution of correlated or entangled initial states is also relevant and studied.^{187,206–212} The state of $S$ at time *t* is the reduced density matrix, cf. Eq. (38),

Taking the time-derivative gives

from which we want to derive an autonomous equation, a *master equation*, for $\rho (t)$. For product initial states (48), Eq. (50) can be cast into the form of an integro-differential equation for $\rho (t)$, which is called the *Nakajima–Zwanzig master equation.*^{59} For an arbitrary initial state $\rho SB(0)$, where $S$ and $B$ are correlated or entangled, the master equation for $\rho (t)$ will depend on the initial $SB$ correlation, see e.g., Refs. 103, 187, 206, 208, 210, 211, and 213–216.

With very few exceptions, a master equation for $\rho (t)$ (in an analytically closed form) cannot be derived without recourse to approximations. In Secs. IV C–IV E, we will discuss various approximations in the well-studied weak coupling limit as well as the ultrastrong and intermediate limit. Before we do so, we will introduce some useful terms and notations and then begin the discussion with the exactly solvable model of the open dynamics of the quantum harmonic oscillator.

Unless otherwise stated, the analysis of the master equations discussed in Secs. IV C–IV E assumes the Hamiltonian $Htot$ defined in Eq. (3) with a single interaction term $X\u2297B$ in Eq. (4). We will make use of decomposition (23) of the operator *X* into a sum of energy eigenoperators *X _{m}*. In the interaction picture (denoted by the tilde), one has $X\u0303m(t)=Xm\u2009ei\omega mt$. These quantities define two time scales—one associated with the Bohr frequencies

*ω*and the other associated with the differences of Bohr frequencies, $\omega mn=\omega m\u2212\omega n$.

_{m}Furthermore, the *bath correlation function G*(*t*) plays a central role in open system dynamics, particularly at weak coupling, as it determines the memory time of the $SB$ interaction. It is defined as the autocorrelation function of the bath operator *B*,

and is the well-known dynamic version of the imaginary-time bath correlation function introduced in Sec. II G. It satisfies the Kubo–Martin–Schwinger (KMS) condition $G(t)=G(\u2212t\u2212i\beta )$, where *β* is the inverse temperature of the bath Gibbs state $\tau B$. *G*(*t*) is assumed to vanish as $t\u2192\u221e$, and this decay defines the time scale, $tB$, of the correlations of the quantum noise. We also define the time dependent coefficients,

which will appear in the master equations below.

### B. Exact dynamics of the damped quantum harmonic oscillator

A general method to solving^{96} the exact dynamics of a damped harmonic oscillator described by Eq. (10) follows a path integral functional integral approach. For general spectral density $J(\omega )$ and coupling strength *λ*, and a broad class of initial states of the global system, including the global Gibbs state $\tau SB$ and entangled states, it is shown^{96} that their steady states are all the same. The initial state $\tau SB$ is clearly a stationary state of the global evolution, and its reduced system state is $\tau MF$. Hence, $\tau MF$ is the steady state for the whole class of initial states considered. This state is given in Eq. (17) in Sec. II.

An alternative method to show the dynamical convergence to $\tau MF$ is obtained in Ref. 217 using Heisenberg–Langevin equation methods. More recent work^{100} addresses general *N*-body quantum Brownian motion, for which the general two-time correlation functions are explicitly evaluated.^{100} This is done for initial states $\rho SB(0)=\rho (0)\u2297\tau B$, cf. Eq. (48), and in the steady state limit of long times $t\u2192\u221e$. Then, the above stationary state argument^{96} for $\tau SB$ is again applied to prove that the steady state of the *N*-body system, for initial states $\rho SB(0)=\rho (0)\u2297\tau B$, must coincide with the MFG state $\tau MF$.

Working from the functional integral description of the system dynamics, it is also possible to construct an exact non-Markovian master equation^{78} and specify its steady state in terms of its Wigner function.^{218} A generalization to the case of the driven damped quantum oscillator has been also obtained.^{219}

Taken together, the studies above firmly establish that the dynamical steady state $\rho ss$ of a quantum oscillator, under dynamics given by Eq. (10) for a broad class of initial global states, is *exactly* the mean force Gibbs state $\tau MF$ at all system-bath coupling strengths *λ* and for general spectral density $J(\omega )$.

### C. Weak coupling dynamics

In most cases, approximate forms for the master equation are obtained on the assumption that the system-bath coupling is weak enough that a second order perturbative treatment suffices. This defines the “weak coupling limit” for master equations, as indicated in Fig. 2. Weak coupling applies in a variety of contexts, including quantum optical systems,^{148–152} nuclear magnetic resonance,^{220–229} solid state and molecular physics,^{153,155,230,231} and in certain biological systems.^{115,232–237} Weak coupling expansions generally lead to master equations that are relatively simple and have useful immediate physical interpretations. For example, they predict basic properties of open quantum systems such as decoherence and thermalization^{6,59,238–242} as well as more complicated properties such as laws of thermodynamics,^{243–248} environment-assisted quantum transport,^{249–253} super-radiance and supertransfer,^{254,255} emergence of dark states,^{256–258} and decoherence-free subspaces.^{259} Even if a weak-coupling approximation is not strictly valid, many important system properties can be qualitatively captured in this approximation. Often, phenomenological equations of the GKSL form^{260} are used,^{152,233} which implicitly assume weak coupling. Moreover, for the example of light-harvesting complexes, the system (representing the electronic DoFs) is often strongly coupled only to a finite number of distinguished vibrational modes.^{261} If these are included into an enlarged system, then the weak coupling approximation can be applied to the enlarged system now interacting with the remainder of the bath.^{237,262}

The weak coupling limit allows for the introduction of two important approximations. The first is the *Born approximation*, which is based on the assumption that, as far as the system dynamics is concerned, any changes in the state of the bath, or any correlations that might develop between the system and bath due to their interaction,^{59,187,214,215} can be neglected. The second approximation is the *Markov approximation* in which the rate of change of the state of the system at time *t* is assumed to depend only on its state at that same time, and not on its history. The Born approximation is justified due to the difference of the sizes of $S$ and $B$ (small influence of $S$ on $B$), and the Markov approximation is justified via a separation-of-time scale argument.^{59} These approximations are almost universally inserted into derivations of master equations, for example, for the Bloch–Redfield ME discussed in Sec. IV C 2.

#### 1. Weak coupling: Davies theory and resonance theory

A rigorous analysis of the weak coupling limit has been provided in pioneering work by Davies.^{140,263,264} The master equation is derived in the *Bogoliubov–van Hove limit*,^{265,266} where $\theta =\lambda 2t$ is taken as constant, but $\lambda \u21920$ and simultaneously $t\u2192\u221e$. The rescaled time *θ* is sometimes called the coarse-grained time. For initial product states $\rho \u2297\tau B$, and under certain conditions on the bath correlation function *G*(*t*), a master equation can then be derived rigorously without making the explicit Born–Markov approximation. Specifically, Davies shows that the reduced system dynamics $\rho \u0303(\theta )$ (interaction picture indicated by $\u0303$) is well approximated by some $\rho \u0303D(\theta )$, i.e.,

where the convergence is uniform on the segment $\theta \u2208[0,\Theta ]$ for an arbitrary finite Θ and $||\xb7||$ denotes the trace norm. Here, $\rho \u0303D(\theta )$ is the solution to the *Davies master equation* in the interaction picture, and for the coarse-grained time $\theta =\lambda 2t$,^{140,263,264}

where $L\u0303D$ does not depend on *θ* and is given by^{267}

The positive damping rates are given by $\gamma m=2\u2009Re[\Gamma m]$, and the term that renormalizes the system energy is

where *X _{m}* and $\Gamma m$ are the energy eigenoperators and the master equation coefficients defined in Eqs. (23) and (52), respectively. Note that this renormalization commutes with the system's bare Hamiltonian $HS$, i.e., $[HS,\Delta H||]=0$.

We now discuss the implications of the convergence (53) at coarse-grained times $\theta \u22650$. If one assumes non-degenerate energy levels of $HS$ and non-degenerate non-zero frequency differences *ω _{mn}*, see Sec. IV A, then Eq. (55) predicts no coupling between system populations (diagonal entries of $\rho \u0303$ in the $HS$ eigenbasis) and the coherences (off-diagonals). The coherences are found to decay to zero as $t\u2192\u221e$.

The above assumptions imply that *γ _{m}* satisfy the detailed balance conditions $\gamma m=e\u2212\beta \omega m\gamma \u2212m$. As a consequence, the steady state of Eq. (54) is then readily shown

^{140,263,264}to be the system Gibbs state

*τ*, cf. Eq. (1). For degenerate energy levels, the steady state is nonunique

^{240,242,255,268–270}and the situation becomes more complex.

Thus, to summarize, working in the Bogoliubov-van Hove limit and assuming the generic nondegenerate case, Davies showed that the steady state is the standard Gibbs state *τ*, confirming the validity of Gibbs statistical physics, see Fig. 2. However, this jars with the static assumption of the steady state being the MFG state $\tau MF$, which includes $O(\lambda 2)$ and higher order corrections in comparison to *τ*, see Eq. (20). There are a number of reasons for this incongruence.

First, in order to obtain an equation in the so-called *GKSL form*,^{6,260,271–274} Davies made the so-called secular approximation (see Sec. IV C 3) in his derivation of Eq. (54). But this removes those dynamical features that could lead to a final state $\tau MF$, and instead forces convergence to *τ*. Second, the Davies master equation is of second order in *λ*, but relaxation also takes place on time scales of order $\lambda \u22122$. Small errors accumulate over long times, and thus, the steady state is trustworthy only up to the zeroth order. Third, the convergence in Eq. (54) is proved for arbitrary large but finite rescaled time intervals $\theta \u2208[0,\Theta ]$. Davies' proof requires *λ* to be smaller and smaller as Θ increases. Hence, strictly speaking, it is not guaranteed that the steady state of the Davies master equation coincides with the true steady state in the weak-coupling limit.

These limitations are overcome by the *quantum resonance theory*, which improves Davies' results. Namely, under the same conditions used to show RtE (Sec. III D), it is shown in Refs. 176, 185–187, and 275 that Eq. (53) can be improved to

valid for all times $t\u22650$ (no coarse graining necessary), where *C* is a constant independent of *λ* and *t*. Equation (57) shows that the solution of the Davies master equation approximates the system dynamics to $O(\lambda 2)$*at all time scales*. In particular, this guarantees that the system Gibbs state *τ*, which is the steady state predicted by the Davies master equation, is the true steady state up to $O(\lambda 2)$. This is entirely consistent with the result of RtE, saying that the stationary state is $\tau MF=\tau +O(\lambda 2)$. It is further shown in Refs. 176 and 275 that by adding to the Davies generator $L\u0303D$ higher order terms,^{276}

the solution $\rho \u0303KM(t)$ of the corresponding master equation $ddt\rho \u0303KM(t)=L\u0303KM\rho KM(t)$ has the following properties. First, it is asymptotically exact, meaning that the stationary state is the exact mean force Gibbs state, $limt\u2192\u221e\rho \u0303KM(t)=\tau MF$. Second, the populations (diagonal density matrix elements in the $HS$ basis) of the system are approximated by those of $\rho \u0303KM(t)$ to $O(\lambda )$ for all times $t\u22650$. The coherences of the system (off-diagonals) are guaranteed to be approximated by those of $\rho \u0303KM(t)$ to $O(\lambda )$ for times *t* in the windows $\lambda 2t<C1$ and $\lambda 2t>C2$ for some constants *C*_{1}, *C*_{2}. The corrections $L\u0303k$ are constructed by an explicit perturbation procedure.

We close this section by mentioning that in the literature, perturbation theory is usually carried out by simply neglecting higher order terms without controlling their size for large times. In contrast, the resonance theory is a rigorous approach in which remainder terms are estimated to be small for all times.

#### 2. Weak coupling: Bloch–Redfield master equation (BRME)

For $Htot$ with a single interaction term $\lambda \u2009X\u2297B$, the general form of the second order ($O(\lambda 2)$) master equation obtained by invoking the Born and Markov approximations is known as the *Bloch–Redfield master equation*,^{59,220–222,277}

The energy renormalization contributions consist of a contribution that commutes with $HS$,

and one that generally does not commute with $HS$,

where *X _{m}* and $\Gamma m(t)$ are the energy eigenoperators and the master equation's damping coefficients defined in Eqs. (23) and (52), respectively.

The second line in Eq. (59) contains the damping matrix

The time dependence of the generator $LtBR$ implies that the master equation is non-Markovian in the sense that the time evolution operator Λ_{t}, defined by $\rho (t)=\Lambda t\rho (0)$, is such that $\Lambda t\Lambda s\u2260\Lambda s+t$, i.e., it does not satisfy the semigroup property criterion for Markovianity.^{278} This time dependence makes solving Eq. (59) much more difficult than solving the Davies ME [Eq. (54)] for which $L\u0303D$ is *θ*-independent. However, due to the bath correlation decay, the functions $\Gamma m(t)$ saturate to constant values for times $t\u226btB$. In this regime, one can replace $\Gamma m(t)$ by $\Gamma m(\u221e)\u2261\Gamma m$ and obtain the time independent generator $L\u221eBR$. At large times, the generators $LtBR$ and $L\u221eBR$ coincide, and since we are interested in steady-state solutions, we will always mean the simpler, time-independent version $L\u221eBR$.

In Sec. IV C 1, we found that the steady state of the Davies master equation is the Gibbs state *τ*, i.e., this steady state does not show any signature of the bath interaction other than the bath's inverse temperature *β*. Having next introduced the Bloch–Redfield equation, we are now in the position to discuss the complex world of its steady state(s).

A first setback is that, contrary to the Davies ME, there is no analytical expression for the steady state of the Bloch–Redfield equation. Two approaches concerning the issue of determining the steady state can be identified, which are discussed below. The first is the one adopted by Geva *et.al.*^{102} and in more detail by Mori and Miyashita^{103} in which the aim is to derive an expression for the time derivative $d\rho /dt$ and confirm to second order in *λ* that the mean force Gibbs state is such that this derivative vanishes, i.e., $\tau MF$ “qualifies” as the dynamical steady state. The second approach adopted, e.g., by Fleming *et al.,*^{279} Thingna *et al.,*^{104} Subaşi *et al.*,^{100} and Purkayastha *et al*^{66} is to construct a master equation for $\rho (t)$ to search for its long time steady state solution to second order in *λ* and to compare this with the mean force Gibbs state also evaluated to second order.

As an example for the first approach, Mori and Miyashita^{103} consider the time derivative $\rho \u0307(t)$ to second order in *λ* (weak coupling), which is given in terms of $\rho (t)$ and the total system state $\rho SB(t0)$ at an earlier time *t*_{0}. Like the quantum oscillator case discussed in Sec. IV B, this initial state is chosen as either (i) $\rho (t0)\u2297\tau B$ or (ii) the global Gibbs state $\tau SB$. Assuming now that $\rho (t)$ is set equal to the mean force Gibbs state $\tau MF$ [also to second order in *λ* see Eq. (20)], Mori and Miyashita confirm^{103} that the derivative $\rho \u0307(t)$ is zero (to second order in *λ*) for both initial state choices (i) and (ii), i.e., $L\u221eBR(\tau MF)=0=\rho \u0307$. Interestingly, the same statement is recovered by letting $t0\u2192\u2212\u221e$ while *no constraint* is imposed on the form of the initial state in that infinite past, i.e., all memory of the initial conditions is lost. It is quite general then to view the above conclusion to hold *independently of the initial state of the total system*, at least to second order in *λ*.

The above result^{103} is applicable for a general quantum systems *S* in contact with a bath at inverse temperature *β* and is valid to second order in *λ*. It confirms that $\tau MF$ at inverse temperature *β* is a steady state of the Bloch–Redfield dynamics $L\u221eBR$ given in (59).

Mori and Miyashita also note that there is, however, an important ambiguity regarding the steady state as follows. The algebraic equation $L\u221eBR(\rho )=O(\lambda 3)$ for the steady state yields a steady state $\rho ss(2)$ accurate to $O(\lambda 2)$. However, for any second order generator $L(2)$, and for *W* a traceless operator that is diagonal in the $HS$ basis, one has $L(0)(W)=\u2212i\u2009[HS,W]=0$ and $L(2)(W)=O(\lambda 2)$. Thus, for states $\rho ss(2)+\lambda 2W$, one finds that they are also steady state solutions to the same order in *λ* as $\rho ss(2)$ itself, i.e., $L\u221eBR(\rho ss(2)+\lambda 2W)=O(\lambda 3)$ due to the linearity of the evolution $L\u221eBR$ in the density matrix *ρ*, cf. Eq. (59). This ambiguity implies that while the steady state $\rho ss(2)$ is correct to second order in $\lambda 2$ for the *off-diagonal elements* in the $HS$ basis, the *diagonal elements* are correct only to zeroth order that is the order where it is correct to simply replace the diagonal elements of $\tau MF$ with those of *τ*, see Eq. (20).

Fleming *et al*^{279} (see also Ref. 280) show that the above ambiguity in the steady state is a consequence of the perturbative approach used to construct the master equation: solving a 2*n* order perturbative master equation will yield the diagonal elements of the steady state accurate only to order $2n\u22122$ (while the off-diagonal elements are determined to order 2*n*). Thus, resolving the second order indeterminacy requires expanding the master equation to *fourth order*, and then using degenerate perturbation theory to find the needed corrections, thus fixing *W* in the above. A slightly different proof of the last conclusion is given in Ref. 280.

Thingna *et al*^{104} reach a similar conclusion to Ref. 279, but bypass the need to calculate fourth order terms in the master equation, typically a very complex task. They do so by making use of analytic continuation methods to derive the corrections to the diagonal elements of the steady state based on its second order off-diagonal elements alone. They confirm that the resulting steady state solution, now correct to second order in *λ* and uniquely defined, is identical to the mean force Gibbs state $\tau MF$ to $O(\lambda 2)$. Subaşi *et al.*^{100} generalize this result to the case where the system–bath interaction is of the more general form (4), while the initial state is restricted to the product state (48). For a double-quantum-dot charge qubit, described by the spin-boson model,^{66} the above results^{100,103,104,279} on the convergence of the dynamics (to second order in *λ*) to the $\tau MF$ (to second order in *λ*) are calculated and illustrated in detail in Ref. 66.

The Bloch–Redfield master equation does have the problem of generating non-positive probabilities, at least for some early times. The trustworthiness of this master equation has been questioned^{60} for this reason, though inclusion of a slippage of initial conditions^{281–283} corrects this without having an impact on the final steady state.^{284,354} However, what we saw here is that its long time steady state, $\rho ss$, has corrections to *τ* that are partially consistent with $\tau MF$, in line with the expectation from the static point of-view discussed in Sec. II. Assuming that the weak coupling limit is justified [here meaning expansions in $\lambda 2$ as well as $\lambda 4$ to get the $O(\lambda 2)$ corrections for the diagonals], it, hence, appears that the BRME is trustworthy when it comes to predicting long time behavior and equilibration. It also gives well-behaved, positive predictions at shorter times when the initial system state is assumed to be close to $\tau MF$.

It is worthwhile to comment on the above results^{100,103} in comparison to those of RtE, cf. Sec. III D. Coincidence of the *perturbative expansion* of the steady state of the BRME with that of the mean force Gibbs state up to a finite order does not, in fact, guarantee that these states are exactly the same. To contrast, the results on RtE prove the exact (i.e., in all orders of *λ*) equality of these states. This is proven for any *λ* smaller than a fixed strictly positive number *λ*_{0}. For this reason, the results of Refs. 100 and 103 are most relevant in the context of discussions about the choice between the Bloch–Redfield and the Davies (secular Bloch–Redfield) master equation: Which one reproduces the true equilibrium state more precisely. If, for a concrete problem to solve, it is important to have a steady state that contains corrections to the Gibbs state *τ*, then the Bloch–Redfield equation or the master equation (58) should be preferred.

#### 3. Weak coupling: Further results

A consequence of taking the long time form $L\u221eBR$ for the Bloch–Redfield generator is that the master equation is not of the GKSL form^{6,271–274} as the damping matrix, *γ _{mn}*, defined in Eq. (62), is not necessarily positive definite. Negative probabilities can now occur for some initial states of the system

^{6,285,286}at times less than $tB$, the bath correlation time. This weakness can be removed by making a further approximation, the secular approximation, already introduced by Davies, cf. Sec. IV C 1, which involves removing those terms $Xm\rho Xn\u2020$ in Eq. (59) that oscillate, in the interaction picture, at frequencies $\omega mn\u226btS\u22121$, where $tS$ is the system relaxation time. The resultant full-secular-approximation form of the Bloch–Redfield master equation can be readily shown

^{59}to be of the GKSL form and, moreover, has the consequence that the steady state is the Gibbs state,

*τ*. However, if nearly degenerate Bohr frequencies are found to occur, then the corresponding non-secular terms in the master equation are important and should not be removed.

It has been noted that the full secular approximation is quite often made indiscriminately, even in circumstances where it is not justified. This appears motivated, at least in part, by it leading to a GKSL master equation that has the Gibbs state as its steady state, thus connecting the ME dynamics to Gibbs state physics, see Fig. 2. If only those non-secular terms that satisfy $\omega mn\u226btS\u22121$ are removed,^{113,287–291} called partial secular approximation, the resultant master equation will still not be of the GKSL form, and moreover, the steady state is no longer guaranteed to be either the Gibbs state or a second order approximation to the mean force Gibbs state. A further approximation, often made in practice, is to ignore the imaginary parts of the master equation coefficients $\Gamma m$, see Eq. (52). It turns out^{103,292} that the Bloch–Redfield equation with $\Gamma m$ replaced by $Re[\Gamma m]$ has the Gibbs state *τ* as its steady state.

Recently, a rigorous derivation of a *unified* GKSL quantum master equation beyond the secular approximation was proposed.^{293} A rigorous procedure leads, in general, to a partial secular approximation followed by certain modifications in the arguments of the spectral density function. Since the derivation is rigorous, the resulting equation is thermodynamically consistent, see also Refs. 294 and 295. In particular, a Gibbs state with respect to a modified system Hamiltonian (in which all nearly degenerate energy levels and Bohr frequencies become exactly degenerate) is its steady state.

### D. Ultrastrong coupling limit and strong decoherence limit

Beyond master equations based on weak coupling expansions, it is possible to develop^{296} a dynamical perturbation theory that extends to the ultrastrong coupling regime, see Sec. II F. Consider again the total Hamiltonian (21). For a non-degenerate *X* with spectral decomposition into $Pn=|xn\u27e9\u27e8xn|$ given in Eq. (26), the trick is to decompose $HS$ in the pointer basis ${|xn\u27e9}$ from the outset,

Here, $\epsilon m$ and Δ_{mn} are real and complex numbers, respectively. One may now assume that Δ_{mn} can be treated as small compared to a very large system–bath coupling strength *λ*. More precisely, one should compare Δ_{mn} with another quantity of the dimensionality of energy, e.g., $\lambda 2J(\omega )$ for a characteristic Bohr frequency *ω* may serve as such quantity. Both large *λ* and small Δ_{mn} imply that the decoherence in the “pointer basis” ${|xm\u27e9}$ takes place on a time scale much smaller than the relaxation of the populations, $pm=\u27e8xm|\rho |xm\u27e9$, see Ref. 205 for a detailed analysis of the spin-boson model. Using this observation, a Pauli master equation for the populations in the pointer basis $pm=pm(t)$ is obtained,^{296}

where the rate constants *k _{mn}* follow detailed balance,

and are proportional to $|\Delta mn|2$ and decrease exponentially^{296} with $\lambda 2$. This is a generalization of the Förster approximation, which is well-known in the theory of excitation energy transfer,^{69} see Sec. II F. The detailed balance condition implies that steady-state populations *p _{m}* are proportional to $e\u2212\beta \epsilon m$. The off-diagonal elements in the pointer basis, $\u27e8xm|\rho (t)|xn\u27e9$, tend to zero for an arbitrary

*t*>

*0 as $\Delta mn\u21920$ (or $\lambda \u2192\u221e$). Thus, the dynamical steady state at $\lambda \u2192\u221e$ is exactly the static MFG state $\tau MF$ given in Eq. (27). The proof presented in Ref. 296 extends to multiple interactions $X(j)$ in Eq. (4), degenerate spectra of the $X(j)$, as well as combinations of weak- and strong-coupling parts of the interaction Hamiltonian. All these cases assume strong decoherence between certain subspaces, thus giving the name “strong decoherence limit” for the corresponding perturbation theory. The steady state at high, but not infinitely large,*

*λ*has now also been derived

^{296}and is found to agree well with the steady state that is numerically solved in Fig. 1 at strong coupling, see also Subsection IV F. Also, numerical calculations of Ref. 125 show excellent correspondence between the MFG state and the steady state at high, but not infinitely large

*λ*.

Pioneering work^{118,124} on the ultrastrong coupling limit has previously conjectured a dynamical steady state based on arguments from einselection theory.^{122} They argued that the state should be diagonal in the pointer basis *P _{m}*, which is now confirmed analytically.

^{296}The steady state diagonals in the pointer basis were conjectured to be $\u27e8xm|\tau |xm\u27e9$ provided that the initial state is the Gibbs state

*τ*. This conjecture differs from the expressions derived above. Moreover, such a state is clearly non-stationary in the case of non-zero Δ

_{mn}. Nevertheless, since Δ

_{mn}are much smaller than $\lambda 2J(\omega )$, decoherence in the pointer basis occurs much faster than the transfer between the different pointer states.

^{296}Hence, the conjectured state is, in fact, rapidly achieved due to strong decoherence. However, thereafter, the system still slowly relaxes to the MFG state $\tau MF$. This resolves the discussion

^{82,118}about the correct form of the steady state at ultrastrong coupling.

Finally, compared to the mathematical rigorous convergence proof in Sec. III D for small enough coupling, the proof outlined here for the ultrastrong coupling is only “physically rigorous,” i.e., based on a microscopic model and physically plausible assumptions. For the spin-boson model, and under certain conditions on the spectral density, the convergence has also been proven mathematically rigorously using a combination of the polaron transformation and quantum resonance theory,^{71,297,298} see Subsection IV E 2.

### E. Intermediate coupling

So far, we have been discussing two limiting cases: weak and ultrastrong coupling limits, see Fig. 2(b). However, the case of intermediate coupling is highly interesting. Since in this case, generally there is no basis for a perturbation theory, non-perturbative techniques need to be developed and employed. We will briefly discuss two known analytical methods that can explore this regime: the reaction coordinate approach and the polaron transformation.

#### 1. Reaction coordinate approach

The reaction coordinate approach is a non-perturbative approach to dealing with intermediate and strongly coupled systems. It was originally developed quite some time ago^{302,303} but has recently been much extended to various applications in the quantum thermodynamic context.^{20,25,29,30,304–315} The basic idea is to redefine that part of the bath $B$ that contributes most strongly to the system–bath coupling into a collective single degree of freedom, the reaction coordinate (RC). This coordinate is then incorporated along with the original system into an enlarged system $S\u2009&RC$. The remnant degrees of freedom of the original bath $B$ now constitute a reduced bath $B\u0303$. The whole point is that the enlarged system $S\u2009&RC$ may interact only weakly with the reduced bath $B\u0303$ with a new small coupling parameter $\lambda \u0303$. Thus, the impact of $B\u0303$ on $S\u2009&RC$ can be treated by the perturbative approaches discussed in the weak coupling, see Sec. IV C. A development of the reaction coordinate method leads to a mapping of a free bosonic bath into a chain of interacting harmonic oscillators coupled to a system at one end.^{203,204} This mapping was mentioned in Sec. III E in relation to many-body physics. It can be used for numerical simulations: the so-called time evolving density matrix using orthogonal polynomials (TEDOPA) algorithm.

The reaction coordinate approach can offer insight into the steady state of the original system when the coupling to the original bath is in the intermediate coupling limit, see Fig. 2. First, one obtains an expression for the steady state of the enlarged system $S\u2009&RC$, which is usually just the Gibbs state $\tau S&RC$ for the corresponding Hamiltonian of $S&RC$. By tracing out the reaction coordinate, the steady state of the original system $S$ is then be found. In Ref. 305, numerical results are presented that strongly indicate that this state is not the Gibbs state *τ*. In Refs. 20 and 307, the steady state of $S$ is shown to be $\tau MF+O(\lambda \u0303)$, where $\tau MF$ is the mean force Gibbs state for $S$, as defined in Eq. (6). In Ref. 316, steady states obtained by the reaction coordinate method are compared with the perturbative result (22) for the MFG state in the weak coupling regime.

#### 2. Polaron transformation

In Sec. II G, we mentioned the use of the polaron transformation for the derivation of approximations to the mean force Gibbs state. A quantum master equation in the polaron-transformed frame was derived in Refs. 138 and 317–320 Since the polaron transformation mixes the system and bath DoFs, originally, the method allowed to evaluate only the populations (in the system eigenbasis). Formulas for the coherences were derived in Ref. 261. As in the static considerations, if the off-diagonal elements of the system Hamiltonian in the pointer basis are small enough with respect to the system–bath coupling, the polaron-transformed Bloch–Redfield master equation correctly describes the dynamics. Again, the variational (instead of full) polaron transformation can be applied; the corresponding master equation and its application to excitation energy transfer and excitonic Rabi rotations in a driven quantum dot were considered in Refs. 321–323.

The steady states of the (full) polaron-transformed master equation were studied in Ref. 324, see also the review article.^{137} They confirm convergence to the mean force Gibbs state.

Combining the polaron transformation and quantum resonance theory, it is shown rigorously in Ref. 297 that the spin-boson model discussed in Subsection II G exhibits RtE (see Sec. III) for arbitrary values of the coupling strength *λ*. The authors show that under certain conditions on the spectral density [in particular, $J(\omega )$ should be proportional to *ω ^{s}* for $s\u22653$ for small

*ω*], if the coupling Δ in the system Hamiltonian (28) is smaller than a certain value $\Delta 0$ (depending on the other parameters of the model), then the coupled $SB$ dynamics converges to the joint $SB$ equilibrium state $\tau SB$ asymptotically in time. In particular, the reduced state of the system converges to the mean force Gibbs state, which rigorously proves the result of Sec. IV D for this particular model. In the subsequent works,

^{71,298}the reduced dynamics of the population and coherences of the two-level system is derived for all times.

The analysis of Ref. 71 is carried out in a more general setting, including lower temperatures, collective and local (independent donor and acceptor) baths, and donors and acceptors coupled with different strengths to the bath(s). It leads to a generalized Marcus formula^{325} with the original one (high temperature, common bath, and same donor and acceptor coupling strength) as a special case.

We stress, however, that, even in the case of proper (super-Ohmic) spectral densities, the polaron transformation followed by weak coupling perturbation theory in the polaron frame cannot serve as a universal tool. It has a certain range of validity, and establishing precise conditions on where the weak coupling perturbation theory in the polaron frame can truly be applied remains largely an open problem.

#### 3. Low density limit

A common assumption in the above discussion is the choice of interaction (11), which is linear in creation and annihilation operators of the bath. However, open quantum systems theory also studies other types of interactions, bath models, and limits, such as the singular coupling limit and the low density limit.^{59,60} While the singular coupling limit is shown to be equivalent to a special kind of the weak coupling limit,^{326,327} which is captured in the unified weak coupling framework,^{293} the low density limits^{328–330} represent significantly different theory.

Now, the interaction $VSB$ is quadratic in the creation and annihilation operators, which represents a collision of a system with a particle of the bath. The state of the bath is a grand canonical (instead of canonical) equilibrium state. The low density limit is the limit $n\u21920$, where *n* is the density of particles of the bath. Thus, the collisions of the system with the particles of the bath are not weak, but rare. A close relation of this model to the collision model (or repeated interactions model), popular in the theory of open quantum systems and quantum thermodynamics,^{331–339} is established in Ref. 340.

In the low density limit, one can derive a Markovian quantum master equation of the GKSL form for the reduced system state $\rho (t)$.^{328–330} The steady state(s) of this equation are much less studied than those for the weak coupling MEs discussed above. Some sufficient conditions for the Gibbs state *τ* to be a unique steady state were derived in Refs. 328 and 341. However, in a recent paper,^{342} it is shown that, under certain (non-generic) conditions on the spectrum of $HS$, the system state converges to a unique steady state that is different from *τ*. It is an open question whether this steady state is the mean force Gibbs state.

### F. Steady states with non-perturbative numerical methods

There is a range of powerful numerical methods able to solve the dynamics of an open quantum system for different coupling limits as well as the intermediate coupling regime. These include the hierarchical equations of motion (HEOM)^{299,300} and TEMPO, which is based on time-evolving matrix product operators.^{301} We will not describe these numerical methods here. We used HEOM [in its high-temperature version (Ref. 62)] to calculate the dynamics of the qubit system shown in Fig. 1. HEOM gives numerically exact results. In both the weak and strong coupling cases, we find excellent agreement of the steady state $\rho ss$ with the high temperature MFG state formula (37) from Ref. 363. This is due to the fact that for the parameters given in Fig. 1, the temperature is sufficiently high that even in the strong coupling case, the condition $\u2113\beta <2$ is still satisfied. For the (moderately) strong coupling case in Fig. 1, we find very good agreement of the steady state $\rho ss$ also with the MFG state formula (27) for the ultrastrong coupling limit ($\lambda \u2192\u221e$) derived in Ref. 82. Inclusion of corrections to this formula for large, but finite, *λ*, obtained in Refs. 125 and 296 gives the most precise match with the numerical steady state.

## V. CONCLUSIONS AND OPEN QUESTIONS

We now conclude on the findings outlined above and provide an extensive list of open questions.

Section II summarized the static point of view. This view includes non-negligible coupling into a modified equilibrium state, the MFG state $\tau MF$, and provides the backbone for much current research on constructing a thermodynamic framework that goes beyond Gibbs state physics.^{8,11–17,19,21,23,24,26–31,33–35} Despite the growing importance of this framework for the assessment of thermodynamic processes of nanoscale and quantum systems, e.g., in terms of heat and work exchanges as well as entropy production on the level of fluctuations and averages, explicit expressions for the concrete functional shape of $\tau MF$, are only known in a handful of cases.

Beyond the exactly solvable quantum oscillator, explicit expressions for general systems are known for weak [neglecting $O(\lambda 4)$ and higher] and ultrastrong [neglecting $O(\lambda \u22121)$ and higher] coupling limits. It remains an open question **(a)** as to whether or not there exist relatively simple, tractable expressions for the MFG state for intermediate coupling regimes, at least for specific system choices. This is the interesting regime where, loosely speaking, the system's interaction with its environment is of comparable scale to the system's bare energy. For example, it may be possible to construct useful expansions for *λ* that neither expand for small or large *λ*, but some intermediate value *λ*_{0}.^{343,362} However, we note that previous considerable analytical effort has not resulted in “simple” expressions even for classical systems. It also remains open **(b)** whether there exists a general criterion, or criteria, that allows one to judge, just looking at the Hamiltonian of a particular problem, whether mean force corrections are going to be important or not in the equilibrium state. For the weak coupling limit, such criteria exist [see, e.g., Eq. (25)], which reveal a complex interplay of system energies, system–bath coupling parameters, and bath temperature. Analogous condition(s) in other coupling regimes remain to be uncovered.

Another open question is the extension of the above results, all valid for coupling to a single (continuum) bosonic bath, to **(c)** simultaneous coupling to multiple bosonic baths and to **(d)** fermionic baths. **(e)** Nonlinear couplings in the bath ladder operators could also be considered.

Furthermore, connections between quantum and classical MFG states are well worth exploring. First, they would **(f)** provide a direct connection^{8} to numerical simulations in (classical) chemistry.^{50,52–54} These routinely include potential of mean force corrections in the calculation of arrangements of molecules in solutions in a vast range of contexts. A timely example are simulations of the catalytic mechanism of proteases in respiratory syndromes.^{55} Second, **(g)** a comparison between quantum and classical cases could also illuminate quantum signatures in the mean force Gibbs state, which are not present in the classical counterpart.^{344} Another intriguing question **(h)** is whether one can derive, or not, the MFG state $\tau MF$ in some manner from entropy-maximization arguments in line with such derivations of the standard Gibbs state as outlined in the beginning of Sec. II. For both the classical or quantum cases, the difficulty here is that the strong-coupling corrected system energy and entropy functionals can depend on temperature.^{8,19,31,36}

Section III outlined how equilibrium arises dynamically, culminating in the mathematical proof that a $SB$ complex does equilibrate exactly to the *global* Gibbs state $\tau SB$, when the initial state is not too far away from it, and the coupling *λ* is finite and small. Under those assumptions, this result immediately dynamically justifies the postulated $\tau MF$ used in the statics section. A major open question is to **(i)** give a detailed quantification of the upper bound on *λ* in terms of physical quantities and **(j)** to prove the global convergence to $\tau SB$ and/or the local convergence^{345} to $\tau MF$ for a larger (or the whole) range of coupling strengths *λ*. In regard to this, all bounds obtained so far shrink with temperature *T*. It remains an open problem to **(k)** show RtE at low (or zero) temperatures^{346} $T\u226a\lambda 2$.

Another issue is dimension: in all proofs of RtE, the system $S$ is assumed to be finite-dimensional. Giving **(l)** a proof of RtE for infinite-dimensional systems (with discrete spectrum), such as a harmonic oscillator, is an important problem to solve. Likewise for an environment with a discrete spectrum, **(m)** estimating the effective dimension $deff$ in Eq. (40) in terms of the volume and/or the number of bath modes is an open problem.

Finally, current studies on quantum thermodynamics often make use of the properties of CPTP maps and define the system's dissipated heat *Q* (with implications for work, entropy, efficiency, etc.) as the energy received by an environment^{347} during a global unitary operation on an $SB$ complex. However, if the environment is small, such as a qubit, then the very notion of it being a thermodynamic context is being abandoned. If the environment is somewhat larger, such as a discrete “bath” of harmonic oscillators, then this is somewhat reasonable given the results discussed in Sec. III B. Yet, it is worth keeping in mind that even in this case, the dynamics does not produce the irreversible character, see Sec. III D, often assumed in thermodynamic arguments. A major task is **(n)** the careful assessment of how the results outlined in Sec. III come to bear on such thermodynamic characterizations.

Section IV outlined a selection of microscopic approaches, in particular master equations, which are capable of giving analytical details of a system's evolution as well as an indication of its steady state. While the application of the secular approximation in the weak coupling MEs of Davies and others predict *τ* as steady state, the Bloch-Redfield ME predicts $\lambda 2$-corrections to *τ*. However, here the issue was that not all $\lambda 2$-corrections are in fact captured by the $\lambda 2$-order BRME: some corrections require higher order expansions.

Most ME treatments assume linear coupling to bath modes. However, in the low density limit, quadratic (collision-like) couplings^{342} are important and give rise **(o)** to open problems concerning their steady states and consistency with $\tau MF$. In the intermediate coupling limit, reaction coordinates or polaron transformation methods can be used, but they cannot serve as universal tools as they have their own limitations. Finally, the ultrastrong coupling regime stands out because here the ME discussed in Sec. IV D leads, without any ambiguities, to the corresponding mean force Gibbs state (27) as detailed in Sec. II F (which differs from the standard Gibbs state *τ*).

To conclude, in this Review, we provided an introduction to several current avenues in the disparate fields of open quantum systems, strong coupling thermodynamics, and beyond. The aim of these fields is to uncover the bath's signature on a nanoscale system's equilibrium state as well as elucidate the system's approach to equilibrium. Impressive results have been achieved addressing this timely challenge—but many key questions remain open. Solving these will provide much needed clarity on how to consistently characterize the thermodynamics of nanoscale and quantum systems.

While some researchers may feel it is obvious that the MFG state $\tau MF$ should be the equilibrium state, we highlight that much current research in quantum thermodynamics still tacitly assumes it to be the Gibbs state *τ*. This includes many master equation derivations as well as the theory of thermal operations and thermodynamic resource theory to name a few. Causes for this scientific mismatch might be the unwieldy formal definition of $\tau MF$, partially resolved in Sec. II, and the conflict between unitary evolution and irreversibility, partially resolved in Sec. III. The results outlined above provide a glimpse of a theory that goes beyond Gibbs statistical physics and will find applications in a variety of fields, where the exchange of energy on the nanoscale is essential, from quantum chemistry and biology, to magnetism and nanoscale heat management.

## ACKNOWLEDGMENTS

We thank Ahsan Nazir, Chris Jarzynski, Dvira Segal, and Jonathan Keeling for valuable comments on a draft of the manuscript. We are indebted to Philipp Strasberg for extensive and thoughtful comments and for bringing up the Gaussian limit argument that underpins much of the success of the Caldeira–Leggett model. M.M. thanks Gennady Berman, Jürg Fröhlich, Alain Joye, Martin Könenberg, and Michael Sigal for guidance and collaborations over many years. A.S.T. is grateful to Alexander Pechen, Alexander Teretenkov, Camille Lombard Latune, and Oleg Lychkovskiy for valuable comments, bibliography suggestions, and fruitful ongoing discussions. J.D.C. and J.A. thank Simon Horsley, Marco Berritta, Stefano Scali, and Federico Cerisola for many stimulating discussions on the subject of this Review. M.M. is supported by a *Discovery Grant* from the *Natural Sciences and Engineering Research Council of Canada* (NSERC). J.A. and J.D.C. acknowledge funding from the Engineering and Physical Sciences Research Council (EPSRC) (No. EP/R045577/1). J.A. further acknowledges EPSRC support in the form of a Doctoral Training Grant and thanks the Royal Society for support.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The dynamical evolution data displayed in Fig. 1 are available from author A. S. Trushechkin upon reasonable request.

## References

^{59}

^{348}and developed in Refs. 349–352.

^{271–273}is $\rho \u0307=\u2212i[H\u2032,\rho ]+\u2211n[Ln\rho Ln\u2020\u221212{Ln\u2020Ln,\rho}],$ where $H\u2032$ and

*L*are system operators and $H\u2032$ is self-adjoint. Often, this form is called the Lindblad form, but the GKSL form is more correct, see Ref. 274.

_{n}^{260}