Kinetic simulations of magnetostatic equilibria

Kinetic simulations of the lowest-order unstable mode of relativistic magnetostatic equilibria

Krzysztof Nalewajko1 2 3 , Jonathan Zrake2 , Yajie Yuan2 , William E. East2 , Roger D. Blandford2
1affiliation: Nicolaus Copernicus Astronomical Center, Bartycka 18, 00-716 Warsaw, Poland;
2affiliation: Kavli Institute for Particle Astrophysics and Cosmology, SLAC National Accelerator Laboratory, Stanford University, 2575 Sand Hill Road M/S 29, Menlo Park, CA 94025, USA
3affiliation: NASA Einstein Postdoctoral Fellow

We present the results of particle-in-cell numerical pair plasma simulations of relativistic 2D magnetostatic equilibria known as the “ABC” fields. In particular, we focus on the lowest-order unstable configuration consisting of two minima and two maxima of the magnetic vector potential. Breaking of the initial symmetry leads to exponential growth of the electric energy and to the formation of two current layers, which is consistent with the picture of “X-point collapse” first described by Syrovatskii. Magnetic reconnection within the layers heats a fraction of particles to very high energies. After the saturation of the linear instability, the current layers are disrupted and the system evolves chaotically, diffusing the particle energies in a stochastic second-order Fermi process leading to the formation of power-law energy distributions. The power-law slopes harden with the increasing mean magnetization, but they are significantly softer than those produced in simulations initiated from Harris-type layers. The maximum particle energy is proportional to the mean magnetization, which is attributed partly to the increase of the effective electric field and partly to the increase of the acceleration time scale. We describe in detail the evolving structure of the dynamical current layers, and report on the conservation of magnetic helicity. These results can be applied to highly magnetized astrophysical environments, where ideal plasma instabilities trigger rapid magnetic dissipation with efficient particle acceleration and flares of high-energy radiation.

Subject headings:
magnetic reconnection — acceleration of particles — plasmas

1. Introduction

Observations of high-energy photons from certain astrophysical sources — such as blazars, gamma-ray bursts, pulsars, etc. — often reveal dramatic energy dissipation and efficient particle acceleration on very short time scales. Examples include rapid gamma-ray variability of blazars (e.g., Aharonian et al., 2007; Aleksić et al., 2011; Saito et al., 2013; Hayashida et al., 2015), radio galaxies (e.g., Aleksić et al., 2014), gamma-ray bursts (e.g., Abdo et al., 2009), and the Crab pulsar wind nebula (Tavani et al., 2011; Abdo et al., 2011; Buehler et al., 2012). In many cases, the environment of these events is thought to be a highly magnetized collisionless plasma — AGN jets, pulsar wind nebulae, etc. This suggests scenarios involving efficient localized dissipation of magnetic energy allowing for rapid particle acceleration (e.g., Begelman et al., 2008; Giannios et al., 2009; Nalewajko et al., 2011; Uzdensky et al., 2011; Clausen-Brown & Lyutikov, 2012).

It is important to consider the astrophysically realistic situations leading to efficient magnetic dissipation. The foremost requirement is the localized reversal of the magnetic field lines. Using the Harris-type current layer as initial condition presumes a highly synchronized global reversal of the magnetic field polarity with a steady supply of magnetized plasma. Such initial condition includes a rather arbitrary structure on kinetic scales, the impact of which can probably be neglected only in sufficiently large simulations (Sironi & Spitkovsky, 2014; Guo et al., 2014; Werner et al., 2016). In the context of pulsar wind nebulae, global reversals can be readily realized in the equatorial striped wind (Coroniti, 1990; Lyubarsky & Kirk, 2001; Kirk & Skjæraasen, 2003), and their consequences for particle acceleration and high-energy emission are actively investigated (Sironi & Spitkovsky, 2011; Baty et al., 2013; Zrake, 2015). In the context of astrophysical jets, it has long been hypothesized that such global reversals can take place (e.g., Lovelace et al., 1997), however, a convincing demonstration of such scenario is yet to be made. The alternative is that magnetic dissipation is triggered by jet instabilities, in particular by the current-driven kink modes (Begelman, 1998). In such a case, interactions of large-scale magnetic structures (eddies) are expected to create many transient localized current sheets that facilitate magnetic dissipation (O’Neill et al., 2012; Mizuno et al., 2012).

Currently, most investigations focus on relativistic magnetic reconnection, and indeed, great progress has been made over the past few years, in large part due to increasingly powerful direct kinetic plasma simulations. When starting from a uniform Harris-type current layer, it has been convincingly demonstrated that relativistic reconnection leads to very efficient particle acceleration. The resulting energy distributions can be characterized as power-laws with broad cut-offs, and the slope depending mainly on the background plasma magnetization . In the limit of , the energy distributions become extremely hard, with (Sironi & Spitkovsky, 2014; Guo et al., 2014; Werner et al., 2016). Such hard energy distributions mean that the maximum particle energy is limited by the initial average magnetic energy per particle. Indeed, for sufficiently large systems, the characteristic cut-off energy scales like (Werner et al., 2016).

However, the dissipation efficiency of relativistic reconnection is limited by reconnection rates being of order (e.g., Liu et al., 2015). Our group proposed an idea of magnetoluminescence, a generic dissipation process that allows for rapid and efficient conversion of electromagnetic energy into radiation (Blandford et al., 2014, 2015). With this end, we started to explore a range of novel magnetostatic equilibria.

Recently, some of us investigated a class of unstable magnetostatic equilibria known as the “Arnold-Beltrami-Childress” (ABC) fields, using relativistic magnetohydrodynamics (RMHD) and force-free electrodynamics (FF) codes (East et al., 2015). Here, we present the first results of the subsequent investigation of these equilibria with the kinetic particle-in-cell (PIC) code Zeltron. The simulations presented here are two-dimensional, with pair plasma, and they focus on the lowest-order unstable mode. Nevertheless, we obtained a very rich physical picture including the evolution of the sheared current layers, regular and stochastic modes of particle acceleration, and rapid dissipation of magnetic energy.

In Section 2, we define the initial configuration and its implementation in the Zeltron code. In Section 3, we present the results including the evolution of total energy components, conservation of magnetic helicity, the evolving structure of the current layers, evolution of the particle energy distribution, and analysis of individual energetic particles. Our results are discussed in Section 4, and our conclusions are presented in Section 5.

2. Simulation setup

We performed two-dimensional PIC numerical simulations using the Zeltron code111 (Cerutti et al., 2013). The core algorithms include a Finite-Difference-Time-Domain (FDTD) method for advancing the electromagnetic fields on staggered grids (Yee, 1966), the “Boris push” for advancing the particles, a smoothing filtering of the electric fields, and a Poisson solver for matching the electric fields with the charge density. The radiation reaction implemented in this code was not used in these simulations.

Our simulations proceed from smooth magnetostatic equilibria. They are kinetic generalizations of the MHD states whose gas pressure is uniform and Lorentz force vanishes. We focus on force-free equilibria satisfying the Beltrami condition , where is a constant. In the present study we use the following magnetic field


where for linear domain size . This state has a wavelength that is smaller by a factor than the domain scale, and resembles the conventional ABC magnetic fields we used in previous studies, but has been rotated by . Translational symmetry of this state along is enforced throughout the evolution by the “2.5D” simulation scheme. This state contains two minima and two maxima of the out-of-plane vector potential . Those extrema are centered on helical flux tubes aligned with the -axis, which are linearly unstable to pairwise coalescence instability as we reported in East et al. (2015). As such coalescence moves plasma toward a lower energy (longer wavelength) magnetic configuration, the initial state must possess a non-zero “magnetostatic free energy”. That is, by our definition, the magnetic energy that could be removed while the frozen-in MHD condition — , hence — is maintained over all but infinitesimal volumes. Such evolution leads to minimization of the magnetic energy as constrained by the global magnetic helicity invariant in the sense of Taylor (1974). For general states having a wavelength , the theoretical free energy fraction is in general , which has the value for our setup. In other words, of the magnetic energy of our initial condition could be used to energize particles if magnetic helicity is conserved.

The current density of the setup in Equation (2) is realized by endowing particles with a locally anisotropic momentum distribution, which we factorize into independent energy and angular parts as


where is the dimensionless particle momentum, is the dimensionless particle energy (Lorentz factor), is the dimensionless particle velocity. Function is the Maxwell-Jüttner energy distribution for dimensionless temperature ; , where is the polar angle measured from the unit vector aligned with the required local current density vector. Parameter is the dipole moment of the local angular distribution of paricle momenta, which determines the local average drifting speed , where is the average speed that is a function of .

s07L400 400 1/4 3.4 0.7 0.33 0.77 0.03 0.11 3.3 105
s07L800 800 1/8 3.4 0.7 0.36 0.78 0.04 0.12 3.7 95
s14L400 400 1/2 6.8 1.4 0.23 0.86 0.07 0.23 3.1 175
s14L800 800 1/4 6.8 1.4 0.25 0.80 0.06 0.21 2.9 180
s14L1600 1600 1/8 6.8 1.4 0.27 0.78 0.06 0.20 3.1 255
s27L800 800 1/2 13.6 2.7 0.20 0.84 0.09 0.30 2.9 380
s27L1600 1600 1/4 13.6 2.7 0.21 0.81 0.10 0.38 2.7 480
s55L1600 1600 1/2 27.2 5.5 0.17 0.83 0.17 0.59 2.5 790
s55L3200 3200 1/4 27.2 5.5 0.18
Table 1Parameters of the simulation runs. Two characteristic magnetization values are reported: and . Parameter is the characteristic constant of the dipole moment of particle momentum distribution. Parameter is the e-folding growth rate of the total electric energy. Parameter describes the global efficiency of magnetic energy conversion, where is the measured conversion efficiency, and is the theoretical expectation for relaxation to the isotropic Taylor state. Parameters and describe the number and energy fractions contained in the high-energy particle distribution tail at . Parameter is the index of the power-law energy distribution tail estimated at . Parameter is the maximum particle energy measured at the level of the normalized energy distribution at .

For pair plasma, the total current density module is , where is the number density of both electrons and positrons. Combining this with previous relations, we obtain an expression for equilibrium value of the pair number density:


where is a constant equal to a characteristic value of the dipole moment (note that ). The upper limit imposed on the dipole moment introduces a lower limit on the particle number density, and hence an upper limit on the magnetization of the simulated system. The characteristic value of the “cold magnetization” can be expressed as:


where is the nominal gyroradius. Hence, for the maximum value of the dipole moment, the magnetization scales linearly with the system size . For isotropic distributions, the cold magnetization is related to the plasma skin depth . We also define the “hot magnetization” as , where is the specific enthalpy. However, one should note that even the initial magnetization is not uniform due to the non-uniformity of .

We have performed several simulations for different values of and , as reported in Table 1. Common parameter values are and . The numerical resolution is , and the total number of particles per cell is 128.

3. Results

3.1. Snapshot maps

Figure 1 shows the distribution of magnetic field component , particle number density , and average particle Lorentz factor at several time steps for Run s55L1600. The initial condition includes two maxima (red) and two minima (blue) of , which correspond to the extrema of the magnetic vector potential . By , the initial symmetry of the setup is broken, as can be seen on the particle number density map. Two current layers form, one in the center of the simulation domain, and one in the corners. By , the current layers grow in length, they accelerate particles to high energies, and they undergo a minor tearing instability. This process can be characterized as the “X-point collapse” first investigated by Syrovatskii (1966). The sharp current layers seen at are disrupted, and the subsequent evolution of the magnetic domains proceeds through a series of oscillations forming complex, chaotic substructures. However, in the final stage, we observe a clear separation between a dense cold plasma located along the magnetic domain boundaries, and a dilute hot plasma partially filling the domain interiors.

Figure 1.— Snapshots from the Run s55L1600. Time is normalized to the light-crossing time scale . Top panels: magnetic field component ; middle panels: number density of electrons and positrons; bottom panels: average Lorentz factor of electrons and positrons.

3.2. Evolution of total energy

Figure 2 shows the time evolution of the main components of the total energy: magnetic, kinetic (including thermal), and electric. The total energy in our simulations is conserved at the level better than , improving with increasing magnetization. All simulations (except for s55L3200) are run for at least 10 light-crossing time scales, at which point the evolution of the total energy components is largely complete. Parameters of the initial configuration determine the characteristic magnetization value, and hence the initial energy shared between the magnetic field and particles.

Figure 2.— Evolution of the total energy components compared for different runs. All values are normalized to the initial total energy. Line colors indicate the mean value of the hot magnetization: (gray), (red), (green), and (blue). Line types indicate the simulation domain size: (dotted), (dashed), and (solid).

The global efficiency of magnetic dissipation is calculated as , where , is the total magnetic energy at , and is the total magnetic energy at . As we report in Table 1, for all simulation runs. This indicates that magnetic relaxation toward the Taylor minimum is effectively complete by that time. We note here that in 2.5D, the energy of the most relaxed state generally exceeds the Taylor minimum energy due to constraints on the change of magnetic topology that were not considered by Taylor, who intended to characterize fully 3D evolution. Nevertheless, Zrake & East (2016) found the additional constraints became significant only when . Since we have , nearly complete relaxation to the Taylor minimum energy state is to be expected in the cases presented here. We also estimate the peak -folding time scale of the total magnetic energy as .

The initial electric energy, at the level of of the total energy, is determined by the residual charge density due to the finite number of particles per cell. At , the electric energy begins to grow exponentially. By , the growth of the electric energy saturates at the level of of the total energy. At about the same time, the magnetic energy decreases sharply, and the kinetic energy begins to increase. There is no simple energy transformation between magnetic and kinetic forms. Instead, the electric energy begins to oscillate, and this oscillation is damping over many light-crossing time scales. These oscillations are reflected more clearly in the evolution of the magnetic energy component, rather than of the kinetic component.

We have also calculated the evolution of non-ideal electric field parallel to the magnetic field . Figure 2 shows that the growth of the non-ideal electric energy is negligible when compared with the growth of the total electric energy. For low magnetizations it is difficult to distinguish the total non-ideal component from the noise. In fact, non-ideal regions form distinct compact spatial structures, while most of the domain volume is consistent with ideal MHD at all times.

We have measured the growth rates of the total electric energy during the linear instability phase as -folding time . Figure 3 shows the dependence of on the mean magnetization . For each value of , the growth rate increases systematically with the domain size, which suggests that we do not achieve a complete convergence. Considering only the growth rates measured for the largest simulation for each magnetization value, we fitted them with the following function , where is the Alfvén speed. For the scaling parameters, we found and .

Figure 3.— Linear instability growth rate, defined as e-folding time scale of the total electric energy normalized to , as function of the mean value of the hot magnetization. The point size indicates the simulation domain size . The blue arrow indicates the asymptotic value measured in FF simulations (East et al., 2015).

3.3. Magnetic helicity conservation

Zrake & East (2016) studied the decay of magnetic turbulence in the force-free limit in both 2D and 3D. They investigated the conservation of magnetic helicity , which is formally broken only by the presence of regions with non-ideal electric field, as (e.g. Brandenburg et al., 2015). In both 2D and in 3D, they reported that total magnetic helicity was indeed conserved to a high precision, which improved with resolution. They also found that in 2D, could be decomposed into a continuous distribution , each value of which is a separate invariant. Here we verify these findings in 2D PIC simulations for the first time.

In Figure 4, we show that total magnetic helicity is conserved at the level better than , improving with the decreasing particle momentum dipole moment . As we will demonstrate later on, smaller dipole moment allows for higher particle density and better screening of the electric fields. Consequently, the volume fraction of the non-ideal field regions decreases, and this may explain better conservation of the magnetic helicity, in line with the results of Zrake & East (2016) in the force-free limit. In Figure 4, we also show that while is well conserved for all values of , magnetic energy is dissipated during the linear instability phase for all values of .

Figure 4.— First panel (from the left): the accuracy of total helicity conservation compared for several simulations. Line colors and types are the same as in Figure 2. Second panel: time evolution of the helicity profile for Run s14L1600. Third panel: time evolution of the magnetic energy profile for Run s14L1600. For the two right panels, the color scale indicates simulation time normalized to the light-crossing time scale .

3.4. Structure of the current layers

Two current layers are formed during the late stages of the instability, at . They appear as thin structures characterized by high density and high average particle energy, they grow in length and undergo tearing instability until they are disrupted. These are also regions where non-ideal parallel electric fields are present, with , and hence they enable efficient particle acceleration.

Since the current layers are not present in the initial configuration, it is interesting to characterize their perpendicular structure and evolution. To this end, we first determined the exact location and orientation of the layers by fitting a two-dimensional Gaussian model to the spatial distribution of the average particle energy. The time range when this can be done robustly is limited — for the current layers are not detectable, and for the current layers bend and eventually disappear. We introduce a local coordinate system with parallel/perpendicular vectors measured along/across the major axis of the current layer. As the location and orientation of the current layers changes slowly, this is repeated at every timestep of interest. Next, we calculated perpendicular profiles of various parameters across the current layer centroid. These profiles are shown in Figure 5 at 4 different times to illustrate a fairly complex evolution of the current layer.

Figure 5.— Perpendicular profiles measured across the centroid of the current layer for Run s27L800 at several times: (dashed lines), and (solid lines). The panels show: (1) number density of electrons and positrons normalized to the initial value, (2) average energy of electrons and positrons, (3) three components of the current density normalized to the initial number density, (4) three components of the mean velocity of positrons, (5) three components of the magnetic field, (6) three components of the electric field and the value. Subscript ‘’ means the component parallel/perpendicular to the major axis of the current layer.

The density profile develops a narrow spike, and the average particle energy profile shows a similar but somewhat broader structure. On the other hand, the current density component shows a more complex profile with a narrow core and broader wings. This appears to be a combination of the narrow density profile with the broad profile of the -component of the drifting velocity. We also observe a small component of alternating parallel current on the scale comparable with the velocity structure. This current is of opposite sign to the background in-plane current that was already present in the initial configuration. Evolution of magnetic field involves a gradual steepening of the parallel component , and a systematic increase of the out-of-plane component . The electric field increases systematically in all components. It is interesting that the component remains very uniform in the process. When we subtracted the ideal field contribution , we found that . The non-ideal field component peaks in the middle of the layer, and has thickness scale comparable to the structure.

We fitted the perpendicular profiles of particle number density , average particle energy , and with Gaussian models with uniform background. In Figure 6, we show the evolution of the amplitudes and perpendicular dispersions of the three parameters. This reveals the complex structure of the current layer, with the three parameters characterized by different thickness scales. The general trend in time is an increase in the amplitudes of all three parameters, the thinning of the density and average energy profiles, and the broadening of the profile. The amplitudes of the average particle energy grow exponentially in time as long as we can measure their profile, with the growth rate increasing systematically with the magnetization . On the other hand, the amplitudes of appear to grow linearly in time, and the growth rate is similar for . The width scales of the density and temperature profiles converge to a value , roughly independent of . The difference between the width scales measured for the density and temperature profiles is probably not significant. However, the width scale of the clearly depends on and increases with . We suggest that the width scale of the density and temperature profiles corresponds roughly to the plasma skin depth , and that the width scale of the region corresponds to the typical gyroradius of particles heated in the layer .

Figure 6.— Evolution of the transverse structure of the current layer in the linear phase for runs: s14L1600 (red lines), s27L1600 (solid green lines), s27L800 (dashed green lines), and s55L1600 (blue lines). We have calculated the transverse profiles of particle number density (left panels), average particle energy (middle panels) and non-ideal electric field (right panels) across the centroid of the current layer. These profiles were fitted with a model consisting of a Gaussian and a constant background. The upper panels show the evolution of the amplitude of the Gaussian model, and the lower panels show the evolution of the width scale in units of . Occasional sharp spikes are due to the passage of plasmoids.

We also calculated the volume-weighted distribution functions for key scalar quantities and . Figure 7 shows the 99 percentile values as functions of simulation time. The values of consistently reach the level of , independently of the setup parameters, at about the time of linear instability saturation. The situation with the values is less clear, with a number of sharp peaks observed at different times. The smallest values are recorded for Run s14L1600, which is qualitatively consistent with the results shown in Figure 6. For other runs, the peak values are in the range .

Figure 7.— Evolution of the 99 percentile volume-weighted values of (left panel) and (right panel) compared for several runs. Line colors and types are the same as in Figure 2.

3.5. Particle acceleration - total spectra

In Figure 8, we show the momentum distributions of electrons and positrons as functions of simulation time. The distributions are normalized to the peak of the initial distribution, which in all cases is the Maxwell-Jüttner distribution for . The initial distribution begins to evolve at (blue) in a period of rapid and regular (linear) energization of a relatively small fraction of all particles. The uniformity of the individual distributions on the log-log plot suggests an exponential growth of both the maximum particle energy and the number fraction of participating particles. For the case of Run s55L1600, this acceleration phase is efficient enough to form a second peak of the momentum distribution. The linear acceleration phase eventually ceases, and the subsequent non-linear evolution of the momentum distribution leads to a gradual formation of a power-law. The “valley” between the initial Maxwellian component and the high-energy bump is promptly leveled. The maximum particle energy increases slowly but systematically in several steps corresponding to the global oscillation of the merged magnetic domains.

Figure 8.— Evolution of the momentum distributions of electrons and positrons for selected runs. The color scale indicates simulation time normalized to the light-crossing time scale . The distributions are sampled at linearly uniform time intervals, however, the intervals are different for each simulation.

We have roughly estimated the index of the power-law tail of the distribution for each simulation at . The values are reported for selected runs in Figure 8. There is a clear trend of increasing with , from for to for . We have not attempted to measure the power-law index rigorously as a function of the simulation time, as accurate automated decomposition of the distribution function requires very complicated modeling (Werner et al., 2016).

Instead, we investigate what is the fraction of all particles contained in the high-energy tail by subtracting the low-energy Maxwellian component. In Figure 9, we plot as functions of simulation time: the fraction of the particle number , the fraction of the particle energy , and the maximum particle energy , defined formally at the level of for the normalized distributions. During the linear acceleration phase, both , and grow roughly exponentially with time. In the non-linear phase, the growth becomes much slower, although it does not fully saturate before , especially for higher values of .

Figure 9.— Upper left panel: time evolution of the number fraction contained in the non-thermal high-energy tail of the electron distribution. Upper right panel: time evolution of the energy fraction contained in the high-energy tail. Lower left panel: time evolution of the maximum particle Lorentz factor measured at the level of for normalized electron distributions. Line colors and types are the same as in Figure 2.

The values of , and at for all simulations are reported in Table 1 and plotted in Figure 10. The number fraction ranges from for to for . The corresponding values of the energy fraction are and . Both and appear to scale with the magnetization like , with the exception of the highest value, although simulations with even greater are needed to verify this scaling. In any case, this trend could not continue to very high values of , as . The ratio , corresponding to the ratio of average particle energies , shows a slight increase from for to for . The maximum particle energy ranges from for to for , and is consistent with a linear dependence on .

Figure 10.— Left panel: dependence of the number (blue squares) and energy (red circles) fractions contained in the high-energy tail of the electron distribution, averaged for , on the mean hot magnetization value. Right panel: dependence of the maximum Lorentz factor measured at the level of for normalized electron distributions, averaged for . The point size indicates the simulation domain size .

3.6. Individual particle histories

At the beginning of each simulation, we randomly selected individual particles in order to record their detailed histories. We will now use this information in order to better understand particle acceleration in the linear stage of the simulation. In Figure 11, we show a sample of particle energy histories for particles that obtain the highest energies at . For every particle, this reveals a sharp acceleration episode during the saturation of the linear instability. Now we ask the question, whether the difference in the final energy values obtained by these most energetic particles for each simulation is due to stronger effective value of the parallel electric field, or rather due to longer acceleration time scale. To this end, we have analyzed the energy histories , identifying episodes of the largest monotonic energy gain . We applied a mild smoothing to , and we excluded episodes with . Then, we calculated the mean electric field accelerating the particle from , where and .

Figure 11.— Energy histories for individual particles selected for their high energy at . Each diagram shows roughly most energetic particles.

In Figure 12, we show the distribution of vs. for three simulations with different magnetization values, indicating also the final particle energy at . For each simulation, there is a considerable spread in the values of and calculated for individual energetic particles. While and determine the energy gain for the main acceleration episode, there is also a substantial scatter between and . The mean values of are for , respectively, and the corresponding mean values of are . With this, we can account for a factor of difference between the maximum particle energy between and , while according to Table 1 the difference between values for these simulations is by factor . The difference appears to be both due to increase in typical electric field strength between and , as well as due to longer acceleration time scale between and .

Figure 12.— Duration vs. mean electric field calculated for the main acceleration episodes for representative samples of energetic individually tracked particles. The area of the markers is proportional to the final particle energy calculated for . Each diagram shows roughly most energetic particles.

4. Discussion

Simulating force-free magnetostatic equilibria offers a new approach for studying relativistic magnetic dissipation, an alternative to simulating the Harris-type current layers. A major conceptual advantage of this approach is that there are no current layers or any structures on the kinetic scale in the initial plasma configuration. The formation of current layers results from the saturation of an ideal plasma instability, and this rapidly creates conditions for efficient particle acceleration. The most energetic particles follow the evolving current layers until they disintegrate on the dynamical time scale, and afterwards they interact chaotically with the slowly decaying turbulence.

The initial configuration consists of two positive and two negative domains of magnetic vector potential . These symmetrically interlocked domains attract each other and hence they remain in an unstable equilibrium until the particle distribution noise breaks the positional symmetry. This triggers an exponential increase in the bulk plasma velocities , and hence an exponential increase in the ideal electric fields . Figure 2 shows that this instability is not seen in the time evolution of the global non-ideal electric energy. The e-folding growth rate of the total electric energy is measured to scale with the mean hot magnetization roughly like . In the limit of , this is in excellent agreement with the value (see Figure 3) measured in the FF simulations (East et al., 2015).

The origin of the current layers can be traced to the converging motion of two magnetic domains resulting in a pile-up of particles at the magnetic null points. The exponential increase in the plasma velocity creates an exponentially increasing flux of particles forced into the current layers. This driven reconnection is in contrast to the undriven reconnection resulting from relaxed Harris-type current layers, where magnetic diffusion regions are characterized by low plasma density. The current layers in our simulations systematically grow in length, they become hotter and thinner, and they are subject to only mild tearing mode instability. We observe that around the saturation moment of the ideal instability the inflow of plasma into the current layer is interrupted, and the current layer is rapidly disintegrated by internal motions of the newly formed magnetic domain encompassing the layer.

Throughout the evolution of the current layers, we observe a systematic increase and smooth structure across the layer of the total electric field. This field originated outside the current layers from the motions of ideal plasma as . In the ideal MHD limit, such field should vanish in the middle of the current layer, and we would strictly expect that . The ideal MHD condition is obviously violated inside the current layer, and is dominated by . The smooth structure of the component suggests that it is advected into the layer from outside, rather than created locally in a separate process. The thickness scale of the region is significantly larger than the thickness scale of the particle density pile-up, moreover, it is systematically increasing with the simulation time. This additional thickness scale is roughly of the order of the gyroradius corresponding to the mean particle energy inside the layer, and it can be seen imprinted into the profile of the current density. The Vlasov momentum equation dictates that for each particle species the following must be satisfied: , where is the momentum density, and is the pressure tensor. We have checked that in general the pressure gradients do not support the , and hence the particle momentum density has to increase, i.e., particles are collectively heated in the regions.

The particle energy distributions evolve from the initial relativistic Maxwell-Jüttner to form a high-energy power-law tail on the time scale of in two stages. In the first, linear stage, particles are heated by regular electric fields in the linearly growing current layers. These particles produce a distinct high-energy bump in the energy distribution, with the number of heated particles and their maximum energy increasing exponentially with time. In the second, non-linear stage, particles energized in the linear stage become detached from the magnetic potential, and they interact chaotically with secondary current layers arising in the slowly damped turbulent plasma motions, diffusing in the energy space. The power-law distributions arise in the non-linear stage, however, their indices appear to be determined already in the linear stage. They can be robustly predicted by the ratio of the two distribution peaks at the saturation moment of the linear instability, or in other words, by the energy fraction of particles interacting with the current layers in the linear stage.

We showed that the acceleration process in the linear stage is very regular with both the maximum particle energy and the number fraction of the high-energy particles increasing exponentially with time. On the other hand, we showed that individual particles are accelerated linearly in time, as if subjected to a roughly constant value of . The exponential growth of the number fraction of high-energy particles can be related to the exponential growth of the flux of particles joining the current layer. The exponential growth of the maximum particle energy in the linear stage is consistent with the growth of the peak value of the non-ideal electric field component . In the non-linear stage, the growth of is close to linear, which would correspond to a constant value of . The final values of the non-thermal particle fractions appear to scale with the mean magnetization as . These scalings are expected to saturate below unity for , and the index may result from a transition between non-relativistic and ultra-relativistic regimes. On the other hand, the maximum particle energy scales roughly like , as has been found for sufficiently large simulations of relativistic reconnection initiated from the Harris-type layers (Werner et al., 2016). Figure 12 shows that the final particle energy is determined both by the non-ideal electric field strength and by the time spent by the particle in the acceleration region.

The emergence of the power-law component can be attributed to the second-order Fermi process. This process can be understood in the following way — (1) high-energy particles emerging from the linear current layers are characterized by large gyroradii, and they follow on chaotic orbits, independently of the small-scale magnetic structures; (2) oscillations of the merged magnetic structures create multiple short-lived current layers with different orientations of the non-ideal electric field; (3) high-energy particles are more likely to pass across these new acceleration regions, however they have comparable chances to be accelerated or decelerated, depending on the sign of . We do not find any clear signatures of the first-order Fermi process in the non-linear acceleration phase (cf. Hoshino, 2012; Beresnyak & Li, 2016).

The power-law indices for the highest mean cold magnetization probed are , which is significantly higher than the values reported for Harris-layer simulations with comparable upstream values of (Sironi & Spitkovsky, 2014; Guo et al., 2014; Werner et al., 2016). This may actually be useful for astrophysical applications, where there are more cases for high particle distribution indices ; in this sense the Harris-type reconnection in the ultra-relativistic limit appears to be too efficient a particle accelerator. The lower acceleration efficiency can be explained by the lower fraction of magnetic free energy — in our case vs. in the case of Harris-type layers. The fraction of free magnetic energy in our 2D simulations is only slightly lower than the theoretical prediction of corresponding to the isotropic Taylor state (see Section 2). In the case of periodic Harris-type setups, the final configuration departures further from the Taylor state, and hence an analytical prediction of is more difficult. Finally, we note that the free energy associated with ABC magnetic fields can be increased indefinitely by considering higher-order harmonics. However, in such case the topological constraints limit the dissipation efficiency in 2D, and hence investigating highly efficient magnetic dissipation requires large 3D simulations (Zrake & East, 2016).

A certain disadvantage of the “ABC” magnetostatic equilibria is that their mean magnetization is numerically limited by the spatial resolution , at least under our conservative requirement that the kinetic scale should always be resolved. This limitation arises from the uniform distribution of magnetic field gradients that have to be supported by volume-filling currents. Effectively, our study is confined to one order of magnitude in plasma magnetization. This is in contrast to the Harris-type configuration, where magnetic field gradients are confined to narrow layers, and hence a wide range of upstream magnetizations, spanning almost three orders of magnitude, can be studied at fixed numerical resolution (Guo et al., 2014; Werner et al., 2016).

Our simulations correspond effectively to guide-field reconnection, and given the Beltrami condition we have no way to vary the level of the guide field. Even though vanishes initially along the separatrices, it is readily advected, roughly incompressibly, into the current layers, so that at all times everywhere in the layer (see Figure 7). This is a potential concern if we would like to apply these results to the gamma-ray flares in the Crab Nebula, as the prospect for breaking the synchrotron photon energy limit in the Harris-layer reconnection scenario depends on the guide field (Cerutti et al., 2013, 2014). However, what we need in order to break the synchrotron limit is , and these components depend on the exact trajectories of energetic particles. Simulations of force-free magnetostatic equilibria taking into account the synchrotron radiation reaction are required to address this problem, and will be presented in a future study (Y. Yuan et al., in preparation).

5. Conclusions

We have investigated magnetic dissipation and particle acceleration in the simplest case of relativistic unstable 2D “ABC” magnetostatic equilibria by means of fully kinetic PIC simulations with pair plasma. An ideal instability leads to the formation of dynamical current layers with complex internal structure, where a fraction of particles is heated collectively by non-ideal electric fields. Saturation of the instability leads to a slowly decaying turbulence which diffuses the energetic particles in the energy space, effectively forming power-law energy distributions. The power-law slopes are significantly softer () when compared with the case of ultra-relativistic reconnection initiated from the Harris-type layers (). This kind of magnetic dissipation corresponds to driven relativistic reconnection with a substantial guide field. However, since most of the dissipation and particle acceleration proceed over a single light-crossing time scale, this is an attractive scenario for the production of rapid high-energy radiation flares observed in blazars, pulsar wind nebulae, GRBs, etc.

We thank Maxim Lyutikov and Lorenzo Sironi for stimulating discussions. K.N. was supported by NASA through Einstein Postdoctoral Fellowship grant number PF3-140130 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. We acknowledge the use of computational resources obtained from XSEDE (Stampede), University of Colorado Research Computing (Janus), and KIPAC/SLAC (Bullet).


  • Abdo et al. (2009) Abdo, A. A., et al., 2009, Nature, 462, 331
  • Abdo et al. (2011) Abdo, A. A., Ackermann, M., Ajello, M., et al., 2011, Science, 331, 739
  • Aharonian et al. (2007) Aharonian, F., et al., 2007, ApJ, 664, L71
  • Aleksić et al. (2011) Aleksić, J., et al., 2011, ApJ, 730, L8
  • Aleksić et al. (2014) Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014, Science, 346, 1080
  • Baty et al. (2013) Baty, H., Petri, J., & Zenitani, S. 2013, MNRAS, 436, L20
  • Begelman (1998) Begelman, M. C., 1998, ApJ, 493, 291
  • Begelman et al. (2008) Begelman, M. C., Fabian, A. C., Rees, M. J., 2008, MNRAS, 384, L19
  • Beresnyak & Li (2016) Beresnyak, A., & Li, H. 2016, ApJ, 819, 90
  • Blandford et al. (2014) Blandford, R., Simeon, P., & Yuan, Y. 2014, NuPhS, 256, 9
  • Blandford et al. (2015) Blandford, R., East, W., Nalewajko, K., Yuan, Y., & Zrake, J. 2015, arXiv:1511.07515
  • Brandenburg et al. (2015) Brandenburg, A., Kahniashvili, T., & Tevzadze, A. G. 2015, PhRvL, 114, 075001
  • Buehler et al. (2012) Buehler, R., Scargle, J. D., Blandford, R. D., et al. 2012, ApJ, 749, 26
  • Cerutti et al. (2013) Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147
  • Cerutti et al. (2014) Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2014, ApJ, 782, 104
  • Clausen-Brown & Lyutikov (2012) Clausen-Brown, E., & Lyutikov, M. 2012, MNRAS, 426, 1374
  • Coroniti (1990) Coroniti, F. V. 1990, ApJ, 349, 538
  • East et al. (2015) East, W. E., Zrake, J., Yuan, Y., & Blandford, R. D. 2015, PhRvL, 115, 095002
  • Giannios et al. (2009) Giannios, D., Uzdensky, D. A., & Begelman, M. C., 2009, MNRAS, 395, L29
  • Guo et al. (2014) Guo, F., Li, H., Daughton, W., & Liu, Y.-H. 2014, PhRvL, 113, 155005
  • Hayashida et al. (2015) Hayashida, M., Nalewajko, K., Madejski, G. M., et al. 2015, ApJ, 807, 79
  • Hoshino (2012) Hoshino, M. 2012, PhRvL, 108, 135003
  • Kirk & Skjæraasen (2003) Kirk, J. G., & Skjæraasen, O. 2003, ApJ, 591, 366
  • Liu et al. (2015) Liu, Y.-H., Guo, F., Daughton, W., Li, H., & Hesse, M. 2015, PhRvL, 114, 095002
  • Lovelace et al. (1997) Lovelace, R. V. E., Newman, W. I., & Romanova, M. M., 1997, ApJ, 484, 628
  • Lyubarsky & Kirk (2001) Lyubarsky, Y., & Kirk, J. G. 2001, ApJ, 547, 437
  • Mizuno et al. (2012) Mizuno, Y., Lyubarsky, Y., Nishikawa, K.-I., & Hardee, P. E. 2012, ApJ, 757, 16
  • Nalewajko et al. (2011) Nalewajko, K., Giannios, D., Begelman, M. C., Uzdensky, D. A., & Sikora, M., 2011, MNRAS, 413, 333
  • O’Neill et al. (2012) O’Neill, S. M., Beckwith, K., & Begelman, M. C., 2012, MNRAS, 422, 1436
  • Saito et al. (2013) Saito, S., Stawarz, Ł., Tanaka, Y. T., et al. 2013, ApJ, 766, L11
  • Sironi & Spitkovsky (2011) Sironi, L., & Spitkovsky, A. 2011, ApJ, 741, 39
  • Sironi & Spitkovsky (2014) Sironi, L., & Spitkovsky, A. 2014, ApJ, 783, L21
  • Syrovatskii (1966) Syrovatskii, S. I. 1966, Soviet Ast., 10, 270
  • Tavani et al. (2011) Tavani, M., Bulgarelli, A., Vittorini, V., et al. 2011, Science, 331, 736
  • Taylor (1974) Taylor, J. B. 1974, Physical Review Letters, 33, 1139
  • Uzdensky et al. (2011) Uzdensky, D. A., Cerutti, B., & Begelman, M. C. 2011, ApJ, 737, L40
  • Werner et al. (2016) Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8
  • Yee (1966) Yee, K. 1966, IEEE Transactions on Antennas and Propagation, 14, 302
  • Zrake & East (2016) Zrake, J., & East, W. E. 2016, ApJ, 817, 89
  • Zrake (2015) Zrake, J. 2015, arXiv:1512.05426
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description