Diffusion Phenomena in a Mixed Phase Space
Matheus S. Palmero Gabriel I. Díaz Peter V. E. McClintock Edson D. Leonel
cond-mat.stat-mech nlin.CD













^1Departamento de Física, UNESP - Univ Estadual Paulista, Av. 24A, 1515, Bela Vista, 13506-900, Rio Claro, SP, Brazil

^2 Department of Physics, Lancaster University, Lancaster, LA1 4YW, United Kingdom

^3 Instituto de Física, IFUSP - Universidade de São Paulo, Rua do Matão, Tr.R 187, Cidade Universitária, 05314-970, São Paulo, SP, Brazil




05.45.-a, 05.45.Jc, 47.52.+j



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.


Diffusion Phenomena in a Mixed Phase Space
    Matheus S. Palmero^1,2,3, Gabriel I. Díaz^3, Peter V. E. McClintock^2 and Edson D. Leonel^1
    
===============================================================================================





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.





§ 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 <cit.>. Dynamical
systems with strong chaotic motion often exhibit diffusive behavior <cit.>.
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 <cit.>, often generating
so-called anomalous diffusion <cit.>.

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
<cit.>. 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
<cit.> 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. <ref> 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 <ref> develops a theoretical
framework yielding analytical results for normal diffusion in a mixed
phase space. In Sec. <ref> we compare these analytical results
with numerical data. Conclusions are drawn in Sec. <ref>.




§ MODEL AND PHASE SPACE

The FUM <cit.> is a version of the Fermi accelerator, which was
originally introduced by Enrico Fermi <cit.> 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 ε, as shown
schematically in Fig.<ref>.





The system is described by a two-dimensional, nonlinear,
area-preserving map T(V_n,ϕ_n)=(V_n+1,ϕ_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 [We 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 ε is relatively small.] of the FUM is

    T:{[ V_n+1 =|V_n - 2εsin(ϕ_n+1)|;  ϕ_n+1 = [ϕ_n + 2/V_n] (2π);                             ]..
The term 2/V_n corresponds to the time between collisions and
- 2εsin(ϕ_n+1) gives the gain or loss of
velocity/energy in each collision.

The phase space V ×ϕ 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 <cit.>, 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 <ref> 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,ϕ_0=π]; 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. (<ref>), the action/velocity is still continuing
its diffusion through phase space. For all panels of Fig. <ref>,
ε=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_ f≈2√(ε).
The localization of such a curve can be
obtained by using a connection with the standard mapping
<cit.>, which is written as

    T:{[    I_n+1 =I_n +Ksin(θ_n),; θ_n+1 = [θ_n + I_n], (2π);                           ].

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 K0; and (ii)
a transition from local chaos when K<K_c to global chaos for K>K_c.
The parameter K_c=0.9716… 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_ 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. <cit.>
for more details in a family of area preserving mappings) is made in
the first equation of mapping (<ref>) by using the fact that the invariant
spanning curve is written as V_n=V_ f+Δ V_n where Δ
V_n≪ V_ f is a small perturbation of the typical value
V_ f. A first order approximation leads to the expression
V_ f≈2√(ε).



§ ANALYTICAL PROCEDURE

Essentially, the action variable V is undergoing a diffusion process
within the bounded space V ∈ [0,V_ f]. This can be described by the
diffusion equation with no flux through its boundaries ∂ρ(0,t)/∂ V=∂ρ(V_ f,t)/∂ V=0,
∀ t>0. Thus, the problem may be reduced to that of solving the
diffusion equation to obtain the probability density function
ρ(V,t); once this has been integrated along the bounded space
⟨ V⟩= ∫_0^V_ f Vρ(V,t) dV, it yields a
theoretical prediction for the average velocity of ⟨ V⟩ 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 <cit.>. 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 ρ(V,0)=δ(V-V_0),
where V_0 is the initial velocity of the particles. The fundamental
solution of the diffusion equation is given by

    ρ(V,t)=1/√(4π Dt)e^-(V-V_0)^2/4Dt .
It is a Gaussian function (normal distribution). Solutions of this type
are well-known and widely applicable in science
<cit.>, especially, in statistical physics. Normal
distributions are characterized by their mean value μ≡⟨
V⟩ and variance σ^2 ≡⟨ V^2⟩ - ⟨
V⟩^2. Likewise, the diffusion coefficient can be written as a
function of the time derivative of the variance
D=1/2dσ^2/dt  →  σ^2=2Dt. The
solution can then be rewritten in terms of μ and σ^2 as

    ρ(V;μ,σ^2)=1/√(2πσ^2)e^-(V-μ)^2/2σ^2 .
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 ∈ [0,V_ f] with ∂ρ(0,t)/∂
V=∂ρ(V_ f,t)/∂ V=0 is given by


    ρ(V;μ,σ^2)=1/√(2πσ^2)∑_m=-∞^∞[.
    exp(-(V-2mV_ f-μ)^2/2 σ^2)
    
    +exp(-(V-2mV_ f+μ)^2/2σ^2).].

As it stands, however, this solution is not normalized for the space interval V ∈ [0,V_ f]. To effect normalization, it is necessary that A ∫_0^V_ fρ(V;μ,σ^2) dV  = 1, with A equal to a normalization constant. Considering the error function property erf(-x)=-erf(x), the normalized solution is given by


    ρ(V;μ,σ^2)=1/2A√(2πσ^2)∑_m=-∞^∞[.
    exp(-(V-2mV_ f-μ)^2/2 σ^2)
    
    +exp(-(V-2mV_ f+μ)^2/2σ^2).] ,

with 
A=∑_m=-∞^∞erf(μ-2mV_ f/√(2σ^2
))-erf(μ-V_ f-2mV_ f/√(2σ^2))-
erf(μ+2mV_ f/√(2σ^2))+erf(μ+V_ f+2mV_ f/√(2σ^2)).



Fig. <ref> shows how the analytical solution for the probability density given by Eq. (<ref>) fits the numerical simulation data for the FUM. The initial conditions for the analytic curve are the same as those used in constructing Fig. <ref>, with V_0=0.01, σ_V_0=0.001 and ε=0.001. This comparison <cit.> 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 ⟨ V⟩= ∫_0^V_ f xρ(V;μ,σ^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 <cit.>, we find that the solution can be written as

    ⟨ V⟩=1/2A√(π)∑_m=-∞^∞√(2σ^2)(Δ^(1)exp+Δ^(2)exp)
    
    + √(π)[(2mV_ f-μ)Δ^(1)erf+(2mV_ f+μ)Δ^(2)erf] ,
where



    Δ^(1)exp   =   exp(-(μ-2mV_ f/√(2σ^2))^2)-exp(-(μ-2mV_ f+V_ f/√(2σ^2))^2),
    Δ^(2)exp   =   exp(-(μ+2mV_ f/√(2σ^2))^2)-exp(-(μ+2mV_ f-V_ f/√(2σ^2))^2),
    Δ^(1)erf   =   erf(μ-2mV_ f-V_ f/√(2σ^2))-erf(μ-2mV_ f/√(2σ^2)),
    Δ^(2)erf   =   erf(μ+2mV_ f/√(2σ^2))-erf(μ+2mV_ f+V_ f/√(2σ^2)).

Defining an important auxiliary variable z=μ/√(2σ^2) and a new parameter ṽ=V_ f/√(2σ^2) making the necessary re-arrangements, the average velocity within the FUM is given analytically by

    ⟨ V⟩=μ/2A∑_m=-∞^∞1/z√(π)(Δ^(1)exp+Δ^(2)exp) 
    
    + 1/μ[(2mV_ f/√(2σ^2)-μ)Δ^(1)erf+(2mV_ f/√(2σ^2)+μ)Δ^(2)erf] ,
with A=∑_m=-∞^∞erf(z-2mV_ f/√(2σ^2))-erf(z-ṽ-2mV_ f/√(2σ^2))-erf(z+2mV_f/√(2σ^2))+erf(z+ṽ+2mV_ f/√(2σ^2)) and



    Δ^(1)exp   =    e^-(z-2mV_ f/√(2σ^2))^2-e^-(z-2mV_ f/√(2σ^2)+ṽ)^2,
    Δ^(2)exp   =    e^-(z+2mV_ f/√(2σ^2))^2-e^-(z+2mV_ f/√(2σ^2)-ṽ)^2,
    Δ^(1)erf   =   erf(z-2mV_f/√(2σ^2)+ṽ)-erf(z-2mV_ f/√(2σ^2)),
    Δ^(2)erf   =   erf(z+2mV_ f/√(2σ^2))-erf(z+2mV_ f/√(2σ^2)+ṽ).



Furthermore, the mean μ and variance σ^2 are calculated, by construction, over the point charge, which is characterised by an unbounded diffusion process. According to our initial mapping, Eq. (<ref>), the point charge mapping is given by
V_n+1=V_n-2εsin(ϕ), where ϕ is an uniform random variable, as observed in Fig. <ref>, it is then
possible to write the mean and variance for the point charge as


    μ_n+1   =⟨ V_n+1⟩=⟨ V_n⟩⇒μ=μ_0=V_0
    σ^2_n+1   =⟨ V^2_n+1⟩ - ⟨ V_n+1⟩^2 = ⟨ V^2_n⟩ +2ε^2 - ⟨ V_n⟩^2 
       ⇒σ^2_n+1=σ^2_n + 2ε^2 .
Following the theory of difference equations <cit.>, assuming a large number of iterations and small values of ε, it is then possible to write σ as a function of n

    σ^2_n+1-σ^2_n=d σ^2/dn ⇒ σ^2(n)=σ^2_0 + 2ε^2 n .
This result is important because it carries the information that the variance is a function of the number of iterations σ^2=σ^2(n), connecting the solution of the diffusion equation to the discrete mapping of the FUM. Moreover, the initial variance σ_0 is zero if the initial profile is considered a perfect Dirac delta function. This also tells us the diffusion is normal, since σ∝√(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=ϵ^2. Then z is also a function of the number of iterations n such that




    z=μ/√(2)σ(n) ⇒  z(n)=V_0/2ε√(n) .


Substituting Eq. (<ref>) into Eq. (<ref>), 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.



§ ANALYTICAL × NUMERICAL RESULTS


Fig. <ref> compares the numerical simulation data with the analytic predictions of Eqs. (<ref>),(<ref>). The expression for ⟨ V⟩, given by Eq. (<ref>), 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 ṽ=1. First, taking z=1

    z=1  ⇒ V_0/2ε√(n)=1  ⇒  n=(V_0/2ε)^2  ⇒  n=V_0^2/4ε^2 ,
but in this case, n=n_x marking the first crossover [The crossover number is defined as the iteration number that marks some dynamical change, for example the change from growth to saturation.]. Thus

    n_x=V_0^2/4ε^2 .
Secondly, taking ṽ=1

    ṽ=1  ⇒ V_ f/2ε√(n)=1  ⇒  n≈(2√(ε)/2ε)^2  ⇒  n≈1/ε ,
but now, n=n'_x marking the second crossover. Thus

    n'_x≈1/ε .
Another important result is the limit

    lim_σ→∞⟨ V⟩ =lim_n →∞⟨ V⟩ = V_ f/2≈√(ε) ,
which provides the saturation value V_sat.
Then, Fig. <ref> shows the average velocity for an ensemble of 10^3 particles, all with initial velocity V_0=2×10^-3, taken within the interval ϕ_0∈[0,2π]. The analytical predictions for the first crossover n_x, Eq. (<ref>), the second crossover n'_x, Eq. (<ref>), and the saturation plateau when n→∞, Eq. (<ref>), are shown by the dashed lines.

The analytical approach, yielding Eq. (<ref>), 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'_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 <cit.> conclude that, for a long enough interval, averaging over the ensemble smooths the observables so that the stickiness can largely be neglected.



Fig. <ref> shows compares numerical data with the corresponding analytical predictions for three different initial velocities V_0 and values of the control parameter ε. 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_ f≈2√(ε). Then, for each value of the parameter ε, 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 ε.

We emphasize that Eqs. (<ref>, <ref>, <ref>) represent the first analytic predictions to be made for the Fermi-Ulam model. They agree well with what was proposed on purely empirical grounds <cit.> more than a decade ago. Three hypotheses were then proposed, based on a scaling analysis: (i) n_x ∝V_0^2/ε^2 which agrees perfectly with Eq. (<ref>), and we now also obtain the proportionality constant 1/4; in addition (ii) n'_x ∝1/ε which agrees with Eq. (<ref>); and finally (iii) V_sat∝ε^α, with α≈1/2, which agrees with Eq. (<ref>).



§ 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 <cit.> is an interesting possibility for future work.


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). 




 99
  lich A. J. Lichtenberg, M. A. Lieberman. Regular and Chaotic Dynamics, Springer (1992).

  ott E. Ott, Phys. Rev. Lett. 42, 1628 (1979).

  chirikov B. V.Chirikov. A universal instability of many-dimensional oscillator systems, Phys. Rep. 52, 263 (1979).

  venegeroles1 R. Venegeroles, Phys. Rev. Lett. 101, 054102 (2008).

  zasla G. M. Zaslasvsky, Physics of Chaos in Hamiltonian Systens, Imperial College Press (2007).

  robnik T. Manos, M. Robnik, Phys. Rev. E 89, 022905 (2014).

  venegeroles2 R. Venegeroles, Phys. Rev. Lett. 102, 064101 (2009).

  meiss O. Alus, S. Fishman and J. D. Meiss, Phys. Rev. E 96, 032204 (2017).

  manchein R. M. da Silva, M. W. Beims, C. Manchein, Phys. Rev. E 92, 022921 (2015).

  dif1 P. J. Basser, D. K. Jones, NMR in Biomedicine 15, 07 (2002).

  dif2 H. A. Kramers, Physica, 7 284 (1940).

  dif3 J. D. Meiss, Chaos 25, 097602 (2015).

  dif4 E. G. Altmann, J. S. E. Portela and T. Tél, Rev. Mod. Phys. 85, 869 (2013).

  dif5 G. M. Zaslavsky, Physics Reports 371, 461 (2002).

  fum M. A. Lieberman and A. J. Lichtenberg, Phys. Rev. A 5, 1852 (1971).

  fermi E. Fermi, Phys. Rev. 75, 1169 (1949).

  simpliver Suppose that the oscillating wall is fixed, however when the particle suffers a collision, it exchanges momentum as if the wall were moving. This simplified version is indeed valid when the nonlinear parameter ε is relatively small.

  tiago T. Pereira, D. Turaev, Phys. Rev. E 91, 010901 (2015).

   vfisc E. D. Leonel, J. A. Oliveira, F. Saif, J. Phys. A
44, 302001 (2011).

  crank J. Crank. The Mathematics of Diffusion, Oxford University Press (1975).

  gauss1 H. Shore, Communications in Statistics – Theory and Methods 34, 507 (2005).

  gauss2 J. M. Rohrbasser, J. Véron, Population 58, 303 (2003).

  gauss3 K. Pearson, Biometrika 13, 25 (1920).

  robnik2 Črt Lozej and Marko Robnik, Phys. Rev. E 97, 012206 (2018).

  jacobi J. M. Borwein and P. B. Borwein, Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity, Wiley (1987).

  aprox C. M. Bender, S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory, Springer (1999)
  
  harsoula M. Harsoula, K. Karamanos and G. Contopoulos, Phys. Rev. E 99, 032203 (2019).

  prl2004 E. D. Leonel, P. V. E. McClintock and J. K. L. da Silva, Phys. Rev. Lett. 93, 014101 (2004).

  gelfreich K. Shah, V. Gelfreich, V. Rom-Kedar and D. Turaev, Phys. Rev. E 91, 062920 (2015).
 


