Kinetic study of radiation-reaction-limited particle acceleration during the relaxation of unstable force-free equilibria
Many powerful and variable gamma-ray sources, including pulsar wind nebulae, active galactic nuclei and gamma-ray bursts, seem capable of accelerating particles to gamma-ray emitting energies efficiently over very short time scales. These are likely due to rapid dissipation of electromagnetic energy in a highly magnetized, relativistic plasma. In order to understand the generic features of such processes, we have investigated simple models based on relaxation of unstable force-free magnetostatic equilibria. In this work, we make the connection between the corresponding plasma dynamics and the expected radiation signal, using 2D particle-in-cell simulations that self-consistently include synchrotron radiation reaction. We focus on the lowest order unstable force-free equilibrium in a 2D periodic box. We find that rapid variability, with modest apparent radiation efficiency as perceived by a fixed observer, can be produced during the evolution of the instability. The “flares” are accompanied by an increased polarization degree in the high energy band, with rapid variation in the polarization angle. Furthermore, the separation between the acceleration sites and the synchrotron radiation sites for the highest energy particles facilitates acceleration beyond the synchrotron radiation reaction limit. We also discuss the dynamical consequences of radiation reaction, and some astrophysical applications of this model. Our current simulations with numerically tractable parameters are not yet able to reproduce the most dramatic gamma-ray flares, e.g., from Crab Nebula. Higher magnetization studies are promising and will be carried out in the future.
Subject headings:magnetic reconnection — acceleration of particles — radiation mechanisms: non-thermal — plasmas
Many powerful gamma-ray sources are found to produce dramatic flares that suggest rather high radiation efficiency. One example is the Crab Nebula (e.g., Abdo et al., 2011; Tavani et al., 2011; Buehler & Blandford, 2014), which flares in the energy range 100 MeV1 GeV roughly once per year. The biggest flare reached an isotropic luminosity of the pulsar’s spin down power, and the flux doubling time scale can be as short as a few hours (Buehler et al., 2012). This suggests an isotropic fluence erg, taking erg s and hr. If we define the radiation efficiency as the ratio between and the magnetic energy contained in the volume , we would get , assuming an average magnetic field of G. The radiation mechanism is generally believed to be synchrotron, but the spectral peak goes beyond the classical synchrotron radiation reaction limit MeV (e.g., Landau & Lifshitz, 1975). This, together with the fast variability, large total energy and non-detection in other wavebands (e.g., Weisskopf et al., 2013; Rudy et al., 2015), presents great challenges to existing particle acceleration theories. Rapid variability is also notably detected in many Active Galactic Nuclei (AGN). For example, the flat spectrum radio quasars (FSRQ) 3C 273 (e.g., Rani et al., 2013) and 3C 279 (Hayashida et al., 2015) produce GeV flares on hourly time scales; the radio galaxy IC310 Aleksić et al. (2014), the FSRQ PKS 1222+21 (Aleksić et al., 2011) and several BL Lac objects, e.g. PKS 2155-304 (e.g., Aharonian et al., 2007) and Markarian 501 (Albert et al., 2007), are observed to flare in TeV band with variability time scales a few minutes. More strikingly, quite a few flares from blazars are found to be accompanied by large polarization angle swings in the optical band (Abdo et al., 2010; Marscher et al., 2008; Blinov et al., 2015; Kiehlmann et al., 2016), suggesting that the dissipation region has a coherent, strong magnetic field.
In many of these cases, the rapid gamma-ray variability is likely to be a consequence of efficient electromagnetic dissipation in a highly magnetized, relativistic outflow from the central engine (neutron star or black hole). Our group has proposed a general idea, called magnetoluminescence, to account for the most dramatic gamma-ray flaring events (Blandford et al., 2014, 2015). Magnetoluminescence refers to large scale, catastrophic conversion of electromagnetic energy into high-energy, non-thermal radiation. It could be triggered by some macroscopic, ideal instabilities in the outflow, which leads to formation of regions with non-ideal electric field that accelerate particles efficiently and process through a large volume rapidly. In the radiation-reaction-limited regime, the energy is quickly removed by radiation and, so, this might conclude with an implosion.
Currently in the literature, the process of electromagnetic dissipation in a relativistic plasma has been primarily studied in terms of magnetic reconnection at a planar current layer (e.g., Kagan et al., 2015). In such a configuration, fast reconnection and energy release proceeds through tearing instability of the current layer. It has been found that in relativistic regime, a fast reconnection rate, , can be reached, especially when the magnetization is high; efficient particle acceleration produces a power law distribution of non-thermal particles , and the power law index hardens toward 1 as magnetization increases to very high values (e.g., Guo et al., 2014; Sironi & Spitkovsky, 2014; Liu et al., 2015; Werner et al., 2016). Magnetic reconnection has been proposed to be the underlying mechanisms for the Crab flares (e.g., Uzdensky et al., 2011; Arons, 2012; Cerutti et al., 2013; Baty et al., 2013), and mini-jet models based on magnetic reconnection are invoked to explain the rapid variability from AGN jets (Giannios et al., 2009; Nalewajko et al., 2011). However, most of the existing kinetic simulations of magnetic reconnection start from a Harris-type current sheet that is already on kinetic scales; this is quite artificial and does not capture the more typically dynamic nature of current sheet formation and evolution.
Our group has been investigating new configurations that could be more generic in revealing the basic properties of electromagnetic dissipation in a highly magnetized, relativistic plasma. One class of examples is the force-free equilibria in a Cartesian periodic box. Previous analytical analysis and force-free/MHD simulations (East et al., 2015; Zrake & East, 2016) have found that the higher order states are generally unstable to ideal modes; they release their free magnetic energy within a single dynamic time scale while keeping the helicity constant. Recently, Nalewajko et al. (2016) carried out 2D particle-in-cell (PIC) simulations of the lowest order unstable configuration, in a mildly relativistic pair plasma. Most interestingly, they show the self-consistent formation of current layers, which evolve on Alfven time scales as determined by the global dynamics. Particles are accelerated by the parallel electric field in the current layers, as well as second order processes during the late stage of relaxation. The non-thermal particle fraction increases with the magnetization. Similar studies have also been carried out independently by Lyutikov et al. (2016) lately.
In this work, we extend our 2D kinetic simulations to include synchrotron radiation reaction self-consistently using the PIC code Zeltron, in a regime where individual particles are ultrarelativistic. We extract the detailed radiation signatures and study the dynamical consequences of radiation reaction systematically. Some aspects of the radiation effect have been investigated before in PIC simulations of Harris-type reconnection (Cerutti et al., 2012, 2014; Kagan et al., 2016); here we see new features as a consequence of the dynamic evolution of the current layers. Also, for the first time, we calculate the polarization signals from the PIC simulations.
The rest of the paper is organized as follows. We introduce our setup in §2, and present the results in §3, including the evolution of the instability (§3.1), mechanisms of particle acceleration and electromagnetic dissipation (§3.2), the radiation signatures (§3.3) and the dynamical consequences of radiation reaction (§3.4). We discuss the astrophysical applications in §4 and present our conclusions in §5.
2. Setup of the problem
2.1. Kinetic description, scalability and numerical method
For a relativistic pair plasma, the evolution of electromagnetic field is governed by the Maxwell equations with the charge and current densities given self-consistently by the distribution function for each species as
where and are the proper velocity and Lorentz factor of individual particles, respectively. The distribution functions satisfy the conservation of particles in phase space
and the acceleration of individual particles is determined by
We use the Abraham-Lorentz-Dirac formalism for the radiation reaction force (Jackson, 1999)
where is the proper time. In the case of synchrotron radiation or Inverse Compton (IC) scattering in the Thomson regime, in the ultrarelativistic limit, the first term (jerk) can be neglected and the radiation reaction force becomes essentially a drag opposing the particle motion. For synchrotron radiation,
where is the component of acceleration that is perpendicular to the particle velocity. For IC in the Thomson regime with isotropic background photon field,
where is Thomson cross section and is the soft photon energy density. In the ultra-relativistic limit (), the distribution function is scalable with particle energy , therefore the equations can be made dimensionless if we measure lengths in terms of the system size , time in units of , and let , , , :
and is the nominal Larmor radius, is the nominal magnetization parameter, is the nominal magnetic energy density, is the background photon energy density, characterizes the strength of synchrotron cooling
Here the characteristic frequency of synchrotron radiation is , particle gyro frequency is , and the synchrotron radiation reaction limit occurs at MeV. Thus, the evolution of the system is determined by the dimensionless parameters , , and , making the problem fully scalable.
We use the electromagnetic, relativistic Particle-In-Cell code Zeltron
2.2. Initial conditions
Note. – The columns from left to right are: run id, the grid resolution , the parameter , the ratio between the light crossing time scale and the cooling time scale , the growth rate of the electric field during the linear evolution stage where , energy partition measured at —including the change of magnetic energy , electric energy , the change of particle kinetic energy , total radiated energy , and the efficiency of high energy radiation , where is the initial magnetic energy contained in the simulation domain, and is the peak synchrotron frequency of the initial distribution.
We consider the lowest-order unstable modes of linear force-free equilibria in a 2D periodic box with size . One particular example is described by
where and the field satisfies
The topology of the equilibrium can be readily seen as two pairs of flux tubes on which the magnetic field lines wind helically with opposite orientation but the same sense of helicity (Fig. 1). The magnetic field magnitude is not uniform in the box: its peak is and the root-mean-square value is .
In the equilibrium and the distribution function for each species should satisfy steady state Maxwell-Vlasov equation (neglecting radiation reaction for the moment):
In practice, we follow the approach of Nalewajko et al. (2016) and approximate the particle distribution function as
where is the cosine of particle pitch angle with respect to local magnetic field, such that equations (22)(23) and the first two moments of equation (21) are satisfied. We assume that the particle distribution is uniform in space; in energy space satisfies the Maxwell-Jüttner distribution
where is related to the temperature as . Since we are working in the ultrarelativistic regime, so the distribution function becomes simply
which is free from the rest mass scale. The angular distribution is assumed to be
which guarantees that the pressure is isotropic and uniform, while the average drift velocity is
Assuming that electrons and positrons have equal number density and drift in opposite directions, we get from Equation (23)
We set the initial number density using the nominal values and ( to make sure ):
With this assumption, the nominal magnetization parameter as defined before becomes
where is the characteristic gyro radius of a thermal particle. We note that the ratio between the skin depth and is
Therefore, the computational grid should resolve or electron skin depth, depending on the value of .
Given the functional form of , the evolution of the system is now fully determined by the following dimensionless parameters: , and two of the other three parameters: , and . Among the latter we choose and : indicates the initial partition between magnetic energy and particle kinetic energy, while is a parameter characterizing the charge multiplicity, or equivalently, level of opposite drifting between electrons and positrons: small implies large multiplicity/small relative drift. We have performed a series of simulations with different values for these dimensionless parameters, to explore the underlying scaling relations. In this paper, we focus on one particular choice of and , consider only synchrotron radiation, and explore the dynamical effects of different as well as the radiation signatures. Specifically, we set , , —the corresponding box-averaged, warm magnetization is , and , where is the enthalpy density of the plasma. Specification of for each run is listed in Table 1. We use a grid resolution , correspondingly . We’ve tested the convergence using both lower and higher resolution grids with and we find that our main results do not depend on the resolution. The number of particles per cell for each species is 144 (slightly larger than Cerutti et al., 2013; Nalewajko et al., 2016). The time step is set to to satisfy the CFL stability condition (Courant et al., 1967). During the evolution, the instability as found by East et al. (2015) does grow spontaneously from the initial conditions we use, but to reduce the computational cost in most of the simulations we introduce a small amplitude, long wavelength perturbation of the form when , and turn on radiation reaction at the same time.
3. Results of the simulations
In this section, we first consider the case where the radiation reaction is not yet dynamically important (run 1). We summarize briefly the evolution of the instability in §3.1, then discuss in §3.3 the radiative signatures of the instability, including variability, spectra and angular distribution, as well as polarization. We analyze the dynamical consequences of radiation reaction systematically in §3.4.
3.1. Evolution of the instability–mildly radiative cases
The overall evolution of the field configuration is qualitatively similar for the range of parameters we explored, and also resembles the mildly relativistic cases studied by Nalewajko et al. (2016). Take Run 1 as an example. Figure 2(a) shows the evolution of different energy components. At the start of the simulation, the random noise in charge density due to finite number of particles cause the electric energy to grow to a plateau times the initial total magnetic energy . The initial small fluctuations roughly settle down after one light crossing time scale, and the instability grows from the perturbation.
During the linear growth of the instability, the two pairs of flux tubes with the same sign of start to merge; each pair produces a layer in between with enhanced current and number density (Figure 3, first column). The electric field grows exponentially as , with in this particular case. Non-ideal regions with start to be produced in the current layer, while remains less than due to the advection and compression of magnetic flux into the current layers. The exhausts of the current layers are the locations where becomes maximal (Figure 4). Some charge density appears surrounding the flux tubes, because the in-plane motion of the flux tubes gives rise to a Lorentz force that moves positive and negative charges in opposite directions along the field lines, accumulating different signs of charges at different sections of the flux tube boundaries. The boundaries separating the flux tubes where relative motion occurs turn out to be tangential discontinuities.
The current layer continues to get thinner and longer, as plasmas are pushed into the layer from both sides and ejected along the exhausts. Eventually the thickness of the current layer is determined by the skin depth of the plasma in the sheet: . When the aspect ratio of the current layer becomes large enough, plasmoids—namely, small islands with local concentration of current in this 2D case—start to form and grow in the current layer (Figures 3 and 4, second and third columns). They, too, get ejected from the ends of the current layer, and when they collide with the ambient magnetic field, secondary current layers are produced, as well as fast waves that propagate into the neighboring magnetic domain.
The initial current sheet only lasts for about one dynamical time scale. It gets destroyed as the overall field structure continues on a large amplitude damped oscillation, with part of the energy going back and forth between magnetic form and electric plus kinetic form. During the oscillation, transient current layers with are still produced. In the end, the system roughly settles into the longest wavelength equilibrium, with a relatively dense, cool plasma near the magnetic separatrices and a hot, dilute plasma near the center of the magnetic domain (Figure 3, last column).
During the whole evolution, deviation from helicity conservation is less than 0.6% (Figure 2b). At the end, of the magnetic energy has been dissipated. As a comparison, theoretically the initial state has an energy times that of the true ground state (the longest wavelength solution allowed in the periodic box), so the maximum amount of magnetic free energy is . The actual released magnetic energy is smaller because in 2D there are additional topological constraints making the final state different from the true ground state (Zrake & East, 2016).
3.2. Mechanisms of particle acceleration and dissipation of electromagnetic energy
The isotropic particle spectra at a series of simulation times are shown in Figure 5. Consistent with Nalewajko et al. (2016), we find that high energy particles are first accelerated in the initial current layers by the nonideal electric field , forming a bump on the tail of the distribution. This appears more clearly after we subtract the thermal Maxwellian component (Figure 5a). The high energy non-thermal tail expands in energy range and the spectrum gets harder as the current sheet continues to stretch. Eventually the tail reaches an extent of about one decade in energy, within which a roughly power-law-like distribution is established. The hardest power-law index is achieved when the main current sheets are in the plasmoid dominated reconnection stage. After the main current layers dissolve, the system enters a turbulent relaxation process, particles diffuse in momentum space, gradually forming a steeper power law connecting the high energy bump with the thermal Maxwellian. There’s also an overall heating of the background plasma (see also, Lyutikov et al., 2016). At the end of the simulation, the power law has a spectral index , and the nonthermal particles comprise 12% in number and 36% in energy. Due to the modest in our simulations, the particle spectrum is relatively soft. It has been shown by Nalewajko et al. (2016) that the spectrum does get harder as increases.
Figure 6 shows the rate of energy transfer from the electromagnetic field to the particles, , as well as that contributed by and —the components of electric field parallel and perpendicular to the magnetic field, respectively. Since the initial magnetic configuration contains null points where , at these locations we let . However, we find that as the evolution starts, is advected into the current layer, so does not vanish in the current layer, is always less than , and , are well defined. It can be seen that although dominates at earlier times when the ideal instability just starts to develop, most of the dissipation happens upon the saturation of the instability and is dominated by . Figure 7 shows maps of and at representative time points. We observe that mostly operates at the primary current sheets, and also the secondary current layers formed when the plasmoids collide with neighboring flux tubes.
3.3. Radiation signatures—mildly radiative case
As a benchmark example, in this section we analyze the radiation signatures of the mildly radiative case, run 1. While the total emitted power as a function of time is quite smooth overall (Figure 2a), the radiation received by a fixed observer can be highly variable, as we now show.
shortest variability time scale
We calculate light curves in different wavebands as received by observers located on the plane. Observers off the plane are not included because the z translational invariance makes it ambiguous to define the emitting/receiving time for the radiation coming into/out of the plane. But fortunately, emission within the plane is already revealing most of the generic behaviors. In Figures 8, 9 and 10 top panel, three examples are shown where the observers are looking from , counterclockwise of , and , respectively. Several remarkable features can be noticed: (1) very sharp, high intensity peaks, with durations (well above the time resolution we use to calculate the light curves), are seen in certain directions. In particular, peaks from positrons are typically observed around directions, and those of electrons close to directions; somewhere in between one may not see any of the sharp peaks. (2) These peaks only appear during the early evolution of the system—between the saturation of the linear stage and the destruction of the first current layers. (3) High frequency radiation is dramatically more variable than the low frequency radiation.
The third point is also made clear by the power spectra of the light curves, shown in Figure 11. The power spectra in different wavebands can be reasonably fitted by power laws, which get shallower as we go to higher energy band. Such behavior is similar to what has been observed in Harris-type reconnection simulations (Cerutti et al., 2013). In particular, we notice that for radiation emitted along and directions, there appears to be some power excess on the time scale .
In the following we investigate carefully the fastest variability seen in the light curves. In Figure 12, we zoom in around the three biggest peaks in the light curve, and compare them directly with the emissivity spacetime diagram (the emissivity has been integrated over y). After identifying the emission event, we pinpoint the responsible features spatially on the 2D emissivity map and total synchrotron power map (Figure 13a). It turns out that these are features across, moving toward the observer with a speed close to c. These are high energy particle beams ejected from the ends of the current layers, and at this particular time point their trajectories are almost tangent to the line of sight. Here compactness, beam sweeping and light travel time effect are among the factors that cause the sharp emission peaks. In particular, if we look at the three bright tracks on the spacetime diagram that are responsible for the three highest peaks in the light curve (Figure 12), their span on the emitting time axis can be associated with the time it takes for the opening angle of the particle beam to turn through the observer receiving angle (), but as the tracks lie almost parallel to the light cone, due to the light travel time effect, the beam sweeping time becomes negligible on the receiving time axis; the actual variability time scale is determined by the instantaneous spatial extent of the emitting structure. We have verified that the measured duration of the spikes does not depend on the size of the receiving angle we use for the light curve calculation, as long as it is small enough ().
Another point to notice is that, these spatially compact emitting structures form following the merging of a plasmoid into the surrounding magnetic field. From Figure 3 row 3 and 4 we see that the plasmoids contain dense concentration of high energy particles. These plasmoids move out at a speed that reaches ; the corresponding induced electric field on the back side of the plasmoid tends to accelerate particles along the current while that on the front side tends to decelerate particles. High energy particles are thus bunched by the plasmoids; they experience a further kick at the secondary current layers that form when the plasmoids collide with the ambient field, getting compressed/elongated in a similar manner as the secondary current layers. The spatial bunching and compactness is important in producing the sharp, high magnitude peaks in the light curves, as these are not observed when the ejecta are more diffuse.
Furthermore, as electrons turn counterclockwise and positrons turn clockwise on the plane when they exit the current layer (for the particular instance we are looking at), observers see emission peaks from electrons mostly near and positrons near . Figure 13(b) shows some aspects of the particle trajectory bending and beam bunching. We plot on top of the magnetic field lines the position of tracked high energy positrons, and compare their energy and synchrotron power. These particles are selected such that at the end of the simulation they all have ; they are only the most energetic of all the positrons within the simulation domain. It can be seen that particles that eventually reach high energies are first accelerated in the biggest current layers, where they do not radiate much; synchrotron radiation only becomes important when particles are ejected from the current layers and their trajectories start to bend significantly. Figure 14 shows the trajectory and energy history of a few particles that are part of the bunch responsible for the highest peak in the light curve. The sharp rise in and near the dashed line is a result of increased curvature in particle trajectory as they get out of the current layer, and the subsequent quick drop in is not because of cooling (particle energy is almost unchanged) but due to reduced curvature, as indicated by the perpendicular acceleration felt by the particles in Figure 14.
The energetics of a single spike can be estimated based on its duration and peak flux . We have seen that is determined by the scale of the ejected plasmoid , which varies depending on the history of the plasmoid. Take the sharpest spike as an example: we have and the peak flux is sr (Figure 12), where is the initial total magnetic energy. An observer who assumes the radiation to be isotropic would deduce . If the observer has a knowledge of the average magnetic field, she could calculate the magnetic energy contained in a volume whose size is determined by the variability time scale, and get . Thus, the observer would conclude that the apparent radiative efficiency is . Of course this depends on as will be shown in §3.4.
The result can be understood from a simple physical argument. Suppose that the average particle number density in a plasmoid is and average Lorentz factor is , upon the destruction of the plasmoid particles start to turn in a magnetic field of strength , so the observed fluence should be where is the beaming angle, determined by the velocity dispersion in the beam. Meanwhile, the inferred magnetic energy content based on the variability time scale is . Assuming the plasmoid gas pressure is times the ambient magnetic pressure before its destruction, namely , we have . In the above example, , , and from Figure 15, the beaming angle is , so we have . could be of order 1 or larger as the pressure is highly anisotropic and it is the momentum that dominates. Thus, the simple estimation is roughly consistent with the above measurement.
Despite the attractive high peak intensity and fast variability produced by the small plasmoids, the total energy involved is small. Depending on the viewing angle, an observer may see emission with longer time scales and larger total energetics albeit their small peak flux. The peaks with longer time scales are produced by more diffuse ejecta or smooth field structures lit up by distributed high energy particles, which evolve on dynamic time scales. Another remark to be made is that, in this particular example the current sheet itself is not rotating, but in reality it can be dynamic and can turn around during the evolution. This could introduce further variability.
Origin and beaming of high energy radiation
During the saturation of the linear instability and the early stage of nonlinear evolution, particles are being efficiently accelerated in the current layers by the parallel electric field. Plasmoids are not very good particle accelerators except for those undergoing rapid acceleration themselves, but they could trap particles that are accelerated in the current layer. High energy particles within the current layer have a fan-like angular distribution spanned around ; they turn toward the ends of the current layers due to the reconnected magnetic field. These fan-like features are readily seen in Figure 16, especially in high energy bands. However, particles within the current layers do not produce a large amount of high energy synchrotron radiation—this is demonstrated by the absence of corresponding fan-like structure in the radiation angular distribution (Figure 15). Plasmoids do radiate a significant amount of power, as shown in the synchrotron power map (Figure 3 last row), but this radiation peaks in relatively low frequency—in Figure 15 we see emission from plasmoids mainly in the intermediate energy band. The main reason for the lack of high energy radiation from the current layer itself (plasmoids included) is that the curvature of the particle trajectory is small.
Significant high energy radiation emerges from just downstream of the exhausts of the current layers, where high energy particles are being dumped onto the surrounding magnetic field. The ejecta can be intermittent and compact, as the spatial distribution of high energy particles is modulated by the spontaneously formed plasmoids. In particular, the unrelaxed, small plasmoids formed at the late stage of current sheet evolution produce the most compact, high intensity emission upon their destruction; as the beams turn around in the magnetic field, they give narrow stripes in the angular distribution of radiation (Figure 15) and we see sharp peaks when the beams sweep across the line of sight. On the other hand, the fully relaxed, large plasmoids produce relatively diffuse emission.
Later on, the initial current layers get destroyed and the system evolves in a more turbulent way. High energy particle beams are gradually dispersed both in configuration space and in momentum space. In this particular example particles do not cool significantly; high energy particles eventually spread over most of the simulation domain. As a result, the high energy radiation becomes more and more diffuse and rapid variability will no longer be observed.
Time dependent spectrum of observed radiation
The quiescent radiation spectrum peaks at Hz in the ultraviolet; as the instability develops and accelerates a significant fraction of particles to energies , a new radiation component also emerges. As an example, we divide the time series of radiation received by the observer located at into equally spaced time windows and calculate the synchrotron spectrum within each of these windows. Figure 17 shows the results both for large time windows throughout the simulation duration and for small time windows around the highest peaks. It can be seen that during the evolution, the peak frequency of the high energy radiation can reach times the quiescent value.
Although the overall particle spectrum as shown in Figure 5 only exhibits a steep power-law tail, the instantaneous distribution of particles moving along a certain direction can have a dominant component on the tail. Figure 18 shows the distribution of particles within around near the time point when the highest peak in the light curve is produced. Evidently the positrons have an almost mono-energetic beam with .
We calculate the linear polarization degree and polarization angle as a function of time, in different wavebands, for observers located on the plane. Examples are shown in Figures 8, 9 and 10. The initial equilibrium produces a high polarization and the polarization angle is aligned with the plane, as one would expect from the symmetry of the configuration. At the start of the instability, there is an overall drop in the polarization degree, which is especially obvious in the low energy band. However, for the high energy radiation, during the “flares” (i.e. sharp peaks in the light curves), polarization degree increases significantly, accompanied by large change in polarization angle. This is because the high energy emitting particle beams are compact and they sample a relatively strong, ordered magnetic field, and their high radiation flux outweighs the contribution from the other parts of the simulation domain. As an example, we looked at the polarization produced by the bunch of particles we tracked in Figure 14. We notice that these particles have gyro radius that is larger than or comparable to the scale length of the field, and the polarization angle is determined by the instantaneous orientation of the orbital plane instead of the local magnetic field, because of the presence of the electric field. The calculated polarization angle of the emission from a typical particle in the bunch is about 64, as shown by the magenta triangle in Figure 8. Since this bunch dominates the total intensity, the resultant polarization angle of the whole box deviates from the nominal value and reaches around . The polarization angle of the low energy component shows only slight change during most of the evolution. When the system is fully evolved to settle into a new equilibrium (the ground state), the all-frequency integrated polarization degree again settles at and the polarization angle also returns to its initial value.
3.4. Dynamical consequences of radiation reaction
In this section, we do a systematic comparison between the runs with different values of .
Total emitted energy
A comparison of the total emitted power as a function of time for runs with different is shown in Figure 19(a). In order to bring the synchrotron power from different runs to the same scale, we have divided by in the plot. It can be seen that when is sufficiently small (), the evolution of is very similar, suggesting that in this regime radiation reaction does not affect the dynamics and the radiated energy simply scales with . However, when is large enough such that accelerated particles get close to the radiation reaction limit, the scaling deviates from a linear relation with : the back reaction of synchrotron radiation suppresses high energy particles, thus reduces the radiative output. We have tracked the energy for different components in the system, and the results are shown in Table 1. The radiated energy is directly drawn from particle kinetic energy; in the high runs, the kinetic energy sees a negative growth due to the strong cooling.
In run 3 and 4, we also see sharp spikes in the received radiation light curves when the viewing directions are close to or , with similar variation time whereas the peak scales with . From an observer’s point of view, these spikes indicate radiative efficiency and 10, respectively.
Evolution of the magnetization
Figure 19(b) shows the evolution of the ratio between the electromagnetic energy and particle kinetic energy as a function of time for runs with different . When , the evolution of the instability heats up the particles, correspondingly we see a drop in . On the other hand, when starts out large such that the cooling time scale starts to be comparable with the dynamic time scale of the system, before the instability saturates the radiative cooling of particles reduces while the field is still largely unchanged (since reducing the Lorentz factor of ultrarelativistic particles does not change the current), as a consequence increases. The saturation of the instability again reduces , but the high level of radiative cooling keeps from dropping to very low values.
Nalewajko et al. (2016) have shown that the growth rate of the instability depends on (and possibly ). An effective increase in due to radiative cooling will affect the growth rate. We measured the growth rate during the linear evolution where the electric field can be written as ; this is shown in Table 1. The slight increase in growth rate can be attributed to increased effective magnetization as a result of radiative cooling. This is also consistent with the results of Cerutti et al. (2013), who reported that the reconnection rate increases in the strong cooling regime.
In Figure 20(a) we plot the relative change of helicity for runs with different . It can be seen that for all the runs the deviation from helicity conservation is smaller than 0.6%. Since the evolution of helicity is determined by , the negligible change in helicity means that the volume integral of is small. We also checked the average and maximum values of as a function of time, shown in Figure 20. Although locally can reach rather high values, the mean is small, indicating that the volume with non-negligible is insignificant. And the evolution of does not have any strong dependence on . This seems to suggest that in the regime we’ve explored, synchrotron radiation reaction cannot support volumetric non-ideal electric field; the scale of the non-ideal region is still determined by the kinetic effects.
Particle spectrum and radiation spectrum
Radiation reaction limits the energy gain of the highest energy particles, as shown in Figure 21 where we plot the isotropic particle spectra at a series of simulation times for run 4. While the initial direct acceleration in the current layers produces a high energy bump in a similar fashion as run 1 (Figure 5), this new component loses its energy within roughly one light crossing time. Different portions of the light curve and time dependent synchrotron spectra, as seen by an observer located at , are shown in Figure 22. Comparing with run 1 (Figure 17), we see that the peak emission frequency at the flux maximum is still efficiently boosted by orders of magnitude in this case, although bumping close to the radiation reaction limit 160 MeV. This is because the direct electric field acceleration is fast enough and the radiative loss in the current layer is not significant, as we will show in the following subsection.
We note that when the cooling is fast, observing directions on the plane may not be optimal for seeing the beamed highest energy emission. This is because the particle beams that are ejected from the current layers are initially pointing off the plane with projections on the plane along (e.g., Figure 16); these particles then turn toward the , directions, and the strongest radiation comes from some point in between. In Figure 23, we show a few examples of instantaneously radiated spectrum along specific directions off the plane, for run 4 and run 5. (This is not the spectrum received by an observer as we cannot define the receiving time unambiguously in this 2D case.) In particular, for the highest run we have (Figure 23 right panel), the peak frequency of synchrotron radiation can reach the radiation reaction limit.
In mildly radiative cases, we have shown in Figures 13 and 14 that the highest energy particles are initially accelerated by the parallel electric field in the current layers, where they follow Speiser orbits (e.g., Uzdensky et al., 2011; Cerutti et al., 2013), while at the same time being deflected toward the ends of the current layers. When they are deep in the current layer, the radiative loss is small due to reduced curvature as a result of acceleration; radiation only becomes significant when they get out of the current layer and start to get bent in stronger ambient magnetic field. Later on, during the large scale oscillation of the fields, some particles continue to gain energy when they bounce off an expanding magnetic domain from outside. In this phase, the acceleration is more stochastic.
When starts off large, the first acceleration phase, namely the action of parallel electric field, is still efficient enough to boost particles to high energies, but the second stochastic phase cannot compete with the cooling. Figure 24 shows a few representative particle trajectories and their energy history. The peaks in (e.g., near the vertical dashed line) are due to sudden increase in curvature when the particles get out of the current layer, similar to that in Figure 14. The difference in this case is that particles are cooling significantly at the same time, so the peak drops off very quickly. After the dissolution of the current layer, all the high energy particles cool rapidly due to radiative loss and no efficient energization by stochastic processes is seen.
In our simulations, the initial configuration is an unstable force-free equilibrium that can get destroyed over a single dynamic time scale, so a natural question one might ask is whether such a structure can form in the first place. In a realistic astrophysical environment, the situation is much more complicated since there’s continuous motion, energy/mass injection and/or loss, so the underlying equilibrium keeps evolving ceaselessly. We imagine that a few processes may build up plasma configurations of high magnetic free energy. The first example is that, in a pulsar wind, there’s an equatorial current sheet, and random reconnection can happen at locations that are causally disconnected. This could lead to the development of highly tangled flux rope structures globally. As the flow expands, these tangled structures get frozen out as different parts stay out of causal connection. When the flow slows down eventually, they get back into causal contact and the configuration may look like some of the higher order equilibria. The second possibility is that, the pulsar wind is initially striped with a wavelength much shorter than the termination shock radius , e.g. cm while cm for the Crab. These stripes—small scale fluctuations—continue to go through inverse cascade, forming larger and larger flux ropes, eventually reaching a large enough scale that is relevant for the acceleration of PeV particles (Zrake, 2016). Another possibility involves the polar jet which has relatively high magnetization and is kink unstable—suggesting that it possesses free magnetic energy. Recently Lyutikov et al. (2016) also proposed that the intermediate latitude post shock flow can produce regions with high magnetization and current carrying filaments, resembling the flux tubes in force-free configurations.
Based on these considerations, we think it instructive to use the unstable force-free equilibria as a simple testbed to study the subsequent particle acceleration and radiation. With a kinetic approach, we are able to self-consistently extract the radiation signatures and study systematically the effect of radiation reaction. We find that the simple model is teaching us a lot about the generic properties of electromagnetic dissipation and radiation in magnetized, relativistic plasmas.
The first remarkable feature is the rapid variability observed during the evolution from an initially smooth configuration. Though no current layers or other singular structures are embedded in the starting equilibrium state, its free energy is released through an ideal instability driven by the flux tubes’ tendency to merge, which inevitably forces the formation of a current layer around the X-point. This is consistent with Syrovatskii’s theory of current sheet formation via X-point collapse (Imshennik & Syrovatskiǐ, 1967; Syrovatskii, 1966; Nalewajko et al., 2016). Kinetic scales become important here: the thickness of the current layer is dictated by the plasma skin depth (Nalewajko et al., 2016); the formation of plasmoids, in the situation of collisionless plasma, is also determined by kinetic processes—the tearing instability growth rate is roughly proportional to the plasma frequency (Zelenyi & Krasnoselskikh, 1979), though the plasmoids can grow to macroscopic scales. The rapid variability in synchrotron radiation is produced by high energy particle beams that are ejected from the ends of the current layers, having been bunched by the tearing instability. This also leads to high energy radiation being much more variable than low energy radiation. The picture is in some sense similar to Cerutti et al. (2013). We have shown that for a single pulse, the apparent radiation efficiency can be relatively high, especially when is large. In particular, Figure 21 for run 4 () shows that the peak power per steradian seen by an observer at reaches where is the total magnetic energy contained in the simulation domain at , and the (isotropic) radiation efficiency is .
This makes it attractive to associate these spikes with the Crab flares, but some caveats need to be taken into account. In the Crab flare case, we have shown that . This is still about two orders of magnitude larger than the maximum we got from the current simulations, and we may need to extrapolate the results using simulations with different magnetization, box sizes and in future studies. In addition, if we do make a literal identification between a single spike from the simulation and a Crab flare event, this would mean that the size of the emitting PeV particle bunch is on the order of —comparable to the Larmor radius of the PeV particles cm, and the simulation box would correspond to a length scale of cm—comparable to the radius of the termination shock. Such a demanding requirement is related to the intrinsic problem of insufficient energy contained in the kinetic scale beams in our current simulations. Recalling that the energy in the beams comes from the volume that has been processed by the current layers, one would hope that going to large limit where both the instability growth rate and the magnetic energy density are larger, the energy deposited into the high energy particle beams would be more promising. At the same time, we should keep in mind that the scale separation in our simulation is still quite far from the astrophysical reality. In the Crab, the pressure is dominated by TeV particles, whose gyro radius is cm, and the plasma skin depth would be much smaller than (unless ). In contrast, our simulation has a box size that is only and the emitting region size only a few . (A larger also means larger charge multiplicity for a fixed . The multiplicity is defined as the ratio between the actual number of pairs and the minimum number needed to support the current. In our setup essentially , see Equation 29.) It is desirable to test the regime of large in the future with additional computational power (c.f. Sironi et al., 2016, but it’s a quite different setup with Harris current layer and cold, nonrelativistic background plasma). One further catch is the relatively soft spectral index for the nonthermal particle distribution (isotropic) we get here. The biggest Crab flare has a spectrum Buehler et al. (2012), which requires a particle distribution with , or mono-energetic if the multi-wavelength constraint is taken into account (Weisskopf et al., 2013). Though the angle dependent particle distribution could have a harder spectrum (§3.3.3), the actual anisotropy needs to be tested using more realistic scale separation. Nalewajko et al. (2016) found that the spectral index decreases as increases for a similar setup, and such a trend is also true in Harris layer reconnection (Sironi & Spitkovsky, 2014; Guo et al., 2014; Werner et al., 2016), so higher simulations are promising.
The second important feature we notice is that as the dominant acceleration mechanism here is parallel electric field acceleration, the curvature of accelerated particles becomes small when they are in the current layer and significant synchrotron radiation only takes place after they get out of the current layer. Such a separation between the acceleration site and radiation site would allow the highest energy particles to get beyond the radiation reaction limit (Uzdensky et al., 2011; Cerutti et al., 2013). The premise is that the available electric potential should be enough to boost particle energy from the thermal sea to the radiation reaction limit. This is challenging to realize in PIC studies due to the limited dynamic range. In our simulations we have chosen a modest sigma and the relaxation process is able to produce highest energy particles with , so in order to test whether we can get over the radiation reaction limit, we tried to artificially increase for the thermal sea, or equivalently, increase the energy of the thermal particles. We find that for as large as , we can get high energy synchrotron radiation instantaneously peaking at the radiation reaction limit. However, when starts out larger than , the thermal population would cool down back to around at the saturation of the instability so it no longer makes sense to further increase at the beginning of the simulation. As a comparison, for the Crab, if the majority of the particles are at TeV energy range and the average magnetic field is 1 mG, then . It remains to be seen whether higher configurations could do a much better job. We expect this to be the case as Nalewajko et al. (2016) find that both the fraction of non-thermal particles and the maximum energy of accelerated particles scale with . This will be tested in the near future.
As we change the cooling parameter , we did a comparison among the various cooling regimes to see the effects of radiation reaction on the dynamics. The strong synchrotron cooling does have a significant impact on the highest energy particles, but since they are not energetically dominant, this hardly has any effect on the global dynamics. We tried to decompose the non-ideal electric field using a generalized Ohm’s law (not shown here), and find that synchrotron cooling is not contributing any noticeable resistivity in the regime we’ve explored. This is partly because the majority of particles are not yet reaching the radiation reaction limit so that the radiation reaction force on any plasma volume remains much smaller compared to the other force terms (inertia, pressure gradient and Lorentz force), partly because the location where non-ideal electric field arises are not locations where synchrotron radiation is most significant. The situation could be quite different if the particles that are counter streaming along z in the current layer excite gyro-resonance instability, which enhances the synchrotron radiation within the current layer, or if the dominant radiation mechanism is not synchrotron but inverse Compton—in these cases radiation reaction force might contribute significantly in supporting non-ideal electric field over some extended volume (e.g., Uzdensky, 2016).
Besides the modest and modest scale separation , another limitation in our current simulations is the 2D constraint. It does not allow for any kink instability or other variations along z direction (cf. Cerutti et al., 2014), and the acceleration length along z is unlimited. In addition, the accelerated particles are counter streaming in the z direction—this could lead to excitation of collective modes (e.g. gyroresonance) that contribute to anomalous resistivity to support non-ideal electric field (e.g. Treumann & Baumjohann, 1997). Also, in the strongly cooling runs, the most interesting directions for fast time variability and efficient synchrotron production might be out of the simulation plane. These will need to be tested in 3D simulations. However, 3D runs are much more computationally expensive at the moment. We plan to carry this out in the near future.
We have performed 2D PIC simulations of a magnetostatic equilibrium that belongs to the lowest order unstable force-free states. We work in the regime where the individual particles are ultrarelativistic, and include synchrotron radiation reaction self-consistently. A benchmark example, where and , is examined in detail to obtain the radiation signatures, and a systematic study of different cooling regimes is performed to understand the dynamical effects of radiation reaction.
We find that the evolution of the system is consistent with previous force-free, MHD (East et al., 2015) and PIC simulations (Nalewajko et al., 2016, mildly relativistic cases). The ideal instability eventually leads to current sheet formation, where most of the particle acceleration and electromagnetic dissipation happens. Regions with develop in the current layers, but the volume they occupy is negligibly small. As a result, helicity is pretty well conserved during the whole evolution. We also do not see regions with develop, because of the advection of guide field into the current layers.
The highest energy particles are first accelerated in the current layers by the parallel electric field, where they do not radiate much due to the small curvature of their trajectory, despite the presence of guide field in the current layer. Most of the radiation is produced when particles are ejected from the current layers—their trajectories start to bend significantly in the ambient magnetic field which changes direction at the end of the current layer. Such a separation between acceleration site and synchrotron radiation site could in principle facilitate acceleration beyond the synchrotron radiation reaction limit.
We find that the fastest variability in synchrotron radiation is produced when compact plasmoids that contain high energy particles are ejected from the ends of the current layer and get destroyed. These give beamed radiation as the particles released from the plasmoids start to turn in the ambient magnetic field. An observer sees high intensity radiation when the beam happens to be aligned with the line of sight. As a result, the high energy radiation is much more variable than the low energy radiation, and these flares are accompanied by an increase in the polarization degree and rapid change of polarization angle in the high energy band. The variability time scale is determined by the spatial extent of the emitting structure. In our simulations, this can be as short as , and the peak flux per steradian can reach in the case , giving an apparent radiation efficiency . However, the total energy involved in these spikes is small. Though this setup, with the parameters that were numerically tractable here, is not yet enough to directly reproduce the feature of the Crab flares, this work suggests that runs at higher are promising. We plan to explore this in future work.
- affiliation: 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
- affiliation: 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
- affiliation: Nicolaus Copernicus Astronomical Center, Bartycka 18, 00-716 Warsaw, Poland
- affiliation: NASA Einstein Postdoctoral Fellow
- affiliation: 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
- Abdo, A. A., et al. 2010, Nature, 463, 919
- Abdo, A. A., et al. 2011, Science, 331, 739
- Aharonian, F., et al. 2007, ApJ, 664, L71
- Albert, J., et al. 2007, ApJ, 669, 862
- Aleksić, J., et al. 2011, ApJ, 730, L8
- —. 2014, Science, 346, 1080
- Arons, J. 2012, Space Science Reviews, 173, 341
- Baty, H., Petri, J., & Zenitani, S. 2013, MNRAS, 436, L20
- Blandford, R., East, W., Nalewajko, K., Yuan, Y., & Zrake, J. 2015, arXiv:1511.07515 [astro-ph], arXiv: 1511.07515
- Blandford, R., Simeon, P., & Yuan, Y. 2014, Nuclear Physics B Proceedings Supplements, 256, 9
- Blinov, D., et al. 2015, MNRAS, 453, 1669
- Buehler, R., & Blandford, R. 2014, Reports on Progress in Physics, 77, 6901
- Buehler, R., et al. 2012, The Astrophysical Journal, 749, 26
- Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2012, The Astrophysical Journal Letters, 754, L33
- Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147
- —. 2014, ApJ, 782, 104
- Courant, R., Friedrichs, K., & Lewy, H. 1967, IBM Journal of Research and Development, 11, 215
- East, W. E., Zrake, J., Yuan, Y., & Blandford, R. D. 2015, Physical Review Letters, 115, 095002
- Giannios, D., Uzdensky, D. A., & Begelman, M. C. 2009, MNRAS, 395, L29
- Guo, F., Li, H., Daughton, W., & Liu, Y.-H. 2014, Physical Review Letters, 113, 155005
- Hayashida, M., et al. 2015, ApJ, 807, 79
- Imshennik, V. S., & Syrovatskiǐ, S. I. 1967, Soviet Journal of Experimental and Theoretical Physics, 25, 656
- Jackson, J. D. 1999, Classical electrodynamics, 3rd edn. (New York: Wiley)
- Kagan, D., Nakar, E., & Piran, T. 2016, arXiv:1601.07349 [astro-ph], arXiv: 1601.07349
- Kagan, D., Sironi, L., Cerutti, B., & Giannios, D. 2015, Space Sci. Rev., 191, 545
- Kiehlmann, S., et al. 2016, arXiv:1603.00249 [astro-ph], arXiv: 1603.00249
- Landau, L. D., & Lifshitz, E. M. 1975, The classical theory of fields, 4th edn. (Oxford ;: Pergamon Press,)
- Liu, Y.-H., Guo, F., Daughton, W., Li, H., & Hesse, M. 2015, Physical Review Letters, 114, 095002
- Lyutikov, M., Sironi, L., Komissarov, S., & Porth, O. 2016, arXiv:1603.05731 [astro-ph, physics:physics], arXiv: 1603.05731
- Marscher, A. P., et al. 2008, Nature, 452, 966
- Nalewajko, K., Giannios, D., Begelman, M. C., Uzdensky, D. A., & Sikora, M. 2011, MNRAS, 413, 333
- Nalewajko, K., Zrake, J., Yuan, Y., East, W. E., & Blandford, R. D. 2016, arXiv:1603.04850 [astro-ph], arXiv: 1603.04850
- Rani, B., Lott, B., Krichbaum, T. P., Fuhrmann, L., & Zensus, J. A. 2013, A&A, 557, A71
- Rudy, A., et al. 2015, ApJ, 811, 24
- Sironi, L., Giannios, D., & Petropoulou, M. 2016, ArXiv e-prints
- Sironi, L., & Spitkovsky, A. 2014, ApJ, 783, L21
- Syrovatskii, S. I. 1966, Soviet Ast., 10, 270
- Tavani, M., et al. 2011, Science, 331, 736
- Treumann, R. A., & Baumjohann, W. 1997, Advanced space plasma physics
- Uzdensky, D. A. 2016, in Astrophysics and Space Science Library, Vol. 427, Astrophysics and Space Science Library, ed. W. Gonzalez & E. Parker, 473
- Uzdensky, D. A., Cerutti, B., & Begelman, M. C. 2011, The Astrophysical Journal Letters, 737, L40
- Weisskopf, M. C., et al. 2013, The Astrophysical Journal, 765, 56
- Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8
- Yee, K. 1966, IEEE Transactions on Antennas and Propagation, 14, 302
- Zelenyi, L. M., & Krasnoselskikh, V. V. 1979, Soviet Ast., 23, 460
- Zrake, J. 2016, ApJ, 823, 39
- Zrake, J., & East, W. E. 2016, ApJ, 817, 89