[ Article ]
Journal of the Korean Astronomical Society - Vol. 52, No. 5, pp.173-180
ISSN: 1225-4614 (Print) 2288-890X (Online)
Print publication date 31 Oct 2019
Received 10 Apr 2019 Accepted 25 Jul 2019

# INVESTIGATING THE PULSAR WIND NEBULA 3C 58 USING EMISSION MODELS

SEUNGJONG KIM ; JAEGEUN PARK ; HONGJUN AN
Department of Astronomy and Space Science, Chungbuk National University, Cheongju 28644 brelkly@gmail.comgeunjaep@gmail.comhjan@chungbuk.ac.kr

Correspondence to: H. An

Published under Creative Commons license CC BY-SA 4.0

## Abstract

We present IR flux density measurements, models of the broadband SED, and results of SED modeling for the Pulsar Wind Nebula (PWN) 3C 58. We find that the Herschel flux density seems to be slightly lower than suggested by interpolation of previous measurements in nearby wavebands, implying that there may be multiple electron populations in 3C 58. We model the SED using a simple stationary one-zone and a more realistic time-evolving multi-zone scenario. The latter includes variations of flow properties in the PWN (injected energy, magnetic field, and bulk speed), radiative energy losses, adiabatic expansion, and diffusion, similar to previous PWN models. From the modeling, we find that a PWN age of 2900–5400 yrs is preferred and that there may be excess emission at ~ 1011 Hz. The latter may imply multiple populations of electrons in the PWN.

## Keywords:

pulsars: general, ISM: individual objects: 3C 58, radiation mechanisms: non-thermal

## 1. INTRODUCTION

High-energy emission in astrophysical objects is produced by non-thermal processes and can help us to understand energetic environments in the Universe. In many cases, the processes responsible for high-energy radiation are synchrotron and inverse-Compton emission which can be studied in blazars, supernova remnants (SNRs), gamma-ray binaries, and pulsar wind nebulae (PWNe). It is believed that particles are accelerated to extremely high energies in these objects (e.g., Drury 1983; Sironi et al. 2015; Kang 2010), resulting in broadband continuum emssion from radio to TeV energies. In particular, PWNe are very important to understand relativistic shocks, magnetohydrodynamic (MHD) flows (Sironi et al. 2015), and high-energy leptons in our Galaxy (e.g., the multi-GeV positron excess; Yüksel et al. 2009).

A PWN is formed near the center of a supernova remnant as the pulsar wind flows outwards (Kennel & Coroniti 1984) and interacts with the ambient medium. Leptons in the pulsar wind are bulk-accelerated beyond the light cylinder (Aharonian et al. 2012) and interact with the surrounding medium to form a termination shock (Kennel & Coroniti 1984). The cold energetic particles are further accelerated in the shock (Drury 1983). The shock-accelerated particles propagate outwards, interact with the magnetic field and the ambient photon fields, and emit radiation across the whole electromagnetic spectrum via synchrotron and inverse-Compton processes, resulting in a characteristic broadband spectral energy distribution (SED). In addition, magnetic reconnection is also possible in PWNe, which may leave distinct features in the electron and photon distributions (Lyutikov et al. 2018).

Broadband SEDs of young PWNe commonly display a double-hump structure, where a low-energy hump is produced by the synchrotron radiation of shock-accelerated particles, and a high-energy hump by inverse-Compton upscattering of low-energy background radiation by the same particles. Hence, the "direct" low-energy spectrum can be used to infer the shock-accelerated particle distribution which is also imprinted in the high-energy spectrum via interaction with soft-photon fields. Therefore, by accurately measuring and modeling the broadband SEDs, one can study the relativistic shock acceleration processes and plasma flow in PWNe and related systems (e.g., gamma-ray binaries; An & Romani 2017).

Recent MHD studies (e.g., Sironi et al. 2015) found that the particle distribution in PWNe has a high-energy power-law tail with an index of p ≈ 2–3 which is responsible for IR-to-X-ray emission, and that the maximum synchrotron photon energy emitted by the shock-accelerated electrons is in the keV range. However, photon indices of 1.8–2.5 measured in some PWNe imply p = 2.5–4, slightly larger than the theoretical prediction, and the maximum synchrotron photon energy can reach MeV levels (Kargaltsev et al. 2017). In addition, photon indices in PWNe vary with distance from the central pulsar in general, and (cooled synchrotron) X-ray spectra of some PWNe exhibit curvature (e.g., Slane et al. 2004; An et al. 2014; Nynka et al. 2014; Madsen et al. 2015; An 2019). These suggest that the particle acceleration and flow in PWNe are highly complex (e.g., diffusion, magnetic reconnection, and particle generation; Reynolds 2009; Tang & Chevalier 2012; Lyutikov et al. 2018). Therefore, further observational and theoretical investigations of PWNe are required, and in this work we study the typical TeV PWN 3C 58 using IR observations and SED models.

3C 58 is a TeV-emitting PWN at a distance of 3.2 kpc (Roberts et al. 1993) which has been studied at multiple wavelengths (e.g., Green & Scheuer 1992; Slane et al. 2008; Abdo et al. 2009; Aleksić et al. 2014; An 2019). It has a relatively flat radio spectrum which gradually curves down at high radio frequencies (Arnaud et al. 2016). Measurements of the SED in the IR band are difficult due to foreground stars (Green & Scheuer 1992; Slane et al. 2008), but limited data suggest that the photon index in this band is ΓIR ≈ 2. Combined with recent X-ray measurements (e.g., Slane et al. 2002, 2004; An 2019) that find a photon index ΓX ≈ 2.3, the IR-to-X-ray SED implies a spectral break between the IR and the X-ray bands. Intriguingly, Planck and Herschel recently obtained high-frequency radio/IR SEDs (Arnaud et al. 2016) that reveal a relatively sharp break at ~ 1011 Hz which was not accurately modeled previously. With high-quality broadband SEDs available and displaying a typical PWN morphology (torus and jet), 3C 58 is a particularly interesting target to study PWN physics.

While the broadband SED of 3C 58 was modeled previously (e.g., Torres et al. 2013; Tanaka & Takahara 2013; Li et al. 2018), some of its physical properties are still uncertain; the age (between the pulsar's characteristic age of ~800 and ~5000 years estimated from the PWN expansion; Stephenson 1971; Slane et al. 2002) and the average magnetic field strength (~20μG and ~80μG; Aleksić et al. 2014; Green & Scheuer 1992) are still controversial. The maximum particle energy was inferred as ~40–140TeV (Tang & Chevalier 2012; An 2019) but models with higher particle energies could also explain the SED. Hence, it is important to review the previous studies, and compile and supplement observational data to improve the SED characterization. Then an SED model must be applied to the refined SED to characterize the physical properties of the PWN.

Here, we supplement the 3C 58 SED with IR measurements made with Herschel (Section 2). We first confirm the previous 500 μm measurement (Arnaud et al. 2016) and extend it to higher IR frequencies. We then apply emission models to the up-to-date broadband SED (Section 3), and discuss and conclude in Section 4.

## 2. CHARACTERIZING THE IR SED OF 3C 58

Because IR SEDs are direct imprints of accelerated particles in PWNe, it is important to measure the IR SEDs accurately. For 3C 58, the highest-frequency Planck, Herschel/500 μm, and Spitzer data (Arnaud et al. 2016; Slane et al. 2008) can provide some information on the particle distribution, but the uncertainties in the Spitzer measurements are large. In addition, the Herschel measurement may suffer from background systematics which require scrutiny. We thus reanalyze the Herschel data to confirm the previous measurement and supplement the IR SED with data taken at other frequencies.

We use archival 3C 58 observations made with the Herschel spectral and photometric imaging receiver (SPIRE) in three bands: 250 μm, 350 μm, and 500 μm (λ/∆λ ≈ 3). The SPIRE observations were taken in "large map" mode, a scan over a 32'×32' field at a speed of 30'' s–1. We process the data using the standard photometer large map pipeline in the Herschel interactive processing environment (HIPE) version 15.0.1.

### 2.1. Flux Density Measurements

We present the Herschel three-band images in Figure 1. A clear diamond-shape morphology similar to the one observed in radio and X-ray images of 3C 58 (Slane et al. 2004) is present at the center, suggesting that 3C 58 is detected in these bands. We measure the source flux densities using elliptical and surrounding annular apertures for the source and the background, respectively (see Figure 1). We estimate measurement uncertainties following Ngyuen et al. (2010); Ciesla et al. (2012); Valtchanov (2017).

Herschel SPIRE three-band images of 3C 58. The center wavelengths are 500μm, 350μm, and 250μm, respectively, from left to right. White solid-line ellipses indicate the source region, dashed lines denotes background regions.

Because of strong foreground and background emission, accurately measuring the flux density of the diffuse and faint PWN is difficult, and hence the estimate for PWN flux density may vary with the selection of specific background regions. Typically, an annular region surrounding the source is used in Herschel photometry; however, in the case of 3C 58 the emission from the surrounding region is highly anisotropic. In particular, regions at R > 425'' include strong clumpy emission and are unlikely to represent the background affecting the source region; we therefore exclude these regions. We test several background regions near the source (R < 425'') and check how sensitive the flux-density measurements are to the background selection. We select background regions with areas similar to the source region, which leaves three independent annular regions (277'' < R < 329'', 329'' < R < 380'', 380'' < R < 425'').

We measure the source flux density with the above three background selections and find that, depending on the background, the inferred source flux density varies by 5%, 45% , and 88% at 500 μm, 350 μm, and 250 μm, respectively. Note that the 500 μm flux density does not depend strongly on the background and is consistent with previous measurements (Arnaud et al. 2016). The fluctuations are large in the higher-frequency bands, arguably because the source is fainter in those bands. To account for this variation, we take (FmaxFmin)/2 in each band and add these values to the measurement uncertainties in quadrature. The results are presented in Table 1 and in Figure 2.

Flux densities of 3C 58 measured with Herschel

SED data and a one-zone model for 3C 58. Data are from previous works. Radio (red): Green & Scheuer (1992), Planck (cyan): Arnaud et al. (2016), Herschel (orange): this work, IRAS (black): Arnaud et al. (2016), Spitzer (red): Slane et al. (2008), X-ray (brown and green): Slane et al. (2004) and An (2019), GeV (purple): Li et al. (2018), and TeV (blue): Aleksić et al. (2014). Note that there are two spectral breaks in the low-energy hump because we used a power law with two breaks for the electron energy distribution.

### 2.2. IR Slope Constraints

In order to infer the Radio-to-IR SED quantitatively, we fit the SED with a broken power law. The inferred IR photon index is ΓIR = 1.97 ± 0.06 (energy index αIR = 0.97), consistent with the estimate without the Herschel data (ΓIR = 1.93). However, we note that the Herschel points lie slightly below the interpolation of the Planck and Spitzer data (Figure 2), suggesting the presence of a small hump at ~ 1011 Hz already implied by previous IR SED measurements (Arnaud et al. 2016). We further check if the Herschel data require another spectral break (i.e., separate particle population) using broken-power-law fits with and without a ~500 μm break. An ftest comparison of the fits finds that such a break does not improve the fit significantly. Because the Herschel data do not significantly alter the measured SED slope, we ignore the Herschel points and use a single power law across the IR band in the modeling described in the following sections.

## 3. MODELING THE BROADBAND SED OF 3C 58

The shape of an energy distribution of accelerated particles in a PWN is not well known because relativistic shock (termination shock) acceleration is not yet fully understood and magnetic reconnection further complicates the system (e.g., Lyutikov et al. 2018). For the sake of simplicity it is often assumed that the particle distribution is a power law and that the power-law index (p) can be inferred from observations assuming synchrotron emission with νFνναIR+1ν–(p–3)/2, where ν is the observed photon frequency and αIR is the energy index of the photon distribution in the IR band (e.g., Longair 2011). These injected electrons then propagate in the PWN which we assume to have homogeneous (one-zone model) or spatially-varying (multi-zone model) properties assuming spherical symmetry. We then compute an SED of synchrotron and inverse-Compton (IC) emission of the electrons (see, e.g., Finke et al. 2008, for reference), and match the model SED with the observed one. Note that for the IC scattering we use the full Klein-Nishina formula. We assume a distance of 3.2 kpc to 3C 58 (Roberts et al. 1993).

### 3.1. One-Zone Model

We first try to match the observed SED (Figure 2) using a simple one-zone model which assumes that the flow properties are homogeneous over the PWN. This is an oversimplification but the parameter values inferred from this approach can be used as guides for the more realistic multi-zone modeling discussed in Section 3.2. We assume a simple broken power-law electron distribution which accounts for injected (emitting IR to optical photons) and cooled electrons (X-rays). We find that the broken power law cannot explain the radio data; if we extend the soft power-law to low electron energies, the model overpredicts the radio fluxes. Hence, we invoke another break at a high radio frequency and use a doubly-broken power law for the emitting particles. We use the CMB and Galactic IR (with T = 20K with an energy density of 0.8 eV cm–3; Torres et al. 2013) radiation as external seeds for inverse-Compton upscattering, as in previous PWN studies.

In addition, we limit the maximum electron energy to γe,max = 4 × 108, which introduces a sharp cutoff at ~ 1019 Hz. This is to accommodate a potential X-ray break found with NuSTAR (at ~25 keV; An 2019); the significance of that detection is relatively low, and so using larger values of γe,max is possible. This changes the cut-off frequency to a higher value and does not con ict significantly with the current X-ray data. The electron indices are determined from the observed SEDs (see Figure 2): ΓRadio = 1, ΓIR ≈ 2 and ΓX = 2.3 corresponding to p1 = 1, p2 = 3 and p3 = 3.5. We then adjust model parameters (γe,b, B and background density) to match the observed broadband SED. We note that upscattering of the Galactic IR background contributes dominantly (over upscattering of CMB) to the GeV-TeV SED. The model and the parameters are presented in Figure 2 and Table 2, respectively.

Parameters of the one-zone SED model used in Figure 2

The model matches the data well in general, but the GeV and TeV data are difficult to fit simultaneously. The highest-energy Fermi-LAT point does not connect smoothly to the TeV SED, which might be due to a cross-calibration issue. In our model, the γ-ray SED can be adjusted by changing the IR background density, meaning that the SED can be matched easily without changing the low-energy SED once the cross-calibration issue is resolved. Keeping this in mind, the one-zone model provides basic properties of 3C 58 phenomenologically. However, the degree of the cooling break p3p2 ≈ 0.5 we estimate from the model is different from theoretical expectation of unity for synchrotron cooling in uniform magnetic fields B. This discrepancy implies that the physical properties of the PWN are not homogeneous (e.g., Reynolds 2009) as is the case in other PWNe. We therefore continue our analysis using a multi-zone model.

### 3.2. Multi-Zone Model

The physical properties of 3C 58 (e.g., B) vary with distance from the central pulsar, and the energy input to the PWN changes with time because of the energy loss of the pulsar. In principle, these effects can be computed from first principles by using MHD simulations with relativistic shock acceleration, but the theory of shock acceleration and MHD simulations still require improvements. We therefore use a phenomenological power-law prescription (e.g., Reynolds 2009) for the spatial variation of magnetic field B = B0 (R/R0)αB and bulk speed V = V0 (R/R0)αV , where R is the distance to the emission zone and R0 is the distance from the pulsar to the termination shock. For diffusion, we use Bohm's law DB–1Ee so D = D0 (R/R0)αB (γe/108). For the value of D0, we are guided by previous studies (e.g., Tang & Chevalier 2012; Porth et al. 2016) and further adjust it to match the observed SED. The model-SED shape is largely insensitive to the value of D0, but in general larger values of D0 produce broad humps in the IR and TeV bands (or deficits in the X-ray and GeV bands). Note that a simple broken-power law electron distribution (Figure 3) is sufficient to match the SED in our multi-zone model because the cooling break is automatically produced in the optical band (~ 1015 – 1016 Hz).

Injected and summed (spatially and temporally) particle distributions for the low-age case (2900 years).

We inject particles at the termination shock (R = R0) and follow their evolution out to the PWN boundary at RPWN = 3.7 pc (assuming a distance of 3.2 kpc; Torres et al. 2013). Note that the size RPWN is not a precise value and scales with the distance. We assume that the injection power decreases with time like ηL0(1 + t/τ0)–(n+1)/(n–1) and τ0 = 2τc/(n–1)–tage due to the pulsar's energy loss (Pacini & Salvati 1973; Torres et al. 2013), where η is the particle conversion efficiency, L0 is the initial luminosity, τc is the pulsar's characteristic age, and n = 1.1 (5400 yr) or n = 2 (2900 and 3800 yrs) is the braking index that quantifies the spin-down of the pulsar. We further assume that toroidal magnetic flux is conserved (αB +αV = –1; Reynolds 2009). Note that time variability of B is implicitly included in the prescription B = B0(R/R0)αB. While a change of B in time may have a different functional form, it is small and unlikely to be important.

For calculating SEDs, we use 1000 radial (linear), 10 000 energy (logarithmic), and 100 000 temporal (linear) bins. At a given time step we calculate the SED of the radiation from the particles in each radial bin. The particles lose energy via radiation (synchrotron and Compton) and expansion, meaning that the distribution at a given position varies. At the subsequent time step, we move the particles in space using advection (V = V0(R/R0)αV) and diffusion (random walk with a Gaussian width of $\sqrt{2Dt}$). At their new location, the particles cool via adiabiatic expansion (V), synchrotron radiation (B), and inverse-Compton upscattering. The latter two processes are used to compute the SED as function of position. Because particle diffusion is implemented as random walks, we run 100 simulations and take the average for each time step. We repeat this process for the age of the PWN, and numerically integrate over all radial and temporal bins (Figure 4) to compute a spatially integrated SED of the PWN. Before comparing the model SED to the observed SED of 3C 58, we verified that our numerical integration is accurate to ~10% by comparing the result with that derived from an analytic expression obtained for a simple power-law electron distribution with synchrotron-energy loss (Longair 2011).

Multi-zone models of 3C 58 for different ages: 2900 yr, 3800 yr, and 5400 yr, respectively, from top to bottom. Data points are the same as those in Figure 2.

Because there are many model parameters, we constrain some of them prior to modeling the SED. X-ray images constrain the location of the termination shock to R0 ≈ 0.1 pc (inner nebula emission; Slane et al. 2004). The bulk-speed parameters V0 and αV can be constrained by observations provided that the size, the age, and the expansion speed of the PWN are known, from

 $\begin{array}{c}{t}_{\mathrm{a}\mathrm{g}\mathrm{e}}={R}_{0}^{{\alpha }_{V}}\left[{R}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(1-{\alpha }_{V}\right)}-{R}_{0}^{\left(1-{\alpha }_{V}\right)}\right]/{V}_{0}\left(1-{\alpha }_{V}\right)\\ and\\ {V}_{0}=0.002c{\left({R}_{\mathrm{m}\mathrm{a}\mathrm{x}}/{R}_{0}\right)}^{-{\alpha }_{V}}.\end{array}$

Assuming that the bulk flow and the magnetic field do not increase with distance (i.e., αV , αB < 0) and that the bulk speed V is higher than the measured expansion speed of 0.002c for 3C 58 at the outer boundary (Bietenholz 2006), we infer that the PWN age is 2900–5600 yr. Note that this range is only approximate since the size (no sharp image boundary) and expansion speed are not accurately known.

We calculate SEDs for three realistic PWN ages (2900 yr, 3800 yr, and 5400 yr) to include limiting cases (αV or αB = 0). For each age, we compute V0 and αV and adjust the other parameters N0, p1,2, γe,b, B0, D0, and Galactic IR background energy density to match the observed SED. The results are presented in Figure 4 and Table 3. At high ages, radiation is inefficient (dominant adiabatic cooling), requiring that the pulsar's spin-down power needs to convert into particle wind more efficiently (e.g., n ≈ 1; Kim & An 2019). The 5400-yr model requires n = 1.1 which appears extreme but is needed to avoid η > 1. It is possible to match the observed SED using n ≈ 2 and η < 1 by varying the other parameters (e.g., αB). However, in our idealized model, these parameters are determined by observations and toroidal magnetic-flux conservation, meaning they cannot vary freely.

Parameters of the multi-zone models used in Figure 4

The agreement between data and model is generally good for all three cases. In the 5400-yr model, the shape of the radio SED is determined mainly by adiabatic cooling and thus depends only weakly on p1. Various values of p1 describe the data equally well; we therefore do not inject any radio-emitting particles below γe = 1.1×104 (no p1, Table 3) in the 5400-yr model presented in Figure 4. We note that the 3800-yr model curves smoothly near the radio break and underpredicts data points in the Planck and the Herschel bands.

## 4. DISCUSSION AND CONCLUSIONS

We investigated the physical properties of the pulsar wind nebula 3C 58 using Herschel data and SED models. We confirmed that the 500 μm flux density is lower than expected from interpolation of the fluxes observed in nearby bands even after considering systematic uncertainties in the estimates of the background emission. We added two more data points to the SED even though the new high-frequency measurements do not significantly change the SED shape due to large uncertainties. We updated the broadband SED by including new X-ray measurements, and used emission models to study flow properties in 3C 58.

We found that the Herschel flux densities were lower than what is expected from interpolation of Planck and Spitzer photometric data. While the measurements at 250 μm and 350 μm are highly uncertain, that at 500 μm is quite robust; an independent analysis of the same data found a consistent value (Arnaud et al. 2016). The highest-frequency Planck and the 500-μm Herschel points suggest that the flux decreases around 500 μm and then turns up, causing a small hump at ~ 1011 Hz. While of limited statistical significance, we speculate that there may be two populations of electrons: one emitting in the radio to Herschel band (subdominant at higher frequencies) and another at higher frequencies. Morphological differences between the Herschel and the Chandra images (internal structure) likewise support the idea that (part of the) IR emission may be produced by a different physical component. Alternatively, the deficit in the 500-μm flux density (or excess in the other IR bands) may merely be a statistical fluctuation, or the excess may be produced by dust emission (e.g., De Looze et al. 2019). Additional IR observations will be necessary for arriving at a robust conclusion.

Using one-zone modeling, we could reliably locate the break frequencies. We found the low-energy break to be at γe = 3.1 × 104, corresponding to an emission frequency of ~ 1011 Hz. Shock-acceleration theory (e.g., Sironi et al. 2015) suggests that pulsar-wind electrons with a relativistic Maxwellian energy distribution are injected into the termination shock, and that the distribution develops a non-thermal high-energy tail via shock acceleration (Sironi et al. 2015). In this case, the SED below the radio break could be a residual of the relativistic Maxwellian distribution injected by the pulsar (radio to ~ 1011 Hz emission) as previously hypothesized (e.g., Dubus et al. 2015); the particle Lorentz factor we infer is similar to that expected in pulsar winds (Aharonian et al. 2012). This may be the additional population of electrons hinted at by the IR SED. Including an additional electron population easily explains the IR SED. In this case IC scattering may leave a distinctive feature in the γ-ray SED; further investigation of this feature is ongoing. Either way, the existence of the IR hump is not beyond doubts and requires observational confirmation.

The one-zone model provided estimates of the emission power of the two SED humps: power emitted via Compton scattering (PIC) and synchrotron radiation (Psync). For a given IR photon energy density (Urad), the magnetic energy is given by UB = B2/8π = UradPsync/PIC, resulting in B ≈ 30 μG. The break energy γe,b2 = 107 suggests the age of the PWN to be ~ 2800 yr which is close to the range inferred from the size and expansion speed (2900–5600 yrs; Section 3.2). Note that one-zone modeling provides phenomenological and spatially "averaged" quantities. The parameters are highly correlated and different sets of parameters can explain the SED equally well.

While our multi-zone model is conceptually similar to previous ones (Torres et al. 2013; Tanaka & Takahara 2013; Li et al. 2018), there are two key differences. First, we assumed a lower FIR spectral energy than previous studies to match the Planck and Herschel data. For a given X-ray SED, this implies that the cooling break is at a higher frequency. Second, the maximum electron energy we assumed is lower than those of previous models, and therefore our SED model cuts off at a lower energy (~25 keV) in the hard X-ray band. This is necessary to fit the tentative spectral cutoff at 25 keV (An 2019). Although the current X-ray SED of 3C 58, without good measurements in the 10 keV–1MeV band, does not show a clear cutoff, future deep observations in this band (e.g., with NuSTAR and AMEGO) will certainly help to measure the cutoff energy more accurately; the cutoff energy can then be compared with those expected in relativistic-shock acceleration.

Assuming a power-law dependence of the flow speed V and a distance of 3.2 kpc, the current measurements of the size and expansion speed constrain the age of 3C 58 to 2900 – 5600 years. Different scaling relations would imply that the bulk flow or the magnetic field increase with distance from the central pulsar, which appears unphysical.

In our multi-zone modeling we find that the SED shape is slightly different for the middle-age model (3800 yr); in the Planck to Herschel band, the middle-age model predicts a smoother transition than the others do (Fig. 4). This is due to the cooling of the low-energy electrons being dominated by adiabatic expansion rather than (inefficient) synchrotron radiation (τad ≈ 200 yr vs τsy ≈ 20 000 yr) (e.g., Longair 2011). The injected electron distribution has a sharp break and thus the emission should likewise display a sharp break when cooling is insignificant. However, the broadband SED is a superposition of emission from all radial zones (cooled for different times) in the entire PWN, and the sharp break may be blurred.

In our model, the effect of adiabatic cooling for 1011-GHz emitting particles increases with the PWN age (long cooling time and smaller αV). The cooling is very weak for the 2900-yr model, and therefore the radio spectrum (p1 ≈ 1) and the shape of the radio–IR break directly re flect the injected particle distribution (sharp break). For the 5400-yr model, the adiabatic cooling is very strong, and particles at the break energy (γe,b) cool down to very low energies and emit sub-GHz photons. In this case, the radio SED shape is mainly determined by the superposition of cooled particles, not by the injection (no injection of radio emitting particles in our model), and the emission from cooled (radio) and uncooled (IR) particles can be distinguished. The break in the SED still remains rather sharp.

The 3800-yr case is located between the two extreme cases; the particles do not adiabatically cool to very low energies, and the transition from cooled to uncooled SED is blurred (smooth). Assuming that the true age of 3C 58 is located near 3800 yr, the multi-zone models also suggest an unmodeled population of particles emitting in the Planck to Herschel band – an observation similar to our discussion of the IR SED shape, although the excess may be due to external dust emission. Conversely, based on the good agreement between model SED and data, one may argue that the most likely ages of 3C 58 are either 2900 yrs or 5400 yrs. In this case, either the flow speed or the magnetic field is expected to be constant across the PWN. Spatially resolved spectroscopy of optical knots (to derive expansion speeds) can help to distinguish the models.

To summarize, our results suggest that there may be multiple populations of electrons in 3C 58. Unfortunately, the Herschel measurements are of limited value and the model parameters we infer are not necessarily unique. Therefore, we cannot yet provide a firm conclusion from our study. More observational and theoretical input is required. IR background density measurements can help to reduce the number of free parameters of the model, and further theoretical studies of relativistic shocks, MHD flows, and improved emission models can help to understand PWN physics better. Finally, diffusion coefficients in PWNe may not follow Bohm's law. The coefficients can be inferred by simultaneously modeling spatially resolved and integrated SEDs, which we plan to do in the near future. With its well-measured broadband SED and spatially resolved X-ray spectra, 3C 58 is certainly one of the best targets for this type of studies.

## Acknowledgments

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (NRF-2017R1C1B2004566).

## References

• Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Discovery of Pulsations from the Pulsar J0205+6449 in SNR 3C 58 with the Fermi Gamma-ray Space Telescope, ApJ, 699, L102.
• Aharonian, F. A., Bogovalov, S. V., and Khangulyan, D. 2012, Abrupt Acceleration of a `Cold' Ultrarelativistic Wind from the Crab Pulsar, Nature, 482, 507. [https://doi.org/10.1038/nature10793]
• Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014, Discovery of TeV γ-ray Emission from the Pulsar Wind Nebula 3C 58 by MAGIC, A&A, 567, L8.
• An, H. 2019, NuSTAR Hard X-ray Studies of the Pulsar Wind Nebula 3C 58, ApJ, 876, 150. [https://doi.org/10.3847/1538-4357/ab18a6]
• An, H., & Romani, R. W. 2017, Light Curve and SED Modeling of the Gamma-ray Binary 1FGL J1018.6–5856: Constraints on the Orbital Geometry and Relativistic Flow, Apj, 838, 145. [https://doi.org/10.3847/1538-4357/aa6623]
• An, H., Madsen, K. K., Reynolds, S. P., et al. 2014, High-energy X-ray Imaging of the Pulsar Wind Nebula MSH 15-52: Constraints on Particle Acceleration and Transport, ApJ, 793, 90. [https://doi.org/10.1088/0004-637X/793/2/90]
• Arnaud, M., Ashdown, M., Atrio-Barandela, F., et al. 2016, Planck Intermediate R esults – XXXI. Microwave Survey of Galactic Supernova Remnants, A&A, 586, A134.
• Bietenholz, M. F. 2006, Radio Images of 3C 58: Expansion and Motion of Its Wisp, ApJ, 645, 1180. [https://doi.org/10.1086/504584]
• Ciesla, L., Boselli, A., Smith, M. W. L., et al. 2012, Submillimetre Photometry of 323 Nearby Galaxies from the Herschel Reference Survey, A&A, 543, A161. [https://doi.org/10.1051/0004-6361/201219216]
• De Looze, I., Barlow, M. J., Bandiera, R. et al. 2019, The Dust Content of the Crab Nebula, MNRAS, 488, 164. [https://doi.org/10.1093/mnras/stz1533]
• Drury, L. O. 1983, An Introduction to the Theory of Diffusive Shock Acceleration of Energetic Particles in Tenuous Plasmas, Rep. Prog. Phys., 46, 973. [https://doi.org/10.1088/0034-4885/46/8/002]
• Dubus, G., Lamberts, A., & Fromang, S. 2015, Modelling the High-energy Emission from Gamma-ray Binaries Using Numerical Relativistic Hydrodynamics, A&A, 581, A27. [https://doi.org/10.1051/0004-6361/201425394]
• Finke, J. D., Dermer, C. D., & Böttcher, M. 2008, Synchrotron Self-Compton Analysis of TeV X-ray-selected BL Lacertae Objects, ApJ, 686, 181. [https://doi.org/10.1086/590900]
• Green, D. A., & Scheuer, P. A. G. 1992, Upper Limits on the Infrared Flux Density of the Filled-centre Supernova Remnant 3C 58, MNRAS, 258, 833. [https://doi.org/10.1093/mnras/258.4.833]
• Kang, H., 2010, Cosmic Ray Spectrum in Supernova Remnant Shocks, JKAS, 43, 25. [https://doi.org/10.5303/JKAS.2010.43.2.025]
• Kargaltsev, O., Klingler, N., Chastain, S., & Pavlov, G. G. 2017, Toward Understanding the Physical Underpinnings of Spatial and Spectral Morphologies of Pulsar Wind Nebulae, J. Phys. Conf. Ser., 012050. [https://doi.org/10.1088/1742-6596/932/1/012050]
• Kennel, C. F., & Coroniti, F. V. 1984, Confinement of the Crab Pulsar's Wind by Its Supernova Remnant, ApJ, 283, 694. [https://doi.org/10.1086/162356]
• Kim, M., & An, H., 2019, Measuring Timing Properties of PSR B0540-69, JKAS, 52, 41.
• Lemoine, M., Kotera, K. & Pétri, J. 2015, On Ultra-high Energy Cosmic Ray Acceleration at the Termination Shock of Young Pulsar Winds, JCAP, 7, 16. [https://doi.org/10.1088/1475-7516/2015/07/016]
• Li, J., Torres, D. F., Lin, T. T., et al. 2018, Observing and Modeling the Gamma-ray Emission from Pulsar/Pulsar Wind Nebula Complex PSR J0205+ 6449/3C 58, ApJ, 858, 84. [https://doi.org/10.3847/1538-4357/aabac9]
• Longair, M. S. 2011, High Energy Astrophysics (Cambridge: Cambridge University Press).
• Lyutikov, M., Temim, T., Komissarov, S., et al. 2018, Interpreting Crab Nebula's Synchrotron Spectrum: Two Acceleration Mechanisms, MNRAS, 489, 2. [https://doi.org/10.1093/mnras/stz2023]
• Madsen, K. K., Reynolds, S., Harrison, F., et al. 2015, Broadband X-ray Imaging and Spectroscopy of the Crab Nebula and Pulsar with NuSTAR, ApJ, 801, 66. [https://doi.org/10.1088/0004-637X/801/1/66]
• Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, HerMES: The SPIRE Confusion Limit, A&A, 518, L5.
• Nynka, M., Hailey, C. J., Reynolds, S. P., et al. 2014, NuSTAR Study of Hard X-ray Morphology and Spectroscopy of PWN G21.5-0.9, ApJ, 789, 72. [https://doi.org/10.1088/0004-637X/789/1/72]
• Pacini, F., & Salvati, M. 1973, On the Evolution of Supernova Remnants. Evolution of the Magnetic Field, Particles, Content, and Luminosity, ApJ, 186, 249. [https://doi.org/10.1086/152495]
• Porth, O., Vorster M. J., Lyutikov, M. and Engelbrecht N. E. 2016, Diffusion in Pulsar Wind Nebulae: An Investigation Using Magnetohydrodynamic and Particle Transport Models, MNRAS, 460, 4135. [https://doi.org/10.1093/mnras/stw1152]
• Reynolds, S. P. 2009, Synchrotron-loss Spectral Breaks in Pulsar-wind Nebulae and Extragalactic Jets, ApJ, 703, 662. [https://doi.org/10.1088/0004-637X/703/1/662]
• Roberts, D. A., Goss, W. M., Kalberla, P. M. W., et al. 1993, High Resolution Hi Observations of 3C 58, A&A, 274, 427.
• Sironi, L., Keshet, U. & Lemoine, M. 2015, Relativistic Shocks: Particle Acceleration and Magnetization, SSRv, 191, 519. [https://doi.org/10.1007/s11214-015-0181-8]
• Slane, P. O., Helfand, D. J., & Murray, S. S. 2002, New Constraints on Neutron Star Cooling from Chandra Observations of 3C 58, ApJ, 571, L45. [https://doi.org/10.1086/341179]
• Slane, P. O., Helfand, D. J., van der Swaluw, E. & Murray, S. S. 2002, New Constraints on the Structure and Evolution of the Pulsar Wind Nebula 3C 58, ApJ, 616, 403. [https://doi.org/10.1086/424814]
• Slane, P. O., Helfand, D. J., Reynolds, S. P., et al. 2008, The Infrared Detection of the Pulsar Wind Nebula in the Galactic Supernova Remnant 3C 58, ApJ, 676, L33. [https://doi.org/10.1086/587031]
• Stephenson, F. R. 1971, A Suspected Supernova in AD 1181. IAU Colloq., 8, 10 (Cambridge: Cambridge University Press). [https://doi.org/10.1017/S0252921100028153]
• Tanaka, S. J., & Takahara, F. 2013, Properties of Young Pulsar Wind Nebulae: TeV Detectability and Pulsar Properties, MNRAS, 429, 2945. [https://doi.org/10.1093/mnras/sts528]
• Tang, X., & Chevalier, R. A. 2012, Particle Transport in Young Pulsar Wind Nebulae, ApJ, 752, 83. [https://doi.org/10.1088/0004-637X/752/2/83]
• Torres, D. F., Cillis, A. N. & Martín Rodriguez, J. 2013, An Energy-conserving, Particle-dominated, Time-dependent Model of 3C 58 and Its Observability at High Energies, ApJ, 763, L4. [https://doi.org/10.1088/2041-8205/763/1/L4]
• Valtchanov, I. 2017, The Spectral and Photometric Imaging Receiver (SPIRE) Handbook, Herschel Explanatory Supplement, https://www.cosmos.esa.int/web/herschel/legacy-documentation-spire.
• Yüksel, H., Kistler, M. D. & Stanev, T. 2009, TeV Gammarays from Geminga and the Origin of the GeV Positron Excess, Phys. Rev. Lett., 103, 051101. [https://doi.org/10.1103/PhysRevLett.103.051101]

### Figure 1.

Herschel SPIRE three-band images of 3C 58. The center wavelengths are 500μm, 350μm, and 250μm, respectively, from left to right. White solid-line ellipses indicate the source region, dashed lines denotes background regions.

### Figure 2.

SED data and a one-zone model for 3C 58. Data are from previous works. Radio (red): Green & Scheuer (1992), Planck (cyan): Arnaud et al. (2016), Herschel (orange): this work, IRAS (black): Arnaud et al. (2016), Spitzer (red): Slane et al. (2008), X-ray (brown and green): Slane et al. (2004) and An (2019), GeV (purple): Li et al. (2018), and TeV (blue): Aleksić et al. (2014). Note that there are two spectral breaks in the low-energy hump because we used a power law with two breaks for the electron energy distribution.

### Figure 3.

Injected and summed (spatially and temporally) particle distributions for the low-age case (2900 years).

### Figure 4.

Multi-zone models of 3C 58 for different ages: 2900 yr, 3800 yr, and 5400 yr, respectively, from top to bottom. Data points are the same as those in Figure 2.

### Table 1

Flux densities of 3C 58 measured with Herschel

Frequency [GHz] Flux [Jy] Uncertainty [Jy]
600 2.2 0.4
857 1.6 0.5
1200 1.4 0.7

### Table 2

Parameters of the one-zone SED model used in Figure 2

Parameter Value
p1, p2, p3 [1, 3, 3.5]
γe,b1, γe,b2, γe,max [3.1 × 104, 107, 4 × 108]
B(magnetic field) 30 μG
Background CMB (2.7K) density 0.25 eV/cm3
Background IR (20K) density 0.8 eV/cm3

### Table 3

Parameters of the multi-zone models used in Figure 4

Parameter 2900 years 3800 years 5400 years
p1 1.0 0.8 · · ·
p2 2.87 2.88 2.79
γe,b 9.3 × 104 7.4 × 104 1.1 × 104
γe,max 5 × 108 3.2 × 108 2.5 × 108
B0 40 μG 140 μG 300 μG
V0 0.077 c 0.013 c 0.002 c
R0 0.1 pc 0.1 pc 0.1 pc
αB, αV 0, –1 –0.5, –0.5 –1, 0
BG IR (20 K) density 1.7 eV/cm3 5 eV/cm3 25 eV/cm3
D0 2 × 1027 cm2/s 1 × 1026 cm2/s 1 × 1026 cm2/s
n 2 2 1.1