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.
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 , 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 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 theory1–4 and is assumed in “thermal operations,” which investigate the properties of CPTP maps5 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 solution8 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 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 be the system density matrix at time t and denote the system steady state as
(Q) Is the system steady state the Gibbs state τ?
If not, how does 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 . One then asks whether the reduced state of the system alone, , stops changing (significantly) at late times t. If this is so, then one may define (one or more) system steady states . It is often assumed that the initial state is uncorrelated, , where is the bath Gibbs state (bath Hamiltonian and temperature T). Furthermore, it is common to make the so-called Born approximation, for all times , 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 and .
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 , where 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, , 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 . 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 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) and the reduced global equilibrium state 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 , consisting of the system and the remaining part, called the bath. Deciding which part of an interacting complex is and which is 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 and are not equal partners: the bath influences the thermodynamic and dynamical properties of the system significantly, while cannot move too far from its initial state. This is modeled by taking bath Hamiltonians , which have a continuum of energies, while system Hamiltonians 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 of position space), continuous bath energies arise in the limit of infinite volume (), see also Sec. III A.
The total Hamiltonian of a general complex has the following form:
where is a dimensionless coupling constant and is the system-bath interaction operator. The latter is generally of the form
where and are operators acting on the system and bath Hilbert spaces, respectively.
Note that throughout the text, we will subindex states for system+bath with , e.g., the global Gibbs state is τSB, and we will subindex states for the bath alone with , but for the system states, we will drop the index , e.g., the system Gibbs state is denoted as τ. We will, however, keep the index S for the system Hamiltonian . Furthermore, unless otherwise stated, we set .
One distinguishes several regimes related to the magnitude of the coupling strength λ. Usually, these are set by a comparison of typical system () energy differences and energies associated with the interaction . 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 may be performed, and all orders 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.40 Gibbs' original derivation41 (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 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 , 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 , becomes . 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 .
A. Mean force Gibbs state
We now consider the system and bath complex, , to be in the global Gibbs state associated with the total Hamiltonian [Eq. (3)]. The emergence of this state can be justified—for now—by considering that the compound has been in very weak thermal contact for a long time with a super-bath14,48 R at inverse temperature . Gibbs state statistical physics for the compound then tells us that the equilibrium state for is the Gibbs state
where 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, differs—sometimes substantially—from the Gibbs state , as we will see below. The naming arises from casting in the Gibbsian (exponential) form
for an effective Hamiltonian 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 chemistry49–55 since the 1930s.
Before proceeding with the discussion of the state , let us first comment on the HMF. Unlike the bare system Hamiltonian 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 in Eq. (3). Note that the is not uniquely defined since, e.g., a constant may be added to it without changing in Eq. (7). This is due to the fact that such a constant would also change the partition function . (The constant also cancels when calculating energetic differences). A common choice is to set with 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 equality13,14,58 and the Clausius inequality15–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 , which is uniquely defined by the formal identity [Eq. (6)]. However, unfortunately, giving explicit expressions for 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 , 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
where ωk are the frequencies of the oscillator (modes) k, the creation and annihilation operators ak obey the commutation relations , and X is an arbitrary system operator. gk are complex numbers that weigh the strength of the interaction between and the oscillator mode at frequency ωk. Even though , Eq. (8), has infinitely many energy levels, those levels do not fill a continuum of values. Thus, technically, according to our bath definition, 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 is often included in the Hamiltonian, which physically arises whenever coupling is introduced via the difference of coordinates, for example, , instead of a product, e.g., . 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., , 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 computers7,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 counterterm74 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 . 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 , giving the Hamiltonian77
where xk and pk are the bath oscillator position and momentum operators, respectively, , and the coupling operator is . The Caldeira–Leggett model is a key model for open quantum systems theory21,78,79 with widespread application to quantum tunnelling80 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 by a continuum of frequencies, for instance, . 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 model9,81,82,173,176,275 and then analyze the full statics and dynamics. This allows one, in particular, to control the perturbation theory for all times, even . 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 , and the continuous mode Hamiltonian associated with Eq. (8), becomes9,81
where is the creation operator smoothed out with gk, which in this context is sometimes called the form factor, a square-integrable complex function of . 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, , 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 and
For a discrete model, an alternative to specifying the coupling constants gk 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 .
If 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 as 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 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 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 ck). The most explicit results can be obtained in the case of the Drude–Lorentz spectral density,
in the limit of continuous modes. Here, 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 moments97 of the oscillator position and momentum operators: , and . The first moments trivially vanish, while the second moments are given in terms of the partition function at inverse temperature β,
where denotes the gamma function and is the first Matsubara frequency. The functions , for j = 1, 2, 3, denote the roots of the cubic polynomial
Reference 95 implies the MFG state as the Gibbs state of an effective harmonic oscillator Hamiltonian (the HMF) . They formally establish the mass and frequency 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 derived21,98 using Heisenberg equations of motion for an arbitrary spectral density , e.g.,
where the mass was set to m = 1 and again. Here, is a Green's function, and is a complex susceptibility given by
where the integral is understood as the principal part integral. Setting 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 energies21 of damped quantum and classical harmonic oscillators as a functional of and inverse temperature β.
E. Expansion of 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 can be considered. Such expansions are based on the Kubo identity99 and are sometimes referred to as “canonical perturbation theory,”100 equivalent to a standard time dependent perturbation expansion with t replaced by . 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 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 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 with system Hamiltonian 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 is used immediately without specifying any coupling constants ck. Here, 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 for an arbitrary system was found to be
Here, the decomposition of the Hermitian system operator X into a sum of energy eigenoperators Xm is used, where Xm are defined by
with ωm being the Bohr frequencies of the system (energy differences of ). Furthermore, the function is defined as
where the integral is understood as a principal part integral. Expression (22) evidences the appearance of coherences in the basis in the system's equilibrium state (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 , see Fig. 2 and Refs. 116–120. Here, it is also possible82 to find an explicit expression for 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 xn are real numbers and Pn are orthogonal projectors of rank one. By expanding the global Gibbs state to orders of and explicitly integrating out the bath oscillators, the mean force Gibbs state simplifies to82
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 , an immediate consequence is that will maintain coherences in the energy basis , i.e., for some m (see Fig. 1 and Sec. IV F). Corrections to Eq. (27) with respect to 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 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 , 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.
and the system coupling operator is , where 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 , where with
The bath part remains unchanged, , and the interaction becomes
where λ has now moved inside . When the integral in the definition of κ converges, it turns out that . 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, is proportional to , but it diverges whenever is proportional to ωs for . This is a restriction of the polaron transformation method.
The factor κ represents the above-mentioned “phonon cloud,” while the Bx and By operators represent fluctuations around this cloud. One may hope that these fluctuations are not large and 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 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 ], 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 are more “localized” superpositions of the pointer basis vectors than the eigenvectors of . 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 . One can show that the diagonal part (in the pointer basis) of the desired MFG state formally coincides with the diagonal part of the reduced system state , i.e., for j = 1, 2. Approximate expressions for were obtained135–137 using second-order perturbation theory with respect to ,
where with is the Gibbs state corresponding to the system Hamiltonian in the polaron frame. The next term in Eq. (32) is
and are the imaginary-time bath correlation functions. The operators in imaginary time, such as and , are defined as . 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 (in pointer basis), does not suffice—the bath DoFs of the total polaron-transformed Gibbs state are also required. Up to the first order in , it is found136,137 to be
where , and , where the functions Sm and Km can be explicitly evaluated and .
In Refs. 135–137, the following (super-Ohmic) spectral density is considered:
where γ (or, more precisely, ) determines the system–bath coupling strength, while is the cutoff frequency and determines the rate of relaxation of the bath correlation functions in time. The above expressions for the elements of are compared with numerically exact simulations. It turns out that the approximation works well for the cases of fast bath and ultrastrong coupling (large λ). It remains an open question to simplify expressions (32) for 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 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 interacts, via , with several baths that all have the same temperature. Instead of (21), the Hamiltonian is
for independent baths,
The system Hamiltonian can be decomposed as , where is the part that is diagonal in the basis and HJ 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, is an operator-valued reorganization term. The result [Eq. (37)] implies that the interstate coupling constants Jmn contained in HJ, describing hopping between states and , get rescaled by the bath interaction into effective coupling constants. The rescaling itself is temperature dependent and vanishes for . For a dimer system with , expression (37) is found to be accurate (Ref. 363) for temperatures satisfying . 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).
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 used153,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 as the equilibrium state. This immediately implied that the MFG state is the equilibrium state for the system alone. Now, we abandon the super-bath and adopt the point of view that 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 (return to equilibrium).
The Hamiltonian of the total, closed complex determines the dynamics in time t from an initial state according to the Schrödinger equation,
also determines the global Gibbs state at inverse temperature β, see Eq. (5). We now discuss two mathematical intricacies relating to equilibration toward , issues that have been kept under the carpet in Sec. II.
There are formal definitions of irreversibility156—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, , with corresponding eigenstates , then the evolution is . 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 with Hamiltonian (kinetic energy),
where is the Laplacian. For finite V, the energies of HV are discrete and the dynamics generated by HV is quasi-periodic (a sum of oscillating terms). Particles are confined to V, and boundary effects cause recurrence. However, for , the spectrum of HV becomes continuous (equal to ). The dynamics is now irreversible in the following sense: Given an arbitrary finite region of observation , 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 equals 100% at all times. This is simply a consequence of the global unitarity of . 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 is trace-class (meaning that its trace is finite). However, for Hamiltonians which have continuous spectrum, we always have159 , 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 having a continuum of energies. In contrast, the system is “small and simple,” typically having a Hamiltonian with discrete (i.e., not continuous) spectrum. The total interacting Hamiltonian for the 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 ) some more degrees of freedom (by coupling it to a small system ), 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 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 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 , where the effective dimension is defined by
in which is the time-averaged state,
The recurrence time grows exponentially with , see Ref. 164. There are estimates for interacting many-body systems,163 indicating that is exponentially large in the joint system plus bath size. Indeed, strong evidence exists that is exponentially large for almost any wavefunction.165,166 One may conjecture that 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 associated with an infinitely extended system proceeds via two steps: First, one takes the “thermodynamic limit” for the bath , resulting in a continuous spectrum of , and a characterization of its own equilibrium state at inverse temperature β. Second, the much “smaller” system is coupled to the bath, which results in a new system plus bath equilibrium state .
We now explain the first step for a bath consisting of free bosons with Hamiltonian HV, see Eq. (39). The procedure was first carried out in Ref. 173 and further explained in Refs. 160 and 174. Consider a finite volume of position space, say a cube of side length L, centered at the origin. The momenta kj and eigenstates of a single particle are explicitly known. [The single particle Hamiltonian is just the Laplacian ; set in Eq. (39)].
By applying a suitable selection of creation operators [cf. definition after Eq. (11)] to the vacuum state (V indicates finite volume), one builds
which is the state describing n1 particles of momentum k1 and n2 particles of momentum k2, and so on, in the volume V. Now, one increases the volume V, keeping the density fixed, in such a way that the discrete distribution μj tends to a pre-selected function of continuous momenta . This is called the (continuous) momentum density distribution because is the number of particles per unit volume (in position space) having momenta in the volume 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 with for some finite (but arbitrary) (the are the Fourier transforms of ). This limiting procedure defines the average of the observable A in the infinite volume state. The values of the expectation functional , for all observables A, define the infinite volume state. Now, the question is how to represent this state as a vector (or density matrix) . One can find a suitable Hilbert space and a normalized state , such that (inner product and trace of ). Here, π is a representation of the observables, mapping each A to an operator acting on . The vector , or equivalently, the density matrix , is often times called the purification of the infinite volume state. The triple is called the Gelfand–Naimark–Segal representation.101,175 It is given explicitly as follows for any prescribed momentum density distribution .160,173,176 The Hilbert space is , where is the usual Fock space for the bosonic gas in which a general N particle state is given by
where is the N-particle wave function in momentum representation. [Note that in Eq. (44), we integrate over all the possible continuous values of the momenta of the N particles; 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 is represented as the vector , where is the vacuum state of the Fock space . The representation map is given by
The desired density is correctly reproduced, as one can check easily . The construction works for all momentum density distributions . Upon choosing Planck's black body distribution, , the state is the (purification of the) infinite volume equilibrium state for the bosonic gas in .
We now discuss the second step—introducing the (much smaller) system. For an uncoupled complex, with the infinitely extended , the global equilibrium state is , where . The full (interacting) equilibrium state, corresponding to an interaction operator , see Eq. (3), is given by101,176,177 . Here, is called the noninteracting Liouville operator with (defined on system density matrices ρ) and (acting on the purification Hilbert space ), where .
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 , where L is the Liouville operator. It plays the role of the Hamiltonian, but now this evolution acts in the infinite volume Hilbert space , 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 . 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 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 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 asymptotics and RtE
We say that the property of RtE holds for a class of initial states and a class of observables if
for all initial states and all observables . In dynamical systems' parlance, Eq. (46) means that the Gibbs state is dynamically attractive and has a basin of attraction containing the class of states . The convergence is measured by limits of expectations of observables . We cannot expect Eq. (46) to hold for all states and all observables. For instance, the initial state , the equilibrium at a different temperature , is stationary, so it does not converge to (nor does any other initial state approaching 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 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). coupling is [see Eq. (3)] with as in Eq. (11). The class of initial states consists of all states that can be obtained by a local modification of , and the observable algebra 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, 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 [ and are well coupled so that relaxation effects are visible at ]. Two remarks are in order: (i) The result [Eq. (46)] holds for small enough λ, but the final state is exactly the coupled equilibrium state , to all orders in λ. (ii) The result [Eq. (46)] is a statement about the full dynamics, not merely the reduced system dynamics.
The class of observables contains , the algebra of all observables acting on the system alone. For the system dynamics, the primary consequence of the RtE relation (46) is that for all system observables ,
where ρ and are an arbitrary system density matrix and the bath equilibrium state, respectively. In other words, under the conditions mentioned above, namely, that , 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 . The same holds for correlated initial states , which are not of the product form (48), as long as they remain in the class 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 ETH194–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 dynamics at leading to Eq. (46), which also implies that the system converges to the stationary state . 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 , as stated in Eq. (3).
A. General dynamical setting
The dynamical point of view for an open system is generally concerned with describing its state evolution when is brought into contact with a heat bath 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 , cf. Eq. (2). One may then address a more refined version of our initial questions (Q):
(Q′) Is the system's dynamical steady state equal to the mean force Gibbs state discussed in Sec. II?
Below, we summarize some key results on steady states of various dynamical systems and make a connection to where possible.
To proceed, consider the general interaction of the form (4). If the operators commute with , 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 does not commute with , then energy exchange processes between the system and bath are enabled.
The initial state, , 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 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 . For product initial states (48), Eq. (50) can be cast into the form of an integro-differential equation for , which is called the Nakajima–Zwanzig master equation.59 For an arbitrary initial state , where and are correlated or entangled, the master equation for will depend on the initial correlation, see e.g., Refs. 103, 187, 206, 208, 210, 211, and 213–216.
With very few exceptions, a master equation for (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 defined in Eq. (3) with a single interaction term in Eq. (4). We will make use of decomposition (23) of the operator X into a sum of energy eigenoperators Xm. In the interaction picture (denoted by the tilde), one has . These quantities define two time scales—one associated with the Bohr frequencies ωm and the other associated with the differences of Bohr frequencies, .
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 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 , where β is the inverse temperature of the bath Gibbs state . G(t) is assumed to vanish as , and this decay defines the time scale, , 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 solving96 the exact dynamics of a damped harmonic oscillator described by Eq. (10) follows a path integral functional integral approach. For general spectral density and coupling strength λ, and a broad class of initial states of the global system, including the global Gibbs state and entangled states, it is shown96 that their steady states are all the same. The initial state is clearly a stationary state of the global evolution, and its reduced system state is . Hence, 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 is obtained in Ref. 217 using Heisenberg–Langevin equation methods. More recent work100 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 , cf. Eq. (48), and in the steady state limit of long times . Then, the above stationary state argument96 for is again applied to prove that the steady state of the N-body system, for initial states , must coincide with the MFG state .
Working from the functional integral description of the system dynamics, it is also possible to construct an exact non-Markovian master equation78 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 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 at all system-bath coupling strengths λ and for general spectral density .
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 thermalization6,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 form260 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 and (small influence of on ), 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 is taken as constant, but and simultaneously . The rescaled time θ is sometimes called the coarse-grained time. For initial product states , 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 (interaction picture indicated by ) is well approximated by some , i.e.,
where the convergence is uniform on the segment for an arbitrary finite Θ and denotes the trace norm. Here, is the solution to the Davies master equation in the interaction picture, and for the coarse-grained time ,140,263,264
where does not depend on θ and is given by267
The positive damping rates are given by , and the term that renormalizes the system energy is
where Xm and 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 , i.e., .
We now discuss the implications of the convergence (53) at coarse-grained times . If one assumes non-degenerate energy levels of and non-degenerate non-zero frequency differences ωmn, see Sec. IV A, then Eq. (55) predicts no coupling between system populations (diagonal entries of in the eigenbasis) and the coherences (off-diagonals). The coherences are found to decay to zero as .
The above assumptions imply that γm satisfy the detailed balance conditions . As a consequence, the steady state of Eq. (54) is then readily shown140,263,264 to be the system Gibbs state τ, cf. Eq. (1). For degenerate energy levels, the steady state is nonunique240,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 , which includes 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 , 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 . 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 . 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 (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 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 . This is entirely consistent with the result of RtE, saying that the stationary state is . It is further shown in Refs. 176 and 275 that by adding to the Davies generator higher order terms,276
the solution of the corresponding master equation has the following properties. First, it is asymptotically exact, meaning that the stationary state is the exact mean force Gibbs state, . Second, the populations (diagonal density matrix elements in the basis) of the system are approximated by those of to for all times . The coherences of the system (off-diagonals) are guaranteed to be approximated by those of to for times t in the windows and for some constants C1, C2. The corrections 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 with a single interaction term , the general form of the second order () 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 ,
and one that generally does not commute with ,
The second line in Eq. (59) contains the damping matrix
The time dependence of the generator implies that the master equation is non-Markovian in the sense that the time evolution operator Λt, defined by , is such that , 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 is θ-independent. However, due to the bath correlation decay, the functions saturate to constant values for times . In this regime, one can replace by and obtain the time independent generator . At large times, the generators and coincide, and since we are interested in steady-state solutions, we will always mean the simpler, time-independent version .
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 Miyashita103 in which the aim is to derive an expression for the time derivative and confirm to second order in λ that the mean force Gibbs state is such that this derivative vanishes, i.e., “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 al66 is to construct a master equation for 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 Miyashita103 consider the time derivative to second order in λ (weak coupling), which is given in terms of and the total system state at an earlier time t0. Like the quantum oscillator case discussed in Sec. IV B, this initial state is chosen as either (i) or (ii) the global Gibbs state . Assuming now that is set equal to the mean force Gibbs state [also to second order in λ see Eq. (20)], Mori and Miyashita confirm103 that the derivative is zero (to second order in λ) for both initial state choices (i) and (ii), i.e., . Interestingly, the same statement is recovered by letting 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 result103 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 at inverse temperature β is a steady state of the Bloch–Redfield dynamics given in (59).
Mori and Miyashita also note that there is, however, an important ambiguity regarding the steady state as follows. The algebraic equation for the steady state yields a steady state accurate to . However, for any second order generator , and for W a traceless operator that is diagonal in the basis, one has and . Thus, for states , one finds that they are also steady state solutions to the same order in λ as itself, i.e., due to the linearity of the evolution in the density matrix ρ, cf. Eq. (59). This ambiguity implies that while the steady state is correct to second order in for the off-diagonal elements in the 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 with those of τ, see Eq. (20).
Fleming et al279 (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 2n order perturbative master equation will yield the diagonal elements of the steady state accurate only to order (while the off-diagonal elements are determined to order 2n). 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 al104 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 to . 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 results100,103,104,279 on the convergence of the dynamics (to second order in λ) to the (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 questioned60 for this reason, though inclusion of a slippage of initial conditions281–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, , has corrections to τ that are partially consistent with , 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 as well as to get the 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 .
It is worthwhile to comment on the above results100,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 for the Bloch–Redfield generator is that the master equation is not of the GKSL form6,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 system6,285,286 at times less than , 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 in Eq. (59) that oscillate, in the interaction picture, at frequencies , where is the system relaxation time. The resultant full-secular-approximation form of the Bloch–Redfield master equation can be readily shown59 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 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 , see Eq. (52). It turns out103,292 that the Bloch–Redfield equation with replaced by 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 develop296 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 given in Eq. (26), the trick is to decompose in the pointer basis from the outset,
Here, 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., for a characteristic Bohr frequency ω may serve as such quantity. Both large λ and small Δmn imply that the decoherence in the “pointer basis” takes place on a time scale much smaller than the relaxation of the populations, , 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 is obtained,296
where the rate constants kmn follow detailed balance,
and are proportional to and decrease exponentially296 with . 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 pm are proportional to . The off-diagonal elements in the pointer basis, , tend to zero for an arbitrary t > 0 as (or ). Thus, the dynamical steady state at is exactly the static MFG state given in Eq. (27). The proof presented in Ref. 296 extends to multiple interactions in Eq. (4), degenerate spectra of the , 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 derived296 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 work118,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 Pm, which is now confirmed analytically.296 The steady state diagonals in the pointer basis were conjectured to be 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 , 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 . This resolves the discussion82,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 ago302,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 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 . The remnant degrees of freedom of the original bath now constitute a reduced bath . The whole point is that the enlarged system may interact only weakly with the reduced bath with a new small coupling parameter . Thus, the impact of on 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 , which is usually just the Gibbs state for the corresponding Hamiltonian of . By tracing out the reaction coordinate, the steady state of the original system 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 is shown to be , where is the mean force Gibbs state for , 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.
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, should be proportional to ωs for for small ω], if the coupling Δ in the system Hamiltonian (28) is smaller than a certain value (depending on the other parameters of the model), then the coupled dynamics converges to the joint equilibrium state 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 formula325 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 limits328–330 represent significantly different theory.
Now, the interaction 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 , 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 .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 , 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 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 is still satisfied. For the (moderately) strong coupling case in Fig. 1, we find very good agreement of the steady state also with the MFG state formula (27) for the ultrastrong coupling limit () 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 , 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 , are only known in a handful of cases.
Beyond the exactly solvable quantum oscillator, explicit expressions for general systems are known for weak [neglecting and higher] and ultrastrong [neglecting 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 connection8 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 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 complex does equilibrate exactly to the global Gibbs state , 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 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 and/or the local convergence345 to 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) temperatures346 .
Another issue is dimension: in all proofs of RtE, the system 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 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 environment347 during a global unitary operation on an 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 -corrections to τ. However, here the issue was that not all -corrections are in fact captured by the -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) couplings342 are important and give rise (o) to open problems concerning their steady states and consistency with . 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 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 , 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.
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.
Conflict of Interest
The authors have no conflicts to disclose.
The dynamical evolution data displayed in Fig. 1 are available from author A. S. Trushechkin upon reasonable request.