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

I. Pope  Columbia Astrophysics Laboratory, 550 West 120th Street, New York, NY
10027, USA  K. Mori  Columbia Astrophysics Laboratory, 550 West 120th Street,
New York, NY 10027, USA  M. Abdelmaguid  NYU Abu Dhabi, PO Box 129188, Abu
Dhabi, United Arab Emirates  J. D. Gelfand  NYU Abu Dhabi, PO Box 129188, Abu
Dhabi, United Arab Emirates  S. P. Reynolds  Physics Department, NC State
University, Raleigh, NC 27695, USA  S. Safi-Harb  Department of Physics and
Astronomy, University of Manitoba, Winnipeg, MB R3T 2N2, Canada  C. J. Hailey
Columbia Astrophysics Laboratory, 550 West 120th Street, New York, NY 10027,
USA  H. An  Chungbuk National University, Chungdae-ro 1, Seowon-Gu, Cheongju,
Chungbuk, 28644 South Korea  P. Bangale  Department of Physics and Astronomy
and the Bartol Research Institute, University of Delaware, Newark, DE 19716,
USA  P. Batista  DESY, Platanenallee 6, 15738 Zeuthen, Germany  W. Benbow
Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA
J. H. Buckley  Department of Physics, Washington University, St. Louis, MO
63130, USA  M. Capasso  Department of Physics and Astronomy, Barnard College,
Columbia University, NY 10027, USA  J. L. Christiansen  Physics Department,
California Polytechnic State University, San Luis Obispo, CA 94307, USA  A. J.
Chromey  Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA
02138, USA  A. Falcone  Department of Astronomy and Astrophysics, 525 Davey
Lab, Pennsylvania State University, University Park, PA 16802, USA  Q. Feng
Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA
J. P. Finley  Department of Physics and Astronomy, Purdue University, West
Lafayette, IN 47907, USA  G. M Foote  Department of Physics and Astronomy and
the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA
G. Gallagher  Department of Physics and Astronomy, Ball State University,
Muncie, IN 47306, USA  W. F Hanlon  Center for Astrophysics $|$ Harvard &
Smithsonian, Cambridge, MA 02138, USA  D. Hanna  Physics Department, McGill
University, Montreal, QC H3A 2T8, Canada  O. Hervet  Santa Cruz Institute for
Particle Physics and Department of Physics, University of California, Santa
Cruz, CA 95064, USA  J. Holder  Department of Physics and Astronomy and the
Bartol Research Institute, University of Delaware, Newark, DE 19716, USA  T.
B. Humensky  Department of Physics, University of Maryland, College Park, MD
20742-4111, USA  NASA GSFC, Greenbelt, MD 20771, USA  W. Jin  Department of
Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA  P.
Kaaret  Department of Physics and Astronomy, University of Iowa, Van Allen
Hall, Iowa City, IA 52242, USA  M. Kertzman  Department of Physics and
Astronomy, DePauw University, Greencastle, IN 46135-0037, USA  D. Kieda
Department of Physics and Astronomy, University of Utah, Salt Lake City, UT
84112, USA  T. K. Kleiner  DESY, Platanenallee 6, 15738 Zeuthen, Germany  N.
Korzoun  Department of Physics and Astronomy and the Bartol Research
Institute, University of Delaware, Newark, DE 19716, USA  F. Krennrich
Department of Physics and Astronomy, Iowa State University, Ames, IA 50011,
USA  S. Kumar  Department of Physics, University of Maryland, College Park, MD
20742-4111, USA  M. J. Lang  School of Natural Sciences, University of Galway,
University Road, Galway, H91 TK33, Ireland  G. Maier  DESY, Platanenallee 6,
15738 Zeuthen, Germany  C. E McGrath  School of Physics, University College
Dublin, Belfield, Dublin 4, Ireland  C. L. Mooney  Department of Physics and
Astronomy and the Bartol Research Institute, University of Delaware, Newark,
DE 19716, USA  P. Moriarty  School of Natural Sciences, University of Galway,
University Road, Galway, H91 TK33, Ireland  R. Mukherjee  Department of
Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA  S.
O’Brien  Physics Department, McGill University, Montreal, QC H3A 2T8, Canada
R. A. Ong  Department of Physics and Astronomy, University of California, Los
Angeles, CA 90095, USA  N. Park  Department of Physics, Engineering Physics &
Astronomy, Queen’s University, Kingston, ON K7L 3N6, Canada  S. R. Patel
DESY, Platanenallee 6, 15738 Zeuthen, Germany  K. Pfrang  DESY, Platanenallee
6, 15738 Zeuthen, Germany  M. Pohl  Institute of Physics and Astronomy,
University of Potsdam, 14476 Potsdam-Golm, Germany and DESY, Platanenallee 6,
15738 Zeuthen, Germany  E. Pueschel  DESY, Platanenallee 6, 15738 Zeuthen,
Germany  J. Quinn  School of Physics, University College Dublin, Belfield,
Dublin 4, Ireland  K. Ragan  Physics Department, McGill University, Montreal,
QC H3A 2T8, Canada  P. T. Reynolds  Department of Physical Sciences, Munster
Technological University, Bishopstown, Cork, T12 P928, Ireland  E. Roache
Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA
I. Sadeh  DESY, Platanenallee 6, 15738 Zeuthen, Germany  L. Saha  Center for
Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA  G. H.
Sembroski  Department of Physics and Astronomy, Purdue University, West
Lafayette, IN 47907, USA  D. Tak  DESY, Platanenallee 6, 15738 Zeuthen,
Germany  J. V. Tucci  Department of Physics, Indiana University-Purdue
University Indianapolis, Indianapolis, IN 46202, USA  A. Weinstein  Department
of Physics and Astronomy, Iowa State University, Ames, IA 50011, USA  D. A.
Williams  Santa Cruz Institute for Particle Physics and Department of Physics,
University of California, Santa Cruz, CA 95064, USA  J. Woo  Columbia
Astrophysics Laboratory, 550 West 120th Street, NY 10027, USA

###### Abstract

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-LAT 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\pm 0.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\sim 3$ $\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.

Boomerang, pulsar wind nebula, NuSTAR, leptonic, X-ray, TeV gamma-rays

††facilities: NuSTAR, Chandra, Fermi-LAT, VERITAS††software: HEAsoft (v6.28),
NuSTARDAS (v2.0.0), XSPEC (v12.11.1; Arnaud (1996)), CIAO (v4.13), Naima
(Zabalza, 2015)††thanks: I. Pope<EMAIL_ADDRESS>M. Capasso,
<EMAIL_ADDRESS>N. Park<EMAIL_ADDRESS>J.
Woo<EMAIL_ADDRESS>
##  1 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_{\gamma}\sim 100$ TeV
(Abeysekara et al., 2020; Cao et al., 2021). Eventually, after $\tau\sim 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.

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.

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\times
10^{37}$ erg s-1 and a characteristic age of $\sim$10 kyr (Halpern et al.,
2001a). A compact PWN with a $r\sim 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 $\sim 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 the bow structure (Kothes et al.,
2001, 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, 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 $\mu$G (Liang et al., 2022).

In the X-ray band, the pulsar, its compact nebula ($r\sim 30$′′) and diffuse
emission over $r\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}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\sim 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^{\degree}$. 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 12CO cloud
$J=1\rm{-}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 $\gamma$ (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.

Table 1: Boomerang distance estimates d | method | citation
---|---|---
0.8 kpc | H I radial velocities | Kothes et al. (2001)
3 kpc | column density ($N_{H}$) | Halpern et al. (2001b)
5-7.5 kpc | dispersion measurement | Yao et al. (2017); Abdo et al. (2009b)
12 kpc | H I radial velocities | Pineault & Joncas (2000)

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\times 10^{21}$ cm-2. Since the
column density through the entire Galaxy is only $8.4\times 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\pm 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^{\circ}$. 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 ${\dot{E}}$. Additionally,
the H I column density from X-ray observations of $6.3\times 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 $\sim 15$′) and
tail (scale $\sim 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.

##  2 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 region varied according to the point spread function (PSF)
of each telescope (§2.2 and 2.3). All errors are given to the $2\sigma$
confidence level unless explicitly stated otherwise.

###  2.1 X-ray observations

Chandra observed the Boomerang region with its CCD ACIS-I array on 2002 March
15 (Obs ID 2787, PI Halpern) for $\sim 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\sim 10$′′; Ng & Romani, 2004), (3) X-ray PWN ($r\sim 30$′′), and
(4) diffuse X-ray emission ($r\sim 100$′′). These X-ray features are resolved
by Chandra as shown in the Chandra image in the 0.5–8.0 keV 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).

####  2.1.1 NuSTAR timing analysis

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:  (top) Chandra 0.5-8 keV image of the Boomerang PWN, smoothed to
$\sigma=\sim 0.5^{\prime\prime}$ (1 pixel). The color scale is in units of
counts per pixel. A $r=1^{\prime\prime}$ 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.

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\times 10^{-6}}_{-3\times
10^{-6}}$ ms ($3\sigma$ 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 $\dot{P}$ 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 $\phi=0.17$ and 0.30, as well as $\phi=0.70$ and
0.85 using extractor for subsequent spectral and imaging analysis of the
nebular emission. The effect of any leftover pulsar component post phase
extraction is negligible, as we determine in Section 2.1.3.

####  2.1.2 NuSTAR and Chandra imaging analysis

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
$\sim$5′′ radial bin widths.

The Chandra observations detected the Boomerang PWN extending to the bounds of
its radio emission ($r\approx 100^{\prime\prime}$; 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 energy-
resolved 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 $\sigma=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_{\rm in}=5^{\prime\prime}$ and $r_{\rm out}=100^{\prime\prime}$ 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.

(a)  (b)

Figure 5:  NuSTAR 3–10 keV (a) and 10–20 keV (b) background subtracted flux
images of the Boomerang PWN [10-5 counts s-1 cm-2]. 1420 MHz radio contours
are shown in green and the position of PSR J2229$+$6114 is marked by the blue
cross.  Figure 6: NuSTAR X-ray radial profiles of the Boomerang PWN in the
soft (red) and hard (blue) bands, compared to the NuSTAR PSF. The vertical
lines indicate the SHERPA best fit 2D-Gaussian FWHM in each energy band. The 1
sigma errors are indicated by the shaded regions. The horizontal bars
designate the 4.75′′ radial bin widths.

The background subtracted source profiles, as shown in Figure 6, are extended
above the NuSTAR PSF up to $r\sim 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\pm 2$′′ for the soft band and $20\pm 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_{\rm in}=1^{\prime\prime}$ (to mask out the pulsar emission)
to $r_{\rm out}=100^{\prime\prime}$ (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 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.

(a)  (b)

Figure 7: NuSTAR-only spectra (a) and NuSTAR–Chandra joint spectra (b) fit to an absorbed power-law model.  Table 2: Boomerang PWN spectral fitting parameters. All errors are given to the 90% confidence level. Parameter | NuSTAR |  NuSTAR \+ Chandra | Chandra
---|---|---|---
Model | tbabs*const*pow | tbabs*const*pow | tbabs*pow
$N_{H}$ [$10^{22}$ cm-2] | 0.89 (frozen) | 0.89 (frozen) | 0.89${}^{+0.15}_{-0.14}$
$\Gamma_{X}$ | 1.52 $\pm$ 0.15 | 1.52 $\pm$ 0.06 | 1.52${}^{+0.13}_{-0.12}$
Flux normalizationaaFlux normalization at 1 keV [$10^{-4}$ photons cm-2 s-1 keV-1]. For the NuSTAR and joint spectral fit, the flux normalization corresponds to the FPMA and Chandra spectra, respectively. | 1.89${}^{+0.63}_{-0.49}$ | 1.95 $\pm$ 0.11 | 1.94${}^{+0.33}_{-0.28}$
Cross-normalization | 0.85${}^{+0.11}_{-0.10}$ | 0.92${}^{+0.09}_{-0.08}$ | —
$\chi^{2}_{\nu}$ (d.o.f) | 1.21 (79) | 0.95 (226) | 0.79 (141)
$F_{X}$ (0.5–10 keV)bbAbsorbed X-ray flux in $10^{-12}$ ergs s-1 cm-2. | 0.98 $\pm$ 0.06 | $1.10\pm 0.03$ | 1.15 $\pm$ 0.04
$F_{X}$ (10–20 keV)bbAbsorbed X-ray flux in $10^{-12}$ ergs s-1 cm-2. | 0.70 $\pm$ 0.04 | — | —

####  2.1.3 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 $\sigma$
over background counts, respectively. In order to determine the hydrogen
column density accurately, we fit the 0.5–8.0 keV 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_{\rm H}=6.2^{+1.0}_{-0.9}\times
10^{21}~{}\rm{cm}^{-2}$. We then repeated the Chandra spectral fit using the
Wilms abundance table (Wilms et al., 2000), obtaining a higher value of
$N_{\rm H}=8.9^{+1.5}_{-1.4}\times 10^{21}~{}\rm{cm}^{-2}$. Radio observations
of PSR J2229$+$6114 measured its DM to be $(204.97\pm 0.02)$ cm-3 pc (Abdo et
al., 2009b). Using a linear relation between $N_{\rm H}$ and DM as well as its
slope errors (He et al., 2013), we derived $N_{\rm H}=(4.3\rm{-}8.8)\times
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_{\rm 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_{\rm H}=8.9\times 10^{21}~{}\rm{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.0
keV 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
$\Gamma=1.52\pm 0.06$ with a reduced chi-squared value of
$\chi_{\nu}^{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^{\prime\prime}$ 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.

###  2.2 Fermi-LAT data selection and analysis

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.

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 selected events with time stamps between MET 239557418 (2008-08-04
15:43:37.000 UTC) and MET 644798253 (2021-06-07 22:37:28.000 UTC), in a
10°-wide region around the 4FGL catalog counterpart to PSR J2229+6114, i.e.
4FGL J2229.0+6114. The analysis was conducted in the energy range 3 GeV - 2
TeV.

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).

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 = $\mathrm{-2ln(L_{max,0}/L_{max,1})}$ where
$\mathrm{L_{max,0}}$ and $\mathrm{L_{max,1}}$ are the maximum likelihoods for
a model excluding and including the source, respectively
111https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/Likelihood_overview.html),
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
$\mathrm{\sim 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
$\mathrm{\sim 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
$\mathrm{2.91\times 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.

###  2.3 VERITAS observation and analysis

Table 3: VERITAS upper limits with 99% confidence level Ethreshold (TeV) | index | Upper limits (10-12 s-1 cm-2)
---|---|---
0.38 | 2.0 | 1.10
0.38 | 2.5 | 1.23
0.35 | 3.0 | 1.52
0.79 | 2.0 | 0.40
0.72 | 2.5 | 0.48
0.72 | 3.0 | 0.49

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$\sigma$
significance within 25 hours. VERITAS has an angular resolution of $\sim$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.

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.

##  3 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
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., 2015, 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.

###  3.1 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_{\rm min}$ and $E_{\rm max}$. While a
hard cutoff at $E_{e}=E_{\rm 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.

Figure 10: NAIMA simple power-law SED fits with B-field set to 2.6 mG (left) and $5~{}\mu$G (right). The complete set of parameters for each fit, as well as for other B-field scenarios, is shown in Table 4.  Table 4: Fit parameters from various leptonic NAIMA model cases. Distance [kpc] |  $B$ [$\mu$G] | $p$ |  $E_{\rm min}$ [GeV] |  $E_{\rm max}$ [TeV] | Electron energy [erg] | Magnetic field energy [erg]aa$\eta_{E}$ and $\eta_{B}$ are calculated by dividing the electron and magnetic field energy by their sum, respectively. | $\eta_{E}$ | $\eta_{B}$
---|---|---|---|---|---|---|---|---
0.8 | 2600 | 2.5 | 0.32 | 22 | $2.3\times\ 10^{41}$ | $2.0\times\ 10^{48}$ | $\sim 0$ | $\sim 1$
0.8 | 100 | 2.5 | 1.64 | 120 | $2.9\times\ 10^{43}$ | $2.9\times\ 10^{45}$ | 0.01 | 0.99
0.8 | 10 | 2.5 | 5.18 | 400 | $1.0\times\ 10^{45}$ | $2.9\times\ 10^{43}$ | 0.97 | 0.03
0.8 | 5 | 2.5 | 7.33 | 500 | $2.8\times\ 10^{45}$ | $7.3\times\ 10^{42}$ | 0.998 | 0.002
7.5 | 2600 | 2.5 | 0.32 | 22 | $2.0\times\ 10^{43}$ | $1.6\times\ 10^{51}$ | $\sim 0$ | $\sim 1$
7.5 | 100 | 2.5 | 1.64 | 120 | $2.5\times\ 10^{45}$ | $2.9\times\ 10^{48}$ | 0.001 | 0.999
7.5 | 10 | 2.5 | 5.18 | 400 | $9.1\times\ 10^{46}$ | $2.9\times\ 10^{46}$ | 0.76 | 0.24
7.5 | 5 | 2.5 | 7.33 | 500 | $2.5\times\ 10^{47}$ | $7.3\times\ 10^{45}$ | 0.97 | 0.03

We first focused on reproducing the radio spectral break at $\sim 5$ GHz and
X-ray spectra. We found that the radio data are adequately fit by various sets
of $B$ and $E_{\rm min}$ values, as listed in Table 4. We have shown that the
radio spectral break can be caused by the selected $E_{\rm 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_{\rm 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=\Gamma_{x}-1$. We therefore infer from $p=2.5$ that
$\Gamma_{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 $\mu$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~{}\mu$G.

NAIMA outputs the total energy of the electron population, $U_{E}$, integrated
between $E_{\rm min}$ and $E_{\rm max}$. Assuming the B-field is uniform
throughout the interior of the PWN ($r<100^{\prime\prime}$), we calculated the
total magnetic field energy as $U_{B}=\frac{B^{2}}{8\pi}\times V_{\rm PWN}$,
where $V_{\rm PWN}=(4/3)\pi R_{\rm PWN}^{3}$ 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 $\sim$1$\rm{-}10$ $\mu$G to B
$\;\buildrel>\over{\sim}\;$ 100 $\mu$G requires decreasing the electron
density dramatically. As a result, a fraction of the total energy allocated to
PWN B-field ($\eta_{B}$) is $\sim 1$ for B $\;\buildrel>\over{\sim}\;$
100$\mu$G. Such a high magnetization parameter is unusual if compared to other
PWN systems – e.g., Martin et al. (2014) found $\eta_{B}=7\times
10^{-4}\rm{-}0.02$ from six PWNe including the Crab nebula (which possesses
the highest $\eta_{B}$ value). If we attain an $\eta_{B}$ value comparable to
those deduced from the other PWNe, it points to a lower B-field range of B
$\;\buildrel<\over{\sim}\;$ 10 $\mu$G. While these arguments as well as the
suggested B-field range (B$\sim$5$\rm{-}10$ $\mu$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.

###  3.2 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_{\rm PWN}$ increased linearly with time over the short lifetime
of the new PWN and matched its current radius (at $t=t_{\rm 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 $\dot{E}=2.2\times 10^{37}$ ergs
s-1. A fraction of the spin-down power ($\eta_{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-\eta_{g}-\eta_{B})\dot{E}$ where $\eta_{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_{\rm PWN}=R_{\rm
PWN}/t_{\rm age}$). We then calculated a radiation SED from the particle
distribution at $t=t_{\rm 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.

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.

We also attempted to constrain our model based on the lack of X-ray
variability over the last $\sim 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_{\rm syn}(t)$ evolution and found solutions without
significant X-ray variability over the last $\sim 20$ years.

Figure 12: Same plot as Figure 11 for $d=7.5$ kpc.  Table 5: Model parameters for the case presented in §3.2 Parameter |  $d=0.8$ kpc |  $d=7.5$ kpc
---|---|---
Age [yr] | 30 | 1000
$p$ | 2.4 | 2.4
$E_{\rm min}$ [GeV] | 7.9 | 10.0
$E_{\rm max}$ [PeV] | 1.0 | 1.0
$\eta_{g}$aaThe gamma-ray efficiency of PSR J2229+6114 from Abdo et al. (2009b). $f_{\Omega}$ (Eq. (3) in the cited work) was assumed to be 1. | $1\times 10^{-3}$ | $7\times 10^{-2}$
$\eta_{B}$ | $7\times 10^{-5}$ | $1\times 10^{-2}$
$B$ [$\mu$G] | 2.1 | 4.5
Electron energy [erg] | $8.9\times 10^{45}$ | $2.9\times 10^{47}$
Magnetic field energy [erg] | $1.3\times 10^{42}$ | $6.3\times 10^{45}$
Expansion velocity [km/s] | $1.3\times 10^{4}$ | $3.9\times 10^{3}$

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_{\rm age}\sim 30$
yr for $d=0.8$ kpc. If $t_{\rm age}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}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 $\tau\sim 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\sim 2~{}\mu$G (with a very low
magnetization of $\eta_{B}=7\times 10^{-5}$). Otherwise, the model over-
predicts synchrotron radiation fluxes in the radio and X-ray band. We also
found that $E_{\rm min}$ is well constrained by the radio spectral break and
the overall flux normalization of electrons (which depends sensitively on
$E_{\rm min}$ for a given $\dot{E}$ and lifetime). The expansion velocity
should be high ($V_{\rm PWN}=1.9\times 10^{4}$ km/s) in order to reach the
observed PWN radius within $\sim 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_{\rm esc}<0.05$
yr), were we able to fit the SED data with higher $\eta_{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\times 10^{47}$ ergs. While the injected electron spectral
parameters such as $p$, $E_{\rm min}$ and $E_{\rm max}$ are similar between
the two distance cases (Table 5), $\eta_{B}$ is higher at $1\times 10^{-2}$
while the current B-field is $4.5~{}\mu$G. The expansion velocity is
$3.9\times 10^{3}$ km/s, which is more in line with the observed PWN
velocities ($v_{\rm PWN}\sim 1,000\rm{-}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\sim 5~{}\mu$G in order to
match the low radiation efficiency ($L_{\rm syn}/\dot{E}$). 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_{\rm 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_{\rm
min}$. In both cases, $E_{\rm min}\sim 10$ GeV fit the radio SED data well,
and thus we attribute the break to the PWN’s intrinsic particle injection
distribution.

###  3.3 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_{\rm 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 $\dot{E}$ ($=2.2\times
10^{37}$ ergs s-1) and characteristic age $t_{\rm ch}$ (= 11 kyr), we first
derive the system’s true age ($t_{\rm age}$) and initial spin-down luminosity
$\dot{E}_{0}$

$\displaystyle t_{\rm age}$ $\displaystyle=$ $\displaystyle\frac{2t_{\rm
ch}}{p-1}-\tau_{\rm sd}$ (1) $\displaystyle\dot{E}_{0}$ $\displaystyle=$
$\displaystyle\dot{E}\left(1+\frac{t_{\rm age}}{\tau_{\rm
sd}}\right)^{\frac{p+1}{p-1}},$ (2)

where $p$ and $\tau_{\rm 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 $\dot{E}(t)$ into
$(1-\eta_{B})\dot{E}(t)$ and $\eta_{B}\dot{E}(t)$, respectively. The allocated
electron energy is distributed between $E_{e}=E_{\rm min}$ and $E_{\rm 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
$\rho_{\rm ej}(r)\propto 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_{\rm
radio}/\dot{E}=3\times 10^{-8}$ and $L_{\rm X}/\dot{E}=7\times 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 $\eta_{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_{\rm 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 cm3) 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.

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 Table 4.  Table 6:
Model parameters for the four cases presented in §3.3

Model parameter | Case A | Case B | Case B’ | Case C
---|---|---|---|---
Source distance [kpc] | 0.8 | 7.5 | 7.5 | 7.5
SN explosion energy [ergs] | $3.0\times 10^{50}$ | $3.6\times 10^{51}$ | $4.0\times 10^{51}$ | $1.9\times 10^{51}$
SN ejecta mass [$M_{\odot}$] | 1.6 | 4.2 | 4.9 | 2.0
ISM density [cm-3] | 160 | 0.3 | 0.3 | 0.9
Magnetic field strength [$\mu$G] | 1.5 | 3.7 | 2.9 | 2.5
PWN radius [pc] | 0.35 | 3.8 | 4.6 | 2.9
Pulsar braking index | 5.6 | 2.9 | 2.8 | 3.1
Pulsar spin-down timescale [kyr] | 3.9 | 8.9 | 9.1 | 7.7
True age [kyr] | 0.64 | 1.8 | 2.2 | 2.1
Initial spin-down power [ergs/s] | $2.7\times 10^{37}$ | $2.2\times 10^{37}$ | $3.4\times 10^{37}$ | $3.5\times 10^{37}$
Wind magnetization ($\eta_{B}$) | $8\times 10^{-8}$ | 0.006 | 0.005 | 0.0007
$E_{\rm min}$ [GeV] | 0.4 | 22.3 | 20.0 | 8.7
$E_{\rm max}$ [PeV] | 0.2 | 4.0 | 3.5 | 1.2
$E_{\rm break}$ [TeV] | 0.14 | 335.77 | 300.00 | 0.01
p1 | 1.5 | 2.6 | 2.6 | 1.6
p2 | 3.1 | 1.3 | 1.3 | 2.3

††footnotetext: 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_{\rm break}$, respectively.

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

(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\sim 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 Table 6 and Figure 13 (upper right 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 $\mu$G while the
PWN expanded to $r=3.8$ pc. Note that the power-law spectral index is softer
below than above $E_{\rm 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_{\rm
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_{\rm 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~{}\mu$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 $\sim 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\sim 2.3$ (i.e., standard Fermi).

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\sim 2\rm{-}4~{}\mu$G.
Alternatively, we found that an extremely magnetized PWN with $\eta_{B}\sim 1$
can fit the SED data well. More specifically, when we adopt $\eta_{B}=0.99$,
only 1% of the pulsar’s rotational energy is allocated to particle injection
and the current B-field is $\sim 100~{}\mu$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 ($\eta_{B}=7\times
10^{-4}\rm{-}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’.

##  4 Discussion

We discuss constraining the properties of the Boomerang PWN and its recent
evolution based on the X-ray and multi-wavelength 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.

###  4.1 Constraining PWN magnetic field

As shown in §3, our SED study strongly suggests that the current B-field
should be $B\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}3\mu$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 $\tau_{\rm syn}=1.2B_{\rm mG}^{-3/2}E_{\rm
keV}^{-1/2}$ 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\sim 0.7$ pc and
$\sim 0.4$ pc, respectively. As expected, these radii are significantly
smaller than the PWN radius assumed for the SED modelling ($\sim$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_{\gamma}=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_{\rm PWN}\sim 1,000$ km s-1 or
$3.3\times 10^{-3}c$ (Reynolds et al., 2018). The Crab nebula is known to
expand at a velocity of $V_{\rm PWN}\sim 1,500$ km s-1 (Bietenholz et al.,
1991) and its radius increases almost linearly with time – $R(t)\propto
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\sim
7\rm{-}10~{}\mu$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_{\rm PWN}\sim 2,400$ km/s and $\sim 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\sim
6\rm{-}14~{}\mu$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 $\mu$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.

###  4.2 PWN evolution

The dynamical SED model fitting §3.3 suggests that the Boomerang PWN is
currently re-expanding to $r\sim 3$ pc after being crushed by a SNR reverse
shock $\sim 1000$ yr ago, much shorter than the pulsar’s $\sim 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 $\sim 1000$ yr. After the PWN compression amplified B-field to $B\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10$ mG, its B-field decreased, as
of present day, to $B\sim 3~{}\mu$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\sim 0.4\rm{-}0.7$ pc)
should be powered by a population of young electrons injected over the last
$\sim 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.

Another key parameter is the proper motion of the pulsar. Although there is
yet no direct measurement of the proper motion, the radio and X-ray morphology
of the PWN suggest that the pulsar is moving in the north–west direction. It
was proposed that the Boomerang-like radio morphology was caused by the PWN
colliding into the high density ISM region in the north–west boundary of SNR
G106.3$+$2.7. In this scenario, the head is more extended because the relic
electrons from the PWN diffuse into a lower density region. In addition, the
X-ray torus–jet morphology data was fit by a PWN emission model (Ng & Romani,
2004) also supporting the proper motion in the north–west direction. Although
a future Chandra observation over $\sim 20$ year baseline may be able to
detect the proper motion, it may be difficult if the source distance is indeed
$\sim 8$ kpc or greater as suggested by our SED and morphology studies
presented in this paper.

##  5 Conclusions

We summarize our findings from our multi-wavelength investigation of the
Boomerang PWN and its surrounding region.

  * (1)

We detected a 51.67 ms pulsation with 3–20 keV NuSTAR data and it allowed us
to separate the on-pulse component from the PWN emission. With NuSTAR we
detected hard X-ray emission from the Boomerang PWN up to 20 keV and found
that its size decreases from $r=33\pm 2^{\prime\prime}$ (3–10 keV) to $20\pm
2^{\prime\prime}$ (10–20 keV).

  * (2)

Our analysis of the 2002 Chandra data, after excising the pulsar emission,
yielded $N_{\rm H}=8.9^{+1.5}_{-1.4}\times 10^{21}$ cm-2 and the PWN photon
index of $\Gamma=1.52^{+0.13}_{-0.12}$. The hydrogen column density is
consistent with the value derived from the pulsar’s DM measurement of
$(4.3\rm{-}8.8)\times 10^{21}$ cm-2. Joint Chandra and NuSTAR spectral
analysis of the PWN measured no X-ray flux variability and the 0.5–20 keV
spectra are consistent with a single power-law model with $\Gamma=1.52\pm
0.06$.

  * (3)

Our analysis of the most updated Fermi-LAT and VERITAS data yielded no
detection of the Boomerang PWN. We set an upper limit on the gamma-ray flux
above 50 GeV to $F\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}10^{-12}$
ergs s-1 cm-2.

  * (4)

Among the previously suggested source distances ranging from 0.8 kpc to 7.5
kpc (Kothes et al., 2001; Abdo et al., 2009b; He et al., 2013), we found that
$d\sim 8$ kpc provides the most plausible solution to fitting the SED data.
The widely-used source distance of 0.8 kpc leads to too low radiation
efficiencies that cannot be produced by any of our SED models with reasonable
parameters. We note, however, that all SED models used in this paper simplify
the complex nature of Boomerang. While the offset between the radio and X-ray
peak emission may suggest a multi-zone system, modelling was done under the
assumption that Boomerang is a single-zone system. Implementing an multi-zone
model may boost our understanding of the emission offset between the radio,
X-ray, and TeV peak emission. The dynamical SED model used in this paper also
assumes that Boomerang is a spherically symmetric PWN, an obvious
simplification of the actual morphology. As can be seen in Figure 13, the
dynamical SED fits predict strong emission at slightly higher radio
frequencies and/or X-ray energies than what has been observed so far. Future
analysis of observations in these energy bands could further test the validity
of these fits.

  * (5)

The radiation efficiencies (i.e. $F(4\pi R^{2})/\dot{E}$ ratios) of Boomerang
are $3\times 10^{-6}(\frac{R}{8~{}kpc})^{2}$, $6\times
10^{-4}(\frac{R}{8~{}kpc})^{2}$ and $4\times 10^{-4}(\frac{R}{8~{}kpc})^{2}$
in the radio, X-ray and TeV bands, respectively. Assuming a distance of 8 kpc,
these values are now more consistent with those of other PWNe with similar
spin-down powers (Kargaltsev et al., 2013).

  * (6)

Our SED modelling requires that the current PWN B-field should be low ($B\sim
2\rm{-}4~{}\mu$G). Our X-ray PWN size measurements suggest a slightly higher
value ($B\sim 6\rm{-}14~{}\mu$G). We ruled out the high B-field value ($B=2.6$
mG) suggested by the radio observations (Kothes et al., 2006), as well as the
B-field value (B = 140 $\mu$G) suggested by Liang et al. (2022).

  * (7)

The Boomerang PWN should currently be in its re-expansion phase after having
been crushed by the SNR reverse shock. The relic electrons injected earlier
before the SNR crush should have produced synchrotron and ICS radiation
elsewhere, e.g. in the head region. We attribute the radio break at 4–5 GHz to
the minimum energy of the injected electrons. While the origin of this break
in the electron spectrum is still unknown, we are confident that the radio
break is not caused by synchrotron cooling as was hypothesized by Kothes et
al. (2006).

  * (8)

The origin of the head–tail morphology could be related to the PWN propagating
in an inhomogeneous density region and interacting with the SNR, as suggested
by HD simulations. Following Kothes et al. (2006), we also suspect that the
PWN is currently propagating in the northwest direction and expanding in high-
density ISM. In order to understand the origin of VHE and UHE gamma rays and
their connection to the SNR-PWN system, it would be important to resolve the
spatial distribution in the TeV emission in the head and tail regions as well
as determine the distance and proper motion of the PSR J2229$+$6114.

Support for this work was provided by NASA through NuSTAR Cycle 6 Guest
Observer Program grant NNH19ZDA001N. VERITAS is supported by grants from the
U.S. Department of Energy Office of Science, the U.S. National Science
Foundation and the Smithsonian Institution, by NSERC in Canada, and by the
Helmholtz Association in Germany. This research used resources provided by the
Open Science Grid, which is supported by the National Science Foundation and
the U.S. Department of Energy’s Office of Science, and resources of the
National Energy Research Scientific Computing Center (NERSC), a U.S.
Department of Energy Office of Science User Facility operated under Contract
No. DE-AC02-05CH11231. We acknowledge the excellent work of the technical
support staff at the Fred Lawrence Whipple Observatory and at the
collaborating institutions in the construction and operation of the
instrument. , GAMERA (Hahn, 2016)

## References

  * Abdelmaguid et al. (2022) Abdelmaguid, M., Gelfand, J., Straal, S., & Gotthelf, E. 2022, X-ray Spectral Analysis and Modeling of the TeV source HESS J1640-465, Zenodo, doi: 10.5281/zenodo.6674275
  * Abdo et al. (2009a) Abdo, A. A., Allen, B. T., Aune, T., et al. 2009a, The Astrophysical Journal, 700, L127, doi: 10.1088/0004-637x/700/2/l127
  * Abdo et al. (2009b) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009b, The Astrophysical Journal, 706, 1331, doi: 10.1088/0004-637x/706/2/1331
  * Abdo et al. (2013) Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, The Astrophysical Journal Supplement Series, 208, 17, doi: 10.1088/0067-0049/208/2/17
  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, Science, 358, 911, doi: 10.1126/science.aan4880
  * Abeysekara et al. (2020) —. 2020, PRL, 124, 021102, doi: 10.1103/PhysRevLett.124.021102
  * Acciari et al. (2009) Acciari, V. A., Aliu, E., Arlen, T., Aune, T., & Bautista, M. 2009, ApJ, 703, L6, doi: 10.1088/0004-637X/703/1/L6
  * Acero et al. (2016) Acero, F., Ackermann, M., Ajello, M., et al. 2016, The Astrophysical Journal Supplement Series, 223, 26, doi: 10.3847/0067-0049/223/2/26
  * Albert et al. (2020) Albert, A., Alfaro, R., Alvarez, C., et al. 2020, The Astrophysical Journal, 896, L29, doi: 10.3847/2041-8213/ab96cc
  * Amenomori et al. (2021) Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2021, Nature Astronomy, 5, 460, doi: 10.1038/s41550-020-01294-9
  * Anders & Grevesse (1989) Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197, doi: 10.1016/0016-7037(89)90286-X
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Arons (2012) Arons, J. 2012, Space Science Reviews, 173, 341, doi: 10.1007/s11214-012-9885-1
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Bachetti et al. (2022) Bachetti, M., Huppenkothen, D., Khan, U., et al. 2022, StingraySoftware/stingray: Version 1.0, v1.0, Zenodo, doi: 10.5281/zenodo.6394742
  * Bao & Chen (2021) Bao, Y., & Chen, Y. 2021, ApJ, 919, 32, doi: 10.3847/1538-4357/ac1581
  * Bietenholz et al. (1991) Bietenholz, M. F., Kronberg, P. P., Hogg, D. E., & Wilson, A. S. 1991, ApJ, 373, L59, doi: 10.1086/186051
  * Bietenholz & Nugent (2015) Bietenholz, M. F., & Nugent, R. L. 2015, MNRAS, 454, 2416, doi: 10.1093/mnras/stv2112
  * Breuhaus et al. (2022) Breuhaus, M., Hinton, J. A., Joshi, V., Reville, B., & Schoorlemmer, H. 2022, arXiv e-prints, arXiv:2201.03984.  https://arxiv.org/abs/2201.03984
  * Burgess et al. (2022) Burgess, D. A., Mori, K., Gelfand, J. D., et al. 2022, arXiv e-prints, arXiv:2204.03185.  https://arxiv.org/abs/2204.03185
  * Cao et al. (2021) Cao, Z., Aharonian, F. A., An, Q., et al. 2021, Nature, 594, 33, doi: 10.1038/s41586-021-03498-z
  * Cholis et al. (2018) Cholis, I., Karwal, T., & Kamionkowski, M. 2018, PRD, 98, 063008, doi: 10.1103/PhysRevD.98.063008
  * Coerver et al. (2019) Coerver, A., Wilcox, P., Zhang, H., et al. 2019, ApJ, 878, 126, doi: 10.3847/1538-4357/ab21d0
  * Cordes & Lazio (2002) Cordes, J., & Lazio, T. 2002
  * Cox et al. (1986) Cox, P., Kruegel, E., & Mezger, P. G. 1986, A&A, 155, 380
  * Fang et al. (2022) Fang, K., Kerr, M., Blandford, R., Fleischhack, H., & Charles, E. 2022, Phys. Rev. Lett., 129, 071101, doi: 10.1103/PhysRevLett.129.071101
  * Ferrand & Safi-Harb (2012) Ferrand, G., & Safi-Harb, S. 2012, Advances in Space Research, 49, 1313, doi: 10.1016/j.asr.2012.02.004
  * Fujita et al. (2021) Fujita, Y., Bamba, A., Nobukawa, K. K., & Matsumoto, H. 2021, ApJ, 912, 133, doi: 10.3847/1538-4357/abf14a
  * Gaensler & Slane (2006) Gaensler, B. M., & Slane, P. O. 2006, Annual Review of Astronomy and Astrophysics, 44, 17, doi: 10.1146/annurev.astro.44.051905.092528
  * Ge et al. (2021) Ge, C., Liu, R.-Y., Niu, S., Chen, Y., & Wang, X.-Y. 2021, The Innovation, 2, 100118, doi: 10.1016/j.xinn.2021.100118
  * Gelfand (2017) Gelfand, J. D. 2017, Radiative Models of Pulsar Wind Nebulae, ed. D. F. Torres, Vol. 446, 161, doi: 10.1007/978-3-319-63031-1_8
  * Gelfand et al. (2009) Gelfand, J. D., Slane, P. O., & Zhang, W. 2009, ApJ, 703, 2051, doi: 10.1088/0004-637X/703/2/2051
  * Giacinti et al. (2020) Giacinti, G., Mitchell, A. M. W., López-Coto, R., et al. 2020, A&A, 636, A113, doi: 10.1051/0004-6361/201936505
  * Hahn (2016) Hahn, J. 2016, in International Cosmic Ray Conference, Vol. 917, 34th International Cosmic Ray Conference (ICRC2015)
  * Halpern et al. (2001a) Halpern, J. P., Camilo, F., Gotthelf, E. V., et al. 2001a, The Astrophysical Journal, 552, L125, doi: 10.1086/320347
  * Halpern et al. (2002) Halpern, J. P., Gotthelf, E. V., Camilo, F., Collins, B., & Helfand, D. J. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 271, Neutron Stars in Supernova Remnants, ed. P. O. Slane & B. M. Gaensler, 199
  * Halpern et al. (2001b) Halpern, J. P., Gotthelf, E. V., Leighly, K. M., & Helfand, D. J. 2001b, The Astrophysical Journal, 547, 323, doi: 10.1086/318361
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103, doi: 10.1088/0004-637X/770/2/103
  * Hartman et al. (1999) Hartman, R. C., Bertsch, D. L., Bloom, S. D., et al. 1999, ApJS, 123, 79
  * Hattori et al. (2020) Hattori, S., Straal, S. M., Zhang, E., et al. 2020, ApJ, 904, 32, doi: 10.3847/1538-4357/abba32
  * He et al. (2013) He, C., Ng, C. Y., & Kaspi, V. M. 2013, ApJ, 768, 64, doi: 10.1088/0004-637X/768/1/64
  * Joncas & Higgs (1990) Joncas, G., & Higgs, L. A. 1990, Astron. Astrophysics. Suppl., 82, 113
  * Kargaltsev et al. (2013) Kargaltsev, O., Rangelov, B., & Pavlov, G. G. 2013, arXiv e-prints, arXiv:1305.2552.  https://arxiv.org/abs/1305.2552
  * Kolb et al. (2017) Kolb, C., Blondin, J., Slane, P., & Temim, T. 2017, ApJ, 844, 1, doi: 10.3847/1538-4357/aa75ce
  * Kothes et al. (2004) Kothes, R., Landecker, T. L., & Wolleben, M. 2004, The Astrophysical Journal, 607, 855, doi: 10.1086/383610
  * Kothes et al. (2006) Kothes, R., Reich, W., & Uyanıker, B. 2006, The Astrophysical Journal, 638, 225, doi: 10.1086/498666
  * Kothes et al. (2001) Kothes, R., Uyaniker, B., & Pineault, S. 2001, The Astrophysical Journal, 560, 236, doi: 10.1086/322511
  * Kuiper & Hermsen (2015) Kuiper, L., & Hermsen, W. 2015, Monthly Notices of the Royal Astronomical Society, 449, 3827, doi: 10.1093/mnras/stv426
  * Lang et al. (2010) Lang, C. C., Wang, Q. D., Lu, F., & Clubb, K. I. 2010, ApJ, 709, 1125, doi: 10.1088/0004-637X/709/2/1125
  * Liang et al. (2022) Liang, X.-H., Li, C.-M., Wu, Q.-Z., Pan, J.-S., & Liu, R.-Y. 2022, A PeVatron Candidate: Modelling the Boomerang Nebula in X-ray Band, arXiv, doi: 10.48550/ARXIV.2209.03809
  * Liu et al. (2022) Liu, Q.-C., Zhou, P., & Chen, Y. 2022, The Astrophysical Journal, 926, 124, doi: 10.3847/1538-4357/ac3fb3
  * López-Coto et al. (2022) López-Coto, R., de Oña Wilhelmi, E., Aharonian, F., Amato, E., & Hinton, J. 2022, arXiv e-prints, arXiv:2202.06899.  https://arxiv.org/abs/2202.06899
  * MAGIC Collaboration et al. (2023) MAGIC Collaboration, Abe, H., Abe, S., et al. 2023, A&A, 671, A12, doi: 10.1051/0004-6361/202244931
  * Martin et al. (2014) Martin, J., Torres, D. F., Cillis, A., & de Oña Wilhelmi, E. 2014, MNRAS, 443, 138, doi: 10.1093/mnras/stu1125
  * Ng & Romani (2004) Ng, C. Y., & Romani, R. W. 2004, ApJ, 601, 479, doi: 10.1086/380486
  * Nynka et al. (2014) Nynka, M., Hailey, C. J., Reynolds, S. P., et al. 2014, ApJ, 789, 72, doi: 10.1088/0004-637X/789/1/72
  * Ohira et al. (2018) Ohira, Y., Kisaka, S., & Yamazaki, R. 2018, MNRAS, 478, 926, doi: 10.1093/mnras/sty1159
  * Oka (2021) Oka, T. 2021, 37th ICRC (ICRC2021), 37, 796
  * Park & VERITAS Collaboration (2015) Park, N., & VERITAS Collaboration. 2015, in International Cosmic Ray Conference, Vol. 34, 34th International Cosmic Ray Conference (ICRC2015), 771\.  https://arxiv.org/abs/1508.07070
  * Perkins et al. (2009) Perkins, J. S., Maier, G., & The VERITAS Collaboration. 2009, arXiv e-prints, arXiv:0912.3841.  https://arxiv.org/abs/0912.3841
  * Pineault & Joncas (2000) Pineault, S., & Joncas, G. 2000, The Astronomical Journal, 120, 3218, doi: 10.1086/316863
  * Reynolds et al. (2018) Reynolds, S. P., Borkowski, K. J., & Gwynne, P. H. 2018, ApJ, 856, 133, doi: 10.3847/1538-4357/aab3d3
  * Straal et al. (2023) Straal, S. M., Gelfand, J. D., & Eagle, J. L. 2023, The Astrophysical Journal, 942, 103, doi: 10.3847/1538-4357/aca1a9
  * Strong et al. (2000) Strong, A. W., Moskalenko, I. V., & Reimer, O. 2000, The Astrophysical Journal, 537, 763, doi: 10.1086/309038
  * Taylor & Cordes (1993) Taylor, J. H., & Cordes, J. M. 1993, ApJ, 411, 674, doi: 10.1086/172870
  * Temim et al. (2013) Temim, T., Slane, P., Castro, D., et al. 2013, ApJ, 768, 61, doi: 10.1088/0004-637X/768/1/61
  * Temim et al. (2015) Temim, T., Slane, P., Kolb, C., et al. 2015, ApJ, 808, 100, doi: 10.1088/0004-637X/808/1/100
  * Torres et al. (2013) Torres, D. F., Martín, J., de Oña Wilhelmi, E., & Cillis, A. 2013, MNRAS, 436, 3112, doi: 10.1093/mnras/stt1793
  * Wik et al. (2014) Wik, D. R., Hornstrup, A., Molendi, S., et al. 2014, ApJ, 792, 48, doi: 10.1088/0004-637X/792/1/48
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016
  * Wood et al. (2017) Wood, M., Caputo, R., Charles, E., et al. 2017, in International Cosmic Ray Conference, Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 824.  https://arxiv.org/abs/1707.09551
  * Xin et al. (2019) Xin, Y., Zeng, H., Liu, S., Fan, Y., & Wei, D. 2019, ApJ, 885, 162, doi: 10.3847/1538-4357/ab48ee
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29, doi: 10.3847/1538-4357/835/1/29
  * Zabalza (2015) Zabalza, V. 2015, in ICRC, Vol. 34, 34th ICRC (ICRC2015), 922.  https://arxiv.org/abs/1509.03319