Open quantum system dynamics and the mean force Gibbs state

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⁵ 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⁸ or a quantum spin with the phononic modes within a material,⁹ 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 $ρ(t)$ be the system density matrix at time t and denote the system steady state as

We ask

• (Q) Is the system steady state $ρss$ the Gibbs state τ?

•      If not, how does $ρ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.

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 $ρSB(0)$⁠. One then asks whether the reduced state of the system alone, $ρ(t):=trB[ρSB(t)]$⁠, stops changing (significantly) at late times t. If this is so, then one may define (one or more) system steady states $ρss$⁠. It is often assumed that the initial $SB$ state is uncorrelated, $ρSB(0)=ρ(0)⊗τB$⁠, where $τB:=e−HB/kBT/trB[e−HB/kBT]$ is the bath Gibbs state (bath Hamiltonian $HB$ and temperature T). Furthermore, it is common to make the so-called Born approximation, $ρSB(t)≈ρ(t)⊗τB$ for all times $t≥0$⁠, 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 $τSB∝e−Htot/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, $τMF=trB[τ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 $ρ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 $τ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) $ρss$ and the reduced global equilibrium state $τ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.

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⊂ℝ3$ of position space), continuous bath energies arise in the limit of infinite volume (⁠$R→ℝ3$⁠), see also Sec. III A.

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

where $λ∈ℝ$ 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 $ℏ=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 $λ−1$ may be performed, and all orders $λ−1$ 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.

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.⁴⁰ Gibbs' original derivation⁴¹ (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 $⟨U⟩$ 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.⁴⁷ 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 $⟨U⟩$⁠, becomes $β=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 $⟨U⟩$⁠.

We now consider the system and bath complex, $SB$⁠, to be in the global Gibbs state $τ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 $β=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−β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, $τMF$ differs—sometimes substantially—from the Gibbs state $τ∝e−βHS$⁠, as we will see below. The naming arises from casting $τ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 $τ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 $τMF$ in Eq. (7). This is due to the fact that such a constant would also change the partition function $ZMF=trS[e−β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”¹¹ 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,³¹ and quantum measurements have been included in a stochastic description of strongly coupled quantum systems.³³ 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.³⁵ 

We now return to the MFG state $τMF$⁠, which is uniquely defined by the formal identity [Eq. (6)]. However, unfortunately, giving explicit expressions for $τ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.

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†$ and annihilation operators a_k obey the commutation relations $[ak,aℓ†]=δk,ℓ$⁠, and X is an arbitrary system operator. g_k 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 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 $λ2∑k|gk|2X2/ωk$ is often included in the Hamiltonian, which physically arises whenever coupling is introduced via the difference of coordinates, for example, $(xk−X)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., $∑kω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.⁷³ Here, the inclusion of the counterterm⁷⁴ 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)−λ2∑k|gk|2X2/ω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ω02x2/2$⁠, giving the Hamiltonian⁷⁷ 

where x_k and p_k are the bath oscillator position and momentum operators, respectively, $gk=−ck ℏ/(m mkωk)$⁠, and the coupling operator is $X=x m/2$⁠. The Caldeira–Leggett model is a key model for open quantum systems theory^21,78,79 with widespread application to quantum tunnelling⁸⁰ and studies of decoherence and the quantum-classical transition.⁷⁸ 

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 ${ωk}k$ by a continuum of frequencies, for instance, $ω∈[0,∞)$⁠. 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→∞$⁠. 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).⁸³ The continuous mode limit then leads to $k∈ℝd$⁠, and the continuous mode Hamiltonian associated with Eq. (8), becomes^9,81

where $a†(g)=∫ℝddk gk ak†$ 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∈ℝ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ℓ†]=δ(k−ℓ)$⁠, 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∝∫0∞dω ω aω†aω$ and $VSB∝X⊗∫0∞dω gω aω†+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 $∫dk$⁠.

If $J(ω)∝ω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(ω)/ω$ as $ω→0$ 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(ω)coth(βω/2)∝J(ω)/ω$ 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.⁹⁴ 

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, $ω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⁹⁷ of the oscillator position and momentum operators: $⟨x⟩, ⟨p⟩, ⟨x2⟩, ⟨p2⟩$⁠, and $⟨px⟩$⁠. The first moments trivially vanish, while the second moments are given in terms of the partition function $Z(ω0,γ)$ at inverse temperature β,

where $Γ(z)$ denotes the gamma function and $ν=2π/β$ is the first Matsubara frequency. The functions $μj(ω0,γ)$⁠, for j = 1, 2, 3, denote the roots of the cubic polynomial

With these expressions, the second moments given in unit-free form [$x̃$ and $p̃$⁠, see Eq. (10)], and including $ℏ$ explicitly, are⁹⁵ 

Reference 95 implies the MFG state as the Gibbs state of an effective harmonic oscillator Hamiltonian (the HMF) $HSeff=p2/(2m(λ))+m(λ)ω02(λ)x2/2$⁠. They formally establish the mass $m(λ)$ and frequency $ω0(λ)$ 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,⁹⁷ 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(ω)$⁠, e.g.,

where the mass was set to m = 1 and $ℏ=1$ again. Here, $G(ω)=−1/(ω2−ω02[1−χ(ω)])$ is a Green's function, and $χ(ω)$ is a complex susceptibility given by

where the integral is understood as the principal part integral. Setting $t=t′$ 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²¹ of damped quantum and classical harmonic oscillators as a functional of $χ(ω)$ and inverse temperature β.

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 $τMF$ can be considered. Such expansions are based on the Kubo identity⁹⁹ and are sometimes referred to as “canonical perturbation theory,”¹⁰⁰ equivalent to a standard time dependent perturbation expansion with t replaced by $−iβ$⁠. In the mathematical literature, the expansion in λ is known as the perturbation theory of Kubo–Martin–Schwinger (KMS) states.¹⁰¹ 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 $τ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(ω)$ is used immediately without specifying any coupling constants c_k. Here, $[qω,pω′]=i δ(ω−ω′)$ 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 $τ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_m are defined by

with ω_m being the Bohr frequencies of the system (energy differences of $HS$⁠). Furthermore, the function $Dβ(ω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 $τMF$ (see Fig. 1 and Sec. IV F). Coherences are often considered a quantum “resource,”¹⁰⁵ 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,⁸² 

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 λ.

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 $λ→∞$⁠, see Fig. 2 and Refs. 116–120. Here, it is also possible⁸² to find an explicit expression for $τ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_n are orthogonal projectors of rank one. By expanding the global Gibbs state $τSB$ to orders of $1/λ$ and explicitly integrating out the bath oscillators, the mean force Gibbs state simplifies to⁸² 

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]≠0$⁠, an immediate consequence is that $τMF$ will maintain coherences in the $HS$ energy basis $|em⟩$⁠, i.e., $⟨em|τMF|em⟩≠0$ for some m (see Fig. 1 and Sec. IV F). Corrections to Eq. (27) with respect to $λ−1$ 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.¹²⁷ In the ultrastrong coupling limit, which corresponds to the Förster theory of excitation energy transfer,⁶⁹ 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.

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).¹²⁹ 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=σz$⁠, where $σ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̃tot=UHtotU†=H̃S+HB+ṼSB$⁠, where $H̃S=ε2σz+κ Δ2σx$ with

The bath part remains unchanged, $HB=∫0∞ dω ω aω†aω$⁠, and the interaction becomes

where λ has now moved inside $ṼSB$⁠. When the integral in the definition of κ converges, it turns out that $trB[τB e±i2λR̂]=κ$⁠. 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(ω)$ is proportional to $ω3$⁠, but it diverges whenever $J(ω)$ is proportional to ω^s for $s≤2$⁠. This is a restriction of the polaron transformation method.

The factor κ represents the above-mentioned “phonon cloud,” while the B_x and B_y operators represent fluctuations around this cloud. One may hope that these fluctuations are not large and $ṼSB$ 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 $λ2J(ω)$], where the polaron transformation is trivially reduced to the identity transformation, and (ii) the limit of small Δ (weak tunneling limit).^137,138 Since $|κΔ|<|Δ|$⁠, the eigenvectors of $H̃S$ 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 $τ̃SB=U τSB U†∝e−βH̃tot$⁠. One can show that the diagonal part (in the pointer basis) of the desired MFG state $τMF$ formally coincides with the diagonal part of the reduced system state $η̃=trB[τ̃SB]$⁠, i.e., $τMF,jj=η̃jj$ for j = 1, 2. Approximate expressions for $η̃$ were obtained^135–137 using second-order perturbation theory with respect to $ṼSB$⁠,

where $η̃(0)=e−βH̃S/Z̃$ with $Z̃=trS[e−βH̃S]$ is the Gibbs state corresponding to the system Hamiltonian in the polaron frame. The next term in Eq. (32) is

where

and $Gmn(β)=trSB[Bm(β)BnτB]$ are the imaginary-time bath correlation functions. The operators in imaginary time, such as $σm(β)$ and $Bm(β)$⁠, are defined as $O(β)=eβ(H̃S+HB)Oe−β(H̃S+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 $τMF,12$ (in pointer basis), $τ̃MF$ does not suffice—the bath DoFs of the total polaron-transformed Gibbs state are also required. Up to the first order in $ṼSB$⁠, it is found^136,137 to be

where $Λ=ε2+(κΔ)2, Sm(β′)=trS[σm(β′) σ− τ̃]$⁠, and $Km(β′)=trB[Vm(β′) cos (2λR̂) τB]$⁠, where the functions S_m and K_m can be explicitly evaluated and $σ−=(σx−iσy)/2$⁠.

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

where γ (or, more precisely, $λ2γ$⁠) determines the system–bath coupling strength, while $ω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 $τMF$ are compared with numerically exact simulations. It turns out that the approximation works well for the cases of fast bath $ωc>Δ$ and ultrastrong coupling (large λ). It remains an open question to simplify expressions (32) for $τ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.

An approximate expression for $τ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⟩}$ interacts, via $Xn=|n⟩⟨n|$⁠, with several baths $n=1,…,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ε+HJ$⁠, where $Hε$ is the part that is diagonal in the $|n⟩$ 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, $Λ=λ2∑n=1N∫0∞dω Jn(ω)ω Xn≡∑n=1Nℓn Xn$ is an operator-valued reorganization term. The result [Eq. (37)] implies that the interstate coupling constants J_mn contained in H_J, describing hopping between states $|m⟩$ and $|n⟩$⁠, get rescaled by the bath interaction into effective coupling constants. The rescaling itself is temperature dependent and vanishes for $β→0$⁠. For a dimer system with $ℓ1=ℓ2=ℓ$⁠, expression (37) is found to be accurate (Ref. 363) for temperatures satisfying $ℓβ≲2$⁠. 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).

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.¹⁵³ 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).⁸⁰ 

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.

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.

In the static approach, see Sec. II A, we have assumed a super-bath to justify the system+bath Gibbs state $τSB$ as the equilibrium state. This immediately implied that the MFG state $τ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 $τSB$ (return to equilibrium).

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

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

There are formal definitions of irreversibility¹⁵⁶—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.¹⁵⁷ Indeed, if the energies are discrete, $E1,E2,…∈ℝ$⁠, with corresponding eigenstates $|ψj⟩$⁠, then the evolution is $e−itH=∑je−itEj|ψj⟩⟨ψ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⊂ℝ3$ with Hamiltonian (kinetic energy),

where $Δj=∂xj2$ is the Laplacian. For finite V, the energies of H_V are discrete and the dynamics generated by H_V is quasi-periodic (a sum of oscillating terms). Particles are confined to V, and boundary effects cause recurrence. However, for $V=ℝ3$⁠, the spectrum of H_V becomes continuous (equal to $[0,∞)$⁠). The dynamics is now irreversible in the following sense: Given an arbitrary finite region of observation $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 $ℝ3$ equals 100% at all times. This is simply a consequence of the global unitarity of $e−itHV$⁠. However, all observations made in any finite volume (such as a laboratory) show dynamical irreversibility.¹⁵⁸ 

The second point to highlight is that in writing Eq. (5), one implicitly assumes that the matrix $e−βHtot$ is trace-class (meaning that its trace is finite). However, for Hamiltonians $Htot$ which have continuous spectrum, we always have¹⁵⁹ $trSB[e−βHtot]=∞$⁠, 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−βHtot$ in Sec. III C, we first comment on the decription of (apparent) irreversibility for noncontinuous systems.

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 $ρ¯SB$ is the time-averaged state,

where $|Ek⟩$ are the eigenvectors of $Htot$⁠, see Refs. 161 and 162. The quantity $trSB[ρ¯SB2]$ is equal to the time average of the Loschmidt echo, or the survival probability¹⁶³ of the initial state, i.e.,

The recurrence time grows exponentially with $deff$⁠, see Ref. 164. There are estimates for interacting many-body systems,¹⁶³ 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 

The mathematical construction of the equilibrium state $τ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 $τB$ at inverse temperature β. Second, the much “smaller” system is coupled to the bath, which results in a new system plus bath equilibrium state $τ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⊂ℝ3$ of position space, say a cube of side length L, centered at the origin. The momenta k_j and eigenstates $|Ψj⟩$ of a single particle are explicitly known. [The single particle Hamiltonian is just the Laplacian $−Δ$⁠; set $mj=1/2$ in Eq. (39)].

By applying a suitable selection of creation operators $a†(Ψj)$ [cf. definition after Eq. (11)] to the vacuum state $|0V⟩$ (V indicates finite volume), one builds

which is the state describing n₁ particles of momentum k₁ and n₂ particles of momentum k₂, and so on, in the volume V. Now, one increases the volume V, keeping the density $μj=nj/|V|$ fixed, in such a way that the discrete distribution μ_j tends to a pre-selected function $μ(k)$ of continuous momenta $k∈ℝ$⁠. This $μ(k)$ is called the (continuous) momentum density distribution because $μ(k)dk$ is the number of particles per unit volume (in position space) having momenta in the volume $dk⊂ℝ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†$⁠, a_x with $x∈V$ for some finite (but arbitrary) $V⊂ℝ3$ (the $ax†,ax$ are the Fourier transforms of $ak†,ak$⁠). This limiting procedure defines the average $⟨A⟩β,∞$ of the observable A in the infinite volume state. The values of the expectation functional $A↦⟨A⟩β,∞$⁠, for all observables A, define the infinite volume state. Now, the question is how to represent this state as a vector (or density matrix) $τB$⁠. One can find a suitable Hilbert space $H$ and a normalized state $|Ω⟩∈H$⁠, such that $⟨A⟩β,∞=⟨Ω|π(A) Ω⟩H=trH[|Ω⟩⟨Ω|π(A)]$ (inner product and trace of $H$⁠). Here, π is a representation of the observables, mapping each A to an operator $π(A)$ acting on $H$⁠. The vector $|Ω⟩$⁠, or equivalently, the density matrix $|Ω⟩⟨Ω|$⁠, is often times called the purification of the infinite volume state. The triple $(H,π,Ω)$ is called the Gelfand–Naimark–Segal representation.^101,175 It is given explicitly as follows for any prescribed momentum density distribution $μ(k)$⁠.^160,173,176 The Hilbert space is $H=F⊗F$⁠, where $F$ is the usual Fock space for the bosonic gas in which a general N particle state is given by

where $Φ(k1,…,kN)$ is the N-particle wave function in momentum representation. [Note that in Eq. (44), we integrate over all the possible continuous values $k1,…,kN∈ℝ3$ of the momenta of the N particles; $k1,…,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 $μ(k)$ is represented as the vector $|Ω⟩:=|0F⟩⊗|0F⟩$⁠, where $|0F⟩$ is the vacuum state of the Fock space $F$⁠. The representation map is given by

The desired density $μ(k)$ is correctly reproduced, as one can check easily $⟨Ω|π(ak†aℓ) Ω⟩H=μ(k)δ(k−ℓ)$⁠. The construction works for all momentum density distributions $μ(k)$⁠. Upon choosing Planck's black body distribution, $μ(k)=1/(eβω(k)−1)$⁠, the state $|Ω⟩$ is the (purification of the) infinite volume equilibrium state $τB$ for the bosonic gas in $ℝ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 $ρ⊗|Ω⟩⟨Ω|$⁠, where $ρ∝e−βHS$⁠. The full (interacting) $SB$ equilibrium state, corresponding to an $SB$ interaction operator $VSB$⁠, see Eq. (3), is given by^101,176,177 $τSB∝e−β(L0+λπ(VSB))/2(ρ⊗|Ω⟩⟨Ω|)e−β(L0+λπ(VSB))/2$⁠. Here, $L0=LS+ LB$ is called the noninteracting Liouville operator with $LSρ=[HS,ρ]$ (defined on system density matrices ρ) and $LB=HB⊗ 1F−1F⊗HB$ (acting on the purification Hilbert space $F⊗F$⁠), where $HB=∫ℝ3dkω(k)ak†ak$⁠.