Discovery of a Pulsar Wind Nebula Candidate
Associated with the Galactic PeVatron 1LHAASO
J0343+5254u

Stephen Dikerby, Shuo Zhang, Tulun Ergin, Naomi Tsuji, Kaya Mori, Fabio

Acero, Samar Safi-Harb, Shunya Takekawa, Jooyun Woo

► To cite this version:

Stephen Dikerby, Shuo Zhang, Tulun Ergin, Naomi Tsuji, Kaya Mori, et al. Discovery of a Pulsar Wind Nebula Candidate Associated with the Galactic PeVatron 1LHAASO J0343+5254u. The Astrophysical Journal, 2025, 983 (1), pp.21. 10.3847/1538-4357/adb7e0. http://arxiv.org/abs/2005.07172

HAL Id: hal-04986172

https://hal.science/hal-04986172v1

Submitted on 22 May 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Discovery of a Pulsar Wind Nebula Candidate Associated with the Galactic PeVatron 1LHAASO J0343+5254u

Stephen DiKerby<fnref n="1" />, Shuo Zhang<fnref n="1" />, Tülün Ergin<fnref n="1" />, Naomi Tsuji<fnref n="2" />, Kaya Mori<fnref n="3" />, Fabio Acero<fnref n="4" />,<fnref n="5" />, Samar Safi-Harb<fnref n="6" />, Shunya Takekawa<fnref n="7" />, and Jooyun Woo<fnref n="3" />

## Abstract

The astronomical origins of the most energetic galactic cosmic rays and gamma rays are still uncertain. X-ray follow-up of candidate “PeVatrons”—systems producing cosmic rays with energies exceeding 1 PeV—can constrain their spatial origin, identify likely counterparts, and test particle emission models. Using ~120 ks of XMM-Newton observations, we report the discovery of a candidate pulsar wind nebula, a possible counterpart for the LHAASO PeVatron J0343+5254u. This extended source has a power-law X-ray spectrum with spectral index Γₓ = 1.9—softer at greater distance from the center—and asymmetric spatial extension out to ≈2'. We conduct leptonic modeling of the X-ray and gamma-ray radiation from this complex system, showing that a fully leptonic model with elevated IR photon fields can explain the multiwavelength emission from this source, similar to other very high-energy pulsar wind nebulas; excess gamma-ray emissivity not explained by a leptonic model may be due to hadronic interactions in nearby molecular cloud regions, which might also produce detectable astroparticle flux.

*Unified Astronomy Thesaurus concepts: X-ray astronomy (1810); Pulsar wind nebulae (2215); High energy astrophysics (739); Ultra-high-energy cosmic radiation (1733); Nebulae (1095); X-ray sources (1822)*

## 1. Introduction and Context

Galactic PeVatrons, astronomical systems producing cosmic rays (CRs) with energies reaching or exceeding 1 PeV = 10¹<fnref n="5" /> eV, are the most energetic astrophysical objects in our galaxy, hosting processes far exceeding the energies reached in terrestrial, artificial accelerators. The exact nature of and relationship of these objects with the galactic environment are still poorly understood; in particular, the “knee” in the CR spectrum around ~3 PeV suggests a cutoff or transition in the astrophysical systems producing CRs. Direct investigation of the source of these CRs at extreme energies is complicated by the redirection of charged CR particles by interstellar magnetic fields over galactic distance. Understanding these astrophysical systems requires very high-energy (VHE; *E* ~ 0.1–100 TeV) and ultrahigh-energy (UHE; *E* > 100 TeV) gamma-ray telescopes and neutrino detectors that detect uncharged particles that travel unimpeded from galactic sources to Earth.

The Large High Altitude Air Shower Observatory (LHAASO; Z. Cao et al. 2019) observes VHE and UHE gamma rays with energies well into the petaelectronvolt range. Over the past few years, LHAASO has detected 90 sources, 43 reaching UHE energies (Z. Cao et al. 2024). LHAASO has two component detector arrays: the Water Cherenkov Detector Array (WCDA), operating at 0.1–30 TeV, and the Square Kilometer Array (KM2A), at 30 TeV–30 PeV.

In this work, we examine one particular region of VHE emission near R.A. = 55°,5, decl. = 53°.1. Originally described in Z. Cao et al. (2021) as a single extended source, named LHAASO J0341+5258, further analysis with additional data and better sensitivity discerned two KM2A and one WCDA sources in this region of the galactic plane (Z. Cao et al. 2024). In the first LHAASO catalog, 1LHAASO J0343+5254u (previously LHAASO J0341+5258) has a reported test statistic (TS) for photons above 100 TeV of 20.2, suggesting that the underlying astrophysical system is a PeVatron (due to the gamma-ray emission above 100 TeV) and leading to substantial interest in follow-up observations and analyses.

In the original discovery paper (Z. Cao et al. 2021) and in subsequent independent analyses (P. Bangale & X. Wang 2023; A. De Sarkar & P. Majumdar 2024), special attention was given to a nearby Fermi Large Area Telescope (Fermi-LAT) point source, 4FGL J0340.4+5302. In the 4FGL-DR4 catalog (S. Abdollahi et al. 2022; J. Ballet et al. 2023), 4FGL J0340.4+5302 is “unassociated,” having no known counterpart at other wavelengths. Because of the sharp cutoff in the Fermi-LAT spectrum around ~2 GeV, these previous papers have treated this Fermi source as a possible gamma-ray pulsar, but without radio pulsations or gamma-ray timing solutions this classification is tentative. A. De Sarkar & P. Majumdar (2024) discussed the Fermi-LAT gamma-ray emission at 4FGL J0340.4+5302 and made a case for that pulsar-like source being the lower-energy counterpart to the LHAASO source, establishing an upper limit on the X-ray flux of approximately 10⁻¹² erg s⁻<fnref n="1" /> cm⁻<fnref n="2" /> at the Fermi source using brief archival Chandra data. Fortunately, there are two Swift X-ray Telescope (D. N. Burrows et al. 2005) observations of this Fermi source from 2012, totaling ~5 ks,

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48820, USA</footnote>

<footnote n="2">²Faculty of Science, Kanagawa University, 3-27-1 Rokukakubashi, Kanagawa-ku, Yokohama-shi, Kanagawa 221-8686, Japan</footnote>

<footnote n="3">³Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027, USA</footnote>

<footnote n="4">⁴FSLAC IRL 2009, CNRS/IAC, La Laguna, Tenerife, Spain</footnote>

<footnote n="5">⁵AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, F-91191 Gif sur Yvette, France</footnote>

<footnote n="6">⁶Department of Physics and Astronomy, University of Manitoba, Winnipeg, MB, R3T 2N2, Canada</footnote>

<footnote n="7">⁷Department of Applied Physics, Faculty of Engineering, Kanagawa University, 3-27-1 Rokkakubashi, Kanagawa-ku, Yokohama, Kanagawa 221-8686, Japan
Received 2024 December 24; revised 2025 February 18; accepted 2025 February 18; published 2025 April 2</footnote>
available in the HEASARC archive, which can give additional context to this Fermi-LAT source.

Summing these two observations with XImage v4.5.1, we find no X-ray point or extended source in the 95% Fermi-LAT uncertainty ellipse, for a $3\sigma$ upper limit for 0.3–10.0 keV absorbed flux of $\approx 1 \times 10^{-13}$ erg s$^{-1}$ cm$^{-2}$ (assuming $\Gamma_X = 2$). Notably, this upper limit is below the X-ray synchrotron flux predicted by A. De Sarkar & P. Majumdar (2024) in both their leptonic and hadronic models for the combined multiwavelength (MWL) spectra. Ongoing work on additional Chandra observations that partially overlap with the 95% ellipse of 4FGL J0340.4+5302 has detected dim X-ray point sources with fluxes of $\sim 10^{-13}$ erg s$^{-1}$ cm$^{-2}$ in the Fermi-LAT region (F. Acero et al. 2025, in preparation, private correspondence), which will certainly increase our understanding of 4FGL J0340.4+5302 as a component of the complex gamma-ray landscape near 1LHAASO J0343+5254u.

It is still uncertain what electromagnetic and CR emission model is most appropriate for UHE sources like 1LHAASO J0343+5254u. Leptonic models—dominated by electron synchrotron inverse-Compton (IC) upscattering of ambient photons by those electrons, or synchrotron self-Compton (SSC)—can explain a wide range of high-energy behaviors using a single population of parent particles in systems like pulsar wind nebulas (PWNs). Alternatively, hadronic processes like pion decay and CR interactions in supernova remnants can account for VHE and UHE emission, and hadronic processes are required to produce the neutrino emission that is suspected to accompany galactic PeVatrons. The discovery of neutrino flux contained in the plane of the Milky Way by the IceCube observatory (R. Abbasi et al. 2023) suggests that hadronic processes are at work in galactic VHE/UHE systems, as neutrinos are produced in copious amounts only by hadronic processes.

A notable challenge for understanding the MWL behavior of galactic PeVatrons is the difficulty in confident association with lower-energy sources. Not only are the uncertainty ellipses or spatial extensions of UHE sources quite large, but their locations in the galactic plane guarantee a bevy of gamma-ray and X-ray sources coincident but unrelated to the VHE and UHE emission. Furthermore, it is entirely possible for a source like a pulsar to be bright enough for detection in gamma rays but too dim for detection in X-rays, or vice versa (D. A. Smith et al. 2023), so an X-ray source may not necessarily be related to a nearby VHE/UHE emission region. Special care should therefore be taken to coordinate low-energy follow-up to gamma-ray emission to search for possible counterparts.

In this work, we report on XMM-Newton (F. Jansen et al. 2001) observations from 2024 February targeting 1LHAASO J0343+5254u (observation IDs 0923400401, 0923400801, and 0923401401). In Section 2, we describe our observations, data reduction, and fitting of X-ray and gamma-ray data. In Section 3, we discuss MWL leptonic modeling; and in Section 4, we discuss our X-ray and MWL modeling results and a related radio project. Finally, in Section 5, we present our conclusions and discuss next steps.

## 2. Observations and Data Reduction

We downloaded the three observations (0923400401, 0923400801, and 0923401401) from the XMM-Newton Science Archive,<fnref n="8" /> and used the standard XMM-Newton

analysis pipeline implemented in Python `pysas v21.0.0` (C. Gabriel et al. 2004) to produce refined event files for each observation. Starting with `epproc` and `emproc`, we produced first-look light curves of the entire XMM-EPIC field of view, to evaluate particle flaring times and establish a desirable rate cut. We found that unfiltered per-second event rates below 40 cps for pn data and below 4 cps for MOS1 and MOS2 data excluded particle flares and times with high background. A manual time interval filter was only necessary for observation 0923400801, due to a low plateau in the event rate surrounding a smaller flare. These rate cuts create 87.7 ks of good time between the three observations.

Applying these rate cuts, we used `tabgtigen` to create good time intervals for each observation with additional filters, using event patterns 0 through 4 and energies $0.2 \text{ keV} < E < 15 \text{ keV}$ for pn data and patterns 0 through 12 and energies $0.2 \text{ keV} < E < 12 \text{ keV}$ for MOS1 and MOS2 data. Finally, we use `evselect` to create cleaned, filtered event lists for each observation.

Examining the XMM-EPIC-MOS2 data by eye, we noticed a substantial extended region of X-ray emission in the southwest corner of the 1LHAASO J0343+5254u WCDA region. By XMM-Newton conventions, this “primary” source is given the name XMMU 034124.2+525720, shown in the insert to Figure 1.

In the XMM-Newton EPIC-MOS2 field of view, we noted several secondary X-ray point sources besides the extended primary region. These secondary sources deserve additional attention, especially as they are scattered across the entire region of interest for 1LHAASO J0343+5254u. Using the exposure-corrected and filtered EPIC-MOS2 data, we use the XIMAGE routines `detect` and `sosta` to identify all sources with a signal-to-noise ratio (S/N) > 4 in the MOS2 data, finding 11 in the field of view upon excluding the extended primary source region described above. We label these sources 2 through 12, in no particular order.

To generate the spectrum of the primary source, we used a circular selection region centered on the maximum count rate in the primary region, at R.A. = 55.3509, decl. = 52.9527. We set the size of the primary region at $0.03 = 1/8$, encompassing the entirety of the extended emission while not extending to either of the nearby point sources. This size is similar to other X-ray PWNs related to teraelectronvolt emission regions, like those linked to G75.2+0.1 (J. Woo et al. 2023, having X-ray size $\approx 1'$) and HAWC J1826-128 (D. A. Burgess et al. 2022, $\approx 3'$), an indication that this new source may be a PWN like those objects.

We select a nearby background region distant from the primary source in order to dodge any contamination from the outer limits of the primary extended emission. The background region was selected to be entirely contained in both the EPIC MOS1 and MOS2 fields of view, located at R.A. = 55.5050, decl. = 52.0523, with radius 0.055. Using `evselect`, we extracted background and source spectra for the primary and secondary sources using the clean, filtered event lists, generating `.arf` and `.rmf` files with `arfgen` and `rmfgen`, respectively.

We used `epicspeccombine` to separately create summed EPIC-MOS and EPIC-pn spectra (it is not feasible to merge EPIC-MOS and EPIC-pn data into a single file). For some of the secondary point sources, their spatial location placed them outside of the XMM-EPIC-MOS1 field of view, so only MOS2

<footnote n="8">⁸ The XMM-Newton Science Archive data become public on 2025 March 25.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The summed XMM-EPIC MOS1 and MOS2 rate map, showing the LHAASO regions (white dashed lines—LHAASO J0341 is the source described in Z. Cao et al. 2021, while the others are the KM2A and WCDA sources from Z. Cao et al. 2024), the primary source region (green line), the selected background region (yellow dashed line), the nearby 4FGL 95% uncertainty ellipse (red line), and the 11 detected secondary point sources (numbered in cyan). The lime insert shows the primary source in greater detail, along with the fitted absorbed power-law slopes in each annular analysis region.</figcaption></figure>

data were used for those. Finally, we used the *ftools* function *ftgrouppha* to group each spectrum with minimum S/N > 3 in each grouped energy bin. This final step reduces the number of bins for fitting but guarantees that each energy bin would individually be a substantial detection and that each bin has approximately Gaussian errors appropriate for χ² fitting.

## 2.1. Primary Source Extension

To better constrain the spatial extension of the primary source, we used the surface brightness analysis tool in *ds9* to extract a surface brightness profile in the four cardinal directions from the center of the primary source in the summed EPIC-MOS count-rate data. We extract the per-square-arcsecond count rate in rectangular tracks going north, east, south, and west from the centroid of the primary source. Each rectangular track has width 0.6, to average out small-scale variations.

Figure 2 shows that the surface brightness of the extended primary source drops down to background levels around ~1' from the centroid position. However, the southern portion of the primary source extends farther, out to ~2', while the northern portion has a steeper drop in surface brightness than any other direction, establishing the asymmetric spatial shape of the primary source. The EPIC-MOS point-spread function

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The average flux density along four tracks extended in the cardinal directions (east—red solid; west—blue dashed; north—lime dotted—dashed; and south—cyan dotted) from the primary source center, along with the EPIC-MOS PSF for a 10' offset and 3 keV energy (black). The insert shows the regions used in the extraction of these count rates. The uncertainties (black for east, cyan for south) assume Poisson statistics and are displayed every 0.36'.</figcaption></figure>
<table><caption>**Table 1**
Spectral Fitting Results for the Primary Source</caption><thead><tr><th></th><th>Exposure<br>(ks)</th><th>χ<sup>2</sup></th><th>D.o.F.</th><th>nh<br>(10<sup>22</sup> cm<sup>-2</sup>)</th><th>Γ<sub>x</sub></th><th>F<sub>x</sub>/10<sup>-12</sup><br>(erg s<sup>-1</sup> cm<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Initial</td><td></td><td></td><td></td><td>0.741</td><td>2</td><td>1</td></tr><tr><td>All Observations</td><td>87.7</td><td>497.76</td><td>491</td><td>1.68 ± 0.12</td><td>1.91 ± 0.08</td><td>1.29 ± 0.05</td></tr><tr><td>0923401401</td><td>43.0</td><td>601.45</td><td>572</td><td>1.65 ± 0.18</td><td>1.90 ± 0.12</td><td>1.24 ± 0.06</td></tr><tr><td>0923400401</td><td>21.6</td><td>287.96</td><td>304</td><td>1.75 ± 0.26</td><td>1.98 ± 0.18</td><td>1.25 ± 0.09</td></tr><tr><td>0923400801</td><td>23.1</td><td>386.61</td><td>342</td><td>1.81 ± 0.31</td><td>1.98 ± 0.20</td><td>1.18 ± 0.10</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** 0.3–8.0 keV XMM-Newton spectrum of the primary source, fitted with an absorbed power law. An instrumental feature of the EPIC-MOS detector at 1.7 keV is masked out for our fitting.</figcaption></figure>

(PSF) is estimated using the XMM-Newton calibration file for the MOS camera (XRT1_XPSF_0016.CCF), using the King profile description at a 10' distance from the center of the field of view at 3 keV. It is much smaller than the extension seen for the primary source, as shown by the black line in Figure 2.

## 2.2. Primary Source Fitting

Using `Xspec v.12.14.0h` (K. A. Arnaud 1996), we conducted an absorbed power-law fit to the background-subtracted primary source spectrum for each individual observation and for the summed spectra. For this fitting, we load both the EPIC-MOS and EPIC-pn data into `xspec` and conduct a simultaneous fit of the EPIC-MOS and EPIC-pn data. An instrumental silicon feature appears in the EPIC-MOS data as a deficit in the spectrum at approximately 1.7 keV and an instrumental 8 keV feature appears in the EPIC-pn data. These features are due to the background spectrum of the cameras; by selecting a background region in a different CCD panel than the source, an instrumental excess feature in the background CCD becomes a deficit in the spectrum. We mask out data points between 1.65 and 1.75 keV for EPIC-MOS and between 7.5 and 8.5 keV for EPIC-pn to avoid any impact on our results.

The fit was of the form `tbabs * cflux * powerlaw`, with `tbabs` giving X-ray absorption by the ISM, `cflux` calculating the 0.2–12.0 keV energy flux, and `powerlaw` being a simple power law with photon index $\Gamma_x$. For this fit, the galactic $n_H$ column density was started at the cataloged value ($n_H = 0.741 \times 10^{22} \text{cm}^{-2}$; HI4PI Collaboration et al. 2016) but

allowed to vary as a free parameter. For `tbabs`, we use cross-sectional data from D. A. Verner et al. (1996) and the abundance calibration from J. Wilms et al. (2000).

The spectrum of the primary source, shown with a power-law fit for the summed observations in Figure 3, appears as a featureless power law with photon index $\Gamma_x = 1.9 \pm 0.1$. The fitted parameters for the primary source are shown in Table 1, with 95% uncertainties generated by the `error` command. The fits from each individual observation are within the errors of the summed fit, so there is no evidence to suggest any interobservation variability in the X-ray spectrum of this source. We find no evidence for a high-energy cutoff or break in the power-law spectrum of the primary source, suggesting that the peak of the $\nu F_\nu$ flux is near or above the top of XMM's energy band.

To test whether $\Gamma_x$ varies based on the distance from the central, brightest part of the primary source, we also subdivided the primary source into annular regions concentric with the center of the region and performed spatially resolved spectral analysis. We created a core region of radius 0.6, an inner annulus beginning there and extending to 1/2, and an outer annulus from 1/2 to 1/8, shown in the insert of Figure 1. We conducted identical power-law fits using all XMM-EPIC data for these three separate regions. We fix $n_H = 1.68 \times 10^{22} \text{cm}^{-2}$ as the value from the overall fit.

For the core, inner, and outer regions, we find photon indices $\Gamma_x$ of $1.80 \pm 0.09$, $1.99 \pm 0.09$, and $2.27 \pm 0.15$, respectively. These fits show a statistically significant difference between the power-law slopes of the core part of the primary source and the annular outer part, demonstrating that the outermost regions of the primary source have notably softer X-ray spectra.

To facilitate MWL fitting of this system, we extract de-absorbed $E^2 \frac{dN}{dE}$ fluxes for the fully fitted all-observation primary model, by setting $n_H = 0$ and using four logarithmic energy bins between 2 keV and 8 keV, with proportional uncertainties obtained via `xspec` for the model flux. These fluxes can be incorporated into a wider spectral energy distribution (SED) spectrum and compared to other fluxes directly.

Conducting timing analysis of the innermost 0.6 region of the primary source, we detect no evidence for X-ray pulsations coming from the region. However, the pn-mode data taken by XMM-Newton are not in the mode most appropriate for superior timing sensitivity, so additional X-ray observations with more discerning timing registrations are warranted.

## 2.3. Secondary Sources

We fit the EPIC-MOS data for each secondary source with an absorbed power-law model to evaluate the general shape of each spectrum. In this case, we kept $n_H$ fixed to the previous fitted value of $1.68 \times 10^{22} \text{cm}^{-2}$. Cataloged $n_H$ values do not
<table><caption>**Table 2**

Observational Parameters, Background-subtracted Count Rate in XMM-EPIC, and Power-law Fits for the Secondary Sources</caption><thead><tr><th>ID</th><th>R.A. (J2000)<br>(HH:MM:SS)</th><th>Decl. (J2000)<br>(dd:mm:ss)</th><th>Count Rate<br>(1 ks<sup>-1</sup>)</th><th>log F<sub>x</sub>, Γ<sub>x</sub> = 2<br>(erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>log F<sub>x</sub><br>(erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>Γ<sub>x</sub></th><th>Notes</th></tr></thead><tbody><tr><td>2</td><td>03:42:23.2</td><td>53:07:34.9</td><td>1.2 ± 0.2</td><td>-13.5 ± 0.1</td><td>-13.5 ± 0.1</td><td>1.8 ± 0.5</td><td></td></tr><tr><td>3</td><td>03:41:28.3</td><td>53:00:45.4</td><td>4.2 ± 0.2</td><td>-12.9 ± 0.1</td><td>-12.9 ± 0.1</td><td>1.7 ± 0.2</td><td></td></tr><tr><td>4*</td><td>03:42:12.1</td><td>53:11:42.4</td><td>1.0 ± 0.1</td><td>-13.1 ± 0.1</td><td>-13.2 ± 0.1</td><td>1.4 ± 0.7</td><td>EPIC-MOS2 only</td></tr><tr><td>5</td><td>03:42:35.5</td><td>53:11:30.2</td><td>1.5 ± 0.2</td><td>-13.2 ± 0.1</td><td>-13.1 ± 0.2</td><td>1.0 ± 0.7</td><td></td></tr><tr><td>6</td><td>03:43:15.7</td><td>53:03:33.4</td><td>0.8 ± 0.2</td><td>-12.8 ± 0.2</td><td>-13.0 ± 0.4</td><td>1.2 ± 2.0</td><td></td></tr><tr><td>7</td><td>03:41:10.8</td><td>53:07:35.1</td><td>1.0 ± 0.2</td><td>-13.4 ± 0.2</td><td>-13.3 ± 0.5</td><td>0.3 ± 1.4</td><td>ROSAT X-ray source</td></tr><tr><td>8</td><td>03:41:36.5</td><td>52:54:49.8</td><td>3.4 ± 0.2</td><td>-13.0 ± 0.1</td><td>-12.0 ± 0.3</td><td>4.0 ± 0.5</td><td>V = 12.5 optical loading</td></tr><tr><td>9</td><td>03:41:01.0</td><td>53:06:54.3</td><td>0.7 ± 0.2</td><td>-13.3 ± 0.2</td><td>-12.5 ± 1.5</td><td>3.21 ± 1.9</td><td></td></tr><tr><td>10</td><td>03:41:05.1</td><td>53:09:26.0</td><td>0.8 ± 0.2</td><td>-13.4 ± 0.2</td><td>-13.0 ± 0.6</td><td>0.2 ± 1.6</td><td></td></tr><tr><td>11*</td><td>03:42:17.8</td><td>53:17:12.6</td><td>0.6 ± 0.1</td><td>-13.0 ± 0.2</td><td>-13.1 ± 0.4</td><td>0.9 ± 0.9</td><td>EPIC-MOS2 only</td></tr><tr><td>12</td><td>03:40:47.4</td><td>53:06:48.7</td><td>2.1 ± 0.2</td><td>-13.1 ± 0.1</td><td>-13.1 ± 0.1</td><td>1.9 ± 0.4</td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** 0.5–500 GeV Fermi-LAT TS map after removing 4FGL J0340.4+5302 from the gamma-ray background model. We include the 1LHAASO J0343+5254u WCDA and KM2A regions (white), the 95% error ellipse of 4FGL J0340.4+5302 (blue ellipse), and the extent of the XMM-Newton primary source (small blue circle). The red points are two X-ray point sources described in upcoming works. NRO CO line contours are shown for velocity ranges of [−20, −4] (yellow) and [−4, +10] (green) km s⁻¹.</figcaption></figure>

substantially vary over the field of our observations, though some of these secondary X-ray sources may be extragalactic, which could impose additional $n_H$ column density. For a first-look flux, we fixed $\Gamma_X = 2$ and subsequently freed that parameter for a full fit. We also note the 0.2–12.0 keV background-subtracted count rate for each source, which is a more observational measurement of the photon flux from each secondary source.

In Table 2, we report the fitted fluxes and power-law slopes for each of the detected secondary sources. Some of the dimmer sources have very soft fitted spectral indices, leading to unrealistic and disproportionate values for their flux; with such a steep spectrum, more photons would be emitted below 1 keV, where X-ray absorption is substantial, so the calculated unabsorbed flux via *cflux* would be unrealistically high. The reported count rate in Table 2 is a reasonable proxy for flux in these cases.

Conducting a cross-reference with the 2 µm All-Sky Survey (M. F. Skrutskie et al. 2006) and with the SIMBAD database, none of the secondary sources are coincident with IR sources down to $m_J = 15$ or SIMBAD optical sources, with the exception of source 8, southeast of the primary. That source in particular is coincident with the star UCAC4 715-027962 of

magnitude $V = 12.5$, listed as a spectroscopic binary in SIMBAD. This star is near the brightness that would incur optical loading in the XMM-EPIC detector in full-frame mode, so its detection as an X-ray source is suspect. Additionally, using a general cross-reference with the SIMBAD catalog, source 7 is noted as a ROSAT X-ray point source. Its reported flux in the ROSAT dim-source catalog (W. Voges et al. 2000) is substantially higher than our observed flux in XMM, suggesting a degree of X-ray variability. These secondary sources may be interesting targets for further follow-up.

## 2.4. Fermi-LAT Data Analysis

We also analyzed Fermi-LAT photon data in the sky region surrounding the XMM-Newton position. The nearby Fermi 4FGL source J0340.4+5302 is not spatially coincident with our XMM-Newton primary source, but could be related to a relic radio PWN related to our candidate X-ray PWN by the movement of a pulsar with a natal kick. We evaluate the spatial extent of the 4FGL source and obtain upper limits at the XMM-Newton position. We performed data reduction, event selection, and time filtering using *fermitools* version 2.2.0 and *fermipy* version 1.1.6 (M. Wood et al. 2017).

We used Fermi-LAT photon data from 2008 August to 2024 August in the 0.5–500 GeV energy range, in a circular region with radius 20° centered at the XMM-Newton primary source. We performed a summed likelihood analysis using the PSF event types, splitting Fermi-LAT events according to their quality of direction reconstruction (PSF0 being the worst reconstruction). We used the following selection as a function of energy: PSF2+3 from 0.5–1 GeV and PSF0+1+2+3 for 1–500 GeV, to keep only the best reconstructed events at low energies. The Fermi-LAT instrument response function version P8R3-SOURCE-V3 was used in the analysis, with the maximum value of the zenith angle of observation being 90°.

We included all 4FGL-DR4 (J. Ballet et al. 2023) gamma-ray point sources and extended sources in the gamma-ray background model. For the diffuse background model, we used the Galactic diffuse emission template (gll_iem_v07.fits) and the relevant extra-Galactic isotropic emission templates for each PSF type (in the form of iso_P8R3_SOURCE_V3_PSF*_v1.txt, where * is the PSF type described above), obtained at the Fermi science center (Goddard Space Flight Center 2024).

To visualize the spatial morphology of 4FGL J0340.4+5302, we deleted that source from the gamma-ray background and reproduced the TS map in Figure 4. Besides the 95% uncertainty ellipse of 4FGL J0340.4+5302, Figure 4 also
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Detections (green dots) and constraining upper limits (black arrows) for our PWN candidate, along with the MCMC best-fit leptonic model (blue). The faint red curve is the leptonic IC spectrum, but with $κ = 1$ photon fields fixed at solar-neighborhood values. The components of the gamma-ray emission are included as thin lines of the same color for each model (SSC is not visible on this scale).</figcaption></figure>

shows additional excess gamma-ray emission around the positional error ellipse of 1LHAASO J0343+5254u in white.

We measured the extension of 4FGL J0340.4+5302 with the extension tool under two different spatial models, Radial-Gaussian and RadialDisk, measuring a likelihood ratio of extension TS<sub>ext</sub> = -2 log L<sub>ext</sub>/L<sub>pt</sub>. If TS<sub>ext</sub> > 16, the tested source is likely extended. The models have TS<sub>ext</sub> = 2.19 and = 2.18, respectively, suggesting that 4FGL J0340.4+5302 is a point source in the energy range of 500 MeV–500 GeV. We did find that TS<sub>ext</sub> increased substantially, depending on the low-energy limit of our event selection, suggesting that there may be substantial energy dependence on the morphology of gamma-ray emission in this region.

The best-fit position of the point-like source with a log-parabola (LP)-type spectrum was calculated using the localize method from the fermipy package. The best-fit position of 4FGL J0340.4+5302 for 500 MeV–500 GeV was found to be R.A. (J2000), decl. (J2000) = 55.1426 ± 0.0174, 53.0656 ± 0.0181, which is compatible with the location reported for 4FGL J0340.4+5302 in the 4FGL-DR4 catalog.

In order to model the gamma-ray flux at the position of our XMM-Newton primary source, we use a two-source model for Fermi-LAT emission in the region, with one source at the position of 4FGL J0340.4+5302 and another at the position of our XMM-Newton primary source. The gamma-ray spectrum of 4FGL J0340.4+5302 was modeled as an LP, the model best fitting its spectrum in the 4FGL-DR4 catalog (J. Ballet et al. 2023), with the form

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-(\alpha + \beta \log(E/E_0))} $$

with $E_0$ being the reference energy and $E_b$ the break energy. The spectral indices are $\alpha$ and $\beta$. The hypothetical additional

source at the XMM-Newton position was modeled with a power law with $Γ_γ = 2$.

In the two-point-sources scenario, 4FGL J0340.4+5302 was detected with a TS value of ~834 (28.9σ), with spectral parameters $α = 3.51 ± 0.27$ and $β = 0.747 ± 0.352$. The total energy flux was found to be $(5.4 ± 0.3) × 10^{-6}$ MeV cm$^{-2}$ s$^{-1}$. No significant gamma-ray emission was detected at the position of the XMM-Newton primary source, so we evaluated the upper limits for that position by assuming a power-law spectrum with $Γ_γ = 2.0$. We find an upper limit on 0.5–500 GeV flux at the position of the XMM-Newton primary source of $6.9 × 10^{-10}$ erg s$^{-1}$ cm$^{-2}$, and we also establish upper limits on smaller energy ranges, as shown in Figure 5.

### 3. MWL Modeling

In this section we combine the X-ray spectrum of the XMM-Newton primary source with the Fermi-LAT upper limits and existing VHE data points from various observatories, creating a MWL SED stretching approximately 12 orders of magnitudes in energy.

#### 3.1. MWL Spectrum Construction

The primary XMM-Newton source has a featureless power-law X-ray spectrum for the overall region, and asymmetrical spatial extensions reaching arcminutes away from a central bright spot. These X-ray properties point toward classification as an X-ray PWN, similar to other systems, like HAWC J1826-128 (the Eel; D. A. Burgess et al. 2022), G106.3+2.7 (the Boomerang; I. Pope et al. 2024), and G75.2+0.1 (the Dragonfly; J. Woo et al. 2023). Given the complex gamma-ray sky around 1LHAASO J0343+5254u, it is possible that the XMM primary source is related to the high-energy emission
<table><caption>**Table 3**
MCMC Best-fit Parameters and $2\sigma$ Uncertainties for the Leptonic Model</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Best Fit</th></tr></thead><tbody><tr><td>E<sub>tot,e</sub></td><td>erg</td><td>(2.1 ± 0.3) × 10<sup>44</sup></td></tr><tr><td>E<sub>cut,e</sub></td><td>TeV</td><td>130 ± 19</td></tr><tr><td>α</td><td></td><td>2.0 ± 0.2</td></tr><tr><td>B</td><td>μG</td><td>5.1 ± 2.3</td></tr><tr><td>κ</td><td></td><td>61 ± 15</td></tr></tbody></table>

observed by LHAASO like the other gamma-ray PWN noted above.

To create an MWL spectrum and test whether the primary source is feasibly related to the LHAASO emission as a PWN, we join our measured X-ray flux points from Section 2.2 with the Swift Burst Alert Telescope (Swift-BAT) upper limits from the nondetections with Swift-BAT (K. Oh et al. 2018), the derived Fermi-LAT upper limits at the position of the primary source from Section 2.4, and the LHAASO spectrum from Z. Cao et al. (2021), to form an SED spanning 9 orders of magnitude in energy. We also incorporate the VERITAS upper limits from P. Bangale & X. Wang (2023) in our model, but those VERITAS upper limits are modeled on a smaller area of the sky, containing only about ~25% of the LHAASO source region. To be as conservative as possible while still including as much information as possible in our spectrum, we multiply the limits reported in P. Bangale & X. Wang (2023) by four, which is a workaround that in all ways treats the limits as conservatively and cautiously as possible; the flux observable with VERITAS is certainly below our recalculated limits. These VERITAS limits, given these caveats, do not in the end affect the leptonic modeling described below.

The LHAASO spectrum reported in Z. Cao et al. (2021) has been slightly updated in the first LHAASO catalog (Z. Cao et al. 2024), with a flux and slope, but not with detailed points as in Z. Cao et al. (2021). The fluxes and slopes reported in Z. Cao et al. (2024) are similar to the overall spectrum described in Z. Cao et al. (2021), a work that included substantially more detailed explanations of the spectrum of this region. We adopt the spectrum from Z. Cao et al. (2021), acknowledging that the gamma-ray emission has been decomposed into subsources in Z. Cao et al. (2024).

We manually searched optical and IR catalogs for low-energy points or extended counterparts to the candidate PWN but found no counterparts. The MWL spectrum constructed from the XMM-Newton and LHAASO detections and other upper limits is shown in Figure 5.

## 3.2. Leptonic Model Fitting

The electromagnetic emissions from systems like HESS J1825-137 (H. E. S. S. Collaboration 2019), the Eel (D. A. Burgess et al. 2022), and the Boomerang (I. Pope et al. 2024) have been successfully modeled as fully leptonic. These models invoke electron synchrotron radiation in a magnetic field of magnitude ~1–10 μG for the X-ray radiation, plus gamma-ray emission via IC scattering of the ambient cosmic microwave background, and near-IR (NIR) and far-IR (FIR) photon fields (plus, to a lesser extent, SSC upscattering). Conducting phenomenological fitting using a synchrotron + IC

model can test whether there are feasible leptonic scenarios to explain the MWL emission from our candidate PWN.

We use the naima Python package (V. Zabalza 2015) to construct a leptonic model for the MWL spectrum and constraining upper limits. We adopt an exponential cutoff power law for the electron population,

$$ \frac{dN}{dE} = N_0(E/E_0)^{-\alpha}e^{-(E/E_{\text{cut,e}})^{\beta}}, $$

between electron energies of 100 MeV and 10 PeV. $N_0$ is an amplitude for the electron distribution, and we set the reference energy $E_0 = 1$ TeV and $\beta = 1$. The electron distribution has three free parameters: the power-law slope $\alpha$ and cutoff energy $E_{\text{cut,e}}$, plus $N_0$. We fix the distance to the system to 1 kpc, a galactic distance similar to the characteristic distances of the spiral arms in the outer galaxy, discussed in Section 4.3, and four of the newly discovered radio sources, discussed in Section 4.4.

The synchrotron emission occurs in a region with an ambient magnetic field $B$ of size 0.5 pc ($\frac{d}{kpc}$) (the physical size corresponding to $2'$ at distance $d$). The IC scattering (D. Khangulyan et al. 2014) occurs in a wider area populated with ambient cosmic microwave background and locally produced synchrotron X-ray photons (which only minimally contribute to upscattering), plus FIR (20 K, $\rho_E = 0.3\kappa$ eV cm$^{-3}$) and NIR/optical (5000 K, $\rho_E = 1.0\kappa$ eV cm$^{-3}$) photon fields. naima integrates the upscattering of these photons by electrons with energy between 1 GeV and 510 TeV. Initially, we fixed the intensities of these ambient photon fields to the solar-neighborhood values. In that case, the red line in Figure 5 shows that only ~10% of the observed LHAASO flux can be explained by IC emission. Many other leptonic models allow the field intensity to be a free parameter, so we multiply the energy densities $\rho_E$ of the NIR and FIR fields by a final fitted parameter $\kappa$. For this IC regime, the FIR field is substantially more relevant than the NIR field, and scaling both densities together ensures that the relative temperature of the overall field stays constant.

As a starting point for fitting, we adopt $\alpha = 2$, $E_{\text{cut,e}} = 100$ TeV, $B = 4.5$ μG, and $\kappa = 50$. While this value for $\kappa$ means a substantially higher FIR field compared to the solar neighborhood, elevated FIR fields up to 2 orders of magnitude higher than local values are not unusual for leptonic models for IC emission from PWNs (D. A. Burgess et al. 2022), and can be achieved in star-forming or gas-/dust-rich regions of the galactic spiral arms.

To more fully explore the parameter space of this model, we use the emcee package to conduct a Markov Chain Monte Carlo (MCMC) ensemble sampling. This simulation uses “walkers,” each starting at a slightly different initial condition, which iteratively explore the multidimensional parameter space of the model. Eventually, the walkers produce an ensemble of parameter vectors that can constrain the best parameters for the model, visualized with a pairs plot. We use an MCMC ensemble sampling routine with 100 individual walkers, running for 200 steps after a 50 step burn-in, to characterize the parameter space of the leptonic model.

For the leptonic fit, the MCMC sampler obtained a best fit of $\alpha = 2.0 \pm 0.2$, $E_{\text{cut,e}} = (130 \pm 19)$ TeV, $B = (5.1 \pm 2.3)$ μG, and $\kappa = 61 \pm 15$ (the amplification of the ambient IR fields by a factor of 61 over the local value in the solar neighborhood). The uncertainty on $B$ is high, due to degeneracies between $B$ and other parameters, including $\alpha$ and $\kappa$. The amplitude of the lepton distribution is $N_0 = 7.0 \times 10^{30}$ eV$^{-1}$, for a total electron
energy of $2.1 \times 10^{44}$ erg. These parameters are summarized in Table 3, and this model is shown in blue in Figure 5.

## 4. Discussion

In terms of X-rays, the primary XMM source has spatial and spectral features much like other PWNs—like the Eel, Boomerang, and Dragonfly, all of which have been linked to nearby >teraelectronvolt emission regions observed in VERITAS (I. Pope et al. 2024), HAWC, and HESS (D. A. Burgess et al. 2022). In these systems, the nearby >teraelectronvolt emission regions are generally more extended than the X-ray emission, which is not always centrally located in the gamma-ray region, similar to our candidate PWN. In this way, our primary source may be linked to the wider LHAASO emission region by the same mechanisms.

A useful comparison is the extended PWN complex HESS J1825-137, described in detail in H. E. S. S. Collaboration (2019). In that work, gigaelectronvolt and teraelectronvolt emission extending out to ~1° around a region of X-ray emission is ascribed to a PWN. The X-ray-emitting nebular region in that work has a compact core and extends asymmetrically to a few arcminutes, similar to our primary X-ray source. In that system, the X-ray emission is ascribed to synchrotron radiation from ~teraelectronvolt leptons in a 3 µG magnetic field, while gamma-ray emission is produced from IC scattering of ambient photons off >gigaelectronvolt leptons farther from the central X-ray region. A similar leptonic model could apply to the joint emission of our candidate PWN and the VHE/UHE emission observed with LHAASO. Though the emission from our candidate PWN does have some substructure, as shown in the insert in Figure 1, our XMM observations do not have sufficient spatial resolution to decisively identify a point source in the candidate PWN that may be the pulsar itself.

In the 2021 paper that originally described gamma-ray emission in this region (Z. Cao et al. 2021), the LHAASO collaboration presented leptonic and hadronic fits to the LHAASO data, using only weak X-ray constraints and Fermi-LAT data relating to the nearby 4FGL J0340.4+5302. To obtain an appropriate fit to the observed gamma-ray cutoff, they used an exponential power-law particle distribution with $\beta = 2$, as opposed to $\beta = 1$ in this work. Given that the 2021 LHAASO data will likely be updated, we do not comment on the wisdom of modifying $\beta = 1$, only noting that an updated LHAASO spectrum will help constrain whether a sharper cutoff in the particle spectrum is warranted. Allowing $\beta > 1$ would better fit the gamma-ray cutoff presented in Z. Cao et al. (2021) but would necessitate an increased B-field magnitude for the synchrotron emission as well.

Extrapolating the models presented in Z. Cao et al. (2021) to kiloelectronvolt energies results in a predicted X-ray spectrum with a peak around ~1 keV. Because our candidate PWN has an X-ray spectrum that must peak substantially higher, modifications to the models proposed in Z. Cao et al. (2021) are warranted.

### 4.1. Leptonic Model Analysis

The blue line in Figure 5 shows the MCMC best-fit leptonic model, along with gamma-ray constituent components, in thinner lines of the same color. The faint red gamma-ray curve shows the IC emission in the leptonic case if $\kappa = 1$, the NIR/FIR fields fixed at the solar value. In the best-fit leptonic model, the values for $\alpha$, $E_{cut,e}$, and B are typical of PWNs related to

VHE emission regions, like the Eel. As mentioned previously, an elevated $\kappa$ would be readily available in a dusty or gas-filled region of space, like a star-forming region, as optical light is reprocessed to lower energies. Low-energy observations near this new candidate PWN will more greatly constrain astrophysical features that might contribute to the ambient IR field.

In the leptonic model, the limits prescribed by Fermi-LAT mean that the electron distribution must have $\alpha \lesssim 2$, a value similar to other VHE PWNs. The featureless X-ray spectrum observed in our XMM-Newton data constrains the magnetic field B producing the synchrotron emission to be high, with $B \approx 5$ µG, akin to other VHE PWNs, though there are substantial uncertainties in our fitted B, due to degeneracies in the model with other parameters. The magnetic field in the leptonic PWN model can be independently constrained by considering the size of the candidate PWN and the cooling time of the leptons diffusing throughout its volume. Equation (6) in S. P. Reynolds et al. (2018) gives the synchrotron-loss timescale $t_{1/2}$, the time for an electron primarily emitting at energy E to lose half its energy in a magnetic field B:

$$t_{1/2} = 1.2 \text{ kyr} \left( \frac{E}{\text{keV}} \right)^{-1/2} \left( \frac{B}{10 \text{ µG}} \right)^{-3/2}$$

Though the great uncertainty in the distance to the candidate PWN prevents accurate measurement of its size, we can parameterize the diffusion timescale of leptons diffusing at velocity $v_{ad} \approx 1000 \text{ km s}^{-1} \approx 1 \text{ pc kyr}^{-1}$ (O. Porth et al. 2016; S. P. Reynolds et al. 2018), using the angular size of the candidate PWN $\theta \approx 2'$ and scaling for some distance d:

$$t_{\text{dif}} = \frac{\theta d}{v_{\text{ad}}} = 0.5 \text{ kyr} \left( \frac{d}{\text{kpc}} \right)$$

Setting these two timescales equal to each other and solving for B gives

$$B = 18 \text{ µG} \left( \frac{d}{\text{kpc}} \right)^{-2/3} \left( \frac{E}{\text{keV}} \right)^{-1/3}$$

This equation presents an upper limit; if B were to exceed the predicted value, leptons would not disperse through the observed angular size of the candidate PWN while cooling. For $d = 1$ kpc and $E = 10$ keV (the XMM spectrum still rising at 8 keV), we obtain an upper limit on B of $B < 8.5$ µG. This qualitative examination would suggest that the leptonic PWN model described above does not require unreasonable magnetic fields for a PWN of this size.

With $\kappa = 1$, no leptonic fit for the entire composite could be obtained, with the gamma-ray emissivity being simply too low to account for more than ≈10% of the observed LHAASO emission. Boosting $\kappa$ even moderately results in IC from the electrons of the PWN being a substantial if not dominant contributor to the observed LHAASO flux. It is also possible that some of the gamma-ray emissivity observed with LHAASO, which is of interest to neutrino observatories like IceCube, is a possible contributor to the galactic neutrino flux (R. Abbasi et al. 2023); neutrino emission is only expected from sources with a hadronic emission component.

In the leptonic model, the total energy in the primary particles is substantially lower than those energies reported for other systems (D. A. Burgess et al. 2022; J. Woo et al. 2023), suggesting that this candidate PWN is a substantially less luminous system than these other objects. In the leptonic
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Reproducing Figure 10 from M. J. Reid et al. (2019), a top-down view of the arms and star-forming regions in the Milky Way galaxy. Given the direction of *l* = 147.0 toward our primary source (green line), it is likely that this candidate PWN is a member of the Local (cyan; *d* < 1 kpc), Perseus (black; *d* < 2 kpc), or perhaps Norma-Outer (red; *d* < 6 kpc) spiral arms.</figcaption></figure>

model, $\log E_{\text{cut}}$/TeV ≈ 2 and 1 µG < *B* < 10 µG—values typical for PWNs that are linked to VHE emission, like the Eel, Dragonfly, or Boomerang.

We only require a single astronomical accommodation to explain most if not all of the observed LHAASO gamma-ray flux (elevated $\kappa$, to increase the ambient photon fields). The candidate PWN could contribute only part of the observed LHAASO gamma-ray emission via IC upscattering, with additional leptonic or hadronic contributions making up the difference. At the moment, we lack the detailed astronomical observations at other wavelengths needed to constrain the many free parameters in a leptonic + hadronic model; further radio, IR, and optical observations will provide further context for our candidate PWN and the broader environment of 1LHAASO J0343+5254u.

## 4.2. Hadronic Models

Though a leptonic model obtains an appropriate fit for this system and others like it, it is also possible that some of the gamma-ray emissivity observed with LHAASO is hadronic in origin, via the decay of neutral pions produced by energetic protons interacting with ambient hadrons. Hadronic emission is of interest to neutrino observatories like IceCube, as a possible contributor to the galactic neutrino flux (R. Abbasi et al. 2023); neutrino emission is only expected from sources with a hadronic emission component.

There are several astronomical unknowns that inhibit a detailed and tightly constrained hadronic model. First, lacking an evolutionary model for a pulsar in the heart of the candidate PWN prevents physical conclusions about the flow of particles into a target molecular cloud, the first step in a hadronic model. Furthermore, there is substantial uncertainty about the molecular cloud that would serve as a target; recent radio observations (discussed below) have only revealed several

tentative cloud regions, with substantial uncertainties in their physical properties and distances. For these reasons, we do not conduct detailed hadronic fitting in this work, leaving more complex models until further MWL data are collected.

## 4.3. The Candidate PWN in the Galactic Environment

For our modeling, we fixed the distance to the primary XMM-Newton source to be *d* = 1 kpc, with its galactic coordinates ($l$ = 147.0, $b$ = -1.8) giving a position on the outskirts of the galactic plane. Supposing that the candidate PWN is of age ≲10 kyr, typical of a >teraelectronvolt PWN (H. E. S. S. Collaboration et al. 2018), it is likely that this system is part of one of the well-mapped spiral arms opposite the galactic center—the homes of star formation in the more rural outer regions of the Milky Way. For comparison, we draw on the work by M. J. Reid et al. (2019), who mapped in detail the star formation regions of the outer arms of our side of the Milky Way, a proxy map of environments where the production of a PWN is likely.

Reproducing Figure 10 from M. J. Reid et al. (2019), and including the vector along *l* = 147.0 for our primary source, we can show in Figure 6 that the distance to our primary source is likely ≲2 kpc, given its likely membership of the Local Arm or Perseus arm, though a distance of ≈4 kpc would be possible if the system is a member of the Norma-Outer arm. 1LHAASO J0343+5254u is an attractive target for further observations, as it is likely closer to us than many other galactic UHE sources and occupies a less busy region of the galactic plane.

## 4.4. Radio Context

We have been conducting observations of radio line emissions using the Nobeyama Radio Observatory (NRO) 45 m radio telescope, to search for associated molecular clouds in the 1LHAASO J0343+5254u region (N. Tsuji et al. 2025). Five molecular clouds are detected within the gamma-ray regions of the LHAASO sources; four have derived distances and $n_{\text{H}}$ in the range 0.3–1 kpc and 200–800 cm⁻³, with the fifth having *d* = 4.0 ± 1.0 kpc and $n_{\text{H}}$ = 2 × 10³ cm⁻³ (related to a nearby optical/IR asymptotic giant star). These radio sources may be molecular clouds, where outflows of energetic protons from the PWN interact with ambient protons to eventually produce gamma rays via neutral pion decay; however, none of them have substantial X-ray flux in our XMM observations, shown in Figure 7. There is also no apparent radio counterpart at the position of the primary X-ray source, nor any extended radio emission pointing to any of the other high-energy sources in this region and relating a kicked pulsar within the candidate PWN back to a radio relic at its birth position.

The upper limits for radio continuum flux at the primary X-ray position are 1.2 × 10⁻¹⁵ erg s⁻¹ cm⁻² at 408 MHz and 4.3 × 10⁻¹⁵ erg s⁻¹ cm⁻² at 1420 MHz. These limits are not particularly constraining, exceeding the fluxes predicted by the leptonic model by several orders of magnitude.

The X-ray emission of our candidate PWN could feasibly be located at the same distance of the detected molecular clouds. In particular, a partial shell-like structure, which was mentioned in Z. Cao et al. (2021), is visible to the east of the primary X-ray source, and the X-ray source seems to be located in the vicinity of the half-shell-like molecular cloud. If they are related to each other, the distance to the primary X-ray source would be ≲0.3 kpc, within the local spiral arm.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The XMM-Newton image at 0.2–12 keV, overlaid with contours of $^{12}$CO (J = 1–0) line emissions within $-40 \text{ km s}^{-1} < V_{\text{LSR}} < 10 \text{ km s}^{-1}$, obtained by the NRO.</figcaption></figure>

**5. Conclusions and Next Steps**

We have reported the discovery of a new candidate PWN in the outer regions of the galactic plane, which may be linked to 1LHAASO J0343+5254u. With its overall featureless X-ray spectrum ($\Gamma_X = 1.9 \pm 0.1$), extended morphology, and spatially varying X-ray hardness, the source is similar to other PWNs related to VHE and UHE emission, like the Eel or Boomerang nebulas. We do not detect a radio point source or substantial X-ray pulsations at this primary source, but further investigations may detect these concrete signs of PWN classification.

We construct an SED linking the primary XMM-Newton source to observed LHAASO emission, using the combined X-ray and VHE emission, plus upper limits established by Swift-BAT, Fermi-LAT, and VERITAS nondetections. We establish a leptonic (synchrotron + IC) model for this combined emission, finding that the leptonic model can explain the MWL emission by invoking elevated ambient IR fields for IC emission. If the IR fields in the environment around our PWN are not high enough to produce all the observed LHAASO emission, then the difference may be made up of gamma rays from hadronic processes in nearby molecular clouds.

In this leptonic model, our primary source may be a VHE PWN, similar to other systems, like the Eel (D. A. Burgess et al. 2022), the Boomerang (I. Pope et al. 2024), HESS J1825-137 (H. E. S. S. Collaboration 2019), or the Dragonfly (J. Woo et al. 2023).

**5.1. Next Steps**

Additional X-ray observations above 10 keV could much more conclusively constrain the magnetic field of the PWN and the total synchrotron emissivity in the PWN. Measuring the X-ray cutoff of the power-law spectrum of the candidate PWN will allow for direct constraints on these physical parameters, so we will propose observations with the Nuclear Spectroscopic Telescope Array (F. A. Harrison et al. 2013) to directly test the model proposed herein and much more tightly constrain *B* in the candidate PWN.

Particularly high-resolution X-ray observations with Chandra could resolve intricate details of the candidate PWN, including the X-ray structure in the central 1′.8 or an X-ray point source in the central region. The detection of a specific X-ray point source would allow for substantially more focused follow-up for a presumptive pulsar, which is not resolved in our XMM-Newton observations.

Finally, radio searches in the region of the primary XMM-Newton source could identify the position of the pulsar itself, allowing for much more detailed physical modeling of an evolving PWN system in an MWL context. In other PWN systems, the radio pulsar is typically located within the X-ray region, suggesting that a radio search for a pulsar within our candidate PWN could be limited to the region immediately around the candidate PWN. The discovery of a radio pulsar would cement the PWN classification and allow for much more detailed evolutionary modeling of the X-ray emission region. Finally, the 2021 spectrum for LHAASO J0343+5254u may be updated with more recent LHAASO data and/or detections with other gamma-ray observatories.

While the candidate PWN is a feasible counterpart to the LHAASO emission, the region around 1LHAASO J0343+5254u may have a nuanced spatial geometry; the two KM2A sources identified in Z. Cao et al. (2024) in particular may mean that the VHE emission in this region is possibly two independent confused or unresolved systems. Continuing analysis of the nearby Fermi-LAT gamma-ray source and radio observations of molecular clouds in the region will constrain the MWL landscape near 1LHAASO J0343+5254u. Detailed in-progress examinations with other gamma-ray telescopes like HAWC and VERITAS will also constrain the environment around 1LHAASO J0343+5254u.

**Acknowledgments**

This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. We are grateful to Alison Mitchell for substantial guidance and comments on this work. Support for K.M. and J.W. was provided by NASA through the XMMC22 grant.

*Software:* FTools (Heasarc 2014), Xspec (K. A. Arnaud 1996), DS9 (W. A. Joye & E. Mandel 2003), FermiPy (M. Wood et al. 2017), XMMSAS (C. Gabriel et al. 2004).

**ORCID iDs**

Stephen DiKerby https://orcid.org/0000-0003-2633-2196
Shuo Zhang https://orcid.org/0000-0002-2967-790X
Tülün Ergin https://orcid.org/0000-0003-2423-4656
Naomi Tsuji https://orcid.org/0000-0001-7209-9204
Kaya Mori https://orcid.org/0000-0002-9709-5389
Fabio Acero https://orcid.org/0000-0002-6606-2816
Samar Safi-Harb https://orcid.org/0000-0001-6189-7665
Shunya Takekawa https://orcid.org/0000-0001-8147-6817
Jooyun Woo https://orcid.org/0009-0001-6471-1405

**References**

Abbasi, R., Ackermann, M., Adams, J., et al. 2023, *Sci*, **380**, 1338

Abdollahi, S., Acero, F., Baldini, L., et al. 2022, *ApJS*, **260**, 53

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17
Ballet, J., Bruel, P., Burnett, T. H., Lott, B. & The Fermi-LAT collaboration 2023, arXiv:2307.12546

Bangale, P., & Wang, X. 2023, arXiv:2308.15643

Burgess, D. A., Mori, K., Gelfand, J. D., et al. 2022, ApJ, 930, 148

Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, SSRv, 120, 165

Cao, Z., Aharonian, F., An, Q., et al. 2021, ApJL, 917, L4

Cao, Z., Aharonian, F., An, Q., et al. 2024, ApJS, 271, 25

Cao, Z., della Volpe, D., Liu, S., et al. 2019, arXiv:1905.02773

De Sarkar, A., & Majumdar, P. 2024, A&A, 681, A34

Gabriel, C., Denby, M., & Fyfe, D. J. 2004, in ASP Conf. Ser. 314, Astronomical Data Analysis Software and Systems (ADASS) XIII, ed. F. Ochsenbein, M. G. Allen, & D. Egret (San Francisco, CA: ASP), 759

Goddard Space Flight Center 2024, Fermi-LAT Background Models, https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html

H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018, A&A, 612, A2

H. E. S. S. Collaboration, Abdalla, H., Aharonian, F., et al. 2019, A&A, 621, A116

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103

Heasarc 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, ascl:1408.004

HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116

Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1

Joye, W. A., & Mandel, E. 2003, in ASP Conf. Ser. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook (San Francisco, CA: ASP), 489

Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014, ApJ, 783, 100

Oh, K., Koss, M., Markwardt, C. B., et al. 2018, ApJS, 235, 4

Pope, I., Mori, K., Abdelmaguid, M., et al. 2024, ApJ, 960, 75

Porth, O., Vorster, M. J., Lyutikov, M., & Engelbrecht, N. E. 2016, MNRAS, 460, 4135

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131

Reynolds, S. P., Borkowski, K. J., & Gwynne, P. H. 2018, ApJ, 856, 133

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191

Tsuji, N., Takekawa, S., Mori, K., et al. 2025, ApJ, 983, 22

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

Voges, W., Aschenbach, B., Boller, T., et al. 2000, IAUC, 7432, 3

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

Woo, J., An, H., Gelfand, J. D., et al. 2023, ApJ, 954, 9

Wood, M., Caputo, R., Charles, E., et al. 2017, ICRC (Bexco), 35, 824

Zabalza, V. 2015, ICRC (The Hague), 34, 922