It is pointed out that the classical phase space distribution in action-angle (a-a) variables obtained from a Wigner function depends on how the calculation is carried out: if one computes the standard Wigner function in Cartesian variables (*p*, *x*), and then replaces *p* and *x* by their expressions in terms of a-a variables, one obtains a different result than if the Wigner function is computed directly in terms of the a-a variables. Furthermore, the latter procedure gives a result more consistent with classical and semiclassical theory—e.g., by incorporating the Bohr-Sommerfeld quantization condition (quantum states defined by integer values of the action variable) as well as the Heisenberg correspondence principle for matrix elements of an operator between such states—and has also been shown to be more accurate when applied to electronically non-adiabatic applications as implemented within the recently developed symmetrical quasi-classical (SQC) Meyer-Miller (MM) approach. Moreover, use of the Wigner function (obtained directly) in a-a variables shows how our standard SQC/MM approach can be used to obtain off-diagonal elements of the electronic density matrix by processing in a different way the same set of trajectories already used (in the SQC/MM methodology) to obtain the diagonal elements.

The traditional Wigner function^{1} *W*(*p*, *x*) corresponding to a quantum mechanical (QM) operator $A\u02c6$ is a function of the Cartesian phase space variables (*p*, *x*) (discussed here for a 1D system, in units for which *ħ* = 1) given by

Here we consider the density (matrix) operator

for a set of discrete states {|*N*〉}, *N* = 0, 1, … , so that Eq. (1) gives the corresponding Wigner functions as

For example, for the two lowest states of a harmonic oscillator—which are of specific interest when using the Meyer-Miller^{2} (MM) approach for representing electronic states—the Wigner function for the 2 × 2 density matrix is easily found (by using the ground and first excited state harmonic oscillator wavefunctions) to be

where dimensionless variables have been used for which *ħ* = *ω* = *m* = 1.

In classical and semiclassical theory it is often useful to use action-angle (“a-a”) variables (*n*, *q*), which are an example of a “generalized” coordinate (*q*) and “generalized” momentum (*n*), rather than the Cartesian variables (*p*, *x*). Though it is not possible to define action and angle operators in a rigorous way quantum mechanically, they are canonically conjugate variables in classical mechanics (i.e., their time evolution is given by Hamilton’s equations when the Hamiltonian is expressed in terms of them) and in semiclassical theory (=“classical-limit quantum mechanics”) one can treat them in a fully consistent fashion.^{3} For a harmonic oscillator, (*n*, *q*) and (*p*, *x*) are related by

and if one uses Eq. (5) to express the density matrix Wigner functions of Eq. (4) in terms of (*n*, *q*), one obtains

There is, however, another way to obtain the Wigner functions for the density matrix in a-a variables, one that is in fact more consistent with semiclassical theory (and also performs much better in many applications), and that is to use Eq. (3) directly in a-a variables,

where the bra-ket 〈*q* | *N*〉 (i.e., the “wavefunction” in angle space for the “state” |*N*〉) is given by^{3}

It is then a simple calculation, using Eq. (8) in Eq. (7), to obtain the following result for the density matrix in a-a variables:

where *N* and *N*′ are integers—the semiclassical (Bohr-Sommerfeld) definition of a quantum state—[they are actually half-integers (because of zero point energy), but one usually takes this into account by replacing *n* by $n+12$ throughout, but this has not been done here to keep the presentation as simple as possible]. For diagonal elements of the density matrix (*N*′ = *N*), for which the operator |*N*〉 〈*N*| is the projection operator onto the state with quantum number *N*, one thus has the simple result

which is essentially a statement of Bohr-Sommerfeld quantization: for state |*N*〉, the phase space distribution in a-a variables is that the classical action *n**is* the quantum integer *N* (or half-integer), and the angle variable *q* is random (i.e., the distribution is independent of *q*).

The result given in Eq. (9) is thus quite different from that given for a harmonic oscillator in Eq. (6), e.g., the distribution of the action variable *n* for the ground state (*N* = *N*′ = 0) in Eq. (6a) is extremely broad, unphysically so,^{4} while Eq. (9b) requires *n* to be the integer *N* for that state, i.e., Bohr-Sommerfeld (semiclassical) quantization. (Note also that the diagonal element for the excited state in Eq. (6b) is not always positive.) For the off-diagonal elements (*N*≠*N*′), one sees that *n* is required by Eq. (9a) to be the average of *N* and *N*′, and in this case the density matrix also depends on the angle variable (i.e., has phase information); note that the angle dependence given by Eq. (9a) is actually the same as that in Eq. (6c).

Our current interest in these Wigner functions in a-a variables is for their use in the symmetrical quasi-classical/Meyer-Miller (SQC/MM) approach^{5,6} for treating electronically non-adiabatic processes. The MM approach characterizes each electronic state as a harmonic oscillator with quantum number *N* = 1 or 0 (for an occupied or unoccupied state, respectively). For the case of 2 electronic states, the relevant electronic space of single excitations thus consists of two direct products of 1D harmonic oscillator states: |1〉 ⊗ |0〉 (for electronic state 1) and |0〉 ⊗ |1〉 (for electronic state 2); the Wigner functions for the 2 × 2 electronic density matrix {*ρ _{ij}*} are therefore the product of two of the 1D Wigner functions of Eq. (9),

As discussed in prior work,^{6–10} the SQC model corresponds to replacing each of the delta functions in Eq. (10) by “pre-limit” delta functions, i.e., the quasi-classical “window” functions centered about the appropriate values of the action variables. This not only makes the calculations much simpler but the smoothing effect this entails also brings the classical results into much better agreement with quantum results than the delta functions themselves.

Fig. 1 shows the 2-dimensional window functions centered about the values (1,0) and (0,1) for the two diagonal elements, and also that centered about $12,12$ (red window) for the off-diagonal element. Since the Wigner functions for the elements of the density matrix must be orthogonal (i.e., the phase space integral of the product of any two elements must vanish), the window functions of the two diagonal elements cannot overlap since they are everywhere non-negative (unlike the Wigner functions in Cartesian variables). However the window function for the off-diagonal element can overlap those of the diagonal elements (as in Fig. 1) because it has a phase factor that enforces the orthogonality.

There is of course no unique pre-limit delta function, and this is where the modelistic aspect of the SQC model enters: one uses intuition and experience to choose these window functions to be able to treat as wide a range of situations as possible. The very simple square histogram window functions in Fig. 1, with our preferred zero point energy (ZPE) *γ*-parameter of 0.366 which sets their widths, have been seen to provide very good agreement with accurate quantum results for a number of model systems, but there will undoubtedly be cases for which other choices may be better. (The window function for the off-diagonal element of the density matrix (centered at $12,12$) is taken to have the same width as those for the diagonal elements.)

Eq. (9) thus provides a practical prescription for computing the full electronic density matrix within the SQC/MM framework: one initializes trajectories as usual—by choosing actions (*n*_{1}, *n*_{2}) within a window about the initial quantum state, i.e., a blue window in Fig. 1, with random angle variables—and then collects the trajectories in the various window functions as they evolve according to the MM Hamiltonian (including the phase-factor for the off-diagonal elements of the electronic density matrix). That is, the entire electronic density matrix is obtainable from this one ensemble of trajectories.

A preliminary result calculated with this approach is provided in Fig. 2 for one of the site-exciton Hamiltonians of Ishizaki and Fleming,^{11} which are often used to model light-harvesting complexes. The figure shows the full 2 × 2 SQC/MM-calculated time-dependent electronic density matrix plotted versus benchmark results calculated using the hierarchical equations of motion (HEOM) methodology of Ishizaki and Fleming (which should essentially be equivalent to exact QM results for this spin-boson type model). One sees that the SQC/MM-calculated matrix elements agree extremely well with the exact HEOM results (using our “standard” histogram window functions with *γ* = 0.366 as described above^{12}).

From these proof-of-concept calculations, we conclude that by augmenting the standard SQC/MM approach with the more general result of Eq. (9), we have a practical (and theoretically sound) prescription for obtaining the complete electronic density matrix within the general framework of a classical MD simulation.

As a final note, it is interesting and illuminating to point out the relation of the above Wigner functions in a-a variables, Eq. (9), to the Heisenberg Correspondence Principle (HCP),^{13} which approximates the matrix elements of any operator $A\u02c6$ as

where *A*(*n*, *q*) is the classical function in a-a variables that corresponds to operator $A\u02c6$. (For example, if operator $A\u02c6$ is the Cartesian coordinate operator $x\u02c6$, with *x*(*n*, *q*) given by Eq. (5), then Eq. (11) gives the exact QM matrix elements of $x\u02c6$ for the two lowest states of a harmonic oscillator.) However, another way to express the matrix element $\u3008N\u2032|\u2009A\u02c6\u2009|N\u3009$ is

i.e., as the trace of $A\u02c6$ with the off-diagonal density operator. Evaluating this trace as a phase space average over a-a variables with the Wigner functions in a-a variables (Eq. (9)) gives

which is easily seen to give the HCP of Eq. (11).

Using Wigner functions obtained directly in terms of action-angle variables, Eq. (9), thus incorporates two venerable semiclassical concepts (for systems which possess a-a variables): (1) the Bohr-Sommerfeld quantization condition (i.e., defining quantum states by integer values of the action variable) and (2) the Heisenberg correspondence principle for obtaining approximate matrix elements of any operator between such states.

The authors thank Mr. Herman Chan of the Whaley Research Group at Berkeley for performing the benchmark HEOM calculations for the density matrix in Fig. 2.

Dr. Stephen Cotton also thanks Dr. Jan Roden, formerly of the Whaley Research Group, for useful discussions surrounding potential applications of the methodologies developed in this paper, as well as Professor Birgitta Whaley herself for her interest and encouragement.

This work was supported by the National Science Foundation under Grant No. CHE-1148645 and by the Director, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

In addition, this research utilized computation resources provided by the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

## REFERENCES

See, e.g., Fig. 4 of Ref. 5 and related discussion.

A complete analysis and further applications will be presented in future work.