Magnetic Field Structure and Faraday Rotation of the Plerionic Supernova Remnant G21.5–0.9

The double-lobed structure in G21.5−0.9 is also observed in other radio PWNe, including Vela (Dodson et al. 2003), G343.1−2.3 (Romani et al. 2005; Y. H. Liu et al. 2022, in preparation), and G76.9+1.0 (Landecker et al. 1993; Arzoumanian et al. 2011; Arumugasamy et al. 2014). 3C58 (Reich 2002; Bietenholz 2006) and DA 495 (Kothes et al. 2008) show only a double-lobed structure in the polarized emission, but not in the total intensity. The positions of the lobes are always symmetric about the pulsar spin axis, except for DA 495, whose pulsar spin orientation is unknown. The symmetry suggests that the X-ray torus might be related to the lobes. The observed double-lobed morphology of the Vela and G76.9+1.0 PWNe was suggested to result from the toroidal emission region, with specific viewing geometry (Chevalier & Reynolds 2011). This model, however, may not be applicable to spherical PWN like G21.5−0.9.

All these lobes are only seen in radio, suggesting the accumulation of less energetic particles in these regions. It is peculiar that they all show a double-lobed morphology, but not a ring-like one, as in the X-ray torus. In Table 2, we summarize the properties of four double-lobed radio PWNe. We notice that the two lobes always show different brightness. Vela has the largest peak brightness difference between the two lobes of a factor of 2. There is also a large range in the physical separation of the lobes, from 0.33 pc in G21.5−0.9 to 8.7 pc in G76.9+1.0. From the table, there seems to be a trend that older PWNe have larger lobe separation. This trend also exists for the ratio of the separation to the size of the X-ray torus. If confirmed, this could imply an outward motion of the lobes.

Table 2. Double-lobed Radio PWNe

	G21.5−0.9	Vela	G343.1−2.3	G76.9+1.0
Age (kyr)	0.87	11	18	9 ∼ 40
Separation between two lobes (pc) (llobes)	0.36	0.36	1.24	8.7
Peak brightness ratio between two lobes	1.07	2	1.14	1.08
Size of X-ray torus (pc) (ltorus)	0.26	0.06	0.1	⋯
llobes/ltorus	1.4	6	12	⋯
${\dot{E}}_{\mathrm{pulsar}}\,({10}^{36}\,\mathrm{erg}\,{{\rm{s}}}^{-1})$	34	6.9	3.4	30

References. The X-ray torus sizes are from Ng & Romani (2008).

Download table as:  ASCIITypeset image

It is intriguing that the RM toward G21.5−0.9 appears to show a symmetric pattern about the pulsar spin axis (Figure 4). Instead of the chance alignment of the foreground RM structure, this more likely results from the Faraday rotation within the PWN. This is supported by the amplitude of the RM variation. Minter & Spangler (1996) derived a structure function for the foreground RM, D_RM(δ θ), in terms of angular scale, due to the fluctuation of both the interstellar electron density and the galactic magnetic field caused by turbulence (see the Appendix for details). The theoretical structure function predicts the average RM difference for a given angular separation. For the case of G21.5−0.9, we obtain ${D}_{\mathrm{RM}}(\delta \theta )\approx (13{B}_{\mathrm{ran}}^{2}+11)\delta {\theta }^{5/3}$ rad² m⁻⁴, where B_ran is the strength of the random B-field component in units of μG and δ θ is the angular separation in units of arcmin. The value of B_ran depends on the observation direction, and its largest empirical value is 9 μG (Haverkorn et al. 2008). This model can be compared with our observations, and the result is shown in Figure 6. The structure function is calculated up to an angular scale of 45'', the same as the PWN size. The observed structure function is larger than the theoretical values at all δ θ for B_ran ≤ 9 μG. This suggests that the foreground Faraday rotation is unable to explain the rapid variation of RM in G21.5−0.9 unless B_ran is exceptionally high. This result indicates that the observed RM variation is unlikely to be due to the foreground.

Zoom In Zoom Out Reset image size

Pulsar wind in general cannot contribute significant Faraday rotation. The main reason is that it consists of electron–positron plasma, such that the Faraday effect is canceled out. Moreover, Faraday rotation is weak for relativistic particles (Jones & Odell 1977), and their number density in pulsar wind is low; for example, less than 10⁻⁵ cm⁻³ for G21.5−0.9 (Hattori et al. 2020). On the other hand, the RM fluctuation we observed (∼80 rad m⁻²) can be explained by electrons originating from the SN ejecta. For a magnetic field strength of 130 μG (Guest et al. 2019) and a nebula size of 2 pc, the observed RM value requires an average electron density of ∼0.8 cm⁻³. We have doubled the density to account for the effect of the internal Faraday rotation (Burn 1966). Assuming that G21.5−0.9 is a sphere, the total number of cold electrons can be contributed by 0.1 M_⊙ of hydrogen or 0.2 M_⊙ of high-Z elements. This is possible for G21.5−0.9, since the swept-up ejecta mass is suggested to be 0.8 M_⊙ (Hattori et al. 2020).

In this picture, the spatial variation of the RM is then caused by the variation of either the line-of-sight magnetic field strength or ejecta electron density. The radial B-field of G21.5−0.9 is expected to produce a spherically symmetric RM pattern rather than an axisymmetric one. Also, nonuniform field strength should result in a correlation between the RM and the surface brightness, which is not observed. We therefore conclude that there is an axisymmetric distribution of ejecta electrons about the pulsar spin axis. This could be due to an asymmetric SN explosion or a latitude-dependent mixing process between the pulsar wind and the SN ejecta. Further simulations are needed to understand the details of these processes.

The RM distributions of the Boomerang PWN (Kothes et al. 2006) and CTB 87 (Guest et al. 2020; Kothes et al. 2020) show similar patterns as G21.5−0.9. They are both symmetric about the pulsar spin axis, and the RM values are lower or more negative along the axis. They have also been suggested to be caused by the internal Faraday rotation (Kothes et al. 2006, 2020). Unlike G21.5−0.9, both of them are evolved PWNe that have been crushed by the SN reverse shock. This could imply that the axisymmetric RM pattern already existed before the crush and was not disturbed by the crush. However, we note that these two sources have different B-field structures than G21.5−0.9. The origins of their RM patterns thus may not be the same as for our case.

Our high-resolution polarization map confirms that G21.5−0.9 has a radial magnetic field overall. A previous infrared study found a toroidal field in the inner PWN with scale of 4'' (i.e., 0.1 pc; Zajczyk et al. 2012). If this extends to larger radii, its influence on the inner PWN may explain the hourglass structure seen in the PL map (see Section 3). It is, however, unclear what the extent of the toroidal field is and how it connects to the global radial field of the radio PWN. Nevertheless, as we will show in Section 5 below, we can reject the existence of any large-scale toroidal field components with comparable strength. This suggests that the toroidal field is only important to the inner nebula.

In addition to G21.5−0.9, the Crab Nebula and DA 495 also show radial magnetic field structures. For the Crab, the radial field is only found near the boundary, and it could be caused by fluid instability between the pulsar wind and the ambient medium (Bietenholz & Kronberg 1990). While numerical simulations suggest that the instability layer can be as thick as 30% of the PWN radius (Bucciantini et al. 2004), this is still insufficient to explain our observations. The radial field of DA 495 has been suggested to be a combination of dipole and toroidal fields (Kothes et al. 2008), which is not the case for G21.5−0.9 (see Section 5 below).

Finally, we note that as shown in Figure 3, the PF is very low near the center and gradually increases toward the edge. This cannot be solely due to beam depolarization, since this effect becomes negligible for angular scales beyond the beam size. As an example, consider a radial B-field emitting at 100% PF locally: the PF is zero at the center, because of beam depolarization, but rises to 85% just one beamwidth outside. In Section 5, we suggest that the radially increasing PF is caused by magnetic field geometry and turbulence. The low PF near the center could also lead to the observed offset between the peaks of the total intensity emission and the polarized emission.

Three-dimensional (3D) simulations have shown that the magnetic field in a PWN could have significant turbulence (Porth et al. 2014). In our model, we consider ordered and disordered (i.e., turbulent) B-field components, and adopt the method of Bandiera & Petruk (2016) to calculate synchrotron emissivity in the presence of small-scale magnetic turbulence. On top of the ordered magnetic field B , a Gaussian random field with variance (∣ B ∣σ)² in each direction is introduced to account for a turbulent environment. The energy ratio between the ordered and disordered fields is thus ${B}_{\mathrm{ordered}}^{2}/{B}_{\mathrm{disordered}}^{2}=1/3{\sigma }^{2}$. A larger σ would decrease the local PF. For simplicity, we assume a constant σ = 0.7 over the entire nebula. The value was chosen so that the PF values are similar to the observations. This gives B_ordered/B_disordered ∼ 0.8, which is close to the previous estimate (Furst et al. 1988).

We considered three models for the ordered magnetic field. The first one is a radial B-field in 3D (the R model) motivated by our polarization map. The second one is a toroidal field plus a radial field (the TR model), because the former was observed in the inner PWN (Zajczyk et al. 2012). The third model is a toroidal field plus a dipole field (the TD model), which was used to explain the observations of DA 495 (Kothes et al. 2008). Mathematically, these can be expressed in spherical coordinates (r, θ, ϕ) as

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{{\boldsymbol{R}}}(r)=-{B}_{0}\hat{r},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{{\boldsymbol{TR}}}(r,\theta )={B}_{0}(\sin \theta \hat{\phi }-\hat{r}),\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{{\boldsymbol{TD}}}(r,\theta )={B}_{0}\left(\cos \theta \hat{r}+\displaystyle \frac{1}{2}\sin \theta \hat{\theta }+\sin \theta \hat{\phi }\right),\end{eqnarray} \tag{ 3 }$$

where B₀ is a constant to characterize the field strength. The form of the toroidal component in the TR and TD models is adopted from Del Zanna et al. (2006). For simplicity, we consider equal strengths for the radial, toroidal, and dipolar field components. The pulsar spin axis points to the + z direction, and θ is the polar angle measured from the spin axis.

We note that the directions of the radial field (inward or outward) and the toroidal field (clockwise or counterclockwise) are not constrained by the surface brightness map because they give the same synchrotron emissivity. However, they induce different signs of internal Faraday rotation. An outward radial field produces a positive RM and an inward field produces a negative RM. The average RM of G21.5−0.9 is ∼65 rad m⁻² near the edge and ∼50 rad m⁻² at the center. This suggests a negative RM contribution by the nebula, hence an inward radial field. The direction of the toroidal field is randomly picked. Except for the magnetic field configuration, the other settings of the models are the same.

The viewing angle between the pulsar spin axis and the line of sight for G21.5−0.9 is 854 (Ng & Romani 2008). We take it as 90° in the modeling. We consider only the region from the TS (R_TS = 0.13 pc) to the outer radius (R_pwn = 1.14 pc), and we assume no emission outside this region.

To compute the synchrotron emissivity, we also need to specify the pulsar wind particle distribution. We adopt the diffusion model deduced from the X-ray photon index and TeV brightness profiles (Tang & Chevalier 2012; Abeysekara et al. 2017). For injected particles following a power-law energy distribution, N(E, r = 0) = KE^−(2α+1), there is an exact formula for the density profile (see Equation (17) in Gratton 1972). When $E\ll 422/({B}_{0}^{2}t)$ and $r\ll \sqrt{{Dt}}$ (in cgs units), this can be approximated by

$$\begin{eqnarray}&&{n}_{\mathrm{pwn}}(r)\propto \displaystyle \frac{K}{{Dr}}\ \ {\rm{for}}\ \ {R}_{\mathrm{TS}}\lt r\lt {R}_{\mathrm{pwn}},\end{eqnarray} \tag{ 4 }$$

where D is the diffusion coefficient and t is the age of the PWN. For our case, the first condition is easily satisfied, given that we are dealing with radio-emitting electrons. The diffusion coefficient D is found to be 3.7 × 10²⁷ cm² s⁻¹ by fitting the X-ray data of G21.5−0.9 (Tang & Chevalier 2012). Although D could be smaller for less relativistic particles, the difference is subtle. Given that the PWN age is 870 yr and the radius is around 1 pc, the second limit is approximately fulfilled. We do not attempt to fit B₀ and K, but arbitrarily set their values and only focus on the relative brightness. The spectral index α is taken to be + 0.16, derived from our observations.

We followed the recipe derived by Bandiera & Petruk (2016) for deriving synchrotron emission in the Stokes I, Q, and U parameters. The final maps in each Stokes parameter are obtained by integrating the emission along each line of sight. We simulate the emission at 5 GHz, the same as the observations. Absorption and scattering are negligible at this frequency. Finally, we smooth the synthetic emission maps by the same beam size as the observations to allow a direct comparison. As we will show in Section 5.3 below, depolarization caused by internal Faraday rotation is negligible at this frequency. We therefore did not account for this effect in our modeling.

Figure 7 shows the simulated surface brightness and polarization maps of different magnetic field configurations. When comparing with the observations, we focus on two prominent features: the radial magnetic field structure and the radially increasing PF. The TR model results in an overall radial field structure. However, there are two local minima in the PF map along the polar axis, rather than a minimum at the center. The low-PF regions are where the toroidal and radial fields are perpendicular to each other. The TD model barely reproduces a radial field structure, and the PF map even has four local minima, for the same reason as above. The R model successfully produces a radial field pattern and radially increasing PF.

Zoom In Zoom Out Reset image size

In Figure 8, we compare the simulated flux density and PF profiles with the observations. It shows that the relative surface brightness is insensitive to the magnetic field structure. The general downward trend is caused by the spherical geometry and the 1/r dependence of the pulsar wind particle density. The differences between the B-field configurations are not significant, and all are roughly consistent with the observations. Note that our model results are very similar to one with constant particle density and the magnetic field strength decreasing as 1/r (Furst et al. 1988). This is due to the degeneracy of these two parameters for synchrotron radiation.

Zoom In Zoom Out Reset image size

Unlike the brightness maps, different magnetic field configurations give very different PF maps. Only the R model shows an increasing trend with radius that generally matches the observations. We argued in Section 4.3 that the angular size of this trend is too large to be attributed to beam depolarization. The increasing trend can be understood as follows. First, considering the case without turbulence (i.e., σ = 0), if we ignore the Faraday effect, the observed polarization direction is perpendicular to the orientation of the B-field projected onto the plane of the sky. For the R model, picking any line of sight, the projected B-field always has the same position angle. Therefore, after projection, the polarized emission from all cells in 3D adds up coherently, such that the observed PF is constant over the entire nebula and equal to the intrinsic PF of each cell. Note that this is not the case for the TR and TD models. They are expected to have lower PFs in the inner nebula, due to the changing magnetic field PA along the line of sight.

Now, considering turbulence, an analogous formula of the local PF in each cell can be written as ${{PF}}_{\mathrm{local}}\sim \tfrac{{PI}\sin \eta }{I(\sin \eta +\sigma )}$, where η is the angle between the line of sight and the magnetic field in 3D. A large turbulence level σ thus leads to a lower PF, and the latter also depends on η. In the R model, the edge of the nebula has η ≈ 90°, and hence the highest PF, while the PF is lowest near the center as η ≈ 0. The same result can be derived rigorously using a full mathematical treatment with the synchrotron recipe (Bandiera & Petruk 2016).

We conclude that the R model is sufficient for explaining the observed radial profiles of brightness and PF. This rules out the existence of a global toroidal or dipolar B-field components with comparable strength to the radial field.

In the following discussion, we focus on the R model to discuss its effect on the internal RM and depolarization. We model the ejecta electron density distribution as ${n}_{\mathrm{sn}}(\theta )={n}_{\mathrm{sn},0}| \cos \theta |$, since it fits our observational results, where ${n}_{\mathrm{sn},0}$ is a constant and θ is the polar angle measured from the spin axis of the pulsar. These electrons are cold, such that their contribution to the synchrotron emission is negligible. We find that for B₀ = 130 μG (Guest et al. 2019), ${n}_{\mathrm{sn},0}=4.2$ cm⁻³ can give an RM fluctuation of ΔRM = 80 rad m⁻² within the PWN. This implies a volume-averaged number density of $\sim 0.5{n}_{\mathrm{sn},0}=2.1\,{\mathrm{cm}}^{-3}$.

Figure 9 shows the simulated RM maps. The internal RM values are all nonpositive because the magnetic field in the near hemisphere points away from the observer. There is no Faraday rotation at the equator, since we assumed ${n}_{\mathrm{sn}}(\theta =\pi /2)=0$ in the model. For a better comparison with the data, we added a constant foreground RM of 105 rad m⁻². The resulting RM map is shown in the middle panel of Figure 9. The RM values are smaller along the pulsar spin axis, which generally matches with the observations. At the PWN boundary, our model suggests minimal internal Faraday rotation due to the short light path. This gives a nearly constant RM value around the boundary, in contrast to what we have observed. The discrepancy could be an observational bias, since some pixels near the edge are clipped in the map due to the low surface brightness.

Zoom In Zoom Out Reset image size

Internal Faraday rotation can lead to depolarization and break the linear relation between PA and λ² at low frequency (Burn 1966). This has been suggested to be the cause of the frequency-dependent depolarization found in the Boomerang and CTB 87 PWN (Kothes et al. 2006, 2020). We illustrate this effect with a line of sight near the strongest internal RM (see the left panel of Figure 9). We calculated the Faraday rotation of the emission from different depths, then summed the rotated Stokes Q and U parameters respectively to determine the resulting PA and PF. This calculation was repeated for different wavelengths, and the results are plotted in Figure 10. We found that the PA varies linearly with λ² at short wavelengths, but starts to deviate significantly at frequencies ≲2 GHz. For the PF, it drops continuously as the wavelength increases, and reaches half of its maximum value at ∼2 GHz. Our result suggests that the internal Faraday effects on the PA and PF are minimal at our observation frequency, but could be significant in the L band. For instance, we cannot simply derive the intrinsic PA using a single RM value, but need to combine it with higher-frequency (e.g., C-band) data to solve for the Faraday dispersion function using RM synthesis (Brentjens & de Bruyn 2005).

Zoom In Zoom Out Reset image size