- Email: [email protected]
Numerical solution of a reduced model of collisionless magnetic reconnection in two and three dimensions
Computer Physics Communications 164 (2004) 23–28 www.elsevier.com/locate/cpc
Numerical solution of a reduced model of collisionless magnetic reconnection in two and three dimensions D. Grasso a,∗ , D. Borgogno a , F. Califano b , D. Farina c , F. Pegoraro b , F. Porcelli a a Burning Plasma Research Group, INFM, Politecnico di Torino, Italy b Dipartimento di Fisica, Università di Pisa, INFM, Italy c Istituto di Fisica del Plasma, CNR, EURATOM-ENEA-CNR Assoc., Milano, Italy
Available online 23 July 2004
Abstract Magnetic reconnection in collisionless regime is investigated in two- and three-dimensional configurations. In the framework of two-dimensional configurations some recent results concerning the difference between dissipative and collisionless reconnection are reviewed. The results of numerical simulation of three-dimensional configurations are presented and analyzed. 2004 Elsevier B.V. All rights reserved. PACS: 52.35.P; 47.65.+a; 94.30.Gm; 52.65
1. Introduction Magnetic field line reconnection in high temperature plasmas is one of the frontier problems in plasma physics due to its relevance to astrophysical and laboratory plasmas. In collisionless fluid regimes, the topology of the magnetic field is broken by the effect of electron inertia, and the reconnection process exhibits Hamiltonian properties. Our analysis is based on a two-fluid Hamiltonian model, derived in Ref. [1]. Numerical integrations of the adopted model in 2D configurations show the formation of current density and vorticity layers aligned with the separatrix of the * Corresponding author. INFM-Politecnico di Torino, Diparti-
mento di Energetica, Corso Duca degli Abruzzi 24, 10129 Torino, Italy. Tel.: +39 0115644435; Fax: +39 0115644499. E-mail address: [email protected] (D. Grasso). 0010-4655/$ – see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2004.06.004
magnetic island. These structures have been interpreted on the basis of the Lagrangian invariants [3]. In the nonlinear phase the magnetic island width saturates at a macroscopic amplitude superimposed to the fine spatial scales arising from the phase mixing of the Lagrangian invariants [4]. The removal of the 2D constraint is expected to show new phenomena. In particular, it leads to coupling between modes with different helicities and with different resonant surfaces and to magnetic field line stochasticity. In addition, 3D perturbations are not constrained by the same topological conservations as in 2D. We will see that the nonlinear interaction between two linearly unstable, resonant modes with different helicities, drives higher order harmonics which modify the nonlinear structure of the reconnection region significantly. In this case the magnetic flux function
24
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
depends on all the three spatial coordinates. Thus, the Hamiltonian function of the equivalent dynamical system that describes the spatial structure of the magnetic field lines is no longer integrable, as in the case of single helicity perturbations, and we observe, by Poincaré maps and Lyapunov exponents, the development of magnetic field line stochasticity.
2. Governing equations We consider here the dissipationless model given in Ref. [1], which describes drift-Alfvèn perturbations in a plasma immersed in a strong, uniform, externally imposed magnetic field, where spatial variations along the field lines are assumed to be small compared to perpendicular variations. In the limit of small ion gyroradii, this system consists of two-fluid, quasi neutral electron and ion equations where small scales effects related to the electron temperature, ρs = (Te /Ti )1/2 i , and to electron inertia, de = c/ωpe , are retained, but magnetic curvature effects are neglected. The √ equations, normalized to the Alfvèn time, τA = 4πnmi Lx /By0 , with By0 the characteristic magnitude of the equilibrium magnetic field, and to the macroscopic scale length, L, are: ∂(ψ + de2 J ) + [ϕ, ψ + de2 J ] − ρs2 [U, ψ] ∂t ∂ϕ ∂U = − ρs2 , ∂z ∂z ∂J ∂U + [ϕ, U ] − [J, ψ] = − ∂t ∂z
(1) (2)
where ψ and ϕ are the magnetic flux and the stream 2 ψ is the current function, respectively. J = −∇⊥ 2 density and U = ∇⊥ ϕ is the vorticity. The Poisson brackets are defined as [A, B] = ez · ∇⊥ A × ∇⊥ B. We adopt here the notation of Ref. [4]. Adding and subtracting the above equations with straigthforward algebra the alternative equations can be obtained: ∂(ϕ± ∓ des G± ) ∂G± + [ϕ± , G± ] = ∂t ∂z
which the functions G± are advected with, are zdependent; at the same time, a z-dependent source term for G± appears at the right-hand side of Eq. (3). After the formation of the magnetic island, very small scales develop due to nonlinear interactions and the “collisionless” simulation becomes computationally challenging. We developed a numerical code based on a finite volumes scheme where filters, based on a FFT algorithm, have been introduced acting only on typical length scales much smaller than any other physical length scale of the system. These filters, which start to progressively smooth out the small scales starting from a characteristic length scale, but leaving unchanged the large scale dynamics even on long times, are described in Ref. [5] (in particular we have used those corresponding to formula C.2.7 in the appendix). The code is based on a finite volume scheme, in which the averaged values of the G± fields are advanced in time using an explicit third order Adams–Bashfort scheme. Since the boundary conditions are periodic, a Fast Fourier Transform method is used to reconstruct the fields at the grid points. In order to address the three-dimensional problem, the code has been parallelized adopting the MPI libraries. In all the numerical investigations of collisionless reconnection we present here, we restrict ourselves to the so-called large ∆ regime, defined by ∆ de > (de /ρs )1./3 , where ∆ is the standard instability parameter. This choice is motivated by the fact that in this regime magnetic reconnection can take place on a fast time scale [2,6].
(3)
where G± = ψ + de J ± de s U , ϕ± = ϕ ± (s /de )ψ. We point out that, in the 2D limit, where ∂/∂z vanishes, G± are Lagrangian invariants. 3D effects come from the fact that the velocities, v± = ez × ∇ϕ± ,
3. 2D results We choose a configuration with a strong magnetic field in the z-direction. The equilibrium magnetic field is B = B0 ez + ∇ψeq × ez , where B0 is constant and ψeq = 1/ cosh2 (x). With this choice, ψeq and its derivatives go to zero at infinity. So, by taking a large integration domain, we can impose periodic boundary conditions also in the inhomogeneity direction. We investigated cases in which ρs is of the order of de . During the evolution of the reconnection process, after a linear phase and an early nonlinear phase in which the process exhibits a quasi-explosive behavior [2,3], a saturated island is reached [4]. In Fig. 1 (top panels) the contour plots of the G+ field
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
25
Fig. 1. Contour plots of the G+ field (top frames) at two different simulation times, show the phase mixing and the filamentation, typical of collisionless reconnection. The formation of small scale is also evident in the contour plots of the current density and vorticity (bottom frames). Superimposed to each plot is the magnetic island.
is drawn at two different simulation times. Superimposed to these contour plots is the magnetic island. We point out that, while the topology of the magnetic flux changes, the topology of the G± fields is preserved. The decoupling between the Lagrangian invariants G± and the reconnected flux ψ is made possible by the formation of current density and vorticity layers. As shown in Fig. 1 (bottom panel), these layers have a cross-shaped structure, aligned to the magnetic island [3]. These layers are accompanied by the formation of increasingly small scales, that generate from the “phase mixing” of the Lagrangian invariant, which are advected by velocity fields v± , along two mirror-symmetric rotation patterns inside the magnetic island. Part of the magnetic energy is transferred to small scale structures in the conserved fields, which
are averaged out in the magnetic flux. This fact allows to access a new equilibrium with a macroscopic island, in spite of the energy conservation property for a Hamiltonian system [4].
4. 3D results We performed our simulations in a 3D periodic slab geometry, starting from a static equilibrium configuration characterized by the following magnetic field, Beq = B0 ez + ∇ψeq × ez ,
(4)
where ψeq = α cos(x) and B0 = −1. All the runs we performed have been carried out with α = 0.48 and s = de = 0.24. The integration domain is defined
26
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
by −Lx /2 < x < Lx /2, −Ly /2 < y < Ly /2 and −Lz /2 < z < Lz /2, with Lx = 2π , Ly = 4π and Lz = 32π , using nx = ny = nz = 64 modes along x, y, z, respectively. The presence of two linearly unstable resonant modes with different helicity modifies significantly the nonlinear behaviour of the system. Indeed, while in the linear phase each mode evolves independently from the other, justifying the single helicity approximation, nonlinearly a strong interaction between the different helicity modes appears. To see that, we performed a numerical simulation starting from the magnetic ˜ ˆ perturbation ψ(x, y, z) = ψ(x)(exp(ik y y + ikz z) + exp(iky y − ikz z)), characterized by the couple of wave numbers ky = 2π/Ly , kz = 2π/Lz . Our numerical results show the presence of perturbations with helicities different from those initially imposed, whose growth rates are in good agreement with the quasi-linear theory predictions, as shown in Fig. 2. Here, as a measure of the growth rate, we plot the time derivative of the logarithm of different modes of the current density evaluated at the rational surfaces. We note that all the modes of order two with respect to the initially imposed perturbed modes (1, ±1), have the same growth rate. This growth rate is twice the one of the modes (1, ±1).
In the linear phase, the resonant magnetic surfaces, where the island’s X and O-points are located, lie at x = xs satisfying
Fig. 2. The time derivative of log |Jm,n (x, t)| dx at the rational surface, are plotted. The second order modes grow, according to quasi linear theory, with rates equal to two times the growth rate of the first order mode (1, 1).
Fig. 3. Current density profiles at different simulation times. During the linear phase (top panels) the current density peaks are localized at the rational surfaces defined by Eq. (5). In the nonlinear phase (bottom panels) the two rational surfaces move towards each-others until the current density peaks merge.
∂ψeq kz =± . ∂x ky
(5)
Nonlinearly, the resonant surfaces move along x-direction carrying along the current density and vorticity structures. This behavior is illustrated in Fig. 3, where the time evolution of the perturbed current density profile as a function of x and at fixed values of y and z, is shown. Initially two peaks are localized at the surfaces defined by Eq. (5). When the nonlinear phase is entered, the two peaks start to move towards each-other until they merge. The strong coupling of modes with different helicities is responsible for the change in the structure of the current density and vorticity layers during the nonlinear evolution. Fig. 4, which shows a density current isosurface and the superimposed contour plots of the same field at the planes x = π , y = 2π , z = −16π , gives a qualitative idea of the topological modification of the current density layer during the transition from
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
27
Fig. 4. Isosurfaces and contour plots of the current density field, J , at two different simulation times.
the linear to the nonlinear phase. The first frame shows that, in the linear phase, the current layer is characterized by the presence of two separate channels. They correspond to the initial double helicities perturbation. The second frame shows that in the nonlinear phase, due to the merging of the two current channels, the structure has changed its topology and looks like a single channel. In the presence of more than one helicity, the Hamiltonian for the magnetic field lines, ψ(x, y, z, t), is no longer integrable. Magnetic field lines satisfy the following set of Hamilton equations, 1 ∂ψ Bx ∂x =− = , (6) ∂z B0 ∂y B0 By 1 ∂ψ ∂y = = (7) ∂z B0 ∂x B0 where x and y play the role of conjugate variables, z the role of time and the true time t is a parameter. Numerical integration of the above equations allowed us to produce Poincaré maps for the magnetic field lines. Fig. 5 shows the maps at t = 50, 60, 90τA and t = 110τA , for the section z = 0. Until t = 50τA the magnetic field lines are well defined over most of all the integration domain. After this time, a sizable stochastic region becomes apparent near the separatrices of the two magnetic islands. At t =
110τA stochasticity has developed so much, that only a few KAM magnetic surfaces are preserved. Most field lines look chaotic and the section consists of a nearly uniform splatter of dots. These figures are in qualitative agreement with the so called “overlap criterion” formulated by Chirikov [7]. Calculating the Lyapunov exponents we evaluated the degree of stochasticity for the magnetic field lines at the different times. At the end of the simulation we observed that field lines close together separate rapidly, with a mean value of the Lyapunov coefficient grater than 1.
5. Conclusions We have presented the numerical solution of a Hamiltonian model for magnetic reconnection. In the two-dimensional limit the process has been followed until a saturated island has been reached. The saturation mechanism is interpreted in terms of the phase mixing of the Lagrangian invariants of the model. In the three-dimensional case we see that, in the early nonlinear phase our results are in good agreement with quasilinear theory. When the full nonlinear regime is entered, a strong modification of the current density and vorticity layers is observed, as a result of the cou-
28
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
Fig. 5. Poincaré maps at different simulation times for a double helicity case.
pling between different helicity modes. The magnetic field structure shows a high degree of stochasticity characterized by the Lyapunov exponents. References [1] T.J. Schep, et al., Phys. Plasmas 1 (1994) 2843. [2] M. Ottaviani, F. Porcelli, Phys. Rev. Lett. 71 (23) (1993) 3802.
[3] [4] [5] [6]
E. Cafaro, et al., Phys. Rev. Lett. 80 (1998) 20. D. Grasso, et al., Phys. Rev. Lett. 86 (2001) 5051. S.K. Lele, J. Comput. Phys. 103 (1992) 16. A. Aydemir, Phys. Fluids B 4 (1992) 3469; R.G. Kleva, et al., Phys. Plasmas 2 (1995) 23; X. Wang, A. Bhattacharjee, Phys. Rev. Lett. 70 (1993) 1627. [7] B.V. Chirikov, Phys. Reports 52 (1979) 265.
Numerical solution of a reduced model of collisionless magnetic reconnection in two and three dimensions D. Grasso a,∗ , D. Borgogno a , F. Califano b , D. Farina c , F. Pegoraro b , F. Porcelli a a Burning Plasma Research Group, INFM, Politecnico di Torino, Italy b Dipartimento di Fisica, Università di Pisa, INFM, Italy c Istituto di Fisica del Plasma, CNR, EURATOM-ENEA-CNR Assoc., Milano, Italy
Available online 23 July 2004
Abstract Magnetic reconnection in collisionless regime is investigated in two- and three-dimensional configurations. In the framework of two-dimensional configurations some recent results concerning the difference between dissipative and collisionless reconnection are reviewed. The results of numerical simulation of three-dimensional configurations are presented and analyzed. 2004 Elsevier B.V. All rights reserved. PACS: 52.35.P; 47.65.+a; 94.30.Gm; 52.65
1. Introduction Magnetic field line reconnection in high temperature plasmas is one of the frontier problems in plasma physics due to its relevance to astrophysical and laboratory plasmas. In collisionless fluid regimes, the topology of the magnetic field is broken by the effect of electron inertia, and the reconnection process exhibits Hamiltonian properties. Our analysis is based on a two-fluid Hamiltonian model, derived in Ref. [1]. Numerical integrations of the adopted model in 2D configurations show the formation of current density and vorticity layers aligned with the separatrix of the * Corresponding author. INFM-Politecnico di Torino, Diparti-
mento di Energetica, Corso Duca degli Abruzzi 24, 10129 Torino, Italy. Tel.: +39 0115644435; Fax: +39 0115644499. E-mail address: [email protected] (D. Grasso). 0010-4655/$ – see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2004.06.004
magnetic island. These structures have been interpreted on the basis of the Lagrangian invariants [3]. In the nonlinear phase the magnetic island width saturates at a macroscopic amplitude superimposed to the fine spatial scales arising from the phase mixing of the Lagrangian invariants [4]. The removal of the 2D constraint is expected to show new phenomena. In particular, it leads to coupling between modes with different helicities and with different resonant surfaces and to magnetic field line stochasticity. In addition, 3D perturbations are not constrained by the same topological conservations as in 2D. We will see that the nonlinear interaction between two linearly unstable, resonant modes with different helicities, drives higher order harmonics which modify the nonlinear structure of the reconnection region significantly. In this case the magnetic flux function
24
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
depends on all the three spatial coordinates. Thus, the Hamiltonian function of the equivalent dynamical system that describes the spatial structure of the magnetic field lines is no longer integrable, as in the case of single helicity perturbations, and we observe, by Poincaré maps and Lyapunov exponents, the development of magnetic field line stochasticity.
2. Governing equations We consider here the dissipationless model given in Ref. [1], which describes drift-Alfvèn perturbations in a plasma immersed in a strong, uniform, externally imposed magnetic field, where spatial variations along the field lines are assumed to be small compared to perpendicular variations. In the limit of small ion gyroradii, this system consists of two-fluid, quasi neutral electron and ion equations where small scales effects related to the electron temperature, ρs = (Te /Ti )1/2 i , and to electron inertia, de = c/ωpe , are retained, but magnetic curvature effects are neglected. The √ equations, normalized to the Alfvèn time, τA = 4πnmi Lx /By0 , with By0 the characteristic magnitude of the equilibrium magnetic field, and to the macroscopic scale length, L, are: ∂(ψ + de2 J ) + [ϕ, ψ + de2 J ] − ρs2 [U, ψ] ∂t ∂ϕ ∂U = − ρs2 , ∂z ∂z ∂J ∂U + [ϕ, U ] − [J, ψ] = − ∂t ∂z
(1) (2)
where ψ and ϕ are the magnetic flux and the stream 2 ψ is the current function, respectively. J = −∇⊥ 2 density and U = ∇⊥ ϕ is the vorticity. The Poisson brackets are defined as [A, B] = ez · ∇⊥ A × ∇⊥ B. We adopt here the notation of Ref. [4]. Adding and subtracting the above equations with straigthforward algebra the alternative equations can be obtained: ∂(ϕ± ∓ des G± ) ∂G± + [ϕ± , G± ] = ∂t ∂z
which the functions G± are advected with, are zdependent; at the same time, a z-dependent source term for G± appears at the right-hand side of Eq. (3). After the formation of the magnetic island, very small scales develop due to nonlinear interactions and the “collisionless” simulation becomes computationally challenging. We developed a numerical code based on a finite volumes scheme where filters, based on a FFT algorithm, have been introduced acting only on typical length scales much smaller than any other physical length scale of the system. These filters, which start to progressively smooth out the small scales starting from a characteristic length scale, but leaving unchanged the large scale dynamics even on long times, are described in Ref. [5] (in particular we have used those corresponding to formula C.2.7 in the appendix). The code is based on a finite volume scheme, in which the averaged values of the G± fields are advanced in time using an explicit third order Adams–Bashfort scheme. Since the boundary conditions are periodic, a Fast Fourier Transform method is used to reconstruct the fields at the grid points. In order to address the three-dimensional problem, the code has been parallelized adopting the MPI libraries. In all the numerical investigations of collisionless reconnection we present here, we restrict ourselves to the so-called large ∆ regime, defined by ∆ de > (de /ρs )1./3 , where ∆ is the standard instability parameter. This choice is motivated by the fact that in this regime magnetic reconnection can take place on a fast time scale [2,6].
(3)
where G± = ψ + de J ± de s U , ϕ± = ϕ ± (s /de )ψ. We point out that, in the 2D limit, where ∂/∂z vanishes, G± are Lagrangian invariants. 3D effects come from the fact that the velocities, v± = ez × ∇ϕ± ,
3. 2D results We choose a configuration with a strong magnetic field in the z-direction. The equilibrium magnetic field is B = B0 ez + ∇ψeq × ez , where B0 is constant and ψeq = 1/ cosh2 (x). With this choice, ψeq and its derivatives go to zero at infinity. So, by taking a large integration domain, we can impose periodic boundary conditions also in the inhomogeneity direction. We investigated cases in which ρs is of the order of de . During the evolution of the reconnection process, after a linear phase and an early nonlinear phase in which the process exhibits a quasi-explosive behavior [2,3], a saturated island is reached [4]. In Fig. 1 (top panels) the contour plots of the G+ field
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
25
Fig. 1. Contour plots of the G+ field (top frames) at two different simulation times, show the phase mixing and the filamentation, typical of collisionless reconnection. The formation of small scale is also evident in the contour plots of the current density and vorticity (bottom frames). Superimposed to each plot is the magnetic island.
is drawn at two different simulation times. Superimposed to these contour plots is the magnetic island. We point out that, while the topology of the magnetic flux changes, the topology of the G± fields is preserved. The decoupling between the Lagrangian invariants G± and the reconnected flux ψ is made possible by the formation of current density and vorticity layers. As shown in Fig. 1 (bottom panel), these layers have a cross-shaped structure, aligned to the magnetic island [3]. These layers are accompanied by the formation of increasingly small scales, that generate from the “phase mixing” of the Lagrangian invariant, which are advected by velocity fields v± , along two mirror-symmetric rotation patterns inside the magnetic island. Part of the magnetic energy is transferred to small scale structures in the conserved fields, which
are averaged out in the magnetic flux. This fact allows to access a new equilibrium with a macroscopic island, in spite of the energy conservation property for a Hamiltonian system [4].
4. 3D results We performed our simulations in a 3D periodic slab geometry, starting from a static equilibrium configuration characterized by the following magnetic field, Beq = B0 ez + ∇ψeq × ez ,
(4)
where ψeq = α cos(x) and B0 = −1. All the runs we performed have been carried out with α = 0.48 and s = de = 0.24. The integration domain is defined
26
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
by −Lx /2 < x < Lx /2, −Ly /2 < y < Ly /2 and −Lz /2 < z < Lz /2, with Lx = 2π , Ly = 4π and Lz = 32π , using nx = ny = nz = 64 modes along x, y, z, respectively. The presence of two linearly unstable resonant modes with different helicity modifies significantly the nonlinear behaviour of the system. Indeed, while in the linear phase each mode evolves independently from the other, justifying the single helicity approximation, nonlinearly a strong interaction between the different helicity modes appears. To see that, we performed a numerical simulation starting from the magnetic ˜ ˆ perturbation ψ(x, y, z) = ψ(x)(exp(ik y y + ikz z) + exp(iky y − ikz z)), characterized by the couple of wave numbers ky = 2π/Ly , kz = 2π/Lz . Our numerical results show the presence of perturbations with helicities different from those initially imposed, whose growth rates are in good agreement with the quasi-linear theory predictions, as shown in Fig. 2. Here, as a measure of the growth rate, we plot the time derivative of the logarithm of different modes of the current density evaluated at the rational surfaces. We note that all the modes of order two with respect to the initially imposed perturbed modes (1, ±1), have the same growth rate. This growth rate is twice the one of the modes (1, ±1).
In the linear phase, the resonant magnetic surfaces, where the island’s X and O-points are located, lie at x = xs satisfying
Fig. 2. The time derivative of log |Jm,n (x, t)| dx at the rational surface, are plotted. The second order modes grow, according to quasi linear theory, with rates equal to two times the growth rate of the first order mode (1, 1).
Fig. 3. Current density profiles at different simulation times. During the linear phase (top panels) the current density peaks are localized at the rational surfaces defined by Eq. (5). In the nonlinear phase (bottom panels) the two rational surfaces move towards each-others until the current density peaks merge.
∂ψeq kz =± . ∂x ky
(5)
Nonlinearly, the resonant surfaces move along x-direction carrying along the current density and vorticity structures. This behavior is illustrated in Fig. 3, where the time evolution of the perturbed current density profile as a function of x and at fixed values of y and z, is shown. Initially two peaks are localized at the surfaces defined by Eq. (5). When the nonlinear phase is entered, the two peaks start to move towards each-other until they merge. The strong coupling of modes with different helicities is responsible for the change in the structure of the current density and vorticity layers during the nonlinear evolution. Fig. 4, which shows a density current isosurface and the superimposed contour plots of the same field at the planes x = π , y = 2π , z = −16π , gives a qualitative idea of the topological modification of the current density layer during the transition from
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
27
Fig. 4. Isosurfaces and contour plots of the current density field, J , at two different simulation times.
the linear to the nonlinear phase. The first frame shows that, in the linear phase, the current layer is characterized by the presence of two separate channels. They correspond to the initial double helicities perturbation. The second frame shows that in the nonlinear phase, due to the merging of the two current channels, the structure has changed its topology and looks like a single channel. In the presence of more than one helicity, the Hamiltonian for the magnetic field lines, ψ(x, y, z, t), is no longer integrable. Magnetic field lines satisfy the following set of Hamilton equations, 1 ∂ψ Bx ∂x =− = , (6) ∂z B0 ∂y B0 By 1 ∂ψ ∂y = = (7) ∂z B0 ∂x B0 where x and y play the role of conjugate variables, z the role of time and the true time t is a parameter. Numerical integration of the above equations allowed us to produce Poincaré maps for the magnetic field lines. Fig. 5 shows the maps at t = 50, 60, 90τA and t = 110τA , for the section z = 0. Until t = 50τA the magnetic field lines are well defined over most of all the integration domain. After this time, a sizable stochastic region becomes apparent near the separatrices of the two magnetic islands. At t =
110τA stochasticity has developed so much, that only a few KAM magnetic surfaces are preserved. Most field lines look chaotic and the section consists of a nearly uniform splatter of dots. These figures are in qualitative agreement with the so called “overlap criterion” formulated by Chirikov [7]. Calculating the Lyapunov exponents we evaluated the degree of stochasticity for the magnetic field lines at the different times. At the end of the simulation we observed that field lines close together separate rapidly, with a mean value of the Lyapunov coefficient grater than 1.
5. Conclusions We have presented the numerical solution of a Hamiltonian model for magnetic reconnection. In the two-dimensional limit the process has been followed until a saturated island has been reached. The saturation mechanism is interpreted in terms of the phase mixing of the Lagrangian invariants of the model. In the three-dimensional case we see that, in the early nonlinear phase our results are in good agreement with quasilinear theory. When the full nonlinear regime is entered, a strong modification of the current density and vorticity layers is observed, as a result of the cou-
28
D. Grasso et al. / Computer Physics Communications 164 (2004) 23–28
Fig. 5. Poincaré maps at different simulation times for a double helicity case.
pling between different helicity modes. The magnetic field structure shows a high degree of stochasticity characterized by the Lyapunov exponents. References [1] T.J. Schep, et al., Phys. Plasmas 1 (1994) 2843. [2] M. Ottaviani, F. Porcelli, Phys. Rev. Lett. 71 (23) (1993) 3802.
[3] [4] [5] [6]
E. Cafaro, et al., Phys. Rev. Lett. 80 (1998) 20. D. Grasso, et al., Phys. Rev. Lett. 86 (2001) 5051. S.K. Lele, J. Comput. Phys. 103 (1992) 16. A. Aydemir, Phys. Fluids B 4 (1992) 3469; R.G. Kleva, et al., Phys. Plasmas 2 (1995) 23; X. Wang, A. Bhattacharjee, Phys. Rev. Lett. 70 (1993) 1627. [7] B.V. Chirikov, Phys. Reports 52 (1979) 265.