High-Resolution Radio Study of the Dragonfly Nebula

The Dragonfly Nebula (G75.2$+$0.1) powered by the young pulsar J2021$+$3651 is a rare pulsar wind nebula (PWN) that shows double tori and polar jets enclosed by a bow-shock structure in X-rays. We present new radio observations of this source taken with the Very Large Array (VLA) at 6 GHz. The radio PWN has an overall size about two times as large as the X-ray counterpart, consisting of a bright main body region in the southwest, a narrow and fainter bridge region in the northeast, and a dark gap in between. The nebula shows a radio spectrum much softer than that of a typical PWN. This could be resulting from compression by the ram pressure as the system travels mildly supersonically in the interstellar medium (ISM). Our polarization maps reveal a highly ordered and complex $B$-field structure. This can be explained by a toroidal field distorted by the pulsar motion.


INTRODUCTION
A pulsar has a strong electromagnetic field that can rip particles off the stellar surface and accelerate them to high energies. This drives a magnetized wind that are composed of mainly electrons and positrons. As the wind streams into the surrounding medium, the particles are shocked and then are further accelerated to ultra-relativistic energies. This results in a shocked wind bubble known as a pulsar wind nebula (PWN) (see Gaensler & Slane 2006, for a review). These systems emit synchrotron radiation from radio to X-rays and inverse Compton (IC) scattering radiation from MeV to TeV gamma-rays. They serve as good nearby laboratories for studying relativistic shock physics and acceleration mechanism of Galactic cosmic rays. They are also suggested to be a main contributor to the detected positron excess in the Galaxy (Blasi & Amato 2011;Della Torre et al. 2015).
Arcsecond resolution X-ray observations revealed small-scale structures in PWNe (Kargaltsev & Pavlov 2008;Kargaltsev et al. 2012) including axisymmetric tori and jets in subsonic moving PWNe (see Ng & Romani 2004, 2008. These are resulting from toroidal magnetic field and anisotropic (latitude-dependent) energy outflow as numerical simulations suggested (e.g., Olmi et al. 2016;Porth et al. 2017). A pulsar inside a supernova remnant (SNR) generally moves subsonically or transonically due to the high sound speed (hundreds of km s −1 ) in the hot environment. Therefore, the PWN structure, such as the jets and tori, are not significantly perturbed. Once it leaves the SNR and enters the interstellar medium (ISM), the sound speed decreases (∼ a few to a few tens of km s −1 ) and the pulsar motion (typically a few hundreds of km s −1 ) becomes supersonic. The pulsar wind is then confined by the ISM ram pressure, resulting in a bow shock. A bow shock PWN has a compact head followed by a long tail that can extend up to a few parsecs (see e.g., Kargaltsev et al. 2017;Reynolds et al. 2017).
One special case is the Dragonfly (G75.2+0.1). This system is powered by PSR J2021+3651, an energetic (spin-down luminosityĖ = 3.4×10 36 erg s −1 ) and young (characteristic age τ c = P/2Ṗ = 17.2 kyr) pulsar with a rotation period of 103.7 ms. The pulsar emission has been detected in radio (Roberts et al. 2002), X-rays (Hessels et al. 2004), and γ-rays (Halpern et al. 2008). The source distance is controversial. While the pulsar dispersion measure implies over 10 kpc (Yao et al. 2017), X-ray spectral study suggests a closer value of 3-4 kpc (Van Etten et al. 2008). A later study obtained a even closer distance D = 1.8 +1.7 −1.4 kpc with 90% confidence interval, using both the X-ray spectrum and red-clump stars as standard candles and employing the extinction-distance relation along the pulsar line of sight (Kirichenko et al. 2015). We will adopt this value throughout our study. This partly overlaps with the estimate by Van Etten et al. (2008) and is also consistent with D ≈ 1 kpc based on the empirical γ-ray "pseudodistance" relation (Saz Parkinson et al. 2010).
On a large scale, the Dragonfly exhibits a cometary tail in X-rays and radio (see Aliu et al. 2014) that points to the general direction of the TeV source HAWC J2019+368 located ∼0.3 • away (Albert et al. 2021). The latter was suggested to be the pulsar birthsite (Mizuno et al. 2017) because TeV emission is originated from IC scattering of less energetic particles. It therefore traces the long-term history of a PWN system. By modeling the spectrum measured with HAWC and Suzaku, it was estimated that the true age of PSR J2021+3651 is ∼7 kyr (Albert et al. 2021).
Previous X-ray observations of the Dragonfly with the Chandra X-ray Observatory found that the PWN has axisymmetric main body of ∼ 20 long. This can be fit with a double tori model. There are also polar jets of ∼ 30 long extending from the pulsar (Hessels et al. 2004;Van Etten et al. 2008). Altogether the overall morphology resembles a dragonfly. The only other system having such a double tori morphology is the Vela PWN (Helfand et al. 2001). While the pulsar has a relatively young characteristic age (17.2 kyr), the host SNR has never been detected. X-ray observation shows a bow-shock like structure surrounding the PWN. This could indicate that the SNR has dissipated and the pulsar is traveling supersonically in the interstellar medium (ISM). This makes the Dragonfly a rare case that shows both bow-shock and torus-jet structure in X-rays (Hessels et al. 2004;Van Etten et al. 2008).
The typical features such as wisps and tori have been found both in radio and X-rays for some PWNe (e.g., the Crab Nebula; Dubner et al. 2017). Our study will supplement in understanding its origin by comparing these structures with simulated results (e.g., Olmi & Bucciantini 2019). The radio PWN can also trace the flow of the particles over longer time scale because of the longer cooling time. In order to look for the radio counterparts of the X-ray structure, we analyzed new VLA 4-8 GHz (C-band) data with arcsecond resolution comparable with Chandra images along with archived 1-2 GHz (L-band) in this study. We will compare the radio bow-shock structure with that observed in X-rays in more details than the previous study in Aliu et al. (2014). The magnetic field is a critical factor in the particle acceleration process. The new data with full polarization measurements allow us to map the magnetic field structure of the nebula in details for the first time. With this new observation in radio, we will constrain the injected particle spectrum to test the different acceleration theoretical models. By combining radio and X-ray spectra we can estimate the B-field strength or age. The observations and data analysis are reported in Section 2 and results are presented in Section 3. Based on these results, we discuss the morphology, injected particle spectrum, polarization properties of the system in Section 4.

C band
We carried out new radio observations of the Dragonfly on 2017 November 2 and December 1 using the Karl G. Jansky Very Large Array (VLA) in the B array configuration. The observations were done in the C-band (4-8 GHz) with full polarization for a total on-source time of 3.5 hr. Data were taken in the continuum mode comprising 32 spectral windows with a total bandwidth of 4 GHz. The shortest array spacing is 243 m, corresponding to a largest angular scale of ∼ 40 that our radio maps are sensitive to. 3C147 was chosen for calibrating the flux density scale as well as the bandpass by the observations. To determine the antenna gains, a secondary calibrator J2015+3710 was observed every 20 minutes (Perley & Butler 2017).
The VLA data reduction was performed using the Common Astronomy Software Applications (CASA; McMullin et al. 2007). We first applied the VLA CASA Calibration Pipeline (Version 5.6.1-8) to flag the shadowed or not-on-source antennas. The edge channels of each spectral window and of each baseband are flagged. The pipeline then removed radio frequency interference (RFI). We solved for the multi-band cross-hand delays of the two basebands (∼2-4 GHz and ∼4-6 GHz) separately. Additional RFI flagging was performed after the pipeline processing.
After flagging and calibration, we employed the task tclean to generate images at 6 GHz in CASA. We adopted a weighting parameter of robust = 1.5 to optimize between the sidelobe suppression and reducing noise levels (Briggs 1995). The threshold for deconvolution is 3 times of the rms flux density. After the first round of cleaning, we applied a clean mask on the image interactively to restrict the region to be cleaned. The final continuum 6 GHz total intensity image has a circular restoring beam of FWHM 3 and rms noise of ∼4 µJy beam −1 , which is similar to the theoretical value.
To generate polarization maps, we first applied similar flagging and calibration procedures as above to the polarization calibrators, then determined the instrumental delay between the two polarization outputs. After that, we solved for the instrumental polarization (the frequency-dependent leakage terms, so-called "Dterms") as well as the polarization position angle using 3C138 and 3C147 as calibrators respectively (Perley & Butler 2013). A natural weighting (robust = 2.0) was used to minimize the noise level. We deconvolved the images in Stokes Q and U simultaneously and restored them with a Gaussian beam of FWHM=4 to boost the signal-to-noise ratio (S/N). The rms noise of the images is compatible with the theoretical levels of ∼3 µJy beam −1 . The polarized intensity and polarization angle maps were created using the Stokes Q and U images with the task immath. We employed the task rmfit to calculate the rotation measure (RM) map, the intrinsic polarization angle (PA0) map after correction for Faraday effect, and the corresponding error maps.
Finally, we divided the whole 4 GHz bandwidth into two equal subbands of 2 GHz each, centering at 5 and 7 GHz, to measure the spectral index. We applied the same imaging and deconvolution procedure as above and restore with a slightly larger circular beam of FWHM 4 to boost the S/N.

L band
We processed an archival L-band observation of the Dragonfly taken with the VLA. The data were obtained as part of a small survey of Galactic EGRET sources in the field around PSR J2021+3651 in the D array on 2000 July 20 and in the C array on 2002 November 2 at a central frequency of 1.47 GHz with a bandwidth of 50 MHz. The flux (primary) calibrator used was J2015+3710 and 3C286 (=J1331+3030) were chosen for phase and bandpass calibration. A three-point mosaic centered on the pulsar was used, with each pointing in the mosaic getting roughly 85 minutes of exposure in C array and 22 minutes in D array. Therefore the shortest spacing of the baseline was 35 m, which correspond to a angular resolution of ∼ 16 in the L-band. The ATNF MIRIAD package (Sault et al. 1995) was used to reduce the data and generate the image. The mossdi clean algorithm was used for deconvolution. Two bright sources near the edge of the field of view were subtracted from the u-v visibilities and multiple rounds of self-calibration were done. Mainly phase self-calibration was used, although a few iterations of amplitude calibration were also done. The final gain restored and gain corrected image has a 12.4 × 11.3 beam with a noise level of roughly 0.13 mJy beam −1 rms near the center of the image.

X-ray
Previous X-ray studies of the Dragonfly show that its spectrum is sensitive to the choice of region (Van Etten et al. 2008;Mizuno et al. 2017). In order to perform multi-wavelength comparison at different regions, we reanalyzed archival Chandra X-ray observations (ObsIDs 7603 and 8502) taken in 2006. The data were reprocessed using the script chandra repro script in Chandra Interactive Analysis of Observations (CIAO) (Fruscione et al. 2006) version 4.13. Both datasets were merged by applying the script merge obs in CIAO after astrometric check. We extracted the source spectra from regions defined by the radio morphology and background spectrum from nearby source free regions. The X-ray spectra were created using the specextract script of CIAO. We grouped the energy bins with at least 20 counts each and fit the spectra from both datasets simultaneously with the package Sherpa (Freeman et al. 2001). We employed an absorbed power-law model with the column density fixed at 6.7×10 21 cm −2 following previous studies. All the X-ray spectral analyses were made in the energy range of 0.5-8 keV.

Radio intensity map
The radio intensity map of the Dragonfly at 6 GHz is shown in Figure 1. The PWN is clearly detected. It has a size of 1. 8×0. 6 and is elongated along a position angle (PA) of ∼ 141 • (north to east). The overall structure consists of four parts: central resolved pulsar region, a western "main-body" region, a fainter eastern "bridge", and a dark gap lying between them. The surface brightness at the pulsar location is ∼ 26 µJy beam −1 . Since the FWHM of the central peak at the position of the pulsar is ∼ 10 , which is larger than its beam size resolution (FWHM=3 ), it should be a compact nebula surrounding the pulsar. The central pulsar region is surrounded by a faint ring-like region with ∼5 thickness. The latter has only half of the average surface brightness of the PWN. There are two brighter lobes southwest of the main PWN and another fainter one further south, which are marked by the dashed contours in Figure 2a.
The western main-body region dominates most of the radio flux density of PWN at 6 GHz. It is ∼100 long and ∼20 wide and has a mean surface brightness of ∼20 µJy beam −1 . It contains two wings that are not perfectly symmetric about the central pulsar. The northwestern (NW) wing has length ∼60 , longer than the southeastern (SE) wing of ∼40 . The radio surface brightness is not uniform in the main-body region. The NW wing has higher peak brightness (∼36 µJy beam −1 ) than the SE one (∼28 µJy beam −1 ). Besides, ∼14% of the area in the NW half has surface brightness ≥ 25 µJy beam −1 , but only < 5% in the SE wing, despite the mean C-band surface brightness in both wings are about the same (20 µJy beam −1 ).
On the eastern side of the PWN there is a relatively faint, long, and narrow bridge region with size ∼50 ×6 . The average surface brightness is ∼16 µJy beam −1 . In between the main-body and the bridge regions, we see a clear 5 wide gap. Its minimum intensity drops to ∼5 µJy beam −1 , which is only one-fifth of the peak emission at the pulsar position and is very close to the rms noise of 4 µJy beam −1 .
The VLA L-band intensity map is shown in Figure 2a with the C-band regions labelled. Overall, the coneshape L-band emission extends from east to west over 10 . It shows two peaks, one at the center of the mainbody and one at 6 west of the pulsar. The brightness is not uniformly distributed but rather patchy. It is unlikely due to the lack of short spacing because the largest angular size that the observation is sensitive to is around 16 . The elliptical C-band emission is located at the head of the cone with ∼50 offset from the sharp edge in the east. The two brighter lobes southwest of the main PWN in the C-band image appear to correlate with the diffused L-band emission and another fainter patch is located near the edge of the L-band emission.

Multi-wavelength comparison
The cone-shape large scale radio structure shown in the L-band is also found in the X-rays (see Figure 2b). The X-ray emission peaks at the pulsar while the peak in the L-band is offset to the south-west for ∼20 . Moving to a smaller scale, the elliptical radio C-band emission is about two times as large compared with the X-ray PWN as shown in Figure 3. The radio PWN's elongation aligns with the X-ray torus and is perpendicular to jets. The wings of the X-ray PWN lie within the radio C-band main-body while the X-ray jet partially overlaps with the radio bridge region. In X-rays, the NW wing, like its radio counterpart, is also slightly larger than its opposite half. We also found that the faint jet seen ahead of the bow shock in the Chandra image (Van Etten et al. 2008) is pretty clear in the L-band (see Figure 2). We plot the cross-section profiles of the radio and Xray intensity maps along the minor and major axes of the PWN in Figure 4. It is clear that the central peak at the pulsar position is wider than the radio beam and the X-ray PSF, implying that there is diffused emission surrounding the pulsar in both radio and X-rays. This can be supported by the fact that the measured radio flux density at the pulsar location is ∼1 mJy at the Lband, much larger than the reported pulsed emission of ∼0.1 mJy (Roberts et al. 2002). The brightness profiles in all three bands are roughly asymmetric about the pulsar position and the intensity falls off more rapidly on the northeast (NE) side along the minor axis.
Along the minor-axis cross-section, the C-band intensity drops from the central peak of ∼ 28 mJy to ∼5%-15% mJy on either side of the pulsar in the ∼10 -wide  . Radio L-band (dotted), radio C-band (solid), and X-ray (dashed) brightness profiles along the (a) minoraxis and (b) the major axis cross-sections. These are measured from a 12 wide region (the large rectangle in the inset plots) for the VLA L-band (1.5 GHz) and a 3 wide region (the small rectangle in the inset plots) for both the VLA Cband (6 GHz) and the Chandra X-ray (0.3-7 keV exposurecorrected and binned to 0.984 ) intensity maps along two axes.
"valleys" and rises to a second peak (see Figure 4a). The two second peaks correspond to bridge and the main-body region as seen in the radio C-band intensity map (see Figure 1). The X-ray profile along the minor axis shows a flat bump (within ∼ 5 west of the pulsar) on the south-west (SW) side of the X-ray profile (Figure 4a). It was fit with a double-torus model in previous X-ray studies (Hessels et al. 2004;Van Etten et al. 2008). In the L-band, the valleys and the second peaks are difficult to identify along the minor axis and are barely seen along the major axis probably due to resolution limit (FWHM ∼12 ). The plateau emission profiles on the western side in both radio bands indicate that the nebula extends to the west, which agrees with the large scale X-ray emission detected with Chan-dra and XMM-Newton (Van Etten et al. 2008;Mizuno et al. 2017).

Radio spectrum
We measure the C-band flux densities of the substructures as shown in Figure 1 at two sub-bands centered at 5 and 7 GHz. For the L-band (1.47 GHz), we only measure the flux densities of the larger regions due to limited resolution. The errors are estimated using the uncertainty of the calibration calibrator and the rms of the emission free region. The spectral results are plotted in Figure 5. The spectral index (α) of different regions are obtained by fitting a power-law to the flux density (S ν ) at different frequencies: S ν ∝ ν α . The results are listed in Table 1. The best-fit spectral index between 5 and 7 GHz α 5−7 GHz is −0.76 ± 0.03 for the overall radio PWN. The main-body region is much flatter (i.e., harder) (α 5−7 GHz = −0.41±0.04) than the bridge region (α 5−7 GHz = −1.33 ± 0.27). If we take the 1.47 GHz flux densities into account, the spectral index of the whole PWN and the main-body region between 1.47 and 7 GHz are α 1.5−7 GHz = −0.74 ± 0.02 and −0.44 ± 0.04, respectively, which are consistent with the result from only the C-band data. Since the angular size of the radio PWN is larger than the largest sensitive scale (∼ 35 at 5.0 GHz and ∼ 25 at 7.0 GHz) of our observation, our data may not be able to recover all the flux densities at large scales.
or 5−7 GHz. Therefore we argue that this is not an issue for the main-body, since the two spectral indexes α 1.5−7 GHz and α 5−7 GHz are consistent. Finally, we calculate the spectral index map between the 5 GHz and 7 GHz images and the result is shown in Figure 6. The spectral index increases (flattens) gradually from the bridge with the α −1 to the main-body region with the mean and median spectral indexes of −0.1 and 0.1 respectively. Figure 6. Spectral index map of the Dragonfly PWN between 5 and 7 GHz. The map is clipped if the total intensity in either map is below 1-σ detection. The circle located at the lower left corner shows the restored beam size and the location of PSR J2021+3651 is marked with an "x". The PWN region is the same as that in Figure 1.

Multi-wavelength spectrum
Our X-ray spectral analysis found a photon index of Γ = 1.55 ± 0.05 for the main body, slightly softer than The solid line indicates the main-body spectrum. For comparison, spectra from the cone-shape nebula region (excluding the X-ray PWN region), which similar to the "outer nebula" region in Van Etten et al. (2008), and the Suzaku observation region from Mizuno et al. (2017) are shown by the dashed and dotted lines, respectively. We extrapolated the lines to the radio band by assuming a spectral break with ∆α = 0.5 due to synchrotron cooling. Γ = 1.47 ± 0.04 for the central X-ray PWN (excluding the inner radius of 2 pulsar region). Moving further out, we found that the X-ray spectrum softens to Γ = 1.77 ± 0.06 in the cone-shape nebula region (excluding the main-body region) and Γ = 2.05 ± 0.12 in an even larger region (15 × 10 ) based on Suzaku observation (Mizuno et al. 2017). The resulting broadband spectral energy distributions (SEDs) are plotted in Figure 7. We find that the X-ray and radio spectra of the main-body region can be connected with a single, unbroken power law. For the cone-shape nebula and the Suzaku observation regions, the C-band data do not have enough sensitivity to provide flux density measurements. We therefore show in the plot only the L-band data and assumed a spectral break with ∆α = 0.5 between radio and X-rays, as expected from synchrotron cooling (see Pacholczyk 1970, for details). Figure 8a shows the fractional polarization map. The entire PWN is highly linearly polarized with polarization fraction (PF) 30% and even > 40% in the two wings of the main-body and the bridge regions. The PF distribution is non-uniform and appears to trace the X-ray PWN morphology, so it is unlikely due to noise. The PF is higher on the X-ray PWN wings but much lower (≤ 30%) in the central part.

Polarization
We attempted to derive a rotation measure (RM) map using the polarization angles in the two subbands at 5 GHz and 7 GHz. However, the resulting RM map has too large uncertainty to be useful. We therefore adopt a constant RM value of −524 rad m −2 of the pulsar (Abdo et al. 2009) to derotate the polarization vectors. The intrinsic magnetic field orientation after Faraday rotation correction is shown in Figure 8b. The magnetic field in the bridge region is predominately along the north-south direction. It then switches over the gap region to mostly east-west in the main-body. The change in direction is unlikely caused by fluctuation of the foreground RM as it shows a good correlation with the overall PWN structure. Finally, we note that there are polarization signals detected in the outer lobes south-west of the PWN.

Pulsar motion
The radio band emission of the Dragonfly in the Lband shows an overall bow-shock morphology, suggesting supersonic motion of the pulsar in the ISM. This is supported by the asymmetry between the NW and SE wings seen in the radio C-band and X-rays, and the distortion of the X-ray jet and the B-field vectors. Nonetheless, we note that the torus-jet structure and the toroidal B-field geometry are not totally disrupted. These imply that the pulsar motion is at most mildly supersonically. For a typical supersonically moving pulsar, its ram pressure is so large that its microstructures would be destroyed (see examples in Kargaltsev et al. 2017).
The apex of the bow shock structure in both radio Lband and X-rays is ∼50 from the pulsar, corresponding to a physical scale r ≈ 0.4 pc at 1.8 kpc. If we take this as the bow shock stand-off distance, the pulsar velocity (v PSR ) can be estimated from the balance between the PWN internal pressure and the ram pressure, such that r = Ė /4πcρv 2 PSR , where ρ is the ambient density. The latter is difficult to determine and the typical values in the cold and hot phases of the ISM give a wide range of v PSR estimates from a few to hundreds of km per second Van Etten et al. 2012).
If the pulsar was born at the centroid of the TeV emission, the ∼0.3 • offset suggests a high transverse velocity of ∼300-2500 km s −1 at a distance of 1.8 +1.7 −1.4 kpc, where the large range of the velocity estimate mostly come from the distance uncertainty. This is possible but the upper bound lies in the high end of the pulsar velocity distribution . In any case, the velocity is clearly supersonic given that the typical ISM sound speed is ∼10-100 km s −1 . We note that the proper motion inferred by the TeV position misaligns with the pulsar spin axis. It would affect the B-field orientation and its interaction with ambient ISM according to the 3-D simulations ) that will be discussed in the later subscetions.

Radio PWN Structure
At a distance of 1.8 kpc, the elliptical morphology of the Dragonfly PWN in radio corresponds to a physical size of 0.9 pc×0.3 pc, comparable to the Vela PWN in radio (Dodson et al. 2003). The radio emission of the Dragonfly has the same overall orientation as the X-ray counterpart, but about two times as large in size. This is a common feature among PWNe due to the longer synchrotron cooling time of the radio-emitting particles.
We can compare the Dragonfly with Vela and B1706 as their central pulsars share very similar properties, including characteristic age, spin-down power, and γ-ray pulsations. On the other hand, their X-ray PWN structures such as torus and jet show different morphology (Dodson et al. 2003;Kargaltsev & Pavlov 2004;Romani et al. 2005;Liu et al. 2022). The radio emission of Vela and B1706 PWNe are found outside the X-ray regions and it is faint at the inner PWN. However, the two wings for the Dragonfly in radio are connected and overlap with the X-ray emission. It is unclear if the difference is due to different viewing geometry or other physical mechanisms.
In the L-band image, a protruding structure ahead of bow shock (see Figure 2) was also found in the Chandra image (Van Etten et al. 2008). The 3-D simulation suggested these kinetic jets are formed due to the leaking PWN particles into the ISM ). These particles escape along the reconnected magnetic field between the PWN and the ISM (Bandiera 2008).
Moving to small scale, our high resolution radio map reveals, for the first time, peculiar features in the Dragonfly. The most intriguing one is the narrow bridge region that has different properties than the main body, including steeper radio spectrum, different B-field orientation, and no X-ray counterpart. It was suggested that the pulsar was born near the TeV source and moves towards east. In this picture, the bridge is then the leading edge of the PWN compressed by the ISM ram pressure. The compression could enhance the B-field and hence the synchrotron emissivity, and perhaps increases the polarization fraction. This is supported by the sharp edge of the bridge in NE in both the C-and L-band images in Figure 4a and also slight deformation of the X-ray PWN in NE as indicated by the white X-ray contours in Figure 3. The pulsar motion could also make the NW wing, which trails the pulsar motion, longer than the SE wing of the main body, in both radio and X-ray bands.
Alternatively, the observed bridge could correspond to the far side of the radio PWN if it has a ring-like structure in 3-D and is viewed near edge-on. Such a ring (with a diameter of ∼0.6 pc in the Dragonfly) could resemble that formed at the termination shock detected in the Crab Nebula in X-rays and radio (e.g., Krassilchtchikov et al. 2017;Dubner et al. 2017) and the putative ring seen in the Vela PWN (Ponomaryov et al. 2018). A ring-like structure is observed in X-rays (called 'double tori' in Van Etten et al. 2008) but is smaller than the radio ring. Perhaps they are formed by different populations of particles. Another possibility is that the X-ray ring is newly produced, just like the X-ray ring of the young Crab Nebula, while the radio ring is the remaining older ring with its X-ray component having cooled down. In this case, the bridge would appear as fainter and narrower than the main-body, which corresponds to the near side, due to Doppler boosting (see Ng & Romani 2004, 2008. The brightness ratio between the two regions is R ≈ 1.15. The inclination angle ζ between the ring axis and the line of sight can be deduced by the relation R = [(1 + β cos ζ)/(1 − β cos ζ)] 2−α . We adopt the best-fit post-shock velocity β ∼ 0.8 inferred by the X-ray torus fitting (Van Etten et al. 2008) and the Cband spectral index α = −0.4 of the main-body. These give ζ ≈ 88 • , which is very close to ∼ 85 • from fitting the X-ray torus Van Etten et al. (2008). This can also explain the rapid fall off of the L-band and C-band intensity profiles seen on the eastern side of the minor-axis cross section in Figure 4a. This model, however, is less preferred, as the ring size is much larger than the X-ray tori. The latter is believed to correspond to the wind termination shock. The formation mechanism of such a radio ring is unclear. This picture also has difficulty explaining the different spectral index in the bridge and the main-body.

Multi-wavelength Spectral Analysis
The radio spectral index of the overall radio PWN is ∼ −0.7 to −0.8, which is much steeper, i.e. softer, than the typical radio PWN spectrum (−0.3 α 0), such as Crab (α = −0. 27 Strom & Greidanus 1992) and G283.2-0.59 (α = −0.13 Ng et al. 2017). As we mentioned, this soft emission is mainly contributed by the bridge region, probably due to compression as the pulsar moves in the ISM. If we consider only the main-body region, which has larger S/N and contributes ∼70% of the PWN flux density at 6 GHz (see Figure 5), its spectral index is ∼ −0.4, more inline with the typical value.
The SED of the main-body region from radio to X-rays can be fit with a single, unbroken powerlaw as shown in Figure 7. This suggests no significant synchrotron burn-off, implying a young age for the particles. For the cone-shape nebula and the Suzaku region, if we assume a break with ∆α = 0.5 between the X-ray and radio spectra, we can infer spectral turn-over (ν b ) at log[ν b (Hz)] = 11.2 +0.7 −1.0 and 14.9 +0.5 −1.0 , respectively. Future direct measurement of their radio spectra can confirm this.

Nebular Magnetic Field Strength
We estimate the so-called equipartition magnetic field strength (B eq ) of the PWN by minimizing the total energy of the particles and the magnetic field. Using standard synchrotrn theory, it can be shown that B eq is related to the synchrotron luminosity L and the emission volume V by where c 12 is a constant weakly dependent on the spectral index and the frequency range (Pacholczyk 1970).
Assuming that the overall PWN has a cone structure in 3D as suggested by the L-band image and taking the main body and the bridge regions as slices of the cone, we estimate their volumes of ∼ 4.1 × 10 54 d 3 1.8 cm 3 and ∼ 3.9 × 10 53 d 3 1.8 cm 3 , respectively, where d 1.8 is the source distance in units of 1.8 kpc. The calculated luminosities (L) are 6.9×10 29 d 2 1.8 erg s −1 for the main-body and 2.3×10 29 d 2 1.8 erg s −1 for the bridge region, based on the integrated flux in the frequency range of 10 7 −10 11 Hz. These give B eq of ∼ 22 µG for the mainbody and a higher value of ∼ 59 µG for the bridge. The result could indicate magnetic field enhancement due to compression by the ISM ram pressure.
For the cone-shape nebula region excluding the inner PWN, similar estimate gives B eq ≈ 21 µG assuming the spectrum shown in Figure 7. This is very close to the estimate of ∼ 20 µG by extrapolating the B-field from the pulsar surface (Van Etten et al. 2008). Similarly, by taking the Suzaku observation region as a cylinder in 3D, we derived B eq ≈ 19 µG, inline with the value for the cone-shape nebula above.
Synchrotron cooling can give rise to a spectral break at frequency where B is the magnetic field strength and t is the age of the system. Using B ∼ 21 µG estimated above and the pulsar's age of 7 kyr, we expect ν b ∼ 10 15 Hz in the cone-shape nebula. This is much higher than the spectral turn-over at ν b ∼ 10 11 Hz inferred from a break of ∆α = 0.5 between the X-ray and radio spectra (see Figure 7). This can be reconciled if the actual B-field is higher than the equipartition estimate or the true age is much larger than 7 kyr. The value ν b ∼ 10 15 Hz obtained from the equation above, on the other hand, is consistent with the inferred cooling break for the Suzaku observation region. X-ray spectral analysis found no obvious photon index variation in the Suzaku region (Mizuno et al. 2017), suggesting that the particles probably have already cooled down and hence the inferred spectral break in Figure 7 could be real.

Magnetic Field Geometry
Our high resolution polarization maps reveal different B-field configurations in different part of the Dragonfly: it is oriented along the north-south direction in the bridge and east-west in the main body. There is no obvious connection with the X-ray torus structure, contrary to a toroidal B-field expected in young PWNe (e.g., Porth et al. 2017). On the other hand, the field geometry seems better correlate with the inferred pulsar motion direction, similar to some bow-shock PWNe, e.g., the Mouse (G359.23−0.82) and the Frying Pan (G315.78−0.23) (Yusef-Zadeh & Gaensler 2005;Ng et al. 2012). In particular, the Mouse has a B-field perpendicular to the pulsar motion ahead of the bow-shock apex. It then switches to parallel after a depolarized region at the shock (Yusef-Zadeh & Gaensler 2005). For the Dragonfly, the B-field also changes direction across the faint gap between the bridge and the main body. Note that the Mouse is believed to be moving highly supersonically (see Klingler et al. 2018) but the Dragonfly is at most mildly supersonically as we argued. It is intriguing that they share some similarities in the B-field geometry. This suggests that the Mach number of the pulsar motion is not the only factor determining the field configuration.
In the main body, the magnetic field orientation can be described by the sum of two components: one follows the torus direction and one follows the pulsar moving direction. This is different than the Mouse, which has Bfield parallel to the pulsar tail. It was suggested that the central pulsar of the Mouse has spin axis perpendicular to its motion (Klingler et al. 2018), such that the two B-field components align.
On the other hand, the magnetic field configuration in the bridge region is different from that of the main body. It fits better with the picture of B-field reconnection between the pulsar wind and the ISM in the head of the PWN, as supported by the 3-D simulation ). This B-field reconnection allow relativistic particles escaping from the pulsar wind and form the kinetic jet as discussed in Section 4.2.
Finally, we note that our polarization result is inconsistent with the idea that the bridge is the far side of a ring-like structure in 3D, since that should give a similar B-field geometry in both the bridge and the main body.

Diffuse Emission Lobes
There are a few diffuse emission lobes located southwest of the pulsar (see Figure 2a). These could be larger scale emission resolved out due to lack of short spacing in our observation. We argue that they are likely real emission rather than sidelobes, as they are located in the pulsar's path and there is also patchy radio emission detected in the L-band. They could be relic pulsar wind particles interacting with inhomogeneous ISM (Toropina et al. 2019) or leftover after the pulsar passed by.
For the southernmost lobe shown in the C-band image in Figure 2a, its emission in the L-band is rather faint. If this extended emission is real, it could indicate an inverted spectrum, different than that of the inner PWN. This requires some re-acceleration mechanism, or it could be caused by instability at the bow-shock interface.

CONCLUSION
In this paper, we reported on a detailed radio study of the Dragonfly PWN using new VLA observations at the C-band (4-8 GHz) together with archival radio Lband VLA data and Chandra X-ray data. The radio emission in the C-band shows elongated structure, with size twice as large as the X-ray counterpart. It consists of a bright, axisymmetric main-body in the west and a narrow bridge in the east. We suggest that the latter is due to compression by the ram pressure when the pulsar travels supersonically in the ISM. This is supported by the high polarization fraction, the soft spectrum (α 1.5−7 GHz = −1.3) and the B-field configuration of the bridge. For the main-body region, the radio spectrum is harder (α 1.5−7 GHz = −0.4), closer to that of a typical PWN, and it can be connected to the X-ray spectrum by a power law without a break. This suggests that the X-ray emitting particles are probably not yet cooled down.
The radio emission of the Dragonfly at 6 GHz has a high linear polarization fraction of over 30% in most regions and even higher in the two wings of the Xray PWN and the bridge. This indicates a highly ordered magnetic field structure. The intrinsic B-field in the bridge predominately are perpendicular to the pulsar motion. It then switches direction across the faint gap and the field configuration in the main body can be described by a toroidal pattern slightly distorted by the pulsar motion. Similar behavior was observed in another PWN that travels highly supersonically in the ISM. However, for the case of the Dragonfly, we argue that the pulsar moves transonically or mildly supersonically since the stand-off distance is large, its jet and torus structures and the B-field are not disrupted.  (Sault et al. 1995) , CIAO (Fruscione et al. 2006), Shepra (Freeman et al. 2001)