CTA sensitivity to a millisecond pulsar population in the center of the Milky Way

(1)

MSc Physics and Astronomy

Gravitation and Astroparticle Physics (GRAPPA)

CTA sensitivity to a millisecond pulsar population in the center of the Milky Way

by

Harm van Leijen

10536299 Master Thesis

60 ECTS

September 2019 - September 2020 Supervisor:

prof. dr. Shin’ichiro Ando

Daily supervisor:

dr Oscar Macias

Second examiner:

(2)

Abstract

Analysis to the Fermi data shows an excess of GeV emission coming from the galactic center. The leading theory suggest that the emission from a population of millisecond pulsars (MSPs) is the cause for this excess. Such a population would, along with the prompt gamma-ray emission, also emit large quantities of electrons and positrons in there environment. These e± _{could potentially inverse-Compton scatter ambient}

pho-tons to gamma-rays that fall within the sensitivity range of the upcoming Cherenkov Telescope Array (CTA). In this thesis we asses the detection potential of CTA to these secondary emissions. To do this, we made a realistic assessment of the different back-ground sources and used these for a likelihood analyses. We predict that CTA will be able to measure the signal if the total IC luminosity from the inner 10° x 10° exceeds the threshold value of 5.4 ∗ 1035 _{erg/s in the most optimistic conditions and 9.4 ∗ 10}36

under the most unfavorable conditions.

Acknowledgements

There are a lot of people that I would like to thank for helping me with this thesis. First and foremost, I would like to sincerely thank my daily supervisor Oscar Macias. He has not only guided me through the project but has also been an academic mentor for me throughout the year. I would also like to thank my supervisor Shin’ichiro Ando for giving me the opportunity to do this project and providing me with very useful feedback. I would especially like to thank both of them for giving me the opportunity to visit a dark matter conference at the Kavli IPMU in Japan. Furthermore, I would like to thank Jacco Vink for being my second supervisor. Finally, this project has made use of the data provided by Deheng Song and Shunsaku Horiuchi, I would also like to thank them for there work.

(3)

1 Introduction

Ever since the Fermi Gamma-ray Space Telescope started surveying the sky in 2008, the available data on gamma rays has sky-rocketed. It has provided us with over 3000 new gamma-ray sources, limits on dark matter annihilation and insights to completely new structures around the galactic center.1 _{Another interesting discovery from the Fermi data is an excess of}

gamma-ray counts originating from the galactic centre.2 This excess is referred to as the galactic center

GeV excess (GCE). For a long time the leading theory was that self-annihilating dark matter particles were the source for this excess3.4 Recently however, studies have shown that a

population of millisecond pulsars could provide a more likely explanation for the excess.1,5–10

Since pulsars have long been speculated to be a source of electrons and positrons,11,12 such

a group of millisecond pulsars could, along with the prompt gamma ray emission, also emit

e± pairs into their stellar environment. These e± could inverse-Compton scatter ambient

photons to very high energies in the environment surrounding the MSPs.5 _{Recent simulations}

to MSPs supported the idea that the MSPs were capable of accelerating e± _{to TeV scale} energies.13 _{The leading observatory for gamma-rays with these energies is the High Energy}

Stereoscopic System (H.E.S.S.). H.E.S.S. surveys to globular clusters showed an extended TeV

γ-ray emission coming from the Terzan 5 Cluster,14 _{this could be explained by the secondary}

radiation from the MSPs. Although observations to 15 other globular clusters could only provide us with upper limits on the flux, the expectation is that future telescopes will give us a better view of these secondary emissions. One such telescope is the Cherenkov telescope array (CTA). CTA is expected to start operating in early 2022, after which it will provide us with unprecedented surveys to gamma-rays from the galactic center.

In this project, we propose that a population of 20-50 thousand MSPs is responsible for the galactic center excess and that this population emits highly energetic e±_{into their surrounding.} We then assess the detection capabilities of CTA to the inverse-Compton emission of such a population. To do this we will model the expected background and signal counts for CTA and use these maps for a likelihood analysis. The thesis is structured as follows: In section II we will give the necessary background information, in section III we will model the expected counts for CTA, we will then use these maps for a likelihood analysis in section IV and finally the results of this analysis will be shown in section V.

(5)

2 Theoretical Background

2.1 Galactic centre GeV excess

For most of the gamma rays that reach the earth, we have a rough understanding of where they come from. The largest part is so called galactic diffusive emission; emission that is produced when charged cosmic ray particles interact with the interstellar gas or when photons are up-scattered by electrons. There will also be a fraction coming from known point sources such as pulsars, supernova remnants or active galactic nuclei.15 _{But if we remove all these}

known gamma ray sources from the Fermi gamma ray sky-maps, we end up with a residual emission at the centre of our universe coming from a yet unknown source. This excess is known as the Galactic Center GeV Excess (GCE). In figure 1 we can see that the GCE has a roughly spherical symmetry around the galactic centre with a emissivity that drops with ∼ r−2.4 _{to ∼} 10 degrees latitude.16 _{The right plot in Fig} ₁_{shows the spectral dependency of the GCE. For}

this plot a lot of different models are used, it becomes visible that even for a wide variety of models a remarkable stable GCE pops up.17

2.1.1 Explanations for the excess

Ever since the discovery of the GCE, there has been a lot of debate about the sources of this excess. When it was first described by Goodenough and Hooper,19 _{they speculated that}

it could be produced by the the annihilation of dark matter particles into standard model particles. This explanation would have a big impact to the understanding of our universe since it would be the first non-gravitational measurements of dark matter.3 _{The main indicator that}

the excess is in fact caused by annihilating dark matter is the similarities in the spectral shape. The flux from annihilating dark matter can be described by the following formula:20

Φ(Eγ, φ) = σν 8πm2 X dNγ dEγ Z los ρ2(r)dl, (1)

Where σ is the dark matter annihilation cross section, mX the mass of the dark matter particle, ρis the dark matter density as a function of the distance from the galactic centre. Gordon and

(6)

Figure 1: Left figure: Figure from Daylan et al. (2016)18_{showing the GCE. All the known gamma}

ray sources are subtracted from the total flux which results in a residual model of the GCE. Right

figure: Figure from Carole et al. (2015)17 _{showing the spectral shape of the GCE.}

particles annihilated to a mixture of τ+_τ− _{and b¯b with an annihilation cross section in the} order of the thermal relic.7 _{This fit is shown in Fig} _2.

Some time after this dark matter interpretation of the GCE, people came with more ’conven-tional’ astrophysical explanations.9 _{They speculated that the excess could also be explained}

by a population of unresolved millisecond pulsars, since such a source would produce a spec-tral shape which is similar to the one we observe from the GCE.7–9 _{Although none are yet}

identified, their long lifetime suggest that MSPs could exist in the galactic center.21

Since the spectral shape of the excess was not enough to distinguish the dark matter from the millisecond pulsar theory, people started looking into other differences. One such difference lies in the statistical distribution of the source: Since the dark matter signal is a diffuse source, we expect a smooth signal with Poisson fluctuations between the pixels. If our signal however comes from MSP-like point sources, we rather expect to see hot- and cold-spots (see figure 3a).10These differences will cause two distinguishable flux probability density function (PDF)

of the photon counts (shown in figure 3b) since high photon counts become more likely when the sources term is clustered.

(7)

Figure 2: Figure from Gordon and Macias (2013),7 showing that dark matter annihilation’s give a reasonable fit to the GCE. The black dots with the red error bars represent the GCE with the 1 sigma error bars. The solid black line represents the self-annihilation of WIMP particles to b¯b 55 % and to leptons 45% of the time.

Figure 3: Figure from Lee et al. (2016)10 showing the difference between the dark matter and millisecond pulsar statistical distribution. The left plot shows that the counts for MSPs are more clustered and the right plot shows that this leads to an increase of the probability of high count rates.

These different flux PDFs were used to apply a Bayesian model comparisons to the data. This technique showed that the data preferred the point sources over a diffusive dark matter one.22

A somewhat similar analysis has been done in a paper by Bartels et al. (2016),1 _{where they}

(8)

Another possible difference lies in the morphology of the excess. The morphology of the GCE has been extensively investigated by Macias et al. (2018),6 _{they used template fitting to}

show that the excess could be statistically better described by the boxy and nuclear bulge, rather than a spherical symmetric model.6 _{Since dark matter is thought to follow a spherical}

symmetric NFW-profile while the stellar distribution is described by this boxy-model, this result shows that the excess is of a stellar nature.

Figure 4: Figure from Bartels et al. (2017),23 _{showing the different models that where used to fit}

the galactic center excess. The left figure is a spherical symmetric NFW-model while the right is a stellar boxy model with a nuclear bulge. They found that the data highly preferred the stellar model.

(9)

2.2 The millisecond pulsars

As mentioned before, the leading theory suggests that the source for the GCE is a population of millisecond pulsars. Pulsars are neutron stars with low rotation times(∼1 sec) and high magnetic fields. Due to the misalignment of the magnetic axis with the rotational axis, we measure a pulsating signal coming from the pulsars. Some of the pulsars are accompanied by another star in a closed binary system, when these stars accrete their mass into the pulsar, the angular momentum also gets transferred and the rotational time of the pulsars can decrease to mere milliseconds.24 _{This is why MSPs are sometimes referred to as recycled pulsars.}25

The strong electromagnetic fields can tear charged particles from the pulsars surface and drag them along the magnetic field lines. This circular motion of electrons and positrons can emit photons through Synchrotron radiation. Some magnetic field lines stretch so far away from the pulsars that the radiation would cause the electrons following these field lines to exceed the speed of light. To avoid this, these magnetic field lines do not reconnect with the pulsar but instead stream into the stellar surrounding. The e± _{that travel through these open field} lines interact with the starlight, the cosmic microwave background (CMB), galactic magnetic fields or the interstellar medium (ISM) to create the secondary radiation that we are interested in for our project. The e± _{luminosity that is produced this way is proportional to the spin} down energy of the pulsars and to the prompt gamma-ray luminosity. It has been estimated that the gamma ray efficiency Lγ/ ˙E is roughly 10 %.26 We therefore write:

Le± = f_e±E = 10f˙ _e±L_γ, (2)

The value of fe± is currently not well defined. Observations of MAGIC to the M15 globular

cluster estimated that the value would fall in the range of ∼ 0.2 − 2.0%.27 _{This however can}

not be simply extrapolated to the galactic center. In our study we investigate the constraining potential of CTA to the value of f± in the galactic center.

(10)

2.3 The Cherenkov Telescope Array

The Cherenkov Telescope Array is a multi national project that aims to detect electromagnetic waves with energies ranging from a couple GeV up to 300 TeV. These high energy photons are produced in the most violent environments of our universe and observations could help us understand the processes that are at play at these places. CTA consists of more than a 100 telescopes in both the northern(at La Palma) and the southern hemisphere(at Chile). These telescope cover a large effective area to provide CTA with a predicted sensitivity which is 10 times larger than the current gamma ray detectors.

Photons with such high energies are not able to reach the earth before reacting with particles in our atmosphere. Direct detection of gamma ray photons is therefore not possible with ground based observatories. CTA is however able to indirectly measure the incoming gamma rays by making use of the Cherenkov effect. Incoming gamma rays will hit particles from the earth atmosphere and induce a cascade of sub-atomic particles. These particles will carry so much energy that they travel faster than the speed of light in the air (99.97 % of c). This excedance of the local speed of light causes the emission of Cherenkov light. The Cherenkov light from all the cascade particles will travel in a cone and overlap on the earths surface. This is measurable by CTA and can be used to reverse calculate the energy of the initial incoming photon.

One of the main goals of CTA will be to make better observations to the galactic center. Since most of these sources are diffuse and CTA is a pointing instrument, this is not always trivial. To map the galactic center, CTA adopts a strategy in which 9 individual pointing positions, centered around the galactic center, are used.28 _{The total observation time is distributed}

equally between these points. The total image is then constructed from these individual measurements. This setup is shown in figure 5.

(11)

Figure 5: Figure from Acharyya et al. (2020)28 showing the observation strategy for CTA . The right figure shows the 9 pointing positions used to measure the galactic center. In the left figure, an extended measurement is shown where 15 additional pointing positions are added.

3 Modelling the data

In order to calculate the sensitivity, we first have to map all the types of radiation that make up the background and the signal that we try to find. We do this by splitting the background into 5 different types and model these types independently. In the sections below, we will explain how these types are modelled.

3.1 Cosmic ray background radiation

By far the most important source of gamma-ray counts are caused by cosmic ray particles that are misidentified as gamma rays. When high energy protons enter the atmosphere, they can interact with particles from the Earth’s atmosphere to create a hadronic shower. The neutral hadrons, created in this process, subsequently decay to a pair of gamma ray photons. Analysis to the shape and the arrival times of these gamma-ray showers allows CTA to efficiently discriminate these showers from gamma-ray showers.29 _{Some showers however will still be}

misidentified as gamma-ray induced. The same goes for electromagnetic showers that are produced by high energy e± _{particles. To model this type of background, we use the data} provided by CTA. They made a realistic assessment to the background rate per energy bin by modelling the incoming cosmic ray flux and subsequently checking the discrimination capability of CTA to these incoming fluxes. We assume a homogeneous exposure in our range of interests, we also assume that the cosmic ray background is a isotropic signal. This means that our

(12)

model for the residual background is just a constant flux varying only with energy. The only variations between pixels are caused by Poisson fluctuations. The rate of these background counts are displayed in figure 12.

3.2 Galactic diffusive emission

An important fraction of the high energy photons that reach the earth are actually produced by secondary emissions from cosmic ray particles. When cosmic ray particles are produced, they can propagate long distances through the galaxy until they encounter photons from the interstellar radiation field (ISRF) or particles from the interstellar medium (ISM). These interactions can produce high energy photons through three processes:

• Pion decay: Neutral pions can be produced when cosmic ray particles interact with the ISM, these neutral pions decay into two photons.

• Bremsstrahlung: When cosmic rays encounter other charged particles, they can get deflected or decelerated. This loss in energy comes in the form of an emitted photon. • Inverse-Compton scattering: photons from the ISRF can get up scattered when they

bounce into cosmic ray particles.

In our project we will combine the pion decay and bremsstrahlung into a single component called the gas correlated emission. To model these components, we use a c++ package called galprop. galprop models the propagation of charged particles by solving the transport equation for a given set of boundary condition and source terms. For our analysis we use the newest version of galprop; v56. This version allows for a more sophisticated 3D analysis of the galaxy with a 3D map of the ISRF and the neutral and molecular hydrogen. Jóhannesson et al. (2018)30 _{showed that this would lead to a significant increase in accuracy compared to}

the more simplified 2D view of galprop v54 . More information about how galprop does this can be found in appendix A. The processes mentioned above will cause the cosmic rays to loss energy during their propagation from the center to the outer regions of the galaxy. These energy losses are hard to predict and assumptions can greatly influence the model. We will therefore implement a new technique in which we divide the GDE in four regions: The first region models the radiation coming from the inner 3.5 kpc, the second region models the radiation coming from 3.5 kpc to 8 kpc, the third from 8 to 10 kpc (this is the region of our

(13)

Figure 6: The morphological maps of the galactic diffusive emission produced by GALPROP.

The maps are divided into two types of radiation (Gas correlated emission and emission caused by inverse Compton scattering ) and four Galactocentric rings. The values represent the flux (s−1sr−1cm−2T eV−1) in a 0.5° x 0.5° pixel for photons centered around 10 TeV.

solar system) and the last region models the radiation from 10 to 50 kpc. This division in rings will make our analysis less model dependant. To obtain the overall normalisation of the flux we fit the IC- and the Gas correlated data independently to the Fermi-Lat observation data of the galactic center. Figure 6 shows the 8 diffusive maps that we end up with; 4 gas correlated maps and 4 IC-maps.

(14)

3.2.1 Alternative GDE model

To prevent theoretical biases in the gas maps, we also create a second set of GDE maps; GDE model 2. In both models, the gas correlated templates are based on the distribution of atomic and molecular hydrogen in the galaxy. The HI templates are based on radio telescope surveys to the 21 cm emission lines of these HI particles.31 _{The distribution of molecular hydrogen can}

not be done directly but is instead based on the distribution of carbon monoxide. The CO is mapped by 115 GHz center for Astrophysics survey of CO by their 2.6 mm emission line.32 _It

is assumed that the distribution of molecular hydrogen closely follows this distribution. These two templates only show the gas column densities of the molecular and atomic hydrogen, the translation of this to a 3D model is where the difference between our model lies:

• Interpolation method: To determine the position of the hydrogen in our line of sight, Galprop assumes a circular motion of the gas clouds around the galactic center and then uses the red-shift of the emission lines to determine the position. The following relation can then be used to decompose the maps into galactocentric rings:

VLSR = V (R)R R − V (R) sin(l) cos(b), (3)

Where VLSR gives the radial velocity, V (R) the rotational velocity at radius R, Rgives the distance of the galactic center to our Sun and l, b are the longitude and latitude of the gas. Because the kinematic resolution is lost for latitudes close to the GC, this approach works only for |l| > 10. Within these boundaries, the densities are calculated by interpolating in each ring from the boundaries.

• The hydrodynamic approach: The boxy-buldge like shape of the galactic center causes the velocities of the gas clouds to be highly non-circular.33 _{The assumptions for the}

in-terpolation method will therefore cause some errors. To prevent this Pohl et al. (2008)33

simulated the non-circular motion of gas around the galactic center and used these ve-locities for the kinematic analysis.

Some molecular hydrogen will not be mixed well with CO and will therefore be missed. Furthermore, possible errors may arise from the assumption that the spin temperature is constant for the whole line of sight (170 K). To compensate for these deficiencies, the alternative model has two additional templates; the negative and the positive dust

(15)

residual map. These maps are obtained by mapping the thermal emission from dust and subtracting the emission already traced by the 21 cm and 2.6 mm emission maps.34

For the first GDE model, these corrections are directly implemented in the hydrogen maps. For GDE model 2 however, these templates are provided separately to provide more freedom to the model. The maps for GDE model 2 are shown in Fig 7.

Comparisons of these two models showed that the hydrodynamic approach gave a much better fit to the available Fermi-data.6,35,36 _{We will however use the maps provided by Galprop}

(16)

4 2 0 2 4 4 2 0 2 4 Latitude [deg] HI ring 0 4 2 0 2 4 4 2 0 2 4 HI ring 1 4 2 0 2 4 4 2 0 2 4 HI ring 2 4 2 0 2 4 4 2 0 2 4 HI ring 3 4 2 0 2 4 Longitude [deg] 4 2 0 2 4 Latitude [deg] H2 ring 0 4 2 0 2 4 4 2 0 2 4 H2 ring 1 4 2 0 2 4 4 2 0 2 4 H2 ring 2 4 2 0 2 4 Longitude [deg] 4 2 0 2 4 H2 ring 3 4 2 0 2 4 Longitude [deg] 4 2 0 2 4 Latitude [deg]

Dust reddening E(B-V) positive

4 2 0 2 4 Longitude [deg] 4 2 0 2 4

Dust reddening E(B-V) negative

Figure 7: The morphological maps of the gas correlated emission constructed by Macias et al.

(2018).6 The maps are divided into a molecular and atomic hydrogen map and the residual dust maps who trace the dark neutral material. The hydrogen maps are also divided into four Galactocentric rings.

(17)

3.3 The signal

For the calculation of the morphology and the spectrum of the signal maps, we use the work done by Song et al. (2019).5 They state that although the prompt gamma-ray radiation from

MSP follows the spatial morphology of the MSP, this does not necessary have to be the case for the secondary radiation since this is a convolution between the cosmic ray source term and all their possible targets (ISRF, ISM and Galactic magnetic fields). This convolution is done by an altered version of galprop in which a possible dark matter term is replaced with a MSP source term with a morphology and spectral shape described in the following sections. All non MSP cosmic ray propagations are then turned off to get a map of the secondary radiation from the MSPs only. We will follow this approach but improve the analysis by using the updated version of galprop (v56).

3.3.1 The spatial distribution of the MSPs

In our analysis we assume that the MSP source term follows the stellar distribution. This distribution is modelled in two parts: The nuclear bulge and the boxy bulge. The boxy bulge describes the inner few kpc of the galaxy which has a boxy/peanut like morphology. To model this, we use model S from Freudenreich(1998)37 which states that the density profile can be

described by: ρbar(R, φ, z) ∝ sech2(Rs) ×      1, R ≤ Rend e−[(R−Rend)/hend]2_, _{R > R} end (4) Where the parameter values are given by table 1 and Rs is described by:

RC⊥ ⊥ = (|X 0_{| /a} x) C⊥ + (|Y0| /ay) C⊥ , (5) RCsk = R C|1 ⊥ + (|Z 0_{| /a} z) C1 , (6)

The nuclear bulge describes a distinctive cluster of stars in the innermost part (∼200 pc) of the Milky Way where the star density and luminosity are relatively high. To model this we follow the work done by Launhardt et al. (2002),38 which divides the nuclear bulge even further into

the nuclear stellar cluster (NSC) and the nuclear stellar disk (NSD) whose density profiles can be given by:

(18)

ρN SC(x, y, z) =                ρ3 1 +_rr c 2−1 , for r/pc < 6 ρ4 1 +r rc 3−1 , for 6 < r/pc ≤ 200 0, for r/pc > 200 (7) ρNSD(x, y, z) =              ρ0 _r rd −0.1 exp−|z|_z d , for r/pc < 120 ρ1 r rd −3.5 exp−|z|_z d , for 120 ≤ r/pc < 220 ρ2 r rd −10 exp−|z|_z d , for r/pc ≥ 220 (8)

All these parts together (boxy bulge + NSC + NSD) have been shown to be a good fit to the galactic center excess extracted from the Fermi data.5,6,39–41 _{In our analysis we also}

look at the possibility of a spherically symmetric profile. We use this as a template for a dark matter self-annihilation signal and check whether our analysis can separate this from the stellar template. For this density profile we use a generalised NFW-profile:

ρ(r) = ρ0 _r R γ _1+r/R s 1+Rs/R (3−γ), (9)

(19)

Stellar density profile parameters Parameter value Rend 3.128 kpc hend 0.461 kpc ax 1.696 kpc ay 0.6426 kpc az 0.4425 kpc Ck 3.501 C⊥ 1.57 rd 1 pc zd 45 pc ρ0 301 M pc−3 rc 0.22 pc ρ3 3.3 x 106 M α 1 β 3 γ 1.2 R 8.25 kpc Rs 23.1 kpc

Table 1: Parameters for all the stellar density profiles from Freudenreich (1998)21and Launhardt et al. (2002)38

3.3.2 MSP e

±

injection term

For the injection spectrum of e± _{from the MSP, we assume a power law with an exponential} cutoff:5

dN dEdt ∝ E

−Γ exp (−E/E

cut) , (10)

Where Ecut gives the cut off energy and Γ the spectral slope. These values are not

perfectly known so we select a few models shown in Table 2. With this equation and the spatial distribution described in the previous section we can now define the source term as:

(20)

Model Γ Ecut (TeV) Baseline 2.0 50 Inj1 1.5 50 Inj2 2.5 50 Inj3 2.0 10 Inj4 2.0 100

Table 2: The parameters for the 5 different injection spectra that we use.

q(~r, E) = c

4πN0

dN

dEdtρ(~r), (11)

Where c/4π is a galprop convention and N0 is a normalisation factor such that the total integral over E and ~r of the equation equals the total injection Luminosity Le±. We do this

for every model from Table 2 and end up with the spectral shapes for the secondary emission that are shown in the left figure of fig 8. Fig 8 also shows the morphology of the secondary emission. 10 2 ₁₀ 1 ₁₀0 ₁₀1 ₁₀2

Energy [TeV]

1015 1014 1013 1012 1011 1010 10 9 10 8

E

2

dN

/d

E[

Te V sc m 2sr

]

IC signal spectra

inj 0 inj 1 inj 2 inj 3 inj 4 dark matter 4 2 0 2 4

Gal. Long. [deg]

4

2 0 2 4

Gal. Lat. [deg]

MSP Signal morphology

4 2 0 2 4

Gal. Long. [deg]

4

2 0 2 4

Gal. Lat. [deg]

DM signal morphology

Figure 8: The left figure shows the IC spectra for the different injection models. The middle and

(21)

3.4 Fermi bubble

Another potential gamma ray source is the Fermi bubbles(FB). The Fermi bubbles consists of two large scale structures that stretch above and below the galactic plane to about 55◦ in latitude. They are characterized by their clear spectral shape that scales with ∼ E−2 for GeV energies with a cutoff around 100 GeV.42 Recent studies have also shown that the

emission towards the base of the bubbles (galactic latitudes(b) < 10◦_{) is brighter and without} a significant energy cutoff before 1 TeV42.43

The relative symmetry between the upper and lower half of the galactic plane suggests an injection mechanism that originates from the galactic center. Several models have been proposed that would explain this injection mechanism and the formation of gamma rays, but none has yet been definitive. One model suggests that the bubbles could be formed by power-ful jets from the active galactic nucleus (AGN).44 _{Other researchers suggest that the bubbles}

could also be explained by a series of accretion events onto the Milky Way’s super massive black hole45 _{or by winds from a series of supernovae explosions.}46

To map the Fermi bubbles we use the work done by Malyshev & Dmitry (2019)47 _{and Macias}

et al. (2018).6 The first group used an image reconstruction technique to get the morphology

of the Fermi bubbles at low latitudes, the later improved these map by filling in the missing data points that were created due to the use of a point source mask template. They did this with the use of a Laplace interpolation technique which solves the Poisson equation on the missing data points with the neighbouring unmasked pixels as boundary condition. The result is shown in the right figure of Fig9. For the spectrum we follow the work from Rinchiuso et al. (2020),48 where they extrapolate the Fermi-data49 to energies above 1 TeV. This gives us two

models: an optimistic (FB min) and a conservative (FB max) one. These models are chosen such that all detection of the FB fall within the two models and the H.E.S.S. measurements of the Galactic ridge are not overshot. The parameters of these two models are shown in Table 3 and the resulting spectral lines are shown in the left figure of Fig 5.

(22)

Model φ0

h

TeV−1cm−2s−1sr−1i Γ Ecut[TeV]

FB max 1 × 10−8 1.9 20

FB min 0.5 × 10−8 1.9 1

Table 3: The parameters of the two models from Rinchiuso et al. (2020)48

Figure 9: The spectral shape and the morphology of the Fermi bubbles that we use for our project.

3.5 High energy point sources

Since our analysis focuses on the center of the galaxy, there are some point sources that could potentially overshine the background models. Modelling these point sources is hard and could lead to inaccuracies. We therefore decided to completely mask the points with known point sources. Since we are interested in high energy gamma-rays, we use the data from the third Fermi-LAT Catalog of High-Energy Sources (3FHL).50 _{We filter this list for all point sources}

that fall into our range of interest and do not have a spectral cut-off, we then create a mask for all pixels within 0.25° of these points. We also decided to mask the inner 0.5° degree of the galactic plane since these have big uncertainties too. Figure 10 shows the mask matrix that we then get for our model.

(23)

Figure 10: Matrix that we use to mask the point sources and the galactic plane.

3.6 Instrument response function

In order to get the expected counts for CTA, the incoming flux has to be convolved with the instrument response function(IRF). For our analysis we use the CTA IRF data that has been made publicly available by the CTA-collaboration. They obtained this data by performing detailed Monte-Carlo simulations in which simulated air-showers are produced which are used for further processing.51 _{This IRF consists of three parts:}

• The effective area: The total collection area of the detectors as a function of the incoming energy. This effective area is shown in the left figure of Fig.11 which shows that the effective area of CTA can be as high as 4.5 million m2_.

• The energy dispersion function: The energy dispersion gives the probability of measuring an energy Ereco as a function of the true incoming energy Etrue. The middle figure of

fig. 11 shows that photons with higher energies have a lower dispersion compared to the lower energies.

• The point spread function: The point spread function describes the probability of mea-suring a position pme as a function of the true incoming position ptrue. The right figure

of fig. 11shows that this PSF is very small in the range of energies that we are interested in (0.01-100 TeV). Since the PSF for CTA is even smaller than the pixel size that we use(0.5◦ _{x 0.5}◦_{), it can be neglected.}

(24)

Figure 11: The three parts of the Instrument response function. These plots are produced with the

python package Gammapy.

Fig. 11shows that all parts off the IRF depend on the offset of the observation. The offset gives the angle between the line of sight of the telescope and the area that is being observed. For our analysis we assume the most favorable observation conditions which are achieved with an on-axis observation(offset = 0). For the optimal conditions, we will also assume an observation time of 500 hours with the southern CTA site at a mean zenith angle of 20°.28

The southern site is chosen because the galactic center is mostly visible from this cite. This data is adopted from the CTA-Performance-prod3bv1-South-20deg-average-50h.root file.

For the range of interest, we decided to focus our analysis on the inner 10° x 10° of the galaxy. This region is divided in 400 (0.5° x 0.5°) pixels, this could not be further decreased since the computational power was limited.

To work with the IRF FITS-files, provided by the CTA-collaboration, a python package called Gammapy is used.52 The Gammapy code is build on Numpy, Scipy and Astropy and

is specifically created to analyse the data from gamma-ray telescopes. Mohrmann et al. (2019)53showed that Gammapy was capable of perfectly reproducing the results from H.E.S.S.

collaboration, indicating that Gammapy is a trustworthy program. We can now calculate the expected count rate for the kthenergy, ithlongitudinal and jthlatitudinal bin with the function:

Nijk = Tobs Z ∆Ek dEγ Z ∞ 0 dE_γ0 φ S ij dEγ Aγ_{ef f}(E0)D(Eγ, E 0 γ), (12) Here Aγ

ef f is the energy dependant effective area, D(Eγ, E

0

(25)

and φS

jk is the incoming flux. The whole function is integrated over the energy bin ∆E and

multiplied by the total observation time Tobs. The two different energies stand for the true

incoming energy: E0

γ and the reconstructed energy:Eγ. If we now put in the fluxes of the

E

₁₀

2

Energy[TeV]

10

IC Gas correlated CR background Fermi bubbles

Figure 12: The left plot shows the flux rates for the high energy background radiations that we use

in this project. The plot also shows the Fermi-LAT data. The spectra produced by Galprop are normalised to coincide with this data. The right plot shows the counts that CTA will measure for these background sources. These are the counts measured in our ROI (10° x 10°) for a 500 hour measurement in the most favorable settings.

(26)

4 The likelihood analysis

Once we have all the count maps, we can get to the likelihood analysis. For our analysis, we define two models; the null hypothesis, which contains only background maps and the alternative hypothesis, which also has a signal component:

θ0,ijk= X n xnBn,ijk, (13) θ1,ijk = X n xnBn,ijk+ xsSijk, (14)

Where xn gives a set of normalisation factors (one for each unique background map),

P

nBn sums over all the types of background and S gives the expected counts from a MSP

source term. These models predict the expected counts for spatial bin i,j and energy bin k. Since CTA is not in operation yet, there is no data to compare these models with. This is solved by creating ’synthetic’ data; all background maps are convolved with the IRF and for each pixel we draw a number from a poisson distribution were the mean is given by the model. The mock data is then created by summing over all these maps and adding a varying MSP source term. An example is given in appendix B. The total likelihood of Model M1 and

M2 to the data can now be calculated by:

L =Y

ijk

pijk, (15)

Where the probability is described by the Poisson distribution function:

pijk =

θijknijke−θijk nijk!

, (16)

Where θijk describes the counts predicted by the model and nijk the counts in the data. To

calculate the maximum likelihood of θ0 and θ1 to the data we use the python package Iminuit. Iminuit finds the values of xnfor us that maximize equation15. This is a very computationally

demanding calculation and to speed this up, we use the quadratic approx of the log-likelihood:

− 2 log L ≈X k X ij (nij − θij) 2 σ2 , (17)

(27)

Where σ gives the standard deviation of the data map. To give this analysis even more freedom, we divide our model and the data into 11 logarithmically spaced energy bins and determine the values of xnthat maximize the likelihood for all these energy bins independently.

This approach gives the model freedom to vary the spectral shape of each background and signal map. To evaluate the significance of the MSP hypothesis, we use a likelihood ratio test. This test is defined as:

T S = −2(ln L0 − ln L1), (18)

Where high values of TS indicate a larger detectable discrepancy between the two models. These values of TS are calculated for each energy bin, to get the total TS, we can simply sum over all the energy bins. For our analysis we apply a Monte Carlo technique where we create 50 independent synthetic data maps and calculate the test statistic for each of these maps. The final TS value is then determinate by the mean of all the values and the error by the standard deviation from the mean. To convert this TS value to a p-value we use the mixture distribution formula:6 p(TS) = 2−n   δ(TS) + n X k=1    n k   χ 2 k(TS)   , (19)

Where δ is the Dirac-delta function, χ2

k is a χ2 function with k degrees of freedom,

_n

k

is the binomial coefficient and n gives the additional degrees of freedom. Since our analysis is divided in 11 energy bins, each with 1 additional normalisation factor, we have 11 additional degrees of freedom for the alternative hypothesis. We can then use this p-value to calculate the n-σ significance in the following way:6

Number of σ ≡√InverseCDFχ2₁, CDF[p(TS), ˆTS], (20)

Where (InverseCDF) gives the inverse cumulative distribution function. From this we can calculate that a 5σ detection, in our case, is acquired with a TS of 41,1. We can now increase the MSP source term in the data until our pipeline reaches this threshold value.

(28)

5 Results

For the analysis, we varied the signal component in the mock data and investigated the impact on the detection capability of CTA. The results are represented in Fig 13and Fig 14. In Fig 13we assume a perfect background model while Fig14assumes a mismodelling scenario. The left plots of these figures display the TS values plotted against the injected signal strength of the five different injection models that we use, the detection threshold of 5σ is shown as a green-dotted line. The 5 different injection spectra produce significantly different results, injection model 1 for example seems to have a detection threshold luminosity that is a factor ∼ 3 lower than that of injection model 2. The most likely explanation for this is that signals that mirror the spectral shape of the background (Inj model 2, with Γ = 2.5) are much harder to detect compared to signals which have a distinguishable spectral shape (Inj model 1, with Γ = 1.5). In the right plots of these figures, the injected luminosity is plotted against the extracted luminosity. this plot shows how capable the pipeline is in extracting the injected signal. It becomes clear that for low injection strengths, the error on the extracted signal is relatively large, but that this error decreases drastically once the injection signal increases. If we compare Fig 13and Fig 14 we see that the mismodelling scenario only increases the error on the TS-values and the extracted luminosities slightly. The detection luminosity is also just slightly larger, indicating the robustness of the pipeline.

Unfortunately, the morphological differences between the four galactrocentric IC emission maps were not significant enough for the pipeline to distinguish them in the data. This was because most of these differences lay in the parts that were further away from the galactic center. These parts were however lost due to our choice of the ROI. We therefore decided to combine the four rings into one single map again and only keep the gas correlated emission maps split into four galactocentric rings.

(29)

1035 ₁₀36 ₁₀37 ₁₀38 100 101 102 103 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 1035 1036 1037 1038 Recovered Luminosity[erg/s] Baseline 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 100 101 102 103 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 1035 1036 1037 1038 Recovered Luminosity[erg/s] Inj1 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 100 101 102 103 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 1035 1036 1037 1038 Recovered Luminosity[erg/s] Inj2 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 100 101 102 103 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 1035 1036 1037 1038 Recovered Luminosity[erg/s] Inj3 5 significance 1035 ₁₀36 ₁₀37 ₁₀38

Injected Luminosity [erg/s]

100 101 102 103 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38

1035 1036 1037 1038 Recovered Luminosity[erg/s] Inj4 5 significance

FB

min

, Perfect GDE

Figure 13: The left figures show the TS values plotted against the injected signal strength for a

perfect background model and the F Bmin model. These results are obtained from 50 Monte Carlo simulations. The dark blue line represents the mean of the result while the light blue region gives the error on this mean. The 5σ detection threshold is shown as a green dotted line. The right plot shows the extracted signal strengths for the injected signals. The diagonal red line represents the perfect case in which the extracted signal matches the injected signal perfectly. The blue patch shows the

(30)

N/

dE

[T

eV

cm

2

s

1

sr

1

]

E

2

d

N/

dE

[T

eV

Energy [TeV]

[T

eV

d

N/

dE

d

N/

dE

[T

eV

Energy [TeV]

[T

eV

cm

2

s

1

sr

1

]

injected signal 4_{extracted signal}

FB

min

, mismodeling of the GDE

(33)

Baseline Inj-1 Inj-2 inj-3 inj-4

F Bmin perfect GDE 1.3*1036 5.4*1035 3.3*1036 2.7*1036 1.0*1036

F Bmin mismodelling 1.8*1036 7.2*1035 3.4*1036 2.8*1036 1.34*1036

F Bmax perfect GDE 7.1*1036 7.5*1036 5.4*1036 6.1*1036 7.2*1036

F Bmax mismodelling 9.0*1036 9.4*1036 6.8*1036 7.8*1036 9.1*1036

Table 4: Summary of the detection threshold luminosities for all the different scenarios we use. The

values are in erg/s.

Baseline Inj-1 Inj-2 inj-3 inj-4

F Bmin perfect GDE 5.0 % 2.0 % 12.7 % 10.4 % 3.9 %

F Bmin mismodelling 6.9 % 2.7 % 13.1 % 10.8 % 5.1 %

F Bmax perfect GDE 27.4 % 28.9 % 20.8 % 23.5 % 27.7 %

F Bmax mismodelling 34.7 % 36.2 % 26.2 % 30.1 % 35.1 %

Table 5: Similar results as Table 4 but here we show the detection lower limits on the electron acceleration fraction f_e± of equation 2.

5.1 Discrimination of dark matter

We also checked whether the pipeline could successfully distinguish a MSP source term from a dark matter one. To do this, we followed the same procedure as before but instead replaced the two models in equation 13 and 14with:

θ0,ijk = X n xnBn,ijk+ xsSijk, (21) θ1,ijk = X n

xnBn,ijk+ xsSijk+ xdmSdm,ijk, (22)

Where Sdm,ijk gives a dark matter source term and xdmits normalisation factor. The mock

data still only contains the background and MSP source term. We can then use this altered pipeline to check whether the additional dark matter template gets extracted from the data. The results of this analysis are represented in Figure 17. In the left figure, the Test statistic is plotted against the injected signal. At no point does the TS value exceed the 5σ detection threshold, indicating that the additional template does not provide a significantly better fit

(34)

to the data. The right plot shows how much of both templates are extracted from the data. Although there seems to be a small overlap, the MSP template is clearly favored. This shows that CTA is capable of distinguishing a dark matter signal template from a MSP one. The same analysis was done with a mismodelled background, these results are shown in figure 18. This plot shows that even when the background is mismodelled, the same results are obtained.

1036 ₁₀37 ₁₀38

Injected Luminosity [erg/s] )

100

101

Test statistic

Detection threshold (5 )

1036 ₁₀37 ₁₀38

Injected Luminosity [erg/s]

1036 1037 1038

Recovered Luminosity[erg/s]

injected signal recovered MSP signal Recovered DM signal

FB

min

, Perfect GDE

Figure 17: Plot showing the discrimination potential of CTA for dark matter and MSP templates.

The left plot shows how much the fit improves when a dark matter template is added and the right plot shows how much of the dark matter and the MSP template is extracted from the data.

1036 ₁₀37 ₁₀38

Injected Luminosity [erg/s] )

100

101

Test statistic

1036 ₁₀37 ₁₀38

Injected Luminosity [erg/s]

1036 1037 1038

Recovered Luminosity[erg/s]

injected signal recovered MSP signal Recovered DM signal

FB

min

, mismodeling of the GDE

(35)

6 Conclusions

In this thesis we aimed to find the detection capability of CTA to a population of MSPs in the GC, the threshold luminosities that we found are shown in table 4. The numbers clearly show the large differences between the models. Especially the decision to use the Fermi max models greatly influences the results. We also translated these results to the detection limits of the electron acceleration fraction (fe±) that we introduced in equation 2. To do this, we

assumed that the prompt gamma ray emission from the MSPs (Lγ) is fully responsible for

the GCE(Lγ = 2.59 * 1037)35 and used the values of Le± for which the detection threshold

value is reached. The threshold values of fe± that we found are shown in table 5. Comparing

these values with limits on the electron efficiency fraction estimated for the MSPs in M15, suggests that CTA will only be able to measure the secondary radiations under very specific conditions. We have however been conservative with our analysis and the limits on fe± are

currently not well known, so even under the assumption of high background uncertainties, CTA has the potential to robustly discover a new population of MSPs in the GC, if the electron acceleration efficiency is in the range of approximately 25% to 35%.

We also investigated whether CTA is capable of distinguishing a spherical symmetric dark matter template from a stellar boxy bulge on. We showed that CTA is very capable of discriminating these two templates. Although there is no theory for self-annihilating dark matter that would both explain the GCE and the emission that is measurable by CTA, these results could potentially be used to further proof the stellar nature of the GCE.

(36)

Appendix A: galprop

To model the secondary emissions, we use a tool called galprop; a publicly available code that models the propagation of charged particles in the galaxy. The code was originally written by Andrew W. Strong and Igor V. Moskalenko in Fortran54 _{but later updated by other}

reasearches to C++. These people are still keeping the package up to date with the latest scientific findings in the galactic parameters. To model the propagation, Galprop solves the transport equation for a given set of boundary conditions and source terms. This transport equation can be given by:

∂ψ ∂t = q(~r, p) + ~∇ · (Dxx~∇ψ − ~V ψ) + ∂ ∂pp 2_D pp ∂ ∂p 1 p2ψ − ∂ ∂p ˙ pψ −p 3(~∇ · ~V )ψ − 1 τf ψ − 1 τr ψ, (23) Where φ(~r, p, t) gives the density per unit of total particle momentum which depends on the following factors:

• The source term q(~r, p): The source term gives the primary injected term of cosmic rays as a function of space and momentum.

• A diffusion and convection term ~∇·(Dxx~∇ψ − ~V ψ): This term describes the convection

(V) and the diffusion term as a function of the spatial diffusion coefficient Dxx = βD0

_ρ

ρ0

δ

. Where ρ is the rigidity of the electrons and β = v/c. This diffusion is caused by random scattering of particles on magnetohydrodynamic waves or discontinuities55

and can be empirically estimated from the boron-to-carbon ratio.? • A reacceleration term ∂_∂_pp2_D

pp_∂∂_p_p12ψ: The scattering of particles can also lead to

reac-celeration of these particles which is represented here as the diffusion in momentum space.

• An energy loss term ∂_∂_ph

˙

pψ − p₃(~∇ · ~V )ψi: Several processes can cause energy losses

(Synchrotron radiation,Bremsstrahlun and Inverse-Compton scattering) which are ac-counted for in this term where ˙p gives the rate of this loss.

(37)

• a particle loss term 1

τfψ −

1

τrψ: The loss of particles is also accounted for in this term

where τf gives the mean lifetime before fragmentation, and τr the mean lifetime before

radioactive decay.

To run the analysis, we first have to add a lot of physical parameters and models. Our choice for most of these parameters is based on the work done by Johannesson (2018),30_where

they find the best fit parameters of the 3D galprop model. A few of the most important parameters are shown in Table 4. For the 3D ISRF model, we chose to use model F98, which is based on the work done by Freudenrech (1998).37 We run these simulations in 200 x 200 x

100 pc3 _{galactic pixels, this was as low as the computational power would allow.}

Parameter value Xh [kpc] ±20.00 Yh [kpc] ±20.00 Zh [kpc] ±6.00 ∆X [kpc] 0.2 ∆Y [kpc] 0.2 ∆Z [kpc] 0.1 D0,xx [1028 cm2s−1] 2.27 δ 0.33 VAlfven [km s−1] 5.26 Vconv [km s−1] 0 dVconv/dz [km s−1 kpc−1] 10.00 R1 [GV] 11.57 γ0 1.90 γ1 2.39 γ0,e 1.60 γ1,e 2.43 γ2,e 4.00 R1,e [GV] 2.18 R2,e [GV] 2.17 × 103

Table 6: Some of the parameters that we use to run the galprop code. Xh, Yh, Zh give the size of the CR halo that we are interested in and ∆X, ∆Y , ∆Z give the pixel size that we use.

(38)

Appendix B: Additional results

The mock data

Figure 19: This figure shows the sum of all the background models. This is what the data would

look like without the Poisson noise.

Figure 20: Two example maps of the mock data. The points that are masked are clearly visible as

(39)

Signal extraction per energy bin

in Fig 13 and Fig 14 the total extracted signal is plotted against the injected one. This was obtained by summing over all the different energy bins. Below we plot the extracted signal strengths for each energy bin independently. It becomes visible that the pipeline works fine except for the highest energy bin. This is to be expected since the fluxes are low here and the signal count rate at these energies is close to zero.

Figure 21: Extracted signal strength plotted against the injected one for each energy bin

(40)

FB

max

results

The results obtained in Section V where done with the F Bmin model, below we present the

result for the same analysis but instead we used the F Bmax model.

1035 ₁₀36 ₁₀37 ₁₀38 100 101 102 103 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 1035 1036 1037 1038 Recovered Luminosity[erg/s] Baseline 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 100 101 102 103 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 1035 1036 1037 1038 Recovered Luminosity[erg/s] Inj1 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 100 101 102 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 1035 1036 1037 1038 Recovered Luminosity[erg/s] Inj2 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 100 101 102 103 Test statistic 5 significance 1035 ₁₀36 ₁₀37 ₁₀38 1035 1036 1037 1038 Recovered Luminosity[erg/s] Inj3 5 significance 1035 ₁₀36 ₁₀37 ₁₀38

FB

10 1 ₁₀0 ₁₀1 ₁₀2 10 13 10 12 10 11 10 10 10 9 10 8 10 7 injected signal 1 extracted signal 10 1 ₁₀0 ₁₀1 ₁₀2 10 13 10 12 10 11 10 10 10 9 10 8 10 7 E 2 d N/ dE [T eV cm 2 s 1 sr 1] injected signal 2 extracted signal 10 1 ₁₀0 ₁₀1 ₁₀2 Energy [TeV] 10 13 10 12 10 11 10 10 10 9 10 8 10 7 injected signal 3 extracted signal 10 1 ₁₀0 ₁₀1 ₁₀2 Energy [TeV] 10 13 10 12 10 11 10 10 10 9 10 8 10 7 E 2 dN/ dE [T eV cm 2 s 1 sr

1] injected signal 4_{extracted signal}

FB

max

, mismodeling of the GDE

Figure 26: The same figures as fig15but for these results we mismodelled the background an used the F Bmax model.

(44)

1036 ₁₀37 ₁₀38 Injected Luminosity [erg/s] )

100

101

Test statistic

1036 ₁₀37 ₁₀38 Injected Luminosity [erg/s]

1036 1037 1038 Recovered Luminosity[erg/s] injected signal recovered MSP signal Recovered DM signal

FB

max

, mismodeling of the GDE

Figure 27: Same figure as Fig. 17 but for these results we mismodelled the background and used the F Bmax model.

(45)

References

[1] R. Bartels, S. Krishnamurthy, and C. Weniger, “Strong support for the millisecond pul-sar origin of the galactic center gev excess,” Physical Review Letters, vol. 116, no. 5, p. 051102, 2016.

[2] D. Hooper and L. Goodenough, “Dark matter annihilation in the galactic center as seen by the fermi gamma ray space telescope,” Physics Letters B, vol. 697, no. 5, pp. 412–428, 2011.

[3] R. K. Leane and T. R. Slatyer, “Dark matter strikes back at the galactic center,” arXiv

preprint arXiv:1904.08430, 2019.

[4] P. Agrawal, B. Batell, D. Hooper, and T. Lin, “Flavored dark matter and the galactic center gamma-ray excess,” Physical Review D, vol. 90, no. 6, p. 063512, 2014.

[5] D. Song, O. Macias, and S. Horiuchi, “Inverse compton emission from millisecond pulsars in the galactic bulge,” Physical Review D, vol. 99, no. 12, p. 123020, 2019.

[6] O. Macias, C. Gordon, R. M. Crocker, B. Coleman, D. Paterson, S. Horiuchi, and M. Pohl, “Galactic bulge preferred over dark matter for the galactic centre gamma-ray excess,”

Nature Astronomy, vol. 2, no. 5, pp. 387–392, 2018.

[7] C. Gordon and O. Macias, “Dark matter and pulsar model constraints from galactic center fermi-lat gamma-ray observations,” Physical Review D, vol. 88, no. 8, p. 083521, 2013. [8] O. Macias and C. Gordon, “Contribution of cosmic rays interacting with molecular clouds

to the galactic center gamma-ray excess,” Physical Review D, vol. 89, no. 6, p. 063515, 2014.

[9] K. N. Abazajian and M. Kaplinghat, “Detection of a gamma-ray source in the galactic center consistent with extended emission from dark matter annihilation and concentrated astrophysical emission,” Physical Review D, vol. 86, no. 8, p. 083511, 2012.

[10] S. K. Lee, M. Lisanti, and B. R. Safdi, “Distinguishing dark matter from unresolved point sources in the inner galaxy with photon statistics,” Journal of Cosmology and Astroparticle

(46)

[11] A. Atoyan, F. Aharonian, and H. Völk, “Electrons and positrons in the galactic cosmic rays,” Physical Review D, vol. 52, no. 6, p. 3265, 1995.

[12] T. Erber, “High-energy electromagnetic conversion processes in intense magnetic fields,”

Reviews of Modern Physics, vol. 38, no. 4, p. 626, 1966.

[13] C. Guépin, B. Cerutti, and K. Kotera, “Proton acceleration in pulsar magnetospheres,”

Astronomy & Astrophysics, vol. 635, p. A138, 2020.

[14] A. Abramowski, F. Acero, F. Aharonian, A. Akhperjanian, G. Anton, S. Balenderan, A. Balzer, A. Barnacka, Y. Becherini, J. B. Tjus, et al., “Search for very-high-energy

γ-ray emission from galactic globular clusters with hess,” Astronomy & Astrophysics,

vol. 551, p. A26, 2013.

[15] S. Abdollahi, F. Acero, M. Ackermann, M. Ajello, W. Atwood, M. Axelsson, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, et al., “Fermi large area telescope fourth source catalog,” The Astrophysical Journal Supplement Series, vol. 247, no. 1, p. 33, 2020. [16] F. Calore, I. Cholis, C. McCabe, and C. Weniger, “A tale of tails: dark matter

interpreta-tions of the fermi gev excess in light of background model systematics,” Physical Review

D, vol. 91, no. 6, p. 063003, 2015.

[17] F. Calore, I. Cholis, and C. Weniger, “Background model systematics for the fermi gev excess,” Journal of Cosmology and Astroparticle Physics, vol. 2015, no. 03, p. 038, 2015. [18] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden, S. K. Portillo, N. L. Rodd, and T. R. Slatyer, “The characterization of the gamma-ray signal from the central milky way: A case for annihilating dark matter,” Physics of the Dark Universe, vol. 12, pp. 1–23, 2016. [19] L. Goodenough and D. Hooper, “Possible evidence for dark matter annihilation in the inner milky way from the fermi gamma ray space telescope,” arXiv preprint arXiv:0910.2998, 2009.

[20] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden, S. K. Portillo, N. L. Rodd, and T. R. Slatyer, “The characterization of the gamma-ray signal from the central milky way: a compelling case for annihilating dark matter,” arXiv preprint arXiv:1402.6703, 2014.

CTA sensitivity to a millisecond pulsar population in the center of the Milky Way

MSc Physics and Astronomy

Gravitation and Astroparticle Physics (GRAPPA)