Could the TeV emission of starburst galaxies originate from pulsar wind nebulae?

While the GeV $\gamma$-ray emission of starburst galaxies (SBG) is commonly thought to arise from hadronic interactions between accelerated cosmic rays and interstellar gas, the origin of the TeV $\gamma$-ray emission is more uncertain. One possibility is that a population of pulsar wind nebulae (PWNe) in these galaxies could be responsible for the TeV $\gamma$-ray emission. In this work, we first synthesize a PWNe population in the Milky Way, and assessed their contribution to the $\gamma$-ray emission of the Galaxy, using a time-dependent model to calculate the evolution of the PWN population. Such synthetic PWN population can reproduce the flux distribution of identified PWNe in the Milky Way given a distribution of the initial state of the pulsar population. We then apply it to starburst galaxies and quantitatively calculate the spectral energy distribution of all PWNe in the SBG NGC 253 and M82. We propose that TeV $\gamma$-ray emission in starburst galaxies can be dominated by PWNe for a wide range of parameter space. The energetic argument requires that $\eta_e \times v_{\rm SN}>0.01 {\rm yr}^{-1}$, where $\eta_e$ is the fraction the spin-down energy going to electrons and $v_{\rm SN}$ is the supernova rate. By requiring the synchrotron emission flux of all PWNe in the galaxy not exceeding the hard X-ray measurement of NGC 253, we constrain the initial magnetic field strength of PWNe to be $<400 \mu$G. Future observations at higher energies with LHAASO or next-generation neutrino observatory IceCube-Gen2 will help us to understand better the origin of the TeV $\gamma$-ray emission in SBGs.


INTRODUCTION
TeV -ray have been detected from the starburst galaxies NGC 253 (Abramowski et al. 2012) and M82 (VERITAS Collaboration et al. 2009).Starburst galaxies have a very high star-formation rate (SFR) in a compact central region, the starburst nucleus (SBN).SBGs are endowed with a high interstellar medium (ISM) gas density ( ISM ∼ 100 − 1000 cm −3 ), and magnetic fields of the order of 100 G (Thompson et al. 2006).The enhanced star-forming activity leads to a high supernova (SN) rate  SN ∼ 0.1 − 1 yr −1 , and subsequently a high injection luminosity of cosmic rays (CRs) including hadrons and leptons.CRs could produce -ray from GeV to TeV by interacting with the local ISM and radiation fields (soft photons) (Persic et al. 2008).
Many studies (Peretti et al. 2019;Wang & Fields 2018;Krumholz et al. 2020;Ha et al. 2021) recently modelled the CR spectrum in SBGs based on their multi-wavelength spectra.They found that electrons typically effectively lose energy within SBN, while it is uncertain whether SBNs are calorimeters for protons.It depends on the properties of CR transport in the ISM such as the diffusion coefficient and the speed of the galactic wind driven from the SBN.For example, Peretti et al. (2019) consider the SBN as a spherical compact region with a leaky-box model, in which the injection of On the other hand, the TeV -ray survey of the Galactic plane reveals multiple classes of sources contributing to the total -ray luminosity.This poses a question that how is the fraction of their contribution to the total -ray flux observed from SBGs.Among various TeV -ray emitters in our Galaxy, PWNe are the dominant class.The energies of pulsar wind electrons and positrons range from ∼ 1 GeV to ∼ 1 PeV, placing their synchrotron and inverse Compton (IC) emission into radio-X-ray and GeV-TeV bands, respectively.According to Kargaltsev & Pavlov (2010), about 60 PWNe associated with known radio or -ray pulsars have been detected, 33 of which have measured TeV fluxes.The birthrate of pulsars is correlated with the star formation rate, and thus PWNe should be more abundant in starbursts galaxies.Based on this, Mannheim et al. (2012) argued that a population of individual PWNe could be responsible for the detected TeV emission from SBGs.They used the properties of Galactic TeV-detected PWNe to estimate the contribution of leptonic -ray emission.Considering the impact of starburst environment on PWNe, Ohm & Hinton (2013) revisited this problem and suggested that PWNe can make a significant contribution to the TeV fluxes, provided that the injection spectrum of particles is sufficiently hard and that the average pulsar birth period is short (∼ 35 ms).However, these works assumed that all PWNe are the same, without considering the difference among individual PWNe and the possible parameter distribution among them.
In this work, we revisit the -ray emission contributed by the PWN population in SBGs, considering the time-dependent evolution of the PWN population in SBGs.We will also discuss the influence of model parameters that are not well determined on the TeV emission of SBGs.The rest of the paper is organized as follows: in Section 2, we describe and discuss the assumptions for the PWNe population model; in Section 3, we use the model in Milky Way and compare with observations available; in Section 4, we apply the model to the starburst region and place constraints on the parameters; Then, we give discussions in Section 5 and conclusions in Section 6.We note that the work in this paper appears in earlier form in the conference report by Chen et al. (2023).

THE PWN MODEL
The time-dependent evolution of relativistic electrons and positrons pairs 1 in PWNe can provide us details about the radiation spectrum as a function of their ages.Cooling of these electrons in the magnetic and radiation fields leads to a multi-wavelength spectrum from radio to ray.H.E.S.S. Collaboration et al. (2018a) introduce a time-dependent model for PWNe.It allows us to trace the evolution of the very-highenergy (VHE;  > 100 GeV) electron population, and hence the radiative output of a PWNe, based on a few general assumptions.We adopt this PWN model, then extend it to the evolution of PWN in SBG.Its essential traits are outlined in the following.
PWNe are powered by the rotational energy of related pulsars, which spin down with time and convert a fraction of their rotational energy into non-thermal particles.Over time the pulsar's rotation slows down, and the energy input rate into the nebula decreases.The spin-down luminosity of the pulsar is the rate at which rotational kinetic energy is dissipated, and is thus given by the equation: where Ω is the angular velocity of the pulsar's spin and  is the neutron star's moment of inertia.For a typical mass of 1.4  ⊙ and a radius of 10 km, we have  ∼ 10 45 g cm −2 .The braking index is assumed to be  = 3 for all pulsars in this work, corresponding to a spin-down due to the dipole radiation only.Then pulsar's the surface magnetic field at the equatorial plane  s can be calculated through the observed spin periods P and the time derivatives  by Gaensler & Slane ( 2006) For the energy spectrum of electrons freshly injected into the nebula we assume the following power-law shape: with a power-law index .Φ 0 () can be calculated imposing where   is the electron conversion efficiency, i.e. the fraction of the spin-down energy going to electrons.The energy distribution of the electrons is assumed to be a powerlaw in the energy range from  min to  max .Varying the boundary 1 Hereafter we do not distinguish positrons from electrons for simplicity.
energies essentially changes the number of particles contained in the inverse Compton (IC) relevant energy range, therefore they are also important parameters that will be discussed in this work.A lowenergy break in the injection spectrum is sometimes assumed (e.g.Torres et al. 2014), but it only impacts the lower ends of the emission spectra.We omit it here because it does not influence the TeV flux.
Generally, the magnetic field in the PWNe is a function of both the distance r from the pulsar and the evolution time of the PWNe.If one only cares about the total emission from a PWN, one can simply consider the average magnetic field in the PWNe which changes only with the age of the PWN.Therefore, following Zhang et al. (2008), the magnetic field evolution is given by where  0 is initial magnetic field of PWNe,   is the B-field parameter and  0 =  0 /2  0 is the characteristic timescale, adopting an index of   = 0.6 which satisfies the conservation of magnetic flux.
As far as energy losses are concerned, we include the possibility that particles can leave the nebula as a result of diffusion, as well as the synchrotron energy loss and IC scattering energy loss.The cooling of the electrons during a time step  is implemented in the model by means of an exponential function: This approach uses an effective cooling timescale  eff , which is adopted from Appendix of H.E.S.S.
The iterative evaluation of Eq. 7 then yields the lepton energy distribution as a function of time.After considering the time-dependent evolution of relativistic particles, electrons in PWNe produce synchrotron and IC emission from radio to VHE -ray energies through interactions with magnetic and radiation fields, respectively.The physics of these processes is described in the comprehensive review article by Blumenthal & Gould (1970), which we follow in the implementation of the radiation mechanisms in our model.

MODELING THE PWN POPULATION IN THE MILKY WAY
Modeling of the PWN population starts with the random generation of a pulsar population.For each pulsar, an initial pulsar spin period  0 is sampled from a normal distribution with average value   0 = 50 ms and the standard deviation   0 = 35 ms, truncated at 10 ms.These distributions are similar to the ones used by both Watters &Romani (2011) andJohnston et al. (2020).We use a typical surface magnetic field strength of pulsar at birth, i.e.  log(B s /G) = 12.65 as an average value with a log-normal distribution of the standard deviation  log(B s /G) = 0.55 with  in Gauss (similar to the values used in Martin et al. 2022;Fiori et al. 2022;Watters & Romani 2011).Assuming typical values of 10 45 g cm −2 and 12 km for the neutron star inertia and radius, these properties determine the spindown history of each pulsar.They set the maximum power available at each time for the non-thermal particle injection into PWNe.We apply the PWN population model to the Galaxy and simulate the expected PWNe TeV flux distribution.The expected generation rate of pulsars in the Galaxy is ∼ 1/100 yr (Faucher-Giguère & Kaspi 2006).We can roughly estimate that around ∼ 1000 -ray-emitting PWNe are generated over a period of  end = 10 5 yr.
By varying the final age, we have verified that  end does not affect the results of the present study, given that the collective luminosity of older PWNe is sub-dominant.

The spatial distribution
We are dealing here with only the young pulsars with ages within 100 kyr.We simply assume that pulsars are born in the Galactic plane and do not move significantly over the first ∼ 100 kyr of their lifetime, so we can ignore the pulsar's proper motion (Arzoumanian et al. 2002;Hobbs et al. 2005).The Galactic radial distribution of pulsars is given by Johnston et al. ( 2020) where   () is the density of pulsars (per kpc 2 ) at radius R (in kpc) from the Galactic Centre and   , i, and   are constants with values of 64.6 kpc −2 , 2.35, and 1.258 kpc, respectively.We give the spatial coordinates of each PWN by this radial distribution function, and randomly assign the angular coordinates.We ignore the spiral arm structure of the Milky Way, which will not affect our results.

The TeV flux distribution
According to the electron energy distribution in the Section 2, the emission arising from synchrotron emission and IC scatterings, which are the most important processes, can be obtained.The target photon fields for IC scattering include cosmic microwave background (CMB), starlight, and infrared photons.The uniform CMB component is modelled as a black-body spectrum with an energy density of 0.26 eV cm −3 and temperature of 2.7 K. The starlight and infrared components can be adopted from the GALPROP code by Porter & Strong (2005).The energy densities of the starlight and infrared fields are 1.92 eV cm −3 and 1.19 eV cm −3 , respectively.The temperatures at the spectral peaks are 7906 K for the starlight field component and 107 K for the infrared.
In Figure 1, we compare the expected number distribution of TeV PWNe from the population synthesis to that of the observation.A summary of the parameters used to generate the population is reported in Table 1, except that the index of lepton injection spectrum is form 1.75 to 2.25 randomly for PWN in Milky Way following H.E.S.S. Collaboration et al. (2018a).The electron conversion efficiency   is fixed to 1 for all PWNe based on previous studies.Indeed, Zhu et al. (2018) modeled the dynamical and radiative evolution of 18 PWNe with a 1D leptonic model.Their results indicate 0.93 <   < 0.99 for six young PWNe ( age < 2400 yr) and five evolved PWNe (2500 yr<  age <4600 yr),   ∼ 0.45 for evolved PWNe MSH 15-52, and 0.5 <   < 0.99 for six mature/old PWNe ( age >6500 yr).Similar conclusions are found in other independent modeling of young PWNe, such as   ∈ [0.6, 0.99] by Torres et al. (2014),   ∈ [0.75, 0.95] by Bucciantini et al. (2011), and   approaching to 1 by de Oña Wilhelmi et al. (2022).
We find that the log  − log  distribution of the synthetic population is well consistent with that of the identified PWNe with TeV fluxes above ∼ 10 −12 photons/cm 2 s.Below this value, the completeness of the observed sample drops due to limited sensitivity of instruments and hence the mocked PWN population outnumbers the observed one.The result is compatible to Fiori et al. (2022), Martin et al. (2022) and Cataldo et al. (2020).To put it shortly, our population model provides a satisfactory description of the currently observed PWN population with flux above ∼ 10 −12 photons/cm 2 s.Among all the generated PWNe, young PWNe with age less than 10 kyr contribute more than 90% of the total TeV luminosity of the entire population.In other words, if all these PWNe are generated in an external galaxy (i.e., at roughly same distance from Earth), most the TeV flux emitted by the population would arise from young PWNe.

TEV EMISSION FROM PWN POPULATION IN SBGS
In the following we will investigate whether a population of PWNe can reproduce the TeV emission in SBGs, and explore the parameter space of some key parameters of PWNe in SBGs.

Input parameters
We explore the influences of some main parameters of our model for each individual PWNe and the entire population.
Supernova rate -the number of pulsars is related to the rate of type-II SN explosions:  SN .It determines the number of PWNe that are generated within  end .The estimation of  SN involves various methods, including stellar population fitting, line and FIR emission analysis, radio source modeling, and direct SN searches, but faces challenges due to star formation history, mass function variations, and uncertainties in direct SN searches (Lacki et al. 2011). SN reported in the literature for M82 and NGC 253 span an order of magnitude, from 0.03 yr −1 to 0.3 yr −1 (H.E.S.S. Collaboration et al. 2018b;Behrens et al. 2022).
Conversion efficiency -the fraction the spin-down energy going to electrons:   .The total spin-down power is usually considered to be divided into electrons injected into PWNe (  ), magnetic fields (  ) and other multi-messenger produced elsewhere from the PWNe ( other ) (Gelfand et al. 2009).The following condition is naturally satisfied:   +   +  other = 1.Since   and  SN are linearly proportional to the amount of energy in relativistic electrons, these two parameters are degenerate.

Photon field composition NGC 253 M82
Energy spectrum of electrons -the differential energy spectrum of the injected electrons can be described by a simple power-law function with index , ranging from  min to  max .Martin & Torres (2022) summarized previous studies, and found that PWNe have very similar electrons injection index  ∈ [2.2, 2.8] in Galaxy.We appropriately enlarge this range in the later simulation to better discuss the allowable parameter space.The energy range is relevant for the normalization of the spectrum and determines the ranges of the synchrotron and IC photon spectra.The maximum energy is related to the accelerations processes at the termination shock.Therefore,  max for each PWN is randomly generated with a uniform distribution from 200 TeV to 800 TeV, following the treatment in Martin et al. (2022).Its value would significantly influence the SED of the PWN for  ≤ 2. On the contrary, the minimum  min will affect the resultant SED of the PWN for  > 2.  min is related to the Lorentz factor of the ultra-relativistic pulsar wind (Chevalier 2000).We consider  and  min as two free and independent parameters.
Pulsar population -a pulsar's initial rotational period  0 and the surface magnetic field at the equatorial plane  s , as these determine the temporal evolution of particle injection.The initial period is particularly relevant because it sets the total available spin-down energy,  rot = 2 2 / 2 0 (with  being the pulsar's moment of inertia).We assume the distribution of parameters among new born pulsars in a SBG is similar to that in our Galaxy, and apply the PWN population generator that has been verified in the previous section to SBGs.
Magnetic field -the strength of the initial magnetic field of PWNe,  0 , affects the synchrotron radiation and cooling of electrons.
For a typical interstellar radiation of our Galaxy, young PWNe are expected to be synchrotron-dominated, with a low IC efficiency.However, the interstellar radiation environment in a SBN region is very different from that in our Galaxy.As shown in the Table 2, the infrared radiation dominates the interstellar radiation of the SBN, and its energy density is more than three orders of magnitude higher than the typical value in our Galaxy, implying that the cooling time of highly-relativistic electrons is much shorter.For an electron with energy of 1 TeV, the cooling timescale is approximately 200 yrs in SBN, while it is 3 × 10 5 yrs in ISM of the Milky Way.The IC cooling in such an environment will dominate over synchrotron losses unless  ≳ 500 G.For typical values of  0 ≈ 100 G, as found in some young PWNe (e.g.Mayer et al. 2012; H.E.S.S. Collaboration et al. 2018a), we expect a large fraction of the injected electron energy to be channeled into the IC radiation at GeV -TeV energies.
In Table 1, we summarize the parameters and relative distribution (or the allowed range) used as initial condition for the simulation of the PWNe population.
The target photon fields in SBG include CMB, infrared, and optical radiation fields, where the latter two in SBG are much stronger than those in our Galaxy.Peretti et al. (2019) obtain energy densities and temperatures of the three IR components arising from dust emission and those of the optical radiation from stars.The parameters of these radiation fields for NGC 253 and M82 are listed in Table 2, which are employed in our following calculations.

NGC 253
NGC 253 is one of the only two starburst galaxies found to emit ray from hundreds of MeV (Abdo et al. 2010) to multi-TeV energies (H.E.S.S. Collaboration et al. 2018b) .Based on the planetary nebula luminosity function, a weighted average of the most reliable distance estimates yields a distance of  = 3.5 ± 0.2 Mpc by Rekola et al. (2005).Melo et al. (2002) derive a star-formation rate of ∼ 3.5 ⊙ yr −1 , based on the far-infrared luminosity in the starburst nucleus of NGC 253.Ohm & Hinton (2013)  pulsar birth rate, and generate 5000 pulsars in our sample since we consider emission of PWNe with ages up to 100 kyr.
For each of the simulated pulsars, we assign the values of  0 and  s to them following the description in Section 3. We assume the same values of   and  min for each of the corresponding PWN, and randomly set the value of their injection spectral index  in the range between  min and  min + 1 with an equal probability.To explore the effect of key parameters ,  min and   ×  SN on emission of their PWNe, we test 140 (= 7 × 20) different combinations of values of  min and  min , where  min = [1.5, 2.1] with a linear increment of 0.1 and  min = [1, 1000] GeV with a logarithmic increment of 0.05 dex.The value of   ×  SN affect the fluxes of PWNe linearly and can be determined by matching the observed flux of the galaxy.
In Figure 2, the red line shows the SED produced by a population of PWNe in NGC 253 with the so-called "baseline" model parameters, i.e.  ∈ U (1.8, 2.8) and  min = 162 GeV.Based on the analysis in Section 4.1, we choose   ×  SN = 0.05 yr −1 as an intermediate parameter to analyze.In this baseline model, the IC emission can explain the TeV emission of NGC 253 while the synchrotron emission does not exceed the X-ray limit.In the figure, the NuSTAR upper limit is obtained by subtracting the two known hard X-ray components from the observed flux: thermal gas and X-ray binaries point sources (Wik et al. 2014) .The remaining represents the unresolved, diffuse non-thermal emission from the galaxy.
We investigate the dependence of the gamma-ray flux on different values of  min and  min .While varying one of these two parameter, we fix other parameters the same as those in the baseline model.By summing up the emission of each PWN, we obtain the SED of NGC 253, as shown in the two panels of Figure 3.The top panel of Fig. 3 illustrates the impact of changing  min , from 1 GeV to 1000 GeV.The figure shows that the gamma flux is sensitive to the minimum energy truncation of accelerated electrons: the larger minimum energy truncation, the higher flux.The bottom panel of Fig. 3 illustrates the dependence of the results on the injection spectral index, from U (1.5, 2.5) to U (2.1, 3.1).As is shown, a steeper 10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 Energy (TeV) spectrum leads to a softer gamma-ray spectrum and a smaller gammaray flux above 0.1 TeV.
The parameter space able to account for the TeV flux of NGC 253 is shown in Figure 4, which is obtained by comparing the model and the observed flux at 1 TeV within 1 statistical uncertainties (H.E.S.S. Collaboration et al. 2018b).Combination of  min and  min in the upper left side of the dashed black curve results in a too high TeV flux while that in the lower right side the dashed white curve results in a too low TeV flux.Therefore, the region between the black and white dashed curves are available parameter space.Three panels are obtained with different values of   ×  SN .We see that the available region is broad and covers the typical values that are usually considered for PWNe, as long as the value of   × SN is not too small.When   × SN ≲ 0.015yr −1 , there will be almost no parameter space that can match the observed TeV flux.Given  SN = 0.05 yr −1 as the typical type-II SN rate in the starburst galaxy, it requires PWNe to have a relatively high pair conversion efficiency of   > 0.3 to make a significant contribution to the TeV radiation of the galaxy.Such a high   is actually common for young PWNe ( age < 10 kyr) in Milky Way according to previous studies, as discussed in Section 3.2.Results of these studies support the PWN-origin of the TeV emission from starburst galaxies.The calculated SED has very similar properties to those of NGC 253, as expected due to the similar target radiation field, magnetic fields, and average particle densities in the SB regions (see Figure 5).
The parameter space of PWNe to account for the TeV flux of M82 is given in Figure 6.The available region in the parameter space for M82 is broader compared to that for NGC 253, because VERITAS observation gives a larger statistical error for the TeV flux.When   ×  SN < ∼ 0.01 yr −1 , we find that there will be very limited parameter space that can match the observed TeV radiation, as can be see from Figure 6.

Constraint on initial magnetic field strength of PWNe
Although the initial magnetic field strength of PWNe dose not affect TeV emission of SBG (see Section 4.1), it can affect the intensity of synchrotron radiation, so the X-ray observations can be used to constrain the initial magnetic field of PWNe.Wik et al. (2014) utilize the NuSTAR and Chandra data to investigate the populations contributing to the galaxy-wide 0.5-30 keV emission from NGC 253.Subtracting the contribution of resolved sources and contribution from diffuse gas thermal emission, they determine the 90% upper limit on the non-thermal flux in the 7-20 keV band.Hard X-ray are mainly produced by synchrotron radiation of TeV electrons, according to   = 3 2  /4  .This limit constrains that the initial magnetic field of PWNe cannot be too high, otherwise the synchrotron emission will exceed the upper limit.Fig. 7 shows the SED under different initial magnetic fields of PWNe.Due to that the energy spectrum is sensitive to the electron injection spectrum, we take  ∈ U (1.6, 2.6) as an example.In each panel, the color lines from dark to light indicate the  min increases from 1 MeV to 100 GeV on the logarithmic interval, respectively.10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 Energy (TeV)  This also proves that the minimum energy truncation  min does not affect the hard X-ray flux.
The initial nebular magnetic field strength in PWNe effects the part of non-thermal electrons energy channeled into high-energy photon for by synchrotron radiation.Based on the allowable parameter space (Fig. 4) given by the limitation of the TeV band, when the magnetic field exceeds 400 G quantitatively, the flux in the hard X-ray band will exceed the upper limit given by the X-ray observations, so we obtain an upper limit of  0 < 400G for the initial magnetic field of PWNe in starburst galaxy.

DISCUSSIONS
We here discussion the influence of some subordinate parameters on the results.
The effect of the maximum age for PWNe -The pulsar age is relevant for the particle accumulation, dynamics and energetic.We have verified, however, that varying the final age of PWNe does not affect the results of the present study, given that the -ray luminosity of the associated PWNe powered by older pulsars is negligible.The number of pulsars with age  age is proportional to  age given a constant pulsar birth rate.On the other hand, the spin-down luminosity decreases with  −2 age for a braking index of 3. Therefore, the total spindown luminosity of pulsars of age  age at the present time roughly scales with  −1 age .Note that due to the intense infrared radiation field in starburst galaxies, the cooling timescale of the high energy electron is quite short.Therefore, electrons injected at early stage, when the spin-down luminosity of those middle-aged pulsars were high, cannot survive at the present time.In addition, for middle-aged and old pulsars, the average electron conversion efficiency   is probably only at the level of 0.1, as constrained by the diffusive -ray emission from the Galactic plane (Yan & Liu 2023).Therefore, the contribution of older pulsars is negligible.
Influence of the magnetic field of ISM -The pulsar may escape the parent SNR shell at a time  cross ≃ 45( SN /10 51 erg) 1/3 ( ISM /1 cm −3 ) −1/3 (  /400 km s −1 ) −5/3 kyr (Draine 2011) due to its high kick velocity.The high-energy electrons can largely escape outer of the PWNe and diffuse into the surrounding ISM, producing the so called "pulsar halos", as found recently in the Milky Way (Abeysekara et al. 2017;Aharonian et al. 2021).Unlike the Milky Way, SBN is characterised by a much higher magnetic field and average density ( ISM ∼ 250 G,  ISM ∼ 250 cm −3 ), so the synchrotron losses may be stronger than the IC loss for these high-energy electrons that have escaped the PWNe.Therefore, the contribution to the TeV emission of SBGs by pulsar halos may be subdominant.
Influence of the magnetic field at the beginning stage -As  <  0 (∼ 0.5kyr), the employed evolution of the magnetic field, i.e., Eq. 2 suggests a constant magnetic field at the beginning.At this stage, the PWN expands freely and its radius increase with time as  PWN ∝  6/5 (Gaensler & Slane 2006).If assuming magnetic flux conservation in the PWN (H.E.S.S. Collaboration et al. 2018a), the magnetic field strength would evolve as () ∝  −12/5 .For a given magnetic field as revealed from the SED of a PWN at the present time, it would imply a much stronger magnetic field strength at the beginning.However, even for the early magnetic field predicted by Eq. 2, relativistic electrons responsible for TeV emission cool very rapidly via the synchrotron radiation.Therefore, those electrons injected at early time cannot survive at the present time in either case.For example, in the PWN of SN 1986J ( ∼ 30 yr,  ∼ 10 39 erg s −1 ,  PWN ∼ 17 mG, Tanaka & Kashiyama 2023) , the cooling timescale of relativistic electrons is about 20(/1TeV) −1 days given the inferred magnetic field.Considering a two-segmented magnetic field evolution scenario, it would affect the resulting gamma-ray spectra at most in the order of ∼ 10%, which is negligible compared to parameters discussed in the previous section.
Effect of braking index -While our calculations are based on the assumption of braking index  = 3 for all generated pulsars, some of detected pulsars present different braking indexes, such as 2.5 for the Crab pulsar.The braking index influences the spin-down history.However, taking a different braking index would not affect our result significantly, mainly due to two reasons.First, the property of the simulated pulsar population need to match that of the observed pulsar sample, for example, in terms of distributions of  and .As a result, the distribution of the spin-down luminosity of pulsars, which is determined by − / 3 , is more or less the same, regardless of the chosen value of the braking index.Consequently, the total electron injection luminosity in their PWNe is insensitive to the braking index.On the other hand, the gamma-ray luminosity at the present day depends on all the cumulative electrons injected in the history, which need be traced back over a time period equal to the cooling timescale of emitting electrons   .Given that pulsar's spindown luminosity evolves with time  as the difference, i.e., the ratio, between the electron injection luminosity at the present day ( =  age ) and a period of time .Given that the magnetic field inside the PWN also evolve with time, we estimate the cooling timescale by only considering the IC cooling for simplicity.In the environment of the starburst nucleus, the typical IC cooling timescale for TeV-emitting electron is about   = 200 yr. 0 of most simulated pulsars range in 800 − 5000 yr.We find that the ratio  is around ∼ 0.5−1.As such, changing  = 3 to  = 2 only alter the value of  by a factor less than 2. Indeed, due to the rapid cooling of TeVemitting electrons in the environment of starburst nucleus, electrons can be only accumulated over a short period of time.As a result, the total amount of emitting electrons does not rely on the braking index which controls the injection history.Instead, it is basically determined by the present-day spin-down luminosity of each pulsar, the distribution of which is calibrated by the observed pulsar sample.In Fig. 8, we show the result with different braking index and we find that the difference in the resulting flux is less than 20%.

CONCLUSION
Starburst galaxies such as M82 and NGC 253 show a harder highenergy -ray spectrum with a higher luminosity than that of Milky Way.While it is generally considered that the TeV -ray emission of the starburst galaxies arise from interactions of injected cosmic ray 10 5 10 4 10 3 10 2 10 1 10 0 10 1 10 2 10 3 Energy (TeV) hadrons with interstellar medium, we found that PWNe in starburst regions may also explain the observed TeV flux.
To model the PWN population and subsequently their emission in the starburst galaxies, we firstly simulate the PWN population in Milky Way for verification of our method.We followed the statistical study of properties of pulsars and PWNe in the previous literature and found that the simulated log  − log  distribution of TeV flux of PWNe in Milky Way is consistent with observations.We then applied the method to the starburst galaxies NGC 253 and M82, and calculate the spectral energy distribution produced by the PWN population in the SBGs.Our main results are as follows: (i) PWNe associated with core-collapse supernovae in starburst regions may explain the observed TeV emission of SBG NGC 253 and M82 with typical parameters which are usually employed in previous literature for PWNe.
(ii) From the perspective of energy budget, we found that it generally requires   ×  SN > 0.01 yr −1 for the PWN population and the galaxy.
(iii) By requiring the synchrotron emission of the PWNe not to exceed the hard X-ray observations of SBGs, we constrain the initial magnetic field strength of PWNe to be less than 400 G.
If GeV and TeV emission of SBGs come from the hadronic radiation of CRs in ISM and the leptonic emission of PWNe respectively, we would expect a spectral break somewhere between GeV and TeV energy.Future observations on SBGs by sensitive gamma-ray instruments such as LHAASO (Cao et al. 2019) and the Cherenkov Telescope Array (CTA) (Shimono et al. 2021) is potential to measure such a feature in the spectra of SBGs, which would serve as a critical test of the scenario.Besides, the next-generation neutrino instruments, such as IceCube-Gen2 (Aartsen et al. 2021), could detect neutrinos from SBGs with a long-term exposure if the TeV emission is dominated by hadronic process (Ha et al. 2021).Thus in combination with future neutrino and gamma-ray observations, the origin of leptonic or hadronic emission from SBGs may be distinguished.

Figure 1 .
Figure 1.Flux distributions above 1 TeV of the full population of sources.The synthetic population (in blue) is directly compared with the firmly identified PWNe by Martin et al. (2022) (in black).The coloured areas represent the errors of the curve.

Figure 3 .
Figure 3. Predicted SED for PWNe populations in the starburst region of NGC 253. -LAT points are dots (blue), and H.E.S.S. points are circle (red).Single case flux calculated by changing the value of the two most relevant parameters, one for each panel.

Figure 4 .
Figure 4.For NGC 253, allowable parameter space for  min ,  min and   ×  SN .Deriving parameter ranges by limiting flux at 1 TeV with 1 statistical uncertainty by H.E.S.S.The blue and white dashed lines represent the upper and lower limits of the 1 respectively, and the middle area is the allowable parameter space.

Figure 5 .
Figure 5. Same as 3, but for the SBG M82.Fermi-LAT points are blue crosses (Ackermann et al. 2012), and VERITAS points are red circles (VERITAS Collaboration et al. 2009).

Figure 6 .Figure 7 .
Figure 6.Allowable parameter space for M82, as Figure 4.The blue and white dashed lines represent the upper and lower limits of the 1 by VERITAS respectively, and the middle area is the allowable parameter space.

Figure 8 .
Figure 8.The impact of different braking index on the final results.The solid lines represent braking index n = 3, the dot-dashed lines represent n = 2.5, and the dashed lines represent n = 2.

Table 1 .
Overview of parameter used in the modeling of PWNe.) indicates a uniform distribution from a to b. N ( ,  2 ) indicates a normal distribution with mean value  and standard deviation .

Table 2 .
Peretti et al. (2019)ad and the temperature  of the three IR components due to dust and the optical one due to stars of NGC 253 and M82 fromPeretti et al. (2019).
Abramowski et al. (2012)II SN rate  SN in the starburst nucleus of NGC 253 to be 0.02 yr −1 , whereas H.E.S.S.Collaboration et al. (2018b)suggest an SN rate within the starburst region of NGC 253 of  SN ≈ 0.05 yr −1 .We take it as a reference Predicted SED for PWNe populations in the starburst region of NGC 253 according to baseline model.Also shown is radio data from the Very Large Array(Carilli 1996), the XMM upper limit is emission from the central source X 34 as given byPietsch et al. (2001), the hard X-ray upper limit is taken fromWik et al. (2014), the Fermi-LAT and H.E.S.S. data is fromAbramowski et al. (2012).Horizontal error bars show the energy band over which a particular observation is made, while vertical error bars show 1 uncertainties; the upper limits are given at 95% confidence level.