Self-consistent chaotic transport in a high-dimensional mean-field Hamiltonian map model D. Mart´ınez-del-R´ıo(cid:5),∗ D. del-Castillo-Negrete(cid:63)†, A. Olvera(cid:5),‡ R. Calleja(cid:5)§ 6 1 0 (cid:5) IIMAS-UNAM, Mexico, D.F. 04510 2 (cid:63) Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831-8071 n a J 5 Abstract Self-consistent chaotic transport is studied in a Hamiltonian mean-field model. The ] S model provides a simplified description of transport in marginally stable systems in- D cluding vorticity mixing in strong shear flows and electron dynamics in plasmas. Self- . h consistencyisincorporatedthroughamean-fieldthatcouplesallthedegrees-of-freedom. t The model is formulated as a large set of N coupled standard-like area-preserving twist a m maps in which the amplitude and phase of the perturbation, rather than being con- [ stant like in the standard map, are dynamical variables. Of particular interest is the study of the impact of periodic orbits on the chaotic transport and coherent structures. 1 v Numerical simulations show that self-consistency leads to the formation of a coherent 2 macro-particle trapped around the elliptic fixed point of the system that appears to- 4 getherwithanasymptoticperiodicbehaviorofthemeanfield. Tomodelthisasymptotic 9 0 state, we introduced a non-autonomous map that allows a detailed study of the onset 0 of global transport. A turnstile-type transport mechanism that allows transport across . 1 instantaneousKAMinvariantcirclesinnon-autonomoussystemsisdiscussed. Asafirst 0 step to understand transport, we study a special type of orbits referred to as sequential 6 1 periodic orbits. Using symmetry properties we show that, through replication, high- : v dimensional sequential periodic orbits can be generated starting from low-dimensional i periodic orbits. We show that sequential periodic orbits in the self-consistent map can X be continued from trivial (uncoupled) periodic orbits of standard-like maps using nu- r a merical and asymptotic methods. Normal forms are used to describe these orbits and to find the values of the map parameters that guarantee their existence. Numerical simulations are used to verify the prediction from the asymptotic methods. ∗[email protected] †[email protected] ‡[email protected] §[email protected] 1 1 Introduction There are many problems in physics and the applied sciences where transport plays an important role. Some examples include the dispersion of pollutants in the atmosphere and the oceans, the efficiency of mixing in chemical reactions, the magnetic confinement of fusion plasmas, and vorticity mixing in fluids among others. The study of transport can be classified in two groups: passive and active transport. In the case of passive transport, the elements which are transported do not modify the velocity field of the flow while in the active case, the advected field determines the velocity field through a self-consistent feedback mechanism. In the present paper, we focus on active transport. The study is based on a simplified model of self-consistent transport in marginally stable systems including vorticity mixing in strong shear flows and electron dynamics in plasmas. The model, originally presented in Ref. [5] and later on studied in Ref. [3], consists of a large set of N mean-field coupled standard-like area-preserving twist maps in which the amplitude and phase of the pertur- bation (rather than being constant) are dynamical variables. An important part of our methodology is based on the use of normal form methods to describe the evolution of the dynamical variables of the system. The rest of this paper is organized as follows. Section 2 presents a brief summary of the model and some results originally contained in Ref. [5] that will be used in the present work. In Sec. 3, we study the asymptotic behavior of numerical simulations of the self-consistent map model and in particular the formation of coherent structures. To study these structures, we introduce a reduced transport map model that mimics the asymptotic behavior and that can be interpreted as an uncoupled system in an external oscillatory field. Large iterations of this map are computed for particular sets of initial conditions to show that global diffusion can occur as a result of the oscillation of the perturbation parameter. In Sec. 4, we study the properties of periodic orbits in the self–consistent map and introduce an asymptotic procedure to find a particular kind of periodic orbits that we name sequential periodic orbits. In Sec. 5 we present a simple example illustrating how these sequential periodic orbits can be computed using normal forms and present the relationships between the parameters of the map that should be satisfied for these orbits to exits. The values of the parameters of the sequential periodic orbits obtained by numerical continuation methods are compared with the normal forms results. It is also shown that these numerical and asymptotic results agree up to machine precision. Section 6 present the conclusions. 2 Hamiltonian mean field model of self-consistent transport The advection-diffusion equation ∂ ζ +V·∇ζ = D∇2ζ, (1) t 2 where ζ is an advected scalar, V is the velocity field, and D is the diffusivity is one of the fundamental models to study transport in the case of incompressible (∇ · V = 0) velocity fields. Restricting attention to 2-dimensions, and introducing the stream-function ψ, V = zˆ×∇ψ, (2) Eq. (1) can be written as ∂ ζ −(∂ ψ)∂ ζ +(∂ ψ)∂ ζ = D∇2ζ. (3) t y x x y Broadly speaking, there are two different classes of transport problems. In the case of pas- sive transport, ζ, evolves without affecting the velocity field, V. A typical example is the transport of low concentration pollutants in the atmosphere and the oceans. In the case of active transport, ζ determines V through a dynamical self-consistency constraint of the form F(ζ,v) = 0 involving an integral and/or differential operator. A paradigmatic example is the 2-dimensional Navier-Stokes equation for an incompressible fluid. In this case, ζ in Eqs. (1) and (3) denotes the fluid vorticity that is self-consistently coupled to the velocity field according to the constrain ζ = (∇×V)·zˆ = ∇2ψ. (4) Active transport problems, are intrinsically nonlinear and as a result, harder to study than passive transport problems. For example, the well-known challenges in understanding fluid turbulence reside in the nonlinearity in Eqs. (3)-(4). It is thus of significant interest to develop simplified models that capture the basic elements of the self-consistent coupling within a relatively simple mathematical setting. One of these simplified descriptions is the Single Wave Model (SWM) originally proposed in plasma physics [15, 4] and used to study self-consistent transport in marginally stable fluids in the presence of strong shear flows [5]. In general, the stream-function ψ has a complicated spatio-temporal dependence. However, in the SWM, ψ is assumed to have the relatively simple form 1 ψ = − y2 +a(t)eix +a∗(t)e−ix, (5) 2 and the vorticity stream-function self-consistent coupling in Eq. (4) reduces to (cid:90) (cid:90) da i −iUa = e−ixdx ζ(x,y,t)dy, (6) dt 2π where a = a(t) is in general complex, and U is a constant real parameter. That is, whereas in the Navier-Stokes equation the stream-function has a general spatio-temporal dependence 3 that is determined at each instant of time by solving the Poisson equation in (4), in the SWM the spatial dependence of the stream-function is given and self-consistency only enters when one determines the temporal dependence of ψ though the amplitude a(t) according to the ordinary differential equation in (6). According to Eq. (2), in the SWM the velocity field consists of a linear component in the x-direction and a traveling-wave component in the y-direction, i.e. V = yi+Re[ia(t)eix]j, where Re denotes the real part. Despite its relative mathematical simplicity, the SWM is able to capture important dy- namics of the full Navier-Stokes equation. In particular, as Fig. 1 from Ref. [5] shows, the SWM exhibit the standard Kelvin-Helmholtz instability leading to the formation of the fa- miliar cats’ eyes vorticity structure found in unstable shear flows. The result in Fig. 1 was obtained from the direct numerical integration of Eqs. (1) and (6) with D = 0.001, U = −1, and initial condition: ζ(x,y,t = 0) = e−y2/2[1−0.2y cos(x)]. (7) Although for the sake of simplicity we have stressed the fluid-dynamics interpretation of the SWM, it is important to reiterate that as discussed in Refs. [4, 5, 6] the SWM has its origins in plasma physics and it can be equally applied to the study of electron dynamics in a Vlasov-Poisson plasma in a neutralizing ion background. In this case, ζ corresponds to the single-particle electron distribution function, ψ is the electrostatic potential, the vorticity equation becomes the Vlasov equation, and the (x,y) variables correspond to the (x,u)-phase space variables. To reformulate the SWM as a finite (but arbitrarily large) degrees-of-freedom Hamil- tonian dynamical system, we assume from now on D = 0 and introduce the point-vortex representation N (cid:88) ζ(x,y,t) = 2π Γ δ[x−x (t)]δ[y −y (t)], (8) j j j j=1 where(x (t),y (t))denotetheLagrangiantrajectoriesofthek = 1,2,...N pointvorticeswith k k intensities Γ . In this case, it can be shown that Eqs. (3), (5) and (6) imply the following set k of Hamiltonian differential equations determining (x (t),y (t)) k k dx ∂H dy ∂H k k = , = − , (9) dt ∂y dt ∂x k k where the Hamiltonian is N (cid:20) (cid:21) (cid:88) 1 H = y2 −a(t)eixj −a∗(t)e−ixj , (10) 2 j j=1 4 and the function a = a(t) is determined from N da i (cid:88) −iUa = Γ e−ixj , (11) j dt N j=1 which is the point-vortex representation of the SWM self-consistent vorticity stream-function relation in Eq. (6). The Hamiltonian model in Eq. (9) is a mean field model in the sense that the dynamics of the particles (x ,y ) is determined by the time dependent field a(t) which k k according to Eq. (11) depends on the mean properties of the particles’ positions. Defining √ a = J e−iθ, p = Γ y , (12) k k k equations (9) and (11) can be written as a full N +1-degrees of freedom Hamiltonian system dx ∂H dp ∂H k k = , = − , (13) dt ∂p dt ∂x k k dθ ∂H dJ ∂H = , = − . (14) dt ∂J dt ∂θ where, (cid:88)N (cid:20) 1 √ (cid:21) H = p2 −2Γ J cos(x −θ) −UJ. (15) 2Γ j j j j j=1 The last step to construct the map model is to perform a time discretization of Eqs. (13) and (14) to obtain the self-consistent standard map model originally proposed in Ref. [5], xn+1 = xn +yn+1 (16a) k k k yn+1 = yn −κn+1 sin(xn −θn) (16b) k k k (cid:112) κn+1 = (κn)2 +(ηn)2 +ηn (16c) 1 ∂ηn θn+1 = θn −Ω+ (16d) κn+1 ∂θn where k = 1,2,...,N, xn,θ ∈ [0,2π), yn,κn ∈ R, Ω = Uτ and the ηn is defined as: k n k N (cid:88) ηn := γ sin(xn −θn), (17) k k k=1 where τ is the discretization time step. The map in Eqs. (16) is a 2N +2 dimensional map with N+1 parameters. The N parameters, γ , represent the intensities of the point vortices, k and the constant parameter Ω is related to the parameter U in the SWM. 5 Figure 1: Cat’s eye formation and vorticity mixing in the single wave model, Eqs. (3), (5) and (6), with initial conditions in Eq. (7). The gray scale denotes the vorticity values with white corresponding to ζ = 1 and dark gray corresponding to ζ =0 [After Ref. [5]]. Note that Eqs. (16a)-(16b) have the structure of the well-known symplectic standard map xn+1 = xn +yn+1 (18a) yn+1 = yn −ε sin(xn −φ) (18b) with constant, fixed perturbation ε, and phase φ (which is usually taken equal to zero). However, in the self-consistent map, the perturbation parameter, κn, as well as the phase, θn, depend on the iteration, n, and their dynamics is dictated by Eqs. (16c)-(16d) which are the discrete map version of the self-consistent SWM coupling in Eq. (11). Although Eqs. (16c)- (16d) do not define a symplectic transformation, they can be rewritten as a symplectic (although implicit) map in the mean field degrees of freedom [5]. The analogy with the standard map allows us to interpret the map in Eqs. (16a)-(16b) as N coupled oscillators through their phase and amplitudes by the mean field in Eqs. (16c)-(16d). We will take advantage of these similarities not only to identify the equations of the map, but also in its perturbation analysis. Note that in the limit γ → 0 in Eq. (17), the oscillators in (16) k depending on θ0 and Ω decoupled and their equations simply correspond to N copies of the standard map. 6 3 Coherent structures and transport The evolution of the self-consistent map in Eqs. (16) has been studied for different initial conditions. Figures 2 and 3 show the results of a simulation of N = 13,440 coupled maps with initial conditions {(x0,y0)} uniformly distributed on the region [0,2π] × [−0.3,0.3] in k k the (x,y) plane and γ = 3×10−6 for k = 1,...N. The initial condition of the mean-field k was κ0 = 10−4 and θ0 = 0, and we assumed Ω = 0. The simulations show that the self-consistent map reproduces the coherent structures observed in the single wave model (Figure 1). In particular, while a group of particles exhibit coherent behavior trapped in the center of the cat’s eye, those particles located in the separatrix exhibit a strong dispersion. Figure 2: Timeevolutionoftheself-consistentmapinEq.(16)withN =13440initialconditionsuniformly distributed in [0,2π]×[−0.3,0.3] with κ0 = 10−4, γ = 3×10−6, θ0 = 0 and Ω = 0. The frames show the k instantaneouscoordinatesoftheN initialconditions,atn=2,6,12,20and66intheregion[0,2π]×[−0.8,0.8] of the (x,y) plane. The colors label the y-coordinate of the initial condition, with red denoting y0 close to k y =0 and blue further away. Theevolutionofthemeanfieldisrepresentedbythevariables(κn,θn), whichareshownin Fig. 3. The behavior of κn starts with a fast growth until it achieves a maximum value, after that, the κn oscillates around a mean value κ¯ and the amplitude of oscillation is bounded by ∆κ. A similar situation is observed with the behavior of ϑn+1 = θn+1−θn. Different types of dynamics, including cases in which the mean field decays to zero or saturates at a constant fixed value can be found in Refs. [1, 3] for similar self-consistent maps. Note that in Fig. 3, κn does not reach the critical value κ = 0.971635406, which is the value of the parameter κ c which corresponds to the destruction of the invariant circle with rotation number equal to the golden mean1 of the standard map [9, 12]. For any κ < κ there is no global diffusion in the c standard map because of the existence of invariant circles, that give rise to barriers in phase space. The existence or non-existence of global diffusion in the self-consistent map depends in a nontrivial way on the dynamics of κn. On a more fundamental level, the observed rapid 1The last invariant circle not homotopic to a point. It must be noted that due the choice of scale of the √ map, the rotation number is ω =2πγ, instead of the golden mean: γ = 5−1. 2 7 growth of κn for a given initial condition is closely related to the linear stability properties of the corresponding initial condition in the single wave model. In particular, Ref. [4] presents the necessary and sufficient conditions for the linear stability (i.e., exponential growth of the mean field amplitude) of a given initial condition in the context of the continuous, N → ∞ limit. These ideas might help understand the conditions for the growth of κn. However, one must be careful before drawing conclusions as the self-consistent map model discussed here is obtained by simplifying drastically the single wave model by approximating the continuous limit when N → ∞. (a) (b) Figure 3: Time evolution of the mean-field variables in the self-consistent map in Eq. (16) with initial conditions taken in a uniform grid in [0,2π]×[−0.3,0.3] with N =13440, κ0 =10−4, γ =3×10−6, θ0 =0 k and Ω=0. κn is shown in (a) and ϑn+1 =θn+1−θn in (b). Perhaps contrary to the intuition, it is observed that global diffusion exists even when κ¯ < κ . It is also worth mentioning that when the instantaneous coordinates (xn,yn) of each c k k degree-of-freedom are plotted on the same plane like in Fig. 2, the amplitude and shape of the cat’s eye structure is in good agreement with the invariant manifolds emanating from the hyperbolic fixed point of the standard map calculated with a perturbation parameter equal to the instantaneous value of κn+1. This gives rise to the following question: What is the mechanism that allows the diffusion across the invariant curves on the self-consistent map? In Figure 4 we observe the formation, for relatively small times, of the macro particle coherent structure trapped around the elliptic fixed point, and at the same time we have the formation of the heteroclinic tangle responsible for the high mixing region around the separatrix of the cat’s eye. 8 Figure 4: Plot of the projected phase space coordinates (x ,y ) on the (x,y)-plane of a simulation of i i the self-consistent map in Eq. (16) in the oscillatory κ regimen. Also shown in black is the heteroclinic tangle generated by the unstable invariant manifold of the hyperbolic fixed point of the standard map with perturbation parameter equal to κn. The initial conditions are the same as those in Fig. 2, except that to enhance the heteroclinic tangle the higher values κ0 =0.005 and γ =0.0005 were used. To address this question, and motivated by the observed oscillation of κn, we propose the following non-autonomous map consisting of two copies of the standard map applied sequentially with alternating values of κ,  yn+1 = yn −κnsinxn κ if n is odd  1 where κn = (19) xn+1 = xn +yn+1  κ if n is even 2 Without loss of generality, we will assume that κ is greater than κ . 2 1 To explore the onset of global diffusion for values of κ and κ less than κ , we considered 1 2 c a set of N initial conditions uniformly distributed in the region [0,2π] × [0,y ] and the max map is iterated n times, with n less than some given maximal value M. We then found the number of initial conditions that reached the semi-space y > y + 2π. In turn, this max means that one orbit could pass through the invariant circles which exist in the standard map with κ and κ less than κ . Figure 5 shows the percentage of initial orbits that reached 1 2 c 9 Figure 5: Onsetofglobaltransportinthenon-autonomousmap(19)fordifferentvaluesofκ¯ =(κ +κ )/2. 1 2 The x-axis corresponds to the difference ∆κ=κ −κ for a fixed value of κ¯. The y-axis corresponds 2 1 to the percentage of initial conditions, (x0,y0), taken from a 100×100 regular partition of the rectangle k k [0,2π]×[0,π/5], that escape due to global transport in the sense that after M = 100,000 iterations satisfy yM ≥π/5+2π. k the semi-space y > y +2π for a given κ and κ . For convenience we show the results as max 1 2 function of ∆κ = κ −κ and κ¯ = κ1+κ2. Note that in most cases κ¯ is lower than κ and yet 2 1 2 c there are cases of global diffusion even when κ < κ < κ . This results provide numerical 1 2 c evidence of the existence of global diffusion. However, as shown in Fig. 5 there are also cases with no diffusion, a result consistent with Moser’s theorem [14] that ensures that for values of κ and κ sufficiently small and for which the twist condition is satisfied, the map (19) 1 2 should have invariant circles that forbid global transport. According to Birkhoff’s theorem, an invariant circle is the graph of a function y = f(x) and the position of this curve in the (x,y) plane is determined by the value of κ [12]. As shown in Fig. 6, the shape and position of invariant curves with a given rotation number for κ (cid:54)= κ in general do not coincide, and the set of points located between these two 1 2 curves can move upward and downward when we iterate the map (19). Like in a revolving 10