Wigner phase space distribution via classical adiabatic switching

Classical molecular dynamics simulations provide simple and computationally efficient ways to follow the time evolution of condensed phase and biological systems using Newtonian trajectories. The major corrections to classical dynamics calculations in large systems typically arise from quantization of the initial density matrix and, in many cases, nonadiabatic effects, usually associated with transitions among Born-Oppenheimer states. Quantum coherence effects, which are extremely important in small systems, are usually washed out in large biological molecules or in processes occurring in liquid environments. Quantization of the initial density matrix is necessary when the zero-point energy (ZPE) is not negligible and quantum dispersion leads to distributions that differ substantially from the Boltzmann function. In order to capture the system’s dynamics via classical trajectories, the quantized density matrix needs to be converted to a phase space function. There are two similar methods for performing such calculations: the linearized semiclassical initial value representation^1–3 (LS-IVR) (which is also known as the Wigner method^4,5 and which has also been derived by linearizing the path integral⁶), where a phase space function is obtained via the Wigner transform of the initial density, and forward-backward semiclassical dynamics^7–9 (FBSD), where the phase space function is given by the coherent state transform¹⁰ with appropriate corrections. The Wigner density is also required in methods that employ quantum-classical Liouville dynamics.^11,12 Recently, our group has developed a quantum-classical path integral (QCPI) methodology,^13–16 which seems promising as a rigorous simulation tool for simulating quantum mechanical process large molecular or condensed phase systems. Unless zero-point energy can be neglected, QCPI calculations also require knowledge of the phase space distribution of the degrees of freedom comprising the system’s environment.

The first step in all the methods mentioned above is the evaluation of the required phase space transform of the density operator, and considerable effort has been devoted to the development of techniques for this task. Several approaches are available for evaluating coherent state representations of the density, including local harmonic approximations,¹⁷ semiclassical propagation in imaginary time,¹⁸ and numerically exact path integral representations.^19,20 Evaluation of the Wigner function presents a more challenging task. Various approximate schemes have been proposed, including local²¹ and variationally optimized⁶ Gaussian approximations, imaginary time semiclassical evolution,²² and the thermal Gaussian approximation^23,24 (which consists of frozen Gaussian dynamics²⁵ in imaginary time) along with extensions that capture quantum corrections.²⁶ These techniques have been successfully applied to condensed phase systems.^27–30 

The Wigner transform⁴ of a density operator $ρ ˆ$ is given by the Fourier-type integral

(We use one-dimensional notation for convenience; the generalization of Eq. (1.1) to many dimensions is straightforward.) The oscillatory nature of the Fourier transform makes numerical evaluation of the Wigner function an extremely challenging task. In particular, the “sign problem” that arises from the oscillatory Fourier factor prohibits the use of Monte Carlo methods,³¹ leaving no practical alternatives for exact evaluation of the integral in systems of many degrees of freedom.

In this paper, we propose a very simple, approximate method for obtaining the Wigner transform of a density operator, which is applicable to a pure state as well as the Boltzmann density. The basic idea is to use the exact Wigner density for a zeroth order Hamiltonian to populate the phase space of the system and to evolve these phase space points by classical trajectories, gradually switching on the perturbation potential. According to the classical adiabatic theorem, a trajectory that lies on a phase space torus that corresponds to a particular eigenstate will maintain a constant action, remaining on the eigenstate of the evolving Hamiltonian as long as the adiabatic switching procedure^32,33 is carried out very slowly. By adjusting the weight of each trajectory to account for the change in the Boltzmann population resulting from the energy change of the trajectory, we are able to adiabatically evolve the zeroth order phase space density to the one that closely approximates the Wigner function of a more complex Hamiltonian. Since most uses of the Wigner function are in connection with classical trajectories, the adiabatic switching step is very easily incorporated in such calculations and requires little additional effort. Perhaps, the most appealing feature of our method is the simplicity of the algorithm, which (unlike other methods, such as the thermal Gaussian approximation^23,24,26) is based exclusively on classical trajectory input. Yet, we find the results to be very accurate, and the resulting Wigner distribution is in all of our test cases nearly indistinguishable from the one obtained by accurate integration of Eq. (1.1). An additional benefit of our method is that (by construction) the Wigner distribution it generates is invariant under classical trajectory propagation, leading to rigorous preservation of thermodynamic averages.

In Sec. II, we motivate the idea, starting with eigenstates of harmonic and anharmonic Hamiltonians and proceeding to the thermal density matrix. Test calculations on various one-dimensional and dissipative systems are presented in Sec. III, and some concluding remarks are given in Sec. IV.

Consider a harmonic oscillator H⁽⁰⁾ of frequency ω₀,

along with its ground state wavefunction $Ψ 0 ( 0 )$ of energy equal to $E 0 ( 0 ) = 1 2 ħ ω 0$⁠. The Wigner function for the density operator $Ψ 0 ( 0 ) Ψ 0 ( 0 )$ is given by the expression

This function can be rewritten as

where

is the energy of the classical system at x₀, p₀. Equation (2.3) indicates that the Wigner density drops to 1/e of its maximum value at phase space points x₀, p₀ whose energy equals the ground state energy of the system. These points form an ellipse, shown schematically in the Figure 1(a), whose intersection with the position axis is the classical turning points of the oscillator.

Next, consider a harmonic oscillator H of frequency ω = αω₀,

for which the ground state energy is $E 0 = 1 2 αħ ω 0$⁠. The phase space points that correspond to this energy are shown in Fig. 1(b). Compared to the ellipse for the oscillator with frequency ω₀, the new ellipse is squeezed by a factor of $α$ along x and stretched by the same factor along the p axis. However, the phase space area remains the same, equal to the action of the ground state, i.e., $S 0 = S 0 ( 0 ) = 1 2 ħ$⁠.

Now, suppose the harmonic oscillator H⁽⁰⁾ is used as the zeroth order Hamiltonian, while the target system is the harmonic oscillator Hamiltonian H of frequency ω = αω₀. Imagine launching a trajectory from a phase space point x₀, p₀ on the 1/e contour of the Wigner function for H⁽⁰⁾, i.e., with energy equal to $1 2 ħ ω 0$⁠, while gradually increasing the oscillator frequency to αω₀. If this process is carried out infinitely slowly, the classical adiabatic theorem guarantees that the action of the trajectory will remain unchanged as its energy slowly changes to $1 2 ħω$⁠. Thus, the trajectory that initially traversed the ellipse of Fig. 1(a) will eventually be found to traverse the ellipse of Fig. 1(b). The evolution of such a trajectory is shown in Figure 2 for a switching time interval equal to 20 oscillation periods. (This switching time, which is not sufficiently long for accuracy, is chosen for clarity of illustration.)

If the above procedure is repeated with many trajectories whose initial conditions trace out the ellipse of Fig. 1(a), these trajectories will be found to trace out the ellipse of Fig. 1(b) at the end of the adiabatic switching process. Thus, the 1/e contour of the Wigner density for the zeroth order Hamiltonian is transformed to the 1/e density contour of the Wigner function that corresponds to the target Hamiltonian H. The above argument can be extended to phase space points selected to correspond to any density contour of the Wigner function. Thus, one sees that the Wigner density of the zeroth order Hamiltonian will be transformed exactly to the correct Wigner density of the target harmonic system if the perturbation potential is switched on infinitely slowly.

The above ideas can be extended to excited states of a harmonic oscillator. The nth eigenfunction of H⁽⁰⁾ is given in terms of Hermite polynomials according to the expression

The corresponding density matrix is

and its Wigner transform is

Rescaling the coordinates as before, i.e., introducing $x= x 0 / α$ (thus, $Δx=Δ x 0 / α$⁠) and $p= α p 0$⁠,

The last expression is recognized as W_n(x, p), the Wigner function for the Hamiltonian of frequency ω = αω₀.

Next, the “primitive” Wentzel-Kramers-Brillouin (WKB) approximation to the wavefunction of a one-dimensional anharmonic Hamiltonian,

Evaluating the Wigner integral, Eq. (1.1), in the limit ħ → 0 via the stationary phase method, and using $p ( x + 1 2 Δ x ) +p ( x + 1 2 Δ x ) ≃2p ( x )$⁠, one finds³⁴ that the Wigner density within the “primitive” WKB approximation has the form

i.e., the Wigner density is localized near the energy shell. For a bound potential, the WKB wavefunction is a linear combination of “primitive” wavefunctions. In this case, the Wigner density has been shown to be oscillatory near the phase space ridge specified by the energy boundary.³⁵ 

Consider again a point x₀, p₀ at the energy eigenvalue of a zeroth order Hamiltonian H⁽⁰⁾ (for example, the harmonic approximation to the potential of H). Under the classical force from H⁽⁰⁾, the trajectory launched from this point traverses the energy boundary $p ( x ) 2 =2m E − V ( x )$⁠, tracing out a closed curve of area equal to the action $S n ( 0 ) = ( n + 1 2 ) ħ$ that specifies the eigenstate. Upon switching on the perturbation H − H⁽⁰⁾ adiabatically, the trajectory deforms to the energy shell of the full Hamiltonian, preserving the value of the action. By virtue of the semiclassical Wigner transform, Eq. (2.11), it follows that the endpoint of this trajectory will again specify a point on the Wigner density ridge.

Based on the above ideas, we suggest that the adiabatic switching process from the Wigner density of a reasonable zeroth order Hamiltonian should yield a good approximation to the Wigner distribution of the corresponding eigenstate of an anharmonic target system. We point out that the adiabatic switching procedure cannot account for subtle quantum mechanical features of the distribution, such as the small shift of its maximum away from the location of the potential minimum which is often observed in the ground state of asymmetric anharmonic systems.

Last, consider the case of finite temperature. Again, we start with a simple zeroth order Hamiltonian H⁽⁰⁾, for which the Boltzmann density is

and consider its spectral expansion in terms of energy eigenstates,

To motivate the procedure, we express the Wigner transform of this density in terms of the Wigner functions of the individual eigenstates,

Imagine carrying out the adiabatic switching procedure separately for each eigenstate. Upon switching on the perturbation H − H⁽⁰⁾ adiabatically, the phase space points x₀, p₀ distributed according to $W n ( 0 )$ will evolve to points x, p with distribution that corresponds (approximately) to W_n. Thus, the adiabatic switching procedure described so far will evolve Eq. (2.14) to the distribution

However, the Wigner transform of the target density,

is given by the expression

One observes that the Wigner densities W_n resulting from the adiabatic switching procedure need to be weighed by the Boltzmann factors corresponding to the energies of the full Hamiltonian, yet, according to Eq. (2.15), the adiabatic switching trajectories carry weights that correspond to the energies of the zeroth order Hamiltonian. To correct this, we need to readjust the density at the point x, p reached by each trajectory.

The easiest procedure for achieving this task is to include the following classical rescaling factor:

If many quantum states are occupied at the given temperature, the density adjustment given in Eq. (2.18) leads to state occupations consistent with the Boltzmann distribution of the target Hamiltonian. One notices that Eq. (2.18) does not account for the ratio of partition functions Z⁽⁰⁾/Z, because the partition functions are not readily available. However, this ratio is a constant and is easily accounted for by normalization.

However, it is easy to see that the classical weight rescaling procedure described above is destined to fail as the temperature approaches zero. In that case, the density operator approaches the ground eigenstate, which is adiabatically transformed correctly without adjustment of trajectory weights. Clearly, the reason for the failure of the classical Boltzmann weight prescription is energy quantization. In the absence of information about the level spacings of the target system, one can approximately remedy this situation by introducing a rescaling factor that depends on the system’s ZPE, which is often available for the zeroth order as well as the target Hamiltonian. The expression for the Wigner density of a harmonic oscillator,

suggests that the proper rescaling factor has the form

Thus, the adiabatic switching procedure with the ZPE-based rescaling factor produces (apart from a normalization factor) the exact Wigner density at all temperatures in the case of a harmonic potential. In the limit of high temperature, βE₀ ≪ 1, this expression reverts to the classical rescaling factor, Eq. (2.18), which is correct for any Hamiltonian. Since most potentials become nearly harmonic at low temperatures, we expect the ZPE-based rescaling procedure to be accurate at low temperatures and also at high temperature for anharmonic systems. Thus, the largest errors are expected at intermediate temperatures for systems where the energy level spacings of the target Hamiltonian deviate significantly from the harmonic structure. Our numerical tests on strongly anharmonic systems indicate that the ZPE-based procedure is quite accurate at all temperatures. Further, these results (presented in Sec. III) suggest that the simpler classical rescaling procedure yields very satisfactory results under most conditions of practical interest, i.e., systems of many degrees of freedom, where full evaluation of the Wigner integral, Eq. (1.1), is impractical. These systems have a high density of states at typical temperatures such that the classically derived scaling factor is sufficiently accurate.

To summarize the procedure, the Wigner density of the target system is obtained from adiabatic transformation of a zeroth order Wigner density with weight adjustment,

where the weight rescaling factor is given by the classical or the ZPE form, and N is a normalization constant. Since the Wigner density is usually generated for the purpose of calculating time-dependent averages of the type

the normalization factor is evaluated concurrently with the dynamical average.

Quasiclassical methods often are concerned about the inconsistency of generating the initial density by quantum mechanical procedures (whenever this task is feasible) and carrying out the dynamics via classical trajectories, as quantized distributions are not invariant under classical propagation.^36,37 By its nature, the adiabatic switching construction of the Wigner distribution guarantees its invariance under classical dynamics. This desirable feature ensures exact preservation of thermodynamic properties.

In this section, we present numerical examples that illustrate the procedure described in Sec. II.

In the first example, we choose the potential

The zeroth order system is chosen as a harmonic oscillator of frequency $ω h = 2$⁠. Figure 3 shows the (renormalized) Wigner function produced via the adiabatic switching method, both with classical and with ZPE rescaling, at various temperatures. At each chosen temperature, we present the obtained phase space distribution, along with the position and momentum distributions,

These results are compared to accurate calculations generated via a basis set calculation and to the classical Boltzmann distribution.

In spite of the very large anharmonicity, it is seen that the adiabatic transform reproduces the Wigner function rather faithfully at all temperatures. It is particularly encouraging that both the classical and the ZPE weight rescaling factors produce accurate results even far from their optimal regimes. The only observable flaw is the absence of a shift in the maximum of the phase space distribution away (toward the right) of the potential minimum, which is seen in the low-temperature distributions obtained by a basis set calculation. This is so because the highest density contours of the Wigner function correspond to trajectories with energies near the potential minimum.

To quantify the respective accuracy attained by the two weight rescaling factors, we show in Figure 4 the variance of the position space distribution as a function of temperature. The simple adiabatic switching with classical rescaling is seen to produce quantitative results at all but the lowest temperatures. The adiabatic switching procedure with ZPE rescaling is quantitatively accurate at all temperatures. The importance of the rescaling factor becomes apparent by presenting the raw adiabatic switching results. It is seen that in the absence of weight rescaling, the method fails to produce accurate results at moderate to high temperatures.

Last, we demonstrate the time invariance of the distribution generated via the adiabatic switching procedure in Figure 5, which shows the distribution initially and also after classical propagation by 10, 100, and 1000 vibrational periods. Since no change of the distribution takes place, thermodynamic properties will be strictly conserved.

Finally, we use the adiabatic switching procedure to generate the Wigner function and dynamics for a harmonic system coupled to a dissipative bath. The Hamiltonian has the form

where m = 1,  Ω = 2. The frequencies and system-bath coupling coefficients are collectively specified by the spectral density.³⁸ We use the Ohmic form,

with the cutoff frequency ω_c = 1.25Ω. The bath was discretized using 60 oscillators with frequency chosen according to the logarithmic discretization of the spectral density³⁹ with ω_max = 4ω_c. The calculations were performed at four parameter sets, for which the system dynamics changes from weakly damped oscillations to near-monotonic decay.

The Wigner function of the separable zeroth order Hamiltonian was generated by using Monte Carlo sampling³¹ to generate phase space points distributed according to the analytic expression for the 122-dimensional Wigner density. The system-bath coupling was switched on adiabatically over a time length of around 32 oscillations. Trajectory weight rescaling was performed using the classical Boltzmann-weighted procedure.

Once the Wigner distribution was constructed, dynamical results were obtained by continuing the trajectories under the forces specified by the system-bath Hamiltonian. We report the real part of the position correlation function of the system,

obtained from the Wigner function according to the classical procedure

Exact results for comparison were obtained analytically.

Figure 6 shows the adiabatic switching Wigner results, along with exact quantum mechanical results for this Hamiltonian. Also shown are results obtained by approximating the Wigner function by the classical Boltzmann phase space density. It is seen that the Boltzmann density severely underestimates the magnitude of the correlation function because of its neglect of zero-point energy. Even though some of the calculations were performed at a very low temperature with respect to the frequency of the system (and the majority of bath oscillators), the classical Boltzmann-weighted procedure produced excellent results. This is so because the large number of degrees of freedom leads to a high density of states. These findings suggest that the simple adiabatic switching procedure is sufficiently accurate for polyatomic systems.

The idea of adiabatic switching, which follows from the adiabatic theorem of classical mechanics, is very old.³² More recently, adiabatic switching was formulated as a trajectory-based procedure for generating semiclassical energy eigenvalues.³³ Even though the idea is strictly valid only when the state of interest is associated with a torus in the full phase space of the system, numerical studies³³ have shown the method to be quite robust even in the presence of chaotic dynamics. However, practical issues are often encountered due to separatrix crossing or when resonant states are present; thus, special care must be taken to choose the initial Hamiltonian in a way that avoids such crossings, whenever possible.^40–42 Adiabatic switching has been applied to calculate vibrational energies in molecules with several degrees of freedom.^43–45 

The adiabatic switching procedure described in this paper is a simple approximate but quite accurate procedure for generating the Wigner transform of the density operator that is valid for pure states or at finite temperature. Because the target density is the Boltzmann operator (or, at zero temperature, the ground state), all degenerate states are to be included, and thus, resonant states do not present a problem. Even though it is conceivable that the Wigner density of a strongly anharmonic system may have small negative parts even at zero temperature, we have not encountered this situation in any of the strongly anharmonic systems we investigated. At finite temperatures, we expect any small-area negative regions of individual eigenstates to be washed out by the Boltzmann averaging procedure; thus, we do not anticipate the Wigner distribution generated by a positive zeroth order density to lack important structure.

The adiabatic switching based method of generating the Wigner distribution is very easy to implement, as it requires only classical trajectory integration without additional information on potential derivatives. Each trajectory employed in the adiabatic switching process can subsequently continue to be propagated in time under the full Hamiltonian to yield dynamical information. The Wigner distribution produced through the adiabatic switching procedure is invariant under classical propagation, preserving thermodynamic averages. Our numerical tests suggest that the adiabatic switching method is also quite accurate under a variety of conditions. Thus, the method is ideally suited for quasiclassical trajectory calculations and also for calculating the phase space distribution in quantum-classical calculations.

This material is based upon work supported by the National Science Foundation under Award No. CHE 13-62826.