Transport in time-dependent random potentials

# Transport in time-dependent random potentials

Yevgeny Krivolapov Physics Department, Technion - Israel Institute of Technology, Haifa 32000, Israel.    Shmuel Fishman Physics Department, Technion - Israel Institute of Technology, Haifa 32000, Israel.
###### Abstract

The classical dynamics in stationary potentials that are random both in space and time is studied. It can be intuitively understood with the help of Chirikov resonances that are central in the theory of Chaos, and explored quantitatively in the framework of the Fokker-Planck equation. In particular, a simple expression for the diffusion coefficient was obtained in terms of the average power density of the potential. The resulting anomalous diffusion in velocity is classified into universality classes. The general theory was applied and numerically tested for specific examples relevant for optics and atom optics.

## I Introduction

Potentials that are random both in space in time are models of noise. Often one seeks for methods of suppressing noise, the present work on the other-hand is focused on interesting properties of dynamics generated by noise. It was subject of many sophisticated studies for nearly 100 years Langevin (1908); Uhlenbeck and Ornstein (1930); Sturrock (1966); Kampen (2007). The response to forces resulting of such potentials typically differs from ordinary diffusion. Specifically, the diffusion coefficients predicted by such mechanisms sensitively depend on the velocity of the particles. Also, if the potential is time-dependent, the energies of the particles will not be constant. The existence of such ‘anomalous diffusion’ has been demonstrated for classical dynamics with spatially and temporally fluctuating potentials Golubović et al. (1991); Rosenbluth (1992); Arvedson et al. (2006); Bezuglyy et al. (2006); Aguer et al. (2009); Bezuglyy et al. (2012). These works claim universal behavior in the sense that for generic random potentials the diffusion coefficient has a universal power-law dependence on velocity , such that as . This in turn implies that asymptotically in time the average velocity satisfies . The average displacement satisfies for one-dimensional systems Golubović et al. (1991); Rosenbluth (1992); Arvedson et al. (2006); Bezuglyy et al. (2006, 2012) (faster than ballistic) and (ballistic transport on average) for systems with dimension higher than one Rosenbluth (1992). In the present paper we will demonstrate that this picture should be extended, and we will introduce such an extension.

The present work is motivated by experiments in optics and in atom optics. The random potentials induced by optical waves in photonic lattices Schwartz et al. (2007) and in atom optics experiments Lye et al. (2005); Sanchez-Palencia et al. (2007), rely on transforming an intensity pattern into an effective potential for the light Efremidis et al. (2002); Fleischer et al. (2003) (the former) or the cold atoms Lye et al. (2005); Sanchez-Palencia et al. (2007) (the latter). Such potentials are naturally described in terms of the Fourier spectrum of the waves inducing them, and their spectral coefficients are assumed to be independent random variables. In the regime where the wave-packet has already acquired large values of velocity, one expects classical mechanics to provide a reasonable approximation. Therefore, we will explore here the classical dynamics for potentials, defined by their spectral content or more precisely their average power spectral density Krivolapov et al. (2012a). In particular, we consider the classical dynamics of a particle in potentials that are random both in space and in time, emphasizing the spreading of the velocity acquired by the particle, as time evolves.

In the present paper we will study dynamics of classical particles in potentials which are defined in terms of their Fourier components. Two types of potentials will be explored:

 V(1)(x,t)=∫^V(k,ω)expi(k⋅x−ωt)dkdω+c.c. (1)

and

 V(2)(x,t)=∣∣V(1)(x,t)∣∣2, (2)

where is a random field chosen such that the distribution of is stationary both in time and space. For time dependent potentials the velocity of the particles may grow, since energy is not conserved. The dominant mechanism of the growth of the velocity is via Chirikov resonances, namely, resonances between the particle dynamics and the external driving. Those resonances occur when the phases in (1) are stationary. For a given instantaneous velocity this happens for

 k⋅v=ω, (3)

and a more complicated expression is found for the potential (2). Indeed, transport takes place only in such regions in phase-space where the density of the resonances is non-zero. The Chirikov resonances provide an intuitive picture for the understanding of transport in phase-space. Since, the potentials are given in terms of their Fourier components it is relatively easy to calculate the spectral content or average power spectral density (PSD), , (see (113)) and with its help also the diffusion coefficient. The spread in velocity is then calculated using the Fokker-Planck (FP) approximation, assuming that the decay of the potential correlations is sufficiently rapid. The random potentials of the form (1) and (2) are natural in optics realizations since these are formed by a superposition of light beams. In recent experiments in optics Schwartz et al. (2007); Levi et al. (2011), light propagates paraxially in a disordered potential, which is produced by utilizing the photo-sensitivity of the medium Efremidis et al. (2002); Fleischer et al. (2003). In the paraxial approximation a probe beam satisfies a Schr dinger like equation where the refractive index plays the role of the potential, the propagation direction plays the role of time and the dispersion relation is approximately . Since the index of refraction is proportional to intensity of the electric field of the writing beam, that varies like of (1), the resulting potential will have the form of (2). Potentials of this type are also relevant for atom optics. There the electric field induced a dipole moment in a neutral atom which interacts with the electric field and results in a potential which is proportional to the intensity of the electric field. The time dependence is induced by a modulation of the intensity of the electric field, such that it is independent of .

Although this work was motivated by experiments in optics, where the dynamics is wave dynamics, the classical dynamics of particles in such potentials is a fundamental question by itself. Therefore, in this work we will study the classical dynamics of particles in potentials of the form (1) and (2). Nevertheless, since it is generally believed that for large velocity (or short wave length) wave dynamics is approximated by classical dynamics, we expect at least qualitatively to describe the original wave problem.

We will classify the various systems into universality classes according to the velocity dependence of the diffusion coefficient in large and small velocity limits. The regimes of the validity of the FP approximation will be studied as well, and the crossover between uniform acceleration and diffusive behavior (in velocity) will be studied. The formalism will be demonstrated for specific examples of the form (1) and (2). Our main result is that although the universality of the form holds always for dimensions higher than one, for one dimensional systems there is a large variety of novel universality classes. A classification scheme for such classes will be presented and it will be demonstrated how these can be realized.

In Sec. II the time required for the particles to ‘discover’ the fluctuations of the potential will be evaluated. After this time we expect the dynamics to be controlled by the FP equation. An expression for the diffusion coefficient in terms of the spectral content or more precisely the average power spectral density (PSD) is derived in Sec. III, it is computed explicitly and its significance is discussed for specific examples. A scaling property of the FP equation is presented in Sec. IV. A detailed calculation of the diffusion coefficient in velocity for dimensions two and higher is presented in Sec. V.1 for potentials of the form of (1), and in Sec. V.2 for potentials of the form (2). In Sec. VI the rich variety of the universality classes found for one-dimensional systems is presented and the asymptotic expansion in powers of is discussed. The results are summarized and discussed in Sec. VII.

## Ii Time domains

In this section we discuss the typical time scale to observe the anomalous diffusion. For this purpose we assume an initial distribution that is narrow in phase-space. Initially the velocity is growing linearly in time (acceleration), since it takes time for the fluctuations of the potential to become effective. First we define the characteristic length and time scales of the potential. Consider the potential given in (1), let be the length over which the potential has an almost constant derivative, resulting in a constant force. It satisfies,

 lx=  ⎷⟨V2⟩−⟨V⟩2⟨|∇V|2⟩. (4)

In a similar way we define the characteristic time scale ,

 lt=   ⎷⟨V2⟩−⟨V⟩2⟨(∂V∂t)2⟩. (5)

Using,

 ⟨V2⟩ = (6) = ∫dω∫dkS(k,ω),

one finds,

 ⟨|∇V|2⟩ = (7) =

where is the average power spectral density (PSD) defined in (16). Therefore the characteristic length scale, is given by,

 l2x=∫dω∫dkS(k,ω)∫dω∫dkk2S(k,ω), (8)

and similarly the characteristic time scale, , is given by,

 l2t=∫dω∫dkS(k,ω)∫dω∫dkω2S(k,ω). (9)

Note, that and are different from the exponential decay rates of correlations. The difference follows from the fact that we are interested in the typical lengths in both time and space where the force is constant, and not in the asymptotic behavior of the correlation function. For times the force resulting from a potential is almost constant, provided that the particle is not displaced more than . For a weak force,

 F0

the particle will travel a distance less than in time , therefore the cross-over out of the uniform acceleration regime takes place at time,

 t∗∼lt. (11)

If the force is strong so that (10) is not satisfied, the particle will reach , at a time shorter than and therefore will be,

 t∗∼√lxF0. (12)

In this work we will consider initial conditions of particles with a narrow distribution both in velocity and position, such that , where is the width of the distribution centered at . For these initial conditions uniform acceleration will take place for time , while for longer time-scales the velocity will exhibit anomalous diffusion. In this work we will assume that the force is weak, such that (10) holds and the crossover to the diffusive regime will take place at (11).

## Iii Anomalous diffusion in velocity

In this section the variation of the velocity of particles by potentials , random both in space and time will be studied on time scales longer than . The equations of motion are,

 dvdt = −∇V(x,t) (13) dxdt = v,

where and are the instantaneous velocity and position. For simplicity we set the mass to be one. We will assume a stationary potential both in space and time and an isotropic distribution of the potential. Therefore, its correlation function is

 ⟨V(x1,t1)V(x2,t2)⟩=C(x1−x2,t1−t2), (14)

and we may take the average of the potentials, that is a constant, to vanish, namely,

 ⟨V(x,t)⟩=0. (15)

The angular brackets denote the average over the ensemble of random potentials. The Weiner-Khinchin theorem implies,

 C(x,t)=∫dω∫dkS(k,ω)expi(k⋅x−ωt), (16)

where is the average power spectral density (PSD) of the potential, that is its spectral content.

The Fokker-Planck (FP) equation for the velocity is given by Kampen (2007),

 ∂P∂t=(∂∂viDij∂∂vj)P, (17)

where is the probability density and is the diffusion tensor, given by the expression,

 (18)

where

 F=−∇V (19)

is the force. Here one assumes that the correlations decay sufficiently fast and the force is sufficiently weak, so that the variation of the velocity on the time-scale of decay of correlations is negligible, and

 x1−x2=v(t1−t2), (20)

this is the usual approximation used to derive the Fokker-Planck equation Kampen (2007). Here and in the rest of the paper summation over repeated indexes is assumed. We will now calculate the force correlation function,

 Kij(x1−x2,t1−t2)=⟨Fi(x1,t1)Fj(x2,t2)⟩, (21)

which can be written using the potential correlation function,

 Kij(x1−x2,t1−t2)=∂2∂x1,i∂x2,j⟨V(x1,t1)V(x2,t2)⟩=−∂2C(x,t)∂xi∂xj|x=x1−x2. (22)

Using (16) one finds,

 Kij(x1−x2,t1−t2)=∫dω∫dkkikjS(k,ω)expi[k⋅(x1−x2)−ω(t1−t2)], (23)

and the diffusion tensor is,

 Dij(v)=12∫∞−∞Kij(vτ,τ)dτ. (24)

Which can be written as,

 Dij(v) = 12∫∞−∞dτ∫dω∫dkkikjS(k,ω)expiτ(k⋅v−ω). (25)

Taking the integral over first, gives a delta function

 Dij(v) = (26)

Then the integral over can be easily performed to give,

 Dij(v)=π∫dkkikjS(k,k⋅v). (27)

Since is integrable by definition, we will define a function , so that

 S(k,ω)=∂2∂ω2F(k,ω). (28)

Then, the diffusion tensor can be written as,

 Dij(v)=∂2∂vi∂vjI(v), (29)

where

 I(v)=π∫dkF(k,k⋅v), (30)

and we have used the fact that the integral is isotropic. Equation (29) can be decomposed to,

 Dij(v) = ∂2∂vi∂vjI(v)=∂∂vi(vjv∂I∂v) = (δijv−vivj/vv2)∂I∂v+vivjv2∂2I∂v2.

For a general isotropic diffusion tensor of the form,

 (32)

the operator on the RHS of (17) is,

 ^L0=∂∂viDij(v)∂∂vj=v−(d−1)∂∂vv(d−1)fP(v)∂∂v. (33)

We turn now to justify this result. It is convenient to use hyper-spherical variables, such that,

 ∂∂vi=∂v∂vi∂∂v+d−1∑j=1∂ϕj∂vi∂∂ϕj=viv∂∂v+d−1∑j=1∂ϕj∂vi∂∂ϕj, (34)

where is the dimension of space and are the angle coordinates. Using (32) and the new variables the operator (33) reduces to,

 ^L0 = ∂∂vi[(δij−vivjv2)fT(v)+vivjv2fP(v)]∂∂vj (35) = ∂∂vi[(δij−vivjv2)fT(v)+vivjv2fP(v)]vjv∂∂v = ∂∂vi[(vi−vi)1vfT(v)+vivfP(v)]∂∂v = ∂∂vivivfP(v)∂∂v.

Which using the chain rule can be written as,

 ^L0 = (δiiv−vivi/vv2)fP(v)∂∂v+∂∂vfP(v)∂∂v (36) = 1v(d−1)fP(v)∂∂v+∂∂vfP(v)∂∂v = v−(d−1)∂∂vv(d−1)fP(v)∂∂v.

Therefore, combining (III) and (33) one identifies . Therefore, the operator is given by,

 (37)

The diffusion coefficient is found to be,

 D(v)=∂2I∂v2. (38)

It can be simplified using the definition (30) and (28),

 D(v)=π∫dk(k⋅^v)2S(k,k⋅v). (39)

Using (39) and aligning the component of with velocity we get,

 D(v)=Sd2∫2π0dθ∫∞0dkkd+1cos2θS(k,kvcosθ), (40)

where is the surface of the dimensional hyper-sphere. Changing variables to, we get

 D(v)=Sdv3∫v−vdyy2√1−(y/v)2∫∞0dkkd+1S(k,ky), (41)

where the factor of two was eliminated due the multiplicity. This expression assumes implicitly that . For one dimensional systems there is no integration over the angle (see Eq. (96) in Sec. VI). The different dependence on the velocity of (41) and (96) is the origin of the richness of the one dimensional behavior. For large velocities (41) has the following asymptotic behavior,

 D(v)∼D3v3, (42)

with

 D3=2Sd∫∞0dy∫∞0dky2kd+1S(k,ky). (43)

Therefore, for dimensions higher than one the asymptotic behavior of the diffusion coefficient is indeed universal for any choice of the PSD and is given by (42). The diffusion coefficient for zero velocity using (40) is given by,

 D(0)=Sd4∫∞0dkkd+1S(k,0). (44)

If , as is the case for , the Fokker-Planck equation for the velocity (17) is asymptotically equivalent to the equation,

 ∂P∂t=(v−(d−1)∂∂vv(d−1)D3v3∂∂v)P, (45)

which has the scaling solution

 P(v,t)=1td/5g(v5t). (46)

The resulting growth of the mean kinetic energy is,

 12⟨v2⟩∼t2/5. (47)

This behavior is considered in the literature as universal Golubović et al. (1991); Arvedson et al. (2006); Bezuglyy et al. (2006, 2012). For dimensions two and higher this is indeed the case, however for one dimensional systems other behaviors are possible, as will be demonstrated in Sec. VI.

## Iv Scaling properties

Let us assume that the average power spectral density, , has natural spatial and temporal frequency scales, and , such that,

 S(k,ω)=V20kd0ω0~S(kk0,ωω0), (48)

where is a constant which determines the strength of the potential. Then (39) takes the form,

 D(v)=πV20kd0ω0∫(k⋅^v)2~S(kk0,k⋅vω0)dk, (49)

where is a dimensionless PSD. Rescaling the variables

 k=k′k0v=v′v0v0≡ω0k0, (50)

gives

 D(v′)=πV20k20ω0∫(k′⋅^v)2~S(k′,k′⋅v′)dk′. (51)

The Fokker-Planck equation for the velocity (17) is invariant under the transformation of variables,

 v → v′ t → t′(πI20k20ω0)−1. (52)

## V Special potentials

In this section the spreading in phase-space is studied for specific potentials. Potentials of the form (1) will be studied in subsection V.1 while potentials of the form (2) will be studied in subsection V.2. In this section the emphasis will be on the behavior for dimensions higher than one, while one dimensional systems will be analyzed in the next section.

### v.1 Potentials which are a superposition of waves

Natural potentials to consider are potentials which are composed of a superposition of standing waves (c.f., Eq. (1)),

 V(x,t)=∫^V(k)expi(k⋅x−ω(k)t)dk+c.c. (53)

with a dispersion relation . We will assume that the amplitudes and the wave numbers are independent random variables. Leading to,

 ⟨^V(k1)^V∗(k2)⟩=V20f(k1)δ(k1−k2), (54)

where and is the distribution of the wave numbers. The correlation function of the potential is given by,

 C(x1−x2,t1−t2)=V20∫f(k)expi(k⋅(x1−x2)−ω(k)(t1−t2))dk+c.c. (55)

Using (16) the PSD is,

 S(q,ω) = 1(2π)d+1∫C(x,t)exp−i(q⋅x−ωt)dxdt = +

By making the integrals over and we get,

 S(q,ω) = V20f(q)δ(ω−ω(q))+V20f(−q)δ(ω+ω(q)). (57)

Starting from (41) and assuming that both and are isotropic gives,

 D(v)=SdV20v3∫v−vdyy2√1−(y/v)2∫∞0dkkd+1f(k)[δ(ky−ω(k))+δ(ky+ω(k))]. (58)

We can now perform the integral over and obtain,

 D(v)=2Sdv3V20∫∞0dkkd−2f(k)ω2(k)√1−ω2(k)/(v2k2)H(v−ω(k)k), (59)

where is a Heaviside step function, which follows from the restriction of the delta function. In the limit of large velocities we have,

 D(v)∼D3v3, (60)

with

 D3=2SdV20∫∞0dkkd−2f(k)ω2(k). (61)

For zero velocity using (44) we have,

 D(0) = Sd2V20∫∞0dkkd+1f(k)δ(ω(k)). = Sd2V20∫∞0dkkd+1f(k)|ω′(0)|δ(k)=0,

where we have assumed that

 ω(k)=0, (63)

has only one solution , which is a general property of physically relevant dispersion relations.

We will proceed by calculating the diffusion coefficient for a particular choice of the dispersion relation, , and a uniform density of wave numbers,

 f(k)={1VdkdR0kR, (64)

where is the volume of the dimensional unit hyper sphere. Therefore, we have

 D(v)=dV202kdRv3⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∫2v0dkkd+2√1−k2/(4v2)vkR/2, (65)

where we have used the fact . Taking the integrals gives,

 D(v)=dV202kdRv3⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩2d+2√πΓ(d+32)Γ(2+d2)vd+3vkR/2, (66)

where is Gauss’s hypergeometric function. For the special case of ,

 D(v)=V20⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩6π(vkR)2vkR/2. (67)

The prediction of the FP equation was compared to the Monte-Carlo simulation and is presented in Fig. 1. It is found that uniform acceleration takes place for times and the FP prediction clearly fail in this regime. This can be expected since from (67) , the assumption that the velocity is constant (20) cannot hold, rendering the FP approximation inconsistent. Nevertheless, for the FP predictions are satisfied. Note that, for any potential discussed in this section the scaling properties of Sec. IV hold.

Numerically, it is found that for large velocities () the average square position is growing ballisticly, , as expected from previous studies Rosenbluth (1992); Golubović et al. (1991); Bezuglyy et al. (2006, 2012), while for small velocities (), , as expected for uniform acceleration.

### v.2 Potentials proportional to the intensity of a superposition of waves

In this subsection the diffusion coefficient for potentials of the form (2) will be calculated. These potentials appear in experiments with neutral atoms and in some experiments in optics, where the refractive index (which plays the role of the potential) of a photosensitive material is proportional to the intensity of light Efremidis et al. (2002); Fleischer et al. (2003). The complex field, , which represents the superposition pattern of waves is given by,

 U(x,t)=∫dk^U(k)expi(k⋅x−ω(k)t), (68)

where are random Fourier coefficients, and is some dispersion relation. Following the analysis of the previous subsection we will assume that is composed of plane-waves with independent phases and amplitudes, leading to

 ⟨^U(k)⟩ = 0 (69) ⟨^U(k1)^U∗(k2)⟩ =I0 f(k1)δ(k1−k2) ⟨^U(k1)^U(k2)⟩ = 0,

where and is the probability density of the wave numbers . The potential is proportional to the intensity of the ,

 V(x,t)=|U(x,t)|2=∫dk1dk2^U(k1)^U∗(k2)expi[(k1−k2)⋅x−(ω(k1)−ω(k2))t], (70)

Using the assumptions (69) we get that the potential is constant on average,

 ⟨V(x,t)⟩ = (71) = I0∫dkf(k)=I0.

The correlation function of the potential is,

 ⟨V(x1,t1)V(x2,t2)⟩ = ∫dk1dk2dk3dk4⟨^U(k1)^U∗(k2)^U(k3)^U∗(k4)⟩ (72) × expi[(k1−k2)⋅x1−(ω(k1)−ω(k2))t1] × expi[(k3−k4)⋅x2−(ω(k3)−ω(k4))t2].

Since we have assumed that the complex amplitudes of the plane-waves are independent random variables we can decompose,

 ⟨^U(k1)^U∗(k2)^U(k3)^U∗(k4)⟩ = ⟨^U(k1)^U∗(k2)⟩⟨^U(k3)^U∗(k4)⟩ (73) + ⟨^U(k1)^U∗(k4)⟩⟨^U∗(k2)^U(k3)⟩ = +

Therefore, the correlation function reduces to,

 C(x1−x2,t1−t2) = I20 + I20∫dk1dk2f(k1)f(k2) × expi[(k1−k2)⋅(x1−x2)−(ω(k1)−ω(k2))(t1−t2)].

The PSD is just,

 S(q,ω) = 1(2π)d+1∫dx∫dtC(x,t)exp−i(q⋅x−ωt) = I20δ(q)δ(ω) + I20(2π)d+1∫dx∫dt∫dk1dk2f(k1)f(k2) × expi[(k1−k2−q)⋅x−[ω(k1)−ω(k2)−ω]t]. (76)

Taking the integral with respect to and we get,

 (77)

or

 S(q,ω) = I20δ(q)δ(ω)+I20∫dkf(k)f(k−q)δ[ω−(ω(k)−ω(|k−q|))]. (78)

The diffusion coefficient is given by (see, (39))

 D(v)=πI20∫dqdk(q⋅^v)2f(k)f(k−q)δ(q⋅v−(ω(k)−ω(|k−q|))). (79)

Changing the variables to,

 q′ = k−q k′ = k, (80)

and suppressing the primes for simplicity gives,

 (81)

Choosing the - components of and such that they are aligned with the velocity yields,

 D(v) = πI20(Sd2π)2∫∞0dq∫2π0dθq∫∞0dk∫2π0dθk(qk)d−1 × (kcosθk−qcosθq)2f(k)f(q)δ(v(kcosθk−qcosθq)−(ω(k)−ω(q))).

Changing variables to, and in (V.2) we have the Jacobian,

 J=1v2√(1−(yk/v)2)(1−(yq/v)2). (83)

In these variables the diffusion coefficient is,

 D(v) = 4πI20(Sd2π)21v4∫∞0dq∫∞0dk∫v−vdyk∫v−vdyq(qk)d−1 ×

For large velocities this gives,

 D(v)∼D3v3, (85)

with

 D3=4πI20S2d∫∞0dqqd−1∫∞0dkkd−2(ω(k)−ω(q))2f(k)f(q). (86)

Since is non-negative the integral for does not vanish and for an appropriate choice of it is convergent. The diffusion coefficient for zero velocity can be calculated from (V.2) by setting ,

 D(0) = πI20(Sd2π)2∫∞0dq∫2π0dθq∫∞0dk∫2π0dθk(qk)d−1 × (kcosθk−qcosθq)2f(k)f(q)δ(ω(k)−ω(q)).

Assuming that is a monotonic function we have,

 δ(ω(k)−ω(q))=1|ω′(q)|δ(k−q), (88)

and therefore,

 D(0) = πI20(Sd2π)2∫2π0dθq∫2π0dθk(cosθk−cosθq)2 (89) × ∫∞0dqq2df2(q)|ω′(q)|.

After integrating over the angles,

 D(0)=πI20S2d∫∞0dqq2df2(q)|ω′(q)|, (90)

that is non-vanishing. Therefore, the potentials described in this section are generic, namely, produce the universal behavior both for small and large velocities. We will now demonstrate the calculation for a particular choice of dispersion relation, and a uniform distribution function on a disk given by (64) the integrals in (V.2) cannot be taken explicitly, however we have calculated them numerically (see Fig. 5) and the resulting dynamics is demonstrated in Fig. 2.