Gamma Rays as Probes of Cosmic-Ray Diffusion Throughout the Galaxy

(1)

GRAPPA

Master Thesis

Gamma Rays as Probes of Cosmic-Ray Diffusion

Throughout the Galaxy

by

Mart Pothast

6079199

July 2018 60 ECTS

Research carried out between September 2017-July 2018

Supervisor Examiner Daily supervisors

Christoph Weniger Shin’ichiro Ando Daniele Gaggero

(2)

(3)

Diffuse gamma-ray emission attributed to decays of neutral pions can be used to probe the diffusion of cosmic rays in the Galaxy. Recent studies have shown that the spectrum of this diffuse pion emission follows a progressive hardening towards the inner Galaxy, which could imply deviations from the standard treatment of cosmic-ray diffu-sion. These results follow from a careful modelling of the gamma-ray sky but should be treated with caution. In this thesis the hardening was studied with a novel tool called

SkyFACT that can account for imperfect gamma-ray models through the use of a large

number of nuisance parameters. A similar trend in the diffuse pion emission was found with respect to previous studies, also while accounting for possible systematics in the inverse Compton emission. Moreover, at high energies the trend in the hardening is also present, possibly distinguishing between models that were proposed to explain this phenomenon. To assess the possibility of contamination by unresolved point sources a dedicated population synthesis code was implemented. Adapting a reasonable amount

of hard spectrum point sources will result in ∼ 10% of the diffuse emission in the inner

Galaxy to originate from unresolved point sources at 10–100 GeV, making it unlikely that this could significantly impact the observed hardening.

(4)

(5)

I am profoundly grateful to my supervisor Christoph Weniger for making this project possible and encouraging me throughout this year. I want to thank Shin’ichiro Ando for agreeing to be the second examiner of this thesis. Also, a special thanks to my daily supervisors: Daniele Gaggero and Emma Storm for all the discussions on cosmic rays, gamma rays and life. I thank SURFsara (www.surfsara.nl) for the use of the Lisa Compute Cluster. And last but not least: thanks to Kim, my friends, family and everyone in C4.172a for their support.

(6)

(7)

CR Cosmic Ray

UHECR Ultra High Energy Cosmic Ray

QLT Quasi-Linear-Theory

ISM InterStellar Medium

Fermi-LAT Fermi Large Area Telescope

DSA Diffusive Shock Acceleration

CMZ Central Molecular Zone

DGE Diffuse Gamma-Ray Emission

ICS Inverse Compton Scattering

ISRF InterStellar Radiation Field

IGRB Isotropic diffuse Gamma-ray Background

DNM Dark Neutral Medium

ROI Region Of Interest

PSC Point Sources

CDF Cumulative Distribution Function

SNR SuperNova Remnant

PWN Pulsar Wind Nebula

GCE Galactic Centre Excess

(8)

(9)

Introduction

The Earth is constantly being bombarded by charged particles from space. So called

cosmic rays have been measured with energies from a few GeV up to 1011GeV, an energy

impossible to achieve with any conceivable particle accelerator. In fact, particle physics started with the study of cosmic rays and particles such as the positron, the muon and the pion were originally discovered through the study of cosmic rays [1–3].

The origin of cosmic rays—where do they come from and how do they reach us?— has been stimulating scientists for more than a century since it was established in the early 1900s that cosmic rays originate from outside the Earth. A discovery that awarded Victor Hess the Nobel Prize in 1936 [4] but was in fact based on work by many scientists [5]. Today, many physics and astrophysics communities are interested in cosmic rays. While there is still much research on their origin (see for example [6]); cosmic rays can additionally serve as a messenger for highly energetic astrophysical events (for example neutrinos form supernova explosions [7]), they can provide clues to dark matter physics [8] and they can serve as a probe for plasma physicists that are interested in the motion of particles in magnetic fields [9].

The population of cosmic rays can be roughly divided as follows: Starting at the low

end of the energy spectrum from 1 GeV up to_{∼ 1 PeV, the spectrum of cosmic rays can}

be approximated by a single power-law. These particles are thought to originate from the Milky Way and supernova remnants are hypothesized to be their main accelerator, although it has not been proven that they can accelerate particles up to PeV energies [10]. At several PeV the spectrum softens. This feature is called the knee and could reflect a cut-off in Galactic accelerators but this section of the cosmic-ray spectrum still has many

open questions [11]. At 1018 eV the spectrum hardens again, creating an ankle shaped

characteristic. Cosmic rays with energies higher than this are accordingly labelled as ultra high energy cosmic rays (UHECR) and while their origin is still unknown [12], it is unlikely that they originate from our own Galaxy as is easy to show: Consider a proton

with an energy E = 1018_{eV, the Larmor radius R}

L= E/qB of such a particle, assuming

the magnetic field in the Galaxy B ∼ 1 µG, is of the order of kiloparsecs. Comparing

this to the height of the Milky Way halo, also of order kpc, it is clear that higher energy cosmic rays cannot be confined to our own Galaxy and that UHECRs are extragalactic.

(12)

smaller percentages of heavier elements [13]. The chemical composition follows roughly the abundances in the interstellar gas with some interesting exceptions. Elements such as Lithium, Beryllium and Boron are relatively much more abundant in cosmic rays than in the interstellar gas. These elements are not generated in stellar nucleosynthesis but do occur as spallation products from cosmic rays. From a careful measurement of these secondary particles it can be shown that cosmic rays traverse a certain amount of matter (grammage) corresponding to traversing the Milky Way several times. This is in fact a sign that cosmic rays undergo a diffusive motion.

Cosmic-ray diffusion throughout the Galaxy is assumed to be due to magnetohy-drodynamic (MHD) turbulence in the interstellar medium (ISM). This can be modelled with quasi-linear theory (QLT) assuming there is a turbulent magnetic field δB that is

much smaller than the large scale Galactic magnetic field B0 [6, 14, 15]. The motion of

cosmic rays can then be described by a transport equation that depends on a diffusion

coefficient that scales with rigidity D_{∝ ρ}δ _{[16]. This phenomenological model has been}

sufficient in reproducing the available data for a long time [17]. There are however many reasons to go beyond this simplistic view (see for example [9, 18, 19]).

The study of cosmic rays by several ground and space based experiments such as the Piere Auger Observatory or the Alpha Magnetic Spectrometer (AMS-02) has provided a great deal of knowledge on the nature of cosmic rays. However, the study of cosmic rays is faced with two problems: (1) the arrival direction of charged particles is randomized by scattering of magnetic fields in the Galaxy so that any spatial information is lost and (2) the only measurements of cosmic rays are taking place here on Earth. One solution for both these problems is gamma-ray astronomy.

When cosmic rays traverse clouds of gas in the Galaxy they can subsequently create

pions. The neutral π0 _{meson decays into two gamma rays which are not deflected from}

their origin. This source of gamma rays is the main constituent of the diffuse Galactic gamma-ray emission as seen by the Fermi Large Area Telescope (Fermi-LAT) [20]. This diffuse emission allows us to study cosmic rays throughout the Galaxy, giving important insight into their diffusive motion.

This thesis deals with the progressive hardening of the cosmic-ray spectrum in the inner Galaxy. Recently, two independent analyses of Fermi-LAT data have found a spectral hardening (a spectrum that is less steep, i.e. relatively more photons with higher energies) in gamma rays resulting from cosmic-ray interactions when moving towards the Galactic centre [21, 22]. This gradient in the cosmic-ray slope cannot be explained by the standard treatment of cosmic-ray propagation with a constant diffusion coefficient and this anomaly could give us a new view on Galactic cosmic-ray physics.

In [23] a non-linear propagation model was proposed where cosmic rays scatter on self-generated turbulence resulting in a harder proton spectrum at locations with higher cosmic-ray densities, i.e. at the sources. Along with an exponentially decaying magnetic field in the Galaxy this model reproduced the progressive trend of the spectral hardening for cosmic rays with an energy of 20 GeV.

(13)

A different approach was put forward by [24]. The authors proposed a phenomeno-logical model where the rigidity dependence of the diffusion coefficient is spatially

depen-dent D _{∝ ρ}δ(R)_{, resulting in a steeper spectrum at large Galactocentric radii R. More}

recently, [25] explained the hardening in the context of anisotropic cosmic-ray transport in the Galactic magnetic field. The harder cosmic-ray spectra in the inner Galaxy are due to parallel escape along the poloidal component of the magnetic field.

It should be noted that modelling the diffuse gamma-ray emission is notoriously dif-ficult and relies on cosmic-ray propagation codes such as DRAGON and GALPROP that in turn rely on gas column densities as measured in other wavelengths and the interstel-lar radiation field (ISRF) model. This makes gamma-ray analyses uncertain and there are several effects that could mimic the spectral hardening of cosmic rays as observed through gamma rays. Most importantly, the inverse Compton emission is not well con-strained and erroneous modelling could mimic a hardening. Furthermore, unresolved point sources could turn up in the diffuse emission and fake a harder spectrum.

Recently a new tool to analyse Fermi-LAT data called SkyFACT was presented [26]. This code accounts for the non-perfect models of diffuse emission by adding a large number of nuisance parameters. This makes it an excellent tool to study the cosmic-ray slope and its uncertainties across the Galaxy. SkyFACT was also recently used to study the Fermi-LAT GeV excess [27].

In this thesis I will discuss the diffuse gamma-ray emission as observed with Fermi-LAT and the implications for cosmic-ray physics. The progressive hardening of the cosmic-ray spectrum in the inner Galaxy was analysed using SkyFACT and its robustness was assessed by carefully analysing the ICS component and estimating the number of unresolved point sources. Moreover, if the hardening is also present at higher energies this would allow to discriminate between the above mentioned models.

This thesis is subdivided as follows: chapter 2 and chapter 3 will provide the necessary theoretical background on, respectively, cosmic rays and gamma rays and introduce the

SkyFACT code; chapter 4 will give an overview of two models that were proposed to

explain the spectral hardening; the actual analyses of the hardening will be presented in chapter 5; chapter 6 will give an assessment of the number of unresolved points sources and whether these could mimic the hardening; finally, chapter 7 contains the conclusion of the results presented in this thesis and the outlook to future studies.

(14)

(15)

Chapter 2

Cosmic Rays

As already mentioned in the introduction of this thesis, cosmic rays cover a wide range of energies and the all-particle spectrum (see Figure 2.4) consist of several broken power laws. Over the last century physicists have been concerned with where these particles came from and how they have reached us. The following sections will cover the basics of particle acceleration at astrophysical sites and the diffusion of cosmic rays throughout our Galaxy.

2.1 Cosmic-Ray Acceleration

One of the first questions one might ask is: where do cosmic rays get their energy from? In other words: how are they accelerated?

The total power of cosmic rays in the Galaxy is approximately LCR ≈ 5 × 1040

erg/s [28]. The nature of this power was for a long time unclear. A connection with supernovae was already proposed by Baade and Zwicky in 1934 [29] but only years later a physical mechanism in the form of diffusive shock acceleration was established. A typical

supernova explosion injects ∼ 1051_{erg into the interstellar medium. With a supernova}

rate of a few per century [30] this gives a total power injected into the Galaxy of about

1041 _{ergs/s [28]. This means that if only a few percent of the power of a supernova}

explosions goes into the acceleration of particles, this would be enough to explain the energy that is present in cosmic rays.

2.1.1 Fermi Acceleration

Enrico Fermi was the first to think about acceleration of particles by magnetic clouds and, doing so, could provide a physical mechanism for the acceleration of cosmic rays [31].

The idea is as follows: suppose a particle increases its energy at every encounter with

a magnetic cloud such that dE_dn = αE, where α is the relative gain in energy at every

encounter and n is the number of encounters. If N (n) is the flux of particles with n collisions and dN

dn =−

N ¯

(16)

for

N (> E)_{∝ E}−1/α¯n, (2.1)

and we obtain a power-law spectrum.

Fermi inspired a detailed theory of acceleration of particles by shocks that was pro-posed in four seminal papers [32–35]. This so called diffusive shock acceleration (DSA) is more relevant, physically, and deals with particles being accelerated by a large shock front (see Figure 2.1). An actual overview of DSA is beyond the scope of this thesis (see for example [6]), but it is fairly easy to show that particles on average gain energy at every shock crossing [28].

The shock front has velocity _−u1 and the shocked gas (downstream) moves away

from the shock with velocity u2 so that the relative velocity in the laboratory frame is

V = u2− u1. Having a particle with initial energy E1 in the lab frame, the energy in the shock frame (denoted with a prime) is E₁0 = γE1(1− β cos θ1), where γ is the Lortenz factor, β = V /c and θ1is the entrance angle and I have assumed the particle is relativistic so that pc = E. The outgoing energy in the lab frame is similarly E2 = γE10(1+β cos θ2). Now the energy change can be calculated as

∆E E1 = 1− β cos θ1+ β cos θ2− β 2_{cos θ} 1cos θ2 1_{− β}2 − 1. (2.2)

And averaging over the angles gives ∆E

E1 ∼

4

3β. (2.3)

We see that there is an average increase in energy at every shock-crossing that scales with velocity to the first order. In traditional Fermi acceleration one can similarly show that

the average energy increase scales with β2. This is why, confusingly, DSA is sometimes

called first order Fermi acceleration, while the original Fermi acceleration is dubbed second order.

2.2 Cosmic-Ray Diffusion in the Galaxy

A second question one might ask is: how do cosmic rays reach us at Earth? The fact that cosmic rays are charged leads one to think that magnetic fields might be involved. It can be shown through the treatment of a particle moving in a background of magnetic waves that the particle undergoes a diffusive motion.

2.2.1 Quasi-Linear Theory

The diffusion of cosmic rays through the Galaxy is commonly treated by so called quasi-linear theory(QLT). This treats diffusion in the Galaxy as scattering of charged particles on magnetic plasma waves called Alfv´en waves. Below I will follow the derivation by [6].

(17)

Figure 2.1Particle acceleration by a shock front as in diffusive shock acceleration. The particle can increase its energy at every shock crossing. Figure taken from [6].

A particle moving in a regular magnetic field B0 = B0z has a simple solution of theˆ Lorentz equation dp_dt = q_cv_{× B, namely:}

vx(t) = v⊥cos (Ωt) (2.4)

vy(t) =−v⊥sin (Ωt) (2.5)

vz(t) = vk = constant, (2.6)

where Ω = qB0

γmc is the gyro-frequency of the particle.

In QLT a perturbation δB is added assumingδB

B0 << 1. The perturbation propagates with the Alfv´en velocity vA = √B_4πρ0 , where ρ is the density of the material. Because

vA v, c, the magnetic perturbations are approximately fixed and we can use the static

limit. Setting

δBx =|δB| sin (kz + ψ) (2.7)

δBy =|δB| cos (kz + ψ), (2.8)

where k is the wave number of the magnetic wave and ψ a phase offset, one can now show there is a so called pitch angle scattering of the particle (see Figure 2.2).

Decomposing p_{→ p}_⊥+ p_k the equation of motion of the particle becomes:

dp_k

dt =

q

(18)

B0zˆ

~p

θ

Figure 2.2 Motion of a particle through a magnetic field. Scattering with pitch angle θ.

Defining µ = cos θ, the angle between p and B (see Figure 2.2), then p_k = pµ(t) and

writing v_⊥= _γmp sin θ = _γmp _{p1 − µ}2 _{so that} pdµ dt = qp γmc p 1_{− µ}2_{[cos (Ωt)δB}_y _{− sin (Ωt)δB}_x_]. _(2.10)

Inserting Equation 2.8 using z = vµt, γ = _mcΩqB and using cos A cos B + sin A sin B =

cos (A_{− B)):} dµ dt = Ω δB B0 p 1_{− µ}2_{cos [(Ω}_{− vµk)t + ψ].} _(2.11)

Integrating to obtain the mean of the pitch angle_{h∆µi one can see that it averages}

to zero for long timescales. But when Ω = vµk the variance_h∆µ∆µiΨ does not vanish:

h∆µ∆µiΨ = πΩ2 δB

B0 2

1_{− µ}2

vµ δ(Ω− vµk)∆t. (2.12)

The power spectrum of magnetic wave modes is defined as P (k) = δB_B22/4π

0/8π, so that

P (k)dk is the energy density in dk waves. Defining kres = Ω/(µv) = 1/(µRL) one can

then write

h∆µ∆µi

∆t = π(1− µ

2_)Ωk

resP (kres). (2.13)

Changing variables back to θ:

h∆θ∆θ

∆t i =

π

4ΩkresP (kres). (2.14)

The average time τ for a particle to change direction by δθ _{∼ 1 is then the inverse of}

above so that λmfp = vτ . We can introduce a diffusion coefficient D = 1₃λmpfv (see

subsection 2.2.2) and D is usually written as

D = 1

3 RLv

(19)

where RL = v/Ω is the Larmor radius and F = kresP(k_B2res)

0 is the power of the magnetic

turbulence at kres.

The diffusion coefficient depends on momentum p and is usually parametrized as D(p) = _ρρ

0 δ

, where ρ = p_q is the rigidity of the particle and δ is a scale factor. The

scale factor δ depends on the power spectrum of the magnetic turbulence (F). Two

often assumed models for this power spectrum are a Kolmogorov [36] or Kraichnan [37] spectrum which predict δ = 1/3 and δ = 1/2 respectively.

The above translates to the following: when 1/kres RL the particle is surfing

along the perturbation and the direction is not changing. When 1/kres RL the

Larmor radius is too big for the magnetic perturbation and it has no effect. So only

when 1/kres∼ RL the particle feels the perturbation and the direction is changing (see

Figure 2.3). A proton with an energy of several GeV has a Larmor radius of about the

size of an astronomical unit (_{∼ 10}11_{m). Thus, the corresponding magnetic perturbations}

have to be of the same order.

k−1 RL

direction is changing k−1∼ RL

Figure 2.3 Particle trajectory (solid line) on the perturbed magnetic field (dashed line). If 1/k is bigger than the Larmor radius the particle is surfing along the perturbation and the direction is unaffected. Oppositely, the particle trajectory is too large for the perturbation and the particle does not feel the perturbation. Only when 1/k has the same size as the particles Larmor radius the particle trajectory is affected by the perturbation.

(20)

2.2.2 Diffusion Equation

Now that we see that cosmic rays undergo diffusion we can set up a diffusion equation. A general diffusion equation can be derived from flux conservation:

∂N

∂t =−∇ · J + Q, (2.16)

where N is the number of particles per unit volume, J is the particle flux and Q is an injection term. A diffusion coefficient D relates the flux to the number of particles

J=−D∇N, so that a general diffusion equation is written as:

∂N

∂t =∇ · (D · ∇N) + Q. (2.17)

D can be interpreted as

D = 1

3λmpfv,¯ (2.18)

where λmpf the average mean free path and ¯v the average velocity of the particle. The full diffusion transport equation for cosmic rays as written down by Ginzburg [15] reads:

∂N

∂t =∇ · (D · ∇N) −

∂

∂E[bi(E)Ni(E)]− ∇ · u Ni(E) + Qi(E, t) − piNi+ vρ m X ki Z _dσ i,k(E, E0) dE Nk(E 0_)dE0_{, (2.19)}

where i, k denote different particle species.

The terms on the right hand side correspond to the following: (1) diffusion, (2)

energy loss with b(E) = dE_dt, (3) convection with velocity u, (4) the source term, (5)

losses by decay and collisions with pi = vρσi/m + 1/(γτi) where σ is the cross-section and γτ the Lorentz dilated lifetime and finally (6) is a cascade term that incorporates gains and losses through fragmentation.

2.2.3 Leaky Box Model

In the leaky box model we assume that cosmic rays are confined to some fiducial volume

where they escape at the boundaries and that there is a mean timescale τdiff which

captures the mean time cosmic rays spent in the Galaxy [28]. For highly energetic particles we can replace the diffusion term by N/τdiff. Going to the steady state limit

∂N

∂t = 0 and neglecting all other terms:

D∇2N = Q⇒ N = τdiff(p)Q(p). (2.20)

(21)

Figure 2.4 Cosmic-ray spectrum as mea-sured on earth by several experiments. Fig-ure taken from [28].

Figure 2.5 Boron to carbon ratio as a function of energy per nucleon as measured by several cosmic-ray experiments. Figure taken from [38].

Secondary cosmic rays are particles that can only be created by spallation processes from primary cosmic rays. Measuring the fraction of the abundance of secondary versus primary cosmic rays one can then estimate the amount of matter transversed and make an estimate of the value of the diffusion coefficient. A common measurement is the Boron to Carbon abundance ratio, where Boron is the secondary and Carbon the primary cosmic-ray species.

Writing down the number density of Boron as in the leaky box model and realizing this should be equal to the decay of other species, dominated by Carbon:

NB= Qτdiff ∝ σC→B NC τdiff, (2.21)

with σC→B the cross section of decay from carbon to Boron. Then the fraction

NB

NC ∝ τdiff ∝ D.

(2.22)

From the measured B/C ratio (see Figure 2.5) one can see that D ∝ ρδ_{, where}

ρ = p/q is the rigidity of the particle. The inferred slope of the diffusion coefficient (which is the slope of the power spectrum of magnetic waves, see subsection 2.2.1) is

about δ = 0.3–0.6. When the injection spectra is a simple power-law Q ∼ ρ−γ_{, then}

from Equation 2.21 one can see that the measured spectrum of primary cosmic rays at earth N _{∼ ρ}−γ−δ _{∼ ρ}−2.7.

(22)

2.2.4 Cosmic-Ray Propagation Codes

The diffusion equation as in Equation 2.19 can be solved analytically in some simple cases and provide important physical insight (see for example [23]), but the full solutions can become complicated very easily. Numerical simulations are necessary to provide accurate predictions for the whole range of cosmic-ray physics and with today’s computers are probably easier and faster than finding an analytic solution.

GALPROP1 was the first pioneering numerical code that solved the cosmic-ray diffusion equation and could provide accurate predictions for measurements of secondary to pri-mary ratios and the diffuse gamma-ray emission from the Galaxy [39, 40]. Since GALPROP several independent codes like DRAGON [41] have been developed that can for example incorporate inhomogeneous diffusion [42].

These numerical propagation codes have a lot of parameters that can be tweaked. To name a few: the height of Galaxy, the normalization of the diffusion coefficient, the diffusion coefficients scaling with rigidity, the injection spectrum, the distribution of sources, etcetera. These parameters are usually fit to the data (direct cosmic-ray data and/or gamma-ray data), but there is still some wiggle room so that several models can be acquired (see for example [43]).

(23)

Chapter 3

Gamma Rays

Gamma rays cover the highest energies of the electromagnetic spectrum. They are produced by non-thermal processes where high energy physics is at play. In this chapter I will focus on the diffuse gamma-ray emission as observed by the Fermi Large Area Telescope and provide an overview of the physical processes responsible for the main components of this emission. I will briefly explain the conventional analysis of diffuse emission and describe the new tool SkyFACT that was used in this thesis.

3.1 Diffuse Gamma-Ray Emission

Since the discovery of the pion it has been speculated that there is intense gamma-ray emission from the Galaxy. This diffuse gamma-gamma-ray emission (DGE) results mainly from (1) π0 _{creation by cosmic rays interacting with interstellar gas, (2) inverse Compton} scattering (ICS) between photons from the interstellar radiation field (ISRF) and cosmic-ray electrons and (3) bremsstrahlung of charged particles moving through a charged medium. See Figure 3.1 for the spectra of these components. I will discuss below the

dominant physical mechanisms at Eγ≥ 1 GeV.

3.1.1 Pion Emission

When cosmic rays traverse the interstellar medium there is a chance that a π0 _{meson is}

created in a proton-proton collision:

p + p→ π0+ X, (3.1)

with X the remaining particles. The π0 subsequently decays into two γ-ray photons that

each have an energy of mπ0/2 = 67.5 MeV in the rest frame of the π0. The gamma-ray

spectrum is symmetric around this pion bump that dominates the spectrum below a GeV. At higher energies we can directly compare the slope of the gamma-ray spectrum with the slope of the cosmic-ray proton spectrum because the production cross section of neutral pions is essentially flat at energies above_{∼ 1 GeV [44, 45].}

(24)

101 ₁₀2 ₁₀3 ₁₀4 ₁₀5 ₁₀6 E [MeV] 10−7 10−6 10−5 10−4 10−3 10−2 E 2dN /dE [MeV cm − 2s − 1sr − 1] Total Pion decay ICS Bremsstrahlung

Figure 3.1Spectra of the main components of diffuse gamma-ray emission according to GALPROP webrun Example 01 (Conventional/2D 4 kpc tuned to agree with ACE).

The interstellar medium is mostly comprised of hydrogen gas. It exists in the form

of neutral atoms (HI), neutral H2 molecules and in ionized form (HII). The HI medium

is traced by the 21-cm line radiation and the molecular medium by 2.6mm 12CO line

emission. This gaseous HI and H2 medium is the main target for cosmic rays creating

π0 _emission.

The distribution of gamma rays resulting from π0 decay closely follows the target

gas so that a template of the gas allows to isolate the π0 component resulting from

cosmic-ray interactions.

The atomic column density (NHI) can be derived form the 21-cm line emission

as-suming a constant spin excitation temperature TS (see e.g. [46]):

NHI(v) =− log 1− TB TS− Tbg TSCi∆v, (3.2)

where Tbg ≈ 2.66 K is the microwave background brightness temperature and Ci =

1.83_{× 10}18 _cm−2_.

The largest uncertainty in the equation above is in the assumed constant value of TS. HI in the ISM is most likely present in different phases so that TS can vary from 40

K to a few thousand K [20]. A typical value that is often used is TS ≈ 125 K [20, 26],

but [21] find a best fit value of TS = 140 K.

(25)

velocity of the gas assuming circular motion: vLSR = R V (R)_R −_RV

sin (l) cos (b), (3.3)

where R = 8.5 kpc and V = 220 km/s (see also Appendix B in [20]). This allows

to compute NHI for different ranges of R by integrating Equation 3.2 over the desired

velocity range.

Molecular H2 gas cannot be detected directly. Instead, it can be inferred by the

2.6-mm line emission of12CO. CO is the second most abundant molecule in the ISM and it

traces H2 because it is its main collisional partner [47]. The H2 column density (NH2) is proportional to WCO, the integrated line intensity of CO, through the XCO = NH2/WCO coefficient. The radial distance of H2 can derived in the same way as for HI.

Because ionized H is sub-dominant the total H column density is approximately:

NH' NHI+ 2XCOWCO. (3.4)

It should be noted that there is no reason the value of XCOis constant across the Galaxy. This is one reason why it is important to know the radial distance of the gas. In this way one can partition the gas into several Galactocentric annuli to create multiple templates that can vary when fitting to the data. This is also important if one wants to study properties of cosmic rays throughout the Galaxy.

The above tracers of gas are not perfect. Especially CO is not a perfect tracer for

H2, leading to an underestimation of the total gas column density. This dark neutral

medium (DNM) can be traced, partly, by dust [48] or it can be inferred from gamma-ray data itself [49].

3.1.2 Inverse Compton Scattering

Inverse Compton scattering (ICS) is the up-scattering of ambient photons by high energy electrons. Instead of the electron gaining energy from an incoming photon with regular Compton scattering, ICS produces a high energy photon that has gained energy from the electron.

Unlike the π0 _{emission which traces gas observed at other wavelengths, there is no}

similar direct observation linked to ICS, so instead it must be calculated. This can be done by using cosmic-ray propagation codes that incorporate a distribution of electron sources and a diffusion model and folds the result with the interstellar radiation fields (ISRF) (see for example [39]).

The ISRF is the electromagnetic radiation that fills the Galaxy (see Figure 3.2). The main contribution is the emission from stars which peaks in the infrared, then there is also a contribution in the far-infrared from dust and a small contribution from the CMB. The ISRF also has to be modelled and relies heavily on the assumed distribution of sources and stellar luminosities. In combination with the uncertainties of cosmic-ray propagation codes, this makes that the ICS template should be treated with care.

(26)

m) µ ( λ -1 10 1 10 102 103 ) -1 m µ -3 m eV cm µ (_λ u λ -2 10 -1 10 1 10

Figure 3.2 ISRF as modeled by [50]. The different lines correspond to different distance from the Galactic center. From top to bottom respectively R = 0, 4, 8, 12 kpc.

3.2 Conventional Diffuse Gamma-Ray Analyses

Conventional analysis of gamma-ray data relies on fitting a linear combination of various templates that correspond to different emission phenomena. These templates can be based on models of cosmic-ray propagation or on observations at other wavelengths,

and can be energy dependent. For example: a template for π0 decay from cosmic rays

interacting with gas (see subsection 3.1.1) will depend on the distribution of gas and on the cosmic-ray distribution as calculated by propagation codes such as GALPROP or

DRAGON.

A typical diffuse gamma-ray model looks like this: Ntotal(l, b, E) = qH(E)NH(l, b) + qICS(E)NICS(l, b, E)+

X

i

NPS(E)δ(l− li, b− bi), (3.5)

where Ntotal, NH(l, b), NICS and NPS indicate the counts expected for the total predic-tion, the gas, ICS and point source components respectively at pixel with longitude l

and latitude b in energy bin E. The values of qH and qICS represent the normalizations

and are determined by the fit to the gamma-ray data in every energy bin. The final sum is over all point sources with l = li and b = bi.

Unfortunately, these templates (Ni in Equation 3.5) come with unknown systematic

(27)

uncertainties (see for example [20]). Most of these models however give a poor fit to the data with strong residuals specifically along the Galactic disk. This can lead to to biased conclusions in indirect dark matter searches or cosmic-ray propagation theory (see for example [51]).

One way to address this problem is to add a large number of nuisance parameters. These parameters allow for variations in the imperfect models and as such improve the fit to the data. This method of image reconstruction is also used in other disciplines such as in positron emission tomography (PET) [52].

3.3 SkyFACT

SkyFACT is a tool that fits Fermi-LAT gamma-ray data with spectral and spatial

tem-plates with a large number of nuisance parameters [26]. It is a hybrid approach between conventional template fitting (see above) and image reconstruction. This allows to ac-count for systematic uncertainties that exist in gas templates and cosmic-ray propagation codes.

The model fitted to the data at pixel p and in energy bin e can be written as a sum over model components k:

Nep = X

k

T_p(k)τ_p(k)_{· S}_e(k)σ_e(k)_{· ν}(k), (3.6)

where T and S are the spatial and spectral templates respectively, τ and σ denote the spatial and spectral nuisance parameters and ν is the overall normalization. For point sources the model is similar but there is no spatial component. When convolved with the Point Spread Function (PSF) the result of Equation 3.6 will give the expected number of counts µ per bin.

3.3.1 Likelihood Fitting

The Fermi-LAT measures number of counts making it natural to use a Poisson likelihood. The likelihood that is minimized, with c the observed and µ the expected number of counts, is the following:

−2 ln LP =−2

X

c_{− µ + c ln (µ/c),} (3.7)

where the sum is over all pixels and energy bins.

Because of the large number of nuisance parameters some sort of regularization is required to reduce shot noise and prevent overfitting. These regularization terms can be seen as priors in a Bayesian framework. They add to the Poisson likelihood so that

L = LP +LR with

−2 ln LR=

X

k

(28)

Here λ, λ0 and λ00 are three hyper-parameters controlling the regularization respectively on the spatial, spectral and overall modulation. Again for the point sources there is no spatial component so also no spatial modulation parameter.

RMEMis the regularization function based on the maximum entropy method (MEM),

it is

λRM EM(x) = 2λ

X

i

1_{− x}i+ xiln xi, (3.9)

where λ controls the strength of the regularization and xi are the model components.

The MEM minimizes the entropy, so that this gives the image that contains the least information while still being consistent with the data [53]. This regularizes the fitted parameters so that no spurious information is added. It is convex with a minimum at xi = 1 and the variance goes as 1_λ for large values of λ. A large value of λ thus constrains the fitted parameters to be close the initial input, while at small λ the minimization will give more weight to the data.

There is also the possibility to add a smoothing regularization in SkyFACT, which constrains the difference between adjacent pixels, effectively smoothing the image. This is controlled by two more hyper-parameters for every component, one for spatial and one for spectral smoothing (see [26]).

The negative of the sum of the Poisson likelihood and the regularization,_{−2 ln L , is} minimized using the L-BFGS-B algorithm [54] to find the best fit parameters for each model component at every pixel and energy bin, respectively, τ and σ and the best fit overall normalization ν.

3.3.2 Parameter Uncertainties

The covariance matrix contains information on the errors of the fitted parameters.

SkyFACTapproximates the covariance matrix by the inverse of the Fisher matrix, which

gives a lower limit on a set of unbiased estimators ˆθ from the Cram´er-Rao bound [55]:

Σij = cov(ˆθiθˆj)≥ I(ˆθ)−1ij . (3.10)

The Fisher matrix is defined as follows: Iij =− ∂2 ∂θi∂θj ln_L , (3.11)

where θi is the i-th model parameter. The average is taken over mock data generated

from the best-fit model. This expression can be calculated analytically (see for example [56]) and the resulting Fisher matrix is, in this case, sparse at the 99% level [26]. This makes the Fisher matrix computationally efficient to work with when dealing with a large number of parameters as is the case here.

Because the inverse of a sparse matrix is not sparse itself, SkyFACT does not calcu-late the covariance matrix directly, but instead samples the errors from a multivariate Gaussian with covariance matrix Σij =Iij−1.

(29)

The Fisher matrix can be decomposed by Cholesky factorization such thatI = LLT_. A vector of standard normal distributed variables with the length of the number of fitted parameters x _{∼ N (0, 1) is created, so that hx}ixji = δij. A multivariate normal

distributed vector y with covariance Σij can then be generated by y = (LT)−1x, which

can be solved for y by back-insertion. The vector y contains the correct errors because hyiyji = (LT)−1xixjL−1 = Σij, so that δˆθ = y contain the variances of the model parameters when averaging over many samples of x.

(30)

(31)

Chapter 4

Physical Interpretations of the

Progressive Hardening of the

Cosmic-Ray Spectrum

For a long time the treatment of cosmic-ray diffusion as explained in chapter 2 has been sufficient in explaining the available data, both from local cosmic-ray measurements and from gamma-ray observations. Today, with the high precision and high statistics of Fermi-LAT this standard treatment can be questioned.

The progressive hardening in the cosmic-ray slope as observed by [21, 22] cannot be explained with the standard treatment of cosmic-ray diffusion (as described in chapter 2), where diffusion is assumed to be constant throughout the Galaxy, and can be sign of a deviation from this standard treatment. Several models have been put forward to explain this anomaly. Below I will briefly discuss the interpretation as put forward by [23] that proposes an explanation in view of non-linear cosmic-ray diffusion and the model by [25] that incorporates anisotropic transport in the complex Galactic magnetic field as a solution to the spectral hardening in the inner Galaxy.

4.1 Non-Linear Effects in Cosmic-Ray Propagation

The relatively simple picture of cosmic rays scattering on a perturbed magnetic field as described in subsection 2.2.1 is likely to be much more complex. There are several reasons that justify non-linear effects in the propagation of cosmic rays (see for exam-ple [6]). One of the reasons is that cosmic rays can be responsible themselves for the turbulent magnetic field on which they scatter, so that the diffusion coefficient depends on the distribution of particles, which in turn depends on the diffusion coefficient, i.e. a non-linear effect is present. Cosmic rays can thus scatter on magnetic waves that are generated by themselves [57].

In [23] the authors use a model where cosmic rays scatter on self-generated waves and show that this can explain harder cosmic-ray spectra near the sources, which are

(32)

more dominant towards the Galactic centre.

They incorporate a simple 1D diffusion-convection equation in phase space (similar to Equation 2.19, but with only diffusion and advection) and assuming only transport along the z-direction:

−∂ ∂z D∂f ∂z + u∂f ∂z − p 3 du dz ∂f ∂p = Q(p)δ(z). (4.1)

In Equation 4.1 f (z, p) is the particle distribution, normalised so that R 4πp2_{f (z, p)dp}

is the particle number density. Furthermore, Q(p) ∝ p−γ _{is the injection spectrum,}

assumed to be a single power-law with index γ, and u = B0/√4πρ is the advection

velocity assumed to be the Alfv´en speed.

They assume D is independent of z and u is constant except for being negative for

negative z and positive for positive z so that du

dz = 2|u|δ(z). One can then integrate

between z = 0₋ and z = 0+ to get

−2D∂f_∂z −2u₃ pdf0

dp = Q(p), (4.2)

where f0(p) = f (0, p).

Solving Equation 4.1 far away from z = 0 at z = H so that ∂f_∂z z=H = 0: D∂f ∂z = uf (z, p), (4.3) is solved to give: f (z, p) = f0(p) 1_{− e}−u(H−z)/D 1_{− e}−uH/D . (4.4)

Using the derivative of Equation 4.4 and putting it into Equation 4.2 gives f0(p) = 3 2u Z ∞ p dp0 p0 Q(p) exp " Z p00 p 3 1_{− e}uH/D # . (4.5)

From subsection 2.2.1 we saw that D = 1₃ RLv

F(k). The power spectrum of the turbulent

magnetic waves, _{F(k), can be approximated by equating the growth rate of magnetic}

waves induced by cosmic rays due to cosmic-ray streaming instability [58]

ΓCR= 16π2 3 u F(k)B2 0 p4v(p)∂f ∂z p=eB0/kc (4.6) and the damping of those waves that occurs at a rate [59]

(33)

with ck= 3.6 and k = 1/RL. This gives F(k) = (2ck)3 16π2p4 B₀2 D ∂f ∂z 2 . (4.8)

With Equation 4.4 we can then solve for the diffusion coefficient

D(z, p) = RLv(p) 3(2ck)3 " B2 0 16π2_p4 1_{− e}−uH/D uf0(p) #2 + 2u(H_{− z).} (4.9)

In the diffusion dominant regime when the first term in Equation 4.9 is dominant over the second term we recover the leaky box model and

f0(p)∝ Q(p) B0

3

∝ p7−3γ. (4.10)

In the opposite limit when advection is dominant the cosmic-ray spectrum is harder f0(p)∝

Q0

B0 ∝ p

1−γ_. _(4.11)

Because Q0 and B0 depend on Galactocentric radius, there is a gradient in the spectral

index of cosmic rays (see Figure 4.1). Intuitively: near the sources the diffusion coefficient is smaller, so that there is more confinement of cosmic rays, resulting in advection to be dominant up to higher energies. Because advection is not dependent on momentum the cosmic-ray spectrum is closer to the injection spectrum, which is harder.

Note that this effect is clearly energy dependent and works only when diffusion is

comparable to advection, this holds roughly for cosmic-ray protons with energy ∼ 10

GeV.

In [23] the authors find a model that describes the gamma-ray data fairly well for cosmic-ray protons of 20 GeV (roughly corresponding to 2 GeV gamma rays). See Figure 4.1.

4.2 Anisotropic Transport

The simplistic picture painted by QLT that diffusion is uniform in the Galaxy is chal-lenged by measurements such as the harder cosmic-ray spectrum in the inner Galaxy. QLT predicts that the motion of cosmic rays is dominated by diffusion parallel to the regular magnetic field:

D_⊥

D_k 1, (4.12)

but this does not necessarily hold in environments with large turbulence, such as injected by supernova remnants [18].

(34)

0 5 10 15 20 25 30 35 R [kpc] 2.2 2.4 2.6 2.8 3.0 Sp ectral index Recchia et al. (2016) Cerri et al. (2017) Acero et al. (2016) Yang et al. (2016)

Figure 4.1 The spectral hardening as observed by [21, 22] with predictions by [23, 24]. The most inner point from [21] is the CMZ. Also note that [22] give the photon index above 2 GeV and [21] give the proton index P2 from fitting a proton spectrum to the gamma-ray data of the

form βP1_ρP2_{, where β = v/c and ρ is the rigidity. The prediction by Cerri et al. (2017) [25] is}

with δ⊥ = 0.7 and D= 0.01 (see text).

The authors of [25] propose a model were the diffusion coefficient is made up of a parallel and a perpendicular part. Simulations of charged particles in magnetic fields

show that the rigidity dependence of D_⊥ can be significantly steeper than that of D_k

[60]. In combination with the recently revealed poloidal component of the Galactic magnetic field; that is mostly dominant in the inner Galaxy, while further away the toroidal component takes over, this could lead to anisotropic diffusion of cosmic rays.

In [25] the authors propose a toy model that illustrates the different diffusion mech-anisms at play in the Galaxy which I will briefly discuss below.

Implementing the magnetic field in the Galaxy as:

BR= 0, Bφ= B0,φ(1− e(−R/R0)), Bz = B0,ze(−R/R0)= BB0,φe(−R/R0), (4.13)

where R0 is the distance where the poloidal component (Bz) takes over and B =

B0,z/B0,φ.

Decomposing the diffusion coefficient as: Dij = D⊥δij+ (D_k− D⊥)BiBj |B|2 , (4.14) with D_k= D₀_kρδk_{, D} ⊥= D0⊥ρδ⊥ = DD0kρδ⊥, (4.15)

(35)

where D = D0⊥/D0k, ρ is the rigidity of the particles and δ⊥ and δk are different scale factors for perpendicular and parallel diffusion respectively.

Considering a simple diffusion equation of the form: ∂N ∂t = ∂ ∂xi Dij ∂N ∂xj + Q, (4.16)

where Q_{∝ ρ}γ _{is the injection spectrum.}

With this setup one can show that the diffusion at R R0 reads:

Dij = Drr Drz Dzr Dzz ! ≈ D⊥ 0 0 D_k ! ≈ Dk Dρδ⊥−δk 0 0 1 ! , (4.17)

where I used Equation 4.13 and Equation 4.14 to calculate Dij. We see that the diffusion

at small R is dominated by the parallel component. For R R0 the diffusion coefficient

is dominated by its perpendicular component Dij ≈

D_⊥ 0

0 D_⊥

!

. (4.18)

Solving for the final particle spectrum N (ρ) one gets:

N (ρ)_∝

(

ρ−(γ+δk) _{for R}_R₀

ρ−(γ+δ⊥) _{for R}_R₀_. (4.19)

So that if the scaling of the diffusion with rigidity is harder for parallel escape this would give a harder spectrum going towards the inner Galaxy.

The authors of [25] implement this different scaling of the perpendicular and parallel

escape in their DRAGON cosmic-ray propagation code and test for different values of D,

(36)

(37)

Chapter 5

Progressive Hardening of the

Cosmic-Ray Spectrum with

SkyFACT

The hadronic component of the cosmic-ray spectrum throughout the Galaxy can be

studied by analysing the diffuse gamma-ray emission that is attributed to π0-decay

(see subsection 3.1.1). Away from the pion bump the gamma-ray spectrum follows the cosmic-ray spectrum, so that the slope of the gamma-ray spectrum at E > 2 GeV is comparable to the slope of the cosmic-ray spectrum.

The hardening as measured by [21, 22, 24] is not a direct measurement and depends on the precise modelling of the gamma-ray sky. SkyFACT can be an excellent tool to robustly characterize this feature because it incorporates a large number of nuisance parameters that can account for imperfect modelling.

A hint of a hardening can already be seen in the results from [26]. In Figure 5.1 the spectrum for Gas ring II looks less steep than for Gas ring III. Furthermore, in Figure 5.2 it can be seen that the rescaling of the spectrum in Gas ring II is harder than in Gas ring III.

In this chapter I will explain how the gamma-ray sky was modelled with SkyFACT and present and discuss the results of the spectral hardening found.

5.1 Method

5.1.1 Data

For this study 9.3 years, from 4 August 2008 to 6 November 2017, of Fermi-LAT Pass 8

data were used. Standard LAT data analysis tools as available from the NASA website1

were used to create counts and exposure maps. Only ULTRACLEAN events (evclass=512)

and FRONT and BACK (evtype=3) events were selected. A zenith angle cut at 90◦was used

(38)

100 ₁₀1 ₁₀2 E [GeV] 10−9 10−8 10−7 10−6 10−5 E 2dN /dE [GeV cm − 2s − 1sr − 1] Full ROI Gas ring I Gas ring II Gas ring III Total Data

Figure 5.1Spectra of the π0_{emission from}

the three gas rings from RUN5 in [26].

10−1 10−1 100 100 101 101

log10(Modulation parameter) Gas ring I

Gas ring II Gas ring III

Figure 5.2 Rescaling of the spectra from RUN5 in [26]. See Figure 5.7 for a detailed description.

to reduce photons from the Earth’s limb and the data quality cuts recommended by the Fermi-LAT collaboration (DATA_QUAL>0) && (LAT_CONFIG==1) were applied. Counts

and exposure maps of 720_{× 81 pixels with a pixel size of 0.5 degrees centered at the}

Galactic center were made. The photon energy was binned in 25 log-spaced energy bins between 0.34–228.65 GeV. Our region of interest spans the whole disk (_{|l| < 180}◦) with latitudes _{|b| < 20.25}◦. The data counts map integrated over all energies is shown in Figure 5.3.f

5.1.2 Model Components

The π0 _{template that was used in the analysis with SkyFACT is based on the gas column}

densities available with the GALPROP public release. The templates are made as described

in subsection 3.1.1. The HI and H2 maps are based on the analytical models in [61],

normalised to, respectively, the Leiden-Argentine-Bonn 21 cm line survey [62] and the

2.6 mm CO line survey [63]. The spin temperature for the HI maps is fixed to TS= 125 K

and XCO = 1.9× 1020 cm−2/(K km/s) (see subsection 3.1.1). Note that no dark gas

correction was made, because it was shown that SkyFACT can recover this through the use of nuisance parameters [26]. Also, ionized gas (HII) (responsible for bremsstrahlung)

180 90 0 −90 −180 l [deg] −20 0 20 b [deg] Data (0.34− 228.65 GeV) 2.0 2.5 3.0 3.5 4.0 log 10 (coun ts)

(39)

was neglected, because it is sub-dominant compared to π0-decay at higher energies (see Figure 3.1).

The resulting gas maps are binned in Galactocentric radii according to their velocity (see subsection 3.1.1). The binning initially used is shown in Table 5.1, this is the same as used by the Fermi collaboration in [21]. See Figure 5.4 for the best-fit templates of the gas in the different radial bins.

As a starting point the measured π0 spectrum from [20] was used. This spectrum is

in agreement with local cosmic-ray observables and is approximately a power-law with an index of 2.7 at energies > 2 GeV.

The ICS model was created with the DRAGON [41] (see subsection 2.2.4) and GammaSky [18, 43] codes. The model is based on the ISRF as documented in [64] with the ‘Ferri`ere’ source model from [65] and the ‘KRA4’ propagation model from [43]. The initial spec-trum was again taken from [20].

The model components and regularization parameters (see section 3.3) used in this thesis are very similar to the settings of the model labeled RUN4 in [26] (see Table 1 therein). The difference between RUN4 and their best fit model RUN5 is that RUN5 includes a template for the GeV excess in the Galactic center. This was not used for this analysis because it did not change the final result. Other differences from RUN4 used here are the larger region of interest; the different binning in the gas maps as described above and the Fermi bubbles are fitted with a flat template, spatially following the region as defined in [66], effectively masking the area. The spectrum of the Fermi bubbles that is used is from [67] fixed to within 1% (λ0= 10000)2. Furthermore, the necessary point sources and extended sources within our ROI were added from the 3FGL catalogue [68]. See Table 5.3 for an overview of the regularization used for each component for each of the different fits to the gamma-ray sky reported in this thesis.

Ring Rmin(kpc) Rmax (kpc)

1 0.0 1.7 2 1.7 4.5 3 4.5 5.5 4 5.5 6.5 5 6.5 7.0 6 7.0 8.0 7 8.0 10.0 8 10.0 16.5 9 16.5 50.0

Table 5.1Binning of gas maps in Galacto-centric radius.

Ring Rmin (kpc) Rmax (kpc)

1 0.0 1.7 2 1.7 4.5 3 4.5 6.5 4 6.5 8.0 5 8.0 10.0 6 10.0 16.5 7 16.5 50.0

Table 5.2Alternative binning of gas maps in Galactocentric radius.

2_{Reminder: λ, λ}0 _{and λ}00 _{denote respectively the spatial, spectral and overall hyper-parameters}

(40)

Components RUN4A RUN4B RUN4Ca RUN4Cb h λ λ0λ00 η η0 · i IGRB [∞ 16 ∞ 0 0 · ] [∞ 16 ∞0 0 · ] [∞ 16 ∞0 0 · ] [∞ 16 ∞0 0 · ] 3FGL PSC [· 25 0 · 0 ·] [· 25 0· 0 ·] [· 25 0· 0 ·] [· 25 0· 0 ·] Extended Sources [0 1∞ 4 0 · ] [0 14 0∞· ] [4 00 1∞· ] [0 14 0∞· ] Fermi bubbles [∞ 10000 ∞ 0 0 · ] [∞ 10000 ∞0 0 · ] [∞ 10000 ∞0 0 · ] [∞ 10000 ∞0 0 · ] ICS [ 1 16 0 100 0 ·] [100 01 16 0·] 0–3 kpc [ 1 1 0 100 0·] [100 01 16 0·] 3–8.3 kpc [ 1 1 0 100 0·] [1001 10000 00 ·] 8.3–50 kpc [ 1 1 0 100 0·] [1001 10000 00 ·] Gas rings 9_{× [}10 16 0 25 0 ·] 7× [10 16 025 0 ·] 9× [10 16 025 0 ·] 9× [10 16 025 0 ·]

Table 5.3Overview of the different model components used in this analysis. The matrix denotes the different regularization hyper-parameters, where λ, λ0 and λ00 are the spatial, spectral and overall modulation parameters respectively and η and η0 the spatial and spectral smoothing parameters respectively. For the binning of the gas rings see Table 5.1 and Table 5.2.

5.1.3 Neutral Pion Spectra

The spectra in each of the gas rings (i.e. the π0 components) are allowed to vary so

that we can reproduce different spectra across the Galaxy. The spectral regularization

parameter is chosen so that the spectra are constrained to within 25% (λ0 = 16) of the

initial spectrum. This allowed the spectra to vary enough to produce a similar spectral hardening as in [21, 22]. Less regularization gave unphysical results, especially in the inner-and outermost radial bins and did not change the spectral index significantly in the other bins.

Away from the pion bump the π0 spectrum follows the cosmic-ray spectrum and

Eγ≈ ECR/10. The resulting π0 spectra for every radial bin (after fitting with SkyFACT)

are fitted with a single power-law in the energy range 2− 228.65 GeV to obtain the

spectral index of the cosmic-ray spectrum. A broken power law of the form dN dE ∝ ( E−a for E < Ebr E−bE_br(b−a) for E≥ Ebr (5.1) with Ebr = 30 GeV was also fitted to assess if the hardening diminishes at higher energies. To get more sensible error bars the spectral index was fitted by sampling from the posterior distribution using the Markov chain Monte Carlo (MCMC) python package

emcee [69]. The flux in each energy bin is assumed to follow a Gaussian distribution

with standard deviation on the differential flux as calculated from the Fisher matrix in

SkyFACT (see subsection 3.3.2). The standard deviation on the flux is affected mostly

by the spectral regularization hyper-parameter λ0 (see section 3.3): larger constraints

(41)

A flat prior from 0 to∞ was used on the normalization to ensure only positive values and a flat prior on the power-law index (or indices) from 0 to 10 so to not get unphysical results. The posterior distribution was sampled in 1000 steps, after removing the first 50, with 500 walkers. The posterior distribution was checked for inconsistencies and the mean value of the slope was chosen as the best fit power-law index. The errors on the spectral index reported are 68% credible intervals.