The nonlinear dynamics of collisionless reconnecting modes is investigated, in the framework of a three-dimensional gyrofluid model. This is the relevant regime of high-temperature plasmas, where reconnection is made possible by electron inertia and has higher growth rates than resistive reconnection. The presence of a strong guide field is assumed, in a background slab model wih Dirichlet boundary conditions in the direction of nonuniformity. Values of ion sound gyro-radius and electron collisionless skin depth much smaller than the current layer width are considered. Strong acceleration of growth is found at the onset to nonlinearity, while at all times the energy functional is well conserved. Nonlinear growth rates more than one order of magnitude higher than linear growth rates are observed when entering into the small- regime.
Observation of explosive collisionless reconnection
in 3D nonlinear gyrofluid simulations.
A. Biancalani, Bruce D. Scott
Max-Planck-Institut für Plasmaphysik, Euratom Association, D-85748 Garching, Germany
Max-Planck-Institut für Sonnensystemforschung, Katlenburg-Lindau, Germany
webaddress of first author: www.ipp.mpg.de/~biancala
Magnetic reconnection occurs in space and laboratory plasmas, when the equilibrium magnetic field lines break due to non ideal effects and reconnect with different magnetic topology . Here, we consider the regime of high-temperature plasmas in the presence of a strong guide magnetic field, where reconnection is made possible by electron inertia and has higher growth rates than resistive reconnection. Such a regime is named collisionless reconnection [2, 3].
Two dimensional studies of collisionless reconnection with gyrofluid models were performed in Ref.  with code REC2, in a periodic configuration with flat temperature and density equilibrium profiles, and more recently in Ref.  and Ref. . Development of very narrow current layers was also found in the nonlinear phase of collisionless reconnection [6, 7, 8, 9]. This suggests the importance of affording high spatial resolution to investigate small scale dynamics, which is crucial in the nonlinear phase of the magnetic island growth. Three dimensional studies with two-fluid model in periodic configuration were also performed in a large- regime  (where is a parameter measuring the ratio of mode wavelength and current layer width). Large- regimes, corresponding to large wavelengths, were often adopted in past investiagations for their apparent property of having higher growth rates, and therefore of explaining better fast reconnection observed in nature. Nevertheless, a comprehensive theoretical model capable of reproducing such observed high growth rates does not exist at present day.
Nonlinear growth acceleration was analytically predicted for the m=1 mode in tokamaks, corresponding to large regimes in slab model . On the contrary, a nonlinear subexponential growth was predicted analytically for the collisionless tearing mode, corrisponding to small regimes in slab model . Numerically however, a nonlinear growth acceleration was found in fluid simulations of collisionless reconnection  even for small regimes . The nonlinear growth rate was shown to increase for smaller values of the ion sound gyro-radius - with the sound speed and the ion cyclotron frequency - and of the electron collisionless skin depth - with the light speed and the electron plasma frequency. The value of the nonlinear growth rate was found up to twice the value of the linear growth rate.
In this paper, we investigate the small- regime in the unexplored limit of very small and , where is the current width: the equilibrium characteristic scale length in the direction of nonuniformity. No gradients of the equilibrium pressure and temperature are set, therefore turbulence is not driven from the beginning of the simulations, but it is allowed to develop selfconsistently with the island evolution during subsequent phases of the simulations.
Numerical experiments are carried out with the GEM-code (Gyrofluid ElectroMagnetic ). This code, originally written for drift wave turbulence in tokamaks, adopts an electromagnetic full-finite-Larmor-radius gyrofluid model of equations for electrons and ions. Temperature fluctuations are included in the model. The GEM-code has been benchmarked with results of Ref.  and with analytical theory in large regime given by Ref. . The gyrofluid model is obtained by taking the first moments of the gyrokinetic equations [15, 16], with the advantage of incorporating finite ion gyroradius effects at arbitrary order. The zeroth and first moments of the gyrokinetic equation for ion and electron species (labeled here by ) are:
where the tilde symbol denotes perturbed gyrocenter densities and velocities , scalar potential and vector potential , and temperatures (cf. Eqs: 34-39, 83, 86, 99-105 of Ref. ). The advective time derivative is given by . The time is normalized in terms of , with the perpendicular scale length normalized to , and the parallel gradient normalized in terms of . Therefore, the size of the contravariant magnetic unit vector component is comparable to , where and and are the scale lengths along and perpendicular to the equilibrium magnetic field. The parameters , and are defined by: , and . The parallel gradient is calculated in the direction parallel to the total (equilibrium plus perturbed) magnetic field, and the gyro-averaged potentials are defined as and . The fields are advected with the the drift velocity and with (where for the electrons we use the bare potential and ). The FLR nonlinearities are represented both by the difference between and , and also by . The electromagnetic fields are linked to the sources by the polarization and induction equations:
where the Padé approximant forms are adopted : , , and is the ion gyro-radius ( here is the ion thermal velocity). Note that in the fluid limit, the gyrofluid equations are equivalent to a two-fluid model. As basic assumptions of this model, we neglect fast magneto-sonic waves, and we consider only frequencies smaller than the ion cyclotron frequency (therefore no whistlers are present in this model). The characteristic scale lengths of this model are the equilibrium current width , the ion gyro-radius , the ion sound gyro-radius , and the collisionless skin depth . The characteristic velocities are the Alfvén velocity and the sound speed , where is the ion particle density and the ion mass.
We consider the plasma regime with , and values of ranging from 0.41 to 4.1. We obtain a very thin reconnection layer with respect of the equilibrium current layer (). This requires a high spatial resolution: we use a grid with up to 2048 points in the perpendicular direction and to 512 in the direction of the sheared field . This also results in a low reconnecting growth rate, requiring very long runs. The parellization of the GEM-code into 64 processors solves the problem of the excessive time length of the runs, bringing the characteristic run time to a few days for the longest runs. A deuterium plasma is considered () in our simulations. The value of the parallel to perpendicular scale ratio has been chosen as . The ion to electron temperature ratio is considered unitary: .
The evolution of a magnetic island is studied in a Harris-pinch equilibrium, where the equilibrium magnetic field is a hyperbolic tangent of the nonuniformity coordinate . We initialize our simulations with a periodic perturbation with wavenumber , where is the box size in the direction of the shear magnetic field. We add to this initial condition a white noise with amplitude ten times smaller than the magnetic island perturbation. The white noise has components in all three directions, and allows the system - whose equilibrium is two-dimensional - to develop structures in the third direction too. Periodic boundary conditions are imposed in and in the direction of the equilibrium magnetic field , and Dirichlet boundary conditions in . The choice of the boundary conditions in differs from those of Ref. , where periodic buondary conditions are chosen also in the direction. The energy of the nonzonal components of the fluctuating fields is measured and the growth rate is calculated as the time derivative of the energy (see Fig. 1): . Such boundary conditions allows free energy to enter from the borders of the simulation box by means of a nonzero gradient of the scalar potential in the direction, corresponding to ExB drifts in the directions. The consistency of our choice of boundary conditions has been checked by doubling the size of the box in the perpendicular direction , and finding that the dynamics under investigation is not changed.
The initial perturbations are let evolve in time, and different phases in the evolution of the magnetic island are found, defined according to the functional behavior of the growth rate. During an initial transient phase, fluctuations with all wave numbers interact with the magnetic island seed. After the transient phase, a clear linear phase is found, where the magnetic island grows exponentially in time and the growth rate is constant (t 100 to t 300 in Fig. 1). At the onset to the nonlinear phase, the magnetic island amplitude is found to grow superexponentially (t 300 to t 350 in Fig. 1). During this acceleration phase, the growth rate increases, reaches a peak value, whose dependence on is shown in Fig. 2, then decreases again. After the nonlinear acceleration phase, the island reaches a saturated amplitude and nonlinearly oscillates around the saturated amplitude. This corresponds to a nonlinear oscillation of the growth rate around the null value. During the saturation phase, the energy flows to secondary instabilities which are formed around the island separatrix.
Simulations are repeated with different plasma parameters, and the value of the peak of the nonlinear growth rate is measured in the different cases. The nonlinear acceleration is found to become more prominent for higher values of . In particular, we repeat the simulation as in the case of Fig. 1, by varying and mantaining the same ratio of , and . In other words, for higher , we are choosing bigger box sizes with bigger equilibrium current sheet’s widths, whereas we keep the same value of . The separation in these different phases is found clearly in all these cases with different beta: we always find an initial transient phase, a linear phase, and a nonlinear acceleration followed by a saturation. The value of the peak of the linear and nonlinear growth rates normalized with the sound time for these different cases is shown in Fig. 2.
The linear growth rate is shown to decrease with a constant slope in the logarithmic graph for increasing , corresponding to a well defined power law. Increasing the value of by one order of magnitude, the linear growth rate decreases also by roughly one order of magnitude. This is consistent with analytical predictions, describing the reconnection process as occurring in a thin “inertial” layer with relative width comparable with . The smaller the value of , the less effective is reconnection. The peak value of the nonlinear growth rate is a measure of how fast the reconnection process occurs during the burst phase. In Fig. 2 we show that, contrary to the linear growth rate, the nonlinear growth rate dependence on becomes weaker and weaker for higher . This means that the nonlinear acceleration becomes more relevant the higher is the value of beta. This also means that reconnection process has a tendency to become independent on the microscopic scale size. Increasing the value of by two orders of magnitude, the nonlinear growth rate decreases by only one order of magnitude. For computational constraints, the maximum value of that we can investigate at present time is . For this value of , the nonlinear growth rate is , to be compared with the linear growth rate . We call such a regime explosive reconnection.
The evolution of the fields during the nonlinear acceleration phase has been investigated, and some markers have been recognized as characteristic of this particular phase. The analysis of their time evolution shows that they evolve with the same characteristic time scales of the magnetic island, with no higher frequency component. At the O-point, secondary vorticity structures are found by taking a cut along the line , . The vorticity structure during the whole linear evolution has definited sign on each side of the plane , consistent with a uniform inflating of the magnetic island. On the other hand, we find an inversion of sign inside the magnetic island starting with the end of the linear phase. In the case of , , the linear vorticity structure is doubled and reproduced inside the separatrix during the nonlinear acceleration phase, while the growth rate passes from the linear value to double of the linear value. At the X-point position, the formation of a Sweet-Parker like current sheet is found, corresponging to a structure with two Y-points. A successive recovering of the X-point structure is found at the beginning of the saturation phase. We find that current layers form during the linear phase along the whole separatrix, and their amplitude and steepness increase strongly during the nonlinear acceleration phase. The end of the acceleration phase is characterized by a saturation of the current spike amplitude, and by the formation of new current structures inside the magnetic island. Acceleration in the magnetic island growth are also found to be linked to the presence of turbulence around the X-point, having the possible effect of multiplying the number of X-points and consequently the rate of reconnection (similarly to results shown in Ref. ).
The ion finite Larmor radius (FLR) effects are studied by varying the ion to electron temperature ratio . We find that in the linear phase the FLR effects yields a slightly faster growth rate, consistently with the analytical predictions . The same is valid also in the nonlinear phase. Nevertheless, no qualitative difference is found in the case of and . The reason is that ion moments don’t enter the physics for -driven processes at low-beta regimes , namely when , and for small . On the other hand, FLR effects were found to be more relevant for 2D models in large regimes .
4 Energy consistency
The check of energy consistency is of crucial importance in reconnection simulations. In fact, reconnection processes can occur in numerical experiments also because of artificial - nonphysical - resistivity, such as numerical dissipation. In general, an artificial resistivity causes that the terms of Ohm’s law calculated separately don’t balance exactly in the subsequent time steps. Therefore, a numerical dissipation can act as a virtual resistivity and the result is a fictitious reconnection, whose growth rate can be in some cases higher than the physical one, depending on the chosen plasma parameters. In nonlinear simulations, smaller and smaller scales are created due to nonlinear coupling of macroscopic modes, and it’s important to avoid the accumulation of energy at the scales of the numerical grid. To this aim, hyperviscosity is added to our model equations, to filter small scales physics out, just above the grid size.
The contribution to the reconnection growth rate of subgrid dissipation, given by hyperviscosity and numerical dissipation, is calculated with a specific diagnostic in the GEM-code. This error growth rate is defined by:
Here is the reconnection growth rate, calculated as the time derivative of the total energy of the system: . The total energy is the volume integral of the sum of the energy densities , , and , being respectively the ExB energy, the thermal state variable energy, the flux variable energy and the magnetic energy . is the sum of the source growth rates: , and . is the sum of the absolute value of the damping rates, namely Landau damping of ion and electrons (where the electron Landau damping is the relevant one for our regime). The contribute of the numerical error is checked to be always negative during the whole run. In this case, the contribute of the subgrid dissipation to the reconnection rate is null.
In summary, we have presented the results of simulations on collisionless reconnection in a 3D Harris-pinch equilibrium, with flat density and temperature initial profiles. The gyrofluid code GEM has been used, and single helicity reconnecting modes have been considered in a small plasma regime. The limit of small , has been investigated. The acceleration phase during the nonlinear growth has been analyzed and the nonlinear growth rate peak scaling has been provided, with respect to the equilibrium plasma parameters. Numerical evidence of explosive reconnection has been found for . The underlying mechanism of our observation is thought to be linked with the presence of the high velocity fields. In fact, a collisionless regime allows the formation of strong velocity fields during the nonlinear acceleration phase. These velocity fields are retained to be responsible for the splitting of the X-point into two Y-points, and as a consequence, the explosive regime is entered. The presence of an acceleration in the growth of the magnetic island, could be important in explaining the fast growth of magnetic islands as seen in nature.
Simulations have been repeated with double or half amplitude of the initial magnetic island seed, and no difference in the value of the nonlinear growth rate has been found. Some markers of the acceleration phase have been pointed out. At the O-point, the formation of smaller vorticity structures inside the magnetic island has been found. At the X-point position, the X-point has been shown to split into two Y-points during the acceleration phase, and then to become again an X-point during the saturation phase. Along the whole separatrix, current spikes at scales narrower than the collisionless skin depths have been found. The importance of these effects in causing the reconnection acceleration is under investigation, and a detailed analysis will be published in a dedicated paper .
It is not quite understood whether the three-dimensionality of the configuration is a key ingredient for nonlinear acceleration. In fact, 2D reconnection is known to be governed by conservation laws: being and functions of the perpendicular variables x and y only, then the current and the vorticity are directed in the direction parallel to the guide field, and there is no interaction between neighbour helicities. On the other hand, in 3D simulations neighbour helicities can interact - for instance J can interact with and - and the nonlinear phase is predicted to have a qualitatively different evolution.
This work was carried out in a collaboration with the Max-Planck institute for solar system research at Katlenburg-Lindau, Germany. This Letter was written in Metz (France). The authors thank J. Buechner and E. Marsch for interesting discussions. Enlightening discussions with F. Pegoraro, F. Califano, D. Del Sarto, T. Ribeiro and F. Zonca are also gratefully acknowledged.
-  H.P. Furth, J. Killeen and M.N. Rosenbluth, Phys. Fluids 6, 459 (1963)
-  A.W. Edwards et al., Phys. Rev. Lett. 57, 210 (1986)
-  F. Porcelli, Phys. Rev. Lett. 66, 425 (1991)
-  B.D. Scott and F. Porcelli, Phys. Plasmas 11(12), 5468 (2004)
-  A. Biancalani and B.D. Scott, 37th EPS Conf. on Plasma Physics (Dublin, Ireland, 21-25 June 2010) http://ocs.ciemat.es/EPS2010PAP/pdf/P4.108.pdf
-  D. Grasso, E. Tassi and F.L. Waelbroeck, Phys. Plasmas 17, 082312 (2010)
-  M. Ottaviani and F. Porcelli, Phys. Rev. Lett. 71, 3802 (1993)
-  D. Grasso, F. Califano, F. Pegoraro and F. Porcelli, Phys. Rev. Lett. 86(22), 5051 (2001)
-  D. Borgogno, D. Grasso, F. Porcelli, F. Califano, F. Pegoraro and D. Farina, Phys. Plasmas 12, 032309 (2005)
-  A.Y. Aydemir, Phys. Fluids B Letter 4(11), 3469 (1992)
-  J.F. Drake and Y.C. Lee, Phys. Rev. Lett. 39(8), 453 (1977)
-  D. Grasso, F. Pegoraro, F. Porcelli, and F. Califano, Plasma Phys. Controlled Fusion 41, 1497 (1999)
-  A. Bhattacharjee, K. Germaschewski, and C. S. Ng, Phys. Plasmas 12, 042305 (2005).
-  B.D. Scott, Phys. Plasmas 12, 102307 (2005)
-  W. Dorland and G.W. Hammett, Phys. Fluids B 5(3), 812 (1993)
-  B.D. Scott, Phys. Plasmas 7(5), 1845 (2000)
-  N.F. Loureiro et al., Mon. Not. R. Astron. Soc. 399, L146 (2009)
-  A. Biancalani and B.D. Scott, to be submitted to Phys. Plasmas (2011)