Synthetic Spectra from PIC Simulations of Relativistic Collisionless Shocks
DDraft version November 13, 2018
Preprint typeset using L A TEX style emulateapj v. 08/22/09
SYNTHETIC SPECTRA FROM PIC SIMULATIONS OF RELATIVISTIC COLLISIONLESS SHOCKS
Lorenzo Sironi and Anatoly Spitkovsky
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544-1001
Draft version November 13, 2018
ABSTRACTWe extract synthetic photon spectra from first-principles particle-in-cell simulations of relativis-tic shocks propagating in unmagnetized pair plasmas. The two basic ingredients for the radiation,namely accelerated particles and magnetic fields, are produced self-consistently as part of the shockevolution. We use the method of Hededal & Nordlund (2005) and compute the photon spectrumvia Fourier transform of the electric far-field from a large number of particles, sampled directly fromthe simulation. We find that the spectrum from relativistic collisionless shocks is entirely consistentwith synchrotron radiation in the magnetic fields generated by Weibel instability. We can recoverthe so-called “jitter” regime only if we artificially reduce the strength of the electromagnetic fields,such that the wiggler parameter K ≡ qBλ/mc becomes much smaller than unity ( B and λ are thestrength and scale of the magnetic turbulence, respectively). These findings may place constraintson the origin of non-thermal emission in astrophysics, especially for the interpretation of the hard(harder than synchrotron) low-frequency spectrum of Gamma-Ray Bursts. Subject headings: acceleration of particles — gamma rays: bursts — radiation mechanisms: non-thermal — shock waves INTRODUCTION
Non-thermal photon spectra from Pulsar Wind Nebu-lae, jets from Active Galactic Nuclei, Gamma-Ray Burstsand Supernova Remnants are usually explained as syn-chrotron radiation from a power-law population of parti-cles, presumably accelerated in collisionless shocks. Themicrophysical details of shock acceleration are still poorlyknown, however, and are the subject of active research.Particle-in-cell (PIC) simulations of colliding plasmashells have shown that Weibel instability (Weibel 1959;Medvedev & Loeb 1999; Gruzinov & Waxman 1999)converts the free energy of counter-streaming flows intosmall scale (skin-depth) magnetic fields (Nishikawa et al.2003, 2005; Silva et al. 2003; Frederiksen et al. 2004;Hededal et al. 2004). The fields grow to sub-equipartitionlevels and deflect and randomize the bulk flow, thus cre-ating a shock (Spitkovsky 2005, 2008a; Chang et al. 2008;Keshet et al. 2009). A few percent of the incoming par-ticles repeatedly scatter off the magnetic turbulence cre-ated by Weibel instability, and eventually populate apower-law high-energy tail in the particle spectrum be-hind the shock (Spitkovsky 2008b; Martins et al. 2009b;Sironi & Spitkovsky 2009).Since most of the magnetic power generated by Weibelinstability is concentrated on scales as small as a fewplasma skin depths, it has been speculated that the emis-sion mechanism in unmagnetized collisionless shocks maybe the so-called “jitter” radiation. Whereas the standardsynchrotron emission applies to large-scale fields, the jit-ter regime is realized if the scale λ of the turbulence issuch that the wiggler parameter K ≡ qBλ/mc (cid:28) Electronic address: [email protected];[email protected] harder than expected from synchrotron radiation (e.g.,Preece 1998). Since PIC simulations self-consistentlyprovide both the strength and the spatial structure ofelectromagnetic fields, as well as the particle distribu-tion, it is possible to calculate the photon spectrum fromfirst principles, thus determining whether the emission ismore synchrotron-like or jitter-like.In this work, we present synthetic spectra extractedfrom PIC simulations of relativistic unmagnetized col-lisionless shocks. In § § § § SHOCK STRUCTURE AND PARTICLE ENERGYSPECTRUM
We use the three-dimensional (3D) electromagneticPIC code TRISTAN-MP (Buneman 1993; Spitkovsky2005) to simulate a relativistic shock propagating into anunmagnetized pair plasma. The shock is triggered by re-flecting an incoming cold “upstream” flow off a conduct-ing wall at x = 0 (e.g., Sironi & Spitkovsky 2009). Thesimulation is performed in the “wall” or “downstream”frame. The incoming flow propagates along − ˆ x withLorentz factor γ = 15, and the shock moves along + ˆ x .To follow the shock evolution for longer times withfixed computational resources, we use a 2D simulationdomain in the xy plane. In the case of an unmagne-tized 2D shock, only the in-plane components of thevelocity, current and electric field, and only the out-of-plane component of the magnetic field are present.Each computational cell is initialized with 16 particlesper species. The relativistic plasma skin depth for theupstream flow ( c/ω p ) is resolved with 10 cells and thesimulation timestep is ∆ t = 0 . ω − . The simulation a r X i v : . [ a s t r o - ph . H E ] A ug Fig. 1.—
Shock structure and particle energy spectra at time ω p t = 2250. a) Number density in the simulation plane, normal-ized to the upstream density. b) Transversely-averaged density.c)-d) Magnetic and electric energy density in the simulation plane,normalized to the upstream kinetic energy density. e) Transversely-averaged magnetic (black) and electric (red) energy density. f) Lon-gitudinal phase space for positrons, shown as a 2D histogram (elec-tron phase space is nearly identical). g) Transverse phase space forpositrons. h) Positron spectrum p d N/ d p versus p (thick solid line; p ≡ γβ = [( γβ x ) + ( γβ y ) ] / is the particle 4-velocity) in theregion x sh − c/ω p < x < x sh + 400 c/ω p (delimited by the ar-rows at the bottom of panel (g)), and the relative contributionsof particles behind (thin solid line) and ahead (dotted line) of theshock. The dashed line corresponds to a power-law distributiond N/ d p ∝ p − α with α = 2 .
5. i) Positron spectrum p x d N/ d p x ver-sus p x (with p x ≡ γβ x ) for particles with p x > p y d N/ d p y versus p y (with p y ≡ γβ y ) for par-ticles with p y > box is 100 c/ω p wide (along y ) and, at the final simula-tion time ω p t = 4500, it is 4500 c/ω p long (along x ).In Fig. 1 we show the internal structure of the shockat time ω p t = 2250. The 2D plots of number density(Fig. 1a), magnetic energy (Fig. 1c) and electric energy(Fig. 1d) show the filaments produced by Weibel insta-bility in the upstream region, with characteristic trans-verse scale ∼ c/ω p . The magnetic filaments are ad-vected with the upstream flow, so in the simulation frametheir magnetic and electric components are comparable(black and red line in Fig. 1e, respectively). At theshock ( x sh (cid:39) c/ω p at time ω p t = 2250), the fila-ments merge and the magnetic energy density peaks at ∼
10% of the upstream kinetic energy density (Fig. 1e).The magnetic field decays farther downstream, where thefield is confined within islands of typical scale ∼ c/ω p (Fig. 1c).The particle energy spectrum behind the shock ( x sh − c/ω p < x < x sh ; thin solid line in Fig. 1h) consists ofa relativistic Maxwellian and a high-energy tail, whichcan be fitted as a power-law of index α = 2 . x sh < x < x sh +400 c/ω p ; dotted line in Fig. 1h), the unshocked beampopulates the low-energy peak at p ≡ γβ (cid:39)
15, whereasthe shock-accelerated returning particles with γβ x > x sh − c/ω p < x < x sh +400 c/ω p ; thick solid line in Fig. 1h)is significantly flatter than in the downstream spectrum. SYNTHETIC PHOTON SPECTRUM
Numerical Technique
We summarize the method introduced by Hededal(2005) and Hededal & Nordlund (2005) (see alsoNishikawa 2009; Martins et al. 2009a) to extract syn-thetic spectra from simulations of collisionless shocks.The electric far-field from a particle with charge q , ve-locity v = β c and acceleration ˙ v = ˙ β c is (Jackson 1999) E ( x , t ) = qc (cid:34) n × { ( n − β ) × ˙ β } (1 − β · n ) R (cid:35) ret , (1)where the unit vector n points toward the observer, atdistance R from the emitting particle. Here, the quan-tity in square brackets is to be evaluated at the retardedtime t (cid:48) = t − R ( t (cid:48) ) /c . The photon spectrum is then com-puted via the Fourier transform of Poynting flux associ-ated with the field in eq. (1). The energy d W receivedper unit solid angle dΩ (around the direction n ) and perunit frequency d ω can be computed as (Jackson 1999)d W dΩd ω = q π c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) + ∞−∞ n × { ( n − β ) × ˙ β } (1 − β · n ) e iω ( t (cid:48) − n · r ( t (cid:48) ) /c ) d t (cid:48) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (2)where r ( t (cid:48) ) is the particle trajectory. Here, we neglectedthe Tsytovich-Razin effect due to the dispersive proper-ties of the plasma (Rybicki & Lightman 1979).In PIC simulations we know the positions, velocitiesand accelerations of simulation particles with time reso-lution ∆ t = 0 . ω − . In order to accurately computethe integral in eq. (2), we interpolate the orbit of the se-lected particles so that to achieve an effective timestepof 0 . t . For a given choice of n , we can then integrateeq. (2) to obtain the photon spectrum from each parti-cle. Assuming that the far-fields by different particles arephase-uncorrelated, the total spectrum will be the sumof the spectra of individual particles.We have implemented eq. (2) and tested it for the casesof synchrotron, bremsstrahlung and wiggler/undulatorradiation, finding good agreement with analytic solu-tions. Following Hededal (2005) and Hededal & Nord- Fig. 2.—
Solid lines: photon spectrum from a power-law popu-lation of particles (d N/ d p ∝ p − α between p min = 50 and p max =2000, with α = 2 .
5) injected behind the shock ( x sh − c/ω p 000 particles moving in the electromagnetic fields self-consistently produced in our PIC simulations. Wetypically follow the particle trajectories over 3000 ∆ t =135 ω − , which is long enough to reach convergence in theshape of the spectrum, but short compared to the char-acteristic time of the shock evolution, such that the spec-trum we obtain may be regarded as “instantaneous.” Wecompute the spectrum for both head-on emission ( n = ˆ x ,which we call “ n x = 1” from now on) and edge-on emis-sion ( n = ˆ y , “ n y = 1” from now on). Our calculationsare performed in the downstream fluid frame. An ad-ditional Lorentz transformation is required if the down-stream medium is moving with respect to the observer. Results We have studied the emission from relativistic shockswith two experiments. First, we have injected a power-law distribution of particles into the downstream region,and we have traced their orbits in a fixed snapshot of thefields of the simulation. The photon spectrum is thencalculated with the technique described above. Second,we compute the spectrum from trajectories of particlesextracted directly from the simulation. In this case, boththe particle distribution and the electromagnetic turbu-lence are provided by the PIC simulation, and the result-ing emission will be the self-consistent spectrum from ourrelativistic shock.Fig. 2 shows the emission spectrum from a power-law population of particles behind the shock ( x sh − c/ω p < x < x sh ). We inject a 2D isotropic distri- bution with d N/ d p ∝ p − α ( α = 2 . p min = 50and p max = 2000. The spectral index and the lower cutoffof the distribution are chosen to mimic the high-energytail of the downstream particle spectrum in the simula-tion (thin solid line in Fig. 1h). The injected particles areevolved in a fixed snapshot of the electromagnetic fieldsfrom the PIC simulation, at time ω p t = 2250.Since electric fields are negligible in the downstreammedium (see Fig. 1d-e), the resulting photon spectrumwill probe the strength and structure of the magneticfields, and it will clarify which regime – synchrotron orjitter – is appropriate to describe the particle emission.The solid lines in Fig. 2 show that the spectrum can bewell approximated by two power-law segments. Regard-less of the observer’s direction n (red line for n x = 1,blue line for n y = 1), the slope at the low frequenciesis remarkably close to 2 / In N dimensions, thehigh-frequency slope should be − [ α − (4 − N )] / 2, whichreduces to − ( α − / − . 25 for N = 2 and α = 2 . 5, inagreement with our spectra (see dotted line). The sim-ilarity between the cases n x = 1 and n y = 1 suggests,given the isotropy of the injected particle distribution,that the downstream magnetic fluctuations are spatiallyisotropic, as seen in Fig. 1c. A transition to the jitter regime should appear whenthe wiggler parameter K becomes significantly smallerthan unity (Medvedev 2000, 2006; Fleishman 2006a,b).We have tried to artificially lower the value of K by de-creasing the strength of the electromagnetic fields, by afactor of 10 (dot-dashed lines in Fig. 2) and 100 (dashedlines in Fig. 2). The high-frequency power-law decreasesin intensity and shifts to lower frequencies, proportion-ally to the average magnetic field. The low-frequencyspectrum becomes softer for decreasing B , approachingthe flat slope expected in the jitter regime when theshock is viewed edge-on. When we impose a magneticfield spectrum of the form B ( k ) ∝ δ (ˆ k − ˆ k ) (here, ˆ k is a fixed direction in k -space) with wiggler parameter K (cid:28) 1, we are able to recover the hard low-frequencyslope ( ∝ ω ) discussed by Medvedev (2000) for head-onemission from shocks. However, we do not observe it forthe downstream turbulence self-consistently generated inthe simulation, suggesting that the magnetic field fluctu-ations are not sufficiently ordered.In Fig. 3 we show the photon spectrum resulting froma sample of particles extracted directly from the PICsimulation, followed near ω p t = 2250 in the time-varying electromagnetic fields of the simulation. The selectedparticles start in the region x sh − c/ω p < x < x sh +400. The photon spectrum (thick solid lines in Fig. 3; redfor n x = 1, blue for n y = 1) confirms that the emissionoccurs in the synchrotron regime, as the 2 / p ≤ 50 (dot-dashed lines in Fig. 3a). At high frequen- For a 3D distribution the low-frequency spectrum is ∝ ω / . The slight difference at high frequencies between n x = 1 and n y = 1 is due to the residual electric fields at x (cid:46) x sh (see Fig. 1d-e), and it disappears if electric fields are neglected while computingthe spectrum. Fig. 3.— Photon spectrum (thick solid lines, in all panels) for par-ticles extracted from the PIC simulation, evolved in time-varyingelectromagnetic fields around ω p t = 2250. Red lines are for head-on emission ( n = ˆ x ), blue lines for edge-on emission ( n = ˆ y ). a)Dot-dashed lines only include the contribution of thermal particles,with 4-velocity p ≤ 50. b) Relative contribution of downstream(thin solid lines) and upstream (dotted lines) particles. c) Timeevolution of the photon spectrum, from ω p t = 2250 (thick solidlines) to ω p t = 4500 (dot-dashed lines). Here, the case n x = 1 hasbeen shifted downward by a factor of 2 for clarity. cies, the emission along ˆ y (blue solid line in Fig. 3a) ismore powerful than along ˆ x (red solid line in Fig. 3a).This reflects the fact that the highest energy particlesare grazing the shock surface (Spitkovsky 2008b) andthey mostly contribute to the emission along ˆ y . Thedifference between n x = 1 and n y = 1 in Fig. 3a canbe directly related to the difference between the particlespectra p x d N/ d p x and p y d N/ d p y shown in Fig. 1i (thicksolid lines, respectively red and blue). In fact, due to rel-ativistic beaming, an observer located along ˆ x will morelikely detect the radiation from particles with p x /p y (cid:29) p x > p y /p x (cid:29) p y > 0) will contribute more to the emission along ˆ y .As shown in Fig. 3b, the spectrum of downstreamparticles (thin solid lines) dominates the total emission(thick solid lines) in the low-frequency bump, whereasthe high frequencies are mostly powered by the spec- trum of upstream particles (dotted lines). This is due toa combination of two effects: first, at the highest particleenergies the upstream particles may dominate by num-ber (compare red dotted and thin solid lines in Fig. 1i);second, the electromagnetic energy in the upstream re-gion is on average larger than in the downstream, sinceupstream electric fields are as strong as magnetic fields,and also because the field is built up ahead of the shockon a length scale larger than the decay scale in the down-stream (Fig. 1e). It is worth pointing out that, at low fre-quencies, the spectra of both downstream and upstreamparticles can be fitted with the synchrotron slope 2/3(black dashed lines).Fig. 3c shows the time evolution of the total pho-ton spectrum, from ω p t = 2250 (thick solid lines) to ω p t = 4500 (dot-dashed lines). The upper cutoff inthe particle energy spectrum grows linearly with time(Spitkovsky 2008b; Sironi & Spitkovsky 2009), and as aresult the photon spectrum extends to higher frequen-cies. Also, the spectral intensity increases, especially athigh frequencies, due to stronger electromagnetic fieldsand a larger number of accelerated particles. DISCUSSION We have computed synthetic photon spectra from 2DPIC simulations of relativistic collisionless shocks byfollowing a sample of simulation particles in the time-varying fields of the simulation. The low-frequency partof the spectrum scales as ∝ ω / , as expected for syn-chrotron emission from a 2D particle distribution. Al-though the electromagnetic fluctuations generated byWeibel instability are on small (skin-depth) scales, theparticle emission does not occur in the jitter regime.In retrospect, this is not surprising. The characteristiclength scale of the magnetic turbulence is λ (cid:38) c/ω p ,and in the shock region, where most of the emissionis produced, the magnetic energy reaches a fraction (cid:15) B (cid:39) . r L / ( c/ω p ) = (cid:15) − / B (cid:39) 3, where r L is the relativisticLarmor radius of a particle moving with the upstreamflow. It follows that K = λ/ ( r L /γ ) (cid:39) γ , so thatthe condition K (cid:28) γ = 3shock is still consistent with synchrotron radiation. Forelectron-ion shocks, the wiggler parameter will be m p /m e times larger, and we are even deeper in the synchrotronregime. Further downstream from the shock, althoughthe magnetic field strength decays, the value of K doesnot significantly change, since short-wavelength modesare progressively damped (Chang et al. 2008), and thecharacteristic scale of the turbulence increases.Although the results presented here apply to a 2D par-ticle distribution, the main conclusions should hold for3D configurations as well. There, we expect the low-frequency slope to be 1 / 3, as in the standard synchrotronradiation. If the GRB emission results from high-energyparticles accelerated in relativistic unmagnetized shocks,it seems that resorting to the jitter radiation is not a vi-able solution for the “line of death” puzzle (Preece 1998).At high frequencies, our results suggest that the contri-bution of upstream particles to the total emission, whichis usually omitted in standard models, is not negligible.It causes the radiation spectrum to be flatter than thecorresponding downstream spectrum, thus partly mask-ing the contribution of downstream thermal particles(Giannios & Spitkovsky 2009). This could potentiallyexplain the absence of clear signatures of downstreamthermal emission in GRB shocks (Band 1993). Simu-lations extending to longer times (and higher particleenergies) will help to clarify these issues.We remark that our calculations do not include ra-diative particle cooling, synchrotron self-absorption and inverse Compton radiation. Still, we have shown thatthe computation of synthetic spectra from self-consistentPIC simulations provides a powerful tool for studyingthe origin of astrophysical non-thermal emission.We thank J. Arons, A. Celotti, J. Kirk, S. Martinsand L. Silva for comments and suggestions, and KITPSanta Barbara for hospitality. This research was sup-ported by NSF grants AST- 0807381 and PHY05-51164.3, as in the standard synchrotronradiation. If the GRB emission results from high-energyparticles accelerated in relativistic unmagnetized shocks,it seems that resorting to the jitter radiation is not a vi-able solution for the “line of death” puzzle (Preece 1998).At high frequencies, our results suggest that the contri-bution of upstream particles to the total emission, whichis usually omitted in standard models, is not negligible.It causes the radiation spectrum to be flatter than thecorresponding downstream spectrum, thus partly mask-ing the contribution of downstream thermal particles(Giannios & Spitkovsky 2009). This could potentiallyexplain the absence of clear signatures of downstreamthermal emission in GRB shocks (Band 1993). Simu-lations extending to longer times (and higher particleenergies) will help to clarify these issues.We remark that our calculations do not include ra-diative particle cooling, synchrotron self-absorption and inverse Compton radiation. Still, we have shown thatthe computation of synthetic spectra from self-consistentPIC simulations provides a powerful tool for studyingthe origin of astrophysical non-thermal emission.We thank J. Arons, A. Celotti, J. Kirk, S. Martinsand L. Silva for comments and suggestions, and KITPSanta Barbara for hospitality. This research was sup-ported by NSF grants AST- 0807381 and PHY05-51164.