Path integral approach to the Wigner representation of canonical density operators for discrete systems coupled to harmonic baths

Practical and accurate representations of fully correlated canonical density operators are essential for the determination of both thermodynamic and dynamic properties of many-body systems. The description of the thermodynamics of a system provides access to quantities like entropy, heat capacity, and various susceptibilities, which provide insight into, for example, the nature of equilibrium phase transitions. On the dynamical side, equilibrium time correlation functions, which require sampling from the full equilibrium Boltzmann operator, lie at the heart of the description of linear and nonlinear spectroscopy,¹ the determination of transport coefficients in condensed phase systems,^2,3 and the calculation of chemical rate constants.^4–7 The development of schemes that accurately represent the canonical density operator has been the objective of a large number of theoretical efforts that have, in turn, produced an impressive spectrum of numerically exact^8–17 and approximate methods.^18–24 However, despite significant progress, the calculation of static and dynamical properties of many-body quantum systems remains a challenging task.

For many complex systems, the phase space formulation of quantum mechanics, as encoded by the Wigner distribution, has provided a particularly convenient platform for the investigation of both dynamics and thermodynamics.^25–29 While the phase space formulation provides a rigorous, if generally impractical, protocol for the evolution of operators via the Moyal bracket,^26,28 its utility lies in its compatibility with the semi-classical hierarchy of techniques. The incorporation of the Wigner approach into these approximate methods not only sidesteps the complications associated with the Moyal bracket expansion but also allows for the choice of the level of sophistication and accuracy necessary for dynamical calculations. Indeed, this is an essential factor as the simple Ehrenfest,^30,31 surface hopping,^32,33 and linearized semiclassical initial value representation^34–36 (LSC-IVR) schemes become the only practical approaches for many complex systems. In addition, the phase space framework has also been an integral component in the development of successful hybrid schemes that combine numerically exact quantum approaches or traditional perturbation theories with classical time evolution.^37–41 

Unfortunately, the Wigner transformation of the canonical density for complex systems can rarely be obtained analytically, and its numerical determination contends with the challenge of the highly oscillatory phase associated with the Fourier transform.³⁴ Nevertheless, a variety of approximations have been developed. These range from the simple replacement of the quantum Boltzmann operator with its classical counterpart, an approximation that is only appropriate at sufficiently high temperatures where the zero-point energy is negligible, to sophisticated path integral-based techniques.^42–57 These approaches have proven useful in the investigation of, for instance, vibrational spectra and relaxation rates,^43,58–60 proton transfer problems,⁵⁶ and quantum diffusion in para-hydrogen⁶¹ and liquid neon.^49,62 The benefits of these approximations notwithstanding, the general accuracy of approximate Wigner transformed density operators in complex systems has been difficult to assess, especially when used in conjunction with dynamical calculations.

Here we show that for impurity-type problems where the system-bath coupling is linear in the bath coordinates and the bath can be approximated as harmonic, the Wigner transform of the canonical density operator can be obtained analytically. For this reason, we focus on the simplest nontrivial model that captures the relaxation and dephasing of generic quantum systems coupled to a quantum bath with arbitrary coupling strength: the spin-boson (SB) model. Of course, even for the SB model, the integration of the bath degrees of freedom required by the Wigner transformation cannot be achieved without a Hamiltonian splitting procedure, achieved in the path integral framework by means of the Trotter approximation. Naturally, the formalism provided here is not restricted to the SB model, but is also applicable to any generalization where the bath remains harmonic and the coupling linear in the bath coordinate. To ensure that the resulting density operator can be used in conjunction with quasi-classical methods such as the Ehrenfest and surface hopping schemes as well as with conventional semiclassical methods, we implement only the partial Wigner transform with respect to the bath degrees of freedom. Extension to the full Wigner transformation can be achieved simply through the use of the mapping variable^63–66 or coherent state⁶⁷ formalisms.

Importantly, the present scheme provides a computationally simple approach to the calculation of thermodynamic data for SB-type systems. By providing a numerically exact representation for the initial conditions used in dynamical simulations, this method also represents an important benchmark for the use of approximate Wigner transformed canonical densities both in the static and dynamic contexts. This property allows us to demonstrate that the current method converges rapidly with respect to the number of path integral slices for a large region of parameter space, and that proper rendering of the canonical density can dramatically influence the accuracy of both thermodynamic and dynamic quantities. Interestingly, the analytical expression derived here reveals the canonical density as a linear superposition of bath distributions with weights determined by the paths allowed in configuration space.

It bears remarking that Moix, Zhao, and Cao have previously developed a related and highly efficient approach based on the influence functional formalism for the calculation of the reduced density matrix of a system coupled linearly to a harmonic bath.⁶⁸ The reduced density matrix, which corresponds to the partial trace over the bath degrees of freedom of the full canonical density operator, permits the calculation of thermodynamic averages of any system operator, but precludes calculation of any non-system property. In contrast, by providing an analytical form for the full canonical density operator, our approach permits the calculation of any thermodynamic average, albeit at a higher computational cost. Another advantage of the present work is that, as stated above, it can be used to efficiently and exactly sample the initial conditions required for quasi- or semi-classical calculations of equilibrium time correlation functions. Similar to the work of Moix et al., our work can be easily generalized to N-level systems coupled to harmonic baths and is not limited to any specific form of the spectral density, $J(ω)$⁠.

The paper is organized as follows. In Sec. II, we introduce the formalism used in the paper. Specifically, in Sec. II A, we present a brief review of the phase space formulation of quantum mechanics. Sec. II B introduces the SB Hamiltonian. In Sec. II C we outline the derivation of the Wigner transformed canonical density for the SB model (the extended derivation can be found in Appendix  A). Sec. III contains the results and in Sec. IV we conclude.

As stated in the Introduction, the phase space formulation of quantum mechanics provides a framework that integrates the use of Monte Carlo sampling of initial conditions coupled with trajectory-based methods associated with quasi- and semi-classical approaches. Within this framework, the trace over two operators can be expressed in phase space as^28,29

where A^W(x, p) and B^W(x, p) are the Wigner transformed versions of operators $A^$ and $B^$⁠, which become the functions of the classical coordinate and conjugate momentum variables x and p, respectively, and f is the number of degrees of freedom with respect to which the Wigner transform is performed. The Wigner transform of an operator, $O^$⁠, is defined as^28,29

As Eqs. (1) and (2) suggest, the phase space formulation can be used to obtain static averages when both A and B in Eq. (1) are independent of time, or correlation functions when at least one of the operators is time evolved. In the following, we will be particularly interested in equilibrium time correlation functions of the form

where $ρ=e−βH/Tr[e−βH]$ is the canonical density operator, $β=[kBT]−1$ is the inverse of the thermal energy, and $B(t)=eiHt/ℏBe−iHt/ℏ$⁠.

The Wigner transform for products of operators (e.g., $[ρA(0)]W$ in Eq. (3)) may be expressed as^28,29

where $Λ↔$ is the Poisson bracket operator

and the arrows above the gradient operators indicate the direction in which they act. For notational simplicity, we henceforth set $ℏ=1$⁠. As a final note, we remark on the fact that it is well-known that the Wigner transform of the density operator need not be positive definite.^28,29 This potential complication presents no difficulties in the calculations that follow.

The formalism we develop here is applicable to Hamiltonians consisting of a finite number of discrete states coupled to a noninteracting harmonic bath, with the coupling assumed to be linear in the bath coordinate. The reason for these restrictions is that the current treatment relies on the influence functional approach, which formally eliminates the bath degrees of freedom in the path integral framework.^69,70

While this restriction may seem severe, it is noteworthy that a wide spectrum of problems in the condensed phase may be mapped to such a Hamiltonian. For instance, the discrete degrees of freedom often correspond to a limited subset of the electronic or excitonic manifold coupled to an environment, often idealized as an infinite set of harmonic oscillators. Such Hamiltonians can be written as a sum of system, bath, and coupling contributions, H = H_S + H_B + H_SB. Perhaps the simplest in this class of models is the SB Hamiltonian.^71,72 In the SB model, the system part consists of two discrete states

where $σi$ corresponds to the i^th Pauli matrix, $2ε$ is the bias energy difference between the two states, and Δ represents the off-diagonal coupling between the two sites and is assumed to be static.

The bath consists of independent harmonic oscillators

where P_k, Q_k, and $ωk$ are the mass-weighted momenta, coordinates, and frequency for the k^th harmonic oscillator, respectively. The last term on the right hand side of Eq. (7) is a constant term added for later convenience. As mentioned above, the system-bath coupling term is assumed to be linear in the bath coordinates and antisymmetric with respect to the system

where c_k is the coupling constant describing the strength of the interaction between the system and the kth oscillator, and $α=±1$⁠. The spectral density, $J(ω)$⁠, fully determines the coupling between the system and the bath and is assumed to take the functional form

where the cutoff frequency $ωc$ determines the correlation time for the bath at finite temperatures, and the Kondo parameter, ξ, is a dimensionless measure of the coupling between the system and bath. The Kondo parameter is also proportional to the reorganization energy of electron transfer theory, $λ=ξωc/π=π−1∫0∞dωJ(ω)/ω$⁠, which represents the energy dissipated after the system undergoes a Frank-Condon transition. The functional form for the spectral density in the second line of Eq. (9) corresponds to the often used Ohmic spectral density⁷¹ with an exponential cutoff. We remark, however, that the approach presented here is not limited to any particular form of the spectral density.

Referring back to Eq. (3), it is clear that an expression for $ρW$ is necessary. Because the system part of the Hamiltonian consists of discrete states, ${|0⟩,|1⟩}$⁠, we focus on deriving an expression for an arbitrary matrix element of the canonical density after a partial Wigner transform with respect to the bath degrees of freedom

where $ρa,b(𝐱+𝐬/2,𝐱−𝐬/2)=⟨𝐱+𝐬/2|⟨a|ρ|b⟩|𝐱−𝐬/2⟩$⁠, $a,b∈{0,1}$⁠, N_ab is a temperature dependent normalization constant, and $Ra,bW(𝐱,𝐩)$ is a bath operator of unit trace, i.e., $∫d𝐱d𝐩Ra,bW(𝐱,𝐩)=1$⁠, which can be interpreted as the bath distribution function. We henceforth drop the dependence of the bath distribution function on the bath coordinates and momenta, (x, p), for notational clarity. We also note that we have included the prefactor $[2π]−f$ in the definition of the Wigner transform of the canonical density so that it obeys the normalization condition $∑a∫d𝐱d𝐩ρa,aW(𝐱,𝐩)=1$⁠.

For systems where the total Hamiltonian can be partitioned into two components that are simple to diagonalize, the path integral framework can provide a convenient route for obtaining the exponentiated form for the Hamiltonian necessary for the calculation of propagators and the Boltzmann factor. In this case, we employ the separation adopted previously by Makri and co-workers in the development of the quasi-adiabatic path integral scheme,^73–76 H = H_ad + H_na, where H_ad = H_S and H_na = H_B + H_SB, which refer to the adiabatic and nonadiabatic components of the Hamiltonian. With this partitioning, we rewrite the Boltzmann factor using the Trotter factorization as an N-membered product of basic path integral units

When N is finite, the above equality ceases to be exact and the error it incurs is of the order $O(N⋅exp{−β[Had,Hna]/2N})$⁠. Also note that the Hermiticity of the Boltzmann factor is maintained by the symmetrical splitting in Eq. (11). Using the Trotter decomposition in Eq. (11), introducing resolutions of the identity in the system and bath subspaces $𝟏S=∑a|a⟩⟨a|$ and $𝟏B=∫d𝐪|𝐪⟩⟨𝐪|$⁠, and performing the integrations over the bath coordinates analytically, it is possible to obtain expressions for temperature-dependent (global) normalization factor and bath distribution in Eq. (10). After integration over the bath degrees of freedom, the sequence of spin-variables that characterize the path integral trajectory in configuration space remains, i.e., the sets ${k0,k1,…,kN}$ where $kj∈{0,1}$⁠. To illustrate this, consider the simpler case of treating the isolated subsystem Boltzmann factor via the path integral procedure (with N = 3) such that

where $|kj⟩∈{|0⟩,|1⟩}$⁠. In the following, we refer to individual realizations of the sequence {k₀, k₁, k₂, k₃} as “paths.” Using this notation

where $N$ is a normalization factor, the ratio $W∼a,b/Wa,b$ corresponds to the weighting factors associated with individual paths, $γ∼p$⁠, and $γ∼p$ (⁠$iκ∼p(l)$ and $κ∼x(l)$⁠) are the path-dependent variances (means) for the Gaussian distributions of the coordinate and momentum of the l^th oscillator, p_l and x_l, respectively. In this notation, the tilde denotes that a quantity is path-dependent. Detailed expressions for these quantities and their derivation can be found in  Appendix A.

The interpretation of Eq. (13) is straightforward. The bath distribution function for the SB- and other impurity-type models where the bath is harmonic and the system-bath coupling linear in the bath coordinate can be expressed as a linear combination of Gaussian distributions in the bath coordinates where each contribution is weighted by a temperature- and path-dependent quantity $W∼a,b/Wa,b$ and for which the average displacements of the bath coordinates and momenta are also path-dependent quantities. Importantly, Eqs. (13) and (14) constitute the main result of the analytical manipulations presented in this work. We emphasize as well that the expressions for the canonical density of the SB model derived here may be used to calculate thermodynamic properties and averages and can be easily incorporated into a quasi- and semi-classical descriptions of the equilibrium time correlation functions. This result is also to be considered in light of related treatments of the density operator, in particular the thermal Gaussian approximation⁷⁷ and the Feynman-Kleinert linearized path integral (FK-LPI) treatment.⁴⁵ In both, the Wigner transformed density operator is expressed as a single function rather than a superposition of Gaussian distributions.

In this section, we present some representative results obtained using Eqs. (13) and (14) for thermodynamic averages of spin variables and dynamic calculations of the correlation function, $Czz(t)=Re⟨σz(0)σz(t)⟩$⁠. The dynamics are calculated using the quasi-classical Ehrenfest method, which propagates the system (bath) variables in the time-dependent mean-field of the bath (system) and can be associated with an expansion of the Moyal operator $eℏΛ↔/2i$ to first order in ℏ.⁷⁸ Via comparison with numerically exact results for $Czz(t)$⁠, we illustrate the sensitivity of the Ehrenfest dynamics to the accurate rendering of the canonical distribution.  Appendix B provides details regarding the implementation of the Ehrenfest method.

Before turning to dynamical calculations, we show some representative calculations of thermodynamic averages of the population difference at equilibrium, $⟨σz⟩$⁠, for different realizations of the SB model. Fig. 1 illustrates the convergence of $⟨σz⟩$ with the number of path integral steps for three cases where β or ξ is increased. As is evident from panel (a), N = 0 path integral slices are sufficient to obtain converged results in the high temperature, weakly coupled case. As panels (b) and (c) indicate, with decreasing temperature and increasing system-bath coupling, the number of path integral slices necessary for the converged calculation of thermodynamic averages increases. This is consistent with the fact that the error associated with the Trotter decomposition is of order $O(N⋅exp{−β[Had,Hna]/2N})$⁠, where the contribution from [H_ad, H_na] generally grows with increasing ξ. Remarkably, even for significantly lower temperature and stronger system-bath coupling (⁠$β=10.0$ and $ξ=5.0$⁠), N = 6 is sufficient to obtain converged results. It is also worth noting that with increasing β and ξ, the polarization of the SB model with net bias becomes more severe, as is indicated by the difference in the magnitude of polarization from panel (a) to (b) and (c), and with the faster onset of full polarization with $|ε|$ between panels (b) and (c).

The current path integral approach to the density operator also permits the facile investigation of the dependence of thermodynamic averages on the continuous variation of parameters. Fig. 2 shows the dependence of the population difference as a function of β with the variation of the applied bias ε, the characteristic frequency of the bath $ωc$⁠, and the coupling between the system and bath ξ. Consistent with physical intuition, panel (a) shows that the system becomes more polarized at equilibrium with increasing bias and favors the polarized state with decreasing temperature. The dependence of the results on the variation of the characteristic frequency of the bath shown in panel (b) indicates that the faster the response of the bath (larger $ωc$⁠), the easier it becomes for the system to reach a stable polarized state, corresponding to the formation of a polaron. Finally, panel (c) shows the dependence of the polarization on the system-bath coupling. The results in panels (b) and (c) also agree with physical intuition which indicates that fast baths and strong system-bath coupling promote polaron formation.

The appropriate representation of the canonical density enabled by the path integral approach presented here also facilitates the calculation of equilibrium time correlation functions. For example, Fig. 3 shows the Ehrenfest results for $Czz(t)$ for the unbiased SB model (⁠$ε=0$⁠) obtained using representations of the canonical density that differ in the number of path integral slices employed. Panel (a), which corresponds to a weak coupling, high temperature case, requires only a minimal number of path integral slices (N = 1) for convergence, indicating that the system and bath are indeed approximately independent. Also consistent with our expectations, the Ehrenfest method, which is most appropriate for systems at high temperature and weak system-bath coupling, is able to recover the exact dynamics easily. This picture changes drastically in panels (b) and (c), which correspond to lower temperatures and greater system-bath coupling. In these cases both the Ehrenfest method and the crude approximation for the density operator that treats the system and bath as approximately independent break down. For these panels, the number of path integral steps necessary for the convergence of the dynamics is N = 5 and 6, respectively. It is noteworthy that the accurate rendering of the equilibrium density operator results in improved accuracy for the dynamics for longer times, correctly capturing the slow relaxation in panels (b) and (c), as well as the short-time behavior up to $t=Δ−1$ quantitatively. Finally, we emphasize again that our scheme for the representation of the canonical distribution can be easily incorporated into other quasi- and semi-classical schemes; we have used the Ehrenfest method to illustrate the advantages of the current approach.

In this work, we have derived an expression for the partial-Wigner transformed canonical density operator of the SB model which can be made arbitrarily accurate with increasing number of path integral slices, N. This approach can be used for the evaluation of thermodynamic averages and in conjunction with quasi- and semi-classical evolution methods for the calculation of equilibrium time correlation functions. Importantly, the current work permits the systematic testing of common approximations to the quantum canonical distribution function (e.g., the thermal Gaussian and FK-LPI approaches). Moreover, the generalization of the procedure presented here to an M-level system coupled linearly to a harmonic bath is straightforward.

We have demonstrated the feasibility of the method in the calculation of thermodynamic averages for the spin and bath variables of the SB model, showing their dependence throughout parameter space. Using the current approach with the Ehrenfest method, we have illustrated the sensitivity of the calculated dynamics to the accuracy of the representation of the canonical density operator, which is especially notable in the low temperature and high system-bath coupling regimes. The compatibility of the expressions provided here with quasi- and semi-classical dynamical schemes opens the door to more accurate semiclassical calculations of, for instance, transport coefficients and rate constants. We reserve the investigation of such properties for future publications.

D.R.R. acknowledges support from NSF Grant No. CHE-1464802. A.M.C. thanks Hsing-Ta Chen for useful conversations and Will Pfalzgraff for helpful comments on the manuscript.

To derive an expression for $ρa,b(𝐱+𝐬/2,𝐱−𝐬/2)$⁠, necessary for the Wigner transformation of the canonical density operator in Eq. (10), we first obtain expressions for the matrix elements of the Boltzmann factor using the path integral procedure outlined in Sec. II C. Specifically, we use the Trotter decomposition in Eq. (11) and introduce resolutions of the identity in the system and bath subspaces, $𝟏S=∑a|a⟩⟨a|$ and $𝟏B=∫d𝐪|𝐪⟩⟨𝐪|$⁠, so that the matrix elements of the Boltzmann factor can be rewritten as

where

In this notation

The path integral unit, $⟨𝐐𝐧|e−βHndkn/2Ne−βHndkm/2N|𝐐m⟩$⁠, in Eq. (A3) takes the following form:^73,74

where $δQn(l)=Q(l)−bn(l)$ is the difference between the coordinate of the l^th harmonic oscillator and its displacement due to the system-bath coupling, $Δbnm(l)=bn(l)−bm(l)$⁠, and $θl=βωl/2N$⁠.

With the previous definitions, it is possible to obtain the following expression:

where A^(l) is a tridiagonal $N−1×N−1$ matrix whose diagonal and off-diagonal entries are equal to 2 and $−sech(2θl)$⁠, respectively. For $N<2$⁠, det[A^(l)] = 1. The path-dependent quantities above (marked by a tilde) take the forms

where

and $𝜹∼(l)=𝐣lT⋅𝐀l−1/cosh(θl)$⁠. Also, when N = 0, $κ∼p(l)=κ∼x(l)=Λ∼(l)=0$⁠.

The path-independent quantities take the following forms:

For $N<2$⁠, $ηl=νl=1$⁠.

Substituting Eqs. (A9) and (A2) into Eq. (A1), setting s = 0, and performing the integration in Eq. (A7) leads to the following expression for the partition function:

where the path-dependent weights take the form $W∼a,b=S∼a,bexp[−∑lΛ∼a,b(l)]$ and $Wa,b=∑pathsW∼a,b$⁠.

One final integration over s in Eq. (10) leads to the following expressions for the partial bath distribution and normalization factor

Clearly, the equations derived above have explicitly used the fact that the bath described by the continuous spectral density, J(ω), can be discretized into individual oscillators. To ensure compatibility with the second line of Eq. (9), we use the approach outlined in Ref. 80 which allows us to decompose the spectral density into f oscillators. The frequency of the k^th oscillator for the Ohmic spectral density takes the form

and the coupling constant,

For the results shown here, we used f = 200–300 oscillators.

The Ehrenfest method^30,31,78,81 is a wavefunction-based approach where the system (bath) evolves in the mean field of the bath (system). In addition, this scheme assumes that the bath dynamics are correctly captured by classical mechanics. One may rigorously formulate the Ehrenfest method by first performing a partial Wigner transform with respect to the bath degrees of freedom of the dynamical object to be calculated, e.g., nonequilibrium average or time correlation function

where X_S (X_B) is a generic system (bath) operator.

The heart of the approximation in the Ehrenfest method lies in the dynamical treatment of the operators. In this scheme, the time-dependence is given by the equations of motion for the system and bath. In the case of the system, the wavefunction is evolved via the quantum Liouville equation under the influence of a modified Hamiltonian,

where

is the modified system Hamiltonian and $λcl(t)=α∑kckQk(t)$ is the classical fluctuation of the bias provided by the classical treatment of the bath. Here, $ρS(t)$ is the density matrix for the system. In Eq. (B1), this corresponds to operator $A(0)$⁠. Since the Ehrenfest approach is a wavefunction based method, initial conditions corresponding to coherences, $ρS(0)=|i⟩⟨j|$ where $i≠j$⁠, must first be sampled correctly for the Ehrenfest method to yield appropriate results. Details regarding the generation may be found in Ref. 82.

The equations of motion for the bath variables are given by the classical Hamilton’s equations subject to the time-dependent bath Hamiltonian

where

and $σ¯z(t)=TrS[ρS(t)σz]$⁠.

Given the previous considerations, $Czz(t)$ takes the form

To calculate $Czz(t)$⁠, a second-order Runge-Kutta scheme was implemented. During individual time steps, $σ¯z(t)$ is kept constant for the evolution of the bath, while $λcl(t)$ is kept constant during the evolution of the system. Over a half time step, the equations for the classical variables take the forms

and

where

Convergence for the correlation functions was achieved using $∼5×104–105$ trajectories.