On the Husimi Version of the Classical Limit of Quantum Correlation Functions

We extend the Husimi (coherent state) based version of linearized semiclassical theories for the calculation of correlation functions to the case of survival probabilities. This is a case that could be dealt with before only by use of the Wigner version of linearized semiclassical theory. Numerical comparisons of the Husimi and the Wigner case with full quantum results as well as with full semiclassical ones will be given for the revival dynamics in a Morse oscillator with and without coupling to an additional harmonic degree of freedom.


Introduction
The quest for a simplified description of quantum dynamics of many body systems in terms of classical dynamical input has recently become ever more prominent [1]. A lot of the progress that has been made recently in the quantum optics community for bosonic systems is based on the use of the so-called truncated Wigner approximation [2][3][4].
Along similar lines, in the chemical physics community, a classical Wigner dynamics was introduced early on by Heller [5], with a precursor in reaction rate theory by Miller [6]. Historically, the heuristic introduction of this Wigner method, based on classical trajectories preceded its derivation from a semiclassical initial value representation of the propagator by Cao and Voth [7] by two decades. Nowadays, the classical Wigner method in chemical physics is also referred to as linearized semiclassical initial value representation (LSC-IVR) and is prominently used by the Miller group [8][9][10]. Technically, the LSC-IVR method is much easier to apply than the full semiclassical initial value representation (SC-IVR) of, e.g., Herman and Kluk [11], as it linearizes the phase difference between interfering classical trajectories and does not contain any oscillating terms in the integrand (no sign problem). In addition to the application to (reactive) single surface dynamics, also the application of LSC-IVR to electronically nonadiabatic processes i.e., those involving transitions between different potential energy surfaces was shown by Miller and co-authors [12,13]. Furthermore, LSC-IVRs are close in spirit to the diagonal approximation that is used to calculate the smooth part of the semiclassical spectral density in chaotic systems [14].
Whereas in the pioneering works of the Heller as well as the Miller groups, Wigner transforms in the LSC-IVR are used, a new semiclassical framework, introduced by Antipov, Ye and Ananth [15], based on Husimi functions, see also [16,17], can be tuned to reproduce existing quantum-limit and classical-limit SC approximations to quantum real-time correlation functions. So far, the applicability of the Husimi

Theory
In the following, we will introduce the general form of correlation functions that are to be studied and then we will compare different implementations thereof based on classical trajectories.

General Correlation Functions
The time correlation function of two arbitrary operatorsÂ andB with a (Heisenberg) time-evolved operatorB(t) is defined as C AB (t) = Tr ÂB (t) = Tr Â e ī hĤ tB e − ī hĤ t (1) is the Hamiltonian of the system under consideration. This system shall have 2N degrees of freedom in phase space, denoted by (p, q).
Dynamic phenomena in complex systems can be described in terms of real-time correlation functions. The general time correlation function has various applications based on the choice of the arbitrary operatorŝ A andB. Various choices for the operators of the correlation function are as follows: is a projection operator on initial state |Ψ α . This leads to the survival probability which is equivalent to the absolute value squared of the auto-correlation function. In the following, we will assume the initial states to be Gaussian wavefunctions.
where Q is the partition function, β = 1/kT inverse temperature andr the position operator. This is the temperature dependent dipole-dipole correlation function used to study IR spectra. Semiclassical investigations for this case have been performed in [16,18,19] in the high temperature limit β → 0.
This choice leads to the reduced density. For the case of a Caldeira-Leggett model, the transition to classicality and the blue shift of the system oscillator have been investigated in [20], see also [21].
Herex is the position operator in one dimension. This choice is used to calculate time-dependent moments.
We will be evaluating the dipole-dipole correlation function and the survival probability using both Wigner LSC-IVR and Husimi LSC-IVR methods in the following subsection.

Comparison between Wigner LSC-IVR and Husimi LSC-IVR
For the following comparison, we restrict the discussion to the case N = 1, i.e., a single degree of freedom. The generalization to arbitrary N is straightforward.
In the classical, so-called LSC-IVR employing Wigner functions, the expression for the general time-correlation function is given by [5,7] where A W (p, q) and B W (p t , q t ) are the Wigner-Weyl transforms of operatorsÂ andB and (p t , q t ) are classical trajectories in phase space that evolve from the initial conditions (p, q) according to Hamilton's equations of motionq A W (p, q) is defined as: Now in the case of the coherent state based LSC-IVR method employing Husimi functions, the general time-correlation function is given by [15] where A H (p, q) and B H (p t , q t ) are the Husimi transforms of the operatorsÂ andB. A H (p, q) is defined as [22] A H (p, q) = z(p, q)|Â|z(p, q) , where |z(p, q) represents a coherent state with width parameterγ, given in position representation by with a phase convention slightly different from Klauder's [23]. We stress that although the Wigner expression is generally valid, the Husimi version of the correlation function only holds for operatorsB that do not contribute to the phase [15]. In passing, we note that a proof that any operator is determined by its expectation in all coherent states is given in [24]. Now we will examine two cases with the help of both LSC-IVR methods (Wigner and Husimi). Firstly this will be the dipole-dipole correlation function and secondly the survival probability.

The Dipole-Dipole Correlation Function
For the case of dipole-dipole correlation function, we choose the operatorsÂ andB as:Â = exp −βĤ q/Q andB =q. Here exp −βĤ /Q is the canonical density-operatorρ. Therefore, the Wigner-Weyl transform [2] is given by (we are neglecting the partition function Q in the following) and the canonical density-operator matrix element is Because we are investigating the classical limit of correlation functions, a high temperature limit seems justified. For high temperatures, i.e, for small β, a short time approximation to the "imaginary time propagator" e −βĤ is possible [25] and it leads to the expression for the canonical density-operator. Applying the short-time approximation to the canonical density-operator and neglecting an additive term proportional to β, the Wigner-Weyl transform thus becomes Similarly, it can be shown that Hence, the dipole-dipole correlation function comes out to be see also [26]. Now, we will calculate the dipole-dipole correlation function using coherent states. The Husimi transform of operatorÂ is given by where we used a stationary phase argument that allows to treat the operatorq like a constant [16]. Again applying the short-time approximation to the canonical density-operator and substituting the explicit expressions for z(p, q)|x and x|z(p, q) we get, Now, we approximate the potential V x +x 2 by its value at x = x = q [16], which is constant and hence can be taken out of the integral. Therefore, we have The total exponent E inside the integral can now be written in the form where and also c = −γq 2 . The value of the Gaussian integral is and det A = a 11 a 22 − a 2 12 =γ Applying the short-time approximation β → 0 to the expression of b T A −1 b and also to det A, we arrive at the final result for the Husimi transform independent of the coherent state width parameterγ. Similarly, The dipole-dipole correlation function comes out to be Hence, we can conclude that the expressions for the Wigner-Weyl (13) and Husimi (25) versions for the dipole-dipole correlation function case in the high temperature limit coincide.

The Survival Probability
For the case of the survival probability, the operatorsÂ andB are projection operators; i.e.,Â =B = |Ψ α Ψ α |, and it is assumed that the wavefunctions in position space are Gaussians. Let us first derive the survival probability using Wigner LSC-IVR. For the Wigner-Weyl transform, we get Defining a = γ

we can use the 1D version of the Gaussian integral (20) and the Wigner-Weyl transform becomes
and, similarly, Therefore, the survival probability becomes Now let us try to get a first guess at the survival probability using the Husimi LSC-IVR method. We stress that the Husimi version for the correlation is not applicable to the present case [15] but we will get a partial answer by using it anyways. OperatorsÂ andB are again projection operators. Their Husimi transforms are thus given by Now, where Ψ α is a Gaussian centered around (p α , q α ) with the width parameter γ. Therefore, To simplify matters, we considerγ = γ and get for the overlap. Similarly, Therefore, the Husimi transforms are given by and Hence, the survival probability follows to be If we compare the Wigner-Weyl transforms and the Husimi transforms, we see that both expressions are quite different. In the Husimi transform (36), there is the factor of 1 2 multiplied to the terms in the exponential, which is not present in the Wigner-Weyl transform (27). In addition, in the Wigner-Weyl transform we have an additional factor of 2 multiplied with the exponential term. Hence, the two expressions for the transforms are different and we find that A H = A w 2 . Also, the final Husimi result (38) does not give unity at initial time as a simple Gaussian integration shows. We note in passing that if we had considered the width parameters of the coherent state (γ) and that of the Gaussian initial state (γ) to be different would not have resolved matters.
We are therefore led to introduce a new approach of calculating the linearized semiclassical survival probability starting from the Herman-Kluk (HK) propagator.

Quasiclassical Staying Probability Using the HK Propagator
Here we develop a consistent approach to the linearized survival (or staying) probability, using the HK propagator. We follow the path used by Cao and Voth as well as by Miller and his co-workers to derive the Wigner LSC-IVR i.e, making a sum and difference coordinate transformation of the integration variables. Let us consider an initial Gaussian wavefunction for N degrees of freedom, centered at (p α , q α ), given by where we have defined (column) vectors and a matrix via Subscripts denote the N degrees of freedom. Later, subscript 1 denotes the Morse oscillator coordinate and the other subscript(s) correspond to the harmonic oscillator(s) acting as the bath for the system. We stress that the HK propagator to be used below is based on real classical trajectories and therefore can account for interference effects but for the inclusion of tunneling effects (not to be studied in the present manuscript) also non-classical trajectories would be needed.
The auto-correlation function is given by Using the HK propagator [11] to evolve the wavefunction in time, this is given by with the classical trajectories (p t , q t ) and the classical action as well as the prefactor where m 11 , m 12 , m 21 and m 22 are the elements of the so-called monodromy (or stability) matrix M [27], The width parameter matrix of the coherent state |z (p, q) isγ, but here for simplicity we will consider the width parameters for the Gaussian initial state and the ones for the coherent state appearing in the HK propagator to be equal.
We look for a quantity with a classical analog, therefore we consider the "probability to stay" which is given as the double phase space integral sometimes referred to as double HK expression. To make progress with the integration, we introduce the sum and difference variablesq with the reverse transformation For general N, the Jacobian for the variable transformation in the double phase space integral is always unity.
We now expand the trajectories aroundq t (p,q),p t (p,q) up to the first order. Using the reverse transformation in (49,50) this leads to Then we expand the action up to the first order aroundS, as the second order difference vanishes. This yields and an analogous formula for S (which denotes the action depending on the primed variables). From classical mechanics we have for the first derivatives and therefore the action difference becomes The prefactor in zeroth order is given by Collecting terms, the total expression for the staying probability is where the exponent E has still to be defined. From using (49,50), we find the contributions as well as The total exponent (including the action difference) thus is Remembering the defining equations of the stability matrix δp t = m 11 ∆p + m 12 ∆q, we can now do the ∆p and ∆q integrations as they are simple Gaussian integrations. To this end, we write the exponent in the form is a 2N element column vector, whereas is a constant term. Hence the Gaussian integration gives The exact analytic calculation of det A using the method of factorization, as shown in the appendix following the lines of Herman's 1986 paper [28], finally gives This cancels the HK-prefactor absolute square.
Therefore, the remaining integral overp andq is: Before evaluating this expression numerically some remarks are in order: • The result is a linearized semiclassical result and we therefore can call it an LSC-IVR result.

•
Because it is based on Gaussian basis functions it is the correct Husimi version in the case of the survival (or staying) probability,Â =B = |Ψ α Ψ α |. • In contrast to a straightforward application of (38), in addition to the Husimi exponent c, there is a in the exponent that contains stability information, similar in spirit to the semiclassical hybrid dynamics, where part of the degrees of freedom are treated using the thawed Gaussian approximation [29].

•
Furthermore also an overall factor of 2 N appears naturally, ensuring normalization at the initial time t = 0.

Numerical Results
For the following numerical investigations, we consider a Morse oscillator as our system of interest. The potential describing the binding of a diatomic molecule is given by We consider the mass of the system to be unity. D is the dissociation energy, β is the stiffness (or range) parameter. The values of the dimensionless potential parameters that we use here are [30]: D = 150 and β = 0.288, leading to ω = 4.988 for the frequency of (harmonic) oscillations around the minimum. We choose the value for the phase-space center of the Gaussian wavefunction |Ψ α to be q α = 3.5 and p α = 0 and we take the width parameter of the Gaussian as γ = 4.
We note that Hamilton's equations are solved with the classical Hamiltonian in all cases. Using instead the Hamiltonian matrix element between coherent states that is sometimes used for the Husimi case would lead to almost identical results for the potential parameters considered here. The phase space coordinates, as well as the stability matrix elements are then determined numerically by using a symplectic leap frog method [27].
The quantity of interest in the following is the survival probability P(t) that will be followed up to the first full revival time, T rev = 2π/β 2 which, in the present case, is around 60 times the classical period of the Morse oscillator, given by T cl = 2π/ω = 2π/( √ 2Dβ) ≈ 2π/5. The number of trajectories used for the full HK as well as for the LSC-IVR results was 10 4 . The time step chosen in all dynamical calculations was ∆t = 2π/(100 × ω) = T cl /100.

Uncoupled Case
In Figure 1, we compare the full numerical solution to the time-dependent Schrödinger equation (TDSE) with the corresponding semiclassical HK result. It can be seen that the two results coincide very nicely, although the HK result due to the loss of norm is not coming back to unity fully at the revival time, which is around t = 75. The fact that the semiclassical HK propagator can reproduce the interference based revival in a Morse oscillator has been investigated in detail by Wang and Heller [30].

Coupling to a Harmonic Bath Degree of Freedom
Now we proceed and ask the question what happens to the revival dynamics, if we couple the Morse oscillator to a very heavy harmonic oscillator of mass m y = 20 such that the total Hamiltonian reads with ω y = 1.17 and the initial width parameter for the ground state wavepacket in the harmonic degree γ y = m y ω y . The very heavy oscillator shall mimic the interaction with a many degree of freedom bath with smaller masses, i.e., a "condensed phase" environment, see also [15]. First, we look at a moderate coupling strength of λ = 1. The corresponding results from the Wigner and the Husimi approach as well as the Husimi without the b T A −1 b 4 term in the exponent are contrasted in Figure 3. It can be seen that the additional term brings the Husimi result closer to the Wigner one, rather similar to the uncoupled case. At longer times there is just a small amplitude wiggling in the signal due to the coupling.
In Figure 4, we turn to a higher coupling strength of λ = 9 and first compare the full solution to the TDSE with the corresponding semiclassical HK result. It can be seen that, as in the uncoupled case, the two results again coincide very nicely, even the small oscillations are reproduced until around t = 75. The revival is not visible any more, due to the strong coupling to the harmonic degree, however.
In Figure 5, the results from the LSC-IVR calculations (Wigner as well as Husimi) in the strong coupling case are displayed. It can be seen that both also show no revivals, as to be expected. For the very strong coupling the underlying classical dynamics becomes strongly chaotic and the numerical inversion of the A matrix turns singular at certain times for a lot of the 10,000 trajectories we propagated. In order to circumvent this singularity problem, we treated the problematic trajectories without the additional term in the exponent from the time on, when the inversion is problematic.

Conclusions and Outlook
The aim of this presentation was two-fold. The first goal was to highlight that the Wigner-Weyl and Husimi transform version of linearized semiclassical theories can lead to the same final formula, whereas they are quite different in the case of the survival probability, where strictly, the simple Husimi version is not applicable. This fact then established a second goal: to find the correct Husimi-like version of the linearized semiclassical IVR. This goal was achieved by straightforward linearization of the double HK expression for the survival (or staying) probability. In the course of the derivation we also slightly generalized a seminal calculation by Herman (see the appendix of [28]) to the case of width parameter matrices that are not proportional to the unit matrix.
The newly developed formalism was then applied to the revival dynamics in a Morse oscillator. The revival, being a quantum interference phenomenon, could, apart from the numerically converged full quantum results, only be observed in the full HK results, however. Both LSC-IVR variants show a monotonically decaying envelope and no revival. Coupling the Morse oscillator to a heavy bath oscillator leads to the absence of the revival in the full quantum calculations also. In this case the LSC-IVR results are predicting the correct behavior to a large degree and full quantum calculations are not necessary, the system dynamics has undergone a quantum to classical transition due to the coupling to the environment, similar to the findings predicted using the semiclassical hybrid dynamics [20].
For the future, it might be worthwhile to look at even more general correlations function, the so-called out of time order correlators [31] in the (linearized) semiclassical IVR approach.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Calculation of det A
In this Appendix we analytically simplify det A from the determinant of a 2N × 2N matrix to the product of determinants of two N × N matrices, so that it cancels out with the HK-prefactor absolute square. We do so, by going along the lines of the appendix of Herman's paper [28]. Our det A can be written as The N × N matrices L, M and N are given by and the matrices X ± and P ± are defined as These matrices obey P T h I, which follow from the properties of the Lagrange bracket [32] and the fact that the transformation (p,q) → (p t ,q t ) is canonical. Therefore, det A can be written as As shown in Appendix B, the determinant on the right-hand side of the above equation can be expressed as the product of three determinants, F, G, and det (γ), where We multiply the transpose of the matrix corresponding to the determinant F from the left by the product matrix whose determinant is det (γ) and from the right by det γ −1 , such that the determinant F is unchanged. Thus we get Now, multiplying the first column in the determinant of G by ī h and adding it to the second column we get In the paper by Herman [28], the width parameter matrix was taken as γ times identity matrix, i.e., the width parameters for all degrees of freedom was chosen to be the same and equal to γ. While, here in this appendix the width parameter matrix can even be taken as a general symmetric and invertible matrix, i.e., the width parameters for all degrees of freedom don't necessarily need to be equal.

Appendix B. Proof of Factorization of det A
For reasons of completeness, here we briefly sketch the factorization of det A, i.e., we prove with F and G given in (A8) and (A9). Thus, with the block matrices Therefore, what we had wished to prove.