# Quantum initial condition sampling for linearized density matrix dynamics: Vibrational pure dephasing of iodine in krypton matrices

###### Abstract

This paper reviews the linearized path integral approach for computing time dependent properties of systems that can be approximated using a mixed quantum-classical description. This approach is applied to studying vibrational pure dephasing of ground state molecular iodine in a rare gas matrix. The Feynman-Kleinert optimized harmonic approximation for the full system density operator is used to sample initial conditions for the bath degrees of freedom. This extremely efficient approach is compared with alternative initial condition sampling techniques at low temperatures where classical initial condition sampling yields dephasing rates that are nearly an order of magnitude too slow compared with quantum initial condition sampling and experimental results.

## I Introduction

Detailed analysis of the vibrational spectroscopy of chromophores in solution can be used to investigate intermolecular interactions in condensed phases. The theory of vibrational pure dephasing and its contribution to spectral line shapes and shifts has been worked out in detail in various limits. Perturbation theory (see references Diestler and Zewail (1979a, b); Zewail and Diestler (1979c); Skinner (1988) and the and the literature cited therein), which assumes that the interaction between vibrational degrees of freedom and the environment is weak, gives an expression that enables the pure dephasing time to be computed from the zero frequency component of the time correlation function of the fluctuations in the energy gap between the vibrational levels. In many experiments, however, there may be strong initial environmental interactions associated with how the system is initially excited so techniques to study dephasing and dissipation beyond the limits of perturbation theory are important, and several non-perturbative ways to study vibrational dynamics have appeared in the literature Skinner and Hsu (1986); Skinner (1988); Hsu and Skinner (1985); Weiss (1999); Makarov and Metiu (1999). One very fruitful way to go beyond perturbation theory, for example, is to employ an idealized model which can be worked out analytically such as a two level system appropriately coupled to a harmonic bath for which the effects of environmental dephasing on lineshape and spectral shifts have been worked out in detail Skinner and Hsu (1986); Skinner (1988); Hsu and Skinner (1985). Several general predictions emerge from this analysis concerning, for example, how lineshapes are affected by characteristics of bath spectral density, and temperature, etc. Experimentalists can use the predictions of this theory to interpret their findings in terms of the nature of the underlying interactions in condensed phase systems. The interactions between the vibrational coordinate and the environment can depend on vibrational excitation in a complex way and this can complicate use of the two level system theory described above. Also employing a harmonic bath when anharmonicity may be significant could make the use of such a theory questionable.

An alternative non-perturbative method that has received considerable attention in the quantum optics literature Dalibard et al. (1992); Dum et al. (1992) and recently has become a focus in molecular science applications Makarov and Metiu (1999) is the so-called stochastic wave function approach Weiss (1999). With this method the evolution of the density matrix is assumed to take a simplified Bloch form parameterized by phenomenological coefficients governing the decay of coherence and the transfer of population. Rather than propagating the density matrix elements, a stochastic wave function represented in terms of relevant basis functions is evolved using statistical rules designed in such a way so that the reduced density constructed from an ensemble of evolved stochastic wave functions reproduces the evolution of the Bloch equations. Since this approach relies on a Bloch model form containing phenomenological parameters it thus provides an efficient, linear scaling approach for fitting experimental data to such models. As these models are often discussed in terms of gas phase collision processes, interpreting such fits in a meaningful way and extracting information about a microscopic mechanism of the relevant decoherence processes operating in condensed phase systems is not straightforward.

A further possible alternative approach is the use of a realistic microscopic model with vibrational state dependent interactions as has been developed in various contexts over the last few years Fredj et al. (1996); Jungwirth et al. (1997); Jungwirth and Gerber (1999). Such an approach requires a quantum dynamical treatment of dephasing and microscopic simulation methods to address this problem have recently been developed Poulsen et al. (2003); Shi and Geva (2003b, c, d); Hernandez and Voth (1998); Miller (2001); Causo et al. (2005b); Martens and Fang (1997). In this paper we extend these methods to include quantum initial condition sampling that should accurately capture the relevant underlying spectral density of the realistic model system. Thus unlike the theories discussed above which are built on use of different models (e.g. Debye or pseudolocal mode spectral densities), with the approach employed here we do not need to assume anything about the spectral density arising from the interactions. As such direct, model free, interpretation of the experimental results is possible.

The specific set of experiments which we will study with this alternative approach come from recent work from Apkarian and coworkers Almy et al. (2000); Karavitis et al. (2001a); Bihary et al. (2001b, 2004a, 2004b); Karavitis and Apkarian (2004c); Karavitis et al. (2005a); Kiviniemi et al. (2005b); Segale et al. (2005c) in which they use Time Resolved Coherent Anti-Stokes Raman Scattering (TR-CARS) to directly probe the dephasing of vibrationally excited I wave-packet components due to interactions of this chromophore with a rare gas matrix. This well studied system is chosen as a benchmark on which to evaluate the approach due to the availability of accurate interaction potential data. Also, as distinct from traditional frequency domain experiments in which lineshapes and frequency shifts are indirectly interpreted in terms of the evolution of the density matrix, these time domain experiments can be directly connected to the theory of decoherence which we outline below.

The phenomenon of vibrational dephasing can be probed when a vibrational subsystem, described by Hamiltonian , with eigenstates given by , is prepared in some coherent superposition state, , for example, in the presence of an environmental subsystem in state which is unaffected by this preparation. Thus the composite system initial state can be described by a separable product state . Important deviations from this separable product initial condition in the limit of strong system - environment coupling have been discussed Pechukas (1994); Laird et al. (1991); Laird and Skinner (1991). If the quantum subsystem and the environment are uncoupled, the initial coherent superposition of vibrational states will be maintained as the whole system evolves, i.e. the amplitudes and relative phases of the component state contributions will be constant, and the composite system wave function will remain separable. In the more general situation, however, the dynamics of the full system is governed by the coupled system-bath Hamiltonian, and the above initially separable wave function will evolve into an entangled state in which the amplitudes and relative phases of the different vibrational state contributions will change with time. In the last equality in the previous result we have projected the composite state onto the diabatic vibrational basis states , and the bath states that appear in this entanglement, labelled according to this basis, are obtained as . The timescale for the variation in the phase of the different vibrational state components is known synonymously as the vibrational dephasing or decoherence time. The variation in the amplitude, on the other hand, reflects the vibrational state population relaxation time. With this description, the bath states contain all the relative phase and amplitude information of the contributions from the different vibrational basis states. With the view outlined here we are assuming that the details of the shapes of the exciting radiation field pulses can be neglected and that we can focus our attention on the evolution of the reduced density matrix. A complete description of the experiments would require considering details of the time dependent interactions of the full system with the radiation field and is beyond the scope of the current investigation of vibrational pure dephasing times.

In the vibrational pure dephasing limit it is assumed that the composite system dynamics results in slow amplitude relaxation of the chosen vibrational basis states on the fast timescale of the fluctuations in the phases of the different component basis states. Thus we suppose that contributions to final state amplitude that originate from initial states , different from , are negligible and we approximate the above expression as . This is equivalent to the vibrationally adiabatic approximation. The timescale for the vibrational pure dephasing process is governed by the distribution of fluctuations in the phase of the state which is determined by two factors: (1) the strength of the interactions between the vibrational subsystem and its environment which will in general depend on the particular vibrational state , and (2) the distribution of initial states of the environment . This paper develops a general approach for computing vibrational pure dephasing rates in condensed phase systems which incorporates both the varying strength of environmental interactions with quantum subsystem state, and the effect of quantum dispersion on the nature of the distribution of the initial environmental states.

The quantum vibrational state dependent intermolecular potential idea developed in the late 1990’s Fredj et al. (1996); Jungwirth et al. (1997); Jungwirth and Gerber (1999) provides a framework for the representation that we employ here to account for the the variation of environmental interaction strength with vibrational state. Specifically the model of Martens and co-workers Riga and Martens (2004, 2005); Riga et al. (2005, 2006) is adopted for this purpose. The main goal of the work here is to explore how thermal and quantum fluctuations of the initial environment influence vibrational pure dephasing in a model system with realistic interactions. Thus we will compare classical and approximate quantum descriptions of the thermal environmental initial distribution and explore their effects on dephasing dynamics. In addition a complete and general derivation of the Feynman-Kleinert-Wigner (FK-Wigner) approach for sampling quantum initial conditions is presented.

The paper is organized as follows: First we outline a general approximate approach for treating the evolution of the vibrational reduced density matrix. With this linearized approximation, quasi-classical trajectories are evolved from quantum initial conditions sampled from the Wigner transform of the initial density. The main methods section of the paper summarizes the FK-Wigner approach used for this sampling. Details of the complete derivation of this approach are given in the appendices. The final methods subsection makes comparisons of this approach with other published techniques. Finally in section III we outline the implementation of these methods for application to computing vibrational pure dephasing rates and compare these data computed using different approximations with available experimental results. Concluding remarks are given in section IV.

## Ii Methods

### ii.1 Vibrational Pure dephasing

The density matrix formulation offers a convenient way to rewrite the wave function description outlined above in a useful form for making further approximations and interpretations. As noted earlier, the experiments of interest involve exciting the vibrational subsystem independently of the environment. The bath degrees of freedom are thus assumed to be initially prepared in thermal equilibrium with the quantum subsystem. The quantum vibrational subsystem is then prepared by rapid pulsed laser excitation, for example, in a non-equilibrium coherent superposition of excited vibrational states as described in the Introduction. We assume that the non-equilibrium composite system is thus initially described by a density operator which is a product form , where is the bath part of the equilibrium density operator for which we will develop approximations, and the non-equilibrium quantum subsystem density for the initially excited coherent superposition state is , which has various component operators, for example, , that depend on what states of the quantum subsystem are coherently excited by the laser pulses.

The coherently excited composite system will evolve from this factored initial state to an entangled state as a result of the coupled full system evolution as discussed above. This entanglement will thus be described by the full system time dependent density matrix with elements in the environmental coordinate and vibrational subsystem state representation given as . The experiments of interest probe only the quantum subsystem states so we will study the reduced density matrix elements obtained by tracing over all the bath degrees of freedom i.e. . Suppose the preparation selects out the component sub-system density operator initially, then the reduced density operator matrix elements at time will have the form

(1) | |||||

and involve forward and backward propagator matrix elements as well as initial bath density operator matrix elements.

Suppose the full Hamiltonian is

(2) |

which can be expressed in our quantum subsystem diabatic vibrational basis as where . Here is the system-bath interaction potential. In the case of vibrational pure dephasing we suppose that for , i.e. the off-diagonal elements are small compared to the diagonal elements, so we can approximate the full system Hamiltonian by the diagonal form for all important . In this case there is no population relaxation between our vibrational basis states so the propagator matrix elements appearing in Eq.(1) can be written in the composite bath subsystem path integral forms as follows:

(3) |

and

(4) |

where the forward path action, for example, is

(5) |

and a similar expression for the action along the backward path is obtained by modifying the vibrational state accordingly to .

Combining these expressions, the reduced density matrix in the pure dephasing limit can be written as

(6) |

In the above expression the bath evolution is still described at the full quantum level. To proceed to a computable expression for the time dependence of the reduced density matrix elements we follow various authors Poulsen et al. (2003); Shi and Geva (2003b, c, d); Hernandez and Voth (1998); Miller (2001) and combine the forward and backward propagators in Eq.(6) and write the product in terms of mean, , and difference, , bath subsystem path coordinates (with a similar transformation for the bath momenta, and ). Next, the phase of the integrand in Eq.(6) is expanded in the path difference variables. In condense phase systems various arguments can be given to justify keeping only low order terms in the path difference Caldeira and Leggett (1983a); Caldeira and Leggett (1983b); Caldeira and Leggett (1985); Poulsen et al. (2003); Causo et al. (2005b). Thus we proceed by truncating the expansion of the phase to linear order in the path difference variables. With this approximation the integrals over the initial difference coordinate can be performed defining the Wigner transform of the initial density operator

(7) |

If we discretize the paths appearing in Eq.(6) by inserting resolutions of the identity in the bath subsystem phase space and take small time steps , thus writing the discrete path variables as , for example, we find that all integrals over the difference coordinates, , and difference momenta, for that appear in the discretized result can also be performed in the linearized approximation since they are integral representations of -functions Bonella and Coker (2005a); Causo et al. (2005b). Thus, the linearized approximation for the evolution of the reduced density operator becomes

(8) | |||||

where

(9) |

and

(10) |

This result indicates that the time evolution of the density operator matrix elements can be approximated by first sampling the initial environment phase space from the Wigner transform of the bath part of the thermal equilibrium density of the full system. In our low temperature calculations we approximate this equilibrium density by assuming that the environment experiences interactions with only the ground state of the vibrational subsystem. At higher temperatures excited vibrational states should be included in this initial sampling. Next, the product of -functions in Eq.(8) indicates that with in the linearized approximation the time dependence of the density matrix elements can be obtained by evolving classical trajectories with the sampled initial conditions and subject to the mean of the forces arising from the quantum states involved in the prepared superposition state. Finally Eq.(8) gives that the contributions from each sampled trajectory to the density are obtained by adding coherently phase factors computed along these trajectories. An approach based on the same approximations outlined above has been used in other work Bonella and Coker (2005a); Causo et al. (2005b) to compute various quantum time correlation functions exploring electron transport, and vibrational energy relaxation, for example. In fact, approaches like this which employ the Wigner transform Hillery et al. (1984) of the equilibrium density as the distribution of initial conditions and evolve the classical dynamics with a mean Hamiltonian have a long history in computing spectroscopic correlation functions, for example, Mukamel (1982); Shemetulskis and Loring (1992). The work of Egorov et al. Egorov et al. (1998) provides an important comparison of this type of approach with alternative classical and mixed quantum-classical methods in the context of computing model condensed phase vibronic spectra.

The rest of the paper describes the computational approach we adopt for sampling the Wigner density and the mean Hamiltionian dynamics calculations that we have conducted to explore the characteristics of the vibrational pure dephasing process in low temperatures crystals. In Fig. 1 we present our dephasing results for I in solid krypton and compare with experiments in order to demonstrate the sensitivity of these type of calculations to the potential field used in the underlying dynamics. The calculation results presented in this figure both use the FK-Wigner initial condition sampling approach detailed in the next section. The different sets of calculated dephasing rates presented here, however, were obtained using different ways of propagating the classical trajectories. The solid curve presents results obtained using the average state force (ASF) to drive the environmental dynamics as suggested by Eq.(8) and discussed above. An alternative approach is to evolve the environmental dynamics using forces arising from the ground vibrational state potential onlyMukamel (1982); Egorov et al. (1998). The curve labeled GSF (Ground vibrational State Force) in Fig. 1 is obtained in this way. We see that both sets of calculated dephasing rates are generally faster than the experimentally observed results. These discrepancies with experiment probably reflect inaccuracy in the model interactions employed in these calculations which we outline later. However, we generally find that the dephasing rates computed using the average surface force are in better agreement with experiment than results computed using the ground vibrational state force. The two calculated dephasing rate curves should coalesce at low vibrational state number where the average potential from states 0 and is essentially that of the ground state. As we increase the average force experienced by the environment in our ASF calculations becomes more and more different from the ground state force employed in the GSF calculations and the deviation between the two curves is apparent. As discussed above, according to Eq. (8), the dephasing rates are obtained by averaging the phase factor in the energy gap between the two relevant vibrational states. If the forces governing the fluctuations in the energy gap arise from both of the states we expect that each state energy will fluctuate in a similar way leading to slower dephasing. If, on the other hand, the solvent responds to only the ground state force, as in the GSF calculations, the excited state energy may fluctuated more strongly as it does not play a role in influencing the environmental dynamics. This explains why the dephasing rates computed with the average force dynamics are slower than those obtained with the ground state force.

The main goal of our work is to study the effect of the distribution of environmental states on vibrational pure dephasing, making particular contact with low temperature experiments. Under these conditions quantum dispersion and tunneling of the environmental degrees of freedom may play a significant role. Thus we will employ various approximations to the Wigner transform of the initial density and compute their effect of the linearized approximation to the dephasing dynamics given in Eq.(8). General methods for computing the Wigner transform of the Boltzmann operator are as yet unavailable. In our work we thus employ an approximation to this operator. First we assume that the temperature is sufficiently high, and that particle masses are large, so that identical particle statistics complications can be ignored. Next, as detailed below, we employ an approach pioneered by Feynman and Kleinert that approximates the many body Boltzmann operator assuming a locally quadratic form for the interactions. This approximate approach has recently been implemented in various condensed phase applications studying quantum effects on transport and spectroscopy by Poulsen and co-workers Poulsen et al. (2003); Poulsen and Nyman (2004); Poulsen et al. (2005). Our study will compare vibrational dephasing results obtained using this approximate Boltzmann operator which can incorporate some quantum initial distribution effects, with classical initial condition sampling techniques which ignore quantum dispersion and tunneling. Our results from all these various calculations will be compared with available experimental results enabling a detailed understanding of the reliability of the different approximations underlying these calculations.

### ii.2 Sampling the quantum initial density for the environment

In this section we review the method developed by Poulsen and co-workers Poulsen et al. (2003); Poulsen and Nyman (2004); Poulsen et al. (2005); Cao and Voth (1993, 1994a, 1994b, 1994c, 1994d, 1994e); Cao and Berne (1990) that adapts the Feynman-Kleinert variational approach for obtaining an optimal local harmonic approximation to the Boltzmann operator which can be easily Wigner transformed. In our application to vibrational dephasing we suppose that the quantum vibration is prepared initially in some superposition state while the environmental degrees of freedom are unaffected by the vibrational excitation and are initially in thermal equilibrium with the vibrational degrees of freedom. As outlined in the previous section this quantity plays a central role as the phase space distribution function that must sampled to provide trajectory initial conditions in linearized approximations for the density matrix dynamics. We will explore the effects of the quantum nature of this initial equilibrium distribution on vibrational pure dephasing of the excited coherence. The development below in terms of the general Boltzmann operator is easily adapted to the case of the environmental degrees of freedom being in equilibrium with the quantum subsystem by supposing that the vibrational degree of freedom acts to provide an external force on the environment. As mentioned above we suppose that the experiments of interest can be interpreted as though the environment was initially in equilibrium with the ground state of the quantum subsystem so the Hamiltonian in the Boltzmann operator below describes interactions of the environment with the external field due to the ground state quantum subsystem.

The principle assumption underlying the Poulsen et al. Poulsen et al. (2003); Poulsen and Nyman (2004); Poulsen et al. (2005) approach to generating an approximate Wigner transform of the Boltzmann operator, is that the off-diagonal elements of the thermal density matrix that enter the Wigner transform expression

(11) |

can be just as well represented as the diagonal elements that appear in the trace expression for the partition function. Accepting this proposal the approach proceeds as follows: First the Feynman-Kleinert method for computing the partition function is employed to obtain a local harmonic approximation to the Euclidian action that is variationally optimized for computing the trace of the density matrix. Next we assume that the off-diagonal elements can be approximated using the optimized local harmonic approximation and this form is employed to compute the Wigner transformation. With this approach a local harmonic approximation is never directly made to the potential, rather an optimal local harmonic approximation to the Euclidean time action is found variationally using the partition function which involves integration over all space and thus contains global information. This local harmonic fitting procedure, however, is not influenced by any off-diagonal information. So it is unclear why it should provide a reliable approximation for generating the Wigner phase space density. One of our goals is to explore this question. In section II.3 we will directly compare density matrix elements for a simple model computed with in this approximation with exact results to explore the reliability of the off-diagonal density obtained with this and other approximate methods.

As outlined above the derivation of the approach starts by following Feynman and Kleinert and considering the partition function, , at inverse temperature, , for the environmental variables, , written as a path integral over cyclic paths i.e.

(12) |

Where, , means integrate over all points in the path, , including the common starting and finishing points, and the Euclidian action is . Here, , is the mass tensor, and , is the inter-particle interaction potential including the effects of external fields (such as, for example, coming from interactions with the vibrator in its ground state). In mass weighted cartesian coordinates, , we can write

(13) |

Since the paths of interest for are real and periodic on their Euclidean time dependence can be written in terms of a Fourier series i.e.

(14) |

where are the Matsubara frequencies, the real zero frequency Fourier coefficient, , is the path centroid, and all other Fourier coefficients, , are in general complex.

If the configuration of the system is of dimension then the cyclic path space integration can be mapped onto an integration over Fourier coefficient vectors of length according to Feynman and Kleinert (1986); Janke and Kleinert (1987); Roepstorff (1994); Kleinert (2004)

(15) |

and the partition function can be shown to be

(16) | |||||

Following Feynman and Kleinert Feynman and Kleinert (1986); Janke and Kleinert (1987) the effective classical potential, , is defined as

(17) | |||||

and the quantum partition function is cast in the suggestive classical form

(18) |

identifying motion of the path centroid over the effective potential as a classical-like way to compute quantum properties using classical statistical mechanics ideas Feynman and Kleinert (1986); Janke and Kleinert (1987); Feynman (1998). Since , at high temperature , the Gaussian factors in the integrand of Eq.(17), for example, will be dominated by small values of , so under these conditions will fluctuate only slightly from and, as Feynman and Kleinert suggested Feynman and Kleinert (1986), the effective potential might be computed using a local harmonic approximation to the full potential in the region around each path centroid location. They used the Gibbs-Bogoliobov-Feynman variational principle to determine the optimal local harmonic approximation at each centroid position.

This approach starts from the exact expression for the partition function obtained by multiplying and dividing by some trial partition function associated with some, yet to be specified, locally harmonic trial action, , thus

(19) | |||||

Due to the concavity of the exponential function, , so we replace the average of the exponential by the exponential of a more simply computed average and obtain the following variational result

(20) |

If we choose to be the Euclidian time action associated with motion in some trial potential , then

(21) |

Further since

(22) |

and the value of the cyclic path integral should be independent of the value of the time at which the function is evaluated since all should be equivalent under the trace we can write Feynman and Hibbs (1965), for example,

(23) |

To proceed we follow Feynman and Kleinert and suppose that the trial action is quadratic in displacements about the path centroid location thus

The term linear in displacement, , multiplying a derivative vector evaluated at vanishes identically when integrated over due to the cyclic nature of the paths on the interval (Eq.(14)), so with this quadratic form only the offset term, , and the local curvature matrix, , need to be determined variationally.

Suppose the matrix, , diagonalizes the curvature matrix giving a diagonal matrix of frequencies, , according to . Transforming to normal mode vectors, , we find that the trial partition function is obtained as

(25) | |||||

The Gaussian integrals over the components of the can be performed analytically provided the frequencies have appropriate values (see detailed discussion below on the restrictions that this places on the matrix ) and, using the fact that

(26) |

the trial partition function becomes

(27) |

so in analogy to Eq.(17) we can write

(28) |

The other quantities that need to be computed to apply the variational result in Eq.(20) according to Eqs.(21) - (23), are averages of the potentials over the trial density. In Appendix A we outline the computation of these quantities obtaining the following general result written in terms of different Gaussian smeared potential forms:

(29) | |||||

where is the Gaussian smeared full potential

(30) |

which can also be computed in Fourier space according to Eq.(A). The Gaussian smeared local harmonic approximate trial potential is

(31) |

which can be used in Eq.(29) to compute . In these results the smearing width matrix is

(32) |

with

(33) |

Thus we find the trial averaged action difference appearing in the variational expression, Eq.(20), can be written as

(34) |

Using the above result we follow Feynman and Kleinert Feynman and Kleinert (1986) and optimize the right hand side of Eq.(20), ,with respect to the variational parameter functions and . The details of this variational calculation are given in Appendix B and it gives the following relationships between the different optimized parameters:

(35) |

Here the optimal curvature matrix is found from the gaussian smeared result

(36) |

Where is the mass weighted Hessian

(37) |

As outlined above the approach requires the computation of Gaussian smeared averages of various functions. The evaluation of such multidimensional integrals is in general complicated and could be performed with a MC procedure. In the application described here, however, the interaction potential functions being smeared are pair decomposable and, as outlined in Appendix C, can be fit to Gaussian forms and the smearing integrals performed analytically.

The structure of the above equations suggests the following iterative procedure for computing the smearing matrix : (1) Choose an initial guess for the matrix. We use the previously converged smearing matrix for the last accepted centroid configuration. (2) Compute the full matrix at the new centroid configuration using Eq.(36) (or Eq.(84)). This can be done efficiently employing the symmetry of the local curvature matrix. (3) Diagonalize to obtain the local harmonic frequencies and the eigenvector matrix as defined above Eq.(25). Finally, (4), use these in Eqs.(33) (or (75)) and (32) to obtain a new estimate of the smearing matrix, , then return to step (2) and iterate this process till convergence.

Using the converged set of parameters, we compute the various Gaussian smeared quantities (using Eq.(30) for example) and obtain the variationally optimized value for the energy offset parameter as in Eq.(35) at the new centroid configuration.

The final steps in the approach involve writing the density operator approximately as the integral over centroid locations of the variationally optimized local harmonic form around each centroid position and then analytically Wigner transforming this locally harmonic result. Detail of these calculations are given in Appendix D. The approximate density has the form

(38) | |||||

Here we define the quantity and Wigner transformation yields

The integral in Eq.(II.2) is easily computed by importance sampled Monte Carlo. Thus a Cartesian centroid configuration, , is sampled. The procedure outlined above Eq.(38) is employed to iterate the width, , and curvature, , matrices to convergence. As indicated, this involves diagonalizing at each iteration to obtain the local harmonic frequencies and normal modes. The converged value of computed from Eq.(28) is used to accept or reject the sampled centroid configuration. Next, to sample the phase space distribution associated with the approximate Wigner transformed Boltzmann operator, a set of normal mode displacements with components are sampled from the component Gaussian distributions with variances and the set of points in Cartesian space is generated to provide a sampling of initial configurations around the point . Similarly, Gaussian random numbers are sampled with variance to provide a set of normal mode momentum vectors which are transformed to cartesian initial momenta according to . The numerical factor in Eq.(II.2) provides an automatic centroid configuration dependent normalization for the sampling of these Gaussian distributions, and so each sample phase space configuration generated in this way carries unit weight.

### ii.3 Comparison of Initial Condition Sampling Methods for Model System

To test the reliability of our implementation of the Feynman-Kleinert approach for sampling the Wigner distribution outlined in the previous subsection we have applied it to several simple exactly solvable one dimensional model systems. In this section we present these results and compare with exact calculations and with another approximate approach due to Shi and Geva (SG) Shi and Geva (2003a, b, c). These authors proceed by multiplying and dividing the Wigner transform expression of the Boltzmann operator by the diagonal elements of the thermal density operator, thus

(40) | |||||

They make the Local Harmonic Approximation (LHA) to the potential, expanding to quadratic order about the position . Note this approximation is made only in the ratio of off-diagonal to diagonal density matrix elements in the integrand, however, the full anharmonic dependence of prefactor diagonal element is included through path integral calculations. With this local harmonic form, the Gaussian integrals can be performed analytically yielding the following result for the case of one dimension:

(41) | |||||

Where depends on the local curvature of the potential about the point and has the form with .

For the purpose of comparison with exact results we computed the density matrix for our 1D model using the numerical matrix multiplication (NMM) approach Klemm and Storer (1973); Thirumalai et al. (1983). It is most convenient to make direct comparison with the full density matrix. The approximate Wigner densities can be inverse Wigner transformed to give approximate coordinate space density matrices using to the following result:

(42) |

Thus the SG approximation to the density matrix becomes

(43) |

with . The corresponding FK approximation to the density matrix takes the form: