Calculation of the equilibrium state of an open quantum system interacting with a bath remains a challenge to this day, mostly due to a huge number of bath degrees of freedom. Here, we present an analytical expression for the reduced density operator in terms of an effective Hamiltonian for a high temperature case. Comparing with numerically exact results, we show that our theory is accurate for slow baths and up to intermediate system–bath coupling strengths. Our results demonstrate that the equilibrium state does not depend on the shape of spectral density in the slow bath regime. The key quantity in our theory is the effective coupling between the states, which depends exponentially on the ratio of the reorganization energy to temperature and, thus, has opposite temperature dependence than could be expected from the small polaron transformation.

The equilibrium state of open quantum systems has attracted significant scientific attention from the early days of quantum theory^{1} and remains of considerable interest in the present day.^{2–7} This issue is relevant in a variety of contexts, from initial states of emission spectra in molecular systems,^{8,9} multichromophoric Förster transfer rates,^{10–12} nonequilibrium steady states in quantum transport problems^{13} or quantum heat engines,^{14} and decoherence problems in quantum information science^{15} to quantum corrections of thermodynamic laws.^{16} Unfortunately, analysis of the equilibrium state is hindered by an enormous amount of degrees of freedom of the bath.

When system–bath coupling is negligible, the equilibrium reduced density operator $\rho ^eq$ of the system is represented by the canonical Boltzmann distribution of the bare quantum system. On the other hand, this is no longer the case if the coupling between the system and the bath becomes stronger. These issues have been investigated both analytically and numerically,^{3–5,17} and there exist straightforward though numerically somewhat expensive methods to calculate $\rho ^eq$ exactly.^{3,6} Unfortunately, no simple analytical formulas are currently available in the literature. This is rectified in this letter, where we derive an analytical expression for the reduced equilibrium density operator that is valid for high temperatures and slow baths. Our result is intuitively understandable, providing a mental picture for this complex issue.

Previously, we have demonstrated non-trivial spectroscopic manifestations of the non-canonical effects and showed numerically that the effects of the bath upon the equilibrium of the system can be accounted for by introducing effective coupling between system states, *J*_{eff}.^{2} It was shown for a molecular dimer that the effective coupling could be well described by the relation

with *J* denoting the original bare coupling, λ_{1,2} being the reorganization energies of the system states, and *a* being a numerical constant, the value of which was determined by the fitting procedure to be ∼0.15. Here, we provide a theoretical justification for these results.

Let us start by defining the total Hamiltonian, which is comprised of the system, bath, and system–bath interaction parts, respectively,

The Hamiltonian of the system can be expressed as *Ĥ*_{S} = *Ĥ*_{ε} + *Ĥ*_{J}, with

denoting the energies of the states (or sites) and

denoting the couplings between the states. *N* numbers the states of the system.

The bath consists of an infinite set of harmonic oscillators

which interact linearly with system states,

Here, *ω*_{j}, $x^j$, $p^j$ are the frequency, dimensionless position, and dimensionless momentum of the *j*th bath oscillator, while *d*_{nj} is a dimensionless coupling between the *j*th bath oscillator and state *n*. We have also introduced bath operators $F^n=\u2212\u2211j\omega jdnjx^j$. Note that we set *ℏ* = 1 throughout the paper. System–bath interaction induces fluctuations of system state energies that can be parameterized by the spectral densities

For simplicity, we assume no correlations between the fluctuations of different sites, $Cmn\u2033\omega =\delta mnCn\u2033\omega $. The overall system–bath coupling strengths are described by the reorganization energies

For convenience, we also define the reorganization energy operator, $\Lambda ^=\u2211nN\lambda n|nn|$. Calculations of the reduced equilibrium density operator require the knowledge of the imaginary time bath correlation function

with

denoting the operator $F^m$ in the imaginary time interaction picture with respect to *Ĥ*_{B} and $\rho ^Beq$ being the equilibrium bath density operator. Tr_{B} denotes the trace over the bath degrees of freedom.

The total equilibrium density operator of the whole system-plus-bath supersystem should correspond to the Boltzmann distribution

with $\beta =kBT\u22121$ and $Z=Tre\u2212\beta \u0124$ being the partition function. The reduced density operator of the system only is obtained by taking the trace over the bath degrees of freedom,

Unfortunately, calculations of this quantity are very difficult. In this letter, however, we derive a simple analytical expression for $\rho ^eq$, requiring only the assumption of high temperature and, thus, keeping correction terms up to the order of *β*^{2}.

This is accomplished as follows. First, we rewrite the relevant exponential operator as

Here, exp_{+} denotes the time-ordered exponential^{18} and $I0$ denotes the imaginary time interaction picture with respect to *Ĥ*_{0} = *Ĥ*_{ε} + *Ĥ*_{SB} + *Ĥ*_{B}. The equilibrium reduced density operator can, therefore, be written as

The trace is calculated by expanding the time ordered exponential and keeping only the terms of the order *β*^{2}. After lengthy but straightforward derivations, we arrive at (see the supplementary material for details)

Here, *Î*_{S} is the identity operator in the system space. The key step is to notice that this expansion can be resummed into a double exponential function. Thus, we finally obtain

with an effective Hamiltonian

The normalization factor is defined as $Zeff=Tre\u2212\beta \u0124eff$. Resummation techniques have been recently applied for a variety of problems in the theory of open quantum systems.^{19–21}

Equations (16) and (17) constitute the main result of this work. They imply that the equilibrium distribution is achieved according to energies *ε*_{n} − λ_{n} and couplings

This result provides justification for our previously obtained expression for *J*_{eff} [Ref. 2 and Eq. (1)], with coefficient $16$ being very close the numerically obtained value of 0.15.

Our expression for *J*_{eff} paints the following physical picture. With increasing system–bath coupling strength λ, excitation sinks deeper in the potential well corresponding to each site; thus, interaction between the sites becomes irrelevant. On the other hand, with increasing temperature, excitation can move more freely on the potential surface and can once again “sense” neighboring sites.

To validate our results, we performed $\rho ^eq$ calculations using the numerically exact stochastic path integral (SPI) approach of Ref. 3, which is based on the unraveling of the influence functional. From the obtained expression for the effective Hamiltonian [Eq. (17)], it follows that $\rho ^eq$ does not depend on the exact form of the spectral density. To confirm this, we performed calculations using three forms of spectral density, namely, the Debye spectral density

where parameter *γ* is the bath relaxation rate, the Ohmic spectral density

and the super-Ohmic spectral density

where in the latter two cases, parameter *ω*_{c} is the cut-off frequency.

For simplicity, we performed simulations for a dimer system, with Δ*ε* = *ε*_{1} − *ε*_{2} = 100 cm^{−1} and *J* = 100 cm^{−1}. We also assumed equal reorganization energies for both sites, λ_{1} = λ_{2} = λ, and performed calculations for several values of λ. In case of the Debye spectral density, we set *γ*^{−1} = 100 fs, and for both Ohmic and super-Ohmic spectral densities, we set *ω*_{c} = 20 cm^{−1}. These values correspond to the slow bath regime (see below). The temperature for calculations was set to 300 K. Interestingly, calculated equilibrium populations in the site basis remained virtually constant even for λ values up to 1000 cm^{−1}. Evaluations of Eq. (16) predict that in case of equal reorganization energies, site basis equilibrium populations are not affected by the non-canonical effects, and this was supported by our numerical calculations. On the other hand, coherences showed significant differences. This is demonstrated in Fig. 1, where we plot the dependence of $\rho 12eq$ values on a dimensionless parameter $2\lambda /kBT$. As this parameter increases, the true value of $\rho 12eq$ goes further away from the value corresponding to the Boltzmann distribution. As predicted by our theory, the results do not depend on the form of spectral density in the slow bath regime. The theoretical values agree extremely well with numerical calculations up to $2\lambda /kBT\u22723$ (λ ≲ 300 cm^{−1} at *T* = 300 K) and reasonably well up to $2\lambda /kBT\u22726$ (λ ≲ 600 cm^{−1} at *T* = 300 K). Application of our theory for larger systems is discussed in the supplementary material.

We have investigated the relationship between the accuracy of our theory and the relaxation time scale of the bath, *τ*_{rel}, which can be well approximated by the inverse of the frequency corresponding to the maximum of the spectral density function. For the Debye spectral density, we have *τ*_{rel} = 1/*γ*, for the Ohmic spectral density, we have *τ*_{rel} = 1/*ω*_{c}, and for the super-Ohmic spectral density, we have $\tau rel=1/3\omega c$. On the basis of extensive calculations (not shown), we conclude that our expression for the equilibrium density operator is accurate for *τ*_{rel}/*β* ≳ 1. Moreover, numerical investigations show that the slower baths result only in a decrease of factor 1/6 in Eq. (18), but the overall form of the expression for *Ĥ*_{eff} remains the same. This is illustrated in Fig. 2 for the Ohmic spectral density. For slow baths (*T* = 300 K, *ω*_{c} = 100 cm^{−1}, resulting in *τ*_{rel}/*β* ≈ 2), the effective resonance coupling, calculated from the equilibrium density operator, follows the prediction of Eq. (2). If the bath is made significantly faster (*ω*_{c} = 800 cm^{−1}, resulting in *τ*_{rel}/*β* ≈ 0.3), we have the same exponential relationship, but with different coefficient instead of 1/6. If we significantly increase the temperature for this fast bath (*T* = 800 K, resulting in *τ*_{rel}/*β* ≈ 0.7), the dependence of *J*_{eff} again comes very close to that predicted by Eq. (2) and would agree completely for even higher temperatures. Therefore, we can conclude that fast, or Markovian, baths do not induce corrections to the equilibrated state of an open quantum system, while the effects are clearly visible for slow baths. This observation is in line with attempts to relate Markovianity with system–bath entanglement.^{22}

Despite the simplicity of our final expressions, in the slow bath regime, they have considerably larger range of validity than the results based on the second-order perturbation theory^{23} with respect to system–bath coupling (the relevant expressions are listed in the supplementary material), even though the latter is more expensive numerically. This is illustrated in Fig. 3, where it can be seen that the perturbation theory approach agrees well with numerically exact results only for $2\lambda /kBT\u22722$, whereas our theory is considerably more accurate throughout the demonstrated parameter regime. The accuracy of our approach can be traced to the resumation of the perturbative series (see the supplementary material).

Our results imply that the eigenbasis of the bare system Hamiltonian is no longer the eigenbasis of the equilibrated system. The true eigenbasis changes in time due to the influence of bath degrees of freedom. This process is equivalent to the basis rotation picture used in Ref. 4. Since for molecular systems, the eigenbasis of the bare system is denoted as the *exciton* basis, following our previous works,^{2,24,25} we denote the true eigenbasis as the *excitonic polaron* basis. To quantify the change of basis numerically, we have calculated the amount of mixing between the two bases, $\theta =\u27e8p1|e1\u27e92$, where |*e*_{1}⟩ denotes the lowest energy excitonic state and |*p*_{1}⟩ denotes the lowest energy eigenstate. The dependence of mixing has been calculated for the previously defined dimer system at 300 K and is shown in Fig. 4. Clearly, the polaronic effects are subtle, but they tend to increase when the dimensionless parameter $2\lambda /kBT$ increases and that might be enough to break some specific symmetries in the system, as demonstrated in Ref. 2. Using the theoretical effective Hamiltonian of Eq. (17) slightly overestimates these effects. As suggested in Ref. 25, the change, or rotation, of the relevant basis during the equilibration might lead to fast transfer processes, that can be observed spectroscopically. The polaronic effects can be expected to become more significant in case the reorganization energies of system states differ substantially, due to the $\u0124\epsilon \u2192\u0124\epsilon \u2212\Lambda ^$ change in the effective Hamiltonian of Eq. (17).

The theoretical approach developed here has straightforward practical implications. The obtained analytical solution for the equilibrium density operator could be used as a benchmark for approximate theories, for both long time limit of real time dynamics or for imaginary time equilibrium solutions, absolving the need to have exact numerical results. This type of strategy for testing accuracy of theoretical approaches for open quantum systems has been employed previously.^{7,26,27} Our approach is also a very inexpensive way to calculate the reduced equilibrium density operator that is a prerequisite for simulations of emission spectra by the hybrid cumulant expansion method proposed in Ref. 9.

At this point, we would like to contrast our result for effective coupling with that coming from the small polaron transformation.^{28–33} Indeed, when applying the small polaron transformation (also in the limiting case of the variational polaron transformation), one obtains that the coupling becomes dressed with bath variables and the thermal averaging results in bath-renormalized coupling

Even though the origins of *J*_{polaron} and *J*_{eff} are quite different, quite often *J*_{polaron} is taken as describing the equilibrium state of the system. For this reason, it is worthwhile to compare it with our theory. There are two essential differences between the expression for *J*_{polaron} and our result of Eq. (18). First, the polaron result shows that the bath-renormalized coupling becomes zero for spectral densities that exhibit behavior $C\u2033\omega \u2009\u223c\u2009\omega $ for small frequencies due to the divergence of the integral in the exponential. This suggests a sharp transition in the dynamical and equilibrium behavior of the system between cases of Ohmic and super-Ohmic spectral densities. On the other hand, our results show that the equilibrium state is completely insensitive to the exact form of the spectral density if the bath can be assumed to be slow. Second, Eq. (22) gives completely different temperature behavior. At high temperatures, $coth\beta \omega /2\u2009\u223c\u20092/\beta \omega $; thus, *J*_{polaron} decreases at *higher* temperature. As follows from Eq. (18), however, *J*_{eff} decreases at *lower* temperature. We suggest that in some cases, *J*_{polaron} might be misrepresented as being related to the equilibrium properties of the system when, instead, it is just a parameter that is useful to decrease the perturbation term. Indeed, polaron transformation leads to an accurate equilibrium state for some parameter regimes.^{4}

In conclusion, we have derived an accurate expression for the equilibrium reduced density operator of an open quantum system that linearly interacts with a harmonic oscillator bath at high temperatures. The equilibrium state can be described in terms of an effective Hamiltonian. Our results, however, should hold also for cases when the system–bath interaction is non-linear, provided that the bath is slow and the interaction term can be written as $\u0124SB=\u2211nF^n|nn|$ with $TrBF^n=0$, because the assumption of linearity was not used in our derivation. By comparison with numerically exact results, we have demonstrated that our theory is accurate for slow baths and up to intermediate values of reorganization energies. In these conditions, the equilibrium state was shown to be the same regardless of the exact form of spectral density. The effective coupling between the quantum states has contrary temperature dependence to what could be expected from the polaron transformation. We have also shown that for slow baths, our result has a significantly larger range of validity than the usual second order perturbation theory with respect to the system–bath coupling. We expect that our insights will prove relevant for analysis of equilibrium fluorescence spectra or non-equilibrium steady states, characterization of the system–bath entanglement and Markovianity, and interpretation of dynamical spectroscopy results.

See the supplementary material for details of the derivation of the expression for the reduced density operator, analysis of the resummation approach, expression of *J*_{eff} from $\rho ^eq$, and summary of the second order perturbation theory for the equilibrium reduced density operator.

Computations were performed using the resources of the High Performance Computing Center “HPC Sauletekis” at the Faculty of Physics, Vilnius University.