Constraining Models of the Pulsar Wind Nebula in SNR G0.9+0.1 via Simulation of its Detection Properties using the Cherenkov Telescope Array
MMNRAS , 1–18 (2020) Preprint 1 October 2020 Compiled using MNRAS L A TEX style file v3.0
Constraining Models of the Pulsar Wind Nebula in SNRG0.9+0.1 via Simulation of its Detection Properties usingthe Cherenkov Telescope Array
M. Fiori, , (cid:63) L. Zampieri, A. Burtovoi, , P. Caraveo, and L. Tibaldo Department of Physics and Astronomy, University of Padova, Via F. Marzolo 8, I-35131, Padova, Italy INAF-Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122, Padova, Italy Centre of Studies and Activities for Space (CISAS) ’G. Colombo’, University of Padova, Via Venezia 15, I-35131 Padova, Italy INAF-IASF Milano, Via A. Corti 12, I-20133 Milano, Italy IRAP, Universit´e de Toulouse, CNRS, CNES, UPS, 9 avenue Colonel Roche, 31028 Toulouse, Cedex 4, France
Accepted 29/09/2020
ABSTRACT
SNR G0.9+0.1 is a well known source in the direction of the Galactic Center composedby a Supernova Remnant (SNR) and a Pulsar Wind Nebula (PWN) in the core. Weinvestigate the potential of the future Cherenkov Telescope Array (CTA), simulatingobservations of SNR G0.9+0.1. We studied the spatial and spectral properties of thissource and estimated the systematic errors of these measurements. The source willbe resolved if the VHE emission region is bigger than ∼ . (cid:48) . It will also be possibleto distinguish between different spectral models and calculate the cut-off energy. Thesystematic errors are dominated by the IRF instrumental uncertainties, especially atlow energies. We computed the evolution of a young PWN inside a SNR using a one-zone time-dependent leptonic model. We applied the model to the simulated CTA dataand found that it will be possible to accurately measure the cut-off energy of the γ -rayspectrum. Fitting of the multiwavelength spectrum will allow us to constrain also themagnetization of the PWN. Conversely, a pure power law spectrum would rule out thismodel. Finally, we checked the impact of the spectral shape and the energy densityof the Inter-Stellar Radiation Fields (ISRFs) on the estimate of the parameters of thePWN, finding that they are not significantly affected. Key words: cosmic rays - gamma-rays - pulsars - supernova remnants
Pulsar Wind Nebulae (PWNe) represent the most numerousclass of identified Galactic Very High Energy (VHE) γ -raysources (de O˜na-Wilhelmi et al. 2013). These objects arehighly magnetized nebulae powered by young and energeticpulsars. Inside these nebulae non-thermal radiation up to ∼ TeV is produced (Rieger et al. 2013).In young PWNe the outer radius of the nebula has notyet started to interact with the reverse shock of the SNR.Therefore, they are particularly interesting objects becausethe uncertainties related to the interaction are not presentand their evolution can be fairly well reproduced by phys-ical models. These models can thus be tested against ob-servations, providing important information on the physicalprocesses at work in these sources (e.g. Gelfand et al. 2009;Mart´ın et al. 2012; Zhu et al. 2015). (cid:63)
E-mail: michele.fi[email protected]
The Cherenkov Telescope array (CTA, Actis et al. 2011)will be capable to study the γ -ray emission of PWNe in greatdetail. With CTA it will be possible to observe PWNe fromfew GeV up to hundreds of TeV, accurately sampling mostof the Inverse-Compton (IC) peak as well as obtaining ameasurement of the spectral cut-off energies where present.In addition, the unprecedented angular resolution will allowus to determine more precisely the γ -ray emission regionsand to investigate the existence of any potential energy-dependent morphology. Thanks to this it will be possibleto test various γ -ray emission models of PWNe and to bet-ter understand their magneto-hydrodynamic structure andevolution.The purpose of this work is testing the capabilities ofCTA in connection with a specific source (SNR G0.9+0.1)while, at the same time, assessing the impact of CTA obser-vations on our understanding of the physical processes oc-curring in PWNe. The source selected is SNR G0.9+0.1 (atTeV energies the source is also referred as HESS J1747-281; © a r X i v : . [ a s t r o - ph . H E ] S e p Fiori, Zampieri, Burtovoi, et al.
H. E. S. S. Collaboration et al. 2018a), a well known com-posite Supernova Remnant (SNR, Helfand & Becker 1987).The bright central core has been unambiguously identified asa PWN through X-ray observations (Gaensler et al. 2001).SNR G0.9+0.1 is composed by a PWN in the core (with adiameter of ∼ (cid:48) ) surrounded by a SNR (with a diameter of ∼ (cid:48) , Dubner et al. 2008). This source has been detected atVHE by HESS (Aharonian et al. 2005), VERITAS (Smith& the VERITAS Collaboration 2015) and MAGIC (Ahnenet al. 2017) only up to ∼ TeV, without any evidence ofa cut-off at TeV energies. Moreover, for all these facilities,the source appears point-like because of the limited angularresolution.SNR G0.9+0.1 is considered to be a young PWN withan estimated age of ∼ − years (Camilo et al. 2009;Sidoli et al. 2000). Due to the projected position of thesource, in the direction of the Galactic Center, and the un-certainties in the electron density model in that direction,the distance is not well determined (between 8 and 16 kpc, assuggested by Camilo et al. 2009). SNR G0.9+0.1 has beenoften adopted as a benchmark to test various theoreticalmodels (e.g. Venter & de Jager 2007; Qiao et al. 2009; Fang& Zhang 2010; Tanaka & Takahara 2011; van Rensburg et al.2018; Torres et al. 2014; Zhu et al. 2018). In the early studiesof Venter & de Jager (2007) and Qiao et al. (2009), only anapproximate treatment of the energy losses was included,while the dynamical evolution of the nebula was not con-sidered. Fang & Zhang (2010) incorporated the dynamicalevolution of the nebula, but assumed an injection spectrumfor the electrons in the form of a Maxwellian plus a power-law tail, instead of the most widely adopted broken powerlaw (as in Tanaka & Takahara 2011, Torres et al. 2014, andZhu et al. 2018). More recently, van Rensburg et al. (2018)presented a more accurate multi-zone time-dependent lep-tonic model to reproduce the spatial properties of the source.In this paper, we did not focus on modelling in detail theenergy-dependent morphology of SNR G0.9+0.1 (the angu-lar resolution at VHE is not sufficient to do it), but adopteda one-zone time-dependent leptonic model, even if it hasbeen shown that lower energy observations with a better an-gular resolution would benefit from multi-zone models (seee.g. van Rensburg et al. 2018; Lu et al. 2019; van Rens-burg et al. 2020). Following Torres et al. (2014) and Zhuet al. (2018) we considered the evolution of a single popu-lation of accelerated electrons inside an expanding uniformmedium in spherical symmetry. This approach turned out tobe sufficiently accurate for reproducing the multiwavelength(MWL) emission of the PWN and allowed us to make pre-dictions on the spectrum of SNR G0.9+0.1 at the highestenergies.Similarly, we used SNR G0.9+0.1 as a test case todemonstrate the improvements that the CTA South arraywill allow us to achieve. The source position, its faintness(only about of the Crab flux) and the small angular sizemake this object a really interesting target for testing thecapabilities of the CTA. Since the extension of the PWN inSNR G0.9+0.1 is comparable to the best angular resolutionachievable with CTA, we expect to be ale to measure its sizeat VHEs. A measurement of the angular size of the source isneeded to better constrain the physical models and to com-pare the source size at different wavelengths. This would helpunderstanding if the VHE emission comes from the central source or if there is some contribution from the SNR shell.In addition, the sensitivity of CTA will be much better upto and above 100 TeV (CTA Consortium 2019), allowing usto measure a possible cut-off at energies higher than 20 TeV(not excluded with the currently available data). Also thismeasurement is important to better constrain the physicalmodels of the nebula, since it will constrain the particle in-jection spectrum, and specifically the maximum energy ofthe electrons (assuming a leptonic model). At such high en-ergies, the inverse Compton emission may be in the Klein-Nishina regime, and thus obtaining such a measurement willbe a very good proxy of the actual maximal electron energy.This in turn may constrain the acceleration process at thePWN termination shock.In this work we present a comprehensive study of thespatial and spectral properties of SNR G0.9+0.1 aiming attesting the observability of specific features in the simulateddata, studying the spatial extension of the TeV emission andthe presence of a VHE cut-off in the spectrum, and compar-ing the data to models of the MWL spectrum. Furthermore,we estimate the systematic uncertainties that may affect ob-servations of SNR G0.9+0.1 carried out with CTA.This paper is organized as follows. In Section 2 we de-scribe the models and the analysis of the spatial and spectralproperties of SNR G0.9+0.1 as seen by CTA. In Section 3we report the results of the simulations. In Section 4 weestimate the systematic uncertainties and in Section 5 wediscuss the results of our analysis. In Section 6 we describethe implementation of a physical model for the emission ofa young PWN inside a SNR. Finally in Section 7 we discussour results and compare the numerical solutions with thesimulations of the CTA observations of SNR G0.9+0.1. To simulate, reduce and analyze the γ -ray data we made useof the software ctools, a software package developed for thescientific analysis of CTA data (Kn¨odlseder et al. 2016).We specified in input: a spatial and a spectral modeldescribing the emission region of SNR G0.9+0.1 and a modelfor the spatial distribution of the cosmic-ray background.For the spectral models we adopted both a power law and apower law with an exponential cut-off (PLEC): dNdE = N (cid:18) EE (cid:19) − Γ , (1) dNdE = N (cid:18) EE (cid:19) − Γ exp (cid:18) EE cut (cid:19) , (2)where N is a normalization factor, Γ the spectral index, E the pivot energy and E cut the cut-off energy. For thespatial model, we use different distributions as described inthe following.SNR G0.9+0.1 is projected in the direction of thecrowded region of the Galactic Center. In order to under-stand which sources can significantly affect the measurementof the flux of SNR G0.9+0.1 and to test the capability ofctools in reproducing the extended emission of the GalacticCenter, we simulate the γ -ray emission in a field of ° x ° around the position of Sgr A*. In doing that, we take into MNRAS000
H. E. S. S. Collaboration et al. 2018a), a well known com-posite Supernova Remnant (SNR, Helfand & Becker 1987).The bright central core has been unambiguously identified asa PWN through X-ray observations (Gaensler et al. 2001).SNR G0.9+0.1 is composed by a PWN in the core (with adiameter of ∼ (cid:48) ) surrounded by a SNR (with a diameter of ∼ (cid:48) , Dubner et al. 2008). This source has been detected atVHE by HESS (Aharonian et al. 2005), VERITAS (Smith& the VERITAS Collaboration 2015) and MAGIC (Ahnenet al. 2017) only up to ∼ TeV, without any evidence ofa cut-off at TeV energies. Moreover, for all these facilities,the source appears point-like because of the limited angularresolution.SNR G0.9+0.1 is considered to be a young PWN withan estimated age of ∼ − years (Camilo et al. 2009;Sidoli et al. 2000). Due to the projected position of thesource, in the direction of the Galactic Center, and the un-certainties in the electron density model in that direction,the distance is not well determined (between 8 and 16 kpc, assuggested by Camilo et al. 2009). SNR G0.9+0.1 has beenoften adopted as a benchmark to test various theoreticalmodels (e.g. Venter & de Jager 2007; Qiao et al. 2009; Fang& Zhang 2010; Tanaka & Takahara 2011; van Rensburg et al.2018; Torres et al. 2014; Zhu et al. 2018). In the early studiesof Venter & de Jager (2007) and Qiao et al. (2009), only anapproximate treatment of the energy losses was included,while the dynamical evolution of the nebula was not con-sidered. Fang & Zhang (2010) incorporated the dynamicalevolution of the nebula, but assumed an injection spectrumfor the electrons in the form of a Maxwellian plus a power-law tail, instead of the most widely adopted broken powerlaw (as in Tanaka & Takahara 2011, Torres et al. 2014, andZhu et al. 2018). More recently, van Rensburg et al. (2018)presented a more accurate multi-zone time-dependent lep-tonic model to reproduce the spatial properties of the source.In this paper, we did not focus on modelling in detail theenergy-dependent morphology of SNR G0.9+0.1 (the angu-lar resolution at VHE is not sufficient to do it), but adopteda one-zone time-dependent leptonic model, even if it hasbeen shown that lower energy observations with a better an-gular resolution would benefit from multi-zone models (seee.g. van Rensburg et al. 2018; Lu et al. 2019; van Rens-burg et al. 2020). Following Torres et al. (2014) and Zhuet al. (2018) we considered the evolution of a single popu-lation of accelerated electrons inside an expanding uniformmedium in spherical symmetry. This approach turned out tobe sufficiently accurate for reproducing the multiwavelength(MWL) emission of the PWN and allowed us to make pre-dictions on the spectrum of SNR G0.9+0.1 at the highestenergies.Similarly, we used SNR G0.9+0.1 as a test case todemonstrate the improvements that the CTA South arraywill allow us to achieve. The source position, its faintness(only about of the Crab flux) and the small angular sizemake this object a really interesting target for testing thecapabilities of the CTA. Since the extension of the PWN inSNR G0.9+0.1 is comparable to the best angular resolutionachievable with CTA, we expect to be ale to measure its sizeat VHEs. A measurement of the angular size of the source isneeded to better constrain the physical models and to com-pare the source size at different wavelengths. This would helpunderstanding if the VHE emission comes from the central source or if there is some contribution from the SNR shell.In addition, the sensitivity of CTA will be much better upto and above 100 TeV (CTA Consortium 2019), allowing usto measure a possible cut-off at energies higher than 20 TeV(not excluded with the currently available data). Also thismeasurement is important to better constrain the physicalmodels of the nebula, since it will constrain the particle in-jection spectrum, and specifically the maximum energy ofthe electrons (assuming a leptonic model). At such high en-ergies, the inverse Compton emission may be in the Klein-Nishina regime, and thus obtaining such a measurement willbe a very good proxy of the actual maximal electron energy.This in turn may constrain the acceleration process at thePWN termination shock.In this work we present a comprehensive study of thespatial and spectral properties of SNR G0.9+0.1 aiming attesting the observability of specific features in the simulateddata, studying the spatial extension of the TeV emission andthe presence of a VHE cut-off in the spectrum, and compar-ing the data to models of the MWL spectrum. Furthermore,we estimate the systematic uncertainties that may affect ob-servations of SNR G0.9+0.1 carried out with CTA.This paper is organized as follows. In Section 2 we de-scribe the models and the analysis of the spatial and spectralproperties of SNR G0.9+0.1 as seen by CTA. In Section 3we report the results of the simulations. In Section 4 weestimate the systematic uncertainties and in Section 5 wediscuss the results of our analysis. In Section 6 we describethe implementation of a physical model for the emission ofa young PWN inside a SNR. Finally in Section 7 we discussour results and compare the numerical solutions with thesimulations of the CTA observations of SNR G0.9+0.1. To simulate, reduce and analyze the γ -ray data we made useof the software ctools, a software package developed for thescientific analysis of CTA data (Kn¨odlseder et al. 2016).We specified in input: a spatial and a spectral modeldescribing the emission region of SNR G0.9+0.1 and a modelfor the spatial distribution of the cosmic-ray background.For the spectral models we adopted both a power law and apower law with an exponential cut-off (PLEC): dNdE = N (cid:18) EE (cid:19) − Γ , (1) dNdE = N (cid:18) EE (cid:19) − Γ exp (cid:18) EE cut (cid:19) , (2)where N is a normalization factor, Γ the spectral index, E the pivot energy and E cut the cut-off energy. For thespatial model, we use different distributions as described inthe following.SNR G0.9+0.1 is projected in the direction of thecrowded region of the Galactic Center. In order to under-stand which sources can significantly affect the measurementof the flux of SNR G0.9+0.1 and to test the capability ofctools in reproducing the extended emission of the GalacticCenter, we simulate the γ -ray emission in a field of ° x ° around the position of Sgr A*. In doing that, we take into MNRAS000 , 1–18 (2020) onstraining Models of the PWN in SNR G0.9+0.1 with CTA account all the known TeV sources and the diffuse emissionin the direction of the Galactic Center, as outlined below. In a box of 3 square degrees around the center of the Galaxy,there are many sources at TeV energies as observed by theHESS, MAGIC and VERITAS collaborations (Aharonianet al. 2006a; Archer et al. 2016; H. E. S. S. Collaborationet al. 2017; Ahnen et al. 2017).We consistently selected all the sources from the HESScatalogue , except for SNR G0.9+0.1 for which we consid-ered all the data included in a joint HESS + VERITAS anal-ysis of the source (Smith & the VERITAS Collaboration2015). The sources considered in our simulation are listedbelow and their spatial and spectral parameters are reportedin Table 1. • HESS J1745-290 (Aharonian et al. 2004): This sourcerepresents the TeV emission coming from the center ofour Galaxy (Acero et al. 2010). It is associated with thesuper-massive black hole Sgr A* or to the candidate PWNG359.95-0.04 (Kistler 2015). It is modelled as a point sourcewith a power law spectrum with an exponential cut-off. Thespectral parameters are taken from Aharonian et al. (2009) • HESS J1741-302 (Tibolla et al. 2008): It is an uniden-tified source detected with HESS at ∼ of the Crab fluxabove 1 TeV. We modelled it as a point source with a powerlaw spectrum. • HESS J1745-303 (Aharonian et al. 2006b): This is anextended and unidentified VHE γ -ray source at a Galacticlongitude of − . ° . The morphology of the source is quitecomplex owing to the presence of 3 major emitting regions.The spatial extension of this source has been modelled usingthe HESS excess map , shown in Figure 1. The spectralmodel is a power law (Aharonian et al. 2008). • Diffuse emission along the Galactic plane (Aharonianet al. 2006a): It is a region of diffuse emission (of approxi-mately ± ° in galactic longitude) probably associated withthe interaction of cosmic-ray particles with molecular cloudsand that contains a number of unidentified sources such asfor example HESS J1746-285 (H. E. S. S. Collaboration et al.2017). This diffuse emission is the only source that can affectour simulation of SNR G0.9+0.1 because the spatial emis-sion regions of these sources overlap. For the spatial modelwe used a section (between . ° < l < . ° and | b | < . ° , inGalactic coordinates) of an image taken from HESS (Figure2) in which the emission coming from HESS J1745-290 andSNR G0.9+0.1 has been previously subtracted. The spectralmodel is a power law. • SNR G0.9+0.1 (Aharonian et al. 2005): The spatialmodel is taken from a radio map at 843 MHz from the veritas.sao.arizona.edu/ The results of the analysis on the sole HESS data (Aharonianet al. 2005) are consistent with the results of the joint analysis. Figure 1.
Excess map of HESS J1745-303 (Aharonian et al.2008), used as spatial model for our simulations.The three dashedcircles indicate the positions of the brightest emitting regions ofthe source.
Figure 2.
HESS excess map of the diffuse emission around theGalactic center (the emission from SNR G0.9+0.1 and HESSJ1745-209 has been previously subtracted) (Aharonian et al.2006a), used as input spatial model for our simulations.
Sydney University Molonglo Sky Survey (radiomap tem-plate hereafter). The map has been prepared for the simula-tion with a technique developed for the analysis of extendedsources in Fermi -LAT . For the spectral model we used asingle power-law for the entire system, as assumed for theHESS and VERITAS observations (Smith & the VERITASCollaboration 2015) since from currently available data itis not possible to discriminate between the emission comingfrom the PWN and the SNR.To simulate observations of the field with the southernCTA facility (CTA-South) we made use of the InstrumentResponse Functions (IRFs) of the baseline array made avail-able by the CTA Consortium (Acharyya et al. 2019). Weprovided as input all the information on the sources listed skyview.gsfc.nasa.gov/surveys/sumss/mosaics/Galactic/J1752M28.FITS fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/extended/extended.html MNRAS , 1–18 (2020)
Fiori, Zampieri, Burtovoi, et al.
Table 1.
Input data used with ctobssim to simulate the VHE emission from a region of ≈ ° x ° around the Galactic center. The reportedpositions are taken from the SIMBAD astronomical database (Wenger et al. 2000), except for the position of the Galactic diffuse emissionfor which we adopt the position of the center of the template map. E is always equal to 1 TeV.Source Spatial Model Position Spectral Model Input ParametersHESS J1745-290 point source RA= ◦· PLEC d N = . × − TeV − cm − s − Dec= − ◦· Γ = . E cut = . TeVHESS J1741-302 point source RA= ◦· power law e N = . × − TeV − cm − s − Dec= − ◦· Γ = . HESS J1745-303 Extended Source RA= ◦· power law f N = . × − TeV − cm − s − HESS excess map a Dec= − ◦· Γ = . Galactic Diffuse Extended Source RA= ◦· power law g N = . × − TeV − cm − s − sr − HESS excess map b Dec= − ◦· Γ = . SNR G0.9+0.1 Extended Source RA= ◦· power law h N = . × − TeV − cm − s − SUMSS radio map (843 MHz) c Dec= − ◦· Γ = . a Figure 1, b Figure 2, c Figure 3 d Aharonian et al. (2009), e Tibolla et al. (2008), f Aharonian et al. (2008), g Aharonian et al. (2006a), h Smith & the VERITAS Collaboration (2015)
Figure 3.
Radio image (843 MHz) of SNR G0.9+0.1 taken fromthe Sydney University Molonglo Sky Survey (SUMSS) and usedas template for some of our simulations. Most of the power in theradio band is coming from the PWN that is surrounded by theless energetic shell of the supernova remnant. above, plus a model for the spatial distribution of the cosmic-ray background (”CTAIrfBackground”). We simulated fourobservations with different observing times centered on theposition of Sgr A*, in the energy range . − TeV: one30 minute observation, one 5 hour observation, one 50 hoursobservation and one 200 hours observation. We simulatedobservations lasting up to 200 hours because we wanted totest the results achievable with CTA under the best assump-tions regarding the observing time. This number is justifiedby the fact that the Galactic Center will be extensively ob-served during the first years of CTA operations (CTA Con- sortium 2019). We then made an unbinned analysis andfitted all the simulated data with the same models given ininput. Applying the maximum likelihood method, we finallycompute the Test Statistics (TS) value for each source . As far as SNR G0.9+0.1 is concerned, we divided the anal-ysis in two parts: first we fixed all the spectral parametersof the source and varied only the spatial model, then wekept fixed the spatial model (one of the previously selectedmodels) and varied the spectral parameters. At this stagewe include in the simulations only the information on SNRG0.9+0.1, the diffuse emission from the Galactic plane andthe cosmic-ray background. The simulated field has a radiusof . ° centered on the source.To understand the capabilities of CTA in resolving thespatial extension of the VHE emission of SNR G0.9+0.1,we perform the simulations using different spatial models inthe energy range 0.2-180 TeV. All the simulated observationslast 200 hours and have fixed spectral parameters (a powerlaw with the parameters reported in Table 1). The spatialmodels used here are: point source (assuming that the VHEemission comes only from the inner part of the remnant), aradio map template (assuming that the VHE emission comesfrom the same region as the radio emission) and various spa-tially uniform radial disk models with different radii, from 1arcsec to 90 arcsec. We then fit all the simulated data withfour different spatial models: a point source model, a ra-dial Gaussian model, a radial disk model, and the radiomaptemplate model. Model fitting has been performed with abinned maximum likelihood analysis (binned cube cen-tered on source position with . ° pixel size bin, 2500 pixel, http://cta.irap.omp.eu/ctools/users/tutorials/quickstart/unbinned.html The square root of the TS value is roughly the Gaussian σ inthe case of one free parameter associated to the source (see e.g.Protassov et al. 2002) http://cta.irap.omp.eu/ctools/users/tutorials/quickstart/fitting.html MNRAS000
Radio image (843 MHz) of SNR G0.9+0.1 taken fromthe Sydney University Molonglo Sky Survey (SUMSS) and usedas template for some of our simulations. Most of the power in theradio band is coming from the PWN that is surrounded by theless energetic shell of the supernova remnant. above, plus a model for the spatial distribution of the cosmic-ray background (”CTAIrfBackground”). We simulated fourobservations with different observing times centered on theposition of Sgr A*, in the energy range . − TeV: one30 minute observation, one 5 hour observation, one 50 hoursobservation and one 200 hours observation. We simulatedobservations lasting up to 200 hours because we wanted totest the results achievable with CTA under the best assump-tions regarding the observing time. This number is justifiedby the fact that the Galactic Center will be extensively ob-served during the first years of CTA operations (CTA Con- sortium 2019). We then made an unbinned analysis andfitted all the simulated data with the same models given ininput. Applying the maximum likelihood method, we finallycompute the Test Statistics (TS) value for each source . As far as SNR G0.9+0.1 is concerned, we divided the anal-ysis in two parts: first we fixed all the spectral parametersof the source and varied only the spatial model, then wekept fixed the spatial model (one of the previously selectedmodels) and varied the spectral parameters. At this stagewe include in the simulations only the information on SNRG0.9+0.1, the diffuse emission from the Galactic plane andthe cosmic-ray background. The simulated field has a radiusof . ° centered on the source.To understand the capabilities of CTA in resolving thespatial extension of the VHE emission of SNR G0.9+0.1,we perform the simulations using different spatial models inthe energy range 0.2-180 TeV. All the simulated observationslast 200 hours and have fixed spectral parameters (a powerlaw with the parameters reported in Table 1). The spatialmodels used here are: point source (assuming that the VHEemission comes only from the inner part of the remnant), aradio map template (assuming that the VHE emission comesfrom the same region as the radio emission) and various spa-tially uniform radial disk models with different radii, from 1arcsec to 90 arcsec. We then fit all the simulated data withfour different spatial models: a point source model, a ra-dial Gaussian model, a radial disk model, and the radiomaptemplate model. Model fitting has been performed with abinned maximum likelihood analysis (binned cube cen-tered on source position with . ° pixel size bin, 2500 pixel, http://cta.irap.omp.eu/ctools/users/tutorials/quickstart/unbinned.html The square root of the TS value is roughly the Gaussian σ inthe case of one free parameter associated to the source (see e.g.Protassov et al. 2002) http://cta.irap.omp.eu/ctools/users/tutorials/quickstart/fitting.html MNRAS000 , 1–18 (2020) onstraining Models of the PWN in SNR G0.9+0.1 with CTA gnomonic projection, and 100 logarithmic energy bins). Atthis stage, we adopted the binned analysis because, for longexposures, the computation time is much shorter than withthe unbinned analysis.After the analysis of the spatial properties of the source,we perform the analysis of the spectral properties fixing allthe spatial parameters. Our goal is to asses the detectabilityof the source in the higher energy range (from 30 TeV upto 180 TeV) and the capability of CTA-South to distinguishbetween different spectral models. We simulate different ob-servations, all lasting 200 hours, with the source spatiallymodelled with the radiomap template and spectrally mod-elled with a power law and various PLEC with different cut-off energies (20 TeV, 30 TeV, 50 TeV, and 100 TeV). Dataare simulated in the energy range between 0.2 TeV and 180TeV. Model fitting has been performed with the binned like-lihood analysis. The spectral energy distribution (SED) ofthe source is extracted using csspec , a specific tool of ctools. In Table 2 we show the results of the unbinned analysisperformed on the four different simulations of the Galacticcenter region mentioned in section 2.1. We report all theTS values and the spectral parameters measured for all thesources in the simulations. These measurements were per-formed to check the detectability of all the simulated sourcesand to determine the needed observing time to reliably re-cover all the parameters of the sources. After 30 minutes ofobservation, all the sources are significantly detected and,as expected, the significance grows increasing the observingtime. Already at 50 hours the inferred parameters are ingood agreement with the input ones. At 200 hours, the in-ferred parameters are very close to the simulated ones andthe associated errors become very small. Therefore, a 200hours observation would lead to the accuracy on the mea-sured parameters of SNR G0.9+0.1 needed for the analysisreported below.We then compared the simulation obtained for an expo-sure of 50 hours with that obtained with HESS in 55 hours (Aharonian et al. 2006a) in a similar energy range (see Fig-ure 4). The images are in good agreement, with the CTAsimulated one having a lower background contamination.With the same observing time, CTA will allow us to ob-tain a wider spectral coverage and a higher signal-to-noiseratio. We performed two different analyses to investigate the re-solving capabilities of CTA. In the first analysis, we carryout different fits of the image simulated using the radiomaptemplate. The fits were performed with four different spa-tial models: point source model, spatially uniform radial disk Figure 4.
Simulation of the Galactic center extended emissionas seen with CTA in an observation of 50 hours in two differentenergy ranges ( . − TeV top panel, − TeV middle panel)and a residual map of the same region from an HESS observationin the energy range ∼ . − TeV (lower panel; Aharonian et al.2006a) . model (with the radius left free during the fit), radial Gaus-sian model (with width left free during the fit) and the ra-diomap model. The results are shown in Figure 5. If theVHE emission follows the radio emission CTA could be ableto detect the source as an extended object because the TSvalue for the point source fit is significantly lower. The ex-tended models have similar TS values, with the radiomaptemplate being slightly more significant, indicating that allthe three models can reproduce well the simulated data andthat the VHE γ -ray emission from outside the PWN (i.e.the emission coming from the SNR shell that can be seen inFigure 3) is almost negligible.In the second analysis, we test the limiting resolvingcapabilities of CTA against the background of the GalacticCenter VHE emission region following the procedure devel-oped to detect an extended source in the Fermi -LAT data(Lande et al. 2012). We have simulated different images as-suming a spatially uniform radial disc with different radii.We then fit all the images with a point source model and aradial disc model with the radius free to vary. This procedureis then repeated for 100 times to account for the statistical
MNRAS , 1–18 (2020)
Fiori, Zampieri, Burtovoi, et al.
Table 2.
Results of the unbinned maximum likelihood analysis on the simulated observations of the Galactic center region. After 30minutes of observation all the sources are significantly detected.Source 0.5 hour observation 5 hour observationSpectral parameters a,b
TS Spectral parameters a,b
TSHESS J1745-290 N = ( . ± . ) × − N = ( . ± . ) × − Γ = . ± . Γ = . ± . E cut = . ± . E cut = . ± . HESS J1741-302 N = ( . ± . ) × − N = ( . ± . ) × − Γ = . ± . Γ = . ± . HESS J1745-303 N = ( . ± . ) × − N = ( . ± . ) × − Γ = . ± . Γ = . ± . Galactic Diffuse N = ( . ± . ) × − sr − N = ( . ± . ) × − sr − Γ = . ± . Γ = . ± . SNR G0.9+0.1 N = ( . ± . ) × − N = ( . ± . ) × − Γ = . ± . Γ = . ± .
50 hour observation 200 hour observationSpectral parameters a,b
TS Spectral parameters a,b
TSHESS J1745-290 N = ( . ± . ) × − N = ( . ± . ) × − Γ = . ± . Γ = . ± . E cut = . ± . E cut = . ± . HESS J1741-302 N = ( . ± . ) × − N = ( . ± . ) × − Γ = . ± . Γ = . ± . HESS J1745-303 N = ( . ± . ) × − N = ( . ± . ) × − Γ = . ± . Γ = . ± . Galactic Diffuse N = ( . ± . ) × − sr − N = ( . ± . ) × − sr − Γ = . ± . Γ = . ± . SNR G0.9+0.1 N = ( . ± . ) × − N = ( . ± . ) × − Γ = . ± . Γ = . ± . a N in unit of TeV − cm − s − and E cut in unit of TeV. E = TeV b Statistical error only fluctuations that can arise from different simulations .Forall the simulated images we compute the significance of de-tecting significant spatial extension for the source by usingthe likelihood ratio test:TS ext = − L RD L PS , (3)where L RD and L PS are the likelihood values of the fits withthe radial disc (RD) and the point source (PS) models. InFigure 6 we show TS ext in function of the simulated sourceradius with the confidence level errors. The value in-creases from very small to large radii, showing that the radialdisc model has a significantly better likelihood (TS ext ≥ )when the source has a radius larger than + − arcsec. Thismeans that if the VHE emission region of SNR G0.9+0.1 isbigger than ∼ . arcmin, and the response of the instru-ment is very well known, CTA will be able to detect it as anextended source even if the PSF of the instrument is larger( ∼ . arcmin). However, it would be difficult to study sub-structures inside the source because the angular size of thesesubstructures would be too small. Different simulations are based on a different random seed forthe Monte Carlo generator that samples the input source modelsto produce observed photon energies and arrival directions. Thisis achieved through the random number generator provided in theGammaLib library(Kn¨odlseder et al. 2016).
Figure 5.
Test Statistics values for different fitting models ap-plied to the simulation in which the VHE emitting region of SNRG0.9+0.1 is modelled with the radio map template. The TS forthe point source fitting model has a lower significance comparedto the other fitting models.
As far as the CTA spectrum of SNR G0.9+0.1 is con-cerned, it is shown in Figure 7. The present analysis aims atunderstanding how well it is possible to recover the expected
MNRAS000
MNRAS000 , 1–18 (2020) onstraining Models of the PWN in SNR G0.9+0.1 with CTA Figure 6.
Significance of the detection of the source as an ex-tended source for the images simulated using a spatially uniformradial disc with different radii. We consider TS ext ≥ as theminimum value for claiming that the source is extended. Theradial disc model has a significantly better likelihood when thesource has a radius larger than + − arcsec. The red area showsthe 95% confidence level error region. Figure 7.
CTA spectral energy distribution of SNR G0.9+0.1,simulated using different cut-off energies. It is clearly possible todistinguish the different spectra. cut-off of this source. This has strong implications for thephysical modeling implemented in section 6 since a differ-ent cut-off energy could lead to inferring different physicalparameters for the nebula. The spectrum has a good statis-tics and therefore the spectral resolution is very good. It isclearly possible to distinguish spectra with different cut-offenergies. This represents a significant improvement in com-parison with currently available data that does not allow todistinguish if the spectral shape of the VHE emission is apower law or a power law with a cut-off at energies higherthan 20 TeV (Figure 8).The maximum cut-off energy detectable in the CTAsimulated spectrum is (cid:38)
TeV while, for the lowest energycut-off considered here (20 TeV), the source is detectableonly up to ∼ TeV.
Figure 8.
Comparison of the spectra simulated in this work withthe data from HESS (black square) and VERITAS (black stars).With the current available data it is not possible to rule out mod-els with cut-off energies higher than 20 TeV.
The spectral analysis of the simulated data returns only thestatistical errors, computed from the covariance matrix ofthe maximum likelihood fitting procedure. But systematicerrors need to be carefully accounted for in order to assessthe accuracy of the results. A fit of simulated data withoutconsidering the systematic errors will lead to overestimat-ing the goodness of the fit and to results that may not berealistic.We considered both the instrumental sources of uncer-tainties and the background related uncertainties. The in-strumental sources of uncertainties are due to the imperfectknowledge of the effective area and the accuracy of the re-constructed energy scale, while the background sources ofuncertainties are due to the cosmic-rays and the Galacticdiffuse emission. In order to translate uncertainties into sys-tematic errors on fluxes and spectral indices we will makesome assumptions on how these uncertainties propagates.In the case of the instrumental uncertainties we startfrom CTA technical requirements and we apply the followingprocedure to measure the associated errors. • Knowledge of the effective area .Uncertainty on the effective area of the system must be < (from the CTA technical requirements). To estimate theeffect of such an uncertainty we followed the method used bythe Fermi-LAT team (Ackermann et al. 2012, Sec. 5.7). Wegenerate perturbed IRFs that represent the worst scenario,extract the spectral parameters and compare them to thoseobtained with the unperturbed IRF. The perturbed effectivearea A (cid:48) ef f is written as: A (cid:48) ef f ( E , θ ) = A ef f ( E , θ ) · ( + ξ Ae f f · B ( E )) , (4)where A ef f is the unperturbed effective area, ξ Ae f f = . isthe uncertainty and B ( E ) a function of the energy (bracket-ing function). Different form for B ( E ) are adopted dependingon the spectral parameter considered. For a simple powerlaw, to maximize the effect on the normalization, the func-tion B ( E ) is written as: B ( E ) = ± , (5)while, to maximize the effects on the spectral index, the MNRAS , 1–18 (2020)
Fiori, Zampieri, Burtovoi, et al. following expression is used: B ( E ) = ± tanh (cid:18) .
13 log (cid:18) EE (cid:19)(cid:19) , (6)where E is the same pivot energy used in Equation 1 and 2.With these two modified IRFs we have reanalyzed the dataand estimated the errors on the spectral parameters fromthe values obtained in the two cases. • Accuracy of the Energy Scale.
The uncertainty on the energy of a photon event candidatemust be < (from the CTA technical requirements). Inorder to estimate the errors on the spectral parameters in-duced by this uncertainty , we took the simulated data andperturbed all the photon energies as: E (cid:48) = E · ( ± ξ Escale ) , (7)where ξ Escale = . . We have then analyzed these data andestimated the errors on the spectral parameters.In the case of the uncertainties related to the knowl-edge of the background we applied a different approach, asdescribed below. • Cosmic-ray Background.
In order to determinate the impact of the uncertainty on thecosmic-ray background we varied its flux of ± from thenominal value. We thus changed the normalization of thebackground according to: N (cid:48) = N · ( ± ξ CRbkg ) , (8)where ξ CRbkg = . . We then analyzed these data and esti-mated the errors on the spectral parameters. Since the devi-ations from the nominal values resulting from this source ofuncertainty seem to be negligible, as discussed in the nextsection, it was not worth considering variations induced bychanges in the photon index of the cosmic-ray background. • Galactic Diffuse Emission.
As mentioned above we modelled the emission from theGalactic plane central region using an HESS observation(Aharonian et al. 2006a). The best fit spectral model forthis observation is a power law with N = . × − TeV − cm − s − sr − and Γ = . with the associated errors σ N = (± . stat ± . syst ) × − TeV − cm − s − sr − and σ Γ = ± . stat ± . syst . Using these errors we calculatean optimistic/pessimistic spectrum from the Galactic cen-ter from: F pess , opt = F ( E ) ± (cid:115)(cid:18) ∂ F ∂ N (cid:19) σ N + (cid:18) ∂ F ∂ Γ (cid:19) σ Γ , (9)where F ( E ) is the best fit value of the flux, the pessimisticcase F pess corresponds to the sign + and the optimistic case F opt to the sign − . This is an approximation of the errorpropagation formula (we lack all the information on the fullcovariance matrix that comes from the analysis made bythe HESS collaboration). The spectrum is shown in Figure9. We have then analyzed these perturbed data and mea-sured the associated errors. We repeated the analysis using In this work we have not taken into account the energy disper-sion since it was computationally too expensive. This value is much bigger then the expected uncertainty onthe residual cosmic-ray background for CTA-South.
Figure 9.
Average spectrum (per steradian) of the Galactic planenear the center region (between . ° < l < . ° and | b | < . ° )as measured by the HESS collaboration (Aharonian et al. 2006a).The shaded area correspond to the error boundary of the HESSmeasurements, in which both statistical and systematic errors aretaken in account. the pessimistic and optimistic estimate of the spectrum andused the spectral parameters of the source inferred in thetwo cases to estimate the errors induced by this systematicuncertainty on it. It is worth to mention that also the un-certainty on the morphology of the Galactic diffuse emissioncan be a source of systematics error. However, at presentwe have not enough information to assess the uncertaintiesrelated to the morphology of diffuse emission. This task isleft for future studies.For all these sources of uncertainty we have repeatedthe simulations one hundred times and we have then takenthe final errors on the average values as representative of theuncertainties induced by the different simulations. In Table 3 we report the values of the systematic errors,computed from the difference between the ”Nominal value”(values computed without perturbing the data) and the val-ues obtained as explained in the previous section.The instrumental systematic uncertainties dominateover the background related sources of error. This is shownin Figure 10 where we plot the errors as a function of energy,assuming a power law spectrum. While the systematics actdifferently at different energies, the background related un-certainties are always small. In the low energy range (wherethe array has the best sensitivity) the instrumental uncer-tainties dominate and are at the same level as the statisticalerrors, while in the higher energy range, the decrease of thesensitivity of CTA-South leads to an increase of the statisti-cal errors. The behavior of the statistical error yields a goodrepresentation of the sensitivity limit of the CTA-South ar-ray. Although the errors reported here are probably over-estimated (especially the instrumental ones), this analysisprovides a good clue on the order of magnitude of the ex-pected systematic uncertainties. According to the results ofour analysis, the background related uncertainties are neg-
MNRAS000
MNRAS000 , 1–18 (2020) onstraining Models of the PWN in SNR G0.9+0.1 with CTA ligible in comparison with the other sources of uncertaintyand have a small impact on the measured spectrum. PWNe are important laboratories to test the processes re-sponsible for the acceleration of charged particles. To thisend, it is crucial to compare real or simulated data withprecise and physically motivated models.Reproducing the broad-band spectrum, from the radioband up to γ -rays, of these sources, requires a dynamicalmodel that describes the evolution of the population of theemitting electrons inside the PWN.A one-zone time-dependent leptonic model is oftenadopted. In this model the main emitting particles are a pop-ulation of electrons that evolves with time and the nebula isapproximated as a sphere where the electrons are uniformlydistributed.This approach has been developed by several authors(e.g. Venter & de Jager 2007; Qiao et al. 2009; Zhang et al.2008; Fang & Zhang 2010; Tanaka & Takahara 2010, 2011;Bucciantini et al. 2011; Mart´ın et al. 2012, 2016; Torres et al.2014; van Rensburg et al. 2018). In this work, we follow theapproach presented by Gelfand et al. (2009). We also test theresult of our implementation for the PWN in SNR G0.9+0.1against those obtained by Zhu et al. (2018) and Torres et al.(2014). The distribution and the evolution of the electronic popu-lation inside the nebula is described by an energy-diffusionequation. The general form of this equation (see equationA1) and the meaning of all the terms of the equation aredescribed in Appendix A. The simplified form used in thiswork is as follow: ∂ N ( E , t ) ∂ t = Q ( E , t ) − ∂∂ E [ b ( E ) N ( E , t )] − N ( E , t ) τ esc ( E , t ) . (10)where N ( E , t ) is the number density of the electrons, Q ( E , t ) the injection rate of electrons at the terminationshock, b ( E ) the variation of the mean energy of the elec-trons per unit time, and τ esc ( E , t ) is a characteristic timescale describing the escape of the electrons from the system.The typical shape adopted for the injection spectrumof the particles is a broken power law. Other types of in-jection spectra have been proposed but all somehow fail toreproduce the observed spectrum or are difficult to motivate(see the discussion in Gelfand 2017). The broken power lawspectrum can reproduce well the different slopes of the syn-chrotron spectrum observed in many PWNe, as the CrabNebula (Atoyan & Aharonian 1996), in the radio and X-raybands. We then assume (Tanaka & Takahara 2010; Buc-ciantini et al. 2011; Mart´ın et al. 2012; Torres et al. 2013;Zhu et al. 2015): Q ( E , t ) = Q ( t ) (cid:40) ( E / E b ) − α for E ≤ E b ( E / E b ) − α for E > E b , (11)where Q ( t ) is a normalization factor determined from the fraction of the spin-down luminosity L ( t ) of the pulsar thatgoes in particles energy and E b is the break energy wherethe slope of the particle spectrum changes. If we write thespin-down luminosity of the pulsar in the form L ( t ) = L (cid:18) + t τ (cid:19) − ( n + )( n − ) , (12)where L is the initial spin-down luminosity, τ the initialspin-down timescale and n the braking index (Gaensler &Slane 2006), we can find the normalization factor Q ( t ) from: ( − η B ) L ( t ) = ∫ E max E min EQ ( E , t ) dE . (13)Here the constant η B , the magnetic fraction of the nebula(Mart´ın et al. 2016), is the fraction of the spin-down lumi-nosity that goes into the electromagnetic field . While − η B is the fraction of the spin-down luminosity that goes in thekinetic energy of the electrons.The escape timescale τ esc ( E , t ) is computed from the as-sumption that particles can escape from the nebula becauseof diffusion. This diffusion inside a PWN arises from the in-teraction of particles with irregularities in the magnetic field(Vorster et al. 2013). Assuming that the diffusion of plasmaacross the magnetic field in the PWN follows Bohm law, τ esc ( E , t ) is given by: τ esc = eB ( t ) R pwn ( t ) Ec , (14)where R pwn is the radius of the PWN.The second term in equation 10 includes the energy vari-ation because of synchrotron radiation, IC scattering, Self-Synchrotron Compton (SSC) and adiabatic losses (Ginzburg& Syrovatskii 1964).The minimum energy E min of the injected electrons is afree parameter in this model and we choose to select a valueequal to the electrons rest mass energy ( . MeV). On theother hand the maximum electron energy E max has to bedetermined because it is strictly related to the accelerationsprocesses at the termination shock. There are different waysto calculate E max . For high magnetic field strengths (forvery young PWNe) one can estimate it by balancing syn-chrotron losses acceleration gains (de Jager et al. 1996). Forlower magnetic field strengths, one needs to consider thatthe highest energy particle must have a gyro-radius compa-rable to the shock radius to participate to the accelerationprocess (de Jager & Djannati-Ata¨ı 2009). Another possibil-ity for estimating E max is to consider the electric potential ofthe neutron star magnetosphere (Bandiera 2008; Bucciantiniet al. 2011; Granot et al. 2017) and determine the maximumenergy that electrons can gain while moving through the po-lar cap potential. We computed E max considering all threedifferent approaches and adopted the second one becausethe other two yield unreasonably high values. The secondcondition is equivalent to impose that the Larmor radius R L must be a fraction (cid:15) < ( (cid:15) containment factor) of the termi-nation shock radius R S . The Larmor radius can be written This is not to be confused with the so called magnetizationparameter σ ( t ) = η B /( − η B ) .MNRAS , 1–18 (2020) Fiori, Zampieri, Burtovoi, et al.
Table 3.
Systematic errors measured using the deviation of the perturbed values from the nominal ones, as explain in the text. Wereport for comparison also the statistical errors computed from the likelihood analysis made with ctools.
Statistical errors N a δ N a δ N / N Γ δ Γ δ Γ / Γ Nominal value . × − . × − . .
306 0 .
006 0 . Systematic errors N (cid:48) a δ N (cid:48) a δ N (cid:48) / N Γ (cid:48) δ Γ (cid:48) δ Γ (cid:48) / Γ A (cid:48) e f f (Eq. 5 -5%) . × − . × − . .
309 0 .
003 0 . A (cid:48) e f f (Eq. 5 +5%) . × − − . × − − . .
310 0 .
004 0 . A (cid:48) e f f (Eq. 6 +5%) . × − . × − . .
333 0 .
027 1 . En. scale ( − ) . × − − . × − − . .
320 0 .
014 0 . En. scale ( + ) . × − . × − . . − . − . Cosmic-ray ( − ) . × − . × − . . − . − . Cosmic-ray ( + ) . × − . × − . .
309 0 .
003 0 . Gal. Diffuse (Opt.) . × − − . × − − . . − . − . Gal. Diffuse (Pess.) . × − − . × − − . .
308 0 .
002 0 . a N and δ N in unit of TeV − cm − s − Figure 10.
Fractional error on the CTA spectrum as function of photon energy, measured assuming a power law model for SNRG0.9+0.1. In the low energy range (where CTA-South will have the best sensitivity) the instrumental uncertainties are dominant, whilein the higher energy range the decrease of the sensitivity of the array leads to an increase of the statistical errors. The background relateduncertainties are very low at all energies. MNRAS000
Fractional error on the CTA spectrum as function of photon energy, measured assuming a power law model for SNRG0.9+0.1. In the low energy range (where CTA-South will have the best sensitivity) the instrumental uncertainties are dominant, whilein the higher energy range the decrease of the sensitivity of the array leads to an increase of the statistical errors. The background relateduncertainties are very low at all energies. MNRAS000 , 1–18 (2020) onstraining Models of the PWN in SNR G0.9+0.1 with CTA as R L = E max eB S , (15)and so the maximum energy becomes: E max = (cid:15) eB S R S . (16)Finally we need an expression for the magnetic field at thetermination shock B S . From Kennel & Coroniti (1984) thepost-shock field is expressed as: B S = κ (cid:114) η B L ( t ) c R S , (17)where κ is the magnetic field compression ratio taken equalto 3 (strong shock condition). The final expression for themaximum electron energy is then: E max = e (cid:15) (cid:114) η B L ( t ) c . (18)To compute the evolution of the magnetic field we con-sider the adiabatic losses due to expansion work done bythe nebula on the surroundings and the energy input fromthe pulsar wind (Pacini & Salvati 1973; Torres et al. 2013;Gelfand et al. 2009): dW B ( t ) dt = η B L ( t ) − W B ( t ) R pwn ( t ) dR pwn dt , (19)where W B = ( π / ) R pwn ( t ) B ( t )/( π ) is the total magneticenergy. The integration over time of this equation leads to B ( t ) = R pwn ( t ) (cid:115) η B ∫ t L ( t (cid:48) ) R pwn ( t (cid:48) ) dt (cid:48) . (20)The last ingredient of the model is the dynamical evo-lution (radius and the expansion velocity) of the PWN. Wecompute it with an iterative approach that is explained inappendix B .The diffusion-loss equation (equation 10) is solved us-ing a freely available code, called GAMERA (Hahn 2015).Once the evolution of the particle spectrum is computed, it ispossible to derive directly the photon spectrum with GAM-ERA. The synchrotron spectrum is computed considering anisotropic pitch angle distribution of the electrons as in Ghis-ellini et al. (1988). The IC emission is computed using thefull Klein-Nishina cross-section (Blumenthal & Gould 1970)on a background radiation field (generally composed by theCMB photons and two Infra-Red components). SynchrotronSelf-Compton (SSC) emission is also included (Atoyan &Aharonian 1996). The model has several parameters that constrain variousphysical properties of the system. Since some of them aresignificantly degenerate, as the distance and the age of thesystem, we decide to fix them by choosing reliable value asreported in the literature (age, distance, energy of the SNexplosion, density of the interstellar medium and photon The caveats of this iterative approach are described at the endof appendix B. libgamera.github.io/GAMERA/docs/main_page.html background, see Table 4). In addition to these parameters,several parameters of the pulsar (spin-down luminosity, pe-riod derivative, characteristic age) are also known and arereported in Table 4. The remaining parameters are those re-lated to the spectrum of the injected electrons population(the break energy and the two indices of the broken powerlaw), the magnetic fraction of the nebula and the contain-ment factor.When fitting the data we leave the injection parame-ters free to vary. The only exceptions are α and E b thatcan be constrained from the radio and X-ray data. As al-ready stated, changing some of the fixed parameters could,in principle, lead to very different values for the fitted pa-rameters. For example, changing the distance of the systemwould lead to different values for the ejected mass of theSN, the age of the system and the densities of the back-ground photon fields for preserving the radius and TeV flux.This would in turn lead to estimating completely differentparameters for the nebula. The distance of the source mustbe estimated accurately to break this degeneracy. However,once the distance is fixed at a certain value, the fitted pa-rameters are fairly well determined. In the following we willnot consider this degeneracy and we will fix the distance ofthe source to 13.3 kpc (as reported by H. E. S. S. Collabora-tion et al. 2018b), since determining it is not the main focusof this paper. The fitting procedure and error estimation ofthe fitted parameters are reported in Appendix C.We tested our implementation against the results pre-sented in Zhu et al. (2018) and Torres et al. (2014), select-ing the same set of data for consistency. The radio data aretaken from Dubner et al. (2008), the X-ray data from Por-quet et al. (2003) and the current VHE data from Aharonianet al. (2005). For the X-ray data, in performing the fit weconsidered only two points, one at the lower and the otherat the higher bound of the energy interval (with the corre-sponding errors). They were computed from the best-fittingpower law reported by Porquet et al. (2003). The rationalbehind this choice was to avoid giving too much weight tothe X-ray data in comparison with the radio data (with onlythree points) , to sample with a similar number of points thesynchrotron and the IC peaks (5 and 7 points, respectively),and to comparatively increase the weight of the TeV data inthe following section. This is crucial to understand to whatextent the better quality of the CTA data will help in esti-mating the parameters of PWNe.The values of the fitted parameters and their compari-son with those found in Zhu et al. (2018) and Torres et al.(2014) are reported in Table 4. Results are consistent. How-ever, the (fixed) value of the ejected mass M ej is slightlydifferent. This difference is likely caused by differences inthe approach adopted to solve Equation 10. However, thediscrepancy does not appear to be particularly relevant con-sidering the actual uncertainty on the knowledge of this pa-rameter.We emphasize that the parameter (cid:15) is loosely con-strained because the data do not cover the part of the spec-trum where the effects of this parameter are more evident(i.e. in the high energy tails of the synchrontron and ICpeaks). Is possible to see this effect in Figure 11 where wevary only (cid:15) between . and . with a constant step of . . This parameter is only constrained to be > . . We MNRAS , 1–18 (2020) Fiori, Zampieri, Burtovoi, et al.
Figure 11.
SED models for SNR G0.9+0.1 computed for α = . and η B = . fixed while (cid:15) is varying from 0.02 to 0.98 with step0.04. We can clearly see that is possible to rule out only reallysmall values of the containment factor ( (cid:15) (cid:46) . ). Figure 12.
Photon (upper panel) and electron (lower panel) SEDmodel for SNR G0.9+0.1 computed with the best fit values re-ported in Table 4. The dataset is the same as in the work of Zhuet al. (2018) (more detail in the text). then took (cid:15) = . as reference value for all the models inthe subsequent analysis.Figure 12 shows the final best fit electron and photonspectra. The reduced chi square of the fit is χ ν = . . The reported value of the reduced chi square is not to be in-tended as an absolute measurement of the goodness of the fit onthe original complete data set (we did not consider all the X-rayspectral points), but only as a reference value useful for compar-ison with the fits of the simulated data reported below.
We applied this model of the PWN evolution to the vari-ous simulated spectra of SNR G0.9+0.1 reported in Section3.2, assuming that the most of the simulated VHE emissioncomes from the central PWN. The spectral range is limitedat 200 GeV to be consistent with the lower limit of the HESSdata and we rebin the spectrum with 10 bins. However, de-pending on the spectral shape, the spectrum can have lessthan 10 bins since at high energies there may be no photons.The results of the model fit are reported in Table 5,while two representative MWL spectra are shown in Figure13. The errors on the γ -ray data-points includes also thesystematic errors computed in the previous section.From Figure 13 we see that in the lower energy part ofthe spectrum (the synchrotron emission peak) the model isalways consistent with the data, while this is not the case athigh energies.The most interesting results is that the value of themagnetization parameters η B is fairly well determined andtend to decrease with increasing cut-off energy, because, forenergy conservation, increasing the maximum energy of theelectrons requires that more power goes in particles ( − η B ) and less in the magnetic field. In general, the MWLspectrum, can constrain it.For a cut-off at − TeV we found a good agreementof the fitted parameters with the values obtained from theHESS data. For a cut-off at a different energy the inferredparameters have significantly different values, which meansthat with the data currently available it is not possible toaccurately constrain them. With the CTA data, which havea higher energy threshold, the estimates will be more accu-rate. The increased sensitivity of CTA will then allow us toobserve this and other PWNe at higher energies and makeaccurate studies on how particles are accelerated at the ter-mination shock.Finally we want to emphasize that the model spectra arenot consistent with a pure power law simulated spectrum forevery value of the parameters (reduced χ ν (cid:39) . ). With thismodel we are not able to reproduce a power law with nomeasured cut-off. Even changing the age and distance of thesource, it is not possible to find a model that has a powerlaw tail up to 180 TeV. The only possibility would probablybe including an hadronic component, but this is beyond thepurpose of this work. We now try to estimate the impact on our results caused bythe uncertainties on the Inter-Stellar Radiation Field (ISRF)at the (unknown) position of SNR G0.9+0.1. In principle adifferent ISRF can affect our measurement of the parameterof the nebula since the shape of the IC component is depen-dent on the background radiation. In the previous analysiswe fixed the parameters of the ISRF. It would have beencomputationally too expensive to let them free.The density and temperature of the Near-Infrared(NIR) and Far-Infrared (FIR) photon field can vary signifi-cantly with the position in the galaxy. Moreover the spectralshape of this emission can be very different from the sim-ple sum of diluted black-bodies (as assumed in the previoussections).
MNRAS000
MNRAS000 , 1–18 (2020) onstraining Models of the PWN in SNR G0.9+0.1 with CTA Table 4.
Fixed and fitted parameters of the model in comparison with those of Zhu et al. (2018) and (Torres et al. 2014, Model 2). Dotsmean that the same value is adopted. All the parameters are computed for the estimated age t a .This work Zhu et al. (2018) Torres et al. (2014) Notes Pulsar and SN parameters (fixed) P [ms] . ... ... from Camilo et al. (2009) (cid:219) P [s s − ] . × − ... ... from Camilo et al. (2009) τ c [yr] ... ... P /( n − ) (cid:219) Pn ... ... fixed at the standard braking index value L ( t a ) [erg/s] . × ... ... from Camilo et al. (2009) t a [yr] ... ... estimated age a τ [yr] ... ... [ τ c /( n − )] − t a L [erg/s] . × ... ... from equation 12 M e j [M (cid:12) ] estimated a E sn [erg] ... ... estimated a d [kpc] . ... . from H. E. S. S. Collaboration et al. (2018b) Environment parameters (fixed) n h [cm ] . ... . from Zhu et al. (2018) T C M B [K] . ... ... from Longair (2008) U C M B [eV/cm ] . ... ... from Longair (2008) T F I R [K] ... ... from Torres et al. (2014) U F I R [eV/cm ] . ... ... from Torres et al. (2014) T N I R [K] ... ... from Torres et al. (2014) U N I R [eV/cm ] ... ... from Torres et al. (2014) Injection parameters E b [TeV] . ... . from Zhu et al. (2018) α . ... . from Zhu et al. (2018) α . ± .
022 2 . ± .
02 2 . fitted η B . ± . . ± .
004 0 . fitted (cid:15) > .
10 0 . ± .
08 0 . fitted PWN parameters b R pwn ( t a ) [pc] . ± .
01 3 .
51 3 . from iterative procedure in Appendix B B ( t a ) [ µ G] . + . − . . + . − . from equation 20 E max ( t a ) c [TeV] >
600 1452 + − from equation 18 a t a , M ej and E SN taken in order to obtain a nebula of ∼ (cid:48) located at 13.3 kpc. b Computed from the PWN dynamics (see Appendix B) c Maximum energy of the electrons in injection at the termination shock of the nebula.
Figure 13.
Photon SED computed with the best-fitting parameters for two out of the six different CTA simulated spectra of SNRG0.9+0.1 reported in Table 5.MNRAS , 1–18 (2020) Fiori, Zampieri, Burtovoi, et al.
Table 5.
Results of the fitting procedure with the PWN model adopted in this work for various CTA simulated spectra of SNR G0.9+0.1.PLEC 20 TeV PLEC 30 TeV PLEC 50 TeV PLEC 80 TeV PLEC 100 TeV PWL
Fitting result χ ν . . . . . . α . ± .
017 2 . ± .
016 2 . ± .
016 2 . ± .
016 2 . ± .
016 2 . ± . η B . ± . . ± . . ± . . ± . . ± . . ± . R pwn (pc) . ± .
01 3 . ± .
01 3 . ± .
01 3 . ± .
01 3 . ± .
01 3 . ± . B ( µ G) . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . In order to estimate the effects of different ISRFs, weperform two different approaches. In the first, we check howmuch the fit differs comparing the case with fixed and freeISRF parameters. To do it we cannot use the full modelsince the computational time would be too large. We thentreated the dynamical evolution in a simplified way, assum-ing a PWN freely expanding in the SNR using just equationB4. We then considered the CTA simulated data with a cut-off at 30 TeV and fitted them leaving α and η B free. Weused a Monte Carlo Markov Chain (MCMC) code ( emcee ,Foreman-Mackey et al. 2013) and made 2500 realizationsof the spectrum. We obtain results similar to those pre-viously found ( α = . + . − . , η B = . + . − . ). Afterthis, we repeated the fit but adding as free parameters theenergy density and temperature for the IR radiation fields( T FI R , U FI R , T N I R , U N I R ). We found in this case a differentISFR, with an higher energy density of the Far IR compo-nent (see Figure 14). However the relevant parameters ofthe PWN did not change significantly, although their errorsincreased ( α = . + . − . , η B = . + . − . ).In the second approach we considered a more realisticradiation field, like the axisymmetric solution for the ISFRof the Milky Way provided by Popescu et al. (2017), anduse it to produce a model with fixed nebula parameters( α = . , η B = . ). We selected the model reportedin the first panel in Figure 9 of Popescu et al. (2017) andrescaled it by a factor ∼ to obtain a similar γ -ray flux as theone of SNR G0.9+0.1. We then used this model to simulatean observation made with CTA, extracted the new spectrumand used it in the MCMC fitting procedure as before. We fitthe usual two parameters α and η B fixing again the valuesfor the ISRF as in the previous analysis and using two di-luted blackbodies to model it. We obtained values that arein very good agreement with the ones used for the prepa-ration of this model ( α = . + . − . , η B = . + . − . ).The results are shown in Figure 15. We also tried to fitthis model leaving all the parameters for the IR radiationfield free to vary and found similar values. While the energydensity of the ISRF is crucial to reproduce the IC compo-nent in the VHE spectrum, its actual spectral distribution isnot, because Comptonized IR photons tends to loose rapidlymemory of their initial energy. In this work we have studied a young PWN inside SNRG0.9+0.1 that is projected near the Galactic Center. De-spite the high background rate, the crowded field and thefaintness of the source we have shown that the CTA-South array enables us to study this region, and in particular thePWN, in great detail.In our analysis of SNR G0.9+0.1 we choose 200 hoursas observing time for the simulations in order to obtain veryaccurate data. This observing time is early achievable be-cause of the projected position of this source, close to theGalactic Center. As reported in the book ”Science with theCherenkov Telescope Array” (CTA Consortium 2019), theGalactic Center is one of the Key Science Project (KSP)for the CTA collaboration. This core programme will runfor the first 3 years of observations with CTA and will pro-duce 525 hours of data from the region of SNR G0.9+0.1.The 200 hours of time needed for our study will be reachedafter ∼ year after the beginning of the observations withCTA-South.Our spatial analysis of SNR G0.9+0.1 shows that if theVHE emission region is larger than ∼ . arcmin CTA willbe able to resolve it, leading to a measurement of the size ofthe nebula in the VHE band. Furthermore our spectral anal-ysis shows that it would be possible to distinguish differentspectral models and calculate the cut-off energy, if present.We could also detect the source at energies higher then 100TeV if the spectrum is a pure power law.We performed also a detailed analysis of the system-atic errors and found that the systematics related to theinstrumental uncertainties dominate, especially at low ener-gies. Despite these errors maybe somewhat overestimated,they provide at least an order of magnitude estimate of theuncertainties that is crucial for our subsequent analysis.We have then implemented a one-zone time dependentleptonic model that computes the evolution of a young PWNinside a SNR in order to obtain some physical informationand to understand what impact on our knowledge of thisPWN CTA may have. We first compared our result withthose obtained by Zhu et al. (2018) and Torres et al. (2014)using the same data-set. We find good agreement, althoughit is difficult to constraint the confinement factor (cid:15) (hencethe maximum energy of the injected electrons in the nebula).Measurements of the flux of SNR G0.9+0.1 at MeV energieswould be needed to obtain a precise value for this parame-ter. However, in the absence of MeV data, an increase of theVHE observing time would help to put constraints on themaximum electron energy because the tail of the IC peak isalso sensitive to it at high energies. From the best fit modelof the currently available data (Table 4) we expect an highenergy cut-off between 20 and 30 TeV. This is a measure-ment that CTA could easily do, as shown in Figure 7, thusallowing us to reduce the uncertainties on the estimated pa-rameters of the PWN (see Table 5).It is worth nothing that the possibility to put a con- MNRAS000
01 3 . ± . B ( µ G) . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . In order to estimate the effects of different ISRFs, weperform two different approaches. In the first, we check howmuch the fit differs comparing the case with fixed and freeISRF parameters. To do it we cannot use the full modelsince the computational time would be too large. We thentreated the dynamical evolution in a simplified way, assum-ing a PWN freely expanding in the SNR using just equationB4. We then considered the CTA simulated data with a cut-off at 30 TeV and fitted them leaving α and η B free. Weused a Monte Carlo Markov Chain (MCMC) code ( emcee ,Foreman-Mackey et al. 2013) and made 2500 realizationsof the spectrum. We obtain results similar to those pre-viously found ( α = . + . − . , η B = . + . − . ). Afterthis, we repeated the fit but adding as free parameters theenergy density and temperature for the IR radiation fields( T FI R , U FI R , T N I R , U N I R ). We found in this case a differentISFR, with an higher energy density of the Far IR compo-nent (see Figure 14). However the relevant parameters ofthe PWN did not change significantly, although their errorsincreased ( α = . + . − . , η B = . + . − . ).In the second approach we considered a more realisticradiation field, like the axisymmetric solution for the ISFRof the Milky Way provided by Popescu et al. (2017), anduse it to produce a model with fixed nebula parameters( α = . , η B = . ). We selected the model reportedin the first panel in Figure 9 of Popescu et al. (2017) andrescaled it by a factor ∼ to obtain a similar γ -ray flux as theone of SNR G0.9+0.1. We then used this model to simulatean observation made with CTA, extracted the new spectrumand used it in the MCMC fitting procedure as before. We fitthe usual two parameters α and η B fixing again the valuesfor the ISRF as in the previous analysis and using two di-luted blackbodies to model it. We obtained values that arein very good agreement with the ones used for the prepa-ration of this model ( α = . + . − . , η B = . + . − . ).The results are shown in Figure 15. We also tried to fitthis model leaving all the parameters for the IR radiationfield free to vary and found similar values. While the energydensity of the ISRF is crucial to reproduce the IC compo-nent in the VHE spectrum, its actual spectral distribution isnot, because Comptonized IR photons tends to loose rapidlymemory of their initial energy. In this work we have studied a young PWN inside SNRG0.9+0.1 that is projected near the Galactic Center. De-spite the high background rate, the crowded field and thefaintness of the source we have shown that the CTA-South array enables us to study this region, and in particular thePWN, in great detail.In our analysis of SNR G0.9+0.1 we choose 200 hoursas observing time for the simulations in order to obtain veryaccurate data. This observing time is early achievable be-cause of the projected position of this source, close to theGalactic Center. As reported in the book ”Science with theCherenkov Telescope Array” (CTA Consortium 2019), theGalactic Center is one of the Key Science Project (KSP)for the CTA collaboration. This core programme will runfor the first 3 years of observations with CTA and will pro-duce 525 hours of data from the region of SNR G0.9+0.1.The 200 hours of time needed for our study will be reachedafter ∼ year after the beginning of the observations withCTA-South.Our spatial analysis of SNR G0.9+0.1 shows that if theVHE emission region is larger than ∼ . arcmin CTA willbe able to resolve it, leading to a measurement of the size ofthe nebula in the VHE band. Furthermore our spectral anal-ysis shows that it would be possible to distinguish differentspectral models and calculate the cut-off energy, if present.We could also detect the source at energies higher then 100TeV if the spectrum is a pure power law.We performed also a detailed analysis of the system-atic errors and found that the systematics related to theinstrumental uncertainties dominate, especially at low ener-gies. Despite these errors maybe somewhat overestimated,they provide at least an order of magnitude estimate of theuncertainties that is crucial for our subsequent analysis.We have then implemented a one-zone time dependentleptonic model that computes the evolution of a young PWNinside a SNR in order to obtain some physical informationand to understand what impact on our knowledge of thisPWN CTA may have. We first compared our result withthose obtained by Zhu et al. (2018) and Torres et al. (2014)using the same data-set. We find good agreement, althoughit is difficult to constraint the confinement factor (cid:15) (hencethe maximum energy of the injected electrons in the nebula).Measurements of the flux of SNR G0.9+0.1 at MeV energieswould be needed to obtain a precise value for this parame-ter. However, in the absence of MeV data, an increase of theVHE observing time would help to put constraints on themaximum electron energy because the tail of the IC peak isalso sensitive to it at high energies. From the best fit modelof the currently available data (Table 4) we expect an highenergy cut-off between 20 and 30 TeV. This is a measure-ment that CTA could easily do, as shown in Figure 7, thusallowing us to reduce the uncertainties on the estimated pa-rameters of the PWN (see Table 5).It is worth nothing that the possibility to put a con- MNRAS000 , 1–18 (2020) onstraining Models of the PWN in SNR G0.9+0.1 with CTA Figure 14.
Photon SED computed with the fixed ISRF background (left panel) and with the free ISRF background (right panel). Thebest fit value are computed with an MCMC procedure.
Figure 15.
In the left panel we show the realistic (Popescu et al. 2017) IR background reprocessed via IC used as input for the simulation.In the right panel the photon SED fitted with just two diluted blackbody. straint on the size of the VHE emission region with CTAwould be crucial to check the goodness of the model, be-cause we could compare it with the model computed radiusand to the size observed at other wavelength.We have shown that MWL data, including CTA data(that will be capable to constrain the cut-off, if present), willlead to a more precise measurement of the magnetization pa-rameter η B of the nebula, that, for simplicity, we consideredto be constant in space and time during the evolution of thenebula. We note also that, with this model, it is not possibleto reproduce a pure power law spectrum. If detected withCTA, this would require a drastic change in the model, suchas the introduction of an hadronic component.We checked also the effects caused by uncertainties onthe ISRF field. A fit leaving the ISRF parameters free leadsonly to small differences in the values of α and η B . Evenapproximating a realistic ISRF with only two diluted black-bodies, the values of α and η B are not significantly affected. ACKNOWLEDGEMENTS
We thank the referee for his useful comments. We wouldalso like to thanks to developer of the software GAMERA,J. Hahn, for the availability at answering very quickly andclearly at our questions on the usage of the code and also formaking available it to all. This paper has gone through inter-nal review by the CTA Consortium. We made use of ctools,a community-developed analysis package for Imaging Air Cherenkov Telescope data. ctools is based on GammaLib,a community-developed toolbox for the high-level analysisof astronomical gamma-ray data. This research made usealso of the following PYTHON packages: MATPLOTLIB(Hunter 2007), NUMPY (van der Walt et al. 2011), AS-TROPY (Astropy Collaboration et al. 2013) and EMCEE(Foreman-Mackey et al. 2013). We acknowledge financialcontribution from INAF through grant ”ASTRI/CTA DataChallenge (ACDC).
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author.
REFERENCES
Acero F., et al., 2010, MNRAS, 402, 1877Acharyya A., et al., 2019, Astroparticle Physics, 111, 35Ackermann M., et al., 2012, ApJS, 203, 4Actis M., et al., 2011, Experimental Astronomy, 32, 193Aharonian F., et al., 2004, A&A, 425, L13Aharonian F., et al., 2005, A&A, 432, L25Aharonian F., et al., 2006a, Nature, 439, 695Aharonian F., et al., 2006b, ApJ, 636, 777Aharonian F., et al., 2008, A&A, 483, 509Aharonian F., et al., 2009, A&A, 503, 817Ahnen M. L., et al., 2017, A&A, 601, A33MNRAS , 1–18 (2020) Fiori, Zampieri, Burtovoi, et al.
Archer A., et al., 2016, ApJ, 821, 129Astropy Collaboration et al., 2013, A&A, 558, A33Atoyan A. M., Aharonian F. A., 1996, MNRAS, 278, 525Bandiera R., 2008, A&A, 490, L3Blondin J. M., Chevalier R. A., Frierson D. M., 2001, ApJ, 563,806Blumenthal G. R., Gould R. J., 1970, Reviews of Modern Physics,42, 237Bucciantini N., Arons J., Amato E., 2011, MNRAS, 410, 381CTA Consortium 2019, Science with the Cherenkov TelescopeArray, doi:10.1142/10986.Camilo F., Ransom S. M., Gaensler B. M., Lorimer D. R., 2009,ApJ, 700, L34Chevalier R. A., 1977, in Schramm D. N., ed., Astrophysicsand Space Science Library Vol. 66, Supernovae. p. 53,doi:10.1007/978-94-010-1229-4 5Chevalier R. A., 2005, ApJ, 619, 839Dubner G., Giacani E., Decourchelle A., 2008, A&A, 487, 1033Fang J., Zhang L., 2010, A&A, 515, A20Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,PASP, 125, 306Gaensler B. M., Slane P. O., 2006, Annual Review of Astronomyand Astrophysics, 44, 17Gaensler B. M., Pivovaroff M. J., Garmire G. P., 2001, ApJ, 556,L107Gelfand J. D., 2017, Radiative Models of Pulsar Wind Nebu-lae. Springer International Publishing, Cham, pp 161–186,doi:10.1007/978-3-319-63031-1 8, https://doi.org/10.1007/978-3-319-63031-1_8
Gelfand J. D., Slane P. O., Zhang W., 2009, ApJ, 703, 2051Ghisellini G., Guilbert P. W., Svensson R., 1988, ApJ, 334, L5Ginzburg V. L., Syrovatskii S. I., 1964, The Origin of CosmicRays. MacmillanGranot J., Gill R., Younes G., Gelfand J., Harding A., Kouve-liotou C., Baring M. G., 2017, MNRAS, 464, 4895H. E. S. S. Collaboration et al., 2017, preprint, 612( arXiv:1706.04535 )H. E. S. S. Collaboration et al., 2018a, A&A, 612, A1H. E. S. S. Collaboration et al., 2018b, A&A, 612, A2Hahn J., 2015, in 34th International Cosmic Ray Conference(ICRC2015). p. 917Helfand D. J., Becker R. H., 1987, ApJ, 314, 203Hunter J. D., 2007, Computing in Science and Engineering, 9, 90Kennel C. F., Coroniti F. V., 1984, ApJ, 283, 694Kistler M. D., 2015, preprint, ( arXiv:1511.01159 )Kn¨odlseder J., et al., 2016, A&A, 593, A1Lande J., et al., 2012, ApJ, 756, 5Longair M. S., 2008, Galaxy FormationLu F.-W., Gao Q.-G., Zhu B.-T., Zhang L., 2019, A&A, 624, A144Mart´ın J., Torres D. F., Rea N., 2012, MNRAS, 427, 415Mart´ın J., Torres D. F., Pedaletti G., 2016, MNRAS, 459, 3868Ostriker J. P., Gunn J. E., 1971, ApJ, 164, L95Pacini F., Salvati M., 1973, ApJ, 186, 249Popescu C. C., Yang R., Tuffs R. J., Natale G., Rushton M.,Aharonian F., 2017, MNRAS, 470, 2539Porquet D., Decourchelle A., Warwick R. S., 2003, A&A, 401, 197Protassov R., van Dyk D. A., Connors A., Kashyap V. L., Siemigi-nowska A., 2002, ApJ, 571, 545Qiao W.-F., Zhang L., Fang J., 2009, Research in Astronomy andAstrophysics, 9, 449Rieger F. M., de O˜na-Wilhelmi E., Aharonian F. A., 2013, Fron-tiers of Physics, 8, 714Sidoli L., Mereghetti S., Israel G. L., Bocchino F., 2000, A&A,361, 719Smith A. W., the VERITAS Collaboration 2015, preprint,( arXiv:1508.06311 )Tanaka S. J., Takahara F., 2010, ApJ, 715, 1248Tanaka S. J., Takahara F., 2011, ApJ, 741, 40 Tibolla O., Komin N., Kosack K., Naumann-Godo M., 2008,in Aharonian F. A., Hofmann W., Rieger F., eds, Ameri-can Institute of Physics Conference Series Vol. 1085, Amer-ican Institute of Physics Conference Series. pp 249–252,doi:10.1063/1.3076652Torres D. F., Cillis A. N., Mart´ın Rodriguez J., 2013, ApJ, 763,L4Torres D. F., Cillis A., Mart´ın J., de O˜na Wilhelmi E., 2014,Journal of High Energy Astrophysics, 1, 31Truelove J. K., McKee C. F., 1999, The Astrophysical JournalSupplement Series, 120, 299Venter C., de Jager O. C., 2007, in Becker W., Huang H. H., eds,WE-Heraeus Seminar on Neutron Stars and Pulsars 40 yearsafter the Discovery. p. 40 ( arXiv:astro-ph/0612652 )Vorster M. J., Tibolla O., Ferreira S. E. S., Kaufmann S., 2013,ApJ, 773, 139Wenger M., et al., 2000, A&AS, 143, 9Zhang L., Chen S. B., Fang J., 2008, ApJ, 676, 1210Zhu B.-T., Fang J., Zhang L., 2015, MNRAS, 451, 3145Zhu B.-T., Zhang L., Fang J., 2018, A&A, 609, A110de Jager O. C., Djannati-Ata¨ı A., 2009, in Becker W., ed.,Astrophysics and Space Science Library Vol. 357, Astro-physics and Space Science Library. p. 451 ( arXiv:0803.0116 ),doi:10.1007/978-3-540-76965-1 17de Jager O. C., Harding A. K., Michelson P. F., Nel H. I., NolanP. L., Sreekumar P., Thompson D. J., 1996, ApJ, 457, 253de O˜na-Wilhelmi E., et al., 2013, Astroparticle Physics, 43, 287van Rensburg C., Kr¨uger P. P., Venter C., 2018, MNRAS, 477,3853van Rensburg C., Venter C., Seyffert A. S., Harding A. K., 2020,MNRAS, 492, 3091van der Walt S., Colbert S. C., Varoquaux G., 2011, Computingin Science and Engineering, 13, 22
APPENDIX A: GENERAL FORM OF ENERGYDIFFUSION EQUATION
Here we describe in detail the energy-diffusion equation usedin this work, starting from its general, non-simplified form(Ginzburg & Syrovatskii 1964): ∂ N i ( E , (cid:174) r , t ) ∂ t = ∇ · [ D i ( E , (cid:174) r , t )∇ N i ( E , (cid:174) r , t )] − ∂∂ E [ b i ( E ) N i ( E , (cid:174) r , t )] + ∂ ∂ E [ d i ( E ) N i ( E , (cid:174) r , t )] + Q i ( E , (cid:174) r , t )− N i ( E , (cid:174) r , t ) τ i ( E , (cid:174) r , t ) + (cid:213) k ∫ P ki ( E (cid:48) , E ) N k ( E (cid:48) , (cid:174) r , t ) dE , (A1) N i ( E , (cid:174) r , t ) is the number density of particles species denotedwith the subscript i . The first term on the right-hand sidedescribes the spatial diffusion of particles inside the nebulaand D i ( e , (cid:174) r , t ) is the diffusion coefficient. The second term de-scribes the continuous energy variation due to accelerationprocesses and energy losses, including adiabatic, synchrotronand IC losses. The function b i ( E ) is the mean energy varia-tion of the particle in unit time. The third term is related tofluctuations in this continuous variation of energy of the par-ticles, whereas the function d i ( E ) is equal to the mean squareof the energy variation per unit time. The term Q i ( E , (cid:174) r , t ) isthe particle injection rate, which in this case originates fromthe acceleration of the particles at the termination shock.The fifth term accounts for the escape of particles from thesystem with the characteristic timescale τ i ( E , (cid:174) r , t ) . Finally, MNRAS000
Here we describe in detail the energy-diffusion equation usedin this work, starting from its general, non-simplified form(Ginzburg & Syrovatskii 1964): ∂ N i ( E , (cid:174) r , t ) ∂ t = ∇ · [ D i ( E , (cid:174) r , t )∇ N i ( E , (cid:174) r , t )] − ∂∂ E [ b i ( E ) N i ( E , (cid:174) r , t )] + ∂ ∂ E [ d i ( E ) N i ( E , (cid:174) r , t )] + Q i ( E , (cid:174) r , t )− N i ( E , (cid:174) r , t ) τ i ( E , (cid:174) r , t ) + (cid:213) k ∫ P ki ( E (cid:48) , E ) N k ( E (cid:48) , (cid:174) r , t ) dE , (A1) N i ( E , (cid:174) r , t ) is the number density of particles species denotedwith the subscript i . The first term on the right-hand sidedescribes the spatial diffusion of particles inside the nebulaand D i ( e , (cid:174) r , t ) is the diffusion coefficient. The second term de-scribes the continuous energy variation due to accelerationprocesses and energy losses, including adiabatic, synchrotronand IC losses. The function b i ( E ) is the mean energy varia-tion of the particle in unit time. The third term is related tofluctuations in this continuous variation of energy of the par-ticles, whereas the function d i ( E ) is equal to the mean squareof the energy variation per unit time. The term Q i ( E , (cid:174) r , t ) isthe particle injection rate, which in this case originates fromthe acceleration of the particles at the termination shock.The fifth term accounts for the escape of particles from thesystem with the characteristic timescale τ i ( E , (cid:174) r , t ) . Finally, MNRAS000 , 1–18 (2020) onstraining Models of the PWN in SNR G0.9+0.1 with CTA the last term accounts for the creation and annihilation ofparticles with a probability distribution P ki ( E (cid:48) , E ) (Ginzburg& Syrovatskii 1964).The equation A1 cannot be easily solved. Suitable ap-proximations are usually made. First of all, we consider onlyone population of particles (electrons), we neglect pair cre-ation or annihilation and we take only the mean value of theenergy losses per unit energy, neglecting any fluctuations inthe continuous energy variation. We also assume an isotropicdistribution of electrons, an isotropic injection term insidethe nebula and a uniform magnetic field (no diffusion effectinside the nebula). With these approximations we can ne-glect the first, the third and the last term in equation A1,that becomes equation 10 from Section 6.1. The escape termin equation A1 is retained, even if we neglect the other dif-fusive terms. Therefore, particles are allowed to escape fromthe nebula, although we do not treat in detail the diffusionprocess. APPENDIX B: RADIUS AND VELOCITYEVOLUTION OF PWN
In this appendix we describe an iterative method similar tothat from Gelfand et al. (2009), which we use to compute theradius R pwn and the expansion velocity v pwn of the PWN ineach time step. For this, it is necessary to take into accountan interaction between the SNR and the PWN expandinginside it.First of all, we determine the properties of the ejectedmaterial between the reverse shock of the remnant and thenebula. Making a standard assumption that an inner corewith initially constant density is surrounded by an outerenvelope with density proportional to r (Blondin et al. 2001;Truelove & McKee 1999; Gelfand et al. 2009), the density ofthe ejecta can be written as: ρ ej ( r , t ) = (cid:40) π E sn v − t t − for r ≤ v t t π E sn v − t t − ( rv t t ) − for r > v t t , (B1)where v t = ( E sn / M ej ) / is the transition velocity be-tween the constant density core and the outer envelope, E sn is an energy of the supernova explosion and M ej is its ejectedmass. The ejecta during this stage is expanding ballisticallyand, therefore, its velocity is equal to v ej = r / t . Since inthis work we study young PWNe, which have not reach thereverse shock of the SNR yet, we are not aiming in furthermodelling of the ejecta.We adopt a thin-shell approximation (Chevalier 2005),considering that the expanding PWN is surrounded by athin shell of swept-up material.Initial condition for our iterative procedure, which esti-mates the radius and the associated expansion velocity, aredetermined as described below. Considering the standardapproximation of an isobaric bubble inside the thin-shell,where the adiabatic losses are dominant, the equation of mo-tion of the mass of the shell M s can be written as (Ostriker& Gunn 1971; Chevalier 1977) M s d Rdt = π R pwn (cid:34) P pwn − P ej − ρ ej (cid:18) dR pwn dt − v ej (cid:19) (cid:35) , (B2)where ρ ej , v ej and P ej are computed at R pwn , and P pwn is the pressure inside the nebula. Since in this phase P pwn (cid:29) P ej , we can simplify neglecting the second term in the righthand side of the equation. From the first law of thermody-namics we can write the following expression: dE pwn dt = L ( t ) − π P pwn R pwn dR pwn dt . (B3)This equation is possible to solve in the approximation of t (cid:28) τ where L ( t ) (cid:39) L . Putting together equations B1,B2 and B3, we obtain the following initial condition for theradius and expansion velocity (Chevalier 1977; Blondin et al.2001): R pwn ( t ) = . (cid:169)(cid:173)(cid:171) E sn L M ej (cid:170)(cid:174)(cid:172) / t / , (B4) v pwn ( t ) ≡ dR pwn dt ( t ) = R pwn ( t ) t . (B5)With this initial condition we can start the iterations,computing new radius of the PWN ( R pwn ( t + ∆ t ) ) togetherwith the magnetic field in the nebula B pwn ( t + ∆ t ) (equation20), the spin-down luminosity L ( t + ∆ t ) (equation 12), themaximum energy of the electrons E max ( t + ∆ t ) (equation 18),and the density and the velocity of the ejecta at R pwn ( t + ∆ t ) . R pwn ( t + ∆ t ) = R pwn ( t ) + v pwn ( t ) ∆ t . (B6)As a second step, we computed the pressure inside thenebula, in order to determine the force acting on the shelland, therefore, a new value of the expansion velocity of thePWN. The net force which affects the shell is proportional tothe difference between the pressure inside P pwn and outsidethe nebula P ej : F pwn ≡ ddt ( M s v pwn ) = π R pwn ( P pwn − P ej ) . (B7)However, the second term of this expression can be neglectedsince it is expected that P pwn (cid:29) P ej .The total pressure inside the nebula is determined as asum of the pressure of the magnetic field P pwn , B and thatof the moving electrons P pwn , e . Calculating the value of themagnetic field B pwn from equation 20, we can determine theenergy stored in the magnetic field: E pwn , B ( t ) = (cid:32) B pwn ( t ) π (cid:33) π R pwn ( t ) . (B8)From equation B8 we obtain P pwn , B as: P pwn , B ( t ) = E pwn , B ( t ) π R pwn ( t ) = B pwn ( t ) π . (B9)The contribution of the second component P pwn , e canbe computed solving equation 10 and extracting the totalenergy from the spectrum of evolved particles: E pwn , e ( t ) = ∫ E max E min E N ( E , t ) dE . (B10)Then, the electron pressure is found as follows: P pwn , e ( t ) = ( γ pwn − ) E pwn , e ( t ) π R pwn ( t ) = E pwn , e ( t ) π R pwn ( t ) , (B11) MNRAS , 1–18 (2020) Fiori, Zampieri, Burtovoi, et al. where γ pwn is equal to 4/3.Finally, we are able to compute new expansion velocityof the nebula. If v pwn ( t ) > v ej ( t ) the new mass of the shellbecomes M s ( t + ∆ t ) = M s ( t ) + π (cid:104) R pwn ( t + ∆ t ) − R pwn ( t ) (cid:105) ρ ej ( t + ∆ t ) . (B12)Otherwise, new mass M s ( t + ∆ t ) is simply equal to M s ( t ) . Thenew velocity v pwn ( t + ∆ t ) , which will be used for calculatingthe radius of PWN in the next iteration, can be found fromthe following expression: v pwn ( t + ∆ t ) = M s ( t ) v pwn ( t ) + ∆ M s v ej ( t ) + F pwn ( t ) ∆ tM s ( t + ∆ t ) , (B13)where ∆ M s = M s ( t + ∆ t ) − M s ( t ) .To compute an evolution of leptons using this iterativeprocedure, we solve advective equation 10 many times. Incase of high energy losses these computations can becometime consuming. To speed up the calculations, we put an up-per limit on the magnetic field inside the nebula during thefirst stages of evolution of the system. We impose that mag-netic field does not exceed µ G during the first 5 yrs andit is < µ G up to 500 yrs of evolution. These constrainsintroduce modest impact to the calculation of the radius ofthe source. Resulting value of the radius is < higher thanthat computed with no upper limits on the magnetic field. Itis worth to mention that this approximation has been testedonly for SNR G0.9+0.1 and may not be valid for youngersources (less than ∼ years), where an higher thresholdfor the magnetic field will be probably needed to better re-produce the observed data. We finally note that once thevalues needed to determine an evolution of the nebula areobtained, we recalculate the particle spectrum without anylimit on the magnetic field. We also checked that the finalphoton spectrum does not differ significantly from that ob-tained using no upper limits on the magnetic field. APPENDIX C: MODEL FITTING
In our fitting procedure we first compute a grid of mod-els spanning a large range of values of free parameters. Wethen compute the chi-square χ statistics for each model ofthe grid and the observational data, and choose the best-fitmodel with the minimal χ . As mentioned in Section 6.2,we leave free to vary only 3 parameters: α , η B and (cid:15) . Othertwo parameters E b and α are fixed to values as in Zhu et al.(2018) in order to perform comparison with their results. Fi-nally, we estimate uncertainties of free parameters using thefollowing procedure: • We produce a three-dimensional (3D) probability gridfrom the χ values obtained for all the models: P D ( α , η B , (cid:15) ) ∝ exp (cid:16) − χ / (cid:17) , (C1) • and normalize it: (cid:213) α ,η B ,(cid:15) P D ( α , η B , (cid:15) ) = , (C2) • We then extract the marginalized (1D) probability dis- tribution for each parameter summing over other two pa-rameters: P D ( α ) = (cid:213) η B ,(cid:15) P D ( α , η B , (cid:15) ) , (C3) P D ( η B ) = (cid:213) α ,(cid:15) P D ( α , η B , (cid:15) ) , (C4) P D ( (cid:15) ) = (cid:213) α ,η B P D ( α , η B , (cid:15) ) . (C5) • Finally, using these marginalized probability distribu-tions, we estimate the confidence interval and σ error foreach parameter, assuming that the distributions are Gaus-sians. MNRAS000