RADIATION SPECTRAL SYNTHESIS OF RELATIVISTIC FILAMENTATION

The relativistic FI is set up by preparing two equal density neutral plasma beams (Figure 1) under fully periodic boundary conditions. Beams were chosen with bulk flows of Γ_Bulk ∈ [2, 4, 6, 10, 15] (Γ_Bulk ∈ [4, 6, 10]) for 3 (2) spatial components and 3 (2) velocity components or 3D3V (2D2V) runs. Simulation volumes were 500Δx × 500Δy × 500Δz = 50³δ³_e (6000Δx × 6000Δz = 300²δ²_e), respectively. Plasma densities were adjusted such that all simulations contained an equal number of relativistic electron skin depths, δ_e = ω⁻¹_pec, independent of Γ_bulk, with ω_pe = [4πn_eq²/m_eΓ_bulk]^1/2, the relativistic plasma frequency. We simulated both baryonic (e⁻ + p⁺) and pair (e⁻ + e⁺) plasmas. A reduced mass ratio was used for the baryonic case; m_i/m_e ≡ 64 suffices in separating ion and electron dynamic timescales.

Zoom In Zoom Out Reset image size

In addition, a subset of the radiation spectra were collected in static (frozen) electromagnetic fields, starting from a well-developed time in the simulation rather than in situ full electromagnetic dynamics. These "frozen" test runs are analyzed and discussed by Medvedev et al. (2010). In summary, they conclude that unlike realistic synthetic spectra collected in situ, spectra obtained in static (yet self-consistent) fields are well (but erroneously) interpreted by a synchrotron template. This along with the obvious advantage of computational speed and ease was a main motivation for employing a completely self-consistent in situ synthesis approach.

We chose the often employed PIC scaling of natural constants, setting c ≡ m_e ≡ q_e ≡ 1. Varying n_e (number or charge density) ensured an equal number of skin depths for all Γ_Bulk for reasons of comparison. Time is henceforth normalized to tω⁻¹_pe ≡ 1. All results were obtained using a sixth-order field solver for Maxwell's equations. The field spectra were only affected qualitatively at high wavenumber (i.e., in the extreme dissipative range). Spectral synthesis was verified to exhibit negligible dependence of field solver order. We used cubic (quadratic) spline interpolation between particles and the field grid in three dimensions (two dimensions).

We incorporate into the PhotonPlasma PIC code, developed in Copenhagen (Haugbølle 2005; Hededal 2005), a discretized version of the standard expression for radiated power from accelerated charges P(ω) ∝  d²W/dΩdω (e.g., Rybicki & Lightman 1979). Using the aforementioned PIC scaling, synthetic radiation frequencies are normalized to ω_pe. Spectra are collected at runtime, and for simulations using N_tot ∼ 10¹⁰ particles, we collect spectra from N_synth ∼ 10⁶ particles without noticeable slow-down in simulation speed. We retain all phase and temporal information of the ensemble signal by integrating the signal as

$$\begin{equation*} P_\omega \propto \bigg|\int ^{t1}_{t0}\sum _j^{N_{\rm synth}}E_{{\rm ret},j}e^{i\omega \phi ^{\prime }}dt^{\prime }\bigg |^2, \end{equation*}$$

with $\phi ^{\prime } \equiv t^{\prime }-\vec{n}\cdot \vec{r_0}(t^{\prime })/c)$. Spectra are recorded every small multiple of a simulation time step. Thus, we may extract "instantaneous" spectral states as P_ω(t₁, 0) − P_ω(t₀, 0). We can specify multiple (far-field) observers at arbitrary orientation, and we can choose multiple regions for spectral collection in the simulation volume. In practice, we have a choice of selecting phase coherent or phase incoherent spectral synthesis. Binning over frequency spectral range was chosen to be logarithmic in Δω, but can in general be either logarithmic or linear in Δω; linear frequency binning is used primarily for test purposes. Particles selected for synthetic radiation contribution are subcycled, typically Δt_synth ⩽ 10⁻¹Δt_sim to resolve high-frequency radiation emitted by particles with β_em ≈ 1 and $\vec{n}\parallel \vec{\beta }$. Six observers (placed at "infinity") in three dimensions (two dimensions) are spread uniformly over angles, $\angle \lbrace \vec{z},\vec{obs}\rbrace \in [0,\pi /2]$ with respect to the streaming direction. Observer "0" designates head-on and observer "5" is edge-on (three-dimensional case).

A number of plasma effects are neglected here. They shall be addressed in future work: plasma frequency cutoff, Razin effect, and self-absorption processes. The spectral synthesis method has been sanity checked exhaustively (Hededal 2005; Hededal & Nordlund 2005).

Growth rates of both linear and nonlinear FI depend heavily on the plasma constituents (Frederiksen et al. 2004); in e⁻ + p plasmas, due to the presence of an inert species, growth rates are reduced (Tzoufras et al. 2006) compared with e⁻ + e⁺ plasmas. Emerging electromagnetic field topologies and power spectra are likewise affected. For our setup, two intervals of instability, and a brief turn-over, are identified for the relativistic FI: (1) linear, t ∈ [0, 12] (t ∈ [0, 14], for e⁻ + e⁺ case), (2) saturation/transition, t = [14, 16] (t ∈ [12, 14]), and (3) nonlinear stage, t ∈ [16,"∞"] (t ∈ [14,"∞"]). We shall use these intervals when comparing spectra from baryonic and pair plasmas. (electro)Magnetic fields develop rapidly during the linear phase. Subsequent to saturation the field spectra converge to a quasistationary configuration. A representative case is plotted in the right-most panel of Figure 1 for late times of the evolution. The perpendicular Fourier spectra of the magnetic field, $\widetilde{B}_{k_\perp }\equiv \langle \widetilde{B}_{\perp }({k_\perp })\widetilde{B}^*_{\perp }({k_\perp })\rangle$, fit a Gaussian + a high-frequency power-law segment on the upper dynamic range, with a dissipative cutoff for k → k_Ny (Nyquist scale). "Parallel" Fourier spectra of the magnetic field, $\widetilde{B}_{k_\Vert }\equiv \langle \widetilde{B}_{\perp }({k_\Vert })\widetilde{B}^*_{\perp }({k_\Vert })\rangle$, are single power laws on k_∥,0 < k_∥ < k_∥,Ny. Except for the earliest linear FI, temporal evolution of $\widetilde{B}_{\perp }$ is modest—see Figure 1 (right panel). Particle distribution functions are well represented by two-shifted Gaussians in the streaming direction and a single Gaussian perpendicular to the streaming direction, f(p_z, p_⊥, Γ_B) ∼ (exp(−(p_z − Γ_B)²) + exp(−(p_z + Γ_B)²))exp(−p²_⊥). The average Lorentz factor of the streaming particles decrease in the simulation (lab) frame during the FI. Temperatures increase with T_⊥>T_∥ for early times and T_⊥ ≈ T_∥ in the late phase. Once the FI has burned all its momentum anisotropy, the total electron ensemble has Γ_bulk ≈ 1 (no bulk flow).

Non-isotropic momentum distributions naturally result from non-isotropic particle acceleration; we expect this to be the case based on our devised toy model of "non-Fermi" acceleration (Hededal et al. 2004; Medvedev 2006b).

We therefore proceeded to compare subsets of the total particle ensemble, in momentum space, specified by the criteria:

$$\begin{equation*} \theta _{C,m}-\delta \theta _C \le \arctan ({|\vec{p_{\perp }}|^{-1}}{\vec{p_\parallel }}) \le \theta _{C,m}+\delta \theta _C, \end{equation*}$$

effectively selecting particles with momenta confined inside a narrow pitch cone of opening angle δθ_C ≡ 1/Γ_Bulk, at varying inclination, θ_C,m ≡ mπ/8, m ∈ {0, 1, ..., 4}, with respect to the streaming direction. Thus, electrons are selected which are significantly beaming their radiative contribution in the direction of the cone; they radiate efficiently for certain observer orientations, θ_C.

Figure 2 (insets 1 and 2) show how, early during the FI, electrons experience a non-isotropic acceleration with directional preference toward forward beaming of 1/Γ_bulk. This "beamed" acceleration continues to drive a small part of the total particle ensemble to high gamma factors ($\gamma _{e,\rm max}\sim 45$). Later, as Γ_bulk decreases (bi-directionally), this anisotropy spreads to include also higher inclination angles (inset 3). This could be either due to early "preference-accelerated" electrons diffusing out of their streaming directed pitch cone (inset 1→inset 4), or due to current filaments going unstable and acquiring a cross-stream component (i.e., they start bending), or both.

Zoom In Zoom Out Reset image size

Figure 3 reveals consistently higher peak frequencies, ω_peak, for increasing observer angles up to a maximum angle 0 < |θ_peak| < π/2. For higher angles, ω_peak stays approximately constant up to θ ∼ π/2. Even at late times (Figure 3, lower panel) peak frequencies ω_peak are observer dependent and apparently increase with streaming axis inclination angle. Still, even though ω_peak(|θ|>0) increases relative to the head-on (θ ≡ 0) observer at a given time, it never exceeds the initial head-on maximum peak frequency, ω_peak(θ ≡ 0), in three dimensions.

Zoom In Zoom Out Reset image size

In the two-dimensional (2D) case viewed in Figure 4 the peak frequency increases for all observers during the FI and $\omega _{{\rm peak}, \rm {2D}}>\omega _{{\rm peak}, \rm {3D}}$ for all cases (open circles of both Figures 3 and 4). This effect is an artifact of the lower-dimensional (two-dimensional) model; particles experience more frequent interactions with the ion filaments, since the "filament hit probability" is considerably smaller in three dimensions due to the reduction in filament filling factor (filaments are actually "sheets" in two dimensions, whereas they are true line-like filaments in three dimensions). Radiation spectra are therefore formed more rapidly by the plasma field structures and with better statistics in two dimensions—albeit artifactually so.

Zoom In Zoom Out Reset image size

Fitting the semi-analytical result for jitter radiation (Medvedev et al. 2010) to three-dimensional (3D) spectra (Figure 5, upper panel) shows markedly better fits than a synchrotron spectrum for the same parameters, during early stage FI, both for baryonic plasmas and pair plasmas. At later stages, also models of synchrotron fit the spectra well, but we emphasize previous findings that the underlying plasma dynamics is not synchrotron in origin Medvedev (2000, 2006a). For a spectrum (and evolution) which has been isotropically averaged over all observer directions—Figure 5, lower panel—making the same qualitative scaled fits of the semi-analytic model produces better agreement if a jitter radiation environment is assumed. Such an averaged spectrum can be seen as the immediate spectral state in the bulk of the plasma two-stream (or locally in a shock), when the simulation frame coincides with the comoving frame of the plasma, as is the case in the FI setup presented here.

Zoom In Zoom Out Reset image size

Our result further agrees quite precisely with the scaling relation of peak frequencies for a semi-analytical jitter radiation model, reported in Medvedev et al. (2010). At early times, we have measured $\omega _{{\rm peak},\rm {3D}}\propto \Gamma ^2_{\rm bulk}$, while at late times $\omega _{{\rm peak}, \rm {3D}}\propto \Gamma ^3_{\rm bulk}$ joint by a smooth transition.

A "by-eye" re-alignment fit (thick gray line, Figure 4) of the same semi-analytic result to the 2D baryonic case shows an underestimatation of high-frequency radiated power. This undershoot artifact is not seen in three dimensions (upper green, solid, curve, Figure 5). It, too, is probably a consequence of the reduced dimensionality again by the above argument concerning 2D result artifacts.

Moreover, in two dimensions, the maximum value corresponding to the abscissae of Figure 2, $\gamma ^2_{{\rm max}, \rm {3D}}\approx 2100$, is $\gamma ^2_{{\rm max}, \rm {2D}}\approx 4200$ which is higher by a factor of $\gamma ^2_{{\rm max}, \rm {2D}}/\gamma ^2_{{\rm max}, \rm {3D}}\sim 2$. This further influences late time spectra (lower panels of Figures 3 and 5). The 2D and 3D cases exhibit gross discrepancy in peak frequencies; they differ by an order of magnitude ($\omega _{{\rm peak},\rm {2D}}/\omega _{{\rm peak},\rm {3D}}\approx 10$), due to a general upward drift of ω_peak in two dimensions. Again, 2D spectra are pushed to higher frequencies, likely due to an enhanced "filament hit probability." The 2D results for e⁺ + e⁻ plasmas exhibit analogous features.