space flow H.López, A.Donoso ∗

A technique for simulating quantum dynamics in phase space i discussed. It makes use of ensembles of classical trajectories to approximate the distribution functions and their derivatives by implementing Adaptive Kernel Density Estimation. It is found to improve the accuracy and stability of the simulations compared to more conventional particle methods. Formulation of the method in higher dimensions is straightforward. Complete analytical description of physical systems exist only for the most simple ideal cases. The vast majority of physical problems require the solution of approximate equations most fre- quently by using numerical calculations. To this end, much effort has been devoted to the devel- opment of suitable techniques to be used in computer simulations. The increasing computational power makes possible the description of more complicated, that is, less ideal systems allowing for more precise predictions of their properties as well as for a better understanding of the underlying phenomena. Among the wealth of numerical approaches to solve dynamical equations by means of discretization, particle methods have gained attention in the last few decades as a conceptu- ally appealing alternative. Conventional numerical solvers rely on a regular discretization of the independent variables to facilitate the evaluation of derivatives to a given order of accuracy. As a consequence, a mesh needs to be defined as the base of the simulations, parts of which remain empty at any given simulation time even with the most sophisticated mesh schemes. The cen- tral idea of particle methods consists of representing functions by sets of particles conveniently located and weighted such that the function is appropriately sampled by the particles. Evolution of the functions is then accomplished by an equivalent evolution of the particles. This meshless or Lagrangian type of discretization has some intrinsic computational efficiency advantages because there is no wasted space. The advantage, however, comes at the expense of having to evaluate derivatives of functions on unstructured arrays of sampling points. The issues about convergence and stability of the resulting algorithms have been the subject of active and fruitfull research. Currently there are some well established techniques to solve dynamical equations with particle methods. Our interest is to provide an alternative ingredient to further improve the performance of numerical solutions of dynamical problems of physical relevance. In the present work we discuss, in the context of particle methods, the implementation of a technique for representing functions and derivatives with adaptive parameters. Adaptiveness require the parameters to be automatically ad- justed to optimal values according to the location of the particles. As shown below, the proposed scheme is conceptually simple and has been tested in a number of different situations. In some simple test models the present work focuses on quantum dynamics, previously treated with more


Introduction
Complete analytical description of physical systems exist only for the most simple ideal cases.The vast majority of physical problems require the solution of approximate equations most frequently by using numerical calculations.To this end, much effort has been devoted to the development of suitable techniques to be used in computer simulations.The increasing computational power makes possible the description of more complicated, that is, less ideal systems allowing for more precise predictions of their properties as well as for a better understanding of the underlying phenomena.Among the wealth of numerical approaches to solve dynamical equations by means of discretization, particle methods have gained attention in the last few decades as a conceptually appealing alternative.Conventional numerical solvers rely on a regular discretization of the independent variables to facilitate the evaluation of derivatives to a given order of accuracy.As a consequence, a mesh needs to be defined as the base of the simulations, parts of which remain empty at any given simulation time even with the most sophisticated mesh schemes.The central idea of particle methods consists of representing functions by sets of particles conveniently located and weighted such that the function is appropriately sampled by the particles.Evolution of the functions is then accomplished by an equivalent evolution of the particles.This meshless or Lagrangian type of discretization has some intrinsic computational efficiency advantages because there is no wasted space.The advantage, however, comes at the expense of having to evaluate derivatives of functions on unstructured arrays of sampling points.The issues about convergence and stability of the resulting algorithms have been the subject of active and fruitfull research.Currently there are some well established techniques to solve dynamical equations with particle methods.Our interest is to provide an alternative ingredient to further improve the performance of numerical solutions of dynamical problems of physical relevance.In the present work we discuss, in the context of particle methods, the implementation of a technique for representing functions and derivatives with adaptive parameters.Adaptiveness require the parameters to be automatically adjusted to optimal values according to the location of the particles.As shown below, the proposed scheme is conceptually simple and has been tested in a number of different situations.In some simple test models the present work focuses on quantum dynamics, previously treated with more conventional particle methods.Dynamical equations correspond to a phase space representation of quantum mechanics, in terms of distribution functions in phase space.Derivation of the equations and a more detailed description of the test models can be found in the references.
We first give a brief review of the implementation of particle methods and present the adaptive kernel approach.Applications to the simulation of quantum phase space dynamics are presented hereinbelow.Finally, as part of the conclusions, we briefly discuss further applications of the method in a different context, where it has been used to simulate compressible and viscous fluid dynamics.

Adaptive particle methods
Given an equation of motion for a distribution function, a conceptually appealing approach to perform numerical calculations is what has become known as particle methods.The idea is quite simple and intuitive.Distribution functions are represented by the density of a set of points (trajectories) in phase space.Then, the time evolution of the system is accomplished by propagating the ensemble of points under their corresponding equations of motion.To evaluate the density, each point contributes as an extended particle by means of a kernel function centered at its location with an appropriate radius of influence.Rigorous proof of convergence and stability can be found, for instance, in reference [1].This idea is the basis of a variety of methods, such as Dissipative Particle Dynamics (DPD) [2] and Smoothed Particle Hydrodynamics (SPH) [3], which have found fruitfull applications in many fields.A crucial step to the implementation of particle methods is the evaluation of the densities and their derivatives.In general, particle methods tend to become unstable in regions of low density (for instance, the tails) due to the lack of particles to approximate the required functions.This is especially important for higher derivatives of density where the numerical errors become amplified.Conventional approaches make use of a fixed radius of influence for all kernels.Instabilities are prevented by using an increasing number of particles, at the expense of raising the computational burden.We have recently implemented an alternative technique that significantly improves the performance of particle methods.It makes use of an idea from mathematical statistics called Adaptive Kernel Density Estimation [4,5] originally proposed for the nonparametric estimation of probability distributions given a set of samples of random variables.Substantial improvement on the approximation of functions and derivatives are obtained by allowing the kernel parameters to dynamically adapt, both globally and locally, to the given data set.Of course, the idea of adaptiveness on particle methods is not new, since there have been attempts for most of them.The scheme we discuss here has the advantages of being quite simple to implement, generalizations to higher dimensions are straightforward and does not increase computational costs too much.
The estimation problem can be stated as follows.For simplicity we assume a two dimensional space, although higher dimensions are treated without complications.Given a set of N points in phase space with locations {q i , p i }, i = 1, . . ., N , we construct an approximate distribution function using Adaptive Kernel Density Estimation [4,5], as where K i (x, y) is a Gaussian kernel function capable of adapting both globally and locally.Globally, the kernels are allowed to be stretched and oriented according to the general shape of the data.Moreover, the kernel widths are locally adjusted, becoming uniformly wider in regions of low density and narrower in regions of high density.The general shape of the data is described by its sample covariance matrix S. In two dimensions it has the elements , where µ is the corresponding sample average.The functional form of the kernels with fixed widths is where D = det S and h = N −1/6 is the optimal window width for two dimensions and Gaussian functional kernels.We have found this form of the kernels to give good results in calculating the densities and first derivatives.Second and higher derivatives, however, require extra effort, as in most of the particle methods.The locally adaptive part of the estimation is accomplished in three steps.First, a pilot density ρ is obtained by equation (1) using fixed width kernels, the same for all particles, as given by equation (2).The second step is to calculate the local bandwidth factors , where ρi is the pilot density evaluated at the location of particle i, g is the geometric mean of the ρi , satisfying log g = N −1 i log ρi , and α = 0.1 is a sensitivity parameter, adjusted by experimentation.The final step is then to include the local bandwidth factors in the kernels, giving the final expression Derivatives of the density function are obtained by directly differentiating the kernel.Thus, the density function inherits the differentiability properties of the kernel.Although being slightly more demanding computationally, having to evaluate the density twice, this refinement substantially improves the quality of the approximation and is essential when derivatives higher than the first one are present in the equations of motion.We found that the adaptive properties of the present scheme lead to accurate representations of the distribution functions and numerically stable propagation in time.There is also the possibility of speeding up the calculations by skipping the adaptive part for a few time steps, depending on the particular problem.Equations of motion for the points should now be derived.Let us assume that the trace of the distribution function is preserved by time evolution.This is the case for quantum phase space flux as well as for any process where mass is conserved, including compressible fluid flow.The case where sources and sinks are present would require a slight modification of the present scheme.Then the equation of motion for the distribution function can be written as a conservation law where f 1 and f 2 are functionals of q, p, ρ(q, p) and its derivatives.Trace preservation implies the existence of a continuity equation ρ t + ∇ • (vρ) = 0, where ∇ is the gradient in phase space coordinates and v = ( q, ṗ) is a velocity field, which needs to be solved.One possible solution comes from identifying the components of the velocity field with the functionals f 1 and f 2 , that is, q = f 1 and ṗ = f 2 .The equations of motion for the points are then obtained by evaluating the velocity field at the particle locations, resulting in a set of coupled ordinary differential equations.Here we note that this approach reports two important advantages.Firstly, the original partial differential equation has been reduced to a set of ordinary equations.Secondly, the differential equations have been reduced one order, due to the use of the continuity equation.
As an example, let us consider the one dimensional diffusion equation, with D a diffusion constant.It can be recast as a conservation law as ρ t + (−Dρ x ) x = 0.A continuity equation for this case is ρ t + (vρ) x = 0, from which we can identify v = −Dρ x /ρ as the velocity field.If the density is constructed, as above, from an ensemble of points, then the equations of motion for the points are given by ẋi with i = 1, . . ., N and where K is the derivative of the kernel.These N coupled equations contain only first derivatives of the kernel, whereas the original partial differential equation is second order.We have solved the diffusion equation in several cases including the phase space Fokker-Planck equation in previous work [6].We now present some applications of the present method.

Quantum phase space flow
Phase space hydrodynamic formulation of quantum mechanics was presented shortly after the early developments of the main theory by the pioneering work of Madelung [7], de Broglie [8] and Wigner [9] as an attempt to rescue the classical perspective from its extravagant predictions.It is in this context where quantum mechanics finds a more intuitive view, as a flow of a compressible fluid in phase space.The celebrated Wigner distribution function has been extraordinarily successfull in describing the similarities and differences between classical and quantum mechanics.It represents the phase space analog of the density matrix, carrying all the physical information about the system.It is used to find expectation values of observables by means of a trace over phase space variables, although its interpretation as a probability distribution is compromised by its lack of positiveness in general.Since the early works of Wigner, alternative phase space distribution functions have been developed [10][11][12], each of them having a different equation of motion.The Husimi function, for instance, has the advantage of being everywhere positive for all times and smoother than the Wigner function, although its equation of motion becomes relatively more complicated.
In previous work [13] we have presented a trajectory-based method for simulating quantum hydrodynamic flow in phase space for a particle escaping from a metastable potential well by tunneling.This test system was selected for being the simplest nontrivial example manifestly exhibiting quantum effects.Exact calculations can be performed by propagating a wavefunction using the Schrödinger equation.For this purpose, we used the Fast Fourier Transform method of Kosloff [14].The problem was formulated in the Wigner representation and a novel technique was presented for evaluating the derivatives of the density.A local Gaussian anzatz was used to fit the required functions to the data, giving qualitatively good results on capturing the tunnel effect.Here we present new results on this system simulated with adaptive kernels.The problem is formulated here in the Husimi representation to ensure positiveness of the distribution function.Moreover, the smoother character of the Husimi function allows for a more accurate representation of the density and its derivatives.This approach leads to a significant improvement of the results compared to more conventional, non-adaptive methods, yielding a nearly quantitative agreement with exact calculations.
We use atomic units ( = 1) throughout.The system [13] consists of a particle of mass m = 2000 (the unit of mass in this system of units is m e , the electron mass) trapped in a metastable potential well V (q) = (mω 2 /2)q 2 − Bq 3 with ω = 0.01 and B = 0.2981.It has a local minimum at q = 0 and a potential barrier at q ‡ = 0.6709 of height V ‡ = 0.015.Under these conditions the dynamics is highly quantum mechanical.The equation of motion for the distribution function is given in the Husimi representation by The initial state of the system corresponds to a Gaussian wavepacket centered near the bottom of the well at q 0 = −0.3 and has the form where σ q = 1/(2mω) and σ p = mω/2 satisfy the minimum uncertainty principle.The fraction of particles with q > q ‡ is calculated as a function of time and compared with the quantum reaction probability obtained from the square of the wavefunction |ψ(q)| 2 integrating from q ‡ to infinity.The results are shown in figure 1.For t = 0 no particle has yet crossed the barrier.An abrupt rise on the curves near t = 500 signals the passage of the fraction of particles having energy higher than the barrier.This is also observed classically.For later times, particles continue escaping by tunneling making the reaction probability grow.More detailed inspections on the results of the simulation show that some particles perform multiple recrossings.Energy is conserved on average but is continuously redistributed among members of the ensemble.This is how some of the particles gain enough energy to surmount the potential barrier.In some cases particles that have crossed the barrier lose energy by transferring it back to the rest of the ensemble thus ending up with the energy lower than the barrier.In this perspective, particles do not mysteriously cross through the barrier.Instead, they act collectively to pass over the barrier by borrowing energy from the ensemble.The energy of a given trajectory changes according to the quantity dH/dt evaluated at the corresponding location.This quantity does not vanish in general for each of the trajectories, as can be verified by substituting the expression for the velocity field on it.On the other hand, the ensemble average dH/dt does vanish for all t for bounded distribution functions.Thus, energy is manifestly conserved only on average.It is also easy to verify that the present scheme obeys Ehrenfest theorem as well, that is, the ensemble average relation ṗ = − V (q) is maintained for all times.It means that, on average, the ensemble moves classically.This property is not shared by other approaches proposed so far, such as the Wigner trajectories methods.The curves in figure 1 show the results of the method proposed here, using the Husimi representation and adaptive estimation, compared to exact results as well as to simulations performed using purely classical motion.The interacting trajectories correctly capture the quantum tunnel effect.Compared to the method previously reported in reference [13] the present scheme yields more accurate results, reducing the overall error by roughly a half (in the root mean square sense).

Conclusions
Phase space formulation of quantum dynamics provides an intuitive picture of the theory allowing for a more direct comparison with classical perspectives.By means of ensembles of interacting classical trajectories, quantum phenomena can be correctly captured in molecular dynamicslike simulations.For this to be done in stable and reliable manner, a crucial ingredient is a faithfull representation of the distribution functions (and derivatives) including regions of lower density.Here we have presented an alternative scheme to perform such approximations.Adaptive kernel density estimation methods provide a simple, yet accurate, means to construct distribution functions from a set of points in phase space giving reasonably well represented low density zones.The formulation goes to higher dimensions unchanged.We are currently pursuing simulations in four and six dimensional phase spaces (corresponding to two and three spatial dimensions) which will be reported elsewhere.
The adaptive kernel density estimation technique is not limited to phase space dynamics.It has also been incorporated into a standard Smoothed Particle Hydrodynamics (SPH) method to simulate conventional fluid dynamics yielding a substantial improvement on the stability of the calculations.In recent work [15], we have reported a variant of SPH method that allows for the simulation of compressible fluid with strong shocks and rarefaction waves in one and two dimensions.The known instabilities of the conventional scheme (the wall heating phenomena) virtually dissapear when adaptive kernels, as presented here, are implemented.The method was found to outperform other approaches, such as Godunov-type methods and SPH formulations based on Riemann solutions.Convergence tests have been performed and reported in reference [16].Another improvement on the SPH method is the elimination of the tensile instability suffered by the regular formulation.Non-adaptive versions of SPH are known to produce such spurious phenomenon on the interface of liquid drop simulations [17].It appears as an artificial clustering of the particles near the surface of the drop.The adaptive kernel method performs only a minimum necessary smoothing to the data making the tensile instability disappear.
We are currently pursuing further refinements to the technique by prescribing kernel functions that can freely adapt not to the general shape of the data but to the close surroundings of each particle.This would allow for the simulation of more general problems such as multiphase fluids and viscous drop oscillations with separation and coalescence.In the quantum phase space dynamics area, we are interested in the simulation of quantum effects of systems in condensed phase.The combined effects of tunneling and quantum dissipation are being simulated in low temperature regimes by solving a Caldeira-Leggett type of equations with adaptive kernel particle methods.The results will be reported elsewhere.

Figure 1 .
Figure 1.Reaction probability vs. time.Solid line represents exact results.Fine dotted line shows classical results.Dashed line presents the results of this work performed in the Husimi representation.Simulations correspond to the initial placement of the center of the wavepacket at q0 = −0.3(in atomic units).