Violation of adiabaticity in magnetic billiards due to separatrix crossings

We consider dynamics of magnetic billiards with curved boundaries and strong inhomogeneous magnetic ﬁeld. We investigate a violation of adiabaticity of charged particle motion in this system. The destruction of the adiabatic invariance is due to the change of type of the particle trajectory: particles can drift along the boundary reﬂecting from it or rotate around the magnetic ﬁeld at some distance from the boundary without collisions with it. Trajectories of these two types are demarcated in the phase space by a separatrix. Crossings of the separatrix result in jumps of the adiabatic invariant. We derive an asymptotic formula for such a jump and demonstrate that an accumulation of these jumps leads to the destruction of the adiabatic invariance. V C 2015 Publishing LLC The magnetic billiard is a well known model in statistical physics and nonlinear dynamics. This model imitates many real complicated systems describing various problems in magnetism, plasma physics, solid state physics, and astrophysics. The destruction of adiabaticity of charged particle motion in magnetic billiards simulates such important phenomena as exponential decay of correlations and breaking of time reversal symmetry in dynamical systems. In this paper, we investigate a new mechanism responsible for the violation of adiabaticity for systems with impacts. This mechanism is related to the separatrix crossing in the phase space. We apply methods of the adiabatic theory of smooth systems to magnetic billiards (i.e., for systems with impacts) and derive the expression determining the rate of adiabatic invariance destruction.


I. INTRODUCTION
Billiards are considered as universal models of many physical processes (see special section devoted to billiard models 27 and review Ref. 19). Although, these are rather simple systems with geometrical laws of particle motion, they are able to simulate various statistical properties, 8,41 statistical relaxation processes, 10 the property of exponential decay of correlations, 12,13 ergodicity of gas models, 6 and breaking of time reversal symmetry of motion. 11 Billiard also serves as a simple model of the Fermi acceleration process. 28,29,36 Among many billiard models, the model of magnetic billiards (i.e., motion of charged particles inside some plane domain with the imposed magnetic field 39 ) plays an important role. The magnetic billiard serves as a useful model for quantum mechanical problems 15,22,30 and modeling of magnetic properties of materials (diamagnetism, paramagnetism, conductance, magnetic edge state, etc.), 1,20,26,31,42 interaction of matter with laser field, 21,38 and molecular/atom dynamics. 2,3,37 In the model of magnetic billiards, the magnetic field amplitude controls a regularity level of the billiard motion. 7,16 There are several investigations of stochastization of motion in the magnetic billiards with and without curvature of the boundaries of domain where charged particles move. 18,23,43,44 One of possible triggers of stochastization of motion is the destruction of adiabatic invariance (for magnetic billiards, the role of this invariant is played by the classical magnetic moment 35 ). It was shown that destruction of the adiabaticity in magnetic billiards can realize in case of an inhomogeneous magnetic field. 34 The primary origin of this process is related to jumps of the adiabatic invariant due to the change of type of particle motion (from drifting along the boundary to gyrorotation at some distance from the boundary without collisions). Different types of particle motion correspond to different domains in the phase space demarcated by a separatrix. There is a general theory of destruction of adiabatic invariance due to the separatrix crossing in smooth systems. 9,32,33,46 Recently, this theory was generalized for systems with elastic impacts. 17 In this paper, we adopt this theory to derive the equation for the rate of the adiabatic invariance destruction in magnetic billiards with an inhomogeneous magnetic field.

II. GENERAL EQUATIONS
We apply the Hamiltonian approach for investigation of charged particle motion in system with impacts (see also Refs. 49 and 50). This motion inside the 2D domain near the curved boundary and under the action of the transverse magnetic field can be described using curved coordinates s and r: s is directed along the boundary, while r is the transverse coordinate (see scheme in Fig. 1).
The corresponding Hamiltonian and sympletic structure have the form 34 Here, p r and p s are components of particle momentum, e and c are particle charge and speed of light, e À1 B is the magnetic field magnitude, k(s) is the boundary curvature, and e is a constant positive parameter. We consider the charged particle motion in the strong magnetic field, i.e., e ( 1 and amplitude of B is of the order of one. The reflection from the boundary r ¼ 0 corresponds to the transformation p s ! p s , p r ! Àp r . We introduce a canonical momentum P s ¼ p s þ e À1 A(r, s) where Thus, Hamiltonian (1) takes the form We expand Eq. (2) up to terms $r 2 where B 0 ¼ B(0, s) and B 0 0 ¼ @B=@rj r¼0 . We introduce new variables r/e ! r, s/e ! s, t/e ! t and expand Hamiltonian (3) as H ¼ H 0 p r ; r; es; P s ð Þþ eH 1 r; es; P s ð Þ; We introduce dimensionless variables and parameters ðp r ; P s Þ ! ðp r ; P s Þ= ffiffiffiffiffiffiffiffi mh 0 p ; and rewrite Eq. (5) as Here, h 0 is a characteristic value of particle energy, while dimensionless (normalized on h 0 ) energy is h. For frozen (s, P s ) Hamiltonian H 0 describes periodic particle oscillations near the boundary r ¼ 0. The phase portrait shown in Fig. 1 demonstrates that there are two types of particle trajectories: with collisions and without collisions with the boundary. The phase region filled by trajectories of these two types are demarcated by the separatrix (shown by the solid grey curve). The region surrounded by the separatrix corresponds to particles that do not collide with the boundary. Particles colliding with the boundary r ¼ 0 are reflected from it. For slowly varying (s, P s ), the position of the separatrix changes due to magnetic field inhomogeneity, and particles can change the type of motion, i.e., particles colliding with the boundary can start moving along the trajectories without collisions. Such separatrix crossing leads to violation of the adiabaticity of the motion because the period of particle oscillations near the separatrix depends on the particle energy singularly (cf. Ref. 17).
To describe this violation of the adiabaticity we introduce the action as the area surrounded by the trajectory 25 and describe the dynamics of this action where the coordinate r max is defined as To introduce I as a new variable we use the generating function Ĥ ¼Ĥ 0 I; es; P s ð Þþ eH 1 I; /; es; P s ð Þ þ eĤ 1 I; /; es; P s ð Þ ; whereĤ 1 is H 1 from Eq. (7) written in new variables and we keep notation I forÎ. Instead of the variable conjugate to I, we use a phase / defined as / is a gyrophase, but it changes from one collision with the boundary to another one. In the new variables, the Hamiltonian H 1 takes the form Fig. 1 shows that for frozen values of (s, P s ) particles move along the periodic trajectories in the phase plane (r, p r ). For slowly changing (s, P s ), the area surrounded by these trajectories and boundary r ¼ 0 can be considered as an adiabatically conserved quantity. This invariant is conserved with an accuracy $ e. However, we can introduce the improved adiabatic invariant J which is conserved everywhere far from the separatrix with the accuracy $e 2 . The adiabatic invariant experiences small jumps in case of change of type of motion. Such change appears when a charged particle moving along the boundary (and reflecting from it) takes off from the boundary and starts drifting along the contours B(r, s) ¼ const. In analogy with smooth systems (without collisions, see Ref. 32), the change of type of motion could result in a jump of I (see Refs. 17 and 34). In case when a particle comes to the boundary again after a certain time, jumps of the adiabatic invariant may result in stochastization of charged particle motion. 34 To estimate the efficiency of this stochastization, one should derive the expression for the jump of the adiabatic invariant.
To derive the expression for the jump, use the improved adia- and brackets denote the averaging over / oscillations. 32 Function u can be calculated at any point along the r-axis. For the point r ¼ 0 (i.e., at the boundary), Eq. (14) takes the form where and T is the period of particle oscillation in the (r, p r ) plane (the time interval between two collisions with the boundary). Thus, to derive the expression for jump DJ, one should consider two terms: DI and Du.

III. JUMP OF THE ADIABATIC INVARIANT
The main variation of the adiabatic invariant corresponds to particle motion near (and taking off from) the boundary. Thus, we consider the close vicinity of the moment of particle taking off from the boundary. This moment corresponds to the equation P s À ffiffiffiffiffi 2h p ¼ 0. Thus, following the general approach (see Ref. 32), we introduce a variable K as a difference of particle energy h and energy calculated at the separa- At the close vicinity of the boundary, K is small (at the boundary H ¼ H sep and K ¼ 0). To derive an expression for change of adiabatic invariant at time-scale of one period of oscillations at the (r, p r ) plane, we need to obtain formula for period of these oscillations T and rate of change eH of area S surrounded by the separatrix.
To obtain formula for T, we use the identity T ¼ 2pð@I=@hÞ ¼ 2pð@I=@KÞ and consider small enough j ¼ K/h (i.e., we consider particle motion near the separatrix). We evaluate integral (8) and expand it around j ¼ 0 Equation (19) provides the following expression for period of particle oscillations: where T sep is regular (not depending on j) part of period (i.e., period calculated along the separatrix), while a ¼ 2/ (3pb).
The value of H can be calculated as H ¼ fS; H 0 g ¼ ð@S=@esÞð@H 0 =@P s Þ À ð@S=P s Þð@H 0 =@esÞ. In the limit j ! 0, we have H % Àð@S=@P s Þð@H 0 =@esÞ where (see Ref. 32) Thus, we have eH ¼ Àpb 0 s ð2hÞ 3=2 =b 2 . Then, the same procedure should be performed with u term given by Eq. (15). For term h H 1 i, we have where / Ã ¼ arcsinð ffiffiffiffiffiffiffiffiffiffiffi 1 À j p Þ. Thus, expansion of Eq. (22) around j ¼ 0 gives For the second term in Eq. (15), one can find Thus, expression for J around j ¼ 0 takes the form Let us estimate effects of variations of terms I and u on jump of the improved adiabatic invariant J. One can show that for j տ e the change of J for one period of the oscillations (i.e., between two collisions) is a value of order e 3 /j 3=2 . The value of j decays by eH between two collisions. Thus, the change of J when j decays from a value $1 to a value $e is $e 2 =e 1=2 $ e 3=2 . The change of J when j decays from a value $1 to a value j Ã ) e is $e 2 =j 1=2 Ã ( e 3=2 . Thus, for calculation of asymptotic of change of J, it is enough to consider motion when j j * . However, the change of eu for this motion is $ej * at most. Take e ( j * ( e 1=2 . Then, the change of eu for this motion is (e 3=2 . Thus, we can neglect the variation of u in our calculation of asymptotic of jump of J and consider only variation of I for motion in the domain of small j. For one period of oscillations in the (r, p r ) plane, we can write the following expression for change of 2pI ¼ S þ h Ð j 0 Tðhj; s; P s Þdj (change between ith and (i þ 1)th collisions with boundary; the index i decreases as the particle moves to the position of taking off): TðH sep þ hj; s; P s Þdj; (26) where T is given by Eq. (20), while the change of the area surrounded by the separatrix DS is and eH is a rate of S variation calculated for the position of taking off, while T i is the period calculated in the ith collisions with boundary. Substituting Eq. (20) to Eq. (26), we get As h(j i -j iþ1 ) ¼ ÀeH (see Ref. 4), Eq. (28) can be rewritten as We calculate the sum of changes (29) starting from the period i ¼ N -1 far from the separatrix and up to period i ¼ 0 (the last period before the separatrix crossing) We also should add the last period of oscillations including the separatrix crossing where T sep is the regular part of the period, i.e., the first term in Eq. (20). The regular part of integral compensates the regular part of DS À1 . Sum of Eq. (30) and (31) gives We introduce the new variable n ¼ hj 0 =ðeHÞ 2 ð0; 1Þ and use expression hj i ¼ hj 0 þ ieH to rewrite Eq. (32) where Thus, the total jump of the adiabatic invariant can be written as The function f(n) can be written in the integral form 17 The final dimensional expression for the jump DJ is where DJ, A, and H are written in dimensional form (we take into account that ffiffiffiffiffi 2h p ¼ 1 when h 0 is equal to the initial value of energy H) Equation (32) shows that DJ $ e 3=2 . The variable n changes significantly in case of small variations of initial particle position far from the position of taking off. Thus, n is traditionally considered as a quasi-random variable with uniform distribution within the range (0, 1) (see Ref. 32).

IV. NUMERICAL SIMULATIONS
To verify Eq. (37), we numerically integrate trajectories of charged particles interacting with the boundary. Equation (37) shows that DJ does not depend on the boundary curvature k(s); thus, we restrict our consideration by a simplified system with the straight boundary. We introduce Cartesian coordinates (x, y) and conjugate momenta (p x , p y ) for the dimensionless Hamiltonian where b(ey) ¼ 1 þ (ey) 2 for x > 0, and energy h ¼ 1/2. The improved adiabatic invariant J coincides with I at the boundary x ¼ 0; thus, we use I defined at x ¼ 0 instead of calculation of function u. An example of the particle trajectory in this system is shown in Fig. 2. The particle moves along the boundary x ¼ 0 reflecting from it and rotating in the magnetic field. This motion corresponds to the increase of jyj, i.e., the particle drifts toward the negative y values. Due to the increase of b with jyj, the particle Larmor radius becomes smaller and the distance between the boundary and the center of the Larmor circle becomes larger during the drift along the boundary. Around y % À100, the particle takes off from the boundary. This moment is shown in detail in the middle panel of Fig. 2. After taking off from the boundary, the FIG. 2. The particle trajectory in the (x, y) plane (left panel) and the moment of particle taking off (middle panel) are shown for system (39) with e ¼ 10 À2 . Right panels show p y and I profiles along the trajectory. particle starts rotating around magnetic field without collisions with the boundary. This rotation is much faster than the drift along the boundary, while the boundary drift is much faster than the gradient drift of particles far from the boundary. Thus, the particle trajectory looks like a sequence of circles slightly shifted one relative to other (see the middle panel of Fig. 2). The moment of taking off corresponds to p y ¼ 1 (i.e., d ¼ 1 À p y = ffiffiffiffiffi 2h p ¼ 0). We show both profiles of p y and the adiabatic invariant I to demonstrate that I starts changing before the p y ¼ 1 point, while for p y > 1 the adiabatic invariant is constant with a high accuracy. This change of the adiabatic invariant between the initial value and final (for p y > 1) value is the jump of J.
To check the dependence of DJ on the small parameter e and on quasi-random variable n, we compare numerical results and the analytical expression for DJ. To this end, we choose three values of e, generate an ensemble of 10 6 particles with different n, and trace these particles backward and forward in time. Then, for each trajectory, we determine the jump of the adiabatic invariant and normalize it on the coefficient ÀA(eH) 3=2 . For trajectories with the same n, we calculate the value of f as a ratio of DJ and ÀA(eH) 3=2 (see Eq. (37)). Comparison of such numerical functions f(n) and integral expression (36) is shown in Fig. 3. One can see that for small enough e the numerical calculations describe analytical expressions quite well.

V. DISCUSSION AND CONCLUSIONS
The violation of adiabaticity described by Eq. (37) is related to a non-smooth dependence of the period of oscillations on energy near the separatrix: T % T sep À 3pa ffiffiffi j p , Eq. (20). For separatrix crossings in smooth Hamiltonian systems considered in Refs. 9, 32, 33, and 46, the period has a logarithmic singularity on a separatrix, T $ lnj. The corresponding jump in adiabatic invariant is DJ $ e. In our case, the singularity is weaker, $ ffiffiffi j p , what results in a smaller jump: DJ $ e 3=2 .
Besides the general theoretical interest, obtained results can be important for description of processes in cosmic plasma systems. Equation (37) shows that the jump of the adiabatic invariance in the system with strong magnetic field and impacts is about $e 3=2 . The amplitude of this jump is significantly larger than the rate of the magnetic moment change in system with strong magnetic field without a change of type of the particle motion, DJ $ expðÀ1=eÞ (see Refs. 14 and 47). There are several space plasma systems where magnetic field configuration assumes the billiard-like motion with the change of type of motion (drift along the boundary and gyrorotation at some distance from the boundary). These are magnetic structures with a local minimum of the magnetic field inside some volume and the strong magnetic field outside it (e.g., plasmoids, 40,45 magnetic holes, 5,24,48 etc.). The particle (electron) motion inside such structures resembles the billiard motion. Thus, our theory predicts that for a time interval Dt $ 1/e a change of the magnetic moment for such particles is about e 3=2 . The averaged jump hDJi n is equal to zero, thus, the rate of the adiabatic invariant diffusion is determined as hðDJÞ 2 i n =Dt $ e 4 . As a result, for a time interval $ e À4 all particles inside such magnetic structures should be significantly mixed due to the destruction of initial magnetic momenta.
To conclude, we have derived the equation for the rate of adiabatic invariance destruction DJ (33, 37) for magnetic billiards. We have shown that DJ does not depend on the boundary curvature k(s) as well as on the magnetic field inhomogeneity across the boundary @B/@r. The scaling of DJ $ e 3=2 and dependence of DJ on n are similar with the scaling and dependence DJ(n) derived before for a potential system with impacts. 17 Thus, there is a certain universality of the process of adiabaticity destructions in the system with impacts. However, in contrast to the system, 17 in the billiard with the inhomogeneous magnetic field, many takes off from the boundary (i.e., many jumps of the adiabatic invariant) result in a stochastization of motion in large phase space domain. 34