A multi-wavelength investigation of PSR J2229+6114 and its pulsar wind nebula in the radio, X-ray, and gamma-ray bands

G106.3$+$2.7, commonly considered a composite supernova remnant (SNR), is characterized by a boomerang-shaped pulsar wind nebula (PWN) and two distinct ("head"&"tail") regions in the radio band. A discovery of very-high-energy (VHE) gamma-ray emission ($E_\gamma>100$ GeV) followed by the recent detection of ultra-high-energy (UHE) gamma-ray emission ($E_\gamma>100$ TeV) from the tail region suggests that G106.3$+$2.7 is a PeVatron candidate. We present a comprehensive multi-wavelength study of the Boomerang PWN (100"around PSR J2229+6114) using archival radio and Chandra data obtained from two decades ago, a new NuSTAR X-ray observation from 2020, and upper limits on gamma-ray fluxes obtained by Fermi and VERITAS observatories. The NuSTAR observation allowed us to detect a 51.67 ms spin period from the pulsar PSR J2229+6114 and the PWN emission characterized by a power-law model with $\Gamma = 1.52\pm0.06$ up to 20 keV. Contrary to the previous radio study by Kothes et al. 2006, we prefer a much lower PWN B-field ($B\sim3$ $\mu$G) and larger distance ($d \sim 8$ kpc) based on (1) the non-varying X-ray flux over the last two decades, (2) the energy-dependent X-ray PWN size resulting from synchrotron burn-off and (3) the multi-wavelength spectral energy distribution (SED) data. Our SED model suggests that the PWN is currently re-expanding after being compressed by the SNR reverse shock $\sim 1000$ years ago. In this case, the head region should be formed by GeV--TeV electrons injected earlier by the pulsar propagating into the low density environment.


INTRODUCTION
Pulsar wind nebulae (PWNe; Cholis et al. 2018) are believed to generate a majority of the energetic leptons in our galaxy.The pulsar's rotating magnetic fields produce a wind of highly relativistic particles that expand out into the shell of the supernova remnant (SNR).High-energy observations of dozens of PWNe detected synchrotron and inverse-Compton upscattering (ICS) of cosmic microwave background photons, ambient infrared (IR) or optical stellar radiation in the X-ray and TeV bands, respectively, suggesting that non-thermal particles are accelerated to TeV or even PeV energies within the PWNe (Arons 2012).The PWN evolution is characterized by three stages: (1) young, termination-shock driven wind nebulae, (2) middle-aged PWNe interacting with their host SNRs, and (3) relic PWNe (Gaensler & Slane 2006;Giacinti et al. 2020).Relativistic winds from the pulsar injected into the SNR center are abruptly decelerated in an inward-facing termination shock, at which particles are accelerated to TeV energies; the post-shock flow further decelerates until it reaches pressure equilibrium with the SNR interior.The bubble of shocked pulsar wind is the observed PWN, which continues to expand until the deceleration of the outer SNR blast wave sends a reverse shock back toward the center, compressing and re-brightening the PWN, at an age of order 1-10 kyr.The PWN continues to interact with the SNR interior until either the SNR dissipates or the pulsar, if born with a substantial kick velocity, escapes the SNR shell, continuing to inflate a PWN.In either case the PWN now interacts directly with the interstellar medium (ISM) (a "relic PWN"; Cholis et al. 2018), often in the shape of a bow-shock nebula.These middle-aged PWNe manifest a vast diversity of highly anisotropic non-thermal emission in multi-wavelength bands.The composite system is formed by its relic PWN interacting with the ambient medium and SNR reverse shock, exhibiting peculiar radio and X-ray morphology (often with nicknames such as Rabbit and Snail).Composite SNRs are of particular interest because they manifest on-going PWN-SNR interaction sites and possibly accelerate particles to TeV-PeV energies (Ohira et al. 2018).Some of the middle-aged PWNe are associated with the PeVatron candidates detected by HAWC and LHAASO above E γ ∼ 100 TeV (Abeysekara et al. 2020;Cao et al. 2021).Eventually, after τ ∼ 100 kyr, electrons and positrons escape from relic PWNe and form extended TeV halos, as revealed around the Geminga and Monogem pulsars (Abeysekara et al. 2017).TeV halos are a new class of gamma-ray sources which are suggested to be the primary source of the positron excess observed at Earth (López-Coto et al. 2022).How and when a PWN evolves through these transitions depends on the progenitor star's characteristic properties and environment within the ISM.Hence, multi-wavelength observations of PWNe in different evolution stages and environments are essential for understanding how particles are injected from the pulsar, diffuse out while cooling, and interact with the ambient gas and their host SNRs.The Boomerang region is one of the most remarkable composite SNRs for its complex multi-wavelength morphology and the recent detection of gamma rays above 100 TeV indicating it to be a PeVatron candidate.Its large-scale radio emission (G106.3+2.7)consists of a compact boomerang-shaped nebula around the radio pulsar PSR J2229+6114 and cometary structure extending toward the southwest.The radio source G106.3+2.7 was first identified as a SNR by Joncas & Higgs (1990) following the Dominion Radio Astrophysical Observatory (DRAO) survey of the northern Galactic plane.Using further DRAO observations in the 408 MHz and 1420 MHz continuum bands, Pineault & Joncas (2000) discerned two distinct regions of SNR G106.3+2.7,labeled the head and the tail (See Figure 1).The head region is characterized by its higher surface brightness and flatter spectral index in comparison to the elongated tail region.Using VLA observations at 20 and 6 cm, as well as ROSAT and ASCA observations, Halpern et al. (2001b) identified a compact radio and X-ray source in the northeast area of the SNR G106.3+2.7 head region and suspected it to be a pulsar with a corresponding PWN.The radio and X-ray detections of a 51.6 ms pulsation from the pulsar, now known as PSR J2229+6114, confirmed this hypothesis (Halpern et al. 2001a).Further radio and X-ray timing studies of the pulsar led to determining a spin-down power of 2.2 × 10 37 erg s −1 and a characteristic age of ∼10 kyr (Halpern et al. 2001a).A compact PWN with a r ∼ 100 ′′ extension was detected in the radio band and was suggested to be associated with SNR G106.3+2.7 based on the subsequent measurement of the same peak H I velocity from the compact Boomerang nebula and the head region (Kothes et al. 2001).While SNR G106.3+2.7 has been labeled as an SNR, no thermal X-ray emission is reported anywhere in the Boomerang complex, and no large-scale radio morphological features are evident that might suggest the supernova blast wave.The larger-scale integrated radio spectral index is −0.61 (Kothes et al. 2006), while that of the PWN alone is ∼ 0 (Halpern et al. 2001a), suggesting a shock acceleration source for the larger scale electrons, but there is no edge brightening apparent in any location.
It has been hypothesized that the Boomerang's shape could be caused by a bow-shock between PSR J2229+6114 and its surrounding medium.However, this was deemed unlikely, as simple modelling of the system under this assumption resulted in a supernova explosion energy far below anything ever recorded; the pulsar also does not lie at the apex of Pope et al.Yao et al. (2017); Abdo et al. (2009b) 12 kpc H I radial velocities Pineault & Joncas (2000) the bow structure (Kothes et al. 2001(Kothes et al. , 2006)).In contrast, based on the boomerang-like radio morphology as well as its proximity to the northeast boundary of SNR G106.3+2.7, it was postulated that the PWN had been crushed by a SNR reverse shock (Kothes et al. 2001(Kothes et al. , 2006)).Further observations of the Boomerang region with the Effelsberg 100-m telescope led to a hypothesis that an interaction with the SNR reverse shock could also account for the low radio luminosity of the PWN with respect to the spin-down power (Kothes et al. 2006).Furthermore, a radio spectral break observed between 4 and 5 GHz was attributed to synchrotron cooling.Under this assumption, Kothes et al. (2006) suggested that the PWN B-field is 2.6 mG and that the PWN was crushed by the SNR reverse shock 3900 years ago.A more recent study based on a model of the diffusion of the relativistic electrons injected into the PWN and X-ray radial profile suggested that the PWN B-field is 140 µG (Liang et al. 2022).
In the X-ray band, the pulsar, its compact nebula (r ∼ 30 ′′ ) and diffuse emission over r < ∼ 100 ′′ were detected by two Chandra observations (17 and 94 ksec) in 2001-2002(Halpern et al. 2001a)).XMM-Newton and Suzaku observations revealed more extended diffuse X-ray emission from the head and tail region (Ge et al. 2021;Fujita et al. 2021).The long Chandra observation in 2002 unveiled a point source at the pulsar position, an incomplete torus of r ∼ 10 ′′ and a jet-like feature.These X-ray features resemble those of the Vela PWN whose pulsar's motion is aligned with its X-ray jet (Halpern et al. 2002).The Chandra ACIS image of the Boomerang PWN was fit by a 3D torus model (Ng & Romani 2004).The brighter side of the torus (west; see Figure 3) is due to Doppler boosting of mildly relativistic magnetic hydrodynamic (MHD) outflow from the termination shock.The best-fit torus model predicts that the pulsar should be moving along the spin-axis (i.e. the jet direction toward the northwest).The prediction that the pulsar is moving toward the northwest does not seemingly agree with the tail morphology of SNR G106.3+2.7.However, the recent numerical studies studying the evolution of PWN while interacting with a host SNR show that the PWN's morphology depends on both the pulsar's proper motion and the region's density gradient (Kolb et al. 2017).
Gamma-ray emission has been observed in the SNR G106.3+2.7 region from GeV up to few hundreds of TeV energy range.The Large Area Telescope on board of Fermi gamma-ray space telescope (Fermi-LAT) detected GeV emission coincident with PSR J2229+6114, which was also associated with EGRET source 3EG J2227+6122 (Abdo et al. 2009a;Hartman et al. 1999).Gamma-ray pulsations were observed above 0.1 GeV, confirming GeV emission originates from PSR J2229+6114 (Abdo et al. 2009b).Using a collection of Fermi-LAT data, Xin et al. (2019) identified emission between 3-500 GeV coincident with the tail region with a source radius 0.25 • .The Very Energetic Radiation Imaging Telescope Array System (VERITAS) detected TeV emission from the tail region and found that the centroid of the TeV source overlaps with 12 CO cloud J = 1−0 emission (Acciari et al. 2009).More recently, the MAGIC collaboration reported TeV detection of the head region as well (Oka 2021).Gamma-ray emission with energies higher than 100 TeV was detected by HAWC (Albert et al. 2020), Tibet AS γ (Amenomori et al. 2021) and LHAASO (Cao et al. 2021); the UHE source is coincident with the VERITAS and Fermi-LAT tail region sources as well as PSR J2229+6114 (Albert et al. 2020).The UHE detection identified the Boomerang region as a PeVatron candidate, but its origin is still debated between the leptonic and hadronic cases associated with the Boomerang PWN and the SNR interaction with molecular clouds, respectively (Ge et al. 2021;Bao & Chen 2021;Fujita et al. 2021;Breuhaus et al. 2022;Liu et al. 2022).Various high energy emission centroids/extents are depicted over the 1420 MHz radio temperature brightness map of SNR G106.3+2.7 in Figure 1.The distance to the Boomerang complex is unusually poorly determined, even among supernova remnants.A list of the various distance measurements is provided in Table 1.Pineault & Joncas (2000) reported radio continuum observations with DRAO at 408 and 1420 MHz, as well as H I observations, in which they identified an absorption feature at −104 km s −1 , giving a kinematic distance of 12 kpc.This would put G106.6+2.9 at a z-height of 607 pc above the Galactic plane, with a linear extent of over 200 pc.However, Halpern et al. (2001b) suggested a much closer distance based on the measurement of an absorbing column density.They observed the PWN region with a radio (VLA) and X-ray (ASCA) telescope in order to search for a counterpart of the unidentified EGRET gamma-ray source 3EG J2227+6122.In this study, they obtained an absorbing column N H of 6.3 × 10 21 cm −2 .Since the column density through the entire Galaxy is only 8.4 × 10 21 cm −3 in that direction, they concluded that the PWN was at least 2 kpc away, perhaps much further, and assumed a fiducial distance of 3 kpc.The pulsar discovery (Halpern et al. 2001a) reported a DM of 200 ± 10 cm −3 pc, which, with the Taylor & Cordes (1993;TC93) electron-density model, implies a distance of 12 kpc (Taylor & Cordes 1993).A revised n e model, NE2001 (Cordes & Lazio 2002), gives 7.5 kpc for the same DM value (Abdo et al. 2009b).Yao et al. (2017) proposed a new n e model (YMW16) to estimate the distance to the pulsars using the same DM value.YMW16 estimated the distance to the Boomerang PWN to be 5.037 kpc with an error of 40 %.This error is larger than 20 %, the threshold YMW16 considered for their model estimation to be satisfactory.Yao et al. (2017) also notes that the distance to the Boomerang PWN showed the largest impact due to the Galactic warp.A distance of 3 kpc would imply, from TC93, a DM of only 75 cm −3 pc.Kothes et al. (2001) suggested a very near distance, 800 pc, based on morphological correspondences between the radio continuum image and channel maps of H I and CO from surveys, at velocities of about −6 km s −1 .Kothes et al. (2004) presented a new technique of H I absorption of polarized emission for distance determinations; they pointed to the absence of an absorption feature in the range of −70 to −55 km s −1 , which they asserted would be present if G106.6+2.9 were further away than the Perseus arm at about 3 kpc.
The ambiguity in the distance of G106.6+2.9 may be rooted in its relatively high Galactic latitude of 2.9 • .Over 85% of the 383 Galactic SNRs in the catalog of Ferrand & Safi-Harb (2012) are closer to the Galactic plane than this.At 3 kpc, for instance, G106.6+2.9 has a z-height of 150 pc, higher than the H I scale height of the Galactic disk of about 100 pc, and perhaps explaining anomalous H I absorption (or its absence).
The properties of G106.6+2.9 are extreme at either distance: 0.8 kpc or 12 kpc.At 0.8 kpc, as pointed out by Kothes et al. (2006), the PWN would have an extremely low ratio of radio luminosity to pulsar Ė.Additionally, the H I column density from X-ray observations of 6.3 × 10 21 cm −2 implies a mean volume density of neutral atomic H of 2.6 cm −3 between Earth and G106.6+2.9 which is unrealistically high.At 12 kpc, Halpern et al. (2001a) state that the pulsar would need to be more efficient than the Crab or Vela pulsars at converting spindown luminosity into > 100 MeV gamma-rays.
In this paper, we present a multi-wavelength analysis of the Boomerang PWN region.We shall refer to the compact radio and X-ray source within 100 ′′ of the pulsar as the PWN, and as distinct from the head (scale ∼ 15 ′ ) and tail (scale ∼ 30 ′ ) regions.Our work is focused on this PWN region, in contrast to the recent publication by Fang et al. (2022) for example, which primarily concerns the larger-scale nebula.Our motivation is to better constrain the characteristic properties and gain further insight on the formation of the Boomerang PWN.In doing so, we gain a better understanding of the Boomerang PWN's relationship to the high-energy emission coincident with SNR G106.3+2.7,mostly confined to the SNR's tail region.We begin by describing the archival Chandra and new NuSTAR X-ray observations and our timing, imaging and spectral analysis ( §2.1).In §2.2 and §2.3 we describe the gamma-ray observations of the Boomerang PWN region and analysis of the corresponding data from Fermi-LAT and VERITAS, respectively.We then combine the multi-wavelength spectral data of the Boomerang PWN and explore various models that could describe the PWN's emission through SED fitting ( §3).For the SED models, we consider the two most extreme source distances, 0.8 and 7.5 kpc, from the Table 1.(From the two distance measurements based on H I radial velocity measurements, we chose 0.8 kpc over 12 kpc as it is more the most recent estimation.)In the end, we determine the source distance from the SED and X-ray morphology Tanalysis of the Boomerang PWN region, within 100 ′′ of the pulsar; We do not consider the emission on larger scales.Finally, we discuss the results from our X-ray and multi-wavelength analysis and constrain the PWN B-field ( §4).We contemplate the current evolution phase of the Boomerang PWN and examine its relation to the high-energy emission coincident with SNR G106.3+2.7.We summarize our results in §5.

OBSERVATIONS AND DATA ANALYSIS
We present X-ray, GeV (Fermi-LAT) and TeV (VERITAS) gamma-ray observations of the Boomerang PWN in the following sections.We performed X-ray analysis of the 2002 Chandra and 2020 NuSTAR observation data ( §2.1).Fermi-LAT and VERITAS data analysis was confined only to the Boomerang PWN region, and the source extraction Pope et al.
region varied according to the point spread function (PSF) of each telescope ( §2.2 and 2.3).All errors are given to the 2σ confidence level unless explicitly stated otherwise.

X-ray observations
Chandra observed the Boomerang region with its CCD ACIS-I array on 2002 March 15 (Obs ID 2787, PI Halpern) for ∼ 94 ks of exposure.The observation files were processed and analyzed using the tools in CIAO v4.13.NuSTAR observed PSR J2229+6114 and its PWN on 2020 September 21 (ObsID 40660001002) for a total of 45 ks of exposure.Data analysis was conducted using the NuSTARDAS v2.0.0 sub-package within HEASOFT v6.28.CIAO v4.13 was also used for its image modelling and fitting application (SHERPA).
The Boomerang PWN region is composed of four components: (1) the pulsar (PSR J2229+6114), (2) a torus-jet feature which represents the termination shock region (r ∼ 10 ′′ ; Ng & Romani 2004), (3) X-ray PWN (r ∼ 30 ′′ ), and (4) diffuse X-ray emission (r ∼ 100 ′′ ).These X-ray features are resolved by Chandra as shown in the Chandra image in the 0.5-8.0keV band (Figure 3).Since NuSTAR (with 58 ′′ half-power diameter) cannot spatially resolve the pulsar from the extended X-ray emission, we first performed timing analysis on the NuSTAR data in order to remove the pulsar emission ( §2.1.1).We then present NuSTAR imaging analysis in different energy bands in comparison with the high-resolution Chandra images below 8 keV ( §2.1.2).In §2.1.3,we analyze X-ray spectral data of the PWN by excising the pulsar emission spatially (for Chandra) and by selecting a phase interval for the off-pulse component (for NuSTAR).

NuSTAR timing analysis
We found that X-ray emission from the Boomerang pulsar + nebula system was detected up to 20 keV.We limited our timing analysis to the 3-20 keV band.The NuSTAR telescope consists of two focal plane modules (FPMA and FPMB) which are described in detail in Harrison et al. (2013).We determined the source centroid for both FPMA and FPMB images using DS9's centroid function.We then applied barycentric correction to the event files using barycorr, and extracted source events from a r = 30 ′′ circular region around the source centroid.Using the Stingray X-ray timing analysis package (Bachetti et al. 2022), we applied the Z 2 algorithm with two harmonics in order to search the combined event times from both focal plane modules for pulsations.A strong periodic signal with P = 51.671495+1×10 −6 −3×10 −6 ms (3σ error bars) was detected with Z 2 2 = 174.This period differs slightly from that reported by Halpern et al. (2001a) after accounting for the Ṗ quoted in the same paper, P = 51.67199357(0),indicating an undetected timing glitch.Folding upon the measured period, we generated a pulse profile in the 3-20 keV band (Figure 2).The on-pulse component is clearly identified as an asymmetric double-peak in the folded lightcurve, consistent with the pulse profiles measured by Halpern et al. (2001a) and Kuiper & Hermsen (2015).We considered our baseline off-pulse component to be around 25 to 30 counts per bin.As a conservative estimate for the pulsed emission, we only excised the clear peaks significantly above the baseline.After calculating a phase value for each photon event, we removed the on-pulse events between ϕ = 0.17 and 0.30, as well as ϕ = 0.70 and 0.85 using extractor for subsequent spectral and imaging analysis of the  Pope et al.
nebular emission.The effect of any leftover pulsar component post phase extraction is negligible, as we determine in Section 2.1.3.

NuSTAR and Chandra imaging analysis
The Chandra observations detected the Boomerang PWN extending to the bounds of its radio emission (r ≈ 100 ′′ ; Halpern et al. 2001a).We thus considered r = 100 ′′ as the outermost boundary of the X-ray nebula for NuSTAR imaging analysis using the phase-resolved event files after excising the pulsed emission.The broad bandwidth of NuSTAR allows us to compare the X-ray image of Boomerang in different energy ranges.We performed energyresolved imaging analysis in a "soft" band (3-10 keV) and "hard" band (10-20 keV).For each of the FPMA and FPMB data, we corrected the positional offsets between the NuSTAR source centroid (measured in the on-pulse images) and the pulsar position measured by Chandra.
The event files for both FPMA and FPMB were split into the soft and hard bands with extractor.Exposure maps were created with nuexpomap for each event file with vignetting effects at 6.5 and 15 keV for the soft and hard band, respectively.The FPMA and FPMB event files were combined for each energy band, and the same was done for the exposure maps.For the purpose of smoothing out spurious features near the detector edges, the summed NuSTAR images were convolved with a Gaussian kernel of σ = 2.46 ′′ (corresponding to the NuSTAR pixel size) before being divided by the corresponding exposure maps (Nynka et al. 2014).The above process produced an exposure corrected mosaic flux image for each energy band.In each energy band, we calculated the background level using a region to the northeast of Boomerang, avoiding the diffuse X-ray emission detected by Ge et al. (2021).The resultant 3-10 keV and 10-20 keV background subtracted flux images of the Boomerang PWN are shown in Figure 5.
In order to characterize the Boomerang PWN emission, we compared radial profiles around the pulsar position in the two energy bands, as well as the NuSTAR PSF for determining a source extension.A set of 20 annuli between r in = 5 ′′ and r out = 100 ′′ were centered on the source centroid of the mosaic image.The radial profiles for each energy band were extracted from these annuli and normalized so that the brightness was set to 1 at r = 0.The same set of annuli was used to create a normalized radial profile of the NuSTAR 8-12 keV PSF to serve as a point source template.Since we found that the NuSTAR PSF varies insignificantly in 3-20 keV, we chose 8-12 keV to produce the radial profile representative for our case.The soft, hard, and PSF radial profiles are shown in Figure 6.
The background subtracted source profiles, as shown in Figure 6, are extended above the NuSTAR PSF up to r ∼ 100 ′′ Ṫhe radial profiles in both the soft and hard bands appear more extended than the NuSTAR PSF profile.Furthermore, the soft band exhibits a slightly wider radial profile than that of the hard band.To determine the size of the nebula in each energy band more quantitatively, we fit the NuSTAR images using SHERPA.We modelled the  X-ray source as a 2D-Gaussian and included a constant background level.The source model was convolved with the NuSTAR PSF and then fit to the NuSTAR image data.After taking into account the telescope dithering, we produced the effective NuSTAR PSF data in 4.5-6 keV and 12-20 keV for the soft and hard bands, respectively (Nynka et al. 2014).The fit yielded a full width at half maximum (FWHM) of 33 ± 2 ′′ for the soft band and 20 ± 2 ′′ for the hard band; the errors represent the 1 sigma confidence intervals (see Figure 6).
For the purpose of illustrating the various regions of interest in the Boomerang PWN and comparing with the NuSTAR images, we created a high-resolution radial profile of the Chandra 0.5-8 keV image shown in Figure 4. To produce this radial profile, we generated a set of 20 annuli around the pulsar position from r in = 1 ′′ (to mask out the pulsar emission) to r out = 100 ′′ (i.e the boundary of the X-ray nebular emission).The profile produced from these annuli was normalized, and the background surface brightness was extracted from the same region used for background spectra extraction.The radius of the diffuse X-ray extent measured by Halpern et al. (2001a) and the NuSTAR 3-10 and 10-20 keV PWN radii are plotted over the resulting radial profile in Figure 4.The FWHM of the soft and hard Pope et al.
X-ray emission detected by NuSTAR (r = 33 ′′ and 20 ′′ ) corresponds to the inner PWN where most of the X-ray nebular emission is concentrated.

Chandra and NuSTAR spectral analysis
In this section, we present our Chandra and NuSTAR spectral analysis of the entire PWN region from r = 100 ′′ .This region is in accordance with the Boomerang radio PWN, thus it is appropriate for subsequent SED studies (e.g., §3.1).Using the CIAO specextract script, the Chandra spectrum was extracted from the r = 100 ′′ circular region around the pulsar position, excluding the pulsar emission in r < 1 ′′ .Increasing the exclusion region to r < 3 ′′ had no significant effect on our spectral fits, suggesting that r < 1 ′′ is sufficient for excising the pulsed emission component.
A source-free region to the immediate northwest of the source extraction region was used for extracting background spectra.In contrast to the Chandra spectral analysis, we performed NuSTAR spectral analysis on the off-pulse events of the Boomerang region.We used a r = 100 ′′ circular source region around the emission centroid (which was previously adjusted to the pulsar position in each module data).The NuSTAR response matrix (RMF) and effective area (ARF) were created using nuproducts for an extended source option.We produced NuSTAR background spectra in two different approaches by modelling the background using nuskybkg (Wik et al. 2014) as well as by extraction from an off-source region with nuproducts.In the former method, we modeled background spectra with nuskybkg by fitting actual background spectra extracted from multiple source-free regions across the NuSTAR detector chips for both modules.In the latter method, a source-free region on the same detector chip as the source was used for generating   background spectra.In both cases, background regions were selected to avoid the additional diffuse non-thermal X-ray emission in the head region (Ge et al. 2021).We found that fitting the source spectra with the background spectra produced by either method yielded statistically identical results.Hereafter we present our NuSTAR spectral analysis results with the standard background extraction using nuproducts.
The NuSTAR and Chandra spectra were adaptively binned to 2.5 and 2.0 σ over background counts, respectively.In order to determine the hydrogen column density accurately, we fit the 0.5-8.0keV Chandra spectrum with an absorbed power-law model in XSPEC v12.12.0.The fit was first carried out using the default abundance data (Anders & Grevesse 1989), resulting in the best-fit column density N H = 6.2 +1.0 −0.9 × 10 21 cm −2 .We then repeated the Chandra spectral fit using the Wilms abundance table (Wilms et al. 2000), obtaining a higher value of N H = 8.9 +1.5 −1.4 ×10 21 cm −2 .Radio observations of PSR J2229+6114 measured its DM to be (204.97± 0.02) cm −3 pc (Abdo et al. 2009b).Using a linear relation between N H and DM as well as its slope errors (He et al. 2013), we derived N H = (4.3−8.8)× 10 21 cm −2 .The hydrogen column densities obtained by the Chandra observation and estimated from the DM measurement are consistent with each other.
We proceeded with spectral fitting using this updated N H value found using the Wilms abundance table.The NuSTAR FPMA and FPMB 3-20 keV spectra were jointly fit using an absorbed power-law model with independently varying FPMA and FPMB flux normalization variables, with the hydrogen column density fixed to N H = 8.9 × 10 21 cm −2 .This allowed us to find the ratio between the FPMA and FPMB flux normalization values (i.e. the cross-normalization factor).We also created a joint fit using the Chandra 0.5-8.0keV and NuSTAR 3-20 keV spectra, again with an absorbed power-law model.In this case, FPMA and FPMB were forced to share a flux normalization value, independent of Chandra's.The resulting joint fit measured the photon index to Γ = 1.52 ± 0.06 with a reduced chi-squared value of χ 2 ν = 0.95 (226 dof), consistent with the individual NuSTAR and Chandra spectral fit photon indices within error (see Table 2).A broken power-law or an additional spectral component was not statistically required given the goodness-of-fit with a single power-law model.The NuSTAR spectra and NuSTAR-Chandra joint spectra with the best-fit models are presented in Figure 7.The parameters from these fits, as well as from the Chandra spectral fit, are listed in Table 2. Given that the cross-normalization factor for the joint fit (C = 0.92 +0.09 −0.08 ) is consistent with 1, no significant X-ray flux change was detected from the PWN between the 2002 Chandra and 2020 NuSTAR observations.
In order to confirm that our NuSTAR PWN spectra was not significantly tainted by pulsar emission post phase extraction, we attempted to characterize the off-pulse pulsar emission.While this was difficult because of the low number of pulsed photon counts, we decided upon the following method.In order to characterize the on-pulse pulsar component, we extracted the on-pulse NuSTAR 3-20 keV spectra from the r = 30 ′′ region around the pulsar position and used the corresponding off-pulse spectra as background.We then jointly fit this spectra with the Chandra 2-8 keV point source spectra, using a tbabs*(pow+const*pow) model.By setting the constant term to zero for the NuSTAR spectra, to one for the Chandra spectra, and allowing the first powerlaw to only fit the NuSTAR spectra, the latter powerlaw characterized the off-pulse pulsar component.We found that the off-pulse pulsar flux component represents only about 10% of the pulsar emission, and only about 5% of the off-pulse emission.We therefore forgo redoing the above spectral analysis with a off-pulse pulsar model component; it is negligible.

Fermi-LAT data selection and analysis
The Fermi Large Area Telescope (Fermi-LAT) is a pair-conversion high-energy (HE) gamma-ray telescope, that can detect gamma rays in the energy range from 20 MeV to above 1 TeV (Atwood et al. 2009).The presented analysis was performed by means of Fermipy, a python-based package that allows to analyse Fermi-LAT data with the Fermi Science Tools (Wood et al. 2017).We used Fermipy version v1.0.1, which is associated to the 2.0.8 version of the Fermi Science Tools.
We used P8R3 SOURCE V3 instrument response functions (IRFs) and event type 3 (front and back conversion type).We binned the data using a spatial size of 0.1 • and 8 energy bins per decade.We modeled all the sources within a box of width 20 • that are included in the second release of the 4FGL catalog (4FGL-DR2; (Acero et al. 2016)), along with the isotropic and Galactic diffuse emission (iso P8R3 SOURCE V3 v1 and gll iem v07).
Pope et al.We first optimized the model by fitting normalization and spectral shape parameters of each source in the region of interest and calculate their Test Statistics (TS = −2ln(L max,0 /L max,1 ) where L max,0 and L max,1 are the maximum likelihoods for a model excluding and including the source, respectively1 ), by using the gta.optimize function.
We then simplified the model by removing sources that have TS below 4 (i.e. with a detection significance below ∼ 2 standard deviations) and number of predicted counts below 1, in order to ease convergence of the fit.Before performing the final fit, we freed sources that have TS above 25 (i.e. with a detection significance of ∼ 5 standard deviations or above), within a radius of 5 • from the target position of 4FGL J2229.0+6114, as well as the isotropic and Galactic diffuse emission components.We also modelled the emission from the tail region of the Boomerang nebula as described by Xin et al. (2019), with a uniform disk with spatial width 0.25 • centered around (R.A., decl.)= (336 • .68,60 • .88)and power-law spectrum.Figure 8 shows the spectral energy distribution (SED) of 4FGL J2229.0+6114.We used the same spectral model as the one reported in the 4FGL-DR2 catalog (PLSuperExpCutoff2 ).The SED was extracted using the gta.sed() method, which fits the flux normalization of the source in each energy bin, using a power-law with a fixed index of -2.We see no emission above 50 GeV (TS of the spectral points is below a threshold value of 4); the differential upper limit in the energy range 50.7 GeV to 2 TeV is 2.91 × 10 −7 MeV cm −2 s −1 at 95% confidence level (CL).As outlined in (Abdo et al. 2013), pulsar spectra in the LAT energy range should cut off exponentially around a few GeV; for this reason, we conservatively only considered the measurements above 50 GeV as PWN emission in order to cut most of the pulsed emission from 4FGL J2229.0+6114.

VERITAS observation and analysis
VERITAS is an array of four imaging atmospheric Cherenkov telescopes (IACTs), located near Tucson, Arizona, designed to measure gamma rays of energies from 100 GeV up to > 30 TeV.Each telescope has a field of view of 3.5 • , and the array can detect a point-like source with 1% of the Crab PWN flux at 5σ significance within 25 hours.VERITAS has an angular resolution of ∼0.1 • (Park & VERITAS Collaboration 2015).The VERITAS Collaboration previously reported the detection of gamma-ray emission from the region of SNR G106.3+2.7 with 33.4 hours (Acciari et al. 2009) based on data collected in the 2008 epoch.The TeV emission was observed near the center extended radio emission (see 1) rather than the location of the Boomerang PWN.From 2009 to 2010, VERITAS accumulated an additional 22.3 hours with a changed array configuration where one telescope was moved to make the array more symmetric, which improved the sensitivity of the array (Perkins et al. 2009).Combined with the previous data set, we used a total exposure of 57.7 hours for the analysis at the location of the Boomerang PWN.As the extension of the PWN measured in radio and X-ray is smaller than the angular resolution of VERITAS, the analysis was performed with an assumption that the emission is a point-like source.Standard VERITAS analysis was performed with two independent analysis methods.Two different event selections were used for the analysis: one selection was optimized to search for emission that was 2-10 % of the Crab Nebula strength and one selection was optimized to search for emission weaker than 2% of the Crab Nebula strength.The event selections optimized to search for the weaker emission reject the largest fraction of background events, resulting in a higher energy threshold.No strong gamma ray emission was detected at the location of the Boomerang PWN.Upper limits at the 99% confidence level for two different energy thresholds were calculated based on the assumption of a power-law spectral energy distribution with a spectral index ranging from 2 to 3. The results are shown in Table 3.These upper limits are shown in the Figure 9 together with the SED of known VHE gamma ray emission in the tail region of SNR G106.3+2.7.

BROADBAND SPECTRAL ENERGY DISTRIBUTION
In this section, fer modelling of the multi-wavelength data from the radio to TeV band.We present analyses using three different leptonic SED models, all one-zone (homogeneous) models with varying degrees of spectral detail.As there was no detection of gamma-ray emission above 3 GeV from the Boomerang PWN, we will not consider the effect of adding hadronic components to the model in this paper.We will include hadronic components in a future paper reviewing the SEDs from the entire SNR.The morphological complexity of the Boomerang PWN and larger-scale system is such that all three models used in this paper are highly simplified.And as can be seen in Figures 3 and 5, the X-ray centroid in both the soft and hard bands is offset from the peak in the radio band, suggesting the possibility that Boomerang is a multi-zone system.By using these one-zone models we assume that both the X-ray and radio

Pope et al.
emission originate from the same source and physical processes as part of a single system.This assumption has been made for other PWNe with similar offsets between X-ray and radio peak emission, such as the PWNe associated with SNRs G54.1+0.3,G327.1+1.1 and MSH 15−56 (Lang et al. 2010;Temim et al. 2015Temim et al. , 2013)).Despite the simplification made in using these aforementioned models, we hope to obtain general constraints on important parameters which could guide more complex models specialized to the unique characteristics of the Boomerang system.For the radio band, we adopted the flux data from Kothes et al. (2006).We used the Chandra and NuSTAR X-ray spectra after correcting for ISM absorption.Since we applied one-zone SED models, both radio and X-ray spectra were extracted from the same r = 100 ′′ region around the pulsar position.The Fermi-LAT flux upper limits were derived from the analyses described in §2.2.As the difference between the VERITAS 0.38-30 TeV and 0.72-30 TeV upper limits found in §2.3 is small, we chose to use the former value.Note that the source extraction regions for Fermi-LAT and VERITAS analyses are subject to the telescope PSF sizes which are larger than r = 100 ′′ .Throughout the SED modelling, we consider two contrasting source distances of 0.8 and 7.5 kpc suggested by the H I velocity (Kothes et al. 2001) and the pulsar's DM measurement (Abdo et al. 2009b), respectively.Below we present three different SED models.NAIMA models radiative SEDs for a given time-independent electron energy distribution ( §3.1).Although this is a simplistic approach, we attempt to constrain several PWN parameters such as B-field and electron spectral index.In §3.2, we applied the time-dependent SED model package GAMERA (Hahn 2016) to the multi-wavelength SED data.This time-dependent approach, which assumes PWN free expansion, revealed several challenges in fitting the PWN SED data and size simultaneously, and further constrained multiple PWN properties.However, both NAIMA and GAMERA proved to be overly simplistic.As can be seen from visual inspection of the radio and X-ray flux data, a power-law fit to the X-ray data will undershoot the observed radio data, a phenomena seen for only a few other PWNe (see Hattori et al. (2020) for example).The best fit results from both NAIMA and GAMERA confirm the difficulty in reproducing both the radio and X-ray data with such simplistic models.We therefore also considered the PWN evolution using the more complex, dynamical SED model developed by Gelfand et al. (2009) in §3.3.The model has been widely used for modelling SED data of various PWNe including the Crab nebula, G21.5−0.9 and composite SNR-PWN systems (e.g., Hattori et al. (2020)).The dynamical PWN model allowed us to track a full evolution path from the free expansion, SNR reverse shock interaction, and re-expansion phases.Both the SED and PWN radius are modeled as a function of time in comparison with the observation data.In this physically motivated approach, we determine the current B-field, pulsar's true age, expansion velocity, and its current evolution stage.

NAIMA SED model
In order to estimate initial PWN parameters from the SED data, we relied on the NAIMA V0.10.0 python package (Zabalza 2015).NAIMA is a time-independent, one-zone SED model used to generate multiple radiative model components from an assumed particle energy distribution.We fit the multi-wavelength SED data assuming that the electron distribution is in the form of a power-law model A(E e /E 0 ) −p between E e = E min and E max .While a hard cutoff at E e = E min is not physically motivated, it is implemented to simplify the model.We generated leptonic radiation models with synchrotron radiation, ICS and synchrotron self Compton (SSC) components.We adopted the cosmic microwave background (CMB) as a seed photon source for the ICS component.No IR seed photons were added so as to maintain the ICS component of the model as a lower limit.Furthermore, including a variable IR seed photon component would only add to the model's degeneracies, an already significant issue which we discuss below.We consider an IR photon field in the more elaborate dynamical model in Section 3.3.While the initial physical parameter estimates of the emitted plasma provided by NAIMA are useful, the model does not directly provide any insight as to where this emitting plasma came from, and -due to the evolution of the PWN inside the SNR, especially once it collides with the reverse shock -the particle spectrum is unlikely to be well described by a single or broken power-law.Our fitting results confirm NAIMA's inadequacy in describing the Boomerang PWN's particle spectrum.
We first focused on reproducing the radio spectral break at ∼ 5 GHz and X-ray spectra.We found that the radio data are adequately fit by various sets of B and E min values, as listed in Table 4.We have shown that the radio spectral break can be caused by the selected E min .The origin of the break in the electron spectrum, used here to fit the radio break, cannot be determined with the model.Given the degeneracy between B and E min , we did not find B = 2.6 mG as a unique solution to reproduce the radio data, contrary to the results of Kothes et al. (2006).On the other hand, the electron spectral index was constrained to p = 2.5 by fitting the radio and X-ray spectra with a synchrotron model component.We expect that (p − 1)/2 = Γ x − 1.We therefore infer from p = 2.5 that Γ x = 1.75, significantly larger than any of the X-ray photon indices found in Table 2, confirming that NAIMA cannot adequately fit both the radio and X-ray data.However, the parameters serve as an initial estimate.For each assumed source distance (0.8 & 7.5 kpc), a list of the model parameters for four representative cases (including B = 2.6 mG) is shown in Table 4 and two representative SEDs are plotted in Figure 10.In the gamma-ray band, the ICS and SSC components remain below the Fermi-LAT and VERITAS flux upper limits as long as B > 5 µG.This is because, as the magnetic field is decreased, a larger number of electrons is required to fit the radio and X-ray data and thus enhances the model GeV-TeV flux through the ICS component.If we consider other seed photon sources than the CMB, the lower limit will be higher than 5 µG.NAIMA outputs the total energy of the electron population, U E , integrated between E min and E max .Assuming the B-field is uniform throughout the interior of the PWN (r < 100 ′′ ), we calculated the total magnetic field energy as PWN is the volume of the PWN.The total electron and magnetic field energies for each fit are provided in Table 4.Given the constraints on the synchrotron radiation fluxes, increasing the magnetic field from B ∼1−10 µG to B > ∼ 100 µG requires decreasing the electron density dramatically.As a result, a fraction of the total energy allocated to PWN B-field (η B ) is ∼ 1 for B > ∼ 100µG.Such a high magnetization parameter is unusual if compared to other PWN systems -e.g., Martin et al. (2014) found η B = 7 × 10 −4 −0.02 from six PWNe including the Crab nebula (which possesses the highest η B value).If we attain an η B value comparable to those deduced from the other PWNe, it points to a lower B-field range of B < ∼ 10 µG.While these arguments as well as the suggested B-field range (B∼5−10 µG) are empirical based on the time-independent one-zone SED model fitting, we will address this issue with more advanced SED models and other X-ray analysis results in the later sections.
Pope et al.

GAMERA SED model
To further characterize the PWN properties, we modelled the time evolution of the particle distribution and radiation SED using GAMERA (Hahn 2016).GAMERA allows us to inject leptons and track their energy distribution as they vary via radiative and adiabatic cooling over time.However, the GAMERA modelling does not account for the interaction between the PWN and the parent SNR which makes the PWN compress and expand, altering the injected particle distribution significantly.The model assumes that the PWN is expanding with a constant velocity, and thus only tracks the recent PWN evolution.The age reflects the time since the PWN started expanding, whether that be from the time of the pulsar's birth or from the time that the reverse shock passed through the PWN.We assumed R PWN increased linearly with time over the short lifetime of the new PWN and matched its current radius (at t = t PWNage ) with the measured radio PWN size.Assuming two different sources distances (0.8 and 7.5 kpc), we adopted the corresponding PWN sizes of r = 0.4 pc and 4 pc.We assumed continuous particle and magnetic field injection into the PWN at the rate of the pulsar's current spin-down power Ė = 2.2 × 10 37 ergs s −1 .A fraction of the spin-down power (η B ) is injected into the PWN in the form of magnetic potential energy U B .The magnetic field B was assumed to be uniform throughout the spherical volume of the PWN.The rest of the spin-down power is injected into particle energy following a power-law distribution A(E/E 0 ) −p .The normalization factor A is determined by the total particle energy of (1 − η g − η B ) Ė where η g is the gamma-ray efficiency of the pulsar.We ignored the fraction of the spin down energy allocated to gamma-ray pulsations; this contribution should not significantly effect the model outputs.The injected particle distribution is evolved by three particle cooling mechanisms, i.e. adiabatic, synchrotron and ICS components.The adiabatic cooling rate is calculated using the expansion velocity of the PWN, which we assumed to be constant (v PWN = R PWN /t age ).We then calculated a radiation SED from the particle distribution at t = t age .We employed a leptonic radiative model comprised of emission from synchrotron radiation and ICS of the CMB and synchrotron-self-Compton scattering emission.Similarly to our NAIMA fitting, we have only considered ICS of CMB and SSC components into account to maintain the ICS component of the model as a lower limit.Still, GAMERA is not able to describe the data very well as summarized below.Adding other potential seed photons, such as IR, would only make fitting the data more challenging for these models, as the IC component would encroach on the Fermi-LAT and VERITAS upper limits.The addition of an IR component is considered for our dynamical model in Section 3.3.
We also attempted to constrain our model based on the lack of X-ray variability over the last ∼ 20 years ( §2.1.3).Since the luminosity of synchrotron radiation is proportional to B 2 , a sizable variation in B-field can lead to detectable X-ray flux variability.As the PWN evolves, decreasing B-field and increasing particle density cancel out with each other and result in a slower synchrotron luminosity evolution, while ICS luminosity keeps increasing as more particles are injected over time.While fitting the GAMERA model to the multi-wavelength SED data, we tracked the synchrotron X-ray luminosity L syn (t) evolution and found solutions without significant X-ray variability over the last ∼ 20 years.7.9 10.0 2.1 4.5 Electron energy [erg] 8.9 × 10 45 2.9 × 10 47 Magnetic field energy [erg] 1.3 × 10 42 6.3 × 10 45 Expansion velocity [km/s] 1.3 × 10 4 3.9 × 10 3 a The gamma-ray efficiency of PSR J2229+6114 from Abdo et al. (2009b).fΩ (Eq.( 3) in the cited work) was assumed to be 1.The best-fit parameters are listed for the short and long distance cases in Table 5. Figure 11 and 12 show a time series of the particle SED, radiation SED, and magnetic field, along with the radiation SED at the current time using the best-fit parameters.Note that both GAMERA models grossly misrepresent the spectrum in the X-ray region, and have other shortcomings as discussed below.
(1) d = 0.8 kpc case: The lifetime of the PWN should be t age ∼ 30 yr for d = 0.8 kpc.If t age < ∼ 20 yr, the model predicts that the 2002 Chandra observation should have detected a much higher X-ray flux than it did.If the lifetime is longer than τ ∼ 30 yr, too many leptons are injected into the PWN to be consistent with the GeV and TeV flux upper limits.Currently, the B-field should be as low as B ∼ 2 µG (with a very low magnetization of η B = 7 × 10 −5 ).Otherwise, the model over-predicts synchrotron radiation fluxes in the radio and X-ray band.We also found that E min is well constrained by the radio spectral break and the overall flux normalization of electrons (which depends sensitively on E min for a given Ė and lifetime).The expansion velocity should be high (V PWN = 1.9 × 10 4 km/s) in order to reach the observed PWN radius within ∼ 30 yr.To relax these stringent constraints, we introduced an extra, non-radiative energy loss term in GAMERA via escaping leptons.Only if we assumed a short particle escaping time (t esc < 0.05 yr), were we able to fit the SED data with higher η B values and older ages.However, such a short particle escape time is unrealistic as it is shorter than the light crossing time of the PWN.
(2) d = 7.5 kpc case: In contrast, the larger distance allows the PWN to radiate away its rotational energy over a more extensive period of time.For example, over 1,000 yr significantly more particle energy can be injected, totalling 3.1 × 10 47 ergs.While the injected electron spectral parameters such as p, E min and E max are similar between the two distance cases (Table 5), η B is higher at 1 × 10 −2 while the current B-field is 4.5 µG.The expansion velocity is 3.9 × 10 3 km/s, which is more in line with the observed PWN velocities (v PWN ∼ 1, 000−1, 500 km/s) from the Crab nebula and Kes-75 (Bietenholz et al. 1991;Reynolds et al. 2018).Overall, we found that the larger distance allows for an older age (i.e.longer particle injection time) and higher magnetization parameter.In either case, the current B-field needs to be as low as B ∼ 5 µG in order to match the low radiation efficiency (L syn / Ė).Consequently, the radio spectral break cannot be caused by synchrotron cooling at B = 2.6 mG as suggested by Kothes et al. (2006).Alternatively, we found that the radio break energy is directly related to the minimum energy (E min ) in the injected particle energy distribution.Unlike the NAIMA model, the added complexities of the GAMERA model broke the degeneracy between B and E min .In both cases, E min ∼ 10 GeV fit the radio SED data well, and thus we attribute the break to the PWN's intrinsic particle injection distribution.

Dynamical PWN evolution model
In this section, we explore the time evolution of the Boomerang PWN using the dynamical PWN model (Gelfand et al. 2009).The SED model takes input parameters for the PWN, SNR and its environment.This model evolves a homogeneous spherical bubble of relativistic electrons and magnetic field, injected according to the pulsar spin-down luminosity and its evolution, following the dynamics of its expansion into first the expanding ejecta of a spherical SNR, and including the eventual compression by the returning reverse shock and subsequent expansion into the interior of a Sedov blast wave.More details on the model description and applications to other PWNe can be found in Gelfand et al. (2009); Gelfand (2017); Hattori et al. (2020); Burgess et al. (2022).The physically motivated model tracks the time evolution of particle energy distribution, radiative SED and PWN properties (e.g., B and R PWN ) by considering particle injection and energy loss due to radiative and adiabatic cooling at each time step.The size and bulk velocity of the PWN are calculated by the pressure balance between the pulsar wind and SNR ejecta.Setting up each model run begins with determining the pulsar's properties.Given the observed current spin-down power Ė (= 2.2×10 37 ergs s −1 ) and characteristic age t ch (= 11 kyr), we first derive the system's true age (t age ) and initial spin-down luminosity Ė0 where p and τ sd are the pulsar's braking index and spin-down timescale, respectively.These input parameters fully characterize the pulsar as a particle injection source.
In the model, the pulsar injects leptons and magnetic field at each time step by partitioning the time-dependent spin-down power Ė(t) into (1 − η B ) Ė(t) and η B Ė(t), respectively.The allocated electron energy is distributed between E e = E min and E max following a broken power-law model.The evolution of both SNR forward and reverse shocks is separately calculated by going through the free expansion and Sedov-Taylor phases.The density profile of SNR ejecta, which follows ρ ej (r) ∝ r −9 until reaching the density of the ISM, is used to calculate the pressure balance between the pulsar wind and ejecta.The ISM density effects the timescale of the SNR evolution and SNR reverse shock.After the SNR reverse shock hits the PWN, the ISM density plays an important role as it imposes additional pressure on the pulsar wind.These pressure factors determine the size and bulk velocity of the PWN at each time step.The injected leptons lose their energy via radiative and adiabatic cooling as the PWN expands over time.The radiative cooling in the model takes into account synchrotron and ICS components which are provided as radiative SED model output.At some point, the SNR reverse shock reaches the PWN and compresses it to the point where the pulsar wind and reverse shock are in a pressure equilibrium before the PWN begins expanding again.
We aimed to reproduce both the multi-wavelength SED data and PWN size with the dynamical model.We note that the flux upper limit obtained by VERITAS changes depending on the assumed spectral index (see Table 3).For the dynamical SED fitting, we used the VERITAS upper limit at the decorrelation energy because the sensitivity of the limit to changes in the spectral index is lowest at this energy.As shown in Figure 13, we considered the upper limit optimized for 2% of the Crab Nebula strength for the model fitting, where the decorrelation energy is 1.12 TeV.As we did for the NAIMA and GAMERA model fitting, we considered the case for both d = 0.8 and 7.5 kpc below.We initially only use CMB as a seed photon source for the ICS component; we then test the effect of adding an IR field to the model.We recognize that the detailed radio and X-ray structure of the PWN (pulsar with small X-ray nebula in the interior of the boomerang-shaped radio arc) is not well represented by a homogeneous sphere, but as with the two previous models, we hope to obtain some general insight into the possible nature and evolution of the PWN with this tool.
(1) d = 0.8 kpc case: Given the very low radiation efficiency (L radio / Ė = 3 × 10 −8 and L X / Ė = 7 × 10 −6 ), it is extremely difficult to allocate the pulsar's expended rotational energy without overshooting the flux data.To suppress the synchrotron radiation in the radio and X-ray bands, we need to minimize the number of injected leptons and PWN B-field.We were able to satisfy both the low radiation efficiency and compact PWN size (r = 0.4 pc at d = 0.8 kpc) with the following (Case A) SED model.For d = 0.8 kpc, any free-expansion phase solution that matched the flux data largely overestimated the PWN size.We therefore had to consider a re-expanding PWN after its interaction with the SNR reverse shock.We were able to roughly reproduce all the radio, X-ray and gamma-ray fluxes by allocating the particle energy to the unobserved energy bands such as MeV and > 100 TeV energies.We arrived at a very low η B and high ISM density.The predicted PWN radius (r = 0.35 pc) is roughly consistent with the Boomerang PWN size.However, the predicted X-ray spectral shape is not consistent with the NuSTAR data and also requires an extremely high braking index (p = 5.6) or very young pulsar age (t age = 640 yr).While the Case A model almost works for reproducing the observation data, we do not consider it compelling due to the unreasonable parameter values.For example, the SNR's interaction with the high ISM density in Case A (160 cm 3 ) at a distance of 0.8 kpc and age of 640 years would produce the brightest thermal SNR ever observed by a large factor.Under the assumption that the PWN confines all of the leptons injected over its entire lifetime, it is nearly impossible to model the observed PWN's faint emission if the source distance is 0.8 kpc.
(2) d = 7.5 kpc case: Despite the larger distance, we still find it constraining to fit the SED data and PWN size simultaneously.Below we present two different cases: (1) a young PWN in the free-expansion phase (Case B) and (2) a re-expanding PWN after SNR crush (Case C).Both cases have the same number of free parameters.In both cases, we need to evolve the PWN size to match r ∼ 4 pc while keeping the number of injected leptons to the level required for fitting the radio, X-ray and gamma-ray SED data.The 3rd column in   4. The SED plots for Case A, Case B, and Case C are shown in Figure 13.The parameters p1 and p2 are the particle indices below and above E break , respectively.
panel) show the model parameters and SED plot, respectively, for Case B. This case assumes that the PWN has been expanding over 1.8 kyr, and it is similar to the GAMERA model for d = 7.5 kpc.As Figure 14 (upper panels) shows, the current B-field is 3.7 µG while the PWN expanded to r = 3.8 pc.Note that the power-law spectral index is softer below than above E break in Case B, which is highly unusual, but has been seen in a few other cases (Hattori et al. 2020;Temim et al. 2015).The spectral break at the highest radio flux point in Model B is caused by the minimum energy of the particle injected at the termination shock.Below this energy, all of the particles were injected earlier and cooled to this regime.Above this energy, there are a mix of particles injected earlier which have cooled, and freshly injected particles that haven't had time to cool yet.The energy spectrum of these two particles populations are different, which leads to a change in slope at the energy separating the two.As can be seen in Figure 13, Case B fits the X-ray data well, while a clear trend away from the model can be seen in the residuals plot in the radio band.Since the emission in the radio band is dominated by low energy particles, this discrepancy between data and model may be due to our simplified assumption that E min stays constant over time.In contrast to Case B, Case C represents the SNR crush scenario, as was suggested by Kothes et al. (2006).We evolved the PWN through the reverse shock crush into re-expansion over t age = 2.1 kyr, as indicated by the PWN radius plot over time in Figure 14 (lower panels).
The current B-field and PWN radius are 2.5 µG and 2.9 pc, respectively.As seen in Figure 13, the model fits the data well.While slightly overshooting the VERITAS upper limit, factoring in the ∼ 20% uncertainty in the upper-limit of the TeV photon density at the decorrelation energy due to the unknown photon index in this band, the model is consistent with the non-detection by VERITAS.Given how close the minimum and break particle energies are, the model favors a single power-law injection spectrum with p ∼ 2.3 (i.e., standard Fermi).

Pope et al.
Reproducing the gamma-ray emission of many PWNe often requires background photon fields in addition to the CMB -either from surrounding dust (e.g., G21.5-0.9;Hattori et al. 2020) and/or nearby stars (e.g., HESS J1640-465 and Kes 75; Abdelmaguid et al. 2022;Straal et al. 2023).Since the gamma-ray emission from the Boomerang has not been detected, the presence of additional background photon fields cannot be directly constrained by the modeling above.Therefore, to determine their possible effect on the parameters derived above, we modeled the SED and dynamical properties of this source assuming an additional photon field with a temperature T = 30 K and energy density 5x that of the CMB, typical values for warm dust in these systems (Cox et al. 1986;Torres et al. 2013).As shown in Table 6, for Case B' we can reproduce the properties of the Boomerang for a set of model parameters similar to that of Case B (where the only background photon field is the CMB).However, in our parameter exploration we were not able to reproduce the observed properties of the Boomerang, assuming this additional background photon field, for a set of parameters similar to Case C in Table 6.This does not necessarily exclude Case C as a reasonable description for the Boomerang, since there are many regions in the Galaxy where the energy density of dust emission is less than that of the CMB (e.g., G21.5-0.9;Strong et al. 2000).
Overall, the PWN evolution model does not reproduce both the SED data and PWN size with reasonable parameters if we assume d = 0.8 kpc.In all of the SED models presented in this section, we found B ∼ 2−4 µG.Alternatively, we found that an extremely magnetized PWN with η B ∼ 1 can fit the SED data well.More specifically, when we adopt η B = 0.99, only 1% of the pulsar's rotational energy is allocated to particle injection and the current B-field is ∼ 100 µG.The smaller number of injected leptons and higher B-field cancel with each other to fit the synchrotron SED in the radio and X-ray bands.However, such a high magnetization parameter is unusual compared to those of six other PWNe (η B = 7 × 10 −4 −0.02) including the Crab nebula (Martin et al. 2014).For the case of d = 7.5 kpc, we consider Case C more plausible since Case B and Case B' suggest that the high-energy particle index is greater than the low energy particle index.Furthermore, in Case C the PWN interacts with the SNR reverse shock, which could explain the offset between the radio and X-ray peak emission.No explanation for this offset can be inferred from Case B or Case B'.

DISCUSSION
We discuss constraining the properties of the Boomerang PWN and its recent evolution based on the X-ray and multiwavelength observations.In the previous sections, we examined the hypothesis proposed from the radio observations (Kothes et al. 2006) that the Boomerang is a highly magnetized PWN (B = 2.6 mG) crushed by an SNR reverse shock 3,900 yr ago, which was made under the assumption that the radio break in Boomerang's spectrum is due to synchrotron cooling.However, the radio break is not necessarily caused by synchrotron cooling, which may explain the discrepancy between the results in this paper and those from Kothes et al. (2006).
Below, we present some implications for this composite SNR-PWN system, especially in the head region, and suggest future observations to further elucidate the origin of the head-tail morphology and UHE emission.

Constraining PWN magnetic field
As shown in §3, our SED study strongly suggests that the current B-field should be B < ∼ 3µG, otherwise the observed synchrotron fluxes will be significantly over-predicted.Combined with the GAMERA SED model results in §3.2, below we consider the energy-dependent X-ray size measurements from the NuSTAR observation and constrain the PWN B-field.The smaller X-ray size (e.g., r = 20 ′′ in 10-20 keV) compared to the radio size (r = 100 ′′ ) is often observed for other PWNe (Coerver et al. 2019).The radio nebula size usually reflects the PWN size determined by the particle flow evolution over the pulsar's age.On the contrary, the smaller X-ray size is determined by synchrotron cooling time.The synchrotron cooling time for electrons emitting X-rays of energy E [keV] is τ syn = 1.2B −3/2 mG E −1/2 keV yr (Reynolds et al. 2018).Since the synchrotron cooling time depends on electron energy, we expect X-ray PWN size to shrink with photon energy.The synchrotron burn-off effect is evident, as we measured different X-ray sizes in two energy bands -33 ′′ (3-10 keV) and 20 ′′ (10-20 keV).Assuming a constant flow velocity over the X-ray synchrotron cooling time, the ratio of the PWN radii between the hard and soft energy bands should be approximately equal to the ratio between the synchrotron cooling times for the respective mid-band energies.Subsequently, adopting the mid energy in each band (6.5 and 15 keV), we estimate that the PWN size ratio between the two energy bands should be 1.5 if the flow velocity is constant, which is close to the observed ratio of 1.7.Hence, we consider the constant flow velocity hypothesis as viable.
Assuming a source distance of 7.5 kpc (as suggested by the pulsar's DM measurement and supported by our SED model fitting), the X-ray angular sizes in the 3-10 and 10-20 keV bands correspond to PWN radii of r ∼ 0.7 pc and ∼ 0.4 pc, respectively.As expected, these radii are significantly smaller than the PWN radius assumed for the SED modelling (∼4 pc at d = 7.5 kpc), which represents the size of the PWN in the radio band.For a given B-field, we can derive an upper limit of the PWN size assuming that leptons moved outward at the speed of light during their lifetime.Using the synchrotron lifetime for the highest energy X-ray emission (E γ = 20 keV), the hard X-ray PWN size of r = 0.4 pc sets an upper limit of B < 0.35 mG, which is lower than the B-field value suggested by Kothes et al. (2006).We can further constrain B-field with a more realistic estimate of the PWN flow velocity.For example, multi-epoch Chandra X-ray observations of Kes-75 measured an expansion of the PWN at V PWN ∼ 1, 000 km s −1 or 3.3 × 10 −3 c (Reynolds et al. 2018).The Crab nebula is known to expand at a velocity of V PWN ∼ 1, 500 km s −1 (Bietenholz et al. 1991) and its radius increases almost linearly with time -R(t) ∝ t 1.264 (Bietenholz & Nugent 2015).If we adopt these flow velocities, the hard X-ray size of the Boomerang PWN in the 10-20 keV band suggests B ∼ 7−10 µG.We also estimated a range of the flow velocity using the output PWN size evolution data (i.e., R(t)) provided by the dynamical PWN model ( §3.3).We found v PWN ∼ 2, 400 km/s and ∼ 400 km/s over the last 60 years of expansion for Case B and C, respectively, and their PWN radius evolution is shown in Figure 14.These PWN velocities yield B ∼ 6−14 µG, which is slightly higher than the B-field suggested by the dynamical model fitting in §3.3.
Our estimates for the PWN B-field, based on the PWN's X-ray spatial properties and SED fitting, are not only significantly below what was suggested by Kothes et al. (2006) but also B = 140 µG, which was suggested by Liang et al. (2022) based on modelling of the diffusion of relativistic electrons in the PWN.Our disagreement with the B-field derived from the analysis done by Liang et al. (2022) is expected given our different assumptions and models.While Liang et al. (2022) incorporated UHE tail emission not necessarily related to the PWN itself within their SED modelling, this paper focuses solely on emission within the general bounds of the PWN region.Furthermore, unlike the dynamical model used in this paper, the PWN model used in Liang et al. (2022) does not include any interaction with the SNR, and the temporal evolution consists only of changes resulting from diffusion and radiative losses.

PWN evolution
The dynamical SED model fitting §3.3 suggests that the Boomerang PWN is currently re-expanding to r ∼ 3 pc after being crushed by a SNR reverse shock ∼ 1000 yr ago, much shorter than the pulsar's ∼ 10 kyr characteristic age.As can be seen in the Case C scenario of the bottom panels of Figure 14, we note that the PWN is undergoing a small second compression due to an overshoot in its re-expansion, which resulted in a negative pressure differential between the PWN's interior and surrounding material.The Boomerang PWN is powered by a population of fresh electrons injected over the last ∼ 1000 yr.After the PWN compression amplified B-field to B > ∼ 10 mG, its B-field decreased, as of present day, to B ∼ 3 µG as a result of the nebula expansion.It is interesting to note the consequence of this temporally varying B-field on the spatial variations of the PWN's spectrum.Particles are generally older further away from their source pulsar (Temim et al. 2015).Subsequently, if the B-field was higher in the past as suggested by the dynamical model, leptons further away from PSR J2229+6114 would generally have experienced more significant synchrotron losses as compared to those found closer to PSR J2229+6114.Therefore, this decrease in B-field over time may explain the significant spectral softening in the X-ray spectrum shown to occur with increasing distance from the pulsar by Ge et al. (2021).
Hydrodynamic (HD) simulation of SNR-PWN interaction exhibited head and tail-like features in the particle density maps (Figure 2 in Kolb et al. 2017).According to this simulation study, there are two factors that determine the composite SNR-PWN morphology -(1) density gradient and (2) pulsar proper motion.
As suggested by Kothes et al. (2006), the density gradient around the head region likely caused their highly asymmetric morphology and Boomerang-like PWN shape in the radio band.While the X-ray bright PWN (r ∼ 0.4−0.7 pc) should be powered by a population of young electrons injected over the last ∼ 1000 yr, the relic electrons injected in prior to the SNR crush should still be producing synchrotron radiation in the radio and X-ray band.We would associate the relic PWN radiation with the head region, while we attribute the lack of TeV emission in the PWN region to a small number of electrons injected by the "refreshed" PWN.Since the total number of relic electrons injected over a few thousand years is much higher, the head region should produce higher TeV emission via ICS than in the PWN.The reported TeV detection in the head region by MAGIC may indicate such ICS emission.The UHE emission in the tail region cannot be caused by relic electrons which should have cooled quickly to GeV-TeV energies.Instead, particle re-acceleration during the SNR-PWN interaction via PWN compression, as proposed by Ohira et al. (2018), may be responsible for producing the UHE emission in the tail.A further study in the entire region, with new VERITAS observation data and more extensive SED study, will be presented in our future paper.

Figure 1 .
Figure 1.CGPS 1420 MHz radio temperature brightness map [K] of the SNR G106.3+2.7 region with the head, tail and PWN indicated by green dashed lines.The pulsar location is marked by the green cross.The white ellipse represents the extent of the gamma-ray emission previously detected by VERITAS.The black plus, yellow cross, and cyan diamond represent the centroids of the gamma-ray emission detected by HAWC, LHASSO, and Fermi-LAT, respectively.

Figure 2 .
Figure 2. NuSTAR 3-20 keV folded lightcurve of PSR J2229+6114.The pulsed phase ranges excised from our PWN imaging and spectral analysis are demarcated by the red regions.

Figure 3 .
Figure 3. (top) Chandra 0.5-8 keV image of the Boomerang PWN, smoothed to σ =∼ 0.5 ′′ (1 pixel).The color scale is in units of counts per pixel.A r = 1 ′′ circular region around the pulsar was excised in order to accentuate the PWN features.The solid green circular region is a r = 100 ′′ circular region around the pulsar position, used for spectral extraction.(bottom) 0.5-8 keV Chandra image (same as from top panel) is shown here in green, overlain with the CGPS 1420 MHz radio temperature brightness map of the Boomerang region in red.

Figure 4 .
Figure 4. Chandra 0.5-8 keV radial profile of the Boomerang PWN.The vertical lines indicate the PWN radius as measured from the diffuse X-ray emission (green), NuSTAR 3-10 (red) and 10-20 keV (blue) energy bands.The 1 sigma errors are indicated by the shaded regions.The error bars designate the ∼5 ′′ radial bin widths.

Figure 8 .
Figure 8. Spectral energy distribution of 4FGL J2229.0+6114.The red points and upper limits were derived from the analysis of LAT data presented in this work.The red solid line indicates the best-fit model to the data points, while the black dashed line is the same spectral model, with the best-fit parameters reported in the 4FGL catalog.The inset reports the residuals as (data-model)/model, for the best-fit model derived from this work and the 4FGL catalog one.The green line denotes the threshold at which we separate the pulsar (left) from the nebula (right) emission.

Figure 9 .
Figure 9. Upper limits of VERITAS at the location of Boomerang PWN with two event selections.VHE gamma-ray emission near the tail region of SNR 106.3+2.7 measured by the IACTs, VERITAS (Acciari et al. 2009) and MAGIC (MAGIC Collaboration et al. 2023), is shown as comparisons.

Figure 11 .
Figure 11.Time evolution of the lepton (top left) and radiation (top right) SEDs, and B-field (bottom left), and the radiation SED at the current time (bottom right) assuming d = 0.8 kpc.

Figure 13 .
Figure 13.Top left panel: PWN evolution SED model for case A assuming d = 0.8 kpc.Top right panel: Case B for d = 7.5 kpc.Bottom panel: Case C for d = 7.5 kpc.The model parameters for these cases can be found in Table4.

Figure 14 .
Figure 14.Top panel: Electron evolution for Case B (left) and Case C (right) of Figure 13 and Table 6.Bottom panel: Rpwn and magnetic field evolution of Case B (left) and Case C (right).

Table 1 .
Boomerang distance estimates

Table 2 .
Boomerang PWN spectral fitting parameters.All errors are given to the 90% confidence level.

Table 4 .
Figure10.NAIMA simple power-law SED fits with B-field set to 2.6 mG (left) and 5 µG (right).The complete set of parameters for each fit, as well as for other B-field scenarios, is shown in Table4.Fit parameters from various leptonic NAIMA model cases.
Distance [kpc] B [µG] p Emin [GeV] Emax [TeV] Electron energy [erg] Magnetic field energy [erg] a a ηE and ηB are calculated by dividing the electron and magnetic field energy by their sum, respectively.

Table 5 .
Model parameters for the case presented in §3.2 Table 6 and Figure 13 (upper right

Table 6 .
Model parameters for the four cases presented in §3.3