text:
# Packing of Compressible Granular Materials

## Abstract

3D Computer simulations and experiments are employed to study random packings of compressible spherical grains under external confining stress. Of particular interest is the rigid ball limit, which we describe as a continuous transition in which the applied stress vanishes as $`(\varphi \varphi _c)^\beta `$, where $`\varphi `$ is the (solid phase) volume density. This transition coincides with the onset of shear rigidity. The value of $`\varphi _c`$ depends, for example, on whether the grains interact via only normal forces (giving rise to random close packings) or by a combination of normal and friction generated transverse forces (producing random loose packings). In both cases, near the transition, the system’s response is controlled by localized force chains. As the stress increases, we characterize the system’s evolution in terms of (1) the participation number, (2) the average force distribution, and (3) visualization techniques.

PACS: 81.06.Rm

Dense packings of spherical particles are an important starting point for the study of physical systems as diverse as simple liquids, metallic glasses, colloidal suspensions, biological systems, and granular matter . In the case of liquids and glasses, finite temperature molecular dynamics (MD) studies of hard sphere models have been particularly important. Here one finds a first order liquid-solid phase transition as the solid phase volume fraction, $`\varphi `$, increases. Above the freezing point, a metastable disordered state can persist until $`\varphi \varphi _{\text{RCP}}^{}`$ , where $`\varphi _{\text{RCP}}`$ is the density of random close packing (RCP)— the densest possible random packing of hard spheres.

This Letter is concerned with the non-linear elastic properties of granular packings. Unlike glasses and amorphous solids, this is a zero temperature system in which the interparticle forces are both non-linear, and path (i.e., history) dependent. \[Because these forces are purely repulsive, mechanical stability is achieved only by imposing external stress.\] The structure of packing depends in detail on the forces acting between the grains during rearrangement of grains; indeed, different rearrangement protocols can lead to either RCP or random loose packed (RLP) systems.

In the conventional continuum approach to this problem, the granular material is treated as an elasto-plastic medium . However, this approach has been challenged by recent authors  who argue that granular packings represents a new kind of fragile matter and that more exotic methods, e.g., the fixed principal axis ansatz, are required to describe their internal stress distributions. These new continuum methods are complemented by microscopic studies based on either contact dynamics simulations of rigid spheres or statistical models, such as the q-model, which makes no attempt to take account of the character of the inter-grain forces .

In our view, a proper description of the stress state in granular systems must take account of the fact that the individual grains are deformable. We report here on a 3D study of deformable spheres interacting via Hertz-Mindlin contact forces. Our simulations cover four decades in the applied pressure and allow us to understand the regimes in which the different theoretical approaches described above are valid. Since the grains in our simulations are deformable, the volume fraction can be increased above the hard sphere limit and we are able to study the approach to the RCP and RLP states from this realistic perspective. Within this framework, the rigid grain limit is described as a continuous phase transition where the order parameter is the applied stress, $`\sigma `$, which vanishes continuously as $`(\varphi \varphi _c)^\beta `$. Here $`\varphi _c`$ is the critical volume density, and $`\beta `$ is the corresponding critical exponent. We emphasize that the fragile state corresponding to rigid grains is reached by looking at the limit $`\varphi \varphi _c^+`$ from above.

Of particular importance is the fact that $`\varphi _c`$ depends on the type of interaction between the grains. If the grains interact via normal forces only , they slide and rotate freely mimicking the rearrangements of grains during shaking in experiments . We then obtain the RCP value $`\varphi _c=0.634(4)(\varphi _{\text{RCP}})`$. By contrast, if the grains interact by combined normal and friction generated transverse forces, we get RLP  at the critical point with $`\varphi _c=0.6284(2)<\varphi _{\text{RCP}}`$. The power-law exponents characterizing the approach to $`\varphi _c`$ are not universal and depend on the strength of friction generated shear forces.

Our results indicate that the transitions at both RCP or RLP are driven by localized force chains. Near the critical density there is a percolative fragile structure which we characterize by the participation number (which measures localization of force chains), the probability distribution of forces, and also by visualization techniques. A subset of our results are experimentally verified using carbon paper measurements to study force distributions in the granular assembly. We also consider in some detail the relationship between our work and recent experiments in 2D Couette geometries .

Numerical Simulations: To better understand the behavior of real granular materials, we perform granular dynamics simulations of unconsolidated packings of deformable spherical glass beads using the Discrete Element Method developed by Cundall and Strack . Unlike previous work on rigid grains, we work with a system of deformable elastic grains interacting via normal and tangential Hertz-Mindlin forces plus viscous dissipative forces . The grains have shear modulus 29 GPa, Poisson’s ratio 0.2 and radius 0.1 mm.

Our simulations employ periodic boundary conditions and begin with a gas of 10000 non-overlapping spheres located at random positions in a cube 4 mm on a side. Generating a mechanically stable packing is not a trivial task . At the outset, a series of strain-controlled isotropic compressions and expansions are applied until a volume fraction slightly below the critical density. At this point the system is at zero pressure and zero coordination number. We then compress along the $`z`$ direction, until the system equilibrates at a desired vertical stress $`\sigma `$ and a non-zero average coordination number $`Z`$.

Figure 1a shows the behavior of the stress as a function of the volume fraction. We find that the pressure vanishes at a critical $`\varphi _c=0.6284(2)`$. Although we cannot rule out a discontinuity in the pressure at $`\varphi _c`$— as we could expect for a system of hard spheres— our results indicates that the transition is continuous and the behavior of the pressure can be fitted to a power law form

$$\sigma (\varphi \varphi _c)^\beta ,$$
(1)
where $`\beta =1.6(2)`$. Our 3D results contrast with recent experiments of slowly sheared grains in 2D Couette geometries  where a faster than exponential approach to $`\varphi _c`$ was found, while they agree qualitative with similar continuous transition found in compressed emulsions and foams .

Figure 1b shows the behavior of the mean coordination number, $`Z`$, as a function of $`\varphi `$. We find

$$ZZ_c(\varphi \varphi _c)^\theta ,$$
(2)
where $`Z_c=4`$ is a minimal coordination number, and $`\theta =0.29(5)`$ is a critical exponent. At criticality the system is very loose and fragile with a very low coordination number. The value of $`Z_c`$ can be understood in term of constraint arguments as discussed in ; in the rigid ball limit, for a disordered system with both normal and transverse forces, we find $`Z_c=D+1=4`$ . As we compress the system more contacts are created, providing more constraints so that the forces become overdetermined.

We notice that $`\varphi _c`$ obtained for this system is considerably lower than the best estimated value at RCP , $`\varphi _{\text{RCP}}=0.6366(4)`$ obtained by Finney  using ball bearings. This latter value is obtained by carefully vibrating the system and letting the grains settle into the most compact packing. Numerically, this is achieved by allowing the grains reach the state of mechanical equilibrium interacting only via normal forces. By removing the transverse forces, grains can slide freely and find most compact packings than with transverse forces. Numerically we confirm this by equilibrating the system at zero transverse force. The critical packing fraction found in this way is $`\varphi _c=0.634(4)`$($`\varphi _{\text{RCP}}`$ within error bars). The stress behaves as in Eq. (1) but with a different exponent $`\beta =2.0(2)`$ (Fig. 1a). At the critical volume fraction the average coordination number is now $`Z_c=6`$ \[and $`\theta =0.94(5)`$, Fig. 1b\], which again can be understood using constraint arguments which would give a minimal coordination number equal to 2D for frictionless rigid balls .

We conclude that the value $`\varphi _c0.6288<\varphi _{\text{RCP}}`$ found with transverse forces corresponds to the RLP limit, experimentally achieved by pouring balls in a container but without allowing for further rearrangements . Experimentally, stable loose packings with $`\varphi `$ as low as 0.60 have been found . In our simulations, $`\varphi _c`$ lower than $`0.6288`$ can be obtained by increasing the strength of the tangential forces. This is in agreement with experiments of Scott and Kilgour  who found that the maximum packing density of spheres decreases with the surface roughness (friction) of the balls.

While previous studies characterized RCP’s and RLP’s by using radial distribution functions and Voronoi constructions , we take a different approach which allow us to compare our results directly with recent work on force transmissions in granular matter. Previous studies of granular media indicate that, for forces greater than the average value, the distribution of inter-grain contact forces is exponential . In addition, photo-elastic visualization experiments and simulations  show that contact forces are strongly localized along “force chains” which carry most of the applied stress. The existence of force chains and exponential force distributions are thought to be intimately related.

Here we analyze this scenario in the entire range of pressures: from the $`\varphi _c`$ limit and up. Figure 2a shows the force distribution obtained in the simulations with friction balls. At low stress, the distribution is exponential in agreement with previous experiments and models. When the system is compressed further, we find a gradual transition to a Gaussian force distribution. We observe a similar transition in our simulations involving frictionless grains under isotropic compression. This suggests that our results are generic, and do not depend, qualitatively, on the preparation history or on the existence of friction generated transverse forces between the grains.

Physically, we find that the transition from Gaussian to exponential force distribution is driven by the localization of force chains as the applied stress is decreased. In granular materials, with particles of similar size, localization is induced by the disorder of the packing arrangement. To quantify the degree of localization, we consider the participation number $`\mathrm{\Pi }`$:

$$\mathrm{\Pi }\left(M\underset{i=1}{\overset{M}{}}q_i^2\right)^1.$$
(3)
Here $`M`$ is the number of contacts between the spheres, $`Z=2M/N`$ is the average coordination number, and $`N`$ is the number of spheres. $`q_if_i/_{j=1}^Mf_j`$, where $`f_i`$ is the magnitude of the total force at every contact. From the definition (3), $`\mathrm{\Pi }=1`$ indicates a limiting state with a spatially homogeneous force distribution ($`q_i=1/M`$, $`i`$). On the other hand, in the limit of complete localization, $`\mathrm{\Pi }1/M0`$ and $`M\mathrm{}`$.

Figure 2c shows our results for $`\mathrm{\Pi }`$ versus $`\sigma `$. Clearly, the system is more localized at low stress than at high stress. Initially, the growth of $`\mathrm{\Pi }`$ is logarithmic, indicating a smooth delocalization transition. This behavior is seen up to $`\sigma `$ 2.1 MPa, after which the participation number saturates to a higher value:

$$\begin{array}{cccc}\hfill \mathrm{\Pi }(\sigma )& & \mathrm{log}(\sigma )\hfill & [\sigma <\text{2.1 MPa}]\hfill \\ \hfill \mathrm{\Pi }(\sigma )& & 0.62\hfill & [\sigma >\text{2.1 MPa}]\hfill \end{array}$$
(4)
This behavior suggests that, near the critical density, the forces are localized in force chains sparsely distributed in space. As the applied stress is increased, the force chains become more dense, and are thus distributed more homogeneously.

How might we expect the participation number to depend upon other system parameters when the forces are transmitted principally by force chains? In an idealized situation, the system has $`N_{FC}`$ force chains, each of which has $`N_z`$ spheres. Each sphere in a force chain has two major load bearing contacts, which loads must be approximately equal. In the lateral directions, roughly four weak contacts are required for stability. These contacts carry a fraction $`\alpha <1`$ of the major vertical load. All other contacts have $`f_i0`$. Under these assumptions,

$$\mathrm{\Pi }=\frac{2}{Z}\frac{(1+2\alpha )^2}{(1+2\alpha ^2)}\frac{N_{FC}N_z}{N}\frac{2}{Z}\frac{(1+2\alpha )^2}{(1+2\alpha ^2)}.$$
(5)
The last inequality becomes an equality iff all the balls are in force chains. From our simulations at large pressure $`\alpha 2/5`$, so at $`Z8`$, $`\mathrm{\Pi }0.62`$, which implies that the system has been completely homogenized. Although Eq. (5) is oversimplified, we believe that the change in slope in Fig. 2c is emblematic of the complete disappearance of well-separated chains.

The localization transition can be understood by studying the behavior of the forces during the loading of the sample. Clearly, visualizing forces in 3D systems is a complicated task. In order to exhibit the rigid structure from the system we visually examine all the forces larger than the average force; these carry most of the stress of the system. The forces smaller than the average are thought to act as an interstitial subset of the system providing stability to the buckling of force chains . We look for force chains by starting from a sphere at the top of the system, and following the path of maximum contact force at every grain. We look only for the paths which percolate, i.e., stress paths spanning the sample from the top to the bottom. In Fig. 3 we show the evolution of the force chains thus obtained for two extreme cases of confining stress. We clearly see localization at low confining stress: the force-bearing network is concentrated in a few percolating chains. At this point the grains are weakly deformed but still well connected. We expect a broad force distribution, as found in this and previous studies. As we compress further, new contacts are created and the density of force chains increases. This in turn gives rise to a more homogeneous spatial distribution of forces, which is consistent with the crossover to a narrow Gaussian distribution.

Experiments: Some of the predictions of our numerical study can be tested using standard carbon paper experiments , which have been used successfully in the past to study the force fluctuations in granular packings. A Plexiglas cylinder, 5 cm diameter and varying height (from 3 cm to 5 cm), is filled with spherical glass beads of diameter $`0.8\pm 0.05`$ mm. At the bottom of the container we place a carbon paper with white glossy paper underneath. We close the system with two pistons and we allow the top piston to slide freely in the vertical direction, while the bottom piston is held fixed to the cylinder. The system is compressed in the vertical direction with an InktronTM press and the beads at the bottom of the cylinder left marks on the glossy paper. We digitize this pattern and calculate the “darkness”  of every mark on the paper. To calibrate the relationship between the marks and the force, a known weight is placed on top of a single bead; for the forces of interest in this study (i.e., from $`0.05`$N to 6 N), there is a roughly linear relation between the darkness of the dot and the force on the bead.

We perform the experiment for different external forces, ranging from 2000 N to 9000 N, and different cylinder heights. The corresponding vertical stress, $`\sigma `$, at the bottom of the cylinder ranges between 100 KPa and 2.3 MPa (as measured from the darkness of the dots). The results of four different measurements are shown in Fig. 2b. For $`\sigma `$ smaller than $`750`$ KPa, the distribution of forces, $`f`$, at the bottom piston decays exponentially:

$$P(f)=f^1\mathrm{exp}[f/f],[\text{ }\sigma <\text{ 750 KPa}],$$
(6)
where $`f`$ is the average force. When the stress is increased above 750 KPa there is a gradual crossover to a Gaussian force distribution as we find in the simulations. For example, at 2.3 MPa we have

$$P(f)\mathrm{exp}\left[k^2\left(ff_o\right)^2\right],[\sigma =\text{ 2.3 MPa}].$$
(7)
where $`kf_o1`$, and therefore $`ff_o`$. Similar results have been found in 2D geometries .

Discussion:  In summary, using both numerical simulations and experiments, we have studied unconsolidated compressible granular media in a range of pressures spanning almost four decades. In the limit of weak compression, the stress vanishes continuously as $`(\varphi \varphi _c)^\beta `$, where $`\varphi _c`$ corresponds to RLP or RCP according to the existence or not of transverse forces between the grains, respectively. At criticality, the coordination number approaches a minimal value $`Z_c`$ (=4 for friction and 6 for frictionless grains) also as a power law. Our result $`Z_c=6`$ agrees with experimental analysis of Bernal packings for close contacts between spheres fixed by means of wax , and our own analysis of the Finney packings  using the actual sphere center coordinates of 8000 steel balls. However, no similar experimental study exists for RLP which could be able to confirm $`Z_c=4`$. A critical slowing down— the time to equilibrate the system increases near $`\varphi _c`$— and the emergency of shear rigidity (to be discussed elsewhere) is also found at criticality. The distribution of forces is found to decay exponentially. The system is dominated by a fragile network of relatively few force chains which span the system.

When the stress is increased away from $`\varphi _c`$ to the point that the number of contacts has significantly increased from its initial value $`Z_c`$ we find: (1) the distribution of forces crosses over to a Gaussian (2) the participation number increases, and then abruptly saturates and (3) the density of force chains increases to the point where it no longer makes sense to describe the system in those terms. Our simulations indicate that the crossover is associated with a loss of localization and the ensuing homogenization of the force-bearing stress paths.


format: text

text:
# Semiclassical description of multiphoton processes

## I Introduction

In the last two decades multiphoton processes have been studied intensively, experimentally as well as theoretically. The inherent time-dependent nature of an atomic or molecular excitation process induced by a short laser pulse renders a theoretical description problematic in two respects. Firstly, a full quantum calculation in three dimensions requires a large computational effort. For this reason, quantum calculations have been restricted to one active electron in most cases . Secondly, an intuitive understanding for an explicitly time dependent process seems to be notoriously difficult, exemplified by pertinent discussions about stabilization in intense laser fields . Many studies have been carried out to gain an intuitive understanding of the two most prominent strong field phenomena, namely High Harmonic Generation (HHG) and Above Threshold Ionization (ATI). In the well established early analytical formulation by Keldysh, Faisal and Reiss the atomic potential is treated as a perturbation for the motion of the electron in a strong laser field .

This picture is still used in more recent models, where the classical dynamics of the electron in the laser field is explicitly considered, e.g. in Corkum’s rescattering model which can explain the cutoff observed in HHG for linearly polarized laser light in one spatial dimension . The corresponding Hamiltonian reads

$$H=H_0+E_0f(t)x\mathrm{sin}(\omega _0t+\delta ),$$
(1)
where $`H_0=\frac{1}{2}p^2+V(x)`$ is the atomic Hamiltonian, $`f(t)`$ is the time-profile of the laser pulse with maximum amplitude $`E_0`$, and $`\omega _0`$ is the laser frequency. The interaction of the electron with the atom is specified by the potential $`V`$.

Lewenstein et al. extended Corkum’s rescattering idea to a quasiclassical model which contains one (relevant) bound state not influenced by the laser field on the one hand side and electrons which only feel the laser field on the other side . This simple model explains qualitatively the features of HHG well. The same is also true for an alternative model, where the electron is bound by a zero-range potential  However, the basic question if and to which extent these multiphoton processes can be understood semiclassically, i.e., by interference of classical trajectories alone, remains unanswered. It is astonishing that no direct semiclassical investigation of the Hamiltonian Eq. (1) has been performed while a number of classical as well as quantum calculations for Eq. (1) have been published. However, only recently, a semiclassical propagation method has been formulated which can be implemented with reasonable numerical effort. This is very important for the seemingly simple Hamiltonian Eq. (1) whose classical dynamics is mixed and in some phase space regions highly chaotic which requires efficient computation to achieve convergence. Equipped with these semiclassical tools we have studied multiphoton phenomena semiclassically in the frame work of Eq. (1). In comparison to the exact quantum solution we will work out those features of the intense field dynamics that can be understood in terms of interference of classical trajectories.

The plan of the paper is as follows. In section II we provide the tools for the calculation of a semiclassical, time-dependent wavefunction. In section III we discuss Above-Threshold-Ionization (ATI) and work out the classical quantities which structure semiclassically the relevant observables. In section IV we use this knowledge for the description of Higher-Harmonic-Generation (HHG). Section V concludes the paper with a comparison of HHG and ATI from a semiclassical perspective and a short summary.

## II Calculation of the semiclassical wave function

A (multi-dimensional) wave function $`\mathrm{\Psi }_\beta (𝐱,t)`$ can be expressed as

$$\mathrm{\Psi }(𝐱,t)=_0^t𝑑𝐱^{}K(𝐱,𝐱^{},t)\mathrm{\Psi }(𝐱^{}).$$
(2)
Here, $`\mathrm{\Psi }(𝐱^{})`$ is the initial wave function at $`t=0`$ and $`K(𝐱,𝐱^{},t)`$ denotes the propagator. We will not use the well-known semiclassical Van Vleck-Gutzwiller (VVG) propagator which is inconvenient for several reasons. Firstly, one has to deal with caustics, i.e. singularities of the propagator, and secondly, it is originally formulated as a boundary value problem. For numerical applications much better suited (and for analytical considerations not worse) is the so called Herman-Kluk (HK) propagator which is a uniformized propagator in initial value representation , formulated in phase space,

$`K^{HK}(𝐱,𝐱^{},t)`$ $`={\displaystyle \frac{1}{(2\pi \mathrm{})^n}}{\displaystyle 𝑑𝐩𝑑𝐪C_{\mathrm{𝐪𝐩}}(t)e^{iS_{\mathrm{𝐪𝐩}}(t)/\mathrm{}}}`$ (3)
$`g_\gamma (𝐱;𝐪(t),𝐩(t))g_\gamma ^{}(𝐱^{};𝐪,𝐩)`$ (4)
with

$$g_\gamma (𝐱;𝐪,𝐩)=\left(\frac{\gamma }{\pi }\right)^{n/4}\mathrm{exp}\left(\frac{\gamma }{2}\left(𝐱𝐪\right)^2+\frac{i}{\mathrm{}}𝐩\left(𝐱𝐪\right)\right)$$
(5)
and

$$C_{\mathrm{𝐪𝐩}}(t)=\left|\frac{1}{2}\left(𝐐_𝐪+𝐏_𝐩i\mathrm{}\gamma 𝐐_𝐩\frac{1}{i\mathrm{}\gamma }𝐏_𝐪\right)\right|^{\frac{1}{2}}.$$
(6)
Each phase space point ($`𝐪,𝐩`$) in the integrand of Eq. (3) is the starting point of a classical trajectory with action $`S_{\mathrm{𝐪𝐩}}(t)`$. The terms $`𝐗_𝐲`$ in the weight factor $`C_{\mathrm{𝐪𝐩}}(t)`$ are the four elements of the monodromy matrix, $`𝐗_𝐲=𝐱_t/𝐲`$. The square root in Eq. (6) has to be calculated in such a manner that $`C_{\mathrm{𝐪𝐩}}(t)`$ is a continuous function of $`t`$. The integrand in Eq. (3) is – depending on the system – highly oscillatory. Although we restrict ourselves to one spatial dimension (see Eq. (1)) the number of trajectories necessary for numerical convergence can reach $`10^7`$. We note in passing that an integration by stationary phase approximation over momentum and coordinate variables reduces the HK-propagator to the VVG-propagator .

In all calculations presented here we have used a Gaussian wave packet as initial wave function,

$$\mathrm{\Psi }_\beta (x^{})=\left(\frac{\beta }{\pi }\right)^{1/4}\mathrm{exp}\left(\frac{\beta }{2}\left(x^{}q_\beta \right)^2\right).$$
(7)
With this choice, the overlap

$$f_{\gamma \beta }(q,p)g_\gamma ^{}(x^{};q,p)\mathrm{\Psi }_\beta (x^{})𝑑x^{}$$
(8)
can be calculated analytically and Eq. (2) reads together with Eq. (3)

$`\mathrm{\Psi }_\beta ^{HK}(x,t)=`$ $`\left({\displaystyle \frac{4\gamma \beta }{\alpha ^2}}\right)^{\frac{1}{4}}{\displaystyle \frac{1}{2\pi \mathrm{}}}{\displaystyle 𝑑p𝑑qe^{iS_{qp}(t)/\mathrm{}}}`$ (9)
$`C_{qp}(t)g_\gamma (x;q(t),p(t))f_{\gamma \beta }(q,p)`$ (10)
with $`\alpha =\gamma +\beta `$. For all results presented here we have taken $`\gamma =\beta `$.

For comparison with our semiclassical calculations we determined the quantum mechanical wave function using standard Fast Fourier Transform split operator methods .

## III Above Threshold Ionization

We start from Eq. (1) with $`\delta =0`$ and use a rectangular pulse shape $`f(t)`$ which lasts for $`4.25`$ optical cycles. This setting is very similar to the one used in.

The energy spectrum of the electrons can be expressed by the Fourier transform of the autocorrelation function after the pulse, i.e. for times $`t>t_f`$,

$$\sigma (\omega )=Re\underset{t_f}{\overset{\mathrm{}}{}}e^{i\omega t}\mathrm{\Psi }(t)|\mathrm{\Psi }_f𝑑t,$$
(11)
where $`\mathrm{\Psi }_f=\mathrm{\Psi }(t_f)`$ is the wave function after the pulse and correspondingly

$$|\mathrm{\Psi }(t)=e^{iH_0(tt_f)/\mathrm{}}|\mathrm{\Psi }_f$$
(12)
is calculated by propagating $`\mathrm{\Psi }_f`$ for some time with the atomic Hamiltonian $`H_0`$ only after the laser has been switched off.

### A Quantum mechanical and semiclassical spectra for ATI

We will present results for two types of potentials to elucidate the dependence of the semiclassical approximation on the form of the potential.

#### 1 Softcore potential

First we apply the widely used softcore potential

$$V(x)=\frac{1}{\sqrt{x^2+a}}$$
(13)
with $`a=1`$ and with an ionization potential $`I_p=0.670`$ a.u.. We have checked that the correlation function differs little if calculated with the exact ground state or with the ground state wave function approximated by the Gaussian of Eq. (7) where $`\beta =0.431`$ a.u. and $`q_\beta =0`$. However, the semiclassical calculation is considerably simplified with a Gaussian as initial state as can be seen from Eqs. (7-10). Therefore we use this initial state and obtain the propagated semiclassical wavefunction in the closed form Eq. (10). In Fig. 1 the quantum and semiclassical results at a frequency $`\omega _0=0.148`$ a.u. and a field strength $`E_0=0.15`$ a.u. are compared. The Keldysh parameter has the value $`1.14`$. The quantum mechanical calculation (dotted line) shows a typical ATI spectrum. Intensity maxima with a separation in energy of $`\mathrm{}\omega _0`$ are clearly visible. The first maximum has the highest intensity while the second maximum is suppressed.

The semiclassical result (solid line) is ambiguous: On the one hand there are clear ATI maxima with a separation of $`\mathrm{}\omega _0`$. All peaks but the first one have roughly the correct magnitude. Again the second maximum is missing. On the other hand we see a constant shift (about $`0.02`$ a.u.) of the spectrum towards higher energies. Therefore, a quantitative semiclassical description is impossible, at least with the present parameters and the softcore potential. Next, we will clarify whether the shift in the spectrum is an inherent problem of a semiclassical ATI calculation or if it can be attributed to properties of the softcore potential.

#### 2 Gaussian potential

To this end we take a potential which has been used to model the “single bound state” situation mentioned in the introduction . It is of Gaussian form

$$V(x)=V_0\mathrm{exp}\left(\sigma x^2\right).$$
(14)
With our choice of parameters $`V_0=0.6`$ a.u. and $`\sigma =0.025`$ a.u., the potential contains six bound states and can be approximated, at least in the lower energy part, by a harmonic potential for which semiclassical calculations are exact. Hence, the semiclassical ATI spectrum with this potential should be more accurate if the discrepancies in Fig. 1 are due to the potential and not due to the laser interaction. The ground state wave function itself is again well approximated by the Gaussian Eq. (7) with $`\beta =0.154`$ a.u. and $`q_\beta =0`$. The laser has a frequency $`\omega _0=0.09`$ a.u., a field strength $`E_0=0.049`$ a.u., and a pulse duration of $`4.25`$ cycles. The Keldysh parameter has the value $`1.87`$.

We obtain a quantum mechanical ATI spectrum (dotted line in Fig. 2) with six distinct maxima. The semiclassical spectrum (solid line) is not shifted, the location of the maxima agrees with quantum mechanics. Hence, one can conclude that the softcore potential is responsible for the shift. The height of the third maximum is clearly underestimated and the details of the spectrum are exaggerated by the semiclassical calculation. Apart from these deviations the agreement is good enough to use this type of calculation as a basis for a semiclassical understanding of ATI.

### B Semiclassical interpretation of the ATI spectrum

#### 1 Classification and coherence of trajectories

With the chosen parameters most of the trajectories ionize during the pulse ($`92\%`$). We consider a trajectory as ionized if the energy of the atom

$$\epsilon (t)=p(t)^2/2+V(q(t))$$
(15)
becomes positive at some time $`t_n`$ and remains positive, i.e. $`\epsilon (t)>0`$ for $`t>t_n`$. Typically, the trajectories ionize around an extremum of the laser field. Tunnelling can not be very important, otherwise the agreement between quantum mechanics and semiclassics would be much worse. The Keldysh parameter of $`1.87`$ suggests that we are in between the tunnelling and the multiphoton regime. Interestingly, the semiclassical description is successful although we are way below energies of the classically allowed over the barrier regime.

An obvious criterion for the classification of the trajectories is the time interval of the laser cycle into which their individual ionization time $`t_n`$ falls, see Fig. 3. Typically ionization of trajectory happens around $`t_n=(2n1)T/4`$ when the force induced by the laser reaches a maximum. Hence, the ionized trajectories can be attached to time intervals $`I_n=[(n1)T/2,nT/2]`$. In Fig. 3 we have plotted four trajectories from the intervals $`I_1`$ to $`I_4`$ which end up with an energy $`E=0.36`$ a.u.. After ionization each trajectory shows a quiver motion around a mean momentum $`p_f`$ . One can distinguish two groups of intervals, namely those with trajectories ionized with positive momentum $`p_f`$ (the intervals $`I_{2k1}`$) and those with trajectories with negative $`p_f`$ (the intervals $`I_{2k}`$). These two groups contribute separately and incoherently to the energy spectrum as one might expect since the electrons are easily distinguishable. One can see this directly from the definition Eq. (11) of the electron energy spectrum. For relative high energies $`\mathrm{}\omega `$ the (short-range) potential may be neglected in the Hamiltonian $`H_0`$ and we get

$`\sigma (\omega )`$ $`=Re{\displaystyle \underset{t_f}{\overset{\mathrm{}}{}}}e^{i\omega t}\mathrm{\Psi }_f|e^{iH_0(tt_f)}|\mathrm{\Psi }_f𝑑t`$ (16)
$`Re{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}e^{i\omega t}\mathrm{\Psi }_f|e^{ip^2t/2\mathrm{}}|\mathrm{\Psi }_f𝑑t`$ (17)
$`={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\delta \left(\omega p^2/2\mathrm{}\right)\left|\mathrm{\Psi }_f(p)\right|^2𝑑p`$ (18)
$`=\left(\left|\mathrm{\Psi }_f(\sqrt{2\mathrm{}\omega })\right|^2+\left|\mathrm{\Psi }_f(\sqrt{2\mathrm{}\omega })\right|^2\right)(\mathrm{}/2\omega )^{1/2}`$ (19)
$`\sigma _{}(\omega )+\sigma _+(\omega ).`$ (20)
Hence, to this approximation, the ATI spectrum is indeed given by the incoherent sum of two terms belonging to different signs of the momenta of electrons ionized in different time intervals as described above.

Figure Fig. 4(a) shows that Eq. (20) is a good approximation.

Only for small $`\omega `$ the spectra do not agree, where the kinetic energy is comparable with the (neglected) potential energy.

Quantum mechanically, all contributions from trajectories which lead to the same momentum $`p_f`$ of the electron are indistinguishable and must be summed coherently. To double check that the interference from different intervals $`I_n`$ is responsible for the ATI peaks, we can artificially create a spectrum by an incoherent superposition $`\stackrel{~}{\sigma }_+=\sigma _2+\sigma _4+\sigma _6+\sigma _8`$ of contributions from trajectories ionized in the intervals $`I_{2j}`$. This artificially incoherent sum (Fig. 4(b)) shows similarity neither with $`\sigma _+(\omega )`$ nor with any kind of ATI spectrum.

#### 2 Classical signature of bound and continuum motion in the laser field

The great advantage of an ab initio semiclassical description lies in the possibility to make dynamical behavior transparent based on classical trajectories, particularly in the case of explicit time dependent problems where our intuition is not as well trained as in the case of conservative Hamiltonian systems. The classical quantities enter semiclassically mostly through the phase factor

$$\mathrm{exp}\left(i[S_{qp}(t)p(t)q(t)]/\mathrm{}\right)\mathrm{exp}[i\mathrm{\Phi }/\mathrm{}]$$
(21)
which each trajectory contributes to the wave function Eq. (10). Although the prefactor $`C_{qp}(t)`$ in Eq. (10) may be complex itself, the major contribution to the phase comes from the effective action $`\mathrm{\Phi }`$ in the exponent of Eq. (21). Figure 5 shows the energy $`\epsilon `$ of the atom and the accumulated phase $`\mathrm{\Phi }`$. One can recognize a clear distinction between a quasi-free oscillation in the laser field after the ionization and the quasi-bound motion in the potential. The latter is characterized by an almost constant averaged bound energy $`\epsilon (t)`$ (Fig. 5(a)) of the individual trajectory giving rise to an averaged linear increase of the phase (Fig. 5(b)). After ionization the phase decreases linearly with an oscillatory modulation superimposed by the laser field. The almost linear increase of $`\mathrm{\Phi }`$ without strong modulation of the laser field during the bound motion of the electron is remarkable, particularly looking at the laser induced modulations of the bound energy seen in Fig. 5(a). The averaged slope of the phase (positive for bound motion, negative for continuum motion) corresponds via $`d\mathrm{\Phi }/dt=E`$ to an averaged energy. The behavior can be understood by a closer inspection of the action

$`\mathrm{\Phi }(t)`$ $``$ $`S_{qp}(t)p(t)q(t)`$ (22)
$`=`$ $`{\displaystyle _0^t}(2TH\dot{p}(\tau )q(\tau )\dot{q}(\tau )p(\tau ))𝑑\tau qp.`$ (23)
Here, $`T=p^2(t)/2`$ refers to the kinetic energy and $`H`$ to the entire Hamiltonian of Eq. (1), the dot indicates a derivative with respect to time, and $`qq(t=0)`$. With the help of Hamilton’s equations and a little algebra $`\mathrm{\Phi }`$ from Eq. (22) can be simplified to

$$\mathrm{\Phi }(t)=_0^t\left(\epsilon (\tau )q(\tau )\frac{dV}{dq}\right)𝑑\tau $$
(24)
where $`\epsilon `$ is the atomic energy Eq. (15). With Eq. (24) we can quantitatively explain the slope of $`\mathrm{\Phi }`$ in Fig. 5(b). For the low energies considered the potential Eq. (14) can be approximated harmonically,

$$V(q)V_0+V_0\sigma q^2$$
(25)
Averaging $`\mathrm{\Phi }`$ over some time yields then $`\mathrm{\Phi }(t)V_0t`$, for any bound energy of a classical trajectory since for an oscillator averaged kinetic and potential energy are equal. Indeed, the numerical value for the positive slope in Fig. 5(b) is $`0.6`$ a.u. in agreement with the value for $`V_0`$.

For the ionized part of the trajectories we may assume that the potential vanishes. The corresponding solutions for electron momentum $`p(t)`$ follows directly from Hamilton’s equation $`\dot{p}=E_0\mathrm{sin}\omega _0t`$,

$$p(t)=\frac{E_0}{\omega _0}\mathrm{cos}(\omega _0t)+p,$$
(26)
where $`p`$ is the mean momentum. Without potential the phase from Eq. (24) reduces to $`\mathrm{\Phi }(t)=p^2(\tau )/2𝑑\tau `$ and we obtain with Eq. (26)

$`\mathrm{\Phi }_c(t)`$ (27)
$`={\displaystyle \frac{U_p}{2\omega _0}}\mathrm{sin}(2\omega _0t){\displaystyle \frac{E_0p}{\omega _0^2}}\mathrm{sin}\omega _0t(U_p+p^2/2)t`$ (28)
with the ponderomotive potential $`U_p=E_0^2/4\omega _0^2`$. We note in passing that Eq. (28) is identical to the time dependent phase in the Volkov state (see the appendix).

#### 3 Semiclassical model for ATI

The clear distinction between classical bound and continuum motion in the laser field as demonstrated by Fig. 5 and illuminated in the last section, allows one to derive easily the peak positions of the ATI spectrum. Moreover, this distinction also supports the so called strong field approximation (e.g. ) where electron dynamics in the laser field is modelled by one bound state and the continuum. While this is postulated in  as an approximation and justified a posteriori by the results the corresponding approximation is suggested in the present context of a semiclassical analysis by the full classical dynamics, i.e., the behavior of the trajectories, as shown in 5. There, we have seen that each classical bound motion leads to the characteristic linear increase of the phase. If the entire phase space corresponding to the initial (ground state) wave function is probed with many trajectories of different energy, the dominant contribution will appear at the bound state energy which implies

$$\mathrm{\Phi }_b(t)I_pt,$$
(29)
where $`I_p`$ is the ionization potential. The time for which a trajectory does not fall into one of the two classes, bound or continuum, is very short (Fig. 5). Hence, we can approximately compose the true phase $`\mathrm{\Phi }=\mathrm{\Phi }_b+\mathrm{\Phi }_c`$. However, we don’t know for an electron with mean momentum $`p`$ when it was ionized. Hence, we have to sum over all trajectories with different ionization times $`\tau `$ but equal final momentum $`p=p_f`$ which leads to the propagated wavefunction

$`\mathrm{\Psi }_f(t,p)`$ $`{\displaystyle _{t_0}^t}𝑑\tau \mathrm{exp}[i/\mathrm{}(\mathrm{\Phi }_b(\tau )+\mathrm{\Phi }_c(t)\mathrm{\Phi }_c(\tau ))]`$ (30)
$``$ $`{\displaystyle \underset{n,m}{}}J_n\left({\displaystyle \frac{E_0p}{\omega _0^2}}\right)J_m\left({\displaystyle \frac{U_p}{2\omega _0}}\right){\displaystyle _{t_0}^t}𝑑\tau e^{i\tau \mathrm{\Delta }_{mn}/\mathrm{}},`$ (31)
where the phase $`\mathrm{\Delta }`$ is given by

$$\mathrm{\Delta }_{mn}=I_p+U_p+p^2/2(n+2m)\mathrm{}\omega _0.$$
(32)
From Eq. (32) and Eq. (30) follows that ATI peaks appear at integer multiples $`n\mathrm{}\omega _0`$ of the laser frequency, when

$$\frac{p^2}{2}=n\mathrm{}\omega _0I_pU_p.$$
(33)
One can also see from Eq. (30) that the ATI maxima become sharper with each optical cycle that supplies ionizing trajectories. Of course, this effect is weakened by the spreading of the wavepacket hidden in the prefactor of each trajectory contribution (see Eq. (10)) not considered here.

Trajectories that are ionized during different laser cycles accumulate a specific mean phase difference. The phase difference depends on the number $`k`$ of laser cycles passed between the two ionization processes:

$$\mathrm{\Delta }\mathrm{\Phi }(p)=kT\left(I_p+\frac{p^2}{2}+U_p\right).$$
(34)
The trajectories interfere constructively if

$$\mathrm{\Delta }\mathrm{\Phi }(p)=\mathrm{\hspace{0.17em}2}\pi l\frac{1}{2}p^2=\frac{l}{k}\omega _0I_pU_p.$$
(35)
If an energy spectrum is calculated exclusively with trajectories from two intervals separated by $`k`$ cycles there should be additional maxima in the ATI spectrum with a distance $`\mathrm{}\omega _0/k`$.

As a test for this semiclassical interpretation of the ATI mechanism we have calculated three spectra with trajectories where the mean time delay between ionizing events is given by $`\mathrm{\Delta }t=T`$, $`\mathrm{\Delta }t=2T`$ and $`\mathrm{\Delta }t=3T`$. For the spectrum Fig. 6 (a) we have used exclusively trajectories from the intervals $`I_2`$ and $`I_4`$ ($`\mathrm{\Delta }t=T`$). One can see broad maxima separated by $`\mathrm{}\omega _0`$ in energy. Trajectories from the intervals $`I_2`$ and $`I_6`$ (see Fig. 6 (b)) form a spectrum where the maxima are separated by $`\mathrm{}\omega _0/2`$ – as predicted for $`\mathrm{\Delta }t=2T`$. In analogy the separation for the ATI maxima in a spectrum with trajectories from the intervals $`I_2`$ and $`I_8`$ is given by $`\mathrm{}\omega _0/3`$ (Abb. 6 (c)). The interference of trajectories ionized in many subsequent cycles suppresses the non-integer maxima according to Eq. (32). If the field strength is high enough the atom is completely ionized during the first cycle. The opportunity for interference gets lost and we end up with an unstructured energy spectrum.

In an extreme semiclassical approximation we would have evaluated the integral in Eq. (30) by stationary phase. The condition

$$d/d\tau [\mathrm{\Phi }_b(\tau )\mathrm{\Phi }_c(\tau )]I_p+p^2(\tau )/2=0$$
(36)
leads to complex ionization times $`t_n`$ whose real part is periodic and allows for two ionizing events per laser cycle, close to the extrema of the laser amplitude. The derivation is simple but technical, therefore we don’t carry it out explicitely here. However, it explains the observation that ionization occurs close to the extrema of the laser field and it also makes contact with the tunnelling process often referred to in the literature since the complex time can be interpreted as tunnelling at a complex ”transition” energy.

Clearly, our semiclassical analysis as described here supports the picture which has been sketched in  interpreting a quantum calculation. The authors assume that wave packets are emitted every time the laser reaches an extremum. The interference of the different wave packets gives rise to the ATI peaks.

In the following we will discuss the process of higher harmonic generation (HHG) which is closely related to ATI. In fact, the separation into a bound and continuum part of the electron description is constitutive for HHG as well, the prominent features, such as cutoff and peak locations, can be derived from the same phase properties Eq. (30) as for ATI. However, there is a characteristic difference, how these phases enter.

## IV High Harmonic Generation

First, we briefly recapitulate the findings of , where we have calculated the harmonic spectrum with the softcore potential Eq. (13). With our choice of $`a=2`$ the ionization potential is given by $`I_p=0.5`$ a.u.. The laser field has a strength $`E_0=0.1`$ a.u., a frequency $`\omega _0=0.0378`$ a.u. and a phase $`\delta =\pi /2`$. The initial wave packet with a width of $`\beta =0.05`$ a.u. is located at $`q_\beta =E_0/\omega _0^2=70`$ a.u.. Note, that the cutoff energy $`E_C`$ in such a symmetric laser scattering experiment is given by

$$E_C=I_p+2U_p.$$
(37)
From the dipole acceleration (see Fig. 7)

$$d(t)=\mathrm{\Psi }(t)\left|\frac{dV(x)}{dx}\right|\mathrm{\Psi }(t),$$
(38)
follows by Fourier transform

$$\sigma (\omega )=d(t)\mathrm{exp}(i\omega t)𝑑t$$
(39)
the harmonic power spectrum (see Fig. 8).

Clearly, our semiclassical approach represents a good approximation. The dipole acceleration shows the characteristic feature that fast oscillations (which are responsible for the high harmonics in Fourier space) show only up after some time, here after $`t=T`$. This is the first time where trajectories are trapped. Trapping can only occur if (i) $`t_n=nT/2`$, (ii) the trajectories reach a turning point (i.e. $`p(t_n)=0`$), and (iii) if at this time the electron is close to the nucleus ($`q(t_n)0`$). The trapped trajectories constitute a partially bound state which can interfere with the main part of the wavepacket (trajectories) still bouncing back and forward over the nucleus driven by the laser. The group of briefly bound (i.e. trapped or stranded) trajectories can be clearly identified, either by their small excursion in space (Fig. 9a) or by the positive slope of their action (Fig. 9b) as it was the case for ATI (compare with Fig. 5). By artificially discarding the initial conditions in the semiclassical propagator which lead to trapped trajectories one can convincingly demonstrate that the plateau in HHG generation is a simple interference effect . Here, we are interested firstly in linking ATI to HHG by using the same separation in bound and continuum parts of the dynamics already worked out for ATI. Secondly, we want to go one step further and construct a wavefunction based on this principle.

Semiclassically, we have to look first at the phases of the observable. Therefore, we define a linear combination for the wavefunction from the respective phase factors for bound and continuum motion. Considering only terms in the exponent the harmonic spectrum Eq. (39) reads simply

$`\sigma (\omega ){\displaystyle 𝑑t\mathrm{exp}(i\omega t)}`$ (40)
$`\left|\mathrm{exp}\left(i\mathrm{\Phi }_c(t)/\mathrm{}\right)+c\mathrm{exp}\left(i\mathrm{\Phi }_b(t)/\mathrm{}\right)\right|^2,`$ (41)
where $`c0`$ is a (so far) arbitrary constant. In principle, $`c=c(t)`$, however its change in time is much slower than that of the optical oscillations of the phases $`\mathrm{\Phi }(t)`$, hence we may approximate $`c`$ by a constant. The bound and continuum phases, $`\mathrm{\Phi }_b`$ and $`\mathrm{\Phi }_c`$, are defined in Eq. (29) and Eq. (28), respectively. For $`\mathrm{\Phi }_c`$ we have $`p=0`$, since this is the dominant contribution from the center of the wavepacket which was initially at rest. The result is shown in Fig. 10. Indeed, the plateau with the harmonics is generated, however, the initial exponential decrease is missing since we have neglected all prefactors of the semiclassical wavefunction which describe the dispersion of the wavepacket.

Consequently, one can evaluate Eq. (40) in stationary phase approximation. The integrand of Eq. (40) becomes stationary if

$$\frac{d}{dt}\left[\mathrm{}\omega t\pm \left(\mathrm{\Phi }_b(t)\mathrm{\Phi }_c(t)\right)\right]=\mathrm{\hspace{0.17em}0}$$
(42)
which happens at

$$\mathrm{}\omega =2U_p\mathrm{sin}^2(\omega t)+I_p.$$
(43)
From Eq. (43) we conclude the cut-off law

$$\omega _{\text{max}}=2U_p+I_p,$$
(44)
as expected for laser assisted electron ion scattering . Using the same expansion into Bessel functions as in Eq. (30) we obtain for the spectrum Eq. (40):

$`{\displaystyle 𝑑t\mathrm{exp}\left(\frac{i}{\mathrm{}}\left[\left(\mathrm{}\omega U_pI_p\right)t+\frac{U_p}{2\omega _0}\mathrm{sin}\left(2\omega _0t\right)\right]\right)}`$ (45)
$`={\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle 𝑑te^{it\left(\mathrm{}\omega U_pI_p+2k\mathrm{}\omega _0\right)/\mathrm{}}\text{J}_k\left(\frac{U_p}{2\mathrm{}\omega _0}\right)}.`$ (46)
Therefore, we see maxima in the harmonic spectrum for

$$\mathrm{}\omega _k=U_p+I_p2k\omega _0.$$
(47)
We can go one step further and construct a full time-dependent wavefunction from this semiclassical approximation, namely

$$\mathrm{\Psi }(x,t)=\mathrm{\Psi }_\beta ^{sc}(x,t)+c\mathrm{\Psi }_0(x)\mathrm{exp}(itI_p/\mathrm{}).$$
(48)
Here, $`\mathrm{\Psi }_0(x)\mathrm{exp}(iI_pt/\mathrm{})`$ is the time dependent ground state wave function (without the laser field) and $`\mathrm{\Psi }_\beta ^{sc}(x,t)`$ is a (semiclassical) wavepacket in the laser field but without potential. Calculating the dipole acceleration and the resulting harmonic spectrum with this wavefunction leads to a remarkably good approximation of the true quantum spectrum (compare Fig. 8 with Fig. 11). The dispersion of the wavepacket leads to the lower plateau compared to Fig. 10.

## V Conclusions

### A Semiclassical comparison between ATI and HHG

Clearly, the main structure such as the plateau, cutoff (HHG) and the occurrence of peaks and their separation in energy (ATI and HHG) is a property of the difference of the classical time-dependent actions $`\mathrm{\Phi }_b(t)\mathrm{\Phi }_c(t)`$ alone.

However, the HHG power spectrum Eq. (39) is an integral over all the time for which the electron wavepacket is exposed to the laser field. In contrast, the ATI spectrum is obtained in the long-time limit $`t\mathrm{}`$ after the laser has been switched off. This difference may explain why the HHG results tend to be better than the ATI results semiclassically: Any semiclassical approximation (which is not exact) become worse for large times.

A second point refers to the fact that the characteristic phase difference $`\mathrm{\Phi }_b(t)\mathrm{\Phi }_c(t)`$ appears already in the wavefunction Eq. (30) for ATI, while for HHG it occurs only in the expectation value Eq. (38). However, this difference is artificial, since the expectation value, or better its Fourier transform the power spectrum, is not the observable of higher harmonic radiation. The correct expression is the dipole-dipole correlation function $`R`$ which can be approximated as $`R|\sigma (\omega )|^2`$ under single atom conditions or in the case of an ensemble of independent atoms which radiate . Hence, in both cases, ATI and HHG, the peak structure appears already on the level of the quantum amplitude (or wavefunction) and is amplified in the true observable.

### B Summary

We have given a time-dependent fully semiclassical description of multiphoton processes. The prominent ATI and HHG features emerge naturally from properties of the classical trajectories whose contributions to the semiclassical wavefunction interfere semiclassically. Any effect of this semiclassical interference can be double-checked by disregarding the phases. This leads (with the same trajectories) to a classical observable. As we have seen, to a good approximation the classical action for an individual trajectory can be composed of one part $`\mathrm{\Phi }_b`$ for the time the electron is bound (disregarding the laser field) and of another part $`\mathrm{\Phi }_c`$ for the time the electron is in the continuum (disregarding the atomic potential). The relevant phase difference $`\mathrm{\Phi }_b\mathrm{\Phi }_c`$ leads in both cases, ATI and HHG, to the prominent harmonic structures in terms of the laser energy $`\mathrm{}\omega _0`$. Finally, we have been able to construct a simple wavefunction for higher harmonics generated in laser assisted scattering. Its key element is an explicitely time-dependent wavepacket of the electron under the influence of the laser field. Starting from an initial Gaussian distribution localized in space the wavepacket disperses in time providing the correct decrease of the intensity of the lower harmonics and in turn the correct height of the plateau.

Financial support from the DFG under the Gerhard Hess-Programm and the SFB 276 is gratefully acknowledged.

##

We want to calculate the semiclassical wave function of a free particle in a laser field according to Eq. (10). A particle in a laser field $`V_L(x,t)=E_0\mathrm{sin}(\omega t)`$ moves with

$`p(t)=p+{\displaystyle \frac{E_0}{\omega }}\mathrm{cos}(\omega t)p+\stackrel{~}{p}(t)`$ (49)
$`q(t)=q+pt+{\displaystyle \frac{E_0}{\omega ^2}}\mathrm{sin}(\omega t)q+pt+\stackrel{~}{q}(t)`$ (50)
The weight factor $`C_{qp}(t)`$ is given by

$$C_{qp}(t)=\left(1\frac{i\mathrm{}\gamma }{2}t\right)^{\frac{1}{2}}.$$
(51)
For the phase factor $`S_{qp}(t)p(t)q(t)`$ we get:

$`S_{qp}(t)p(t)q(t)=`$ $`{\displaystyle \frac{U_p}{2\omega }}\mathrm{sin}(2\omega t)U_pt`$ (52)
$`{\displaystyle \frac{p^2}{2}}t\stackrel{~}{q}(t)pqp`$ (53)
Evaluating Eq. (10) with the stationary phase approximation, which is exact for quadratic potentials, leads to the condition that

$`f(q,p)=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\left(xp(t){\displaystyle \frac{p^2}{2}}t\stackrel{~}{q}(t)p{\displaystyle \frac{\gamma }{\alpha }}qp{\displaystyle \frac{\beta }{\alpha }}q_\beta p\right)`$ (54)
$`{\displaystyle \frac{\gamma }{2}}\left(xq(t)\right)^2{\displaystyle \frac{\gamma \beta }{2\alpha }}\left(qq_\beta \right)^2{\displaystyle \frac{1}{2\mathrm{}^2\alpha }}p^2`$ (55)
must have an extremum. With

$`{\displaystyle \frac{f}{q}}=\mathrm{\hspace{0.17em}0}=\gamma \left[xq(t)\right]{\displaystyle \frac{\gamma \beta }{\alpha }}(qq_\beta ){\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \frac{\gamma }{\alpha }}p`$ (56)
$`{\displaystyle \frac{f}{p}}=\mathrm{\hspace{0.17em}0}=\gamma \left[xq(t)\right]t{\displaystyle \frac{1}{\mathrm{}^2\alpha }}p`$ (57)
$`+{\displaystyle \frac{i}{\mathrm{}}}\left(xpt\stackrel{~}{q}(t){\displaystyle \frac{\gamma }{\alpha }}q{\displaystyle \frac{\beta }{\alpha }}q_\beta \right)`$ (58)
we find

$`q_s={\displaystyle \frac{x\stackrel{~}{q}(t)+i\mathrm{}\beta tq_\beta }{1+i\mathrm{}\beta t}}`$ (59)
$`p_s={\displaystyle \frac{i\mathrm{}\beta }{1+i\mathrm{}\beta t}}\left(x\stackrel{~}{q}(t)q_\beta \right).`$ (60)
After some algebra we arrive at the stationary exponent

$`f(q_s,p_s)`$ $`={\displaystyle \frac{i}{\mathrm{}}}x\stackrel{~}{p}(t){\displaystyle \frac{\beta }{2\left(1+i\mathrm{}\beta t\right)}}\left(x\stackrel{~}{q}(t)q_\beta \right)^2`$ (61)
$`={\displaystyle \frac{i}{\mathrm{}}}x\stackrel{~}{p}(t){\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \frac{\mathrm{}^2\beta ^2t}{2\sigma (t)}}\left(x\stackrel{~}{q}(t)q_\beta \right)^2`$ (62)
$`{\displaystyle \frac{\beta }{2\sigma (t)}}\left(x\stackrel{~}{q}(t)q_\beta \right)^2,`$ (63)
where $`\sigma (t)`$ is given by

$$\sigma (t)=\mathrm{\hspace{0.17em}1}+\beta ^2\mathrm{}^2t^2.$$
(64)
The determinant of the second derivatives of $`f`$ still has to be calculated. With

$`{\displaystyle \frac{^2f}{q^2}}={\displaystyle \frac{\gamma ^4+2\gamma \beta }{\alpha }}{\displaystyle \frac{^2f}{p^2}}={\displaystyle \frac{i}{\mathrm{}}}t\gamma t^2{\displaystyle \frac{1}{\mathrm{}^2\alpha }}`$ (65)
$`{\displaystyle \frac{^2f}{qp}}={\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \frac{\gamma }{\alpha }}\gamma t`$ (66)
we get

$$det\left(\begin{array}{cc}\frac{^2f}{q^2}& \frac{^2f}{qp}\\ \frac{^2f}{pq}& \frac{^2f}{p^2}\end{array}\right)=\frac{2\gamma }{\mathrm{}^2\alpha }\left(\left[1i\gamma \mathrm{}t/2\right]\left[1+i\beta \mathrm{}t\right]\right).$$
(67)
The factor $`\gamma `$ cancels as it should be and we are left with

$`\mathrm{\Psi }_\beta ^{sc}(x,t)=`$ $`\left({\displaystyle \frac{\beta }{\pi }}\right)^{1/4}\sqrt{{\displaystyle \frac{1}{1+i\mathrm{}\beta t}}}`$ (68)
$`\mathrm{exp}\left({\displaystyle \frac{i}{\mathrm{}}}\left[\stackrel{~}{p}(t)x{\displaystyle \frac{U_p}{2\omega }}\mathrm{sin}(2\omega t)U_pt\right]\right)`$ (69)
$`\mathrm{exp}\left({\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \frac{\mathrm{}^2\beta ^2}{2\sigma (t)}}\left(x\stackrel{~}{q}(t)q_\beta \right)^2t\right)`$ (70)
$`\mathrm{exp}\left({\displaystyle \frac{\beta }{2\sigma (t)}}\left(x\stackrel{~}{q}(t)q_\beta \right)^2\right).`$ (71)
This semiclassical time dependent wavepacket is quantum mechanically exact and corresponds to a superposition of Volkov solutions according to a Gaussian distribution at time $`t=0`$ . The fact that the semiclassical wavefunction is exact is a direct consequence of the Ehrenfest theorem which implies that interactions $`Vx^n`$, $`n=0,1,2`$ have quantum mechanically exact semiclassical solutions.


text:
# Stochastic models of exotic transport

## 1 Prerequisites

A simple Brownian motion, in its canonical (model) manifestations, does not seem to hide any surprises. It is our main goal in the present paper to reconsider those ingredients of the standard formalism which, while relaxed or slightly modified, would lead to conceptually new and possibly exotic (in the sense of being non-typical) features.

One of possible ”defects” of the standard theory is that the Brownian motion is incapable to originate spatiotemporal patterns (structures) and rather washes them out, even if initially in existence. In this particular context, ”active” Brownian particles were introduced, -. Most generally, they are supposed to remain in a feedback relationship with the environment (heat reservoir, thermostat, random medium, deterministic driving system or whatever else) through which they propagate. If regarded as a thermal bath, the environment is basically out of equilibrium. Obviously, a concrete meaning of being active relies on a specific choice of the model for the thermostat (deterministic or random, linear or nonlinear etc.) and the detailed particle-medium coupling mechanism.

Another ”defective” feature of the standard theory is rooted in its possible molecular (chaos) foundations where the major simplification is needed while passing from the Boltzmann collision theory to the novel kinetic framework set by the Kramers equation, cf. . In the process, microscopic (molecular) energy-momentum conservation laws completely evaporate from the formalism. Then, there is no obvious way how to reconcile random exchanges of energy and momentum between the Brownian particle and the thermostat with the manifest isothermality requirement. An issue of the heat production/removal and resulting thermal inhomogeneities is normally ignored (cf. however - for a discussion of various thermostat models, and the notion of thermostatting in non-equilibrium dynamics).

That was the starting point of our discussion of the origins of isothermal flows and the ultimate usage of the third Newton law in the mean to justify the concept of the Brownian motion with a recoil (via feedback relationship with the bath), .

In conformity with the metaphor, : ”everything depends on everything else”, we relax the habitual sharp distinction between an entirely passive particle (nonetheless performing a continual erratic motion)and its exclusively active (perpetually at rest and in equlibrium) driving agency - the thermostat.

In the present paper we take the following point of view on major signatures of the particle activity while moving at random due to external (environmental) reasons: any action upon the particle exerted by the environment induces (on suitable space and time scales) a compensating reaction in the carrier random medium. In loose terms of Refs. - we attribute to any (even single) Brownian particle an ability to generate perturbations of the medium (named a selfconsistent field in Ref.  and perturbations of noise in Ref. ) which in turn influence its further movement. That is a non-liner feedback relationship mentioned before: while inducing a random dynamics of a particle (action), the medium suffers alterations (reaction), which modify next stages of the particle motion (feedback). Such features are obviously alien to the standard formalism.

There is no experimentally reliable way to watch Brownian particles individually on time scales below say $`\frac{1}{100}s`$, nor get any insight into the associated fine detailed particle-bath interaction phenomena. The only realistic procedure to quantify observable features of the Brownian motion is to exploit its familiar Janus face (cf. ). Namely, the same mathematical formalism essentially applies to a single particle and to a statistical ensemble of identical non-interacting Brownian particles. Then a hydrodynamical analogy can be exploited on both levels of description, , once we invoke a standard kinetic reasoning (used in passing from the Boltzmann equation to gas or liquid dynamics equations).

Effectively, we need to evaluate (conditional) averages over a statistical ensemble of particles, which (in loose terminology that follows Ref. , see also ) constitutes one component (noninteracting Brownian ”gas”) of a coupled two-component continuous system. Another component is the thermostat.

What are the ultimate kinetic features of the Brownian motion (diffusion process) critically relies on the thermostat model i.e. specific mechanisms of the energy/heat exchange due to the coupling between the thermostat and the Brownian ”gas”. This particular issue we shall investigate in some detail in below, by resorting to somewhat nontypical (exotic in the Brownian motion context) methods.

## 2 Local conservation laws for the Brownian motion and Smoluchowski-type diffusion

It is useful to exploit a standard phase-space argument that is valid, under isothermal conditions, for a Markovian diffusion process taking place in (or relative to) a driving flow $`\stackrel{}{w}(\stackrel{}{x},t)`$ with as yet unspecified dynamics nor concrete physical origin. (In particular, such a flow can receive an interpretation of a selfconsistent field generated in the environment by Brownian particles themselves, cf. Ref. .)

We account for an explicit external force (here, acceleration $`\stackrel{}{K}=\stackrel{}{F}/m`$) exerted upon diffusing particles, while not directly affecting the driving flow itself. Then, infinitesimal increments of phase-space random variables read:

$$d\stackrel{}{X}(t)=\stackrel{}{V}(t)dt$$
(1)
$$d\stackrel{}{V}(t)=\beta [\stackrel{}{w}(\stackrel{}{x},t)\stackrel{}{V}(t)]dt+\stackrel{}{K}(\stackrel{}{x})dt+\beta \sqrt{2D}d\stackrel{}{W}(t).$$
(2)
Following the leading idea of the Smoluchowski approximation, we assume that $`\beta `$ is large, and consider the process on time scales significantly exceeding $`\beta ^1`$. Then, an appropriate choice of the velocity field $`\stackrel{}{w}(\stackrel{}{x},t)`$ may in principle guarantee the convergence of the spatial part $`\stackrel{}{X}(t)`$ of the process to the Itô diffusion process with infinitesimal increments :

$$d\stackrel{}{X}(t)=\stackrel{}{b}(\stackrel{}{x},t)dt+\sqrt{2D}d\stackrel{}{W}(t).$$
(3)
In this case, the forward drift of the process reads $`\stackrel{}{b}(\stackrel{}{x},t)=\stackrel{}{w}(\stackrel{}{x},t)+\frac{1}{\beta }\stackrel{}{K}(\stackrel{}{x})`$. Notice that the $`\beta ^1\stackrel{}{K}`$ contribution can be safely ignored if we are interested in the dominant driving motion.

Throughout the paper we are interested in Markovian diffusion processes, which propagate respectively the phase-space or configuration space probability densities (weak solutions of stochastic differential equations are thus involved). In the configuration space variant, we deal with a Markovian stochastic process whose probability density $`\rho (\stackrel{}{x},t)`$ evolves according to the standard Fokker-Planck equation

$$_t\rho =D\mathrm{}\rho \stackrel{}{}(\stackrel{}{b}\rho )$$
(4)
and the forward drift is not necessarily of the form dictated by Eqs. (1) and (2). We admit here more general forward drift functions, cf. Ref. , which do not allow for a simple additive decomposition and thus are capable to give account of nonlinearities in the particle-bath coupling.

One can easily transform the Fokker-Planck equation to the familiar form of the continuity equation (hydrodynamic mass conservation law) $`_t\rho =\stackrel{}{}(\stackrel{}{v}\rho )`$ by defining $`\stackrel{}{v}=\stackrel{}{b}D\frac{\stackrel{}{}\rho }{\rho }`$. The current velocity $`\stackrel{}{v}`$ obeys a local (momentum per unit of mass) conservation law which directly originates from the rules of the Itô calculus for Markovian diffusion processes, and from the first moment equation in the diffusion approximation of the Kramers theory, :

$$_t\stackrel{}{v}+(\stackrel{}{v}\stackrel{}{})\stackrel{}{v}=\stackrel{}{}(\mathrm{\Omega }Q).$$
(5)
While looking similar to the standard Euler equation appropriate for the lowest order hydrodynamical description of gases and liquids, this equation conveys an entirely different physical message.

First of all, for a class of forward drifts that are gradient fields the most general admissible form of an auxiliary potential $`\mathrm{\Omega }(\stackrel{}{x},t)`$ reads:

$$\mathrm{\Omega }(\stackrel{}{x},t)=2D[_t\varphi +\frac{1}{2}(\frac{\stackrel{}{b}^2}{2D}+\stackrel{}{}\stackrel{}{b})].$$
(6)
Here $`\stackrel{}{b}(\stackrel{}{x},t)=2D\stackrel{}{}\varphi (\stackrel{}{x},t)`$. In reverse, by choosing a bounded from below continuous function to represent conservative force fields (after taking the gradient) i. e. otherwise arbitrary $`\mathrm{\Omega }`$, we can always disentangle the above (Riccatti-type) identity with respect to the drift field.

Moreover, instead of the standard pressure term (consider a state equation $`P\rho ^\alpha ,\alpha >0`$), there appears a contribution from more complicated (derivatives !) $`\rho `$-dependent potential $`Q(\stackrel{}{x},t)`$. It is is given in terms of the so-called osmotic velocity field $`\stackrel{}{u}(\stackrel{}{x},t)=D\stackrel{}{}ln\rho (\stackrel{}{x},t)`$:

$$Q(\stackrel{}{x},t)=\frac{1}{2}\stackrel{}{u}^2+D\stackrel{}{}\stackrel{}{u}.$$
(7)
An equivalent form of the enthalpy-related potential $`Q`$ is $`Q=2D^2\frac{\mathrm{}\rho ^{1/2}}{\rho ^{1/2}}`$.

A general expression for the local diffusion current is $`\stackrel{}{j}=\rho \stackrel{}{v}=\rho (\stackrel{}{b}D\frac{\stackrel{}{}\rho }{\rho })`$. This local flow in principle may be experimentally observed for a cloud of suspended particles in a liquid. The current $`\stackrel{}{j}`$ is nonzero in non-equilibrium situations and a non-negligible matter transport occurs as a consequence of the Brownian motion, on the ensemble average. We thus cannot avoid local heating/cooling phenomena that need to push the environment out of equilibrium. That leads to obvious temperature inhomogeneities, which are normally disregarded, cf. Ref. .

If the forward drift is interpreted as a gradient of a suitable function ($`b=2D\varphi =\mathrm{\Phi }`$) and we take $`\mathrm{\Phi }(\stackrel{}{x},0)`$ as the initial data for the $`t0`$ evolution), then we have:

$$\mathrm{\Omega }=_t\mathrm{\Phi }+\frac{1}{2}|\stackrel{}{}\mathrm{\Phi }|^2+D\mathrm{}\mathrm{\Phi }.$$
(8)
If we decide that the above Hamilton-Jacobi-type equation is to be solved with respect to the field $`\mathrm{\Phi }(\stackrel{}{x},t)`$, its solution (and general solvability issue) relies on the choice of a bounded from below, continuous function $`\mathrm{\Omega }(\stackrel{}{x},t)`$, without any a priori knowledge of forward drifts. Viewed that way, Eq. (7) sets limitations on admissible forms of the space-time dependence of any conceivable self-consistent field to be generated in the bath by Brownian particles (cf. Ref. ). There is no freedom at all for postulating various partial differential equations, if their solutions are to be interpreted as forward drifts of Markovian diffusion processes, .

There is also interesting to observe that a gradient field ansatz for the diffusion current velocity $`\stackrel{}{v}=\stackrel{}{}S`$ allows to transform the momentum conservation law of a Markovian diffusion process to the universal Hamilton-Jacobi form:

$$\mathrm{\Omega }=_tS+\frac{1}{2}|\stackrel{}{}S|^2+Q$$
(9)
where $`Q(\stackrel{}{x},t)`$ was defined before as the ”pressure/enthalpy”-type function. (That form looks deceivingly similar to the standard hydrodynamic conservation law, valid for liquids and gases at thermal equilibrium, where instead of $`Q`$ an enthaply function normally appears.) By performing the gradient operation we recover the previous hydrodynamical form of the law.

In the above, the contribution due to $`Q`$ is a direct consequence of an initial probability measure choice for the diffusion process, while $`\mathrm{\Omega }`$ alone does account for an appropriate forward drift of the process, playing at the same time the role of the volume force potential ($`Q=\frac{P}{\rho }`$ contributes to energy-momentum transfer effects through the boundaries of any volume).

## 3 Moment equations for the free Brownian motion: Getting out of conventions

The derivation of a hierarchy of local conservation laws (moment equations) for the Kramers equation can be patterned after the standard procedure for the Boltzmann equation, . Those laws do not form a closed system and additional specifications (like the familiar thermodynamic equation of state) are needed to that end. In case of the isothermal Brownian motion, when considered in the large friction regime (e.g. Smoluchowski diffusion approximation), it suffices to supplement the Fokker-Planck equation by one more conservation law only to arrive at a closed system.

To give a deeper insight into what really happens on the way from the phase-space theory of the Brownian motion to its approximate configuration-space (Smoluchowski) version, let us consider the familiar Ornstein-Uhlenbeck process (in velocity/momentum) in its extended phase-space form. For clarity of discussion, we discus random dynamics for one degree of freedom only.

In the absence of external forces, the kinetic (Kramers-Fokker-Planck equation) reads:

$$_tW+u_xW=\beta _u(Wu)+q\mathrm{}_uW$$
(10)
where $`q=D\beta ^2`$. Here $`\beta `$ is the friction coefficient, $`D`$ will be identified later with the spatial diffusion constant, and provisionally we set $`D=k_BT/m\beta `$ in conformity with the Einstein fluctuation-dissipation identity.

The joint probability distribution (in fact, density) $`W(x,u,t)`$ for a freely moving Brownian particle which at $`t=0`$ initiates its motion at $`x_0`$ with an arbitrary inital velocity $`u_0`$ can be given in the form of the maximally symmetric displacement probability law:

$$W(x,u,t)=W(R,S)=[4\pi ^2(FGH^2)]^{1/2}exp\{\frac{GR^2HRS+FS^2}{2(FGH^2)}\}$$
(11)
where $`R=xu_0(1e^{\beta t})\beta ^1`$, $`S=uu_0e^{\beta t}`$ while $`F=\frac{D}{\beta }(2\beta t3+4e^{\beta t}e^{2\beta t})`$$`G=D\beta (1e^{2\beta t})`$ and $`H=D(1e^{\beta t})^2`$.

Marginal probablity densities, in the Smoluchowski regime (we take for granted that time scales $`\beta ^1`$ and space scales $`(D\beta ^1)^{1/2}`$ are irrelevant), take familiar forms of the Maxwell-Boltzmann $`w(u,t)=(\frac{m}{2\pi kT})^{1/2}exp(\frac{mu^2}{2k_BT})`$ and the diffusion kernel $`w(x,t)=(4\pi Dt)^{1/2}exp(\frac{x^2}{4Dt})`$ respectively.

A direct evaluation of the first and second local moment of the phase-space probability density:

$$<u>=𝑑uuW(x,u,t)=w(R)[(H/F)R+u_0e^{\beta t}]$$
(12)
$$<u^2>=𝑑uu^2W(x,u,t)=(\frac{FGH^2}{F}+\frac{H^2}{F^2}R^2)(2\pi F)^{1/2}exp(\frac{R^2}{2F})$$
(13)
after passing to the diffusion (Smoluchowski) regime, , allows to recover the local (configuration space conditioned) moment $`<u>_x=\frac{1}{w}<u>`$ which reads

$$<u>_x=\frac{x}{2t}=D\frac{w(x,t)}{w(x,t)}$$
(14)
while for the second local moment $`<u^2>_x=\frac{1}{w}<u^2>`$ we arrive at

$$<u^2>_x=(D\beta D/2t)+<u>_x^2.$$
(15)
By inspection one verifies that the transport (Kramers) equation for $`W(x,u,t)`$ implies local conservation laws:

$$_tw+(<u>_xw)=0$$
(16)
and

$$_t(<u>_xw)+_x(<u^2>_xw)=\beta <u>_xw.$$
(17)
By introducing (we strictly follow the moment equations strategy of the traditional kinetic theory of gases and liquids) the notion of pressure function $`P_{kin}`$ (we choose another notation to make a difference with the previous notion $`P`$ of the pressure function, cf. $`Q=\frac{P}{\rho }`$):

$$P_{kin}(x,t)=(<u^2>_x<u>_x^2)w(x,t)$$
(18)
we can analyze the local momentum conservation law

$$(_t+<u>_x)<u>_x=\beta <u>_x\frac{P_{kin}}{w}$$
(19)
In the Smoluchowski regime the friction term is cancelled away by a counterterm coming from $`\frac{1}{w}P_{kin}`$ so that

$$(_t+<u>_x)<u>_x=\frac{D}{2t}\frac{w}{w}=\frac{P}{w}=Q$$
(20)
where $`P=D^2w\mathrm{}lnw`$.

For comparison with the notations (local conservation laws) of the previous section one needs to replace $`w(x,t)`$ by $`\rho (x,t)`$ and $`<u>_x`$ by $`v(x,t)`$ in all formulas that pertain to the Smoluchowski regime.

Further exploiting the kinetic lore, we can tell few words about the temperature of Brownian particles as opposed to the (posibly equilibrium) temperature of the thermal bath. Namely, in view of (we stay in the Smoluchowski regime) $`P_{kin}(D\beta \frac{D}{2t})w`$ where $`D=\frac{k_BT}{mm\beta }`$, we can formally set:

$$\mathrm{\Theta }=\frac{P_{kin}}{w}(k_BT\frac{D}{2t})<k_BT.$$
(21)
That quantifies the degree of thermal agitation (temperature) of Brownian particles to be less than the thermostat temperature. Heat is continually pumped from the thermostat to the Brownian ”gas”, until asymptotically both temperatures equalize. This may be called a ”thermalization” of Brownian particles. In the process of thermalization the Brownian ”gas” temperature monotonically grows up until the mean kinetic energy of particles and that of mean flows asymptotically approach the familiar kinetic relationship:

$$\frac{w}{2}(<u^2>_x<u>_x^2)𝑑x=k_BT$$
(22)
In view of this medium $``$ particles heat transfer issue, one must be really careful while associating habitual thermal equilibrium conditions with essentially non-equilibrium phenomena, cf. Ref.  for more extended discussion.

## 4 Hydrodynamical reasoning: In search for exotic

Once local conservation laws were introduced, it seems instructive to comment on the essentially hydrodynamical features (compressible fluid/gas case) of the problem. Specifically, the ”pressure” term $`Q`$ is here quite annoying from the traditional kinetic theory perspective. That is quite apart from the fact that our local conservation laws have a conspicuous Euler form appropriate for the standard hydrodynamics of gases and liquids. Let us stress that in case of normal liquids the pressure is exerted upon any control volume (droplet) by the surrounding fluid. We may interpret that as a compression of a droplet. In case of Brownian motion, we deal with a definite decompression: particles are driven away from areas of higher concentration (probability of occurence). Hence, typically the Brownian ”pressure” is exerted by the droplet upon its surrounding.

Following the hydrodynamic tradition let us analyze that ”pressure” issue in more detail. We consider a reference volume (control interval, finite droplet) $`[\alpha ,\beta ]`$ in $`R^1`$ (or $`\mathrm{\Lambda }R^1`$ ) which at time $`t[0,T]`$ comprises a certain fraction of particles (Brownian ”fluid” constituents). The time rate of particles loss or gain by the volume $`[\alpha ,\beta ]`$ at time $`t`$, is equal to the flow outgoing through the boundaries i.e.

$$_t_\alpha ^\beta \rho (x,t)𝑑x=\rho (\beta ,t)v(\beta ,t)\rho (\alpha ,t)v(\alpha ,t)$$
(23)
which is a consequence of the continuity equation.

To analyze the momentum balance, let us allow for an infinitesimal deformation of the boundaries of $`[\alpha ,\beta ]`$ to have entirely compensated the mass (particle) loss or gain:

$$[\alpha ,\beta ][\alpha +v(\alpha ,t)\mathrm{}t,\beta +v(\beta ,t)\mathrm{}t]$$
(24)
Effectively, we pass then to the locally co-moving frame. That implies

$$lim_{\mathrm{}t0}\frac{1}{\mathrm{}t}\left[_{\alpha +v_\alpha \mathrm{}t}^{\beta +v_\beta \mathrm{}t}\rho (x,t+\mathrm{}t)𝑑x_\alpha ^\beta \rho (x,t)𝑑x\right]=$$
(25)
$$=lim_{\mathrm{}t0}\frac{1}{\mathrm{}t}\left[_{\alpha +v_\alpha \mathrm{}t}^\alpha \rho (x,t)𝑑x+\mathrm{}t_\alpha ^\beta (_t\rho )𝑑x+_\beta ^{\beta +v_\beta \mathrm{}t}\rho (x,t)𝑑x\right]=0$$

Let us investigate what happens to the local matter flows $`(\rho v)(x,t)`$, if we proceed in the same way (leading terms only are retained):

$$_{\alpha +v_\alpha \mathrm{}t}^{\beta +v_\beta \mathrm{}t}(\rho v)(x,t+\mathrm{}t)𝑑x_\alpha ^\beta (\rho v)(x,t)𝑑t$$
(26)
$$(\rho v^2)(\alpha ,t)\mathrm{}t+(\rho v^2)(\beta ,t)\mathrm{}t+\mathrm{}t_\alpha ^\beta [_t(\rho v)]𝑑x$$

In view of local conservation laws we have $`_t(\rho v)=(\rho v^2)+\rho (\mathrm{\Omega }Q)`$ and the rate of change of momentum associated with the control volume $`[\alpha ,\beta ]`$ is (here per unit of mass)

$$lim_{\mathrm{}t0}\frac{1}{\mathrm{}t}\left[_{\alpha +v_\alpha \mathrm{}t}^{\beta +v_\beta \mathrm{}t}(\rho v)(x,t+\mathrm{}t)_\alpha ^\beta (\rho v)(x,t)\right]=_\alpha ^\beta \rho (\mathrm{\Omega }Q)dx$$
(27)
However, $`Q=\frac{P}{\rho }`$ and $`P=D^2\rho \mathrm{}ln\rho `$. Therefore:

$$_\alpha ^\beta \rho (\mathrm{\Omega }Q)dx=_\alpha ^\beta \rho \mathrm{\Omega }dx_\alpha ^\beta Pdx=$$
(28)
$$=E[\mathrm{\Omega }]_\alpha ^\beta +P(\alpha ,t)P(\beta ,t)$$

Clearly, $`\mathrm{\Omega }`$ refers to the Euler-type volume force, while $`Q`$ (or more correctly, $`P`$) refers to the ”pressure” effects entirely due to the particle transfer rate through the boundaries of the considered volume.

The missing ingredient of our discussion is the time developement of the kinetic energy of the matter flow $`\frac{1}{2}(\rho v^2)`$ (per unit of mass) transported through the chosen volume. Let us therefore evaluate:

$$_{\alpha +v_\alpha \mathrm{}t}^{\beta +v_\beta \mathrm{}t}\frac{1}{2}(\rho v^2)(x,t+\mathrm{}t)𝑑x_\alpha ^\beta \frac{1}{2}(\rho v^2)(x,t)𝑑t$$
(29)
$$\frac{1}{2}(\rho v^3)(\alpha ,t)\mathrm{}t+\frac{1}{2}(\rho v^3)(\beta ,t)\mathrm{}t+\mathrm{}t_\alpha ^\beta [_t\frac{1}{2}(\rho v^2)]𝑑x$$

We have $`_t(\frac{1}{2}\rho v^2)=\frac{1}{2}v^2(\rho v)\rho v(Q\mathrm{\Omega }+\frac{1}{2}v^2)=(\frac{1}{2}\rho v^3)\rho v(Q\mathrm{\Omega })`$. Consequently, the time rate of the kinetic energy (of the flow) loss/gain by the volume reads:

$$lim_{\mathrm{}t0}\frac{1}{\mathrm{}t}\left[_{\alpha +v_\alpha \mathrm{}t}^{\beta +v_\beta \mathrm{}t}\frac{1}{2}(\rho v^2)(x,t+\mathrm{}t)_\alpha ^\beta \frac{1}{2}(\rho v^2)(x,t)\right]=_\alpha ^\beta (\rho v)(\mathrm{\Omega }Q)dx$$
(30)
In the integrand one immediately recognizes an expression (up to the mass parameter) for the power released/absorbed by the control volume at time $`t`$.

In particular, by taking advantage of the identity $`Q=\frac{P}{\rho }`$ we can rewrite the pure ”pressure” contribution as follows: $`_\alpha ^\beta \rho vQdx=_\alpha ^\beta vPdx`$ which clearly is a direct analog of the standard mechanical expression for the power release ($`\frac{dE}{dt}=Fv`$).

Surely, one cannot interpret the above local averages as an outcome of an innocent operation of the random medium upon particles. There is a non-negligible transfer of energy and momentum to be accounted for in association with the Brownian motion, and a number of problems (when can we disregard temperature gradients ?) pertaining to local heating and cooling phenomena suffered by the environment (in terms of local averages) should be consistently resolved.

## 5 Implementing the feedback

In contrast to molecular chaos derivations based on the Boltzmann equation, , in case of the Brownian motion we have no access into the microscopic fine details of particle propagation. Instead, we need to postulate the mathematically reliable form of the dynamics, even if ending up with obvious artifacts of the formalism (like e.g. non-differentiable paths, or their infinite spatial variation on arbitrarily short time intervals, in case of the Wiener process or, even worse, white-noise input).

Our major hypothesis about an exotic particle-bath coupling pertains to local averages and is motivated by a fairly intuitive picture of the behaviour of suspended particles in thermally inhomogeneous media. Namely, it is well known that hot areas support much lower density of suspended particles (can even be free of any dust admixtures) than the lower temperature areas. Dynamically we can interpret this phenomenon as a repulsion of suspended particles by warm areas and an attraction by the cool ones, cf. Ref. .

In the course of our discussion we have spent quite a while on demonstrating that a non-trivial energy/heat exchange is completely ignored in the standard approach to the Brownian motion (even, if under suitable circumstances one has good reasons to do that). Effectively, the model dramatically violates basic conservation laws of physics. That derives from the assumption that the thermostat is perpetually in the state of rest (in the mean), in thermal equilibrium and thus free of any thermal currents.

The force/acceleration term (we turn back to the three-dimensional notation) $`\stackrel{}{}(\mathrm{\Omega }Q)`$ appears in all formulas that refer to energy/momentum flows supported by local averages of the Brownian motion. Their sole reason is the action of the random medium upon particles which causes an expansion of the Brownian ”swarm” out of areas of higher concentration. That however needs a local cooling of the medium (Brownian particles are being ”thermalized” by increasing their mobility e.g. temperature). On sufficiently low time scales that should amount to an instantaneous reaction of the medium in terms of thermally induced currents: Brownian particles are driven back (attraction !) to the cooler areas i.e. float down the temperature gradients as long as they are non-vanishing, .

We shall bypass the thermal inhomogeneity issue by resorting to an explicit energy-momentum balance information available through local conservation laws for the Brownian motion.

If we regard the term $`\stackrel{}{}(\mathrm{\Omega }Q)`$ as a quantification of the sole medium action upon particles, then the most likely quantification of the medium reaction should be exactly the opposite i.e. $`\stackrel{}{}(Q\mathrm{\Omega })`$.

Told otherwise, medium ”effort” to release momentum and energy from the Brownian ”droplet” (control volume), at a time rate determined by the functional form of $`\stackrel{}{}(\mathrm{\Omega }Q)(x,t)`$, induces a compensating (in view of the heat deficit) energy and momentum delivery to that volume needed to remove the thermal inhomogeneity of the thermostat. As a consequence, Brownian particles propagate through the medium which is non longer in the state of rest and develops intrinsic (mean, on the ensemble average) flows.

A mathematical encoding of this Brownian recoil principle , or third Newton law in the mean  hypothesis is rather well established.

The momentum conservation law for the process with a recoil (the reaction term replaces the decompressive action term) will read:

$$_t\stackrel{}{v}+(\stackrel{}{v}\stackrel{}{})\stackrel{}{v}=\stackrel{}{}(Q\mathrm{\Omega })$$
(31)
implying that

$$_tS+\frac{1}{2}|\stackrel{}{}S|^2Q=\mathrm{\Omega }$$
(32)
stands for the corresponding Hamilton-Jacobi equation, cf. , instead of ”normal” one. A suitable adjustment (re-setting) of the initial data is here necessary, cf. .

In the coarse-grained picture of motion we shall deal with a sequence of repeatable feedback scenarios realized on the Smoluchowski process time scale: the Brownian ”swarm” expansion build-up is accompanied by the parallel counter-flow build-up, which in turn modifies the subsequent stage of the Brownian ”swarm” migration (being interpreted to modify the forward drift of the process) and the corresponding, built-up anew counter-flow.

Perhaps surprisingly, we are still dealing with Markovian diffusion-type processes, . The link is particularly obvious  if we observe that the new Hamilton-Jacobi equation can be formally rewritten in the previous form by introducing:

$$\mathrm{\Omega }_r=_tS+\frac{1}{2}|\stackrel{}{}S|^2+Q$$
(33)
where $`\mathrm{\Omega }_r=2Q\mathrm{\Omega }`$ and $`\mathrm{\Omega }`$ represents the previously defined potential function characterizing any Smoluchowski (or more general) diffusion process. It is however $`\mathrm{\Omega }_r`$ which would determine forward drifts of the Markovian diffusion process with a recoil. Those must come out from the Cameron-Martin-Girsanov identity $`\mathrm{\Omega }_r=2Q\mathrm{\Omega }=2D[_t\varphi +\frac{1}{2}(\frac{\stackrel{}{b}^2}{2D}+\stackrel{}{}\stackrel{}{b})]`$.

After complementing the Hamilton-Jacobi-type equation by the continuity equation, we again end up with a closed system of conservation laws. The system is badly nonlinear and coupled, but its linearisation can be immediately given in terms of an adjoint pair of Schrödinger equations with a potential $`\mathrm{\Omega }`$ (the imaginary unit $`i`$ on the left-hand-side in below is not an error !), . Indeed,

$$i_t\psi =D\mathrm{}\psi +\frac{\mathrm{\Omega }}{2D}\psi $$
(34)
with a solution represented in the polar form $`\psi =\rho ^{1/2}exp(iS)`$ and its complex adjoint makes the job. The choice of $`\psi (\stackrel{}{x},0)`$ sets here a solvable Cauchy problem. Notice that, in view of the Schrödinger-type linearization, for time-indepedent $`\mathrm{\Omega }`$ (conservative forces) the total energy $`_{R^3}(\frac{v^2}{2}Q+\mathrm{\Omega })\rho d^3x`$ of the system is a conserved finite quantity. We thus reside within the framework of so-called finite energy diffusion processes, whose mathematical features received some attention in the literature. In particular, it is known that in the absence of volume forces, a superdiffusion appears, while harmonic volume forces allow for non-dispersive diffusion-type processes, .