Diffusion Phenomena in a Mixed Phase Space

TL;DR

This paper analytically derives average particle velocities in strongly chaotic systems using Brownian dynamics and diffusion equations, validated on a simplified Fermi-Ulam model with mixed phase space.

Contribution

It introduces a novel analytical method for calculating average velocities in chaotic systems with mixed phase space, extending to time-dependent billiards.

Findings

Analytical velocities match numerical simulations.

Method applicable to various time-dependent billiard systems.

Provides a new approach for understanding diffusion in chaotic dynamics.

Abstract

We show that, in strongly chaotic dynamical systems, the average particle velocity can be calculated analytically by consideration of Brownian dynamics in phase space, the method of images and use of the classical diffusion equation. The method is demonstrated on the simplified Fermi-Ulam accelerator model, which has a mixed phase space with chaotic seas, invariant tori and Kolmogorov-Arnold-Moser (KAM) islands. The calculated average velocities agree well with numerical simulations and with an earlier empirical theory. The procedure can readily be extended to other systems including time-dependent billiards.

Equations

$$T:\left\{\begin{array}[]{ll}V_{n+1}=|V_{n}-2\varepsilon\sin(\phi_{n+1})|{}{}\\ \phi_{n+1}=[\phi_{n}+\frac{2}{V_{n}}]\mod(2\pi)\\ \end{array}\right..$$

$$T:\left\{\begin{array}[]{ll}I_{n+1}=I_{n}+K\sin(\theta_{n}),{}{}\\ \theta_{n+1}=[\theta_{n}+I_{n}],\mod(2\pi)\\ \end{array}\right.$$

$$\rho (V,t)=\frac{1}{\sqrt{4\pi Dt}}{e}^{\frac{-(V-V_0{)}^2}{4Dt}}.$$

$$\rho (V;\mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{\frac{-(V-\mu {)}^2}{2{\sigma }^2}}.$$

$$\begin{gathered} \rho (V;\mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\sum _{m=-\infty }^{\infty }[\exp \left(\frac{-(V-2m{V}_{\mathrm{f}}-\mu {)}^2}{2{\sigma }^2}\right) \\ +\exp \left(\frac{-(V-2m{V}_{\mathrm{f}}+\mu {)}^2}{2{\sigma }^2}\right)]. \end{gathered}$$

$$\begin{gathered} \rho (V;\mu ,{\sigma }^2)=\frac{1}{2A\sqrt{2\pi {\sigma }^2}}\sum _{m=-\infty }^{\infty }[\exp \left(\frac{-(V-2m{V}_{\mathrm{f}}-\mu {)}^2}{2{\sigma }^2}\right) \\ +\exp \left(\frac{-(V-2m{V}_{\mathrm{f}}+\mu {)}^2}{2{\sigma }^2}\right)], \end{gathered}$$

$$\begin{gathered} ⟨V⟩=\frac{1}{2A\sqrt{\pi }}\sum _{m=-\infty }^{\infty }\sqrt{2{\sigma }^2}({\Delta }^{(1)}\text{exp}+{\Delta }^{(2)}\text{exp}) \\ +\sqrt{\pi }\left[(2m{V}_{\mathrm{f}}-\mu ){\Delta }^{(1)}\text{erf}+(2m{V}_{\mathrm{f}}+\mu ){\Delta }^{(2)}\text{erf}\right], \end{gathered}$$

$${\Delta }^{(1)}\text{exp}$$

$${\Delta }^{(2)}\text{exp}$$

$${\Delta }^{(1)}\text{erf}$$

$${\Delta }^{(2)}\text{erf}$$

$$\begin{gathered} ⟨V⟩=\frac{\mu }{2A}\sum _{m=-\infty }^{\infty }\frac{1}{z\sqrt{\pi }}({\Delta }^{(1)}\text{exp}+{\Delta }^{(2)}\text{exp}) \\ +\frac{1}{\mu }\left[\left(\frac{2m{V}_{\mathrm{f}}}{\sqrt{2{\sigma }^2}}-\mu \right){\Delta }^{(1)}\text{erf}+\left(\frac{2m{V}_{\mathrm{f}}}{\sqrt{2{\sigma }^2}}+\mu \right){\Delta }^{(2)}\text{erf}\right], \end{gathered}$$

$${\Delta }^{(1)}\text{exp}$$

$${\Delta }^{(2)}\text{exp}$$

$${\Delta }^{(1)}\text{erf}$$

$${\Delta }^{(2)}\text{erf}$$

$$\begin{array}{rl}{\mu }_{n+1}& =⟨V_{n+1}⟩=⟨V_n⟩\Rightarrow \mu ={\mu }_0=V_0\\ {\sigma }_{n+1}^2& =⟨V_{n+1}^2⟩-⟨V_{n+1}{⟩}^2=⟨V_n^2⟩+2{\epsilon }^2-⟨V_n{⟩}^2\\ & \Rightarrow {\sigma }_{n+1}^2={\sigma }_n^2+2{\epsilon }^2.\end{array}$$

$${\sigma }_{n+1}^2-{\sigma }_n^2=\frac{d{\sigma }^2}{dn}\Rightarrow {\sigma }^2(n)={\sigma }_0^2+2{\epsilon }^{2}n.$$

$$z=\frac{\mu }{\sqrt{2}\sigma (n)}\Rightarrow z(n)=\frac{V_0}{2\epsilon \sqrt{n}}.$$

$$z=1\Rightarrow \frac{V_0}{2\epsilon \sqrt{n}}=1\Rightarrow n={\left(\frac{V_0}{2\epsilon }\right)}^2\Rightarrow n=\frac{V_0^2}{4{\epsilon }^2},$$

$$n_x=\frac{V_0^2}{4{\epsilon }^2}.$$

$$\tilde{v}=1\Rightarrow \frac{V_{\mathrm{f}}}{2\epsilon \sqrt{n}}=1\Rightarrow n\approx {\left(\frac{2\sqrt{\epsilon }}{2\epsilon }\right)}^2\Rightarrow n\approx \frac{1}{\epsilon },$$

$$n_x^{\prime }\approx \frac{1}{\epsilon }.$$

$$\underset{\sigma \to \infty }{\lim }⟨V⟩=\underset{n\to \infty }{\lim }⟨V⟩=\frac{V_{\mathrm{f}}}{2}\approx \sqrt{\epsilon },$$

Full text

We present an approach to the analysis of complicated chaotic dynamical systems that are difficult to treat in other ways. It relies on the counter-intuitive application of a fundamental idea from classical continuum physics (probabilistic diffusion) to chaotic systems that are, of course, inherently deterministic. In particular, we consider diffusion and Brownian dynamics in the phase space of the chaotic system and show how the diffusion equation, applied in this unusual context, can provide an accurate description of the average velocity and its evolution. To demonstrate and validate the formalism, we take a well-known example from astrophysics - the Fermi-Ulam model.

I Introduction

The evolution of systems described by Hamiltonians with nonlinear terms in their dynamical equations may exhibit either regularity or chaos. The result is often a mixed phase space containing chaotic seas, invariant tori and Kolmogorov-Arnold-Moser (KAM) islands lich . Dynamical systems with strong chaotic motion often exhibit diffusive behavior ott ; chirikov . An intuitive example of this is to drop colored ink into water, observing how the particles of ink move away from each other, spreading out into the liquid. For a mixed phase space, however, an initial condition e.g. around a KAM island may lead to very complicated behavior. The stability structures influence directly the transport properties of chaotic orbits venegeroles1 , often generating so-called anomalous diffusion zasla ; robnik .

There are many scenarios where, rather than analysing the individual behaviour of a single particle starting from a particular initial condition, it is more interesting to consider the average properties of the system, taking into account an ensemble of particles. Statistical methods can then be used to describe the dynamical phenomena venegeroles2 ; meiss ; manchein . Correspondingly, the properties and construction of the phase space can lead to what are effectively diffusion processes: as the dynamics evolves, there is diffusion of the action, usually associated with the velocity of the particles, through the phase space.

In this work we show that the classical diffusion equation dif1 ; dif2 ; dif3 ; dif4 ; dif5 can be solved via a procedure well-known in electrostatics, namely the method of images, and used to describe the evolution of the average velocity for a system characterised by a mixed phase space. We will demonstrate the effectiveness and utility of this idea by applying it to the well-known and widely-studied Fermi-Ulam model (FUM).

This paper is organised as follows. In Sec. II we describe the FUM, showing the nonlinear map associated with the dynamics and introducing a picture of a diffusion process occurring within its characteristic phase space. Section III develops a theoretical framework yielding analytical results for normal diffusion in a mixed phase space. In Sec. IV we compare these analytical results with numerical data. Conclusions are drawn in Sec. V.

II Model and phase space

The FUM fum is a version of the Fermi accelerator, which was originally introduced by Enrico Fermi fermi as a possible explanation for the production of very high energy cosmic rays. Its acceleration mechanism involves the repulsion of an electrically charged particle by strong oscillatory magnetic fields, a process that is analogous to a classical particle colliding with an oscillating physical boundary. The model consists of a particle bouncing back and forth between two rigid walls, one of which is fixed, whereas the other moves periodically in time with a normalized amplitude $\varepsilon$, as shown schematically in Fig.1.

The system is described by a two-dimensional, nonlinear, area-preserving map $T(V_{n},\phi_{n})=(V_{n+1},\phi_{n+1})$. The velocity of the particle is the action variable and the phase, related to the time-dependent boundary, is the angle variable. Taking into account that the absolute value of velocity changes at the moment of each collision, the mapping for the simplified version 111We approximate the oscillating wall as fixed. In reality, of course, when the particle suffers a collision, it exchanges momentum appropriate to a moving wall, but this simplified version is valid when the nonlinear parameter $\varepsilon$ is relatively small. of the FUM is

[TABLE]

The term $\frac{2}{V_{n}}$ corresponds to the time between collisions and $-2\varepsilon\sin(\phi_{n+1})$ gives the gain or loss of velocity/energy in each collision.

The phase space $V\times\phi$ for the FUM is composed of chaotic seas and KAM islands, and is accordingly classified as a mixed phase space. In addition, it is bounded by an invariant spanning curve which plays the role of a boundary: trajectories of lower velocity will never visit a region above this curve, no matter how many times the trajectory is iterated.

The average velocity of an ensemble of particles inside the FUM grows initially tiago , and then flattens off towards a plateau. This velocity growth and saturation can be interpreted as involving a diffusion process, albeit diffusion not in the physical space of the FUM, but rather in its phase space. Figure 2 shows how this phase-space-diffusion behaves for different numbers of iterations $n$. At $n=0$ we have the initial Gaussian-shaped distribution centered at $[V_{0}=0.01,\phi_{0}=\pi]$; then, one iteration latter, the distribution seems to have become spread out uniformly along the phase axis, a fact that will be used later in the analytic approach. However, diffusion is also starting on the action axis. After $10$ and then $100$ iterations of the mapping Eq. (1), the action/velocity is still continuing its diffusion through phase space. For all panels of Fig. 2, $\varepsilon=0.001$.

It is important to bear in mind that the phase space diffusion is limited down by null velocity and up by the first invariant spanning curve. Its position is approximated by $V_{\rm f}\approx 2\sqrt{\varepsilon}$. The localization of such a curve can be obtained by using a connection with the standard mapping lich ; vfisc , which is written as

[TABLE]

where the parameter $K$ controls the intensity of the nonlinearity of the mapping. There are two transitions in the standard mapping: (i) integrability when $K=0$ to non-integrability for any $K\neq 0$; and (ii) a transition from local chaos when $K<K_{c}$ to global chaos for $K>K_{c}$. The parameter $K_{c}=0.9716\ldots$ identifies the critical value of control parameter where all of the invariant spanning curves are destroyed, letting the dynamics diffuse unbounded in the $I$ direction. This is exactly the transition we want to use in connection with the FUM as an attempt to describe the localization of the first invariant spanning curve. Above the curve in the FUM, one observes local chaos, an infinity of other invariant spanning curves and eventually periodic orbits. Below the first invariant spanning curve only chaos, periodic and quasi periodic dynamics coexist, each one of them being visited as determined by the initial conditions. The procedure to obtain $V_{\rm f}$ consists of describing the position of the first invariant spanning curve in the FUM through a local description of the standard mapping. Then a Taylor expansion (see Ref. vfisc for more details in a family of area preserving mappings) is made in the first equation of mapping (1) by using the fact that the invariant spanning curve is written as $V_{n}=V_{\rm f}+\Delta V_{n}$ where $\Delta V_{n}\ll V_{\rm f}$ is a small perturbation of the typical value $V_{\rm f}$. A first order approximation leads to the expression $V_{\rm f}\approx 2\sqrt{\varepsilon}$.

III Analytical procedure

Essentially, the action variable $V$ is undergoing a diffusion process within the bounded space $V\in[0,V_{\rm f}]$. This can be described by the diffusion equation with no flux through its boundaries $\frac{\partial\rho(0,t)}{\partial V}=\frac{\partial\rho(V_{\rm f},t)}{\partial V}=0,\forall t>0$. Thus, the problem may be reduced to that of solving the diffusion equation to obtain the probability density function $\rho(V,t)$; once this has been integrated along the bounded space $\langle V\rangle=\int_{0}^{V_{\rm f}}V\rho(V,t)dV$, it yields a theoretical prediction for the average velocity of $\langle V\rangle$ of particles inside the FUM.

The solution of the diffusion equation with no flux through the boundaries can be obtained analytically by the method of images, as in electrostatics crank . Basically the idea is to treat the initial Gaussian distribution as a point charge and the boundaries as conducting planes. The solution will then be an infinite sum of Gaussian functions centered at $V_{0}$, due to the infinity of images of the initial profile.

First, we consider a normal diffusion process in one dimension, with no boundaries and with the initial condition $\rho(V,0)=\delta(V-V_{0})$, where $V_{0}$ is the initial velocity of the particles. The fundamental solution of the diffusion equation is given by

[TABLE]

It is a Gaussian function (normal distribution). Solutions of this type are well-known and widely applicable in science gauss1 ; gauss2 ; gauss3 , especially, in statistical physics. Normal distributions are characterized by their mean value $\mu\equiv\langle V\rangle$ and variance $\sigma^{2}\equiv\langle V^{2}\rangle-\langle V\rangle^{2}$. Likewise, the diffusion coefficient can be written as a function of the time derivative of the variance $D=\frac{1}{2}\frac{d\sigma^{2}}{dt}~{}~{}\rightarrow~{}~{}\sigma^{2}=2Dt$. The solution can then be rewritten in terms of $\mu$ and $\sigma^{2}$ as

[TABLE]

Knowing the fundamental solution, and applying the principle of superposition as in the method of images, we may assume that a sum of Gaussian functions is still a solution to the problem. Hence the solution when $V\in[0,V_{\rm f}]$ with $\frac{\partial\rho(0,t)}{\partial V}=\frac{\partial\rho(V_{\rm f},t)}{\partial V}=0$ is given by

[TABLE]

As it stands, however, this solution is not normalized for the space interval $V\in[0,V_{\rm f}]$. To effect normalization, it is necessary that $A\int_{0}^{V_{\rm f}}\rho(V;\mu,\sigma^{2})dV~{}=1$, with $A$ equal to a normalization constant. Considering the error function property $\text{erf}(-x)=-\text{erf}(x)$, the normalized solution is given by

[TABLE]

with $A=\sum_{m=-\infty}^{\infty}\text{erf}\left(\frac{\mu-2mV_{\rm f}}{\sqrt{2\sigma^{2}}}\right)-\text{erf}\left(\frac{\mu-V_{\rm f}-2mV_{\rm f}}{\sqrt{2\sigma^{2}}}\right)-\text{ erf}\left(\frac{\mu+2mV_{\rm f}}{\sqrt{2\sigma^{2}}}\right)+\text{erf}\left(\frac{\mu+V_{\rm f}+2mV_{\rm f}}{\sqrt{2\sigma^{2}}}\right)$.

Fig. 3 shows how the analytical solution for the probability density given by Eq. (6) fits the numerical simulation data for the FUM. The initial conditions for the analytic curve are the same as those used in constructing Fig. 2, with $V_{0}=0.01$, $\sigma_{V_{0}}=0.001$ and $\varepsilon=0.001$. This comparison robnik2 provides a convincing verification of the analytical solution. The good fit indicates that the solution is suitable when considering an initial profile in a chaotic region and neglecting the anomalous diffusion phenomena around KAM islands.

Having obtained this solution, we need to calculate the average as $\langle V\rangle=\int_{0}^{V_{\rm f}}x\rho(V;\mu,\sigma^{2})dV$ in order to be able to predict analytically the average behaviour of the velocity. Because of the lack of symmetry, this calculation is non-trivial but, integrating between the upper and lower limits using the Jacobi Theta function representation jacobi , we find that the solution can be written as

[TABLE]

where

[TABLE]

Defining an important auxiliary variable $z=\frac{\mu}{\sqrt{2\sigma^{2}}}$ and a new parameter $\tilde{v}=\frac{V_{\rm f}}{\sqrt{2\sigma^{2}}}$ making the necessary re-arrangements, the average velocity within the FUM is given analytically by

[TABLE]

with $A=\sum_{m=-\infty}^{\infty}\text{erf}\left(z-\frac{2mV_{\rm f}}{\sqrt{2\sigma^{2}}}\right)-\text{erf}\left(z-\tilde{v}-\frac{2mV_{\rm f}}{\sqrt{2\sigma^{2}}}\right)-\text{erf}\left(z+\frac{2mV_{f}}{\sqrt{2\sigma^{2}}}\right)+\text{erf}\left(z+\tilde{v}+\frac{2mV_{\rm f}}{\sqrt{2\sigma^{2}}}\right)$ and

[TABLE]

Furthermore, the mean $\mu$ and variance $\sigma^{2}$ are calculated, by construction, over the point charge, which is characterised by an unbounded diffusion process. According to our initial mapping, Eq. (1), the point charge mapping is given by $V_{n+1}=V_{n}-2\varepsilon\sin(\phi)$, where $\phi$ is an uniform random variable, as observed in Fig. 2, it is then possible to write the mean and variance for the point charge as

[TABLE]

Following the theory of difference equations aprox , assuming a large number of iterations and small values of $\varepsilon$, it is then possible to write $\sigma$ as a function of $n$

[TABLE]

This result is important because it carries the information that the variance is a function of the number of iterations $\sigma^{2}=\sigma^{2}(n)$, connecting the solution of the diffusion equation to the discrete mapping of the FUM. Moreover, the initial variance $\sigma_{0}$ is zero if the initial profile is considered a perfect Dirac delta function. This also tells us the diffusion is normal, since $\sigma\propto\sqrt{n}$. In addition, it is also possible to calculate the diffusion coefficient, which is a constant and quite intuitive with our suppositions for this case, so that $D=\epsilon^{2}$. Then $z$ is also a function of the number of iterations $n$ such that

[TABLE]

Substituting Eq. (9) into Eq. (8), we can calculate how the average velocity behaves as a function of the number of iterations for the dynamics of the FUM. But of course we now need to check whether, or not, this theory really describes the actual behaviour of the average velocity.

IV Analytical $\times$ Numerical results

Fig. 4 compares the numerical simulation data with the analytic predictions of Eqs. (8),(9). The expression for $\langle V\rangle$, given by Eq. (8), represents a continuous competition between the exponential and error functions, so it is interesting to study their arguments. Based on a graphical analysis, we conclude that there are two changes of behaviour: at $z=1$; and at $\tilde{v}=1$. First, taking $z=1$

[TABLE]

but in this case, $n=n_{x}$ marking the first crossover 222The crossover number is defined as the iteration number that marks some dynamical change, for example the change from growth to saturation.. Thus

[TABLE]

Secondly, taking $\tilde{v}=1$

[TABLE]

but now, $n=n^{\prime}_{x}$ marking the second crossover. Thus

[TABLE]

Another important result is the limit

[TABLE]

which provides the saturation value $V_{sat}$. Then, Fig. 4 shows the average velocity for an ensemble of $10^{3}$ particles, all with initial velocity $V_{0}=2\times 10^{-3}$, taken within the interval $\phi_{0}\in[0,2\pi]$. The analytical predictions for the first crossover $n_{x}$, Eq. (10), the second crossover $n^{\prime}_{x}$, Eq. (11), and the saturation plateau when $n\rightarrow\infty$, Eq. (12), are shown by the dashed lines.

The analytical approach, yielding Eq. (8), clearly agrees well with the numerical simulation data. The correspondence might have been even closer were if not for the fact that the diffusion is not ideal for higher values of $V$, due to the configuration of the phase space. This also explains the fluctuation for $n>n^{\prime}_{x}$. The diffusion around stability structures like KAM islands leads to the very complicated behaviour known as anomalous diffusion. However, the associated stickiness of the dynamics near the islands, though real, is a relatively minor effect given the size of the whole phase space: Harsoula et al harsoula conclude that, for a long enough interval, averaging over the ensemble smooths the observables so that the stickiness can largely be neglected.

Fig. 5 shows compares numerical data with the corresponding analytical predictions for three different initial velocities $V_{0}$ and values of the control parameter $\varepsilon$. It is important to remember that the position of the upper boundary in the phase space, which is the first invariant spanning curve, is approximated by $V_{\rm f}\approx 2\sqrt{\varepsilon}$. Then, for each value of the parameter $\varepsilon$, a different bounded phase space is considered. Again, it is evident that the analytic curves provide an excellent fit to the numerical data, even for relatively large values of $\varepsilon$.

We emphasize that Eqs. (10, 11, 12) represent the first analytic predictions to be made for the Fermi-Ulam model. They agree well with what was proposed on purely empirical grounds prl2004 more than a decade ago. Three hypotheses were then proposed, based on a scaling analysis: (i) $n_{x}\propto\frac{V_{0}^{2}}{\varepsilon^{2}}$ which agrees perfectly with Eq. (10), and we now also obtain the proportionality constant $\frac{1}{4}$; in addition (ii) $n^{\prime}_{x}\propto\frac{1}{\varepsilon}$ which agrees with Eq. (11); and finally (iii) $V_{sat}\propto\varepsilon^{\alpha}$, with $\alpha\approx\frac{1}{2}$, which agrees with Eq. (12).

V Conclusion

We conclude that a combination of the theory of diffusive processes with dynamical systems theory, plus the method of images from electrostatics, provides a powerful method for treating systems described by nonlinear mappings. The method can be expected to work for mixed phase spaces that are delimited by boundaries through which there are no fluxes. Application to the Fermi-Ulam model, taken as an example, has yielded some interesting features and excellent agreement both with numerical simulations and with earlier empirically-based theoretical considerations. Extension of the procedure discussed here to time-dependent billiards gelfreich is an interesting possibility for future work.

Acknowledgements.

We gratefully acknowledge valuable discussions with Professor Roberto Lagos. M. S. Palmero was supported by Fundação de Amparo à Pesquisa do Estado de São Paulo FAPESP, from Brazil, processes number 2014/27260-5 and 2016/15713-0. G. D. Iturry thanks to Brazilian agency CNPq. The research was supported by the Engineering and Physical Sciences Research Council (United Kingdom) under grants GR/R03631 and EP/M015831/1. E. D. Leonel acknowledges support from CNPq (303707/2015-1) and FAPESP (2017/14414-2).