The Effect of Pulsar Geometry on the Observed Gamma-ray Spectrum of Millisecond Pulsars

We analyse 13 yrs of $\textit{Fermi}$-LAT PASS 8 events from 127 gamma-ray emitting millisecond pulsars (MSPs) in the energy range 0.1$-$100 GeV and significantly detect 118 MSPs. We fit the stacked emission with a log parabola (LP) spectral model which we show is preferred to two previously published models. We consider the influence of pulsar properties and observer geometric effects on spectral features by defining energy flux colours for both the individual MSPs, and our stacked model as a baseline. There is no correlation of colours with pulsar luminosity, $\dot{E}$, surface magnetic field or magnetic impact angle. We also find that pulsar geometry has little effect on the observed gamma-ray spectrum which is in tension with previous modelling of gamma-ray emission with respect to pulsar geometry. Our LP MSP model is applicable to problems where an ensemble of gamma-ray MSPs is considered, such as that of the Galactic centre excess or in the case of emission from globular clusters.


INTRODUCTION
In 2009 the Fermi-LAT established 8 millisecond pulsars (MSPs) as pulsed gamma-ray emitters using just 8 months of observations combined with radio timing data Abdo et al. (2009b).The initial discoveries of Abdo et al. (2009b) showed that the pulse peaks of radio and gamma-ray emission in MSPs need not be co-incident, favouring different emission regions for radio and gamma-rays and supporting the prevailing slot-gap and outer gap (OG) models of gamma-ray emission (Muslimov & Harding (2004) and Zhang & Cheng (2003) respectively).
The most recent fully published catalogue of gamma-ray pulsars is "The Second Fermi Large Area Telescope Catalog of Gamma-Ray Pulsars" or 2PC, compiled by the Fermi-LAT collaboration Abdo et al. (2013) 1 .The 2PC uses Fermi-LAT observations to identify 40 MSPs above 0.1 GeV.
At higher energies of 10 GeV−2 TeV, the 3FHL catalog of Ajello et al. (2017) identifies 15 MSPs at Galactic latitude || ≥ 10°with 13 previously identified as pulsed gamma-ray emitters at lower energies using the LAT and the remaining 2 radio pulsars previously undetected in gamma-rays.
In a more recent study than the 2PC, an analysis of 19 MSPs with 7 years of pass 8 Fermi-LAT event data in the range 0.1−300 GeV shows that the empirical MSP death-line (or  below which gamma-ray emission cannot be detected) is 8×10 32 erg s -1 Guillemot et al. (2016).This analysis also demonstrates that   is uncorrelated with , showing a nearly 2 orders of magnitude of scatter for MSPs with an  of 10 34 erg s -1 .For  > 5×10 34 erg s -1 ,   is directly proportional to √ .Similarly the wide variation of the relationship between ,   and  makes it difficult to deduce gamma-ray luminosities a priori from MSP timing information alone.
In contrast, the geometry of individual MSPs is increasingly well understood, through the fitting of models to gamma-ray and radio light curves, with the angles for both the line of sight and magnetic axis, with respect to the pulsar spin axis determined, for all 40 MSPs in the 2PC Johnson et al. (2014) and more recently for 10 MSPs with high signal-to-noise observations Benli et al. (2021).
We note that there is recent modelling and interest in the effects of pulsar inclination (among other variables) on the observed gamma-ray spectra of MSPs ranging from a general one Vigano et al. (2015a,b), to a spectral softening below 400 MeV, to being a lesser influence on the high-energy spectrum Torres et al. (2019); Giraud & Petri (2021).Thus, the effect of pulsar magnetosphere geometry and magnetic inclination angle on the observed gammaray spectra is an open question.Therefore, inclination effects could well affect the the observed gamma-ray spectrum, although to our knowledge this has not been confirmed by observation.
A demonstrable link between pulsar inclination and observed spectra could serve as a constraint on spectral emission models where the pulsar geometry is known, whilst conversely the geometry of the pulsar system may be constrained from the observed gammaray spectrum.We are therefore motivated in this paper to determine for the first time if there is a link between pulsar geometry and observed gamma-ray spectral features.
In this work, we seek to construct the best general spectral model for 118 MSPs and use it to examine the effect of geometry on gamma-ray spectra, whilst excluding any systematic effect of other pulsar properties (surface magnetic field, luminosity , ).We analyse a larger sample of 127 MSPs, rather than the 39 MSPs of the 2PC Abdo et al. (2013), as used by the previous stacked MSP analyses of Xing &Wang (2016) andMcCann (2015).Furthermore, we note that the MSP spectral model we produce is very useful in its own right, and could be applied to any problem where gamma-ray emission is presumed to originate from an ensemble of MSPs plus another component such as emission from dark matter, as relevant to globular clusters or the Galactic centre excess .
This paper is structured as follows.In Section 2 we describe the selection criteria of MSPs for analysis.In Section 3 we describe our analysis method for the detection of MSPs, and in Section 4 we present our results.In Section 5 we derive a spectral model for the stacked emission of detected MSPs.In Appendix A we compare our model to those of Xing &Wang (2016) andMcCann (2015) and determine whether there is any relationship between colours constructed from low and medium energy MSP fluxes, and pulsar , luminosity and surface magnetic field.We also consider the spectral predictions of synchro-curvature emission models incorporating pulsar magnetic inclination and determine if they are suitable for comparison with the observed gamma-ray spectra.Finally in Section 7 we summarise our findings and make suggestions for future work.

MILLISECOND PULSAR SELECTION
We select all 127 MSPs (defined as any pulsar with period ≤ 30 ms) from the Public List of LAT Detected Gamma-Ray Pulsars2 .Their co-ordinates, periods and  values are shown in Table D1.

Initial Photon Event Data Selection
The data in this analysis were collected by Fermi-LAT between 4th Aug 2008 (15:43 h) to 26th August 2021 (00:52 h), (Mission Elapsed Time (MET) 2395574147[s] to 651631981 [s]).We select all pass 8 events which are source class photons (evclass = 128), and PSF3 events with the best quartile direction reconstruction, (evtype = 32), spanning the energy range 0.1−100 GeV.Throughout our analysis, the Fermipy software package3 with version 2.0.8 of the Fermitools is used, in conjunction with the p8r3_source_v3 instrument response functions.We apply the standard pass8 cuts to the data, including a zenith angle 90°cut to exclude photons from the Earth limb, and good-time-interval cuts of DATA_QUAL > 0 and LAT_CONFIG = 1.The energy binning used is 4 bins per decade in energy and spatial binning is 0.1°per image pixel.

Spectral Models
The differential flux, dN/dE, (photon flux per energy bin) of an individual MSP, is described by a spectral model which is either an exponential cut-off power law 2 (Eqn. 1, PLSuperExpCutoff2 ), a log parabola (Eqn.2, LP), or a power law ( Eqn. 3, PL) 4 .
where normalisation (also known as prefactor) =  0 , index1 =  1 ,  0 is the scale,  is the exponential factor and index2 =  2 (which controls the sharpness of the exponential cut-off).
where norm =  0 , and  b is a scale parameter.

MSP Likelihood Analysis
We perform a full likelihood analysis on all 127 MSP targets (Table D1), in the energy range 0.1−100 GeV, using a 25°Radius of Interest (ROI) centred on the nominal MSP co-ordinates and a 40°S ource Radius of Interest (SROI) for each MSP target.The model we use in our likelihood analysis consists of a point source population seeded from the Fermi-LAT's 10 yr fourth point source catalog, data release 2, (4FGL-DR2, gll_psc_v27.fit), extended gamma-ray sources and diffuse gamma-ray emission.The diffuse emission detected by the Fermi-LAT consists of two components: the Galactic diffuse flux, and the isotropic diffuse flux.The Galactic component is modelled with Fermi-LAT's gll_iem_v07.fitspatial map with the normalisation free to vary.The isotropic diffuse emission is defined by Fermi's iso_P8R3_SOURCE_V3_PSF3_v1.txt tabulated spectral data.The normalisation of the isotropic emission is left free to vary.
In addition, all known sources (including MSPs) take initial spectral parameters and position from the 4FGL.
Initially, the setup and optimize methods are run to create count and photon exposure maps and to compute the test statistic (TS) values of all 4FGL sources in the model.The TS value is defined as Eqn. 4 where  0 is the maximum likelihood for a model without a source observations (i.e. the null hypothesis) and  1 is the maximum likelihood for an additional source observation at a given location.
The fit method is then run: fit is a likelihood optimisation method which executes a fit of all parameters that are currently free in the the model and updates the TS and predicted count ( ) values of all sources.From this fit, all point sources with a TS < 4, or with a predicted number of photons < 4 are removed from the model.
The normalisation of all sources within 10°of the MSP are then freed using the free_source method to allow for the Point Spread Function (PSF) of PSF3 events down to 100 MeV (95 % containment at 8°) .The source nearest to the catalogue position of the MSP has prefactor and index spectral parameters (Eqn.3) freed for power law sources and prefactor, index1 and expfactor freed for PLSuperExpCutoff2 sources (Eqn.1).
The shape and normalisation parameters of all sources with a TS > 25 are then individually fitted using the optimize method.Then the fit method is run twice more with an intervening find_sources step.Find_sources is a peak detection algorithm which analyses the test statistic (TS) map to find new sources over and above those defined in the 4FGL model by placing a test point source, defined as a power law (PL) with spectral index 2.0, at each pixel on the TS map and recomputing likelihood.Finally, the sed method generates a spectral energy distribution, with energy dispersion disabled for all MSPs which are known 4FGL sources, and a 2 confidence limit on the determination of instrument upper limits.
We also list the best fit spectral models for MSPs with a PLSu-perExpCutoff2, LP or PL spectral model, in Tables C1, C2 and C3, respectively.For these spectral models we provide a selection of representative spectral energy distributions in Figs.C1 and C2.
We determine the average differential energy flux per energy bin (for flux points of ≥ 2 significance) and the standard error of the mean for the detected MSPs (Table 1).The energy bin centres and lower and upper bin energy (extent) arise from a binning of 4 bins per decade of energy as specified in the analysis.Note that the error (energy dispersion) of the bin extent is < 10 % between 1−100 GeV increasing to 20 % at 0.1 GeV.
The majority of MSPs display emission across the whole energy range between 0.2−10 GeV, while emission between 10−18 GeV is seen in only 41 % of MSPs, and 15 % of MSPs above 18 GeV.

GEOMETRIC EFFECTS ON THE OBSERVED GAMMA-RAY SPECTRA OF MILLISECOND PULSARS
The geometry of a pulsar system is described by (1) the angle between the observer line-of-sight (LOS) and the pulsar spin axis(), (2) the angle between magnetic axis and the pulsar spin axis, also known as the magnetic obliquity () and (3) the orbital inclination angle between the pulsar orbital plane and the observer LOS ().
An impact parameter () describing the closest approach of the magnetic axis and the LOS is defined as  =  −  as in Johnson et al. (2014).
In Guillemot & Tauris (2014),  is derived for Galactic field MSPs, with helium white dwarfs as binary companions, from the relationship of pulsar and white dwarf mass, and the orbital period.Over the life of the system, the values of  and  with respect to the pulsar rotation axis are shown to align, with the detection of gamma-ray emission in Galactic MSPs appearing to be somewhat favoured by smaller values of .
In Vigano et al. (2015a) the effect of parameters and assumptions in the OG model are examined.The magnetic inclination angle of a pulsar is shown to affect the magnetosphere geometry, position of the OG and radius of curvature (r c ). r c is identified as the most important parameter having an order of magnitude impact on the parallel electric field, (responsible for lepton acceleration and consequent curvature radiation), which is noted to affect the observed gamma-ray spectrum qualitatively.
In their follow-up paper, Vigano et al. (2015b) model exemplar pulsar gamma-ray spectra obtained by integrating the energy loss of a single particle with distance travelled, convolved with an effective observed particle distribution which incorporates geometry and beaming effects.They show that photons from the initial part of the particle trajectory exhibit softer low energy spectra below 200 MeV, with an index of 0.68 vs an index of 1.14 for emission from the whole trajectory.They note that such spectra could be observed for pulsars with a favourable viewing geometry of the outer gap.
The effect of inclination on the observed gamma-ray spectrum is also considered by Cerutti et al. (2016) who model the combined low-energy synchro-curvature radiation from the polar cap and higher-energy pulsar synchrotron radiation (HES) from the equatorial current sheet (ECS) using a particle-in-cell (PIC) simulation.They produce phase-averaged observed pulsar energy spectra between 400 keV − 4 GeV for varying  and .These spectra show a trend of decreasing energy flux and a softening of the spectrum between 40−400 MeV, with the disappearance of the HES spectral peak at ∼ 200 MeV as  decreases from 90 to 0°for  = 90°.In contrast, varying  from 0 to 90°for  angles of 45 and 60°generally results in few spectral shape changes.
In Petri (2019), a minimalist model determines lepton particle trajectories and velocities depending only on local  field using the vacuum rotator dipole model (VRDM) of pulsars as first described by Deutsch (1955).The gamma-ray spectrum arising from the resulting curvature radiation is determined and used to produce sky-maps and light curves for varying pulsar phase, angle  and magnetic inclination angle ( in that paper rather than ), with respect to the pulsar rotation axis.This modelling is able to reproduce spectral cut-offs at a few GeV and single or double peaked light curves broadly consistent with Fermi-LAT observations of pulsars.It is also shown that there is a slight hardening of the spectrum be-  1.The count and percentage of significant (> 2) flux points for detected MSPs along with mean flux and associated standard error of the mean.
tween 1−10 GeV, when the annulus of the emitting region increases in extent, as it moves closer to the surface of the star, from 0.5−1 of the light cylinder radius (  ) to 0.1−1   .Beyond the light cylinder there is no change in emission spectrum with increasing extent (1−5   ).Overall the spectral shape is noted to be insensitive to geometry.
In Torres et al. (2019), a synchrocurvature spectral emission model is able to reproduce well the observed spectral energy distribution of 18 out of 32 non-MSPs 5 considered, from X-ray to gamma-ray energies.The free parameters of the model are properties intrinsic to the pulsar, namely, the parallel electric field and the magnetic field gradient, and properties related to the viewing geometry of the observer: a length scale for particle emission relative to the radius of the light cylinder and a normalisation of the accelerated particle distribution with respect to distance ( 0 /  ,  0 respectively).This model can reproduce the gamma-ray spectrum of the 4 MSPs considered but not the X-ray spectrum, whilst 3 other MSPs can be well fitted simultaneously in X-rays and gamma-rays.Extending the model so that there are two emission regions (each with own values of ( 0 /  ,  0 , but same intrinsic properties) and visible to a varying extent depending on observer line-of-sight, is shown to fit both the X-ray and gamma-ray observations.This suggests that, whilst broad-band emission from X-rays to gamma-rays can be sensitive to geometric effects and the visibility of different emission regions, gamma-ray emission is less so.The authors also note that the precise location and extent of the accelerating region has previously being shown not to dominate the predicted high-energy spectral shape.
Finally, Giraud & Petri (2021) model the gamma-ray emission arising from synchrotron and curvature radiation in high-altitude slot gaps at the separatrix between open and closed field lines in the pulsar magnetosphere for pulsars of period 2 ms.They use a VRDM with a mono-energetic lepton distribution to produce phaseintegrated and LOS averaged photon count spectra between 1−30 GeV, with spectral shape and energy limits that vary with magnetic obliquity, , and LOS, , whilst retaining a spectral peak at around 4 GeV.Specifically, increasing the magnetic obliquity (for phase and line of sight averaged spectra) from 30 to 90 °increases spectral energy limits from 7 to 30 GeV.However, they conclude that, while 5 Although their selection criteria for a non-MSP is a pulsar with period > 10 ms, their 18 pulsars all have period > 30 ms and thus also meet our definition of a non-MSP.
the pulsar radio spectra significantly depend on the magnetic obliquity, the high-energy part of the spectrum is much less sensitive to geometry.
The models of Cerutti et al. (2016); Petri (2019); Torres et al. ( 2019); Giraud & Petri (2021) are broadly similar in that they (1) derive the curvature radiation spectrum using a Bessel function of order 5/3, (2) are force-free allowing leptons to travel along magnetic field lines, because the screening of the pulsar electric field by plasma is excluded and (3) balance lepton acceleration and braking due to radiative friction.Petri (2019) and Giraud & Petri (2021) employ a mono-energetic particle distribution whereas Torres et al. (2019) spatially vary the number density of particles.
For the models above, we summarise the effect of geometry on the gamma-ray spectrum in Table 2, varying from no effect to changes in energy limits and spectral shape with varying .However, the models which present explicit spectra for specific geometries are difficult to compare directly with gamma-ray observations of MSPs, as they either do not provide a full co-varying range of  and  (Giraud & Petri (2021)), or are designed to test a model concept rather than to be fitted directly to observations (Cerutti et al. (2016)), and moreover they do not reproduce gamma-ray spectral peaks of a few GeV as seen in MSPs.
As no model spectra are suitable for comparison with observations, we instead consider the observed MSP gamma-ray spectrum to investigate if geometry has any systematic effect.
Specifically, we would like to examine the most rapidly varying part of the MSP spectrum, which is generally away from the cut-off values of the the exponential power law in the 2PC (1.2−5.3GeV), where the spectrum peaks and is often almost flat.Moreover, we wish to probe energy ranges which are noted as being sensitive to emission extent and geometric effects as in Table 2 above.
We therefore define two colour band ratios to examine the requisite energy ranges, a low-energy ratio, (Eqn.5, LE) and a medium energy ratio, (Eqn.6, ME).

𝐿𝐸 =
133 MeV   237 MeV (5) The precise energy centres and limits for LE and ME arise from the binning employed, but the energy bins are also chosen to examine the varying and potentially inclination-sensitive part of the Giraud & Petri (2021) Increasing magnetic inclination increases upper energy limit of  averaged photon spectra from 7 to 30 GeV but high-energy spectra noted to be much less sensitive to geometry than radio bands.If a decreasing value of  results in an increasingly soft spectrum below 400 MeV as in Cerutti et al. (2016) then we should expect a corresponding increase in the LE colour as  decreases from 90 to -45 °.Conversely, if geometric effects systematically affect the spectrum between 1−10 GeV then this should appear as a changing ME colour as  changes, with the caveat that increasing emission extent from within the light cylinder could result in a decreasing ME colour.

THE STACKED MSP SPECTRUM
Whilst we consider the effect of geometry on the spectra of individual MSPs, these spectra can vary between MSPs, and so we must derive a stacked MSP model to serve as a baseline to show if there are systematic spectral changes depending on geometry across the MSP sample.In order to derive this model for an ensemble of Galactic field MSPs, we sum the energy fluxes ( 2    ) in each bin (ignoring ULs) between 0.1−56.2GeV for MSPs detected in the energy range 0.1−100 GeV, with the upper energy range of 56.2 GeV arising solely from the binning employed.We exclude the 1 flux point above 56.2GeV (for the single MSP detection J1903-7051) from the stacked model, as this flux is only marginally significant (TS 5.5).
We define two models for comparison; a super exponential  cut-off power law model spectral model as Eqn. 1, "PLSEC" and a log parabola model as Eqn. 2, "LP".
In the PLSEC model, the parameters index1, Normalisation and exponential factor are left free.The scale and index2 parameters are frozen to 1×10 -3 TeV and 1 respectively.In the LP model, the Norm,  and  parameters are left free, while  b is frozen at 1 TeV.
We then use gammapy version 0.18.26 software to fit the flux summation of 118 MSPs with the LP and PLSEC models in the range 100 MeV−56.2GeV.
In Appendix A we show that the LP model is the best model to describe the stacked spectrum of the 118 MSPs and preferred to the previously published models of Xing &Wang (2016) andMcCann (2015).The parameters of our best fit LP model are presented in Table 3 and all models are presented alongside the stacked MSP spectrum in In Fig. 1.
The LE and ME colours for this LP model are 0.50±0.02and 1.75±0.05respectively.These serve as a baseline model to determine spectral trends for individual MSPs in our selection.We show the energy range used to form the LE and ME colour ratios on Fig. 1.In considering the effect of inclination on spectra, we firstly need to exclude any correlation of LE and ME colours with intrinsic pulsar properties such as , luminosity or characteristic age (/2 ).A plot of MSP spectral colours against luminosity, , age and surface magnetic field indicates that there is no detectable correlation (Fig. 2, Pearson correlation coefficient for ME colours of -0.23, -0.13, -0.30 and -0.16 respectively) We next consider whether the inclination (geometric) properties of the MSP, such as the LOS, , and the magnetic inclination angle, , with respect to the pulsar rotation axis might affect spectral features.Johnson et al. (2014) use a VRDM to derive the best fit  and  values, for the MSPs of the 2PC.They divide the MSPs into 3 classes, Class I, II and III, defined as the MSP's gamma-ray pulse trailing, aligned or leading the MSP's radio pulse, respectively.They then derive inclination angles from a two-pole caustic (TPC) and outer-gap (OG) model for Class I MSPs, a TPC, OG and slot-gap model for Class II MSPs and a pair-starved polar cap model for Class III.We use these magnetic inclination angles to determine the magnetic impact angle, , defined as  =  −  for the three classes of MSPs (Tables 4, 5 and 6).We then plot  for individual MSP colours in Fig. 3, 4, 5.We observe that there is no dependence of ME and LE colours on , implying that LE and ME gamma-ray spectral features of MSPs are not correlated with magnetic inclination and line-of-sight effects.
More recently, Benli et al. (2021) determine  and  for 10 of the same Class I MSPs as Johnson et al. (2014), using a force-free magnetospheric model for gamma-ray emission at the light cylinder and into the striped wind, and radio emission from the polar-cap to reproduce gamma-ray and radio light curves from the 2PC.They then fit the time-lag between radio and first gamma-ray pulse along with the gamma-ray pulse separation and ratio of the gamma-ray pulse heights of light curves in the 2PC for varying inclination and magnetic angle to obtain a best fit.We again determine  from their inclination and magnetic angles (Table 7).
We then plot ,  and  for MSPs in our selection in Fig. 6.The values of  and  are uncorrelated with the model colours we have calculated.However, for , there is some indication of harder ME emission for absolute  values of 5−10°(ME colour ratio ∼1) and generally softer emission outside this range, (ME colour ratio > 2), suggesting that ME spectral features are influenced by the offset of the emission region from the pulsar spin axis, once geometric viewing effects are allowed for.However, the inclination angles of more MSPs would need to be determined with the same force-free model to draw a definitive conclusion, especially as no such effect is seen when using the inclination angles of Johnson et al. (2014).
In conclusion, it appears that pulsar line-of-sight and the magnetic axis inclination seem to have little effect on the observed gamma-ray spectrum at low and medium energies.Therefore the observed gamma-ray spectrum does not allow us to constrain pulsar geometry or to further validate the existing synchro-curvature models of pulsar emission which incorporate geometric effects.However,  and  are determined only for a relatively small pulsar population and so we cannot completely exclude spectral differences arising from pulsar geometry.
Our standalone LP model is also useful in its own right to model ensemble emission from a population of MSPs.Such an application of a stacked MSP model is demonstrated in Brown et al. (2018), wherein, the gamma-ray emission of the globular cluster 47 Tuc is modeled using the MSP model of Xing & Wang (2016) combined with an annihilating dark matter (DM) model to show that a "MSP+DM" model is preferred over MSP emission alone.In this regard we expect our LP model comprising a larger sample   57 % of the photon exposure of our study.The caveat applies that any stacked spectrum should be applied to systems only where the total spectral variation of individual MSPs is less than the model hypothesis under consideration, as demonstrated by the comment on Brown et al. (2018) in Bartels & Edwards (2019) and its reply in Brown et al. (2019).Specifically, our model will be applicable to considering the ensemble gamma-ray emission from MSPs in globular clusters (GCs, Abdo et al. (2009a)), as the characteristic ages 7 of 76 GC MSPs (2.6 × 10 7 − 6.8 × 10 10 yrs) are consistent with those of the model MSPs in Fig. 2. Furthermore, the spectra of disc and bulge MSPs are consistent within uncertainties (Ploeg et al. (2020) ) and so we expect our model based on local disc MSPs to be useful in consideration of the Galactic centre excess problem which could arise from an unresolved population of MSPs (Abazajian ( 2011)).

CONCLUSION
We analyse 127 MSPs from the Public List of LAT Detected Gamma-Ray Pulsars and detect 118 MSPs in the range 100 MeV−100 GeV.We sum the 100 MeV−56.2GeV fluxes of these MSPs and find the best fit to the resulting spectral energy distribution is an LP model and that this model is superior to other published models of stacked MSP emission.
Most previous models of gamma-ray emission from MSPs suggest that pulsar geometry may affect the observed gamma-ray spectrum.These models are difficult to compare with observations as they are either qualitative, do not cover the full range of geometries or are a proof-of-concept of modelling technique which does not reproduce known features of MSP spectra.
In the absence of quantitative spectral emission models for comparison, we use our LP model as a baseline to determine if 7 Calculated for the list of MSPs in GCs (with period and positive period derivative) maintained at http://www.naic.edu/~pfreire/GCpsr.html accessed 29/10/2022 pulsar properties affect spectral shape by using energy flux ratios (or colours) for the low-energy (133/237 MeV) and medium energy (4.2/7.5 GeV) bins.We find that pulsar , surface magnetic field and gamma-ray luminosity are uncorrelated with these colours and hence do not systematically affect the spectral shape in the lowenergy and medium energy bins.Similarly we find that LE and ME colours, are in the main, uncorrelated with magnetic impact angle, indicating that pulsar spin axis inclination and magnetic inclination have little effect on the LE and ME features of the observed gammaray spectrum.There is a hint of symmetry in ME colours for harder emission from magnetic impact angles between 5 and 10 °, but more inclination determinations for a larger pulsar sample would be required to confirm this.
In general, our results show little influence on the low and medium energy gamma-ray spectrum arising from pulsar geometry.
Finally, we note that our MSP spectral model (LP) is useful more generally in the problems of the Galactic centre excess and the ensemble emission in globular clusters arising from a population of MSPs.
The 2PC is a published survey of pulsars observed by Fermi-LAT.It uses 3 years of pass 7 event data in the energy range between 100 MeV and 100 GeV with the 2FGL source catalog as a source model and lists the spectral models and fluxes of 117 pulsars, evenly divided between MSPs, young radio-quiet and young radio-loud pulsars.The survey uses three search strategies for pulsar detection to overcome the difficulty of only one photon being detected in a few million pulsar rotations.Firstly, the known rotation ephemerides of pulsars, obtained mostly through radio and in some cases X-ray observations, are used to fit a timing model with tempo/tempo2 software, to tag the gamma-ray event data with a pulsar phase.The gamma-ray data are then phase-folded to identify any emission peaks.Secondly, blind periodicity searches are used on unassociated sources classed as candidate pulsars because they show no variability and have spectra that can be fitted by an exponential cut-off in the GeV band.This method is challenging because the event data are sparse with only a few photons detected per hour and in addition pulsars may have been missed by being in binary systems, which tends to smear the signal through Doppler shifts arising from orbital motion.Finally, the detection of pulsed radio emission in unassociated sources and the construction of timing models can lead to the detection of gamma-ray pulsations through phase folding methods as above.The 2PC lists 40 MSPs, 20 of which have been detected using this final method.The 2PC increased the then known MSP sample from 8 to 40 MSPs with heliocentric distances up to 2 kpc and a uniform distribution in the sky.The MSPs exhibit between 1−3 gamma-ray peaks and their differential flux spectrum, dN/dE, (photon flux per energy bin) is an exponential cut-off power law as described by Eqn.A1.This equation is functionally equivalent to Eqn. 1 with the exponential factor  replaced by 1/ cut and index  1 having negative sign.For consistency with the 2PC and Xing & Wang (2016) we continue to use symbols , Γ and  which are equivalent to the  0 ,  1 and  2 in Eqn. 1 The MSPs listed in the 2PC also provide a pulsar sample to use in determining models of stacked gamma-ray emission.Xing & Wang (2016) re-analyse 39 of the 40 MSPs in the 2PC (excluding one, J1939+2134, which has a detection significance of ≃ 3 ) using 7.5 years of pass 8 event data in the energy range 100 MeV−300 GeV in 15 energy bins, using the 3FGL as a source catalog model.
In the 3FGL, 33 MSPs are described by an exponential cut-off power law model (Eqn.A1), whereas 6 MSPs are best fitted with a simple power law (PL) (Eqn.3).However, in their analysis they find that an exponential cut-off can be detected in the 6 PL MSPs at > 3  significance.They therefore use an exponential cut-off power law model throughout their analysis.They then stack all flux points from the 39 MSPs with a TS > 9 (equivalent to > 3  significance) to obtain a functional form described by an exponential cut-off power law as Eqn.A1 with Γ = 1.54 +0.10 −0.11 and   = 3.70 +0.95 −0.70 but with  0 and  equal to 1 (hereafter the "Xing and Wang" model).Finally, they recommend that this functional form can be used as a model to find candidate MSPs in unidentified Fermi-LAT sources at high Galactic latitudes.
Alternatively, McCann (2015) constructs a stacked MSP gamma-ray spectrum using an aperture photometry (AP) method rather than likelihood analysis.The AP approach has the advantage over the likelihood approach in that it is model independent and less computationally intensive.However, it does require timing information for the pulsars analysed.McCann chooses 39 MSPs from the 2PC (excluding a different one, J2215+5135 because its off-phase, where emission is at a minimum, is undefined) and consider 4.2 years of pass 7 event data per MSP in the energy range 100 MeV−1 TeV binned at 4 bins per decade of energy.McCann then uses the tempo 2 software to barycentre and phase fold the photon events.McCann then obtains the energy excess counts of all events outside the off-phase (i.e the on-phase), distributed by energy, corrects for exposure and produces a spectral energy distribution from the stacked fluxes.Finally McCann fits the differential flux  2    (as opposed    ) with an exponential cut-off power law, with a functional form as Eqn.A2.This exponential cut-off power law has Γ = 0.54 ± 0.05,  cut = 3.60 ± 0.21 GeV and  = 0.7 ± 0.15 (hereafter the "McCann" model).McCann also makes a check on the performance of the AP method vs the likelihood method of the 2PC by defining a flux ratio for the MSPs of    2   which varies between 0.8 and 0.9 for energies of 250 MeV−8 GeV.
We assess which MSP model(s) (LP, PLSEC, McCann or Xing and Wang) are a preferred description of our stacked MSP spectrum through a likelihood analysis.We allow the normalisation of the models to vary and determine for each normalisation the residual In Fig. 1, we show the likelihood best fit models and the stacked spectrum of 118 MSPs with the LP model fitting the greatest number of flux points.
The preferred MSP model can also be determined using the minimum value of the Akaike Information Criterion (AIC) statistic Akaike (1974), (Eqn.A3).The AIC ranks how well a model fits a data set (compared to other models) and penalises the over-fitting which results from the model having more free parameters.The AIC is a relative measure in that it allows a set of models to be compared with the model exhibiting the lowest AIC score considered superior in that set, but it does not allow a determination of whether any model is best in an absolute sense.The AIC is defined in Eqn.A3 where  is the number of free model parameters and ℓ is the likelihood of the best fit model.
For the purposes of spectral model comparison, a more convenient definition of AIC is Eqn.A4 where  is the number of flux data points or energy bins and  is the residual sum of squares as defined in Eqn.A5 with   being the observed flux and  (  ) the flux predicted by the model for at an energy   for an energy bin .
The AIC statistic for the fit of our models, LP, PLSEC, and the models of Xing and Wang,and McCann, respectively.Our model, LP, having the minimum AIC statistic, is thus the preferred one.Furthermore the evidential significance Δ  of any model, , can be determined as Δ  = (AIC value for model  -AIC value for our model).A Δ  value ≤ 2 indicates a model has substantial support whereas Δ  ≥ 10 indicates a model has essentially no support (Burnham & Anderson (2004)).The Δ  values for Xing and Wang,and McCann are 9.6 and 28.3 respectively and hence these models are disfavoured by the AIC statistic.
The absolute goodness-of-fit of models to the stacked MSPs can also be determined by the  2 statistic which we calculate across the energy range of each bin and take the minimum (best)  2 value in each bin for the best fit likelihood model.The LP, PLSEC, Xing and Wang and McCann models have a  2 statistic of 6.6, 192.0, 1225.5 and 2472.8 respectively between 100 MeV−56.2GeV (Table A2).
The Xing and Wang model is an acceptable fit from 178 MeV to 56.2 GeV whereas the PLSEC model is an acceptable fit from 562 MeV to 17.8 GeV.The McCann model is not a good fit, with  2 values exceeding the critical value in most bins (7 d.o.f./ critical value 24.3 for  = 0.001).In contrast, only the LP model provides an acceptable fit across the whole energy range (8 d.o.f./ critical value 26.1 for  = 0.001) whilst minimising  2 compared to the other models, The LP model is therefore the preferred spectral model overall.
During the final preparation of this paper we became aware of a new MSP spectral model based on 104 MSPs (Wu et al. (2022)) prepared using the same method as Xing & Wang (2016).However, we find that this model is essentially the same as that of Xing & Wang (2016), having log-likelihood = -142.1 and AIC = -221.5,and similar  2 per bin with respect to the critical value as Xing & Wang (2016), so we do not consider it further.
spectral shape is insensitive to geometry.Torres et al. (2019) Gamma-ray spectrum insensitive to geometry but two visible emission regions needed to account for both X-ray and gamma-ray emission.Vigano et al. (2015a)Qualitatively, magnetic inclination affects spectrum via influence on radius of curvature and the consequent parallel electric field.Vigano et al. (2015b)Favourable viewing geometry sampling initial particle trajectories produces softer low-energy spectra.Cerutti et al. (2016)A proof of concept model showing spectral softening below 400 MeV for  =90°as  increases, in contrast with few spectral shape changes for  = 45 and 60°.

Figure 1 .
Figure 1.The likelihood best-fit models for LP and PLSEC (both this work), McCann, and Xing and Wang for the stacked spectrum of 118 MSPs.The LP model is the best fit overall and has just visible 1 uncertainty indicated by grey shading.The blue and red bands (each comprising two bins), indicate the overall energy range used to form the low and medium energy colours respectively as in Eqns. 5 and 6.For clarity the uncertainty of the McCann (1) and the Xing and Wang (3) models is not shown.

Figure 2 .
Figure2.The LE and ME colours for the 118 MSPs plotted against pulsar luminosity, , characteristic age and surface magnetic () field.There is no apparent correlation with LE and ME colours distributed evenly above and below the LP colour values (red and blue bands).Colour ratios are shown either as points with errors (where flux points exist in all bins for the colour), or as upper or lower limits where there are a mixture of upper limit and flux point observations in the energy bands used to form the colour ratios.

Figure 3 .
Figure 3.The LE and ME colours for Class I MSPs vs a magnetic impact angle () from inclination angles in Johnson et al. (2014) for a two-pole caustic model (TPC), and an outer-gap (OG) model.There is no apparent correlation of colour and impact angle.The LP colour values are shown as red and blue bands.Colour ratios are shown either as points with errors, or as upper or lower limits.

Figure 4 .
Figure 4.The LE and ME colours for Class II MSPs vs a magnetic impact angle () determined from inclination angles in Johnson et al. (2014) for a two-pole caustic (TPC), an outer-gap (OG) and a slot-gap (SG) model.There is no apparent correlation of colour and impact angle.The LP colour values are shown as red and blue bands.Colour ratios are shown either as points with errors, or as upper or lower limits.

Figure 5 .
Figure 5.The LE and ME colours for Class III MSPs vs a magnetic impact angle () determined from inclination angles in Johnson et al. (2014), for a pair-starved polar cap (PSPC) model.There is no apparent correlation of colour and impact angle.The LP colour values are shown as red and blue bands.

Figure 6 .
Figure 6.The LE and ME colours for ten MSPs plotted against line of sight inclination ( ), magnetic inclination () from Benli et al. (2021) and magnetic impact angle ( =Line of sight inclination − magnetic inclination).All inclinations are with respect to the pulsar rotational axis.There is no obvious correlation of colour with line of sight or magnetic inclination whilst the magnetic impact angle shows colours with symmetry between 5 and 10°w here the colour ratios are close to 1 indicating a harder spectrum.The LP colour values are shown as red and blue bands.Colour ratios are shown either as points with errors, or as upper or lower limits.

Figure C1 .
Figure C1.Individual MSP spectra fitted with a PLSuperExp2 spectral model.The examples here range from the lowest to highest detection significance (TS).

Figure C2 .
Figure C2.Individual MSP spectra fitted with a log parabola (left) and power law model (right).The examples here are those with the highest detection significance.
Bin Centre Lower Bin Energy Upper Bin Energy MSP Count and Percentage of Sample Mean Energy Flux

Table 2 .
Summary of inclination effects on gamma-ray spectra as indicated by current models, varying from no effect, to changes in model spectra as magnetic obliquity, , varies MSP spectrum as noted above.If an upper limit (UL) is present at 133 MeV or 7498 MeV, then the LE and ME colours represent upper and lower limits respectively.

Table 3 .
The parameters of the best-fit LP model (this work) using Eqn.2, for the stacked differential energy flux of 118 significantly detected MSPs in the energy range 100 MeV−56.2GeV.The LP model is the best fit overall compared to all other models considered.

Table 4 .
Magnetic impact angle () for Class I MSPs determined from inclination angles in Johnson et al. (2014) provided for a two-pole caustic model (TPC), and an outer-gap (OG) model.

Table 6 .
Magnetic Xing & Wang (2016)14)lass III MSPs determined for inclination angles inJohnson et al. (2014)provided for a pair-starved polar cap (PSPC) model. of 118 MSPs with well constrained 1 errors (Fig.1) to be more representative of an MSP ensemble than that ofXing & Wang (2016)which uses just 39 MSPs, a much greater 3 uncertainty and only

Table A1 .
The parameters of the best-fitted PLSEC model (this work) using Eqn. 1, for the stacked differential energy flux of 118 significantly detected MSPs in the energy range 100 MeV−56.2GeV.This model is inferior to the LP model above.foreach energy bin which is the difference between the stacked MSP flux and the model flux evaluated at the centre of each energy bin in the spectrum.We sum the log of the absolute value of each bin residual to obtain a set of log likelihood values, one for each normalisation.We then determine a minimum log likelihood value of −135.2, −131.9, −137.6 and −123.9 for the LP, PLSEC, Xing and Wang, and McCann models, respectively.Using Eqn. 4 and the minimum log likelihood value of each model as  0 and  1 , we determine that LP, and Xing and Wang, are significantly preferred over that of McCann at TS 22.6 (4.7) and 27.4 (5.2) respectively, whilst LP, and Xing and Wang are an equally good fit to the MSP spectrum using a likelihood analysis.

Table A2 .
A breakdown of  2 test statistic by energy bin, for stacked models of MSP emission fitted to the 118 MSP stacked flux, ranked in order of total  2 .The LP model (this work) is best in minimising  2 overall and the only model which provides an acceptable fit across the whole energy range (0.1−56.2 GeV).The McCann model is not a good fit, with  2 values exceeding the critical value in most bins (7 d.o.f./ critical value 24.3 for  = 0.001).All other models 8 d.o.f./ critical value 26.1 for  = 0.001.