X-Ray Characterization of the Pulsar PSR J1849−0001 and Its Wind Nebula G32.64+0.53 Associated with TeV Sources Detected by H.E.S.S., HAWC, Tibet ASγ, and LHAASO

We analyzed the X-ray observation data of J1849 obtained by Chandra, XMM-Newton, NICER, Swift, and NuSTAR. We reprocessed the Chandra ACIS-S data acquired on 2012 November 16 for 23 ks (ObsID 13291) using the chandra_repro tool of CIAO 4.14. The XMM-Newton observation (2011 March 23; ObsID 0651930201) was carried out with the full-frame and small-window mode for the MOS and PN detectors, respectively. We processed these data with the emproc and epproc tools of SAS 20211130_0941, and we further removed particle flares following the standard procedure. Net XMM-Newton exposures after the data screening are 47, 50, and 37 ks for MOS1, MOS2, and PN, respectively. The NICER data were collected between 2018 February 13 and 2022 July 14 (143 observations with ObsIDs 1020660101–5505050501). We processed the NICER observations using the nicerl2 script integrated in HEASOFT v6.31. We also used the Swift windowed-timing mode data taken on 2017 March 19 for 12.5 ks (ObsID 00034978002) to extend the baseline for our pulsar timing study (Section 2.2). The Swift data were processed with the xrtpipeline tool. The source was observed by NuSTAR on 2020 November 24, and we processed the data using nupipeline along with the strict SAA filter. The net exposure of the NuSTAR observation is 51 ks for each of FPMA and FPMB.

A timing solution for J1849, valid between MJD 57830.9 and 58391.0, was reported by Bogdanov et al. (2019; see also Ho et al. 2022). We extend this solution to cover the 5.4 yr period spanning the epochs of the Swift, NICER, and NuSTAR observations. In the Swift data, we extracted events within an R = 16'' circle centered at the pulsar in the 1–6 keV band, while for the NuSTAR data, we used an R = 60'' extraction circle in the 3–60 keV band. Because the source pulsations were not well detected by NICER below 2 keV, we employed the 2–8 keV band for the NICER analysis. Our analysis included 108 NICER observations (for both timing and spectral analyses) in which the pulse profiles were confidently measured. We applied barycenter corrections to the photon arrival times, utilizing the source position (R.A., decl.) = (2822568023, −00216153), and carried out a semi-phase-coherent timing analysis (e.g., An & Archibald 2019).

We constructed a pulse-profile template using the Swift and NICER data for which the existing solution is valid (MJD 57830.9–58391.0). We then fit the template with a double-Gaussian function and used the function in our analysis. We progressively folded the later observations and measured a phase shift of each observation by fitting the observed pulse profile with the function. In the case that notable drifts occurred in the pulse-arrival phases, we updated the timing solution by introducing a higher-frequency derivative to ensure phase alignment. Following this, we created a new pulse-profile template using the revised timing solution and updated the Gaussian function. We repeated this procedure for all Swift, NICER, and NuSTAR observations. To ensure phase coherence across the Swift, NuSTAR, and 4.5 yr NICER data, it was necessary to incorporate four frequency derivatives. The results of this semi-phase-coherent analysis are presented in Table 1 and Figure 1. The source pulsations were significantly detected up to the 50–60 keV band by NuSTAR with a chance probability of 6 × 10⁻⁴.

Zoom In Zoom Out Reset image size

Because our timing solution is not valid at the much earlier XMM-Newton observation epoch due to the uncertainties in the timing parameters, we performed H tests (de Jager et al. 1989) to measure the XMM-Newton pulse profile. For this analysis, we used the 1–10 keV band and an R = 16'' circular region to extract source events from the PN data. We conducted a search for pulsations around the previously reported period while holding $\dot{P}$ fixed at the measured value of 1.4224 × 10⁻¹⁴ (Gotthelf et al. 2011; Kuiper & Hermsen 2015). The pulsations were detected with high significance with an H statistic of ∼4000. The XMM-Newton–measured profile is displayed in Figure 1. The on-pulse profiles (off-pulse emission was subtracted using the phase interval ϕ = 0.75–0.95) measured with the three observatories agree with each other very well, and we further verified that the on+off profiles measured with XMM-Newton and NuSTAR, obtained by subtracting the background taken from circular regions (R = 32'' for XMM-Newton and R = 60'' for NuSTAR) at R ∼ 150'' from the pulsar, also agreed well with each other. It is important to note that the pulsations were significantly detected in the <2 keV band of the PN data. This means that the nondetection of <2 keV pulsations in the NICER data likely results from elevated background levels (e.g., flares) present in the data, rather than originating from intrinsic emission properties of J1849 or the influence of strong Galactic absorption.

While the extended PWN G32.6 has been previously identified (Terrier et al. 2008; Gotthelf et al. 2011; Kuiper & Hermsen 2015; Vleeschower Calas et al. 2018), a careful image analysis of the Chandra, XMM-Newton, and NuSTAR data is required to assess contamination by the pulsar in the PWN emission region. This is crucial for analyzing the low imaging resolution data collected by XMM-Newton and NuSTAR. The optics vignetting effect is not a concern for our investigation of the central region (e.g., <30''). At R = 150'', this effect was estimated to be ∼3%–8% for XMM-Newton and Chandra and ∼14% for NuSTAR.

2.3.1. Analysis of the Chandra and XMM-Newton Images

We present the 2–7 keV count image measured by Chandra in panels (a) and (b) of Figure 2. The central R ≲ 1'' region is significantly affected by the pileup effect caused by the bright pulsar. The R ≲ 40'' region surrounding the pulsar appears brighter than other regions and displays some possible substructures. These include a jetlike feature in the northeast direction, a counterjet in the southwest direction (distinct from the readout trail), and an arc to the south of the pulsar. However, each of these structures is statistically insignificant (≤2.5σ).

Zoom In Zoom Out Reset image size

On a larger scale, the extended PWN with a radius of R ≈ 150'' (Figure 2(c)) is detected by XMM-Newton, as previously reported (e.g., Terrier et al. 2008). This X-ray PWN overlaps with the regions of TeV emission measured by H.E.S.S and LHAASO (H.E.S.S. Collaboration et al. 2018; Cao et al. 2023). HAWC has also identified a source with a 016 offset in the northeast direction (Albert et al. 2020).

We further characterized the PWN emission using radial profiles around the pulsar position. The Chandra profile (Figure 3(a)) exhibits J1849, a narrow core, and a broad wing. The latter two components extend to large distances, representing the extended PWN. The profile also shows spiky features that arise from contamination by point sources within the field. To model the Chandra profile, we considered a uniform background (exposure map) and a PWN model while excluding the piled-up pulsar (R ≤ 6'') and other contaminating sources. However, we found that a single-component model for the PWN emission failed to adequately fit the data. Consequently, we adopted a two-component function, a Gaussian combined with an exponential tail: $A\exp (-{r}^{2}/2{\sigma }^{2})+B\exp (-r/l)$. This function successfully fit the data, as illustrated in Figure 3(a), with the best-fit parameters of σ = 10.8 ± 2.2'' and l = 49.5 ± 5.8''.

Zoom In Zoom Out Reset image size

We also carried out an analysis of the 1D radial profile derived from the XMM-Newton/MOS data (Figure 3(b)). For the analysis of the profile, we utilized the MOS point-spread functions (PSFs) generated using the psfgen tool of SAS to represent the pulsar emission. To construct a PWN model, we convolved the Chandra PWN model with the XMM-Newton PSFs. In this analysis, we adopted a flat background, i.e., the exposure map, and fit the profile by optimizing the normalization factors of the PSF, PWN, and background components. In the fitting procedure, we excluded regions that were contaminated by point sources, such as the small bump at R ≈ 60''. The results are shown in Figure 3(b). Notably, within a radius of R ≤ 50'', the pulsar appears brighter than the PWN, and in the central 16'' region, the estimated 2–10 keV PWN counts amount to <2% of the pulsar's count.

The XMM-Newton/PN small-window image, despite its limited coverage of the PWN, allowed the distinction between the on- and off-pulse emission from the pulsar. We generated radial profiles of the PN data utilizing the on- and off-pulse intervals defined in Figure 1 and then fit these profiles with the same model applied to the MOS data. Our analysis revealed that the 2–7 keV off-pulse emission from the pulsar (i.e., PSF) is 22% ± 1% of its on−off emission (per spin cycle).

2.3.2. Analysis of the NuSTAR Image

As previously mentioned, the prominent pulsar emission dominates over the PWN emission in the XMM-Newton and NuSTAR data. This section aims to precisely characterize the pulsar's spectrum by using the NICER, XMM-Newton, and NuSTAR data (Section 2.4.1). This is crucial to the NuSTAR measurement of the PWN spectrum because contamination of the PWN region by the pulsar's emission is a concern in the NuSTAR analysis. While the strong contamination from the on-pulse emission of the pulsar can be minimized by using the off-pulse interval (Figure 1) in the measurement of the PWN spectrum, there is still some contamination from the pulsar's off-pulse emission. We characterize this off-pulse spectrum using the XMM-Newton data by comparing the "on+off" spectrum, in which contamination from the PWN is negligible, with the "on−off" spectrum, as shown below (Section 2.4.2). Additionally, we assess contamination by other point sources within the PWN (Section 2.4.3). Subsequently, we take into account emissions from the pulsar and other point sources when measuring the PWN spectra (Sections 2.4.3 and 2.4.4).

We performed spectral fitting in XSPEC v12.13.0 c by employing the χ² statistic. To account for Galactic absorption, we adopted the tbabs model with the wilm abundances (Wilms et al. 2000) and the vern cross section (Verner et al. 1996).

2.4.1. On−Off Spectrum of the Pulsar J1849

To measure the on−off spectrum of J1849, we selected events within the on- and off-pulse intervals (Figure 1) from the XMM-Newton/PN, NICER, and NuSTAR data. For extracting on-pulse (source) spectra, circular regions with radii of R = 16'' and 60'' were employed for XMM-Newton and NuSTAR, respectively. The off-pulse (background) spectra were extracted within the same regions used for the source regions. As for the NICER data, source and background spectra were extracted from each of the 108 observations. We then combined these NICER spectra using the addascaspec script to construct high-quality source and background spectra. Corresponding spectral response files were generated following the standard procedure suitable for each observatory. We then grouped the source spectra to ensure at least 20 counts per spectral bin.

We jointly fit the on−off spectra measured by XMM-Newton, NICER, and NuSTAR with an absorbed PL or the log parabolic PL (logpar; Massaro et al. 2004) model. For each spectrum, a cross-normalization factor was incorporated, with the factor for the NICER spectrum set to 1. The energy bands used were 0.3–10, 2–8, and 5–60 keV for the XMM-Newton, NICER, and NuSTAR spectra, respectively. The selection of these energy bands aimed to minimize contamination from background and (cross-)calibration uncertainties (e.g., Madsen et al. 2022).

While the PL fit was acceptable, adopting the logpar model resulted in an improved fit (Table 2 and Figure 4). The latter model was favored with an F-test probability of 5 × 10⁻⁵. The cross-normalization factors for the NuSTAR spectra were consistent with 1 at ≲1σ levels. However, the factor for the XMM-Newton/PN spectrum was measured to be 0.94 ± 0.02, exhibiting a deviation from 1 at a ∼3σ significance level. This discrepancy might indicate some systematic effects, potentially stemming from cross-calibration issues (see Madsen et al. 2017). These cross-normalization factors are taken into account below while measuring the PWN spectrum.

Zoom In Zoom Out Reset image size

Table 2. Spectral Analysis Results

Data	Instrument a	Model b	Energy	NH	Γ/a	b	F2−10 keV c	χ2/dof	Comment
			(keV)	(1022 cm−2)
PSR	XP+Ni+Nu	logpar	0.3–60	6.4 ± 0.4	1.40 ± 0.03	0.38 ± 0.09	4.09 ± 0.07 d	1325/1296	On−off
PSR	XP+Ni+Nu	PL	0.3–60	8.1 ± 0.2	1.42 ± 0.03	⋯	4.29 ± 0.07 d	1342/1297	On−off
PSR	XM+XP	PL	0.3–10	6.7 ± 0.2	1.21 ± 0.03	⋯	4.75 ± 0.06	820/811	On+off
PSR	XP	PL	0.3–10	6.7 e	1.20 ± 0.04	⋯	3.85 ± 0.06 d	456/410	On−off
PWN	CXO	PL	0.3–10	6.4 f	1.96 ± 0.33 ± 0.12 g	⋯	1.44 ± 0.18 ± 0.02 g	43/45
PWN	Nu	PL	5–20	6.4 f	2.64 ± 0.41 ± 0.36 g	⋯	1.94 ± 0.68 ± 0.11 g	56/61
PWN	CXO+Nu	PL	0.3–20	6.4 f	2.25 ± 0.29 ± 0.15 g	⋯	1.37 ± 0.15 ± 0.05 g	104/107

Notes.

^a XP: XMM-Newton/PN; XM: XMM-Newton/MOS; Ni: NICER; Nu: NuSTAR; CXO: Chandra. ^b $K{(E/{E}_{p})}^{[-a-b\mathrm{log}(E/{E}_{p})]}$, with E_p = 10 keV for logpar and K(E/1 keV)^−Γ for PL. ^c Absorption-corrected 2–10 keV flux in units of 10⁻¹² erg s⁻¹ cm⁻². ^d Spin-cycle averaged flux. ^e Fixed at the value obtained from the fit of the XM on+off spectrum. ^f Fixed at the value obtained from the logpar model of the XP+Ni+Nu on−off spectra. ^g The second error denotes the systematic uncertainty (see Section 2.4.3).

Download table as:  ASCIITypeset image

Our estimation of N_H = (6.4 ± 0.4) × 10²² cm⁻² toward the source is larger than the previous results of (4.0–4.5) × 10²² cm⁻² (Terrier et al. 2008; Gotthelf et al. 2011; Kuiper & Hermsen 2015; Vleeschower Calas et al. 2018). We believe that the discrepancy could be ascribed to the different abundance model that these authors used. Although they did not report the abundance model, we were able to reproduce their N_H values by using the angr abundances (default in XSPEC).

2.4.2. Estimations of the Off-pulse Spectrum of the Pulsar

2.4.3. Spatially Integrated Spectrum of the PWN G32.6

2.4.4. Spatially Resolved Spectra of the PWN

We performed a spatially resolved spectral analysis of the Chandra data using three annular regions with sizes of R = 5''–20'', 20''–50'', and 50''–150'' to investigate possible spectral variation within the PWN. We extracted ACIS spectra from these regions and fit them along with background spectra constructed from the blank-sky data (Section 2.4.2). We jointly fit the spectra, taking into account the contamination by the pulsar (i.e., logpar model; see Section 2.4.3). Our analysis revealed an increase in Γ with increasing distance from the pulsar (Figure 5); however, this change is statistically insignificant due to large uncertainties.

Zoom In Zoom Out Reset image size

We also used the XMM-Newton/MOS data to investigate the spatial variation of the PWN spectrum. We extracted spectra from annular regions of R = 30''–75'' and 75''–150'', excluding contaminating sources within the regions. For background estimation, we followed the same procedure employed in the Chandra analysis (Section 2.4.3) using the quiescence particle background data generated by the evqpb of SAS instead of the blank-sky data used for Chandra (see Park et al. 2023a, for more detail). We then fit the spectra with the PL+logpar model, where the logpar model accounts for the contamination from the pulsar. As in the Chandra analysis, we introduced a normalization factor to accommodate the off-pulse emission and employed point-source response functions for the logpar model. The PL+logpar model provided a satisfactory fit to the data, and the best-fit parameters for the PL component are presented in Figure 5.

These XMM-Newton measurements are consistent with the Chandra results, albeit yielding a larger Γ (≈2.1) for the R = 75''–150'' region compared to previous analyses of the same region (Γ=1.7–1.8; Terrier et al. 2008; Kuiper & Hermsen 2015; Vleeschower Calas et al. 2018). Although the difference is statistically insignificant, given the substantial uncertainties, we speculate that the contamination from the pulsar emission in the previous analyses may have contributed to the discrepancy. By omitting the pulsar model from our analysis, we were able to reproduce the previous results. It is worth noting that Gotthelf et al. (2011) also considered the contamination by the pulsar and inferred Γ = 2.1 ± 0.3 for the PWN spectrum within an R = 30''–150'' annular region.

We also estimated systematic uncertainties for these spatially resolved spectra as we did in Section 2.4.3. The systematic uncertainties were added to the statistical ones in quadrature, and Figure 5 displays the summed uncertainties.

Since our multizone PWN model is extensively described in Park et al. (2023a), we provide a concise overview hereafter. We assume that the pulsar injects electrons with a PL distribution

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{e}}{d{\gamma }_{e}{dt}}={\dot{N}}_{0}{\gamma }_{e}^{-{p}_{1}}\end{eqnarray} \tag{ 1 }$$

between ${\gamma }_{e,\min }$ and ${\gamma }_{e,\max }$ into the TS at a distance of R_TS from the pulsar such that the particle energy corresponds to a fraction η_e of the pulsar's spin-down power. The remaining power is contained in B (fraction of η_B ) or radiated as gamma-ray emission from the pulsar (fraction of η_γ ). As the magnetic energy is thought to be almost fully dissipated at the shock (Kennel & Coroniti 1984b), we assume η_B ≪ 1 in the PWN. Additionally, the pulsar (J1849) was not detected by the Fermi Large Area Telescope (LAT; Atwood et al. 2009), indicating η_γ ≪ 1. Therefore, in the following study, we assume η_e ≈ 1. These parameters affect the energy of the electrons injected by the pulsar into the PWN, and thus the predicted flux of the PWN could vary depending on the assumed η_e . A smaller value of η_e (i.e., larger η_γ ) would cause the model to underpredict the measured fluxes (but see Section 3.3).

The electrons propagate within the PWN via advection and diffusion, and they lose energy due to radiation (synchrotron and IC) and adiabatic expansion. Regarding the bulk flow (V_flow) of the electrons and B within the PWN, we adopt a PL dependence on the distance (r) between the pulsar and an emission zone:

$$\begin{eqnarray}&&{V}_{\mathrm{flow}}(r)={V}_{0}{\left(\displaystyle \frac{r}{{R}_{\mathrm{TS}}}\right)}^{{\alpha }_{V}}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&B(r)={B}_{0}{\left(\displaystyle \frac{r}{{R}_{\mathrm{TS}}}\right)}^{{\alpha }_{B}},\end{eqnarray} \tag{ 3 }$$

where we assume α_B + α_V = −1, which holds if the B structure is toroidal and the magnetic flux is conserved (e.g., Reynolds 2009). We assume the diffusion coefficient to be (e.g., Abeysekara et al. 2017)

$$\begin{eqnarray}&&D={D}_{0}{\left(\displaystyle \frac{B}{100\,\mu {\rm{G}}}\right)}^{-1}{\left(\displaystyle \frac{{\gamma }_{e}}{{10}^{9}}\right)}^{1/3}.\end{eqnarray} \tag{ 4 }$$

We let the injected electrons propagate over an assumed age of the PWN, and then we project their radiation onto the tangent plane of the observer to compute both spatially resolved and integrated SEDs.

While the model employs multiple parameters, some of them can be estimated first based on the SED measurements and one-zone model calculations. These estimations serve as inputs for the detailed modeling (Section 3.3) and help with obtaining converged solutions more quickly.

The similar values of ν F_ν for the X-ray (synchrotron) and TeV (gamma-ray) emission indicate that the magnetic energy density (${U}_{B}\equiv \tfrac{{B}^{2}}{8\pi }$) is comparable to the energy density of the cosmic microwave background (CMB; u_CMB) seed photons for IC upscattering. From this, we estimated B ∼ 3 μG. Other complications, such as the contribution of the Galactic IR field and spatially varying B, are considered by the multizone model. Additionally, the measured TeV–to–X-ray flux ratio can provide information on the true age of the PWN. According to the correlations reported by Zhu et al. (2018), the flux ratio of ∼1 for G32.6 suggests that its true age is 4–9 kyr, much smaller than the spin-down age (τ_c = 43 kyr) of J1849. We assume a t_age of 9 kyr in this work, but models with different values of t_age (and other covarying parameters, e.g., D₀) can also fit the data well (see Park et al. 2023a, for parameter covariance).

Given the estimated B strength (B ≈ 3 μG), the synchrotron emission at ∼20 keV detected by NuSTAR suggests a ${\gamma }_{e,\max }$ of ≥6 × 10⁸. The X-ray photon index of Γ ≈ 1.5 in the innermost region corresponds to the PL index of the uncooled particle distribution of p₁ = 2Γ − 1 ≈ 2. These estimations of ${\gamma }_{e,\max }$ and p₁ allow for the determination of ${\gamma }_{e,\min }$, since η_e ≈ 1.

Although the TS of G32.6 has not been resolved, likely due to its compact size obscured by the bright pulsar emission, we assumed R_TS = 0.1 pc, as has been measured for several PNWe (e.g., Ng & Romani 2004; Kargaltsev & Pavlov 2008). This assumption is also in accord with the correlation reported by Bamba et al. (2010). Note that the exact value of R_TS does not significantly impact the model outputs as long as it is sufficiently smaller than the PWN size. Although the X-ray PWN does not exhibit a sharp boundary, we assumed an X-ray PWN size of R_PWN = 150''; this radial position is where we halt the advection flow and lower B to a slightly smaller value (2.3 μG in the PWN to 2 μG outside). When r > 150'' (outside the X-ray PWN region), we do not invoke α_B + α_V = −1; both B and V_flow are assumed to be constant. The large region for the TeV emission (R = 009) corresponds to 11 pc for an assumed distance of 7 kpc (H.E.S.S. Collaboration et al. 2018) and implies that the TeV-emitting electrons (γ_e ≈ 10⁷) can go beyond the X-ray PWN due to diffusion. Considering the diffusion only, $2\sqrt{{{Dt}}_{\mathrm{age}}}\gtrsim 10$ pc gives an estimate of D₀ ≳ 1 × 10²⁶ cm² s⁻¹ for the assumed age of t_age = 9 kyr (e.g., see Tang & Chevalier 2012; Porth et al. 2016; Guest et al. 2019, for values of diffusion coefficients for various PWNe).

We used the aforementioned estimations of the parameters (Section 3.2) as initial values for our multizone model and further adjusted them to fit the measured emission properties of G32.6. For the seeds of the IC emission, we adopted the interstellar radiation field (ISRF) estimated by GALPROP (Porter et al. 2022), along with the CMB radiation. The results are presented in Figure 5 and Table 3. Our model broadly matches the broadband SED, as well as the radial profiles of brightness and Γ.

Table 3. Parameters for the SED Model in Figure 5

Parameter	Symbol	Value
Spin-down power	LSD	9.8 × 1036 erg s−1
Characteristic age of the pulsar	τc	43 kyr
Age of the PWN	tage	9000 yr
Size of the PWN	Rpwn	5 pc
Radius of the TS	RTS	0.1 pc
Distance to the PWN	d	7.0 kpc
Index for the particle distribution	p1	2.12
Minimum Lorentz factor	${\gamma }_{e,\min }$	10
Maximum Lorentz factor	${\gamma }_{e,\max }$	108.9
Magnetic field	B0	7.5 μG
Magnetic index	αB	−0.3
Flow speed	V0	0.04c
Speed index	αV	−0.7
Diffusion coefficient	D0	1.0 × 1026 cm2 s−1
Energy fraction injected into particles	ηe	0.95
Energy fraction injected into B field	ηB	0.004

Download table as:  ASCIITypeset image

The amplitude of the model-predicted SED is determined by a complex interplay between the spectrum of the injected electrons (e.g., η_e and p₁), the time the electrons spend in the emission region (e.g., t_age and V_flow), and the environmental factors (B and external IR seed density u_IR). In this work, we assumed a large value of η_e , indicating weak gamma-ray emission from the pulsar, since J1849 was not detected by the LAT. We verified that different values of η_e could be easily accommodated in our model by adjusting other parameters, such as t_age, B, u_IR, and/or p₁ (see Park et al. 2023a, for the parameter covariance).

An electron injection spectrum with p₁ ≈ 2 was necessary to match the X-ray SED as well as the <10 TeV data. For this p₁ value, a Γ = 1.5 PL emission from uncooled electrons is expected. These uncooled electrons are responsible for the Γ ≈ 1.5 X-ray spectrum in the inner R ≲ 20'' region. On the other hand, the Γ ≈ 2 spectra in the outer regions are produced by cooled electrons. This implies that the cooling break (i.e., the peak of the curved SED) should be visible in the spatially integrated X-ray SED. The computed SED model predicts a peak at 5 keV. Hence, the model predicts that the <5 keV SED exhibits a harder spectrum than above the energy, as indicated by the separate fits of the Chandra and NuSTAR data; however, the broadband SED fits well with a single PL, probably because of the large measurement uncertainties.

In this model, TeV emission is primarily produced by IC scattering of electrons with γ_e ≈ 10⁷–10⁸ off of the CMB photons. Given the inferred B value, the assumed t_age of 9 kyr corresponds to the synchrotron+IC cooling timescale of γ_e ∼ 8 × 10⁷ electrons. Consequently, a spectral break of the IC SED is expected at ∼40 TeV. Moreover, the Klein–Nishina effect becomes significant for ≥20 TeV (IR seeds) and ≥100 TeV (CMB seeds) gamma rays. These factors explain the steep LHAASO SED at ≳10 TeV.

The radial profiles of brightness and Γ provide constraints on α_B (and thus α_V = −1 − α_B ) and D₀ values. We used a slowly varying B, but different sets of α_B and D₀ can also reproduce the profiles. For the measured radial profiles and the assumed PL trend for V_flow (Equation (2)), our model constrains α_V to be between −0.8 and −0.6. However, it is important to emphasize that the choice of a PL for V_flow was made on a phenomenological basis and lacks a physical motivation and that the actual functional form of the flow in PWNe is not well known. Our model has the flexibility to accommodate various trends; e.g., the radial profiles of V_flow and B predicted for ideal magnetohydrodynamic flow by Kennel & Coroniti (1984a; see also Reynolds 2003) can also fit the data when we freely adjust the other parameters (e.g., V₀, B₀, and D₀). These parameter values cannot be determined uniquely because of the covariance between them. The degeneracy can be broken to a certain extent if we can constrain some of the parameters; e.g., measuring the expansion speed of the PWN can help determine V₀ and α_V (Equation (2)), and B₀ can be well constrained if u_CMB + u_IR are known. Nonetheless, the values presented in Table 3, which are similar to the values inferred for other PWNe (e.g., Park et al. 2023b), can account for the measurements.

A curvature in the on−off spectrum of J1849 has been suggested by Terrier et al. (2008) and Kuiper & Hermsen (2015) based on comparisons of a low-energy on−off spectrum and a high-energy on+off spectrum. While these comparisons hinted at the presence of curvature, the definite establishment was hindered as the comparisons involved distinct quantities, on−off versus on+off emissions. Furthermore, issues such as contamination from the PWN in the INTEGRAL data and cross-calibration problems were not addressed in the previous studies. By combining the high-quality X-ray data acquired by XMM-Newton, NICER, and NuSTAR, we were able to conclusively demonstrate evidence of a spectral curvature in the on−off pulse emission of J1849.

The pulsar J1849 is among the four "MeV pulsars" listed in Harding & Kalapotharakos (2017; see Kuiper & Hermsen 2015, for more sources). These pulsars are young and energetic RPPs that exhibit strong emission of pulsed nonthermal X-rays while lacking significant radio or GeV emission. Harding & Kalapotharakos (2017) proposed that accurately characterizing the SEDs of MeV pulsars can aid in understanding the mechanism of pair production in RPPs. According to their study, the frequency of the SED peak (E_SR) of the pulsar emission is predicted to be ∝B_S B_LC, where B_LC represents B at the light cylinder, if the pairs are generated by polar-cap cascades. On the other hand, if the pairs are produced in the outer gap, ${E}_{\mathrm{SR}}\propto {B}_{\mathrm{LC}}^{7/2}$ is expected (Zhang & Cheng 2000; Harding & Kalapotharakos 2017).

Measurements of the SED peaks for PSR B1509−58 (E_SR = 2.6 ± 0.4 MeV; Chen et al. 2016) and PSR J1846−0257 (E_SR = 3.5 ± 1.1 MeV; Kuiper et al. 2018) indicate that the outer-gap scenario is improbable. This is evident, as these two pulsars possess substantially different B_LC values, yet their E_SR values are very similar (Figure 6, top). However, testing the polar-cap scenario with these two pulsars is challenging, since their B_S B_LC values are almost the same (Figure 6, bottom).

Zoom In Zoom Out Reset image size

The NuSTAR data of J1849 enabled us to measure ${E}_{\mathrm{SR}}={58}_{-17}^{+43}$ keV for its SED, thereby confirming the E_SR–B_S B_LC correlation (Figure 6, bottom). The inclusion of the J1849 measurement in the E_SR–${B}_{\mathrm{LC}}^{7/2}$ trend argues against the outer-gap scenario. Measuring E_SR for more MeV pulsars would greatly contribute to a more detailed understanding of the pair-production mechanism.

G32.6 displays a spectrally hard X-ray emission (Γ ≈ 1.5) in the inner region (R < 20'' from the pulsar). We note that Kuiper & Hermsen (2015) reported smaller values of Γ in the inner regions. We speculate that the discrepancy in the Γ measurements, albeit within the uncertainties, is caused by the contamination from the pulsar emission in the analysis of Kuiper & Hermsen (2015). The Γ ≈ 1.5 spectrum indicates that the energy index of the electron distribution is p₁ ≈ 2, although the large uncertainty of ΔΓ ≈ 0.5 in this region poses challenges in accurately constraining p₁. Spatially integrated SEDs can also offer an alternative method for determining the p₁ of the uncooled electron distribution, as it is reflected in the synchrotron SED in the IR-to-optical band or in the IC SED at energies ≲10 TeV. However, measurements in these bands are either lacking or of poor quality, hampering precise estimations of p₁. Moreover, the observed SED is a superposition of multizone emission components from various radial regions projected onto the observer's tangent plane. Thus, the intrinsic spectrum in the inner region might have a lower value of p₁. Nevertheless, if the inferred value of p₁ ≈ 2 is indeed accurate, it may suggest that magnetic reconnection or shock drift acceleration plays a role in particle acceleration at the TS of PWNe (e.g., Summerlin & Baring 2012; Sironi & Cerutti 2017).

Alternatively, the Γ ≈ 1.5 spectrum in the innermost region could be attributed to putative substructures such as the jet and counterjet (Figure 2(b)), as these features can exhibit spectrally hard emission (e.g., Safi-Harb & Kumar 2008; An et al. 2014a). These substructures might potentially contain a distinct electron population. A deeper Chandra observation will facilitate more definitive detections and spectral characterization of these substructures, elucidating the origin of their hard spectra. This will, in turn, provide valuable clues about the mechanisms responsible for particle acceleration at the TS of middle-aged PWNe.

At large distances from the pulsar (e.g., R > 20''), the X-ray spectrum follows a Γ ≈ 2 PL, similar to other PWNe. In comparison to the Γ ≈ 1.5 spectrum in the inner region, this indicates rapid particle cooling. As electrons travel outward within the inner region (R ≈ 20''), they progressively lose kinetic energies. As a result, the emission peak of the highest-energy electrons will appear in the X-ray band and enable determination of the flow speed and B in the PWN. However, the spectral variation over the PWN was not measured with high significance with the current X-ray data, and further confirmation of this spectral softening is necessary with future observations.

We applied a multizone PWN model to infer the properties of the PWN G32.6 (Table 3). While our model broadly reproduces the measurements (Figure 5), the rapid change in Γ in the inner regions, if real, is not very well explained. At face value, the ΔΓ ≈ 0.5 can be explained by synchrotron cooling, which causes a change of the spectral slope of 1 (in ideal cases) between the distributions of the injected and cooled electrons. However, the "observed" spectrum is a superposition of emissions from many different radial regions, and thus the "observed" difference in Γ would be diluted if observed by the current X-ray telescopes and appear smaller than the emitted one (the inner versus outer radial zones). This implies that the difference in the slopes of the electron distributions (relevant to the X-ray band) between R ≲ 20'' (injected) and ≳20'' (cooled) needs to be greater than 1. This is difficult to reproduce within the model unless we employ a BPL distribution for the injected electrons (such BPLs may be required to explain the large amounts of spectral steepening seen in several PWNe; Chevalier 2005).

The parameter values in Table 3 are not well constrained or uniquely determined, as the parameters are interdependent (e.g., see Park et al. 2023a), and the quality of the measurements is rather poor; e.g., various values of η_e can reproduce the observed results, as mentioned in Section 3.3. However, if the energy densities of the seed photons estimated with the GALPROP model, despite its low angular resolution, are accurate, B₀ can be determined relatively well independently of the other parameters based on the shapes and amplitudes of the synchrotron and IC SEDs (Section 3.3). This, in turn, helps estimate ${\gamma }_{e,\max }$ in combination with the NuSTAR spectrum. Moreover, the estimation of B₀ along with our measurements of the brightness and Γ profiles can help pin down α_V and V₀ if the PWN expansion rate is measured.

By fitting the radial profiles of the X-ray brightness and Γ, we inferred the index of the bulk-flow speed to be α_V = −0.7. With this speed profile, electrons can be transported to 8.5 pc, which is sufficient to explain the observed X-ray emission zone with a size of 5 pc. The softening of the X-ray spectrum with increasing distance from the pulsar requires diffusion, as pure advection models predict a constant Γ profile in the inner regions and a rapid increase in the outer regions (Reynolds 2003). In our model, we used the diffusion coefficient of D₀ = 10²⁶ cm² s⁻¹. In this case, electrons responsible for the TeV emission (e.g., off of IR seeds), with γ_e ≈ 10⁷, can move out to a distance of ∼10 pc through diffusion and advection, which can explain the TeV emission size of 0.09°. In this region, located outside the X-ray PWN, the electrons will cool primarily through IC emission and propagate further outward by diffusion. Assuming a spatially uniform ISRF with u_IR = 0.5 eV cm⁻³ and an ISM B (B_ISM) of 1 μG, the cooling timescale of electrons with γ_e ≈ 10⁷ is estimated to be 12 kyr. Then, our model predicts that the TeV-emitting electrons will reach a distance of ∼30 pc from the pulsar in 12 kyr and form an extensive TeV halo based on the assumed diffusion law (Equation (4)). However, this estimate is subject to substantial uncertainty, since the values of u_IR and/or B_ISM may differ significantly.

Our estimate of ${\gamma }_{e,\max }$ corresponds to a maximum electron energy of ${E}_{e,\max }\approx$ 400 TeV, which is comparable to the 740 TeV inferred from modeling the TeV emission of G32.6 (Amenomori et al. 2023). With the lack of measurements at >20 keV, our estimation of ${E}_{e,\max }$ should be regarded as a lower limit. Additionally, if the Γ ≈ 1.5 spectrum in the innermost region represents the intrinsic spectrum of the PWN rather than putative substructures, our model predicts a spectral curvature at ∼5 keV.