Multiband nonthermal radiative properties of pulsar wind nebulae

Aims. The nonthermal radiative properties of 18 pulsar wind nebulae (PWNe) are studied in the 1D leptonic model. Methods. The dynamical and radiative evolution of a PWN in a nonradiative supernova remnant are self-consistently investigated in this model. The leptons (electrons/positrons) are injected with a broken power-law form, and nonthermal emission from a PWN is mainly produced by time-dependent relativistic leptons through synchrotron radiation and inverse Compton process. Results. Observed spectral energy distributions (SEDs) of all 18 PWNe are reproduced well, where the indexes of low-energy electron components lie in the range of 1.0–1.8 and those of high-energy electron components in the range of 2.1–3.1. Our results show that FX/Fγ > 10 for young PWNe; 1 < FX/Fγ ≤ 10 for evolved PWNe, except for G292.0+1.8; and FX/Fγ ≤ 1 for mature/old PWNe, except for CTA 1. Moreover, most PWNe are particle-dominated. Statistical analysis for the sample of 14 PWNe further indicate that (1) not all pulsar parameters have correlations with electron injection parameters, but electron maximum energy and PWN magnetic field correlate with the magnetic field at the light cylinder, the potential difference at the polar cap, and the spindown power; (2) the spin-down power positively correlates with radio, X-ray, bolometric, and synchrotron luminosities, but does not correlate with gamma-ray luminosity; (3) the spin-down power positively correlates with radio, X-ray, and γ-band surface brightness; and (4) the PWN radius and the PWN age negatively correlate with X-ray luminosity, the ratio of X-ray to gamma-ray luminosities, and the synchrotron luminosity.


Introduction
A pulsar wind nebula (PWN) is formed when the wind of a pulsar interacts with the ambient medium, either the supernova ejecta or the interstellar medium (e.g., Rees & Gunn 1974;Kennel & Coroniti 1984;Chevalier 2004;Gaensler & Slane 2006;Bucciantini 2008). It is generally believed that a PWN is mainly composed of a relativistic nonthermal lepton (electron/positron) plasma and magnetic field and can emit nonthermal photons from radio to γ-rays via synchrotron radiation and the inverse Compton (IC) process. Observations have shown that some PWNe, such as the Crab nebula, Kes 75, MSH 15-52, can emit nonthermal emission from radio to very high-energy (VHE) γ-ray bands (e.g., Bühler & Blandford 2014;Reynolds et al. 2017). Currently, about 37 TeV PWNe are firmly identified 1 . More recently, the observations and relevant physical analysis of 19 TeV PWNe have been presented in Abdalla et al. (2017), where 14 TeV PWNe are firmly identified as PWNe by HESS observations. These observations provide an essential reason for studying PWNe.
Various models have been proposed to explain the nonthermal properties of PWNe and a brief review of current models and differences can be found in Torres et al. (2014). Here we focus on the one-dimensional (1D) model of the dynamical and radiative evolution of a PWN in a nonradiative SNR. Here we focus on the 1D model of the dynamical and radiative evolution of a PWN in a nonradiative supernova remnant (SNR) presented in Gelfand et al. (2009). In the frame of the model, the radiative properties during different phases of the PWN evolution 1 http://tevcat.uchicago.edu/ are investigated with a single power-law injection spectrum for the electrons/positrons (Gelfand et al. 2009), with the relativistic Maxwellian and a high-energy power-law tail injection spectrum for the electrons/positrons (Fang & Zhang 2010a), and with two possible forms of injected electron spectra: a broken power law and the sum of a power law at low energy, and a relativistic Maxwellian plus a high-energy power-law tail . It has been shown that the broken power-law form for the electron injection is required when the model applies to the Crab nebula . Recently, a model for describing PWN radiative properties during the dynamical evolutions of the PWN and the SNR is presented in Martin et al. (2016), the radius of the PWN during the free expansion phase and compression in this model is calculated through solving the equations given by Chevalier (2005), where the prescription is similar to that in Gelfand et al. (2009). The model is applied to a TeV PWN CTA 1 (Martin et al. 2016). We note that there is a kind of model in which a broken power-law injection spectrum for the electrons/positrons is assumed but the dynamics beyond reverberation is not included (e.g., Zhang et al. 2008;Bucciantini et al. 2011;Tanaka & Takahara 2011).
To investigate nonthermal radiative properties of each PWN and statistical features of PWNe, observed multiwaveband data for each PWN are required. However, not all PWNe are observed at different bands, due to observation limits, so 18 PWNe from known PWNe are selected here according to the following criteria: (1) the period and period derivative of central pulsar in each PWN are known, and (2) nonthermal emission at radio, X-ray, and TeV bands are detected. In this selection, N 158A is not detected at TeV band, G310.6-1.6 is detected with an upper limit at TeV band, G292.0+1.8 is detected with upper limits at GeV

Model description
For completeness, we briefly review the model for dynamical and radiative evolution of a PWN inside a SNR given in Zhu et al. (2015), which is closely based on that developed by Gelfand et al. (2009; see also Fang & Zhang 2010a). In this model, the large-scale evolution of a composite SNR depends on the mechanical energy E sn of the explosion, the density n H of the ambient medium, the mass M ej of the supernova ejecta, and the spin-down power L(t) of the pulsar. For a given pulsar with a period P, a period derivativeṖ, a braking index n, and an initial spin-down power L 0 , the spin-down power L(t) at time t is given by (Gaensler & Slane 2006) where τ 0 = (2τ c )/(n − 1) − T age is the initial spin-down timescale of the pulsar, τ c = P/2Ṗ is the characteristic age of the pulsar, and T age is the age of the PWN. Since L(t) can be estimated by L(t) = 4π 2 IṖ/P 3 , where I is the pulsar moment of inertia (here I = 10 45 g cm 2 is used), the evolution of the spin-down power can be determined (i.e., L 0 and τ 0 can be estimated) if T age is known. Therefore, the main parameters of a given pulsar are P, P, and n. At present, the vaules of the braking index n for some pulsars have been measured, and we use the measured values if available, otherwise we assume n = 3. We call these parameters the pulsar and ejecta parameters. The evolution of isotropic electron distribution N(E, t) in the PWN is calculated by using whereĖ is the energy-loss rate of the particles with an energy E which includes the contributions of synchrotron radiation, inverse Compton scattering, and adiabatic losses and τ(E, t) is the escape time (for the details of their processes, see Zhang et al. 2008). The last term on the right-hand side in Eq.
(2), Q(E, t), is a source term, i.e., the injected electron number per unit energy per unit time, and is assumed to be a broken power-law form where Q 0 (t) is a time-dependent parameter, E b is the break energy, α 1 < 2 and α 2 > 2 are respectively the spectral indexes of the injection rate at E ≤ E b and E > E b , and E max is the maximum energy of the electrons.
Here the spin-down power of the pulsar is assumed to be distributed between electrons (Ė e (t) = η e L(t)) and magnetic fields (Gelfand et al. 2009) On the other hand, the maximum energy of the electrons can be estimated by the condition that the Larmor radius of the electrons inside the PWN is smaller than the termination shock radius of the PWNe by the containment factor < 1, which is given by where e is the electron charge (e.g., Zhu et al. 2015). We note that E max can also be estimated by the balance between the synchrotron loses and acceleration, which gives E syn max ≈ [3(m e c 2 ) 2 /4e] π/eB pwn . For the parameters of PWNe used here (see Tables 1-3), E max < E syn max is always satisfied. Hence, the injection parameters involving electron injection are E b , α 1 , α 2 , η B , and .
In such a model, the dynamical and radiative properties for a given PWN can be self-consistently studied Gelfand et al. (2009; see also Fang & Zhang 2010a;. On the one hand, the time evolutions of the SNR radius R snr (t), reverse shock radius R rs (t), PWN radius R pwn (t), the position of the neutron star R psr (t), the termination shock radius R ts (t), and the magnetic field B pwn (t) of the PWN can be calculated. On the other hand, the time evolutions of electron spectrum and the spectral energy distribution (SED) of nonthermal photons can be obtained. It should be pointed out that different electron injection forms and electron energy losses will lead to the change in dynamical features of the PWN. In this paper, the broken power-law injection of the electrons (see Eq. (3)) and escape term of the electrons (see the second term on the left-hand side of Eq. (2)) are considered, so the dynamical features for a given PWN will be different from those given in Gelfand et al. (2009), who assumed a single power-law injection without electron escape (e.g., Crab nebulae, see Zhu et al. 2015). We note that both R pwn (t) and B pwn (t) play important roles in the evolution of the energy losses and the SED for a given PWN. Therefore, in our calculations the allowed ranges of R pwn (t) are limited by the values given in Abdalla et al. (2017; see their Tables 1 and 3).
In our calculations, nonthermal photons are produced through synchrotron radiation and inverse Compton (IC) scattering of relativistic electrons (e.g., Zhu et al. 2015). For the synchrotron radiation, the emissivity given in Zhang et al. (2008) is used, which includes the effect of electron pitch angle. The magnetic field B pwn of a PWN evolves with time and can be estimated by Eq. (A.6) in Appendix A. We note that Torres et al. (2014) did not consider the effect of electron pitch angle (see also Tanaka & Takahara 2010;Martin et al. 2012), which will result in parameter differences between our model and the model in Torres et al. (2014) for a given PWN.
For IC scattering, soft photon fields consist of four components: the cosmic microwave background (CMB) in our calculations, the galactic far-infrared (FIR) background, the near-infrared (NIR) and optical photon field due to the stars, and synchrotron radiation. The energy density and temperature of the CMB photons are U CMB = 0.25 eV cm −3 and T CMB = 2.7 K, and the energy densities and temperatures of FIR (U FIR and T FIR ) and NIR (U NIR and T NIR ) can be different for PWNe. We refer to these parameters as soft photon parameters. The emissivity of IC scattering used here is given by Zhang et al. (2008).
To fit the observed spectral energy distribution (SED) for each PWN, pulsar and ejecta parameters and soft photon parameters are fixed, and injection parameters are considered as fitting parameters (see Tables 1-3). The Levenberg-Marquardt (LM) method of the χ 2 minimization fitting procedure is used to find the best-fitting values of injection parameters and their uncertainties. However, because of the lack of data (from radio to optical) for most PWNe, the best-fitting values and their uncertainties of α 1 and E b cannot be easily obtained. Therefore, at first, visual fitting is used to determine the values of α 1 and E b , and then the LM method is used to obtain the values of α 2 , η B , and . With this procedure, we fit 14 PWNe, except for N 158A, G310.6-1.6, G292.0+1.8, and HESS J1303-631. Since the observed data at GeV-TeV and from radio to X-ray band are upper limits for these four PWNe, the LM method can only determine the values of α 2 and η B for N 158A, G310.6-1.6, and G292.0+1.8, and α 2 for HESS J1303-631; other injection parameters are estimated via visual fitting. These values are shown in Tables 1-3.
In addition to the above parameters, six derived parameters (E max , B pwn , R pwn , R rs , R snr , and R ts ) are also given (see Tables 1-3). We note that E max and B pwn , which relate to and η B , are calculated. Because the slight difference in the injection spectrum has little effect on the dynamical properties of PWNe, the uncertainties of R pwn , R rs , R snr , and R ts are not calculated here.

Applications to individual PWNe
The model described in Sect. 2 is now applied to explain the observed SEDs of nonthermal photons for 18 PWNe. The detailed descriptions of observed and derived properties of each PWN are summarized in Appendix B. The PWN sample is divided into three groups based on possible ages T age of PWNe. The first group consists of the PWNe with T age ≤ 2400 yr (we call them young PWNe), the second group with 2400 < T age < 5000 yr (called evolved PWNe), and the third group with T age ≥ 5000 yr (called mature/old PWNe). Although the division of the three groups is arbitrary, it is convenient for us to study the properties of PWNe that lie in different age ranges.   Martin et al. 2014;Zhu et al. 2015;Martin et al. 2016). Two important quantities that influence the electron energy loss and photon SED of a PWN are R pwn (t) and B pwn (t). In Torres et al. (2014), R pwn (t) ∝ t 6/5 is calculated during the SNR's free expansion phase (van der Swaluw et al. 2001(van der Swaluw et al. , 2003 and B pwn (t) is estimated with a method similar to that used in Gelfand et al. (2009). In Tanaka & Takahara (2011), B pwn (t) is calculated by using the magnetic field energy conservation, i.e., (4π/3)R 3 pwn (t)B 2 pwn (t)/(8π) = ηL(t )dt , and R pwn is estimated by reproducing photon SED.

Group 1: young PWNe
There are six young PWNe in this group, N 158A, G21.5-0.9, Crab nebula, Kes 75, G310.6-1.6, and 3C 58, whose ages range from 700 yr to 2400 yr. The parameters of pulsar and ejecta, electron injections, and soft photon fields for these young PWNe are listed in Table 1, and the details for each PWN is given in Appendix B. Using the parameters in Table 1, the SEDs of these young PWNe are calculated. The comparisons of the calculated and observed SEDs are shown in Fig. 1 and derived parameters are listed in Table 1. The multiband SEDs of these young PWNe have been studied (e.g., Tanaka   expansion phases and the dynamical evolution in the free expansion phase used here is similar that used in Torres et al. (2014). The fluxes, F X ≡ F 1−10 keV and F γ ≡ F 1−10 TeV , in X-ray (1-10 keV) and γ-ray (1-10 TeV) bands for each PWN are calculated.
At an age of 760 yr N 158A was not detected in GeV and TeV γ-ray bands, so the sensitivity curves of the Cherenkov Telescope Array (CTA) is used as an upper limit to restrict soft photon parameters (see Table 1). Our results are shown in Table 1 and Fig. 1, and give that R pwn = 0.68 pc and B pwn = 45.22 µG, which is roughly consistent with the values given in Martin et al. (2014, R pwn = 0.7 pc and B pwn = 32 µG). The observed SED can be reproduced and predicted as F X /F γ ∼ 118.
Our results for G21.5-0.9 are shown in Table 1 and Fig. 1. Here the age is assumed to be T age = 900 yr, which is close to T age = 870 yr in Torres et al. (2014). Our results give R pwn = 0.86 pc and B pwn ≈ 84 µG, which are consistent with those in Torres et al. (2014; R pwn = 0.9 pc and B pwn = 71 µG) and those in Tanaka & Takahara (2011; R pwn = 1.0 pc and B pwn = 64, and 47 µG). We note that R pwn is less than that given in Abdalla et al. (2017;<4 pc). The predicted F X /F γ ∼ 130.
For Crab nebulae, the results are the same as those given in Zhu et al. (2015). As mentioned above, the main differences between our results and those in Tanaka & Takahara (2011) and Torres et al. (2014) come from the use of different synchrotron emissivity. Our results give R pwn ≈ 1.95 pc (<3 pc in Abdalla et al. 2017) and B pwn ≈ 116 µG. The predicted For Kes 75, its age is assumed to be 1000 yr, which is slightly older than the ages given in Tanaka & Takahara (2011) and Torres et al. (2014), and its SED is calculated with the same values of pulsar and ejecta parameters and soft photon parameters as those in Torres et al. (2014) and can reproduce the observed SED well with F X /F γ ∼ 13 (see Fig. 1). However, our injection parameters are different from those given in Torres et al. (2014) because the formula of synchrotron emissivity used here (see Eq. (12) of Zhang et al. 2008) is different from that used in Torres et al. (2014, see Eqs. (27) and (28) Table 1 in comparison with those (R pwn = 0.9 pc and B pwn = 19 µG) in Torres et al. (2014) and those (R pwn = 0.29 pc and B pwn = 20 µG) in Tanaka & Takahara (2011).

Group 2: evolved PWNe
The group of evolved PWNe consists of HESS J1813-178, G54.1+0.3, G292.0+1.8, G0.9+0.1, MSH 15-52, and N 157B whose ages range from 2500 yr to 4600 yr. Since the SN explosion energy E sn are unknown for these PWNe, here we assume E sn = 1.0 × 10 51 erg. The related parameters for the evolved PWNe are listed in Table 2 and the comparisons of the modeled and observed SEDs are shown in Fig. 2. HESS J1813-178 is observed to have a distance d ≈ 4.7 kpc Halpern et al. 2012) and R pwn = 4.0 ± 0.3 pc . With T age = 2500 yr, R pwn ≈ 3.7 pc and B pwn = 8.11 µG are obtained here. The predicted F X /F γ ∼ 6. We note that Fang & Zhang (2010b) investigated the PWN in detail, and obtained R pwn = 1.7 pc and B pwn = 5 µG with T age = 1200 yr.
Although the age of G292.0+1.8 is uncertain, R pwn ≈ 3.0 pc is given in Bhalerao et al. (2015). With T age = 2700 yr , which is larger than 1600 yr given in Murdin & Clark (1979) and 2500 yr in Martin et al. (2014), our results show that R pwn = 3.12 pc and B pwn ≈ 24 µG, which are consistent with R pwn = 3.5 pc and B pwn = 16 µG obtained by  and R pwn = 3.5 pc and B pwn = 21 µG by Martin et al. (2014). The predicted F X /F γ ∼ 14.

Estimate of pair multiplicity and wind Lorentz factor
Pair multiplicity, κ multi , and the Lorentz factor, γ w , are calculated for a PWN based on our model. The pair multiplicity can be expressed as where Q is the injection rate by integrating Eq. (2), andṄ is the injection rate with electrodynamic minimum suggested by Goldreich & Julian (1969), i.e.,Ṅ = ( cIΩΩ e 2 ) 1/2 . The values are shown in Table 5. The pair multiplicity ranges from ∼10 5 to ∼10 8 . The Lorentz factor of the pulsar wind can be estimated by where E is the average energy for the SED of each PWN. The values of γ w are listed in Table 4. Clearly, E b is larger than E by up to several orders of magnitude.

Brief summary of calculated results
Here, the main results in the above calculations are briefly summarized. The observed SEDs of 18 PWNe can be reproduced well in the frame of the model. For the electron injection, five free parameters (α 1 , α 2 , E b , η B , ) are used. Our results give α 1 ≈ 1.0-1.8, α 2 ≈ 2.0-3.1, E b ≈ 10 4 -10 7 MeV, η B ≈ 0.1%-11% for the PWNe except for MSH 15-52 (η B ≈ 45%) and CTA 1 (η B ≈ 49%), and ≈ 0.1-0.70. The results for spectral indices are consistent with previous studies (e.g., Bucciantini et al. 2011;Tanaka & Takahara 2011;Torres et al. 2014). According to Sironi & Spitkovsky (2011, 2012, these low-energy electrons are accelerated by relativistic magnetic reconnection in a PWN; its spectral slope in the range at 1.0-2.0. Therefore, the low-energy electron component possibly originates from relativistic magnetic reconnection. The high-energy electron component with spectral slope from 2.0-3.0 are consistent with that produced by the Fermi-type process (e.g., Achterberg et al. 2001); thus, high-energy electrons are from the Fermi-type process in the termination shock. The PWN magnetization is weak, which means that most PWNe are particle-dominated. This conclusion is consistent with result of Torres et al. (2014). The values of are in the range 0.1-0.45, except for HESS J1119-614, and of the break energy E b are in the range 10 4 -10 7 MeV.

Correlation analysis of the PWNe
In this section, the correlations of various physical quantities in our model are studied. To this end, the best linear fit y = p 1 x+ p 0 to the data with the minimum χ 2 technique is used. For each fit, the values of p 1 and p 0 are given, and the Pearson's correlation coefficient r and the probability of the null hypothesis P null are calculated. Here, it is assumed that the correlation between two physical quantities is noneffective if P null ≥ 0.05. In correlation analysis, the sample of 14 PWNe is used in which α 2 , η B , and are obtained through the LM method. All results have been recalculated and various luminosities with uncertainties are listed in Table 6. The uncertainties of the various physical quantities are included. We note that although four PWNe (N 158A, G310.6-1.6, G292.0+1.8, and HESS J1303-631) are not included in the best linear fits, the values of relevant physical quantities of these four PWNe are listed in Tables 4-6 and shown in Figs. 4-10.
Here, the 2σ confidence band for each panel shown in Figs. 4−10 are from the dispersion of the fit to the points.

Correlations between pulsar parameters with injection spectrum parameters and with derived parameters
For a given central pulsar with a period P (s) and a period deriva-tiveṖ (s s −1 ) in a PWN, it is commonly assumed that the magnetic field around it can be approximated as a dipole, so the surface magnetic field B s is given by where c is the light speed, I = 10 45 g cm 2 is the inertia moment, and R = 10 6 cm is the surface radius of the pulsar. At the light cylinder with a radius R LC = (cP)/(2π) ≈ 4.77 × 10 9 P cm, the magnetic field is described as A110, page 10 of 20  For such a pulsar, the potential difference at the polar cap, Φ, can be expressed as Φ = 2π 2 B s R 3 (cP) 2 = 4.2 × 10 20 (P −3Ṗ ) 1/2 stat volts.
We note that the spin-down power L(t) at time t (see Eq. (3)) can be also expressed in terms of P andṖ as For each pulsar, the L(t) and τ c are shown in Tables 1-3, and B s , B LC , and Φ are listed in Table 5. The correlation analysis indicate that (1) there are no correlations between E b , η B , α 1 , α 2 , and the pulsar parameters such as B s , B BL , Φ, L(t), and τ c ; (2) two derived parameters  L(t) ∼ (PṖ) 1/2 , B LC ∼ P 5/2Ṗ1/2 , and Φ ∼ (P −3Ṗ ) 1/2 , these correlations are only consistent with the logical conclusions.

Correlations between luminosities at various bands and L(t)
The possible correlations between the luminosities at radio (1.4 GHz), L r ; X-ray (1-10 keV), L X ; VHE γ-ray (1-10 TeV), L γ ; bolometric luminosity at all wave bands, L bol (which can be obtained from our model results, the corresponding values and errors are listed in Table 5); and the pulsar's spin-down power, L(t), were analyzed. The results show that there are significant positive correlations between L r and L(t), between L X and L(t), and between L bol and L(t) (see Fig. 5), and the best linear fits yield L r ∼ L(t) 1.21±0.17 , L X ∼ L(t) 1.54±0.12 , and L bol ∼ L(t) 1.24±0.12 , respectively. The correlation between L X and L(t) is in agreement with the recent conclusions of Mattana et al. (2009), Kargaltsev & Pavlov (2010), and Kargaltsev et al. (2013). As expected, there is no correlation between L γ and L(t), which is consistent with Mattana et al. (2009), Kargaltsev & Pavlov (2010), and Kargaltsev et al. (2013). This phenomenon may be due to the differences in the IR background and uncertain distance.
Since L r /L γ and L X /L γ are distance-independent, the correlations between L r /L γ and L(t) and between L X /L γ and L(t) are analyzed. Our results are shown in Fig. 6 where the best linear fits give L R /L γ ∼ L(t) 1.28±0.28 , and L X /L γ ∼ L(t) 1.41±0.22 . In fact, the correlation between L X /L γ and L(t) is found in Mattana et al. (2009).

Correlations between synchrotron luminosity and L(t)
Because synchrotron radiation dominates over the radiation from radio to X-ray bands, it is important to study the correlations of synchrotron cooling break energy (E s ) and synchrotron radiation luminosity (L s ) with L(t). Following Tanaka Fig. 9. Correlations between R pwn and L X , L X /L γ , and L s . The solid lines are the best linear fits, which are log L X = −(1.36 ± 0.68) log R pwn + (37.03 ± 0.38), log L X /L γ = −(2.54 ± 0.52) log R pwn + (2.69 ± 0.31), and log L s = −(1.90 ± 0.75) log R PWNe + (37.48 ± 0.24); the correlation coefficients r = 0.45, 0.77, and 0.54; and the probability of the null hypothesis P null = 6.40 × 10 −2 , 1.73 × 10 −4 , and 2.18 × 10 −2 (from left to right). The dashed lines are the 2σ confidence bands for the sample.
the inverse Compton cooling is assumed to be ineffective in most of the evolutionary phases, and the synchrotron cooling break frequency ν s is given by Hz.
Using Eq. (13) and based on our model calculation, E s = hν s and L s are calculated for each PWN in our sample. The results show that there is a strong positive correlation between L s and L(t), the best linear fit yields L s ∼ L(t) 1.64±0.09 , and that there are strong negative correlation between E s and L(t), where E s ∼ L(t) −0.78±0.11 , and E s and L s , where L s ∼ E −1.15±0.11 s . The best linear fit results are shown in Fig. 7.

Correlations between surface brightness S and L(t)
The surface brightness is defined as (e.g., Abdalla et al. 2017) where L is the luminosity at a given band, R PWN is the radius of the PWN, F is the energy flux at a given band, and σ is its angular extent as seen from the Earth. Using our model and Eq. (14), the surface brightness, S R , S X , and S γ , at radio (1.4 GHz), X-ray (1-10 keV), and γ-ray (1-10 TeV) bands are calculated. The results are shown in Fig. 8 and show that there are correlations between L(t) and S R , S X , and S γ . The best linear fits give S r ∼ L(t) 1.78±0.20 , S X ∼ L(t) 1.89±0.32 , and S γ ∼ L(t) 0.71±0.17 . We note that S γ ∼ L(t) 0.81 has been given in Abdalla et al. (2017).

Correlations between luminosity and R pwn and T age
The correlations between the luminosities at various bands and R pwn are analyzed here. In our model, the PWN radius for each PWN can be calculated (see Tables 1-3). The results show that there are no correlations between R pwn and L r , L γ , and L bol , but there are correlations between R pwn and L X , L s , and L X /L γ . The best linear fits give that L X ∼ R −1.36±0.68 pwn , L s ∼ R −1.90±0.75 pwn , and L X /L γ ∼ R −2.54±0.52 pwn , respectively. We note that the correlation A110, page 13 of 20  Fig. 10. Correlations between T age and L X , L X /L γ , and L s . The solid lines are the best linear fits, which are log L X = −(2.44 ± 0.59) log T age + (44.23±1.89), log L X /L γ = −(3.37±0.29) log T age +(12.47±0.94), and log L s = −(2.83±0.53) log T age +(45.37±1.59); the correlation coefficients r = 0.72, 0.95, and 0.80; and the probability of the null hypothesis P null = 7.81 × 10 −4 , 3.15 × 10 −9 , and 7.22 × 10 −5 . The dashed lines are the 2σ confidence bands for the sample.
between L X and R pwn is marginal here. The results are shown in Fig. 9. The same procedure is performed for the correlations between T age and the luminosities at various bands, and almost the same correlations are found. We found that there are correlations between T age and L X , L s , and L X /L γ . These correlations mean that L X ∼ T −2.44±0.59 age , L s ∼ T −2.83±0.53 age , and L X /L γ ∼ T −3.37±0.29 age .
The best linear fit results are shown in Fig. 10.

Discussions and conclusions
In this paper, the leptonic model with a broken power-law injection for the electrons/positrons is applied to the sample of 18 PWNe. The sample is divided into three groups: young, evolved, and mature/old PWNe. Observed SEDs of all 18 PWNe can be reproduced well in this model (see Figs. 1-3). The model parameters obtained in our calculations are listed in the Tables 1−3, and the relevant discussion are described in Sect. 3.
Using a time-dependent modeling of PWNe given by Martin et al. (2012), Torres et al. (2014) modeled the SEDs of nine PWNe and studied their statistical properties. These nine PWNe are included in our sample. As mentioned in Sect. 1, the model given in Torres et al. (2014) does not include the dynamical evolution of a PWN, but our model does. Therefore, some model parameters obtained in these two models are different. Meanwhile, our sample is larger than some previous works. So it is worth comparing our correlation results with the results of Torres et al. (2014) and some previous works.
For the relation between the derived parameters and pulsar's parameter, our results indicate that the maximum electron energy has positive correlations with the magnetic field at the light cylinder, the potential difference at the polar cap, and the pulsar's spin-down power and that the magnetic field in the PWN is positively correlated with the magnetic field at the light cylinder, the potential difference at the polar cap, and the pulsar's spin-down power (see Fig. 4). Other parameters that describe pulsar properties and electron injection have no correlations. These results are consistent with those of Torres et al. (2014).
The results presented in this paper indicate that the spindown power L(t) is correlated with L r , L X , and L bol (see Fig. 5). The results of L(t) versus L r and L(t) versus L X are consistent with those in Torres et al. (2014). Meanwhile, our results show that the L(t) correlates with the ratio L r /L γ and the ratio L X /L γ (see Fig. 6). These correlations are consistent with those given in previous papers (e.g., Mattana et al. 2009;Kargaltsev & Pavlov 2010;Torres et al. 2014).
In this paper, the correlations between L(t) and surface brightness at different bands are presented here (see Fig. 8). Moreover, the correlations between R pwn (T age ) and L X , L X /L γ , and L s are analyzed (see Figs. 9 and 10). For example, the results show that L X /L γ ∼ R −2.54 pwn and L X /L γ ∼ T −3.37 age , which means that old PWNe have smaller values of L X /L γ than do young PWNe. In fact, our results show that the L X /L γ ≥ 10 for young PWNe, 1 < L X /L γ 10 for evolved PWNe, and L X /L γ ≤ 1 for mature/old PWNe. This result may provide a new tool for classifying the evolution state of PWNe. It should be pointed out that the results presented in this paper require confirmation from further observations.