Hunting for PeV proton accelerators
in Galactic supernova remnants by
γ-ray observations

Tomohiko Oka

Department of Physics, Graduate School of Science, Kyoto University
Kitashirakawa Oiwake-cho, Sakyo, Kyoto, Kyoto, 606-8502, Japanoka.tomohiko.25n@st.kyoto-u.ac.jp

January 2023

Abstract
Although cosmic rays (CRs) are an important actor in the universe, their origin is not wellunderstood. In particular, CRs up to PeV energies are believed to be produced by supernova
remnants (SNRs) in our Galaxy, but there is no conclusive observational evidence for this. In
this work, we analyze γ-ray observation data for the two SNRs, G106.3+2.7 and HB9, to reveal
proton acceleration up to PeV in the SNRs.
SNR G106.3+2.7 is a promising candidate as the PeV protons accelerator since previous studies revealed γ rays above 100 TeV around this SNR. However, this contains a pulsar wind nebula(PWN) producing high-energy electrons, and the γ-ray origin is thus unknown. To identify theorigin by spatially resolving the emission region, we performed the TeV γ-ray observations with
better angular resolution (0.07–0.1◦ ) than previous studies in the band by using the MAGIC
telescope. We found significant emissions above 0.2 TeV spatially coinciding with the SNR radio
shell, while emissions above 6 TeV are detected only far from the PWN. This fact suggests thatthe 100-TeV γ ray originates in the SNR rather than the PWN. Furthermore, spectral modeling
results have shown that leptonic models cannot explain the energy spectrum, while a hadronic
scenario works well and suggests proton acceleration up to ∼ PeV in this region.
Next, we tried to measure the time evolution of particle acceleration in SNRs with an aim toverify why γ-ray observations have found not many PeV proton accelerators. For this purpose,
we proposed a new measurement method and applied it to the γ-ray data around SNR HB9. In
this method, we expect that runaway CRs from the SNR illuminate the cloud and generate γ raysreflecting the particle distribution in the SNR at a specific epoch in the past, while the spectrumat the shell reflects the present distribution of particles. We can therefore quantify the evolutionof the acceleration in the SNR by comparing the spectra at the SNR shell and nearby clouds
observed simultaneously. Our Fermi-LAT analysis has revealed new γ-ray emissions spatiallycorrelated with the molecular cloud as well as the SNR shell. We then succeeded in explaining
these spectra with the scenario above, giving evidence for the CR acceleration by the SNR up
to higher energies (≳ 10 TeV) in the past than that produced in the present shell (300 GeV).
We also found that this kind of method allows for estimating the diffusion coefficient around
SNRs with few uncertainties.
The model developed in this method is applied to the observations for SNRs G106.3+2.7
and G335.2+0.1. Using the obtained model parameters, we confirmed that the luminosity of
CRs from SNRs is roughly consistent with that of the measurement at the Earth. Finally, we
discuss the future prospect with the next-generation γ-ray observatories.

iContents
Abstracti

1 Introduction

1

1.1

1.2

1.3

1.4

1.5

Candidate for CR origin: Supernova Remnants . . . . . . . . . . . . . . . . . . .

1

1.1.1

The SNR paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.1.2

SNR as a CR accelerator . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

CR acceleration and diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.2.1

CR acceleration in an SNR . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.2.2

CR propagation in the Galaxy . . . . . . . . . . . . . . . . . . . . . . . .

11

Lights from CRs in SNRs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

1.3.1

Neutral-pion decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

1.3.2

Synchrotoron . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

1.3.3

Inverse Compton scattering . . . . . . . . . . . . . . . . . . . . . . . . . .

17

1.3.4

Bremssthlung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

1.3.5

Contributions to γ-ray flux . . . . . . . . . . . . . . . . . . . . . . . . . .

17

Overview of the SNR observations and PeVatron search . . . . . . . . . . . . . .

20

1.4.1

Observations for SNRs . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

1.4.2

Sub-PeV γ-ray observations . . . . . . . . . . . . . . . . . . . . . . . . . .

23

Research objectives and thesis structure . . . . . . . . . . . . . . . . . . . . . . .

23

2 The MAGIC γ-ray telescope
2.1

2.2

26

Imaging Atmospheric Cherenkov Telescope (IACT) . . . . . . . . . . . . . . . . .

26

2.1.1

Air shower phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

2.1.2

Cherenkov emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

The MAGIC system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

2.2.1

Optics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

2.2.2

Photo-sensor and the readout system . . . . . . . . . . . . . . . . . . . . .

32

2.2.3

Atmospheric monitoring . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

ii

CONTENTS
2.2.4
2.3

2.4

2.5

iii
Drive system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

Event reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

2.3.1

γ-ray energy and flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

2.3.2

Arrival direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

2.3.3

Background estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

Systematic uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

2.4.1

Flux normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

2.4.2

Energy scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

2.4.3

Spectral slope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

Comparison in the performance of IACTs . . . . . . . . . . . . . . . . . . . . . .

47

3 MAGIC observation in the vicinity of SNR G106.3+2.7

49

3.1

PeVatron candidate: SNR G106.3+2.7 . . . . . . . . . . . . . . . . . . . . . . . .

49

3.2

Observation and data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

3.3

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

3.3.1

Morphological study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

3.3.2

Energy spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

Systematic uncertainties in the extended source analysis . . . . . . . . . . . . . .

56

3.4.1

Estimation of the effective area . . . . . . . . . . . . . . . . . . . . . . . .

56

3.4.2

Effect due to θ2 cut . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

4 Spectral modeling and implication to the γ-ray origin of SNR G106.3+2.7

66

3.4

4.1

4.2

Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

4.1.1

Bremsstrahlung-dominated model . . . . . . . . . . . . . . . . . . . . . .

67

4.1.2

IC-dominated model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

4.1.3

Hadronic-dominated model . . . . . . . . . . . . . . . . . . . . . . . . . .

68

Interpretation of the γ-ray origin . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

4.2.1

head region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

4.2.2

tail region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

4.2.3

Remarks on the modeling . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

5 γ-ray observations of SNR HB9 toward the measurement of the DSA history 74
5.1

5.2

A measurement method of the maximum energy evolution in an SNR . . . . . .

74

5.1.1

Previous observations of delayed γ-ray emission . . . . . . . . . . . . . . .

76

5.1.2

Introduction of SNR HB9 . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

Fermi-LAT observations in the vicinity of HB9 . . . . . . . . . . . . . . . . . . .

79

CONTENTS

5.3

5.4

iv

5.2.1

Fermi -LAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

5.2.2

Data selection and analysis . . . . . . . . . . . . . . . . . . . . . . . . . .

82

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

5.3.1

Morphology study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

5.3.2

Spectral results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86

Spectral modelling and its implication . . . . . . . . . . . . . . . . . . . . . . . .

86

5.4.1

Modeling of delayed γ-ray emission . . . . . . . . . . . . . . . . . . . . . .

88

5.4.2

Application to the cloud regions around HB9 . . . . . . . . . . . . . . . .

93

5.4.3

Accuracy of parameter determination . . . . . . . . . . . . . . . . . . . .

94

6 Discussion

98

6.1

Application of the delayed γ-ray model to SNR G106.3+2.7 . . . . . . . . . . . .

98

6.2

Delayed γ-ray candidate: HESS J1626−490 . . . . . . . . . . . . . . . . . . . . . 100

6.3

Discussion on the CR luminosity . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

6.4

Future prospect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.4.1

Future γ-ray observations . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

6.4.2

Neutrino observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

7 Summary110

Acknowledgements

112

Appendices113

A Gas density in the vicinity of SNR G106.3+2.7

114

B Electrons in the delayed γ-ray emission model116

B.1 Derivation of the spectrum of CR electrons . . . . . . . . . . . . . . . . . . . . . 116
B.2 Impact of electron parameter to the delayed γ-ray model . . . . . . . . . . . . . . 117
Bibliography

119Chapter 1

Introduction
Many mysteries remain in the evolution of the universe and the formation of stars and galaxies.
The origin of Cosmic Rays (CRs), which are the key to these mysteries (e.g., Padovani et al.,2020; Ohira & Murase, 2019), is, however, still unclear. CRs are charged particles (∼90% of
these are protons) that travel through space, and were first discovered by Hess (1912). More thana century of observational research since then has revealed the following interesting features:i. The highest energy of CR reaches > 1020 eV (Bird et al., 1995), which is many orders ofmagnitude higher than the energy that can be produced on Earth (13.6 × 1012 eV; Large
Hadron Collider in CERN1 ). Once we can identify an environment that produces suchhigh energy, it could provide a venue for testing the standard theory (e.g., Mattingly,
2005; Aartsen et al., 2021a; Kephart & Weiler, 1996).
ii. With sufficient interaction with its surroundings, particle energies should follow a Maxwellian(thermal) distribution. However, this is not the case for CRs, which follow a power-law
distribution (i.e., a non-thermal distribution), as shown in Fig. 1.1. This indicates that
the CR “explosively” acquire energy and reach the Earth without sufficient interaction.
iii. The energy density of CRs is ∼ 1 eV cm−3 , which agrees well with that of other components
for the universe, as listed in Table 1.1. This confirms that, contrary to (ii), there are someinteractions between CRs and the other components.
Still, the full picture of the CR origin - where and how they originate and are accelerated- remains unclear. An elucidation of their origin is thus one of the most essential issues inastrophysics.

1.1

Candidate for CR origin: Supernova Remnants

As shown in Fig. 1.1, the CR spectrum has multiple breaks, suggesting the existence of multiple
origins. Interestingly, the shape resembles a human foot, and thus certain parts are called by
the name of the foot (e.g., “Knee” at ∼ PeV, “Ankle” at ∼ EeV). As for candidates of the CR1

https://home.cern/science/accelerators/large-hadron-collider

1CHAPTER 1. INTRODUCTION

Intensity × E2 [GeV m 2s 1sr 1]

103

2

E 2.7

101
10 1
10 3
10 5
108

Knee
(~ PeV)

E 3.1

Auger
TA
KASCADE-Grande
IceTopTibet III
Ankle
HAWC
(~ EeV)
CREAM I+III (H+He)
AMS02 (H+He)
BESS-TeV (H+He)
AMS02 (e++e )
1014
1010
1012
1016
1018
Kinetic energy [eV]1020

Figure 1.1: The CR spectrum observed in the Earth. The data for the all-particle flux are
taken from Verzi (2020), Ivanov (2020), Apel et al. (2013), Aartsen et al. (2013), Amenomori
et al. (2008), and Alfaro et al. (2017). In the lower band, the total flux of proton and helium (Yoon et al., 2017; Aguilar et al., 2015b,a; Shikaze et al., 2007) as well as electron andpositron flux (Aguilar et al., 2018) are also shown.
Table 1.1: Energy densities in the Local interstellar medium (Draine, 2011).
Component
Energy density [eV cm−3 ]Cosmic microwave background (TCMB = 2.725 K)
0.26
Far-infrared radiation from dust
0.31
Starlight (hν < 13.6 eV)
0.54
Thermal kinetic energy (3/2)nkT
0.49
Turbulent kinetic energy (1/2)ρv 2
0.22
2
Magnetic energy B /8π
0.89
Cosmic rays
1.39

CHAPTER 1. INTRODUCTION

3

Figure 1.2: The predicted maximum energies as a function of the size and magnetic-field strengthof systems (Hillas, 1984). This is called Hillas plot after the first person who created it.
origin, the bottom-up scenario (i.e., accelerations in multiple celestial species) has been expected,
while the super heavy dark matter is also considered for the CRs at ∼ 1020 eV (Chung et al.,1998). In the bottom-up scenario, the Larmor radius formula can roughly predict the maximum
acceleration energy for each object in terms of how high energy the system can confine. The
maximum energy is hence evaluated as a function of the magnetic-field strength (B) and thesize (RL ) of the system, expressed as:
ZeRL B= 1 PeV Z
pmax ≲
c



B
1 µG



RL
1 pc


(1.1)

where Z is the atomic number, e is the elementary charge, and c is the speed of light. Fig. 1.2
shows a summary of the maximum acceleration energies for each object species based on Eq. 1.1.

CHAPTER 1. INTRODUCTION

1.1.1

The SNR paradigm

Among them, supernova remnants (SNRs) are believed to be a promising candidate for CR
origin, especially up to PeV energies from the point of the Hillas plot. The following discussions regarding the explosion energy budget and elemental composition have also justified thispossibility.
Energy budget in SNRs
Baade & Zwicky (1934) first proposed the possibility of SNRs as the CR source by considering theenergy budget. The explosion energy of SN is typically ESN ∼ 1051 erg (e.g., Utrobin & Chugai,2017), while the energy supply required to sustain CRs in the Galaxy is LCR ∼ 1041 erg s−1 . Ifassuming that SNs occur two or three times per 100 yr (∼0.03 yr−1 ; Tammann et al., 1994) and
10% of the energy can be used for the CR acceleration, the energy supply by SNR is found tobe LSN ∼ 0.1ESN /0.03 yr ∼ 1041 erg s−1 ∼ LCR , which is compatible to the CR energy.
Elemental composition in SNRs
Hayakawa (1956) showed that the elemental abundances of CRs (especially from H to Ni) produced by SN explosions roughly reproduce observations assuming that CR origins are uniformly
distributed in the Galaxy and travel through the interstellar medium (ISM) with a uniform density. By optimizing the contributions from SN Type I, Type II, and ISM, the observed elemental
abundances were successfully reproduced within a deviation by a factor of 2.5 (Yanagita et al.,1990). Furthermore, the deviation can be reduced to a factor of 1.6 by additionally consideringa contribution of Wolf-Rayet stars (Yanagita & Nomoto, 1999).
For these perspectives, the SNR is the leading candidate for the CR origin. Fig. 1.3 shows anexample of the modelling for the CR spectrum assuming the multiple origins (Thoudam et al.,2016), which is based on Eq. 1.1. The modeling assumes that dominant CR origins are regular
SNRs for ≲1 PeV, Wolf–Rayet SNe/SNRs for 1 PeV–1 EeV, and active galactic nuclei (AGNs)
for ≳1 EeV. It should be noted that the model can also reproduce the spectra of heavier nuclei,
as well as protons.

1.1.2

SNR as a CR accelerator

Since SNRs change with age during their lifetimes, the period that an SNR can accelerate CRs
is limited. SNRs are born when the blast of an explosion at the end of a star’s life (supernova)
sweeps up the ISM. The blast then compresses the ISM, creating a shock wave. The shock waves
are expected to accelerate CRs, as we will describe the detail in the next section. The evolutionof an SNR is separated into the following periods (e.g., Draine, 2011, and references therein).

4

H
He
C
O
Ne
Mg
Si
Fe

3

10

102

5

SNR−CRs
WR−CRsEG−RSB93
Total

C/He=0.1

101

3

E × Intensity [cm−2 sr−1s−1GeV2 ]CHAPTER 1. INTRODUCTION

10−1
1

10

102

3

10

104

5

6

10
10
107
Energy E (GeV)

8

10

9

10

10

10

1011

Figure 1.3: An example of the modelling result to the CR spectrum, taken from Thoudam et al.
(2016). The thick lines represent the model emissions from SNRs (blue solid), Wolf-Rayet stars(blue dashed), extra-galactic sources (blue dash-dotted), and total (red solid).

CHAPTER 1. INTRODUCTION

6(i) Free expansion phase
When the SNR age (tage ) is younger than ∼1 kyr, the mass of the swept-up ISM (MISM ) issmaller than that of the ejecta (Mej ), and thus the ejecta freely expands without a deceleration.
The shock speed can be expressed as:
s
vsh =

2ESN
∼ 10000 km s−1
Mej



ESN
1051 erg

1/2 

Mej
M⊙

−1/2
,

(1.2)

where ESN is the explosion energy.
(ii) Adiabatic expansion (Sedov) phase
Once MISM surpassed Mej , the shock decelerates due to the ISM. Assuming a point-like explosion
in a uniform ISM and then using the self-similar solution (Sedov-Taylor solution), the shock
velocity during this period is expressed as:
vsh = 490 km s

−1



tage
10 kyr−3/5 ESN
1051 erg

1/5 

nISM −1/5,
1 cm−3(1.3)

where nISM is the density of the ISM.
(iii) Radiative cooling (Snowplow) phase
When tage is older than 10 kyr, energy loss due to radiative cooling becomes significant. The
shock velocity gradually slows with the relation:vsh ∝ t−5/7
age ,

(1.4)

and finally, down to the speed of sound, by tage reaches 100 kyr. The ejecta merges with the
ISM, and then the shock compression declines. As a result, CR acceleration is no longer possible.

1.2

CR acceleration and diffusion

In this section, we describe the CR acceleration by shock waves (Sect. 1.2.1) and the diffusionof CRs (Sect. 1.2.2) toward an interpretation of the power-law form of the CR spectrum.

1.2.1

CR acceleration in an SNR

Diffusive Shock Acceleration (DSA) mechanism
Diffusive Shock Acceleration (DSA), developed by Blandford & Ostriker (1978) and Bell (1978),is the most fundamental theory for CR acceleration in SNRs. In this theory, we consider the casewhere CRs travel back and forth through the shock wave of the SNR many times, as shown in theconceptual diagram in Fig. 1.4. Here, the following conditions are imposed as the environmentof the shock wave of SNR:CHAPTER 1. INTRODUCTION

7

oc
Sh
kf
nt
ro

Supernova
explosion

post-shock(downstream)

pre-shock(upstream)

vd, ρdin shock rest framevu, ρuin shock rest frame

particle
MHD
waves
MHD
waves

Figure 1.4: Schematic view of the DSA in an SNR. The physical quantities are discontinuous
through the shock front.

CHAPTER 1. INTRODUCTION

8a) The shock compresses the fluid downstream with a compression ratio of r (= vu /vd , where
vu and vd are upstream and downstream velocities, respectively, in the shock rest frame).
b) There are no collisions among particles (commonly known as collision-less shock waves).
c) Particles are isotropically scattered by the Alfvén waves.
d) Conditions b and c allow particles to cross the shock front many times.
e) The shock front can be treated as a plane.
f ) The velocity of the particle can be approximated to the speed of light, and sufficiently fasterthan the shock velocity.
Suppose first that a particle with an energy of E1 crosses the shock wave plane at an angle ofθ1 from upstream to downstream. The energy of this particle in the rest system of downstreamis:E1′ = ΓE1 (1 − βcosθ1 ),(1.5)
where Γ is the Lorentz factor and β is the shock velocity normalized the speed of light (β =vu −vdc ). Once the particle is scattered by a Alfvén wave and crosses again the shock front at anangle of θ2 , the particle energy in the rest system of upstream is:
E2 = ΓE1′ (1 + βcosθ2 ).

(1.6)

Hence, the energy gained during one round trip on the shock wave is:∆E =E2= Γ2 (1 − βcosθ1 + βcosθ2 − β 2 cosθ1 cosθ2 ).
E1

(1.7)

Next, we find the mean values for the two angles, θ1 and θ2 . Eq. 1.7 is written by:∆E = Γ2 (1 − β⟨cosθ1 ⟩ + β⟨cosθ2 ⟩ − β 2 ⟨cosθ1 ⟩⟨cosθ2 ⟩).

(1.8)

If a particle moves at an angle of θ with a velocity of v, the particle approach the shock wave
surface at a speed of v × cosθ, which corresponds to the probability of crossing the shock wave
surface. Therefore, we can obtain the mean values by of weighting with cosθ. The mean values
in Eq. 1.8 are found to be:Rπ
R0
2 × π/2 cos2 θ1 sinθ1 dθ
−2 × −1 x2 dx
2
Rπ
=
=−
⟨cosθ1 ⟩ =R0
3
2 × π/2 cosθ1 sinθ1 dθ
−2 ×
xdx

(1.9)

−1

2
⟨cosθ2 ⟩ = .
3(1.10)

where the integration ranges for θ1 and θ2 are from +π/2 to +3π/2 and from −π/2 to +π/2,
respectively. Using Eq. 1.9, Eq. 1.8 can be rewritten by444∆E = Γ2 (1 + β + β 2 ) ∼ 1 + β.
393

(1.11)

where we assume β ≪ 1 and thus the second-order term and higher are ignored. If a particle
with initial energy of E1 makes n round trips on the shock wave, the total energy acquired is

CHAPTER 1. INTRODUCTION
therefore:

9En∼ (1 + ϵ)n ∼ exp(ϵn) (for n → ∞).
E1(1.12)

where ϵ = ∆E − 1 = 4β/3.
In other words, the energy of the particles is amplified (reduced) due to the Lorentz boost
when particles cross the shock wave from upstream (downstream) to downstream (upstream),
resulting in the energy balance of the round trip being positive. This acceleration mechanism iscalled “Fermi first-order acceleration” because Fermi pioneered this concept (Fermi, 1949) andbecause the energy gained is proportional to the first order of the shock velocity.
This mechanism can explain the CR spectrum with the non-thermal distribution. To obtain
the energy spectrum of particles escaping from this acceleration region, we evaluate the probability of particles escaping from downstream. Since the particle flux (Fd ) escaping from downstream
can be regarded as the amount of particles flowing along with the fluid of downstream, the flux
can be written as:
Fd ∼ Nd vd ,(1.13)
where Nd is the number density of particles in downstream. Regarding the flux of particles thatback upstream, the number of particles is first reduced by half because the particles can backonly with the absolute value of the angle between the direction and the shock wave should beless than 90 degrees. We also take a weighted average with cosθ, as in Equation 1.9. As a result,the flux is obtained as:
R π/2
2 0 cosθsinθdθ
N0 c
1
=
.
Fu ∼ Nd × × c
R π/2
2
4
sinθdθ
2(1.14)

0where we assume the velocity of particles can be approximate to the speed of light. Hence, theprobability that particles escape from the acceleration region (Pesc ) is found to be
Pesc = Fu Fd =

Nd vd
4vd
=
.
Nd c/4c(1.15)

The fraction of the particles that do not escape from the shock region during n round tripon the shock surface is:S(n) = (1 − Pesc )n ∼ exp(−Pesc n) (for n → ∞)(1.16)

Since n can be written by ln(∆E)/ln(1 + ϵ) using Eq. 1.12, the S(n) is rewritten by:
S(n) = exp(−Pesc ×

Pescln(∆E)) ∼ ∆E− ϵ ,ln(1 + ϵ)(1.17)

where we use the approximation of ln(1 − ϵ) ∼ ϵ. Finally, the differential flux is written as
follows:
PescdS∝ ∆E − ϵ −1 .
(1.18)
dE
In this way, we can reproduce the power-law form as in the CR spectrum. With the compression

CHAPTER 1. INTRODUCTION

10

ratio (r = vu /vd ), the spectral index of particles (α) is written as:
α=−

Pesc
4ud
3c
r+2
−1=−
×
−1=−
.
ϵ
c
4(vu − vd )r−1

(1.19)

Rankine–Hugoniot equation
At the discontinuous boundary of the perpendicular shock, as shown in Fig. 1.4, the conservationof mass, momentum, and energy leads to the following relation (called “Rankine–Hugoniot
relation”):
ρ u vu = ρ d vd

(1.20)

ρu vu2 + pu = ρd vd2 + pd
11
wu + vu2 = wd + vd2 ,2
2(1.21)(1.22)

where ρ, v, p and w represent the density, velocity, pressure and enthalpy per unit volume,
respectively. Assuming an ideal gas, the enthalpy can be written by
w = Cp T =γp,(γ − 1)ρ

(1.23)

where γ = cp /cv and cp and cv are the specific heat at constant pressure and the specific heatat constant volume, respectively. Assuming that the shock wave has a large Mach number, i.e.,
the upstream pressure is negligible relative to that of downstream (pu /pd ≪ 1), the following
relations is obtained:
γ+1
ρu
γ−1
γ−1
vd =
vu
γ+1

ρd =

(1.24)(1.25)

The compression ratio (r) is then obtained as:
r≡

vu
γ+1
=
= 4 (for γ = 5/3)vd
γ−1(1.26)

where we consider the case of a monatomic gas and γ is thus 5/3. Substituting this into Eq 1.18,
the spectral index of particles (α) is found to be −2.0.
Maximum acceleration energy by DSA
In Sect. 1.1, we evaluated the maximum acceleration energy determined by the Lorentz force
using Eq. 1.40. Here, we again evaluate it by considering the timescale of CR escape under theenvironment of the shock region. The timescale of particle acceleration is written as:
tacc = ∆T × ϵ,

(1.27)

CHAPTER 1. INTRODUCTION

11

where ∆T is the time taken for one round trip of the shock wave and can be estimated by adding
up the time that particles stay in the upstream and downstream. Hence,∆T =4Du 4Dd
+
,
vu c
vd c

(1.28)

where D is the diffusion coefficient and can be expressed as:D=

λmfp c,3(1.29)

using the mean free path of the particle (λmfp ). The mean free path is determined by themagnetic field and can be expressed using the Larmor radius (RL ) as:
λmfp = ηRL , with η ≡

B
δB

2
≥ 1,

(1.30)

where η, called gyro-factor, is the parameter of strength of the magnetic field turbulence. Smallerη indicates stronger turbulence, especially the limit of η → 1 is called the Bohm limit. Notethat the gyro-factor of SNRs using X-ray observations is estimated to be 1–10 (Tanaka et al.,2020; Tsuji et al., 2021), indicating the particle acceleration in SNR is close to the Bohm limit.
By assuming that r = 4 as in Eq. 1.26 and Du = Dd , tacc can be rewritten by:tacc =

20 cηRL
2 .
3 vsh(1.31)

If the SNR shock radius can be approximated as R ∼ tacc × vsh , then the maximum accelerationenergy is obtained as:
pmax ∼ 0.5 PeV


Z
vsh
η 104 km s−1



B
10 µG



R
10 pc


.

(1.32)The maximum energy is proportional to the shock velocity as well as the magnetic field strengthand the source radius.

1.2.2

CR propagation in the Galaxy

The spectral index of CRs at SNRs is expected to be −2.0 as derived in Sect. 1.2.1, while thatof CR spectrum observed on the Earth is −2.7. The difference in indices would be explained byconsidering the propagation effect in the Galaxy. In this section, we introduce the observationof energy dependency on the diffusion.
The steady-state transport equation can be written by:∂Ni= (Injection) + (Decay from heavier nuclei)∂t− (Loss due to collision) − (Loss due to decay)
X
vρ
vρσi
1
= Qi +
σk→i Nk −
Ni −
Ni
mp
mp
Γτlifek>i(1.33)(1.34)(1.35)

CHAPTER 1. INTRODUCTION

12

where each parameter means the follows.
i, k: the particle species
N : the number density of CR particle
v: the velocity of CR particle
ρ: the matter density in the Galaxy
mp : the rest mass of proton
σk→i : the decay probability from heavier nuclei
σ: the cross section between the particle and the matter in the Galaxy
Γ: the Lorentz factor
τlife : the lifetime of particle
Here, we assume that boron (B) is not elementally synthesized in space (QB = 0) and generated
only by decay from carbon (C) and oxygen (O). We also assume that the lifetime of boron is
sufficiently long to survive without decay during propagation (i.e. 1/(Γτlife )Ni ∼ 1). Once thefraction of CR escaping from the Galaxy is written as exp(−t/τesc ), Eq. 1.33 becomes:vρ
NB=(σC→B NC + σO→B NO )
τescmp

(1.36)

Furthermore, we can transform this equation using NO ∼ NC :λesc
σC→B + σO→B
NB=
NC
1 − λesc /λB
mp(1.37)

where λB (≡ mp /σB ) is the mean free path of the boron in space and λesc (≡ vρτesc ) is thepath length of escaping CR particles. In Eq. 1.37, λB , σC→B , and σO→B can be measured withground-based experiments. λesc is expected to depend on the kinetic energy (velocity) of theparticles, and follows a relation: λesc ∝ E δ . Once we can evaluate the energy dependence of λesc ,we can estimate the effect of CR diffusion on the observed CR spectrum. The observed energy
dependence of the CR boron and carbon abundance ratios is shown in Fig. 1.5. The fitted result(λesc ∝ E −0.6 ) indicates that the index is softened by ∼0.6 due to the CR diffusion (GarciaMunoz et al., 1987). Consequently, the spectral index of CRs arriving from SNR to the Earthcan be estimated to be: α + δ ∼ −2.6, which is comparable to the observed spectrum.
The recent result from more precise observation of secondary-to-primary ratio (e.g., Boronto-Carbon ratio) indicates a flatter slope of δ = −0.333 ± 0.014stat ± 0.005syst (Aguilar et al.,2016) than −0.6. This agrees well with expectation from the Kolmogorov turbulence theory,in which δ is asymptotically equal to −1/3 (Kolmogorov, 1991). However, this also suggeststhat the effect due to CR diffusion alone cannot explain the discrepancy in the spectral slope(∼ 0.7) between the CR spectrum measured at the Earth (∼ 2.7) and the injection spectrumat the CR source (∼ 2.0; expected by DSA) as above stated. This discrepancy would imply
an injection spectrum ∝ E −2.4 , which is preferable from the point of view of anisotropy (Blasi& Amato, 2012). On the other hand, if considering spectral hardening observed in the theCHAPTER 1. INTRODUCTION

13

Figure 1.5: B/C ratio spectrum, taken from Garcia-Munoz et al. (1987). The black lines show
the model curves (an exponential with X0 ∝ E −0.6 for > 1 GeV) using different levels of solar
modulation.
elemental spectra (Ahn et al., 2010) using the broken power-law function:
D(E) = D0Eδ,(1 + (E/Eb )∆δ/s )s(1.38)

instead of
D(E) = D0 E δ ,(1.39)

δ can be estimated to be 0.43–0.53 (Génolini et al., 2017, 2019), which relaxes the disagreement
and is close to a Kraichnan turbulence spectrum (δ = −1/2; Iroshnikov, 1964; Kraichnan,1967). Such spectral hardening could be reproduced if assuming an additional CR source near
the Earth (Kawanaka & Lee, 2021). However, it is still controversial and requires further investigation of CR diffusion, especially around SNRs.

1.3

Lights from CRs in SNRs

The interstellar magnetic field bends the pass of accelerated particles. The Larmor radius of a
proton (Z = 1) with an energy of E is expressed with:
RL ∼ 1 pc

E
3 × 1015 eV



B
3 µG

−1
,

(1.40)

CHAPTER 1. INTRODUCTION

14

Magnetic field

Protonspp→π0CR
sourceπ0-decay γ-ray

Gases

Bremsstrahlung

ObserverElectrons
Photon field (e.g., CMB) Inverse Compton scattering

Synchrotron

Magnetic field

Figure 1.6: Illustration of the γ-ray emission processes around the CR source. Although CRs
lose the information of their acceleration site, the γ-ray emissions would tell us them.
where B is the magnetic-field strength, normalized by the average interstellar value of 3 µG.
This indicates that a PeV proton cannot travel straight for >1 pc without affections of themagnetic field and thus lose the information where they were accelerated. However, CRs also
emit light via several processes as described in Sections 1.3.1 through 1.3.4 (see Fig. 1.6 forschematic illustration). Since light travels straight through space without affection by magnetic
fields, its observation allows us to identify the acceleration site.

1.3.1

Neutral-pion decay

The protons collide with nuclei in the ISM, then produce π 0 as follows:p + p → p + p + π0π+π−.

(1.41)

98% of π 0 decay accompany emitting a pair of γ rays as follows:
π 0 → 2γ.

(1.42)

CHAPTER 1. INTRODUCTION
80
70

σinel [mb]60

15

PDG data
Fit−formula
Kelner(2006)
Kamae(2006)

50
40
30
20
10
0 −1
0
10
10

1

10

2

10

3

4

10 10
Tp [GeV]

5

10

6

10

7

10

8

10

Figure 1.7: Left: Cross section for the proton-proton collisions (Kafexhiu et al., 2014). The
abscissa is the kinetic energy of the proton in the laboratory frame. Black points indicate
measured data and the red line shows the fitted curve for the points. Right: Example of theenergy spectrum of γ-ray produced via the π 0 -decay process (Funk, 2015). α means the assumed
index of the differential proton spectrum.
The energy of γ-ray emission by this process can be written as:
1Eγ = mπ Γπ (1 + βπ cosθ)2
1
1
Eγ, min = Eπ (1 − βπ ), Eγ, max = Eπ (1 + βπ ),22(1.43)(1.44)

where mπ is the rest mass of π 0 (∼ 135 MeV), βπ is the velocity of π 0 , Γπ is the Lorentz factor,and θ is the angle between the direction of π 0 and the γ ray. The fraction of proton kinetic
energy transferring to π 0 is κπ0 ∼ 0.17 (e.g., Funk, 2015), and then the energy fraction of γ-ray
to proton is ∼ 0.09 because the γ-ray energy should be half of π 0 .
The energy loss rate of the π 0 -decay process is given by


dE
dt


π0

∼

κπ0 Ep∆Einel,== κπ0 Ep cnσpp
inel )∆t1/(cnσpp(1.45)

inel is the inelastic proton-proton cross section, and n is the average density of target
where σpp
protons. Hence, the timescale of energy loss due to the neutral pion decay is n −1
1
E
8
τπ0 =
∼
∼ 1.6 × 10 yr
inel
dE/dt
1 cm−3κπ0 cnσpp

inel
σpp
40 mb

!−1
.

(1.46)

inel as a function of energy. The energy threshold of thisThe left panel of Fig. 1.7 shows σpp
γ-ray emission process is determined by the rest mass of π 0 . Hence, the γ-ray spectrum has acharacteristic cutoff at ∼ 100 MeV. The right panel of Fig. 1.7 shows the energy spectrum ofγ rays produced due to the proton-proton interaction. This spectrum has a “bump” structurein the energy band from 100 MeV to 1 GeV and its index after the bump is equal to that ofproton.

CHAPTER 1. INTRODUCTION

16

Note that, CR protons can make not only π 0 and also the charged pions (π + , π − ) by theinteraction with nuclei in the ISM.
p + p → p + n + π+

(1.47)

p + p → n + n + 2π +(1.48)

p + He → 4p + n + π .

(1.49)−

4

π + and π − decay with a lifetime of 2.6 × 10−8 s, and the muon pairs (µ+ , µ− ) are produced.
Furthermore, these muons then decay with a lifetime of 2.2 × 10−6 s, producing electron pairs(e+ , e− ). The Synchrotron emission (the detail of which will be described in 1.3.2) from theseis shown in the right panel of Fig. 1.7.

1.3.2

Synchrotoron

Charged particles moving at relativistic speeds radiate electromagnetic waves when accelerated
by magnetic filed. This radiation is called synchrotron radiation. The energy loss rate of charged
particles, is expressed by4dE= σT cβ 2 Γ2 UB ,(1.50)
dt sync 3
where Γ is the Lorentz factor, σT = 8πr02 /3 = 0.665 × 10−28 m2 (r0 ≡ e2 /(me c2 ); classicalelectron radius) is the Thomson scattering cross section, and UB = B 2 /(8π) is the energy
density of the magnetic field (B). Assuming they have the same energy (E) and substituting
σT ∝ 1/m2 and Γ = E/(mc2 ), the synchrotron radiation intensity ratio of electrons to protons
is found to be


Pp
me 4
=
= 9 × 10−14 .
(1.51)
Pe
mp
Since the synchrotron radiation of protons is much lower than that of electrons, the energy loss
due to its radiative cooling is small, indicating the protons can be accelerated to higher energy
than electrons. Using Eq. 1.50, the timescale for synchrotron cooling can be written as5τsync = 1.25 × 10 yr

Ee
1 TeV

−1 

B
10 µG

−2
(1.52)

where Ee is the energy of electrons. In addition, the average energy of a photon from synchrotron
radiation from a monochromatic electron can be written byϵ ∼ 2 keV

Ee
100 TeV

2 

B
10 µG


.

(1.53)

If the energy distribution of particles is represented by a power-law function N (E) ∝ E α , thenthe photon spectrum can also be represented by a power-law function ofϵdn∝ ϵ(α−1)/2 .
dϵ(1.54)

CHAPTER 1. INTRODUCTION

17

Substituting α = 2 expected by the DSA into Eq. 1.54, we can obtain the photon index of 0.5.
This value agrees well with the observations of SNRs in the radio band (Reynolds et al., 2012).

1.3.3

Inverse Compton scattering

High-energy electrons produce γ rays by passing energy when they collide with light. This
process is called inverse Compton scattering (ICs). The energy loss of charged particles due toICs, is expressed by


dEe
4
= σT cΓ2 β 2 Uphoton ,(1.55)
dt IC 3where Uphoton represents the energy density of the radiation field. Thus, the cooling time scaledue to IC is written by5τIC = 3.1 × 10 yr

1.3.4

Ee
1 TeV

−1 

Uphoton
1 eV cm−3

−1.

(1.56)

Bremssthlung

Bremsstrahlung is the process by which γ rays are emitted when a charged particle is accelerated
by the electric field of an atomic nucleus instead of magnetic field considered in the synchrotron
process. The emission of bremsstrahlung by relativistic electrons (E ≫ me c2 ) is expressed by


dEe
dt


=
Brems

3
σT cαZ(Z + 1)n(lnΓ + 0.36)Ee ,2π(1.57)

where α is the fine-structure constant (1/137.036), Z is the atomic number of the target particle, and n is its number density. Thus, assuming Z = 1, the cooling time scale due to thebremsstrahlung is found to beτBrems ∼ 2.8 × 108 yr

1.3.5



n −1
1
.
1 cm−3(lnΓ + 0.36)(1.58)

Contributions to γ-ray flux

Fig. 1.8 shows the comparison of the radiative cooling timescale. Uphoton and B are set to0.26 eV cm−3 and 3 µG, respectively, the former corresponding to the energy density of CMB
and the latter to the average magnetic-field strength in the Galaxy. The dominant coolingprocess for TeV electrons in SNRs is expected to be synchrotron losses because the magnetic
field should be amplified due to the compression by the shock wave.
The contribution to the γ-ray flux from each radiation process can be also compared using
these cooling times scales. We first compare the fluxes at 1 TeV due to neutral pion decay (Fπ0 )and IC (FIC ). γ rays with an energy of ∼ 1 TeV are produced by electrons (protons) with anenergy of Ee ∼ 20 TeV (Ep ∼ 10 TeV) in the IC process with CMB (neutral pion decay) (Gabici& Aharonian, 2016). The flux ratio can be then written as:
1/τIC × Kep
FIC=∼ 0.27
Fπ0
1/τπ0



Uphoton0.265 eV cm−3



−1
n
100 cm−3



Kep
0.01


,

(1.59)

CHAPTER 1. INTRODUCTION

18

Figure 1.8: The typical cooling timescale due to π 0 -decay (red dotted), Synchrotron (blue solid),
IC (green dashed), and Bremssthlung (magenta dot-dashed).

CHAPTER 1. INTRODUCTION

E2dN/dE [erg cm 2s 1]

10 10

19

Assumed spectrum: Np (E) = N0 E −2.0 e−E/0.1 PeV , Ne (E) = Kep × Np (E)π 0 decay (n = 10 cm−3)
Bremsstrahlung (n = 10 cm−3)
IC (Uphoton = 0.26 eVcm−3)
Synchrotron (B = 10 µG)

10 11
10 12
10 13 3
10

100

103

106 109
Energy [eV]

1012

1015

Figure 1.9: An example of the non-thermal emissions from an SNR. The red dotted, blue solid,
green dashed, and magenta dot-dashed lines represent the emissions from π 0 -decay, Synchrotron,
IC, and Bremssthlung, respectively.
where Kep is the electron-to-proton flux ratio and set to 1%, corresponding to the measured ratioof CRs at the Earth (Aguilar et al., 2015b, 2018). Eq. 1.59 implies that the π 0 -decay radiation
is dominant compared to IC at dense gas regions (e.g., molecular cloud). We also compare theflux from Bremssthlung with π 0 -decay emission. As with Eq. 1.59, using Eq. 1.46 and 1.58, theBremssthlung-to-π 0 flux ratio can be written as:
FBrems
∼ 0.1
Fπ0



Kep
0.01


.

(1.60)

This does not depend on the gas density but only on the electron-to-proton flux ratio (Kep ).
We calculate possible non-thermal emission spectra by primary particles accelerated at an
SNR using the naima package (Zabalza, 2015), as shown in Fig. 1.9. Here we assume thatthe proton spectrum follows the power law with exponential cutoff function, and the electron
spectrum has the same index and cutoff as the proton, with only normalization scaled by Kep
of 1%.

Radio (red) + Optical (yellow) + X-ray (blue)

20
Dec (deg)

CHAPTER 1. INTRODUCTION

TeV γ

50
40

-41.6

30
20

-41.8

10
0
-42

-10
-20

PSF
-42.2

X-ray synchrotron-42.4-226.3815h04m15h02m-225.03
RA (hours)

Figure 1.10: The sky image of SN 1006. Left: Composite view of the remnant of supernova
SN 1006, including X-ray data in blue from the Chandra Observatory, optical data in yellowish
hues, and radio image data in red. Image Credit: NASA, ESA, Zolt Levay (STScI). Right:
TeV γ-ray image of SN 1006 (Acero et al., 2010). The white contours correspond to a constant
2–4.5 keV X-ray intensity as derived from the XMM-Newton flux map and smoothed to the
H.E.S.S. point spread function, enclosing respectively 80%, 60%, 40% and 20% of the X-ray
emission.

1.4

Overview of the SNR observations and PeVatron search

To elucidate the CR origin, search for the PeV proton accelerator “PeVatron” has been conductedthus far (see, e.g., Cristofari, 2021, for a review). This section describes the overview of SNR
observations (Sect. 1.4.1) and 100-TeV γ-ray observations (Sect. 1.4.2) for the PeVatron search.

1.4.1

Observations for SNRs

In the history of SNR observations, Koyama et al. (1995) first discovered the synchrotron X-ray(non-thermal) emission in an SNR named SN1006. The follow-up observation detected TeV
γ rays emitted from the same region as the synchrotron emission (Acero et al., 2010). Thesesuggest an electron acceleration up to ≳ 10 TeV in an SNR. Fig. 1.10 shows the skymaps
of SN1006. The non-thermal emission (i.e., γ rays and synchrotron X-rays) regions spatiallycoincide with the shock front of the SNR, indicating that an efficient CR acceleration occurs
there. In addition, the detailed observations of SN 1006 have revealed a filament structure witha width of < pc of SNR shock waves (Bamba et al., 2003). Using the filament size and themagnetic-field strength derived by the X-ray observation, the maximum energy of electrons isestimated to be 37+4
−7 TeV.
The GeV γ-ray observations of SNRs IC 443, W44, and W51C with AGILE/Fermi-LAT
revealed a bump structure that appears only in the π 0 -decay spectrum (Fig. 1.11), providing

CHAPTER 1. INTRODUCTION21

IC 443W44
10-10

E2 dN/dE (erg cm-2 s-1)

E2 dN/dE (erg cm-2 s-1)

10-10

10-11

Best-fit broken power law
Fermi-LAT
VERITAS (Acciari et al. 2009)
MAGIC (Albert et al. 2008)
AGILE (Tavani et al. 2010)
π0-decay
Bremsstrahlung
Bremsstrahlung with Break

10-12

108

109

1010

Energy (eV)10-11Best-fit broken power law
Fermi-LAT
AGILE (Giuliani et al. 2011)
π0-decay
Bremsstrahlung
Bremsstrahlung with Break

10-12

1011

1012

108

109

1010

1011

1012

Energy (eV)Figure 1.11: Gamma-ray spectra of IC443 (Left) and W44 (Right). Taken from Ackermann et al.
(2013). In each panel, the solid line represents the π 0 -decay emission, while the dashed and dashdotted lines represent the electron bremsstrahlung emissions without/with the low-energy breakat 300 MeV in the electron spectrum.
evidence for proton acceleration in the SNR (Ackermann et al., 2013; Jogler & Funk, 2016;
Giuliani et al., 2011; Cardillo et al., 2016).
TeV γ rays have already been detected at ≳ 20 SNRs2 , and these observations can give alower or upper limit to maximum particle energy at SNRs. However, as mentioned in Sect. 1.3,
there are multiple mechanisms for γ rays, thus we first need to identify the parent particle(i.e., electron or proton) and its emission mechanism. This identification can be approached
by using γ-ray morphology. The fluxes of hadronic and electron bremsstrahlung emissions areproportional to the proton density, while IC does not depend on it (see Sect. 1.3). Hence, We
can distinguish whether the radiation process is IC with the spatial correlation between thedense gas and the γ-ray emission. Moreover, once we assume Kep = 1% (corresponding to
the measurement near Earth), we can consider that hadronic emission is dominant rather thanBremsstrahlung, as in Eq. 1.60.
Such spatial correlation between dense gas and γ-ray emission has been reported in several SNRs. For example, TeV emissions have been found from molecular clouds around SNR
W28 (Aharonian, 2008a), as shown in Fig. 1.12. Since such middle-aged SNR with an age of
20–30 kyr is not expected to accelerate CRs efficiently, the TeV emissions are interpreted byconsidering re-acceleration due to shock-cloud interaction and/or CRs escaped from SNRs in
the past (Uchiyama et al., 2010; Cui et al., 2018). Spatial correlations between dense gas andTeV γ rays have been found also in the young SNRs (e.g., Fukui et al., 2012; Sano et al., 2019),
but the angular resolution of γ-ray observations is still insufficient to spatially separate theemissions from the shell and the dense gas. (e.g., Abdalla et al., 2018a). The latest researchusing the three-dimensional likelihood analysis of energy and morphology found the possibilityof separating them (Fukui et al., 2021), and the next-generation TeV observatory is expected togive more reliable results (Acero et al., 2017). However, its maximum particle energy is at most≲ 100 TeV in any cases and the existence of PeV protons in the SNR has not been confirmed.
2

http://tevcat.uchicago.edu/

Dec J2000.0 (deg)-23oCHAPTER 1. INTRODUCTIONNANTEN

22

12

CO(J=1-0) 10-20 km/s
100
W 28 (Radio Boundary)

HESS J1801-233

80o-23 30´

l=+7.0o60

GRO J1801-2320
l=+6.0o40

G6.1-0.6-24

o6.225-0.569

b=0.0o20b=-1.0oA
HESS J1800-24018h03mW 28A2BC

RA J2000.0 (hrs)18h0Figure 1.12: Skymap of molecular clouds measured with 12 CO (J = 1 − 0) line emission in thevicinity of SNR W28, overlaid with green counters of the TeV emissions (Aharonian, 2008a).
The thin-dashed black circle represents the approximate radio boundary of the SNR.

dJ/dE × E2.7 (G

o

o

o

(a) 25 <l <100 , |b |<5
CHAPTER 1. INTRODUCTION

dJ/dE × E2.7 (GeV1.7 cm-2 s-1 sr-1)
dJ/dE × E2.7 (GeV1.7 cm-2 s-1 sr-1)

10-5 -2
10

10-1

23

1

10

102

10 10

10210

10

3

104

Tibet
AS+MD
Tibet
AS+MD
ARGO-YBJCASA-MIASpace independent CR
Space
IndependentSpace
dependentCR
CR
Pulsarhalo modelSpace
Dependent
CR
10
10-4

-4

o

o

o
o
o o, |b |<5
l <100(b)(a)5025<l<<200, |b |<5
-5

10 -2
10-5
10
10-2

10
10-1

-1

1 1

dJ/dE × E2.7 (GeV1.7 cm-2 s-1 sr-1)

Tibet AS+MD

2

10 10
3

3

10410

4

Energy (TeV)CASA-MIASpace Independent CR

Figure 1.13: Energy spectrum of the diffuse γ rays from the Galactic plane in the region
Space Dependent CR
of |b| < 5◦ , 25◦ < l < 100◦ (Amenomori et al., 2021a). The solid and dashed curves show
the predicted hadronic
emissions from local CR sources with the space-independent and space10-4
dependent models (Lipari & Vernetto, 2018). The dotted curve shows the predicted leptonic
emission from PWNe (Linden & Buckman, 2018).

1.4.2

Sub-PeV γ-ray observationsooo(b) 50 <l <200 , |b |<5

An observation of γ rays above 100 TeV is a reliable (actually only reliable when emission is10-5 -2
3
proven to be hadronic)
methodPeVatron.
Theof the1010-1 to confirm
102
10 recent detection
104
1 the potential
10
Energy (TeV)
diffuse sub-PeV γ rays from the Galactic diffuse plane suggests the presence of PeVatron in the
Galaxy (Amenomori et al., 2021a). Fig. 1.13 shows the energy spectrum of the Galactic diffuse
emission, which suggests the sub-PeV flux can be reproduced with the hadronic emission of local
CR sources rather than leptonic emissions from PWNe3 . Another recent observations with the
LHAASO experiment have revealed 12 sources emitting a sub-PeV γ-ray emission (Cao et al.,2021), as tabulated in Table 1.2. However, the almost sources are possibly originated in a pulsar
wind nebula (PWN) rather than an SNR.
PWNe are generated by the interaction between the circumstellar material and the electronpositron pair wind driven by the pulsar, which composes of a high density, a strong magnetic
field, and high rotation speeds. Its shock wave is considered to accelerate mainly electronsand positrons served by the pulsar, whereas it can also accelerate protons theoretically (e.g.,
Nagataki, 2004). At present, no hint of proton acceleration has yet been obtained, and the latest
Crab Nebula observation with LHAASO did not see it (Cao et al., 2021).

1.5

Research objectives and thesis structure

From the standpoint of believing the paradigm that SNRs are PeVatrons, we can interpret why
a PeVatron has not been discovered yet, with the following two scenarios:
3

Fang & Murase (2021) argued that the leptonic emissions from PWNe still can reproduce the observed subPeV γ-ray spectrum.

CHAPTER 1. INTRODUCTION

Table 1.2: The 100-TeV γ-ray sources detected with LHAASO (Cao et al., 2021).
Name
Possible Origin
Type
Emax [PeV]J0534+2202 PSR J0534+2200
PSR
0.88+0.11−0.11
PSR J1826−1334
PSR
J1825−1326
0.42+0.16−0.16
PSR J1826−1256
PSR
PSR J1837−0604
PSR
J1839−0545
0.21+0.05−0.05
PSR J1838−0537
PSR
J1843−0338 SNR G28.6−0.1
SNR0.26+0.16−0.10
PSR J1849−0001
PSR
J1849−0003
0.35+0.07
−0.07
W43
YMC
SNR G40.5−0.5
SNR
PSR
J1908+0621 PSR J1907+0602
0.44+0.05−0.05
PSR J1907+0631
PSR
PSR J1928+1746
PSR
PSR
J1929+1745 PSR J1930+1852
0.71+0.16−0.07
SNR G54.1+0.3
SNR
PSR J1958+2846
PSR
J1956+2845
0.42+0.03−0.03
SNR G66.0−0.0
SNR
PSR J2021+3651
PSR
J2018+3651
0.27+0.02−0.02
Sh 2−104
HII/YMC
Cygnus OB2
YMC
J2032+4102 PSR J2032+4127
PSR
1.42+0.13−0.13
SNR G79.8+1.2 SNR candidate
J2108+5157––
0.43+0.05−0.05
SNR G106.3+2.7
SNR
J2226+6057
0.57+0.19−0.19
PSR J2229+6114
PSR24

CHAPTER 1. INTRODUCTION

25(1) Some sources in Table 1.2 have not yet been identified well with a counterpart of the 100TeV emission, and they thus may be SNR PeVatrons. Precise γ-ray observations will play
an essential role in identifying γ-ray origins, leading to the discovery of SNR PeVatrons.
(2) We can consider the scenario as follows. An SNR is capable of accelerating particles up to
PeV only when it is younger than 103 yr (Ptuskin & Zirakashvili, 2003, 2005), while suchhigh-energy particles will escape from the SNR at an early stage (Gabici & Aharonian,2007), reducing the opportunity of discovery of PeVatron. To validate this scenario, weneed to measure the time evolution of the maximum acceleration energy at SNR.
In this thesis we explore the possibilities of the two scenarios described above. For the
first point, we focus on LHAASO J2226+6057, which is expected to associate with the SNR
G106.3+2.7 and thus one of the most promising PeVatron candidates. This SNR has a cometary
shape which can be divided into a head and a tail region with different physical conditions. However, it is not identified in which region the 100 TeV emission is produced due to the limited
position accuracy and/or angular resolution of existing observational data. Additionally, it remains unclear whether the origin of the γ-ray emission is leptonic or hadronic. We therefore
carry out new observation for this source using the MAGIC telescopes, which is capable of observing γ rays with a better angular resolution than previous studies. The MAGIC telescopeused for the observations is introduced in Chapter 2, and the MAGIC data analysis for SNR
G106.3+2.7 are described in Chapter 3. The origin of the γ-ray emission is discussed in Chapter 4. For the second point, we propose a new way to measure the time evolution of the SNR, as
described in Chapter 5, and apply this method to observed data. We select SNR HB9 as the firstsample since this SNR has some advantages for this test, which we will mention in Chapter 5.
We finally discuss the feasibility of the SNR paradigm in Chapter 6 based on these results, andsummarize this work in Chapter 7.

Chapter 2

The MAGIC γ-ray telescope
To identify the γ-ray origin of PeVatron candidate SNR G106.3+2.7, further γ-ray observations
with better angular resolution than the previous studies are required, as mentioned in theprevious chapter. In this work, we observe this source with the MAGIC (Major AtmosphericGamma Imaging Cherenkov) telescope, which is one of the Imaging Atmospheric Cherenkov
Telescopes (IACTs) and achieves the highest sensitivity and better angular resolution than 0.1◦
at the TeV energy band. This chapter first introduces the Imaging Atmospheric Cherenkov
technique used in the MAGIC telescope and describes in detail the configuration of the MAGIC
telescope and the analysis method.

2.1

Imaging Atmospheric Cherenkov Telescope (IACT)

2.1.1

Air shower phenomenaLight in the ultraviolet and higher bands cannot be directly observed on the ground due toabsorption by the atmosphere of the Earth. Basically, the observations in their bands are
therefore carried out using telescopes onboard satellites. However, detecting high-energy γ rays
requires a large effective area (because of the infrequent arrival), which is challenging to achieve

Figure 2.1: The two MAGIC telescopes. Image credits: Daniel Lopez/IAC

26

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

27

(b)

(a)
Atmosphere

γ

E=E0

R=0

X0
e-

e+

E0/2

Pair production

2X0
E0/4

e+

γ

γ

e-

Bremsstrahlung

3X0
E0/8
4X0
…

E→EcritFigure 2.2: (a) Sketch of the electromagnetic shower phenomenon. E represents mean energies
of particle or photon. X0 is the radiation length of electron and positron, or a mean free pathof a γ ray, with a dimension of the length. (b) CORSICA simulation of an electromagnetic air
shower, initiated by a γ ray of 100 GeV (https://www.iap.kit.edu/corsika/). Image credit:
Fabian Schmidt (University of Leeds, UK).
with a detector small enough onboard satellites. To observe such high-energy light, an innovative
method was then proposed to use the atmosphere as a detector.
Once γ rays enter the atmosphere, the electron pair production occurs. The electrons thencause bremsstrahlung, producing γ rays, which causes again the pair production. The sequenceof electron pair production and bremsstrahlung is repeated with the energy losses until theelectron energy reaches the critical energy of Ecrit ∼ 83 MeV, where electrons cannot cause
another bremsstrahlung due to ionisation losses (see Fig. 2.2 for a sketch). This phenomenon
is called an electromagnetic air shower. The longitudinal development of electromagnetic air
showers approximately follows the semi-empirical model proposed by Greisen (1956), and isexpressed with the so-called shower age (s) parameter. The shower age is defined as:
s=

E0
3t
with y = ln,
t + 2y
Ecrit

(2.1)

where E0 is the primary energy and t is the shower depth in units of radiation length (X0 ∼36.7 g cm−2 for a dry air). The shower development begins at s = 0, reaches a maximum at
s = 1 (i.e., t = y), and then decays. The average number of electron pairs above the critical

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

28

Figure 2.3: Longitudinal development of an electromagnetic air shower in a hydrostatic, isothermal atmosphere, following the Greisen semi-empirical model (Greisen, 1956). The shower size,
characterized by the number of secondary electrons in the shower (Ne ), is plotted as functionof depth. The dashed lines corresponds to equal-age curves. The dashed lines corresponds toequal-age curves.
energy (Ecrit ) in the shower can be approximated by
0.31
Ne (t, E0 ) = √ et−1.5tln(s) .
y

(2.2)

Figure 2.3 shows the longitudinal development of the γ-ray initiated air shower calculated using
Eq. 2.2. For the IACT observations, it is necessary to capture the longitudinal development in
the FoV. Therefore, an altitude of several thousand kilometers is suitable for the location of thetelescope observing TeV γ rays.

2.1.2

Cherenkov emission

Once a charged particle travels faster than the speed of light in that medium, the medium ispolarized and makes dipole radiation, as shown in the left panel in Figure 2.4. By Huygens’
principle, the phases of the electromagnetic wave are aligned, resulting in a forward-beamed
emission called “Cherenkov Radiation”. Following this, electron pairs in an air shower emit the

+

+
-

+-

-

+ + - ++ - ++
- +-

+

-

+-

+

+ - +- +- + ++

+-

-

++-

-

+

-

-

-

-

v<c

+

- +

-

+

- +
+ ++ - +-

++

+-

+-

-

+- +-

+

+

+-

-

Θ

+-

++

-

-

+
+-

+

+-

+

+

- +- +-

-

k

+

+

29

-

++ + -+ + - +-- +-CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

v>c

v<c

v>c

Figure 2.4: Illustration of the polarization of the medium induced by the crossing of a relativistic particle (Left) and Construction of Cherenkov wave-front (Right), which are taken fromde Naurois & Mazin (2015).
Cherenkov light.
The radiation angle of the Cherenkov light (θ) depends on the velocity of the charged particle
and the refractive index of the medium (n), and can be written by:cosθ =

1(for βn ≥ 1).
βn

(2.3)

The typical refractive index of the atmosphere is n ∼ 1.0003, and then θ ∼ 1.4◦ . With theoverlapping of the Cherenkov lights from an air shower, its shape reproduces the longitudinaldevelopment of the shower.
The number of Cherenkov photons (Nph ) per unit pass length and unit wavelength can be
expressed as:
d2 Nph2παZ 21=(1 − 2 2 ),
(2.4)
2
dxdλ
λ
β n
where α is the fine structure constant and Z is the atomic number of the incident particle. On
the other hand, the photon spectrum will be affected by absorption in the air. Absorption ismainly due to O3 , O2 , and N2 for the ultraviolet band (< 300 nm) and rotational or vibrationaltransition lines of H2 O and CO2 for the infrared band. The resultant observable spectrum isshown in Figure 2.5 (Benn & Ellison, 1998; Bonardi et al., 2014), which peaks at 300–400 nm.
The density distribution of the Cherenkov light is expected to be uniform within the so-called
light pool, the radius (RC ) of which can be expressed as a product of the maximum height of air
shower and θ. Thus, the radius for a γ ray of 100 GeV is estimated to be RC = 10 km × 1.0◦ ∼110 m. In reality, however, the structure is broader than RC due to the multiple scattering effectin the atmosphere. Fig. 2.6 shows simulations of Cherenkov photon density. Cherenkov photon
densities are almost flat within ∼ 100 m, decreasing at large core distances. In addition, thedensity of Cherenkov photons increases at high altitudes, as with the air shower (Fig. 2.3). At
higher altitudes, the telescopes are closer to the emission region, resulting in a larger density ofphotons over a smaller area on the ground. Hence, the energy thresholds of IACTs would belower. On the other hand, there are several concerns with observations at high altitudes (e.g.,CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

30

Figure 2.5: Cherenkov light and night sky background (NSB) spectra in La Palma (CanaryIslands, Spain) at an altitude of 2200 m (Benn & Ellison, 1998; Bonardi et al., 2014).

Figure 2.6: Lateral distributions of Cherenkov photon density for wavelengths from 300 to
600 nm (Hassan et al., 2017). The solid, dotted, and dashed lines represent results at analtitude of 3500, 2500, and 1500 m, respectively. Left: Cherenkov photon density in linear scaleclose to the core position. Right: Broader core distance ranges and logarithmic scale.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPEMAGIC Telescope

Camera (Front view)

PMT module

Back of Camera

31

Winston Corn

Optical fiber
162 mElectronics & PC

Figure 2.7: The focal plane camera of MAGIC telescope. The bottom-left and bottom-right
pictures are taken from Aleksić et al. (2016). The other pictures are taken from https://magic.mpp.mpg.de.
Hassan et al., 2017). First, IACTs at high altitude will see larger fluctuations in the shower
images because of random particles penetrating down to or close to ground level, which worsens
the determination accuracy of the shower parameters. Also, the shower image tends to leakfrom the FoV of the telescope due to the larger angular difference between the shower axis andthe pointing axis. Thus, the angular resolution and energy estimation accuracy worsen in thiscase. As a consequence of these effect, it is better to install the telescopes at the high altitude toachieve the lower energy threshold, or at the low altitude to improve the telescope performance
regardless of the energy threshold.

2.2

The MAGIC system

In this section, we introduce the MAGIC telescopes used in this study. The MAGIC telescopeis located at 2200 m altitude above sea level at the Observatorio del Roque de los Muchachoson the Canary island La Palma, Spain (28.76◦ N; 17.89◦ W). As shown in Figure 2.7 for theschematic view, the flow of the Cherenkov signal at the MAGIC telescope is as follows: (i) the
MAGIC system focuses Cherenkov lights using parabolic mirrors with a diameter of 17 m to
the focal plane camera, (ii) which converts into charge information and transmits them to thereadout boards in the electronic room via optical fiber with lengths of 162 m. We describe thedetail for (i) and for (ii) in sub-sections 2.2.1 and 2.2.2, respectively.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

2.2.1

32

OpticsIn order to detect low-energy γ rays as well as high energies, the MAGIC system focuses
Cherenkov lights using a parabolic mirror with a relatively large diameter of 17 m. A parabolic
mirror can suppress the arrival time dispersion of reflected photons and thus can provide a
better signal-to-noise ratio, resulting in the achievement of the lower energy threshold1 . The
mirror dishes consist of 250- and 247-segmented mirrors for MAGIC-I and II, the area of whichis ∼ 240 m2 . The small deformations of the mirror dish due to a change of observation zenith
angle can be corrected using the Active Mirror Control system (Biland et al., 2007). The totalreflectivity of the mirrors is maintained at ∼ 85% (Doro et al., 2008).
In addition, the hexagonal shape Winston cones are attached to the front of each photomultiplier tubes (PMTs) of the focal plane camera to minimize dead space between PMTs andto exclude NSB from wide angles to the focal plane.

2.2.2

Photo-sensor and the readout system

Photo-multiplier tubes (PMTs)
Each MAGIC focal plane camera uses 1039 PMTs (R10408 from Hamamatsu Photonics), which
has a hemispherical photo-cathode with a diameter of 1 inch. Fig. 2.8 shows the quantum
efficiencies (QEs) of PMTs installed in the cameras. The maximum QEs are 30–35% in the
wavelength range of 300–450 nm, corresponding to the peak of the Cherenkov light spectrum
shown in Fig 2.5.
The average of the PMT gains is 3.0 × 104 with the high voltage (HV) of 850 V, but the
distribution has a spread within (1.0–6.0) ×104 . Such a spread in the gains is unavoidable
during the manufacturing process, but would make biases in the trigger production and theanalysis. We, therefore, align the gains of each pixel by adjusting the applied voltage of each
PMT and using attenuators for those with high gains. This process is called flat-fielding. The
left and right panels in Fig. 2.9 show the distribution of the HVs applied to PMTs and thecharge distribution of the calibration pulses in the MAGIC cameras after the HV flat-fielding
procedure, respectively. The latter distribution is in good agreement between the MAGIC-I andII, leading the less systematic uncertainties on the energy and flux estimation than 8%. We alsofind two pixels in MAGIC-II that could not be flat-fielded well, because the gain is too low evenat the maximum HV. The signal-to-noise ratio of these two pixels are worse than others but
they can still be used in the analysis.
Readout system
The amplified analog signals from the PMTs are converted into optical signal with VCSELs,
and then transmitted to the electronics room via optical fibers. The high sampling frequencyof ∼1 GHz for the analog-to-digital conversion is required to acquire the waveform data of1

One disadvantage of parabolic mirrors is that they have a large coma aberration, yielding a narrow FoV oftelescopes. For example, VERITAS and H.E.S.S. adopt the Davies-Cotton telescope (which can reduce coma
aberration) to achieve a wide FoV instead of a low energy threshold.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

33

120

MAGIC-I

100

MAGIC-II

Number of pixels

Figure 2.8: Quantum efficiency of PMTs for the MAGIC telescopes as function of wavelength (Nakajima, 2013).

MAGIC-I
102

Entries 1039
7775
Mean
RMS
209.2

10

Entries 1039
Mean
7756
RMS
296.9

MAGIC-II

80
60
40
20
0
500

1
600

700

800

900 1000 1100 1200 1300 1400 1500
High Voltage (V)01000 2000 3000 4000 5000 6000 7000 8000 9000
Mean Signal (ADC counts)

Figure 2.9: The distributions of the applied HVs for the MAGIC PMTs and of charges of
the calibration pulses after the flat-fielding procedure (Aleksić et al., 2016). Left: The solidand dashed lines show the distributions in MAGIC-I and MAGIC-II cameras, respectively. The
highest voltage that can be applied to the PMTs is 1250 V. Right: The red and blue data show
the results of MAGIC-I and MAGIC-II cameras, respectively.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

34

Cherenkov light with a time width of a few ns. Such high-speed sampling can also allow to
reduce the NSB contamination with an arrival frequency of 10–100 MHz. For this requirement,
the MAGIC system uses the Domino Ring Sampler version 4 (DRS4)2 as well as the ADC(14 bit resolution), achieving a sampling frequency of 2.05 GHz. The DRS is an ASIC usinga switched capacitor array, which can store charge information in thousands of capacitors insequence and allow the relatively low-speed (32.5 MHz used in MAGIC) digitization with ADC
at the latter part. The DRS was originally developed for the MEG experiment (Ritt et al.,2010), but is currently used in many experiments (e.g., Anderson et al., 2016; Nozaki et al.,2020) because DRS enable us to operate with lower power consumption compared to FlashADC. The calibration method for DRS4 in the MAGIC telescope was established by Sitarek
et al. (2013).
Trigger system
In the trigger generation process at the DAQ system, the MAGIC telescopes make signals inthree stages (Aleksić et al., 2016), aiming to avoid accidental noise data and reduce the deadtime. First, discriminators find pixels in which the amount of charge exceeds the criteria. Oncethe three or more pixel-wise triggers are generated simultaneously at neighboring pixels, theindividual telescope trigger is generated. Finally, the stereo trigger is generated by finding thecoincidence event between the individual telescope triggers of MAGIC-I and II within a 100 ns
time window. The typical rate of the stereo trigger is set to be 250–300 Hz for normal operations.
Figure 2.10 shows the block diagram of the readout and trigger chain.
In-Situ calibration of the PMT camera
Since the PMT gain varies with temperature, it is desirable to calibrate their response during observations. We perform in-situ calibration using a technique called the “F-Factor” method (Mirzoyan, 1997; Schweizer et al., 2002). This principle is as follows (Bencheikh et al., 1992). The
relation between the average of the output (Qout ) and input charges (Qin ) can be written as:⟨Qout ⟩ = ⟨G⟩⟨Qin ⟩ + P,

(2.5)

where G is the PMT gain and P is the pedestal noise. In the case of a perfect detector,
i.e., where the detector response is one-to-one correspondence between input and output, thestandard deviation of the output charge is determined only by the statistical uncertainty of the
input charge, and thus:p(2.6)
σQout = G Qin .
By considering the uncertainty on the detector response, the observed variance will be2σQ= ⟨G⟩2 ⟨Qin ⟩ + ⟨G⟩2 ⟨Qin ⟩out2https://www.psi.ch/en/drs



σQ1p.e.
Q1p.e.

22+ σelec,(2.7)

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

35

Figure 2.10: Schematic view of the readout and trigger chain of the MAGIC telescopes (Sitarek
et al., 2013). The pixel-wise trigger, the individual telescope trigger, and the stereo trigger are
named L0, L1, and L3, respectively.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

36

where the Q1p.e. is the output charge for a single electron in PMT, σQ1p.e. is its standard deviation,and σelec is the combination of the detector, amplifiers and ADC electronics noise. Once ⟨G⟩⟨Qin ⟩is large enough than the pedestal noise, the terms of P and σelec can be neglected from Eq. 2.5
and 2.7, respectively. Hence, we obtain the following relation:2
σQ= ⟨G⟩(1 +
out

2
σQ
1p.e.

Q21p.e.

)⟨Qout ⟩.

(2.8)

Finally, the PMT gain can be written as:
v

u
2
2
u
σQ
σQ
out
⟨G⟩ = 2, with F = t1 + 21p.e. ,F ⟨Qout ⟩Q1p.e.

(2.9)

where F is the so-called F factor (also called excess-noise factor). In Eq. 2.9, Q1p.e. and σQ1p.e.
can be measured in the laboratory before inserted in the telescopes. As a consequence, we canestimate the PMT gain using the pulse data uniformly illuminated over the entire camera witha higher amplitude than the noise. In the MAGIC operation, the pulse evens generated by thecalibration box, installed the center part of the mirror dish, are taken simultaneously during theobservations of air shower events. The calibration box has a Nd:YAG laser that produces pulsesof 0.4 ns FWHM with a wavelength of 355 nm. This can illuminate light pulses onto the PMT
camera uniformly with a fluctuation of less than 2% (Aleksić et al., 2016). This system are also
used to perform the flat-fielding procedure.

2.2.3

Atmospheric monitoring

As mentioned in Sect. 2.1, the air shower depends on the thickness of the atmosphere, i.e., notonly altitude and zenith angles and also the atmospheric transmission. Since the atmospherechanges from moment to moment, it is necessary to monitor it during observations. For themonitoring, we use Light Detection And Ranging (LIDAR; Fruck et al., 2014). The LIDAR
system shoots a laser signal (wavelength: 532 nm, pulse width: 0.5 ns) to the air and measures
the reflected laser light. The atmospheric transmission can then be estimated by measuring thetime interval between the emission of the laser and the arrival of the reflected light. If a measured
transmission is worse than ∼ 80%, we need to consider the additional systematic uncertainty ofthe energy estimation and/or the correction factor (Fruck et al., 2014; Hahn et al., 2014).

2.2.4

Drive system

The mount of the MAGIC telescope consists of a space frame of carbon-fiber reinforced plastic
tubes, achieving the comparatively low weight of 72 tons per telescope. This allows a fast
movement of 180◦ (360◦ ) in azimuth within 20 s (33 s) (Bretz et al., 2009). The systematic
uncertainty on the reconstructed source position is estimated to be less than 0.02◦ (Aleksić
et al., 2016), and thus would not affect analysis results.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

Name
size
width
lengthwl
CoG
CONCtime gradientleakage
MaxHeight
impactNumber of islands

2.3

37

Table 2.1: The main Hillas parameters.
Description
Total charge of selected pixels
length in the minor axis of the shower ellipse
length in the major axis of the shower ellipse
width over length
Center of Gravity in the shower image
Fraction of the image concentrated in the brightest pixels.
slope of the arrival time projection along the major axis
Charge contained in the outer edge pixel of the camera
Height of the shower maximum point from the ground
Crossing point of shower axis on the ground.
Number of non-connected pixels that survived the image cleaning.

Event reconstruction

As the name implies, the data analysis of IACT uses Cherenkov images for event reconstruction.
In the imaging process, the signal is extracted using a sliding window algorithm (Albert et al.,2008a) by finding six consecutive cells (corresponding to 3 ns) selected to be maximized within
the readout window. Also, NSB is removed by taking coincidence within 2–4 neighboring pixels.
We then parameterize features of the images with the so-called Hillas parameters (Hillas, 1985).
The main Hillas parameters are summarized in Table 2.1, some of which are shown in Fig. 2.11,
overlaid on the electromagnetic shower image. Using these parameters, we estimate the energy,
flux, and arrival direction of incident γ rays (Sect. 2.3.1 and 2.3.2) and remove background
events (Sect. 2.3.3). These data analyses are performed with the MAGIC standard analysis
package (Zanin, 2013).

2.3.1

γ-ray energy and flux

As shown in Figure 2.3, γ rays with higher energies make a large air shower, yielding a larger
amount of charge in the shower (size parameter). Fig. 2.12 shows the correlation between
the mean size of a γ-ray event and its energy, obtained with MC simulations. However, theenergy estimation also depends on the incident angle of γ rays and the distance from telescope(corresponding to impact), which thus should be taken into account. We, therefore, estimate theenergy of observed γ rays by evaluating the Hillas parameters with MC data. The evaluationis performed using machine learning, specifically, the Random Forest (RF) method (Breiman,
2001), which has been implemented into the MAGIC analysis by (Albert et al., 2008b).
Figure 2.13 shows the energy resolution and energy bias evaluated with MC simulations. Theenergy resolution has been estimated to be 15–25% for the range of 60 GeV–20 TeV, indicatingthat the estimated energy spectrum may be biased. We therefore compute the response function
using MC simulation data and then derive a true energy by unfolding with the results. The
relation between the distribution of the estimated energy (S) and the true energy (Y ) is expressedCHAPTER 2. THE MAGIC γ-RAY TELESCOPE

38

th

g
len
th

wid

CoG

log( E/GeV )Figure 2.11: Cherenkov image of an electromagnetic air shower, overlaid with the several Hillas
parameters (see Tab. 2.1 for detail).
5000

4.5

4500
4

4000
3500

3.5

3000
3

2500
2000

2.5
1500
1000

2

500
1.5
1.5

2

2.5

3

3.5

4
4.5
5
log( (Size +Size2)/2 / phe )1Figure 2.12: Correlation between the mean size of a γ-ray event and its energy obtained with
MC simulations (Aleksic et al., 2012).

Energy bias and resolution

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

39°

Bias (2010, < 30 )0.35

Resolution°

Bias (2013, < 30 )°0.3

Resolution°

Bias (2013, 30 -45 )Resolution

2000

10000 20000
Etrue / GeV

0.25
0.2
0.15
0.1
0.05
0

-0.05

60

100

200 300

1000

Figure 2.13: Energy resolution and energy bias of the MAGIC telescope (Aleksić et al., 2016).
using the transfer matrix (M ) as follows
Yi =

X
j

Mij · Sj (i = 1, ..., n).

(2.10)

Furthermore, we perform a regularization to suppress the effect of the systematic uncertainty
in the unfolding. Although there are several regularization methods (Albert et al., 2007), thedifference in the results of the spectral index among the methods is estimated to be less than
0.1 (Aleksic et al., 2012). The final results in this thesis are those obtained with the Tikonov
method.
The differential γ-ray flux is defined as:
dNγ
dF
=
dE
dEdtobs dAeff(2.11)

where Nγ is the number of γ-ray events, tobs is the observation time, and Aeff is the effective area.
The effective area is actually obtained by integrating the detection efficiency with azimuthal
direction (Az) and the impact parameter (Pimpact ), and is defined as:
Z 2π
Aeff =ZdAz

0dPimpact ϵ(E, Zd, Az, Pimpact ) × Pimpact ,(2.12)

Note that ϵ depends on the energy and zenith angle (Zd) as well as Pimpact and Az. Using theMC simulation data of γ rays, the effective area can be also computed as follows:
Aeff =

Nsurvived
AMC ,Ninput(2.13)

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

40

×103
170

106

160

104

Zd = 5-15 deg
Zd = 15-25 deg

103

Zd = 25-35 deg
102

Zd = 35-50 deg

10Effective area at 1 TeV [m2]Effective area [m2]

105

150
140
130

∝ 1/cos(Zd)

120
110
100

1

90
102

103
Energy [GeV]

104

10

20

30

40

50

Zenith distance [deg]Figure 2.14: Left: Effective area of the MAGIC telescope as a function of energy. Right:
Effective area at 1 TeV as a function of zenith angle. The red curve shows the fit function of
1/cos(Zd).
where Ninput and Nsurvived are the number of a total simulation events and the survived events
after the analysis cut, respectively, and AMC is the geometrical are used in the simulation.
Fig. 2.14 shows the effective area for point sources derived with the MC simulation. The
zenith dependence of the effective area follows a relation between relative air mass and zenith
and then Aeff ∝ 1/cos(Zd) (the right panel of Fig. 2.14). In observations at the larger zenith
angle, the energy threshold become higher3 but the effective area can be increased, enabling us
to observe higher energy photons.

2.3.2

Arrival direction

In the case of stereo observations (observations with two or more telescopes), the arrival directionof incoming γ rays can be geometrically determined from the crossing point of the major axesof the Hillas ellipse. Although this estimation method has been used conventionally, it may
fail if the Hillas parameters are not determined well (e.g., in the case of low-energy events).
In the MAGIC data analysis, the arrival direction is therefore estimated with the DISP RF
method (Aleksić et al., 2016), which takes into account the time gradient parameter. As shownin Fig. 2.15, the DISP RF method estimates the arrival position as follows:
1 Two estimated points per telescope are obtained using the time gradient parameter.
2 We then find the combination of the estimated positions between MAGIC-I and II so thatthe distance is the shortest
3 We take the average of the two points by weighting the number of pixels in the image.
The estimated angular resolution can be evaluated with MC simulation data or the real
observation data of a point-like source. Fig. 2.16 shows the angular resolution of the MAGIC
3
Fig. 2.14 indicates that an effective area in larger zenith observation is smaller than the lower ones atthe low energy band. We usually set a safety energy threshold to consider the deterioration of the signalto-background ratio due to this effect. The analysis threshold can be approximated by an empirical formula:
74 × cos(Zd)−2.3 GeV (Aleksić et al., 2016).

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

41

B1

B2

MAGIC-I

MAGIC-II

A2

A1

Reconstructed position
with DISP RF method

Reconstructed position
with conventional method

Angular resolution [degree]

Figure 2.15: Schematic view of the DISP RF method. The main axes of the Hillas ellipses
are shown with solid black lines. The red squares represent the estimated positions using thetime gradient parameter. The black and red stars represent the reconstructed position with theconventional and DISP RF methods, respectively.

0.18

68% cont., Crab, 2010, < 30deg

MC

0.16

68% cont., Crab, 2013, < 30deg

MC

0.14

68% cont., Crab, 2013, 30-45deg

MC

0.12
0.1
0.08
0.06
0.04
0.02
0

102

10

3

104
Estimated Energy [GeV]

Figure 2.16: Angular resolution of the MAGIC telescopes as a function of the estimated energy
obtained with the Crab Nebula data sample (points) and MC simulations (solid lines) (Aleksić
et al., 2016). Red (blue) points show the low (medium) zenith angle sample. For comparison the
low zenith angle pre-upgrade angular resolution is shown as gray points (Aleksic et al., 2012).

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

42

T3
T4
T1

T1
T2

Figure 2.17:
Left: Angular resolution of the VERITAS telescopes as a function of
the energy (https://veritas.sao.arizona.edu/about-veritas/veritas-specifications).
Right: Schematic of the VERITAS telescopes as viewed from directly overhead (Perkins &Maier, 2009). The blue lines are for the original array layout, while the red lines are for the newarray layout. The black arrow indicates the relocation of Telescope-1. The small grid squares
are 5 m on a side.
telescopes as a function of an estimated energy. In the lower energy band than ∼ 1 TeV, theangular resolution is determined by the accuracy of the Hillas parameters and, thus, improves
with energy since γ rays with higher energies produce more explicit shower images. At higher
energies above 1 TeV, the pixel size of PMTs determine the resolution limit.
In addition, the resolution is also correlated with the number of telescopes (Ntel ) that can−0.5
observe simultaneously, and roughly proportional to Ntel
(Funk & Hinton, 2009). While Ntelof the MAGIC system is half of VERITAS, the used pixel size is small compared to VERITAS.
The resultant current angular resolution of both telescopes is comparable. Note that the angularresolution in the previous study with the VERITAS telescopes for SNR G106.3+2.7 (Acciari
et al., 2009) is 10–25% worse than the current performance due to the constructed location of
the telescopes. The left panel in Fig. 2.17 shows the angular resolution of the VERITAS telescope
for the several periods. After the previous study, the angular resolution has been improved bymoving the telescope (Perkins & Maier, 2009), as shown in the right panel of Fig. 2.17.

2.3.3

Background estimation

The hadron-initiated air shower
A main background in the shower image analysis is the hadron-initiated shower. CRs arriving
at the Earth produce pions due to their strong interaction with atmospheric nuclei, according
to Eq. 1.41:
p + p → p + p + π0π+π−.
(1.41)

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE
(a)

43

(b)

primary CR(p, α, Fe, ...)

(c)

Atmosphericnucleus

ππ+

µ!µ
"
µ+

π0
nucleons

e+

e+

π0

π±

γ

γ

e-

e-

EM shower

νµ
e+

e+
e-

e-

EM shower

Figure 2.18: (a) Schematic view of a hadron-initiated air shower. (b) Longitudinal developmentof the hadron-initiated showers obtained with CORSICA simulations. (Image credit: Fabian
Schmidt, University of Leeds; https://www.iap.kit.edu/corsika/). (c) Cherenkov image ofa hadron-initiated shower obtained with the MAGIC telescope.
The cross section can be written as:
1/3

1/3

σAA [mb] ∼ 65(Ain + Aair + 1.12)2 ,(2.14)

where Ain and Aair indicate the mass number of the incident CR and the atmosphere, respectively, the latter of which is Aair ∼ 14.5. In this process, π 0 , π + , and π − are generated at aratio of 1 : 1 : 1, respectively. π + and π − then decay due to the weak interaction as follows:
π + → µ+ + νµπ − → µ− + νµ ,(2.15)

and 98% of π 0 decays according to Eq. 1.42:
π 0 → 2γ.

(1.42)The γ rays may produce electromagnetic air showers as described in Sect. 2.1.1. As a result, theCR hadrons can cause air shower phenomena similar to the γ-initiated showers (see Fig. 2.18).
Since the number of the γ-initiated events is much less than the hadron ones (even for brightsources, e.g., Crab Nebula, the fraction is ∼ 10−3 ), it is necessary to distinguish whether theobserved air shower is initiated by hadrons or γ rays.
The main differences between showers initiated by γ rays and hadrons are as follows:
1. In the case of a γ-initiated shower, all the energy of incoming γ rays is converted into
shower development, while most energy of a hadron-initiated shower is lost due to µ and ν
produced by Eq. 2.15. Consequently, the charge in the Cherenkov image (size) of hadroninitiated showers will be less than that of γ rays with the same energy.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

20

44

18

γ -ray (simulation)

16

Hadron (real)# of Events

14
12
10
8
6
4
2
0
0

0.1

0.2

0.3

0.4

0.5
0.6
Hadronness

0.7

0.8

0.9

1

Figure 2.19: The hadronness distributions of γ-ray data generated with the MC simulation andthe hadron data obtained with the real observation.
2 Because of the momentum of the π 0 at the production, the shower spreads significantly not
only in the longitude but also in the lateral direction. The structure of the shower may beasymmetric or have sub-cores within it.
Simulated image of the hadron-initiated shower are shown in the right panel of Fig. 2.18, in which
the difference from the γ-initiated shower (Fig. 2.11) appears well. This difference is also evident
in the Hillas parameters. For example, the hadron-initiated showers have a larger width and wl(width over length) than those of γ-ray showers. Furthermore, the CONC parameter tends to belarge for γ-initiated showers. Thus, we can distinguish between hadron and γ-initiated shower,by evaluating the Hillas parameters. γ/hadron separation is made by the so-called “hadronness”
parameter, which takes a value close to 1 (0) when identified as a hadron(γ)-initiated shower. As
with the estimation of γ-ray energy (Sect. 2.3.1) and arrival direction (Sect. 2.3.2), we use theRandom Forest method (Albert et al., 2007) to obtain the hadronness parameter. The training
data for the machine learning are γ-ray data obtained with MC simulations and hadron dataobtained from actual observations. Figure 2.19 shows the hadronness distributions of γ-ray data
generated with MC simulation and hadron data obtained with actual observations. It indicatesthat the hadron and γ-ray events can be distinguished well, and thus the event selection using
hadronness will be efficient.
Note that, there is a trade-off between reducing the hadron background events and increasing
the gamma-ray statistics. In this work, we select hadronness cut value so that 90% of gamma-ray
events remain for each energy bin by finding with MC simulations.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

45

Camera FoV
OFF

OFF
OFF

ON
Wobble offset

Cameracenter

~0.4 deg
γ-raysourceOFF

OFF

Figure 2.20: The ON (red circle) and OFF (black circle) regions in the Camera FoV.
The wobble-mode observation
As mentioned above, the hadron-initiated shower event is a main background in the IACT
observation. The background reduction using the Hadronness parameter is very useful, but
as shown in Figure 2.19, it cannot eliminate them completely. In addition, electron-initiatedshowers can also be a background event, which also cannot be classified with the Hadronness
cut. Therefore, we estimate the background rate after the Hadronness cut by selecting the OFF
region without including γ-ray sources from the Camera FoV. However, since the sensitivity in
the FoV depends on the distance from the camera center (Aleksić et al., 2016), it is necessary
to choose a location from which the distance is the same. In addition, it is desirable to observe
the ON and OFF regions simultaneously to consider time variations in weather and NSB. For
this reason, we perform IACT observations with the wobble mode (Fomin et al., 1994), in which
the center of the γ-ray source is shifted from the camera center. Fig. 2.20 shows the schematic
view of the ON and OFF regions in the wobble mode. The distance between the camera centerand the γ-ray source is named wobble offset and is typically set to 0.4◦ (a larger wobble offset
may be adapted for extended sources).

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

46

Table 2.2: Summary of Systematic uncertainties in the MAGIC observations (Aleksic et al.,2012). ES, FN, and SL are abbreviations for the energy scale, flux normalization, and spectral
slope, respectively.
Part
Uncertainty
F-Factor
10% ES
atmospheric transmission≲ 10% ESmirror reflectively
8% ES
PMT electron collection efficiency
5% ESLight collection in a Winston Cone
5% ES
PMT quantum efficiency
4% ES
Signal extraction
3% ES
Temperature dependence of gains
2% ES
Charge flat-fielding
2–8% ES FN
Analysis and MC discrepancies≲10–15% FN
Background subtraction
1–8% FN
Broken channels/pixels
3% FN
Mispointing
1–4% FN
NSB
1–4% FN
Trigger
1% FN
Unfolding of energy spectra
0.1 SL
Non-linearity of readout
0.04 SL

2.4

Systematic uncertainties

The systematic uncertainties due to each component are summarized in Table 2.2. The total
effect on the energy spectrum parameters will be also described in Sect. 2.4.1 through 2.4.3.

2.4.1

Flux normalization

The systematic uncertainty for the flux normalization are mainly due to the charge flat-fieldingof the camera, analysis and MC discrepancy, and background subtraction (Aleksic et al., 2012).
The uncertainty can be described as:q
2 + σ2(2.16)
σnorm = σA
SBR ,where σA and σSBR are the uncertainties dut to the effective area and signal-to-background ratio(SBR), respectively. σA has been estimated to be 18% at the low energies (< 100 GeV), 11%
at the medium energies, and 16% at the high energies (≳ 1 TeV) (Aleksić et al., 2016). σSBR is
described as:
0.015
σSBR =.
(2.17)SBR
If the SBR is worse than 1.5%, the uncertainty exceeds 100%.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

47

Table 2.3: Parameters of the existing IACTs.
Telescope
MAGIC
VERITAS
H.E.S.S

location
La Palma (Spain)
Arizona (USA)Namibia

latitude28.76◦ N
31.68◦ N
23.27◦ S

Height [m]2200 m1268 m1800 mMirror Diameter [m]
17 m × 2
12 m × 4
28 m × 1 + 12 m × 4

FoV [deg]
3.5
3.5
3.2, 5.0†† The former is of the 28-m telescope, while the latter is of the 12-m telescope.

2.4.2

Energy scale

The estimate of an energy scale highly depends on the accuracy of the PMT calibration (F-Factor
method, see Sect. 2.2.2) and the atmospheric transmission (see Sect. 2.2.3). This uncertainty is
estimated at 17% at low energies and 11% on the high energies.

2.4.3

Spectral slope

The uncertainty for the energy slope mainly comes from the unfolding accuracy of energy spectraand non-linearity of readout. This uncertainty is estimated to be ±0.15, while in the faint source(e.g., SNR < 0.25), this increases with the SBR as follows:p(1%/SBRLE )2 + (1%/SBRHE )2∆αSBR = 2 ×.
ln(Emax /Emin )2.5(2.18)

Comparison in the performance of IACTs

In 1989, the Whipple telescope first succeed the detection of γ-ray emission from the Crab
Nebula (Weekes et al., 1989). Subsequently, HEGRA made the first successful observation
with a stereo system (Daum et al., 1997), while TeV γ-ray observations of the southern sky
were initiated by CANGAROO (Tanimori et al., 1994). The MAGIC, VERITAS, and H.E.S.S
telescopes are currently in operation. IACTs have many upgrades over the 30-year operation.
The comparison in the performance of existing IACTs are summarized in Table 2.5. With theend of the operation of CANGAROO in 2011, only H.E.S.S. is currently installed in the Southern
Hemisphere.
The sensitivity curves of the gamma-ray telescopes are shown in Fig 2.21. MAGIC, H.E.S.S.,
and VERITAS have comparable sensitivity in the TeV band, while the MAGIC can observe thelower energies since the stereo system has telescopes with a larger diameter. Note that thesensitivity of MAGIC depends on the zenith angle, but there is little difference in sensitivity
between the low and medium zenith observations. Furthermore, the next generation project
CTA, which will have sensitivity about one order of magnitude better than the existing IACTs,
is underway.
The highest sensitivity in the GeV band is provided by Fermi -LAT, which directly detects
primary γ rays with a telescope onboard a satellite (the detail will be described in Sect. 5.2). In
the energy band above 10 TeV, LHAASO provides the highest sensitivity. LHAASO observes
γ rays by detecting the secondary particles rather than the Cherenkov radiation from them.

CHAPTER 2. THE MAGIC γ-RAY TELESCOPE

48

(0 ,

HAW

(l, b
)=

10 11

C (5

0)

yr)

10 12

b)
=

(1 2

0,

45

)

hr)

10 14108

(l,

(50

10 13

Fermi/LAT (10 yr)

VERITAS,
HESS (50 hr)

th
ou
A-S
CT

Differential sensitivity ×E2 [erg cm 2 s 1]

10 10

SWGO,
MAGIC at low zenith (50 hr)
MAGIC at medium zenith (50 hr) LHAASO (5 yr)
1014
1010
1012
Energy [eV]HiSCORE(1000 hr)

1016

Figure 2.21: Differential 5σ sensitivity of gamma-ray telescopes. Black open and filledsquares represent the post-upgraded MAGIC sensitivity for 50 hr at low (< 30◦ ) or medium(30◦ –45◦ ) zenith angle (Aleksić et al., 2016), respectively.
Future instruments (CTA:
Bernlöhr et al. (2013), SWGO: Albert et al. (2019), and HiSCORE: Tluczykont et al.
(2014)) are shown as dotted lines and existing instruments (Fermi-LAT: https://www.
slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm, VERITAS: https://veritas.sao.arizona.edu/about-veritas/veritas-specifications, H.E.S.S.: Holler et al.
(2016), HAWC: Abeysekara et al. (2017), and LHAASO: Addazi et al. (2022)) are shown as solidlines.
Although the resultant angular resolution is inferior to IACTs, the system can achieve a large
effective area. While the IACT can work only at night, this type of detector, the so-called air
shower array, can observe gamma rays throughout the day, benefiting the high-duty cycle.

Chapter 3

MAGIC observation in the vicinityof SNR G106.3+2.7
3.1

PeVatron candidate: SNR G106.3+2.7

One of the 100-TeV γ-ray sources, LHAASO J2226+6057, is spatially consistent with SNR
G106.3+2.7 (Cao et al., 2021), and thus this SNR is the most promising PeVatron candidate.
However, previous studies have shown that this region is very complex.
This region consists of the PWN named “Boomerang” as well as the SNR G106.3+2.7, thelatter of which was first discovered by the radio observations (Joncas & Higgs, 1990). As shownin Fig. 3.1 (a), the SNR has a comet-shaped morphology, where the bright circle and fainterextended parts were labeled “head” and “tail”, respectively. The spectral indices of the head andtail regions are α = 0.49±0.05 and 0.70±0.07 (Pineault & Joncas, 2000), respectively, which are
in rough agreement with the typical value (∼ 0.5) in SNRs (Reynolds et al., 2012). Although theorigin of this comet-shaped structure is not well understood, the HI line observations suggested
that it was formed due to the surrounding gas distribution (Kothes et al., 2001). Under theassumption that the SN occurred in the head region, its shock wave can be expected to havegrown in the southeast direction, where the HI density is relatively low (see Fig. 3.1 b). The
Boomerang PWN is located at the north of the head region. This is driven by an energetic pulsar
PSR J2229+6114, whose spin-down luminosity and the characteristic age are 2.2 × 1037 erg s−1
and ∼ 10 kyr, respectively (Halpern et al., 2001a). The spectrum of the PWN shows a spectral
break at 4.3 GHz attributed to synchrotron cooling, suggesting that its age is ∼ 3.9 kyr (Kothes
et al., 2006). Association of HI and molecular clouds with SNR G106.3+2.7 suggests that the
distance is 800 pc (Kothes et al., 2001), while the estimation from X-ray absorption indicatesthat it is 3 kpc (Halpern et al., 2001b).
In this region, non-thermal X-ray emissions have been found not only in the compact regionof the PWN and also in the entire SNR, i.e., both the head and tail regions (Ge et al., 2021;
Fujita et al., 2021). Fujita et al. (2021) claims that the emission in both regions is generated
by electrons originating in the PWN, while Ge et al. (2021) argue that the tail emission is more
likely due to the electrons accelerated in the shock of the SNR.

49

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

(a)

50

(b)

(c)

Figure 3.1: Skymaps in the vicinity of SNR G106.3+2.7. (a) Radio continuum intensity map
at 1420 MHz. (b) HI line intensity map overlaid with the contours of 12 CO (J = 1 − 0) line
(thick black line) as well as radio continuum emission at 1420 MHz (thick white line). Panels(a, b) are taken from Kothes et al. (2001). (c) TeV γ-ray excess event map obtained with
the VERITAS observation (Acciari et al., 2009). The black and magenta contours represent
the radio continuum and CO line emissions, respectively. The yellow open cross represents thelocation of PSR J2229+6114.

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

51

Many observations in γ-ray bands have also been performed for this region. The FermiLAT and EGRET observations have revealed a pulsed GeV emission from the pulsar (Abdoet al., 2009a; Hartman et al., 1999) and a steady emission above 3 GeV (Xin et al., 2019; Fanget al., 2022). The emission region of the latter is close to the SNR tail ( ∼ 0.4◦ away from theposition of PSR J2229+6114) and was better described with a 0.25◦ radius disk than a point-likesource (Xin et al., 2019). The steady emission was also detected with the VERITAS telescopein the TeV band (Acciari et al., 2009) and labeled VER J2227+608 (Fig. 3.1 c). The shapeof the TeV emission region can be characterized with an elongated two-dimensional Gaussianwith 0.27 ± 0.05 (0.18 ± 0.03)◦ extent in the major (minor) axis. In addition, γ rays above
10 TeV from this region were first detected with the Milagro experiment (Abdo et al., 2007,2009b) and also observed with HAWC (Albert et al., 2020), Tibet ASγ (Amenomori et al.,2021b), and LHAASO experiments (Cao et al., 2021). The HAWC and Tibet ASγ results
suggest a power law spectrum without a cutoff and the spectral indices are 2.25 ± 0.23stat and3.17 ± 0.63stat , respectively. This emission above tens of TeV provides a lower limit on the
maximum energy of the particles accelerated in this object. If the emission process is leptonic,
an exponential cutoff energy of the electron must be higher than 270 TeV (Albert et al., 2020) or
190 TeV (Amenomori et al., 2021b), while if it is hadronic, the maximum proton energy should
be higher than 800 TeV (Albert et al., 2020) or 500 TeV (Amenomori et al., 2021b).
However, it is still inclusive whether parent particles are accelerated in the SNR blast wave
or the PWN complex due to the limited angular resolution of the γ-ray observations. The spatial
coincidence between the molecular clouds and the γ-ray emission supports a hadronic origin of γrays (e.g., Acciari et al., 2009), while the leptonic scenario for the origin is not ruled out (Fujita
et al., 2021). Therefore, further precise γ-ray observations with a better angular resolution thanthat of previous studies are needed to resolve its origin.

3.2

Observation and data analysis

For this purpose, we observed SNR G106.3+2.7 using the MAGIC telescopes between May 2017
and August 2019 for 183.7 hr. To estimate the background simultaneously, all observations were
performed in wobble mode (see Sect. 2.3.3) at three positions (RA = 336.31◦ , DEC = 61.40◦ ;RA = 338.25◦ , DEC = 61.06◦ ; RA = 336.66◦ , DEC = 60.42◦ ) with an offset of 0.57◦ from theposition (RA = 337.05◦ , DEC = 60.96◦ ), which is close to the center of VER J2227+608 (RA= 337.0◦ , DEC = 60.8◦ ). The left panel of Fig. 3.2 shows the zenith angle distribution during
the observations. For the event reconstruction, we also use hadron (real observations) and γray (MC simulation) data at approximately the same zenith angle as the SNR observations, as
shown on the left of Fig. 3.2. We estimate the analysis energy threshold with the MC simulation
data by assuming the energy distribution follows a simple power-law function with an index of
2.0 (the right panel in Fig. 3.2). Since the peak of this histogram can be defined as the analysis
energy threshold for IACT observations, we conservatively set to 200 GeV for this work.
The MAGIC angular resolution, characterized by the point spread function (PSF), for thisanalysis was estimated to be 0.084◦ (68% containment radius) at E > 0.2 TeV and 0.072◦ at

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

52

108
22

107

20
18
Number of events (× 10-3)

Number of events

106
105
104
3

10

SNR G106.3+2.7

102

Hadron1

14
12

Entries
2011230
Mean
181.3
RMS
65.12
Gauss Mean
159.05Gauss Sigma
90.73

8
6

2
0
0

10

Threshold

104

MCs data

10

16

20

30
40
Zenith angle [deg]

50

60

50

100

150
Energy [GeV]

200

250

300

Figure 3.2: Left: Zenith angle distribution in the MAGIC observations for SNR G106.3+2.7(red). Also shown the blue and magenta histograms represent the hadron (real observations)and the γ-ray (MC simulation) data. Right: Energy distributions of the MC simulation data,
the peak energy of which corresponds the analysis threshold.
Table 3.1: Comparison of the conditions in the observations for SNR G106.3+2.7 between
MAGIC (this work) and VERITAS (Acciari et al., 2009)
VERITAS (Acciari et al., 2009) MAGIC (This work)
Observation period
2009
2017–2019
Observation time
33 hr
122 hrEnergy threshold
0.63 TeV
0.20 TeV◦ (> 0.2 TeV)
Angular resolution
0.084
0.11◦
(68%-containment radius)
0.070◦ (> 1.0 TeV)
E > 1 TeV, which is the best angular resolution among the previous γ-ray observations for thisobject (e.g. 68% containment radius of the observation with the VERITAS telescope performed
in 2009 is 0.11◦ ).
The data analysis was performed with the MAGIC standard analysis package (Zanin, 2013).
The data selection was based mainly on the transmission of the atmosphere monitored with a
LIDAR system (Fruck et al., 2014). In this analysis we only selected data with an atmospherictransmission above 85% of the optimum. After quality cuts, the total dead time corrected
observation time is 121.7 hours. We used the wobble-map method (e.g., Vovk et al., 2018) for
estimating backgrounds. To cross-check the results obtained with the MAGIC standard analysis
package, we used the SkyPrism package (Vovk et al., 2018), which includes independent methods
to compute the instrument response functions and estimate the energy spectra using a spatial,
maximum likelihood fit. Both results are in good agreement, as we will mention in Section 3.4.2.
The comparison of the observation conditions between this work and the previous study with
the VERITAS telescope (Acciari et al., 2009) is summarized in Table 3.1.

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

3.3

Results

3.3.1

Morphological study

53

Fig. 3.3 shows the pre-trial significance maps around SNR G106.3+2.7 in different energy bands.
Panel (a) shows the morphology of the γ-ray emission above 0.2 TeV overlaid with the radio
emission contours at 408 MHz measured by DRAO (Pineault & Joncas, 2000) and 12 CO (J =1 − 0) emission contours (Taylor et al., 2003). γ-ray emission above 0.2 TeV from the directionof VER 2227+608 (shown in Fig. 3.1 c and Fig. 3.3 c) is clearly detected. We evaluate thestatistical significance using Eq. 17 of Li & Ma (1983):
s 
S=


Non
1+α
Noff
2 Non ln
+ Noff ln(1 + α), ,α Non + Noff
Non + Noff

(3.1)

where Non (Noff ) is the number of signal (background) events and α is a ratio between theobservation time of the signal region and background regions1 . Integrating the same area as
VERITAS and using Eq. 3.1, the resultant statistical significance is found to be 8.9σ. It is
extended and spatially coincident with the radio shell of the SNR, i.e., the emission region isextending from the SNR head region to the tail region. The emission at the tail coincides withstrong 12 CO (J = 1 − 0) emission, but the overall emission profile does not follow well the CO
distribution. The emission at the head is in fact seen where 12 CO (J = 1 − 0) emission is notobserved. It should be noted that 12 CO (J = 1 − 0) does not trace all existing interstellar gas
as will be discussed in Sect. 4.2.
The panels (b), (c) and (d) of Fig. 3.3 show the maps at 0.2 to 1.1 TeV, 1.1 to 6.0 TeV and
6.0 to 30 TeV, respectively. The morphology of the detected γ-ray emission clearly changes withenergy. By fitting with a symmetric Gaussian function, the center position of the γ-ray emission
in the highest energy band of 6.0–30 TeV is estimated to be (RA, DEC) = (336.66 ± 0.05◦ ,+60.87 ± 0.02◦ ) (J2000), which is offset from the location of PSR J2229+6114 by 0.47 ± 0.03◦(Panel d). On the other hand, the lower energy emission extends close to the pulsar position(Panels b and c). The centroid of the low energy emission for 0.2–1.1 TeV and its distance fromthe pulsar position are found to be (RA, DEC) = (336.99 ± 0.04◦ , +61.04 ± 0.02◦ ) (J2000) and
0.24 ± 0.03◦ . The 1σ extension at 6.0–30 TeV, after removing the effect of PSF is 0.14 ± 0.09◦ ,which is consistent with the value (0.24 ± 0.14◦ ) reported by Tibet ASγ Amenomori et al.
(2021b).
The spatial distribution in Fig. 3.3 (a) appears to have a more complex shape than the
double symmetric Gaussian function, but current statistics allow fitting data with this function.
The best-fit parameters are as in Table 3.2 for the center position, and 0.083◦ (0.087◦ ) for the 1σ
extension of head (tail). Fig. 3.4 shows the residual map after subtracting two Gaussian sources
and its significance distribution of the residuals. The distribution is consistent with the null
hypothesis, which indicates that, with the current statistics, the double Gaussian assumption isvalid (i.e., the reduced chi-squared statistics of the fit is reasonably smaller than 1), though the1

The evaluation with this equation is valid even if the events is not followed to the Gaussian distribution
and/or if the observation times in the signal and background regions are different.

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

54

Figure 3.3: Energy-dependent pre-trial significance maps of SNR G106.3+2.7 observed with
the MAGIC telescopes. (a) The map above 0.2 TeV. The white circle labeled ”PSF” represents
the 0.075◦ size of a Gaussian kernel (corresponding to the MAGIC γ-ray point spread function)
for this analysis. The position of PSR J2229.0+6114 is marked with the open yellow cross. The
cyan contours show the radio emission of SNR G106.3+2.7 at 408 MHz by DRAO (Pineault &Joncas, 2000). The green contours represent 12 CO (J = 1 − 0) line intensity integrated overthe velocity range from −6.41 to −3.94 km s−1 . The white dotted circles show θ2 cut regionsof the head and tail regions, respectively, as shown in Table 3.2. Also shown by white squares
are the pointing positions used in the observations. (b) The map at 0.2–1.1 TeV. The whitecircle labeled ”PSF” represents the 0.100◦ size of a Gaussian kernel as the panel (a). The yellow
solid and dotted circle represent the extension and location of the Fermi -LAT source (Xin et al.,2019) and the analysis region for the head region used in Liu et al. (2020), respectively. (c) The
map at 1.1-6.0 TeV. The white circle labeled ”PSF” represents the 0.065◦ size of a Gaussian
kernel as the panel (b). The green ellipse and dotted circle represent the extended TeV γ-ray
emission of VER J2227+608 and θ2 cut region used in the VERITAS paper (Acciari et al., 2009),
respectively. (d) The map at 6.0 - 30 TeV. The white circle labeled ”PSF” represents the 0.065◦
size of a Gaussian kernel as the panel (c). The magenta solid and dotted circle represent theextended γ-ray emission above 10 TeV observed with Tibet ASγ (Amenomori et al., 2021b) andthe upper limit at 90% confidence level of the Gaussian extension of HAWC J2227+610 (Albert
et al., 2020). The open square and diamond show the centroid of the VHE γ-ray emission
detected with Milagro (Abdo et al., 2009b) and LHAASO (Cao et al., 2021), respectively.

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

55

Number of entries per bin

Data
2

10

Null hypothesis

10

1

10−1

−4

−2

0

2
4
Significance [σ]

6

8

10

Figure 3.4: Top: The residual map after subtracting two Gaussian sources in the energy range
above 0.2 TeV. As for the Gaussian parameters, the locations are the same as the definition ofhead and tail and each 1σ extension radius is 0.085◦ . Bottom: Pre-trial significance distributionof the residual map.
Table 3.2: Regions considered in this work for the analysis of MAGIC data and their modelling.
Sourcehead region
tail regionRA337.◦ 13336.◦ 72

DEC
61.◦ 1060.◦ 84

Radius0.◦ 160.◦ 16

true γ-ray source morphology may be more complex.
To understand the emission mechanism better, we study the γ-ray spectra at the head andthe tail regions. The parameters for the two regions are defined in Table 3.2 and also shown in
Fig. 3.3 (a). The centers of these regions are adapted from a fit to the γ-ray map above 0.2 TeV
with a double symmetric Gaussian. The position of the tail emission is in good agreement withthe peak position observed with VERITAS/Tibet (Acciari et al., 2009; Amenomori et al., 2021b)
and included within the upper limit at 90% confidence level of the Gaussian extension of HAWC
J2227+610 (Albert et al., 2020). The radii of these areas are chosen to be the same for bothregions and of maximum length without the regions overlapping.
Fig. 3.5 shows the so-called θ2 distributions of the two regions, where θ is the opening angle
between the center of the region and the event arrival direction. The excesses are detectedfrom the head and tail regions above 0.2 TeV with statistical significance of 6.2σ and 6.9σ,
respectively, evaluated using Eq. 17 of Li & Ma (1983). Fig. 3.6 shows the energy-dependent θ2plots, in which the energy range is the same as in Fig. 3.3. The significances for 0.2–1.1 TeV are
4.8σ at head and 2.8σ at tail, while for 6.0–30 TeV they are 6.5σ at the tail, and only 2.4σ atthe head, indicating that the magnitude ratio of the head and the tail emissions flips between
the low and high energy bands.

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7
1300

1300
Significance (Li&Ma) = 6.2σ

1200

Significance (Li&Ma) = 6.9σ

1200

ON events

1100

ON events

1000

OFF events

1000

OFF events

Nevents

1100
Nevents

56

900
800

900
800

Head

700
600
0

0.1

Tail

700
600
0

0.2

0.1

θ2 [ deg2 ]0.2
θ2 [ deg2 ]Figure 3.5: θ2 distributions of ON (blue circles) and OFF (black line) events above 0.2 TeV
toward the center of the head region (left) and that of the tail region (right). The region
between zero and the vertical dashed line (at θ2 = 0.0256 deg2 ) has been used to estimate ON
and OFF events. The OFF data represents the average of six regions rotated by 120 and 240 deg
with respect to each wobble center from the ON region.

3.3.2

Energy spectrum

Fig. 3.7 and 3.8 show the γ-ray spectra of the two regions defined in Table 3.2 and the extraction region of VER J2227+608 (Acciari et al., 2009), respectively. Using the forward-folding
method (Aleksić et al., 2016), the spectra are fitted with a power-law function:
dN
= N0
dE



E
3 TeV

−Γ
.

(3.2)

The best-fit parameters are summarized in Table 3.3. The γ-ray spectrum in the tail region
has a higher flux and a marginally harder index than that of the head region. For the VER
J2227+608, using the same integration region as VERITAS, our results are consistent with theirs(Acciari et al., 2009) within the statistical uncertainties in both the index and the normalization
at 3 TeV. The apparent discrepancy seen in Fig. 3.8 between the MAGIC results and the
Tibet ASγ measurement at the 6–20 TeV range, amounts to only 1.4σ statistical uncertainty.
Considering the source extension of VER J2227+608 and the MAGIC PSF, the flux derived inthis work may correspond to ∼ 60% of the whole region estimated with the other experiments.
If this loss is considered, the discrepancy between MAGIC and Tibet ASγ relaxes from 1.4σ
to 1.1σ. In addition, if the systematic uncertainties are taken into account, both results agree
within 1σ.

3.4

Systematic uncertainties in the extended source analysis

3.4.1

Estimation of the effective area

The effective area used in the spectral analysis is estimated by assuming a source radius. In
Sect. 3.3.2, we assumed 0.16◦ as the radii of head and tail. However, as mentioned in Sect. 3.3.1,
the actual radii and morphology of the γ-ray emissions are unknown. Here, we evaluate theCHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7
6500

57

6500
6000

ON events

5500

OFF events

5500

OFF events

Nevents

ON events

Nevents

6000

5000

5000

4500
4000
3500
0

4500

Head at 0.2 - 1.1 TeV: 4.8σ
0.1

4000
3500
0

0.2

Tail at 0.2 - 1.1 TeV: 2.8σ
0.1

θ2 [ deg2 ]
300

300

250

ON events

250

ON events

200

OFF events

200

OFF events

Nevents

Nevents

0.2
θ2 [ deg2 ]

150
100
50
0
0

150
100
50

Head at 1.1 - 6.0 TeV: 7.3σ
0.1

0
0

0.2

Tail at 1.1 - 6.0 TeV: 7.2σ
0.1

θ2 [ deg2 ]0.2
θ2 [ deg2 ]

70

70

60

60

ON events

50

ON events

50
OFF events
Nevents

Nevents

OFF events
40
30

40
30

20

20

10

10

0
−10
0

Head at 6.0 - 30 TeV: 2.4σ
0.1

0.2

0
−10
0

θ2 [ deg2 ]Tail at 6.0 - 30 TeV: 6.5σ
0.1

0.2
θ2 [ deg2 ]Figure 3.6: The same as in Fig. 3.5 but for the energy of 0.2–1.1 TeV (top), 1.1–6.0 TeV (middle),
and 6.0–30 TeV (bottom).
Table 3.3: Comparison of the spectral parameters between this MAGIC results reported here andthe VERITAS ones (Acciari et al., 2009). All sources were fitted with the power-law function of
Eq. 3.2, using a forward-folding method (Aleksić et al., 2016). The systematic uncertainties are
estimated based on as in Sect. 2.4.

SourceN0 at 3 TeV−14(10
cm−2 s−1 TeV−1 )Γχ2 /ndf

head
tail
VER J2227+608 (MAGIC)VER J2227+608 (VERITAS†)

3.8 ± 0.7stat ± 0.7sys
6.0 ± 0.7stat ± 1.0sys
13.1 ± 1.1stat ± 2.1sys
11.5 ± 2.7stat ± 3.5sys

2.12 ± 0.12stat ± 0.15sys
1.83 ± 0.10stat ± 0.15sys
1.91 ± 0.07stat ± 0.15sys
2.3 ± 0.33stat ± 0.30sys

5.5/6
2.6/6
7.1/6
-

† Ref: Acciari et al. (2009).

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

58

E2 dN/dE [TeV cm-2 s-1]

MAGIC head
MAGIC tail
−12

10

10−13
1011

1012
Energy [eV]1013

Figure 3.7: Energy spectra of the head and tail regions. Red and blue data represent the spectraof the head and tail, respectively. The color bow-tie areas show the result of fitting with a simple
power-law function and 1σ statistical uncertainties.
Table 3.4: Details of the Crab nebula data.
Observation period
Oct. 2016 – Mar. 2018
Observation time
23.0 hrZenith angle
30–50◦
Wobble angle
35◦ , 215◦
Offset angle
0.4◦
Analysis energy threshold
100 GeV
impact of different assumptions in the source model on estimating the effective area, as anadditional systematic uncertainty. We investigate the effects of (I) assumed source extension and(II) assumed center position, using the observation data of the Crab Nebula, whose morphology
is well known. The data set used in this study is summarized in Table 3.4.
(I) Uncertainty of the source extension
We examine the difference of the computed effective area and the energy spectrum due to theassumed radius of the γ-ray source (see Fig. 3.9 a for a schematic view). For this, we prepare
the MC simulation data shown in Figure 3.10, in which the radii are set 0.05, 0.1, 0.15, and 0.2◦ .
Figure 3.11 shows the dependence of the effective area on the assumed source radius. The
ratio to the point source data is also shown in the lower part of Fig. 3.11. The resultant
discrepancy between the assumptions is found to be at most 5% at the energy band higher
than the analysis threshold of 100 GeV. The energy spectrum of Crab Nebula is fitted with aCHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

59

MAGIC VER J2227+608 (this work)
MAGIC head (this work)

10−11

2

-2

E dN/dE [TeV cm s-1]

MAGIC tail (this work)

10−12

Fermi-LAT (2019)
VERITAS (2009)
Milagro (2007, 2009)
HAWC (2020)
Tibet ASγ (2021)
LHAASO (2021)

10−13

10−14
1010

1011

1013

1012
Energy [eV]

1014

1015

Figure 3.8: Spectral energy distribution of the whole region of SNR G106.3+2.7. Green data
represents the spectrum of VER J2227+608 region as measured with the MAGIC telescopes.
Also shown as the shaded blue and red regions are the same as in Fig. 3.7. The open triangles,
open stars, open squares, x marks, and open diamonds show the Fermi -LAT (Xin et al., 2019),
VERITAS (Acciari et al., 2009), Milagro Abdo et al. (2007, 2009b), Tibet ASγ (Amenomori
et al., 2021b), and LHAASO measurements (Cao et al., 2021), respectively. The black bow-tie
area shows a power-law fit and 1σ statistical errors measured by HAWC (Albert et al., 2020).

(I) Assumption of the extension(II) Assumption of the position

CameraCenterCameraCenter0

Assumed Source radius

0.1

0.2 0.3 0.4 0.5 0.6

Distance from Camera Center [deg]0.70Assumed Source position

0.1

0.2 0.3 0.4 0.5 0.6

0.7

Distance from Camera Center [deg]Figure 3.9: Schematic diagram in the calculation of systematic uncertainties due to the assumption of (I) the extension and (II) the center position of a γ-ray source.

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

1

1

70

(a) Point source(c) R = 0.10 deg

0.8
3

0.6

0.6
50

0.2

40

deg

deg

0.2
0

30

−0.2
−0.4

2.5

0.4

20

−0.6

0
1.5

−0.2
−0.4

−0.8
0

0.2

0.4

0.6

0.8

1

0.5

−0.8
−1
−1 −0.8 −0.6 −0.4 −0.2

0

deg

0

0.2

0.4

0.6

0.8

1

3.5

0.4

3
2.5

0

2

−0.2
1.5

−0.4

1

−0.6
10

0.6

0.2

2

deg

0.4

4

0.8

60

−1
−1 −0.8 −0.6 −0.4 −0.2

13.5(b) R = 0.05 deg

0.8

60

−0.6

1

−0.8

0.5

−1
−1 −0.8 −0.6 −0.4 −0.2

0

deg

0

0.2

0.4

0.6

0.8

1

0

deg

1

1

(d) R = 0.15 deg(e) R = 0.20 deg

0.8

6

0.8

5

0.6

0.6
5
4

0.4

0.4
4

0.2
3

0

−0.2

deg

deg

0.2

0
3

−0.2
2

−0.4

−0.4

−0.6

−0.6

1

−0.8
−1
−1 −0.8 −0.6 −0.4 −0.2

2

1
−0.8

0

deg

0.2

0.4

0.6

0.8

1

0

−1
−1 −0.8 −0.6 −0.4 −0.2

0

0.2

0.4

0.6

0.8

1

0

deg

Figure 3.10: Distributions of MC simulation data used for the estimation of the systematic
uncertainties (in the camera coordinate). Panel a shows the data distribution for a point source,
while Panels b–e show those for extended sources with a radius of 0.05–0.2◦ .

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

61

106

2

Effective area [m ]

105
104

Point source
R = 0.05 deg
R = 0.10 deg
R = 0.15 deg
R = 0.20 deg

103
102
10
1

A/APoint

1.4

102

103

104

102

103
Energy [GeV]

104

1
0.6

Figure 3.11: The effective areas as a function of an estimated energy. The black data show theresult of the point source MC data. The blue, green, magenta, and red data show the result of
assuming the source radius of 0.05◦ , 0.10◦ , 0.15◦ , and 0.20◦ , respectively. The lower panel shows
the ratio to that of the point source model. The gray shaded area indicates the lower energy
band than the analysis threshold.
single-power law function as follows:
dF
= N0
dE



E
1 TeV

−Γ
.

(3.3)

The best-fit parameters for the case of using the point source MC are N0 = (2.78 ± 0.04stat ) ×
10−11 cm−2 s−1 TeV−1 and Γ = 2.48 ± 0.01stat , respectively. We then compare the fitting parameters with those derived using the other MC data and show the result in Fig. 3.12. The
systematic differences are at most 3% for the normalization flux and 1% for the spectral index,
respectively.
(II) Uncertainty of the source location
Next, we investigate the effect due to the assumed source position, as shown in the panel b of
Fig. 3.9. While the wobble offset (D) of 0.4◦ is correct, we prepare the MC data with the wobble
offset from 0.1 to 0.7◦ shown in Fig. 3.13. Note that the source radius is assumed to be 0.05◦ .
Fig. 3.11 and Fig. 3.12 show the effective areas computed with MC data of Fig. 3.13 andthe comparison of the spectral fitting results, respectively. The difference from the case with

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

1.1
1.08

1.1

(a) Normalization Flux

1.08

1.04

1.04

1.02

1.02

Γ/ ΓPoint

1.06

N0/N0, Point

1.06

1

0.98

0.96

0.96

0.94

0.94

0.92

0.92
0

0.05
0.1
0.15
Assumed source radius [deg](b) Photon index

1

0.98

0.9

62

0.9

0.2

0

0.05
0.1
0.15
Assumed source radius [deg]0.2

Figure 3.12: Assumed source radius dependence of normalization flux (a) and spectral index(b). A source radius of 0◦ corresponds to the point source of MC simulation data.

1

1

3

0.6

18

0.4

0.4

16

0.2

14

−0.2

2

−0.2

−0.4

1.5 −0.4

−0.6

1

−0.8

0.5 −0.8−1−1 −0.8 −0.6 −0.4 −0.2

0

0.2

0.4

0.6

0.8

1

8

0.2

2.5

0

20

0.6

deg

deg

0.2

0.8
10

3.5

0.4

6

0

−0.6

2

10
8

−0.4

6

−0.6

4

−0.8

−1
−1 −0.8 −0.6 −0.4 −0.2

0

12

0

−0.2
4

deg

0

0.2

0.4

0.6

0.8

1

2

−1
−1 −0.8 −0.6 −0.4 −0.2

0

deg

1

0.8

0.6

0.6
40

0.4

0.2

0.4

0.6

0.8

1

0

1

(e) D = 0.6 deg
50

0

deg

1

(d) D = 0.5 deg
0.822(c) D = 0.3 deg

0.8

4

0.6

1

(b) D = 0.2 deg

4.5

deg

(a) D = 0.1 deg
0.8(f) D = 0.7 deg

70

0.8

90

60

0.6

80

0.4

70

0.2

60

0.4
50

30
0

−0.2

20

−0.4

−0.2
−0.4

−0.6

10

−0.8
−1
−1 −0.8 −0.6 −0.4 −0.2

0

deg

0.2

0.4

0.6

0.8

1

0

0

50

30

−0.2

40

−0.4

30

−0.6

20

20

−0.6
10

−0.8
0

40

deg

0.2

deg

deg

0.2

−1
−1 −0.8 −0.6 −0.4 −0.2

0

deg

0.2

0.4

0.6

0.8

1

0

−0.8
−1
−1 −0.8 −0.6 −0.4 −0.2

10
0

0.2

0.4

0.6

0.8

1

0

deg

Figure 3.13: Distributions of MC simulation data, which is used to investigate the systematic
uncertainties due to the assumption of the source position (in the camera coordinate). In Panels
a–f, the assumed distance from the camera center is changed from D = 0.1◦ to 0.7◦ in 0.1◦
increments.

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

63

106

2

Effective area [m ]105104
D = 0.1 deg
D = 0.2 deg
D = 0.3 deg
D = 0.4 deg
D = 0.5 deg
D = 0.6 deg
D = 0.7 deg

103
102
10

A/AD = 0.4 deg

1
1.4

102

103

104

102

103
Energy [GeV]

104

1
0.6

Figure 3.14: Same as Fig. 3.11, but the results with the other MC data of Fig. 3.13.
the correct wobble offset (D = 0.4◦ ) appears to increase in correlation with the absolute valueof D − 0.4◦ . Even if the assumed center position is wrong by 0.2◦ (which is larger than PSF),
the systematic discrepancies in the fitting parameters (flux and index) are still less than 1%.
(III) Dependency on the offset angle
As mentioned in Sect. 2.3.3, the camera sensitivity depends on the wobble offset (Aleksić et al.,2016). This may affect the investigation of this systematic uncertainty on the effective area. In
particular, the distance between the assumed source position and the pointing position in the
SNR G106.3+2.7 observation is 0.42–0.78◦ , which is different from the Crab Nebula observation(D = 0.4◦ ). Fig. 3.16 shows the offset angle dependence of the effective area at 1 TeV. The
computed effective area decreases as the offset angle increases, and its slope steepens after 0.55◦ .
The systematic uncertainty of the effective area would increase following the change in slope,i.e., a factor of ∼ 1.9.
The total effect due the uncertainty of the source assumption can be evaluated by:(Systematic uncertainty due to the MC source model)q= (Offset angle factor) × (extension uncertainty)2 + (position uncertainty)2 . (3.4)

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

1.1
1.08

1.1

(a) Normalization flux

1.08
1.06

1.04

1.04

1.02

1.02

Γ/ ΓPoint

1.06

N0/N0, Point

64

1

1

0.98

0.98

0.96

0.96

0.94

0.94

0.92

0.92

0.9
0

0.1

0.2
0.3
0.4
0.5
0.6
0.7
Assumed distance from Camera Center [deg](b) Photon index

0.9
0

0.8

0.1

0.2
0.3
0.4
0.5
0.6
0.7
Assumed distance from Camera Center [deg]Figure 3.15: Assumed source location dependence of the fitting parameters, normalization flux(a) and spectral index (b). Also shown with the red line is the fit result with a quadratic
function.

3

×10
170

Slope

0.30-0.55 deg= -5.7 × 104

Effective area at 1 TeV [m2]160
1500.55-1.00 deg

140: -1.1 × 10

5

Slope

130
120
110
100

Slope

/Slope0.30-0.55 deg

0.3

0.4

0.5

0.55-1.00 deg

~ 1.9

0.6
0.7
Offset angle [deg]

0.8

0.9

1

Figure 3.16: Effective area at 1 TeV as a function of the offset angle. The solid lines show thefit result with a linear function to the data of 0.3–0.55◦ (red) and 0.55–1.0◦ (blue).

0.8

-2

10−13

10−12

10−13

Tail
3

102

10
Energy [GeV]

104

102

103
Energy [GeV]104

ratio to SkyPrism

ratio to SkyPrism

Head
2
1.5
1
0.5
0

65

Forward-unfolding
Tikhonov
SkyPrism2

10−12

E dN/dE [TeV cm s-1]

-2

Forward-unfolding
Tikhonov
SkyPrism

2

E dN/dE [TeV cm s-1]

CHAPTER 3. MAGIC OBSERVATION IN THE VICINITY OF SNR G106.3+2.7

2
1.5
1
0.5
0

3

102

10
Energy [GeV]

104

102

103
Energy [GeV]

104

Figure 3.17: Comparison in energy spectra between the results of Sect. 3.3.2 and the SkyPrism
analysis.
Hence, the effects on the normalization and slope are estimated as:
σnorm ≲ 6.0%

(3.5)

σindex ≲ 2.7%,(3.6)

respectively. As a result, the systematic uncertainty due to differences in the models used in
the extended source analysis are found to be smaller than the systematic uncertainties in thestandard MAGIC analysis (15–20%) and thus be negligible.

3.4.2

Effect due to θ2 cut

Although we cannot claim the proper source shape of the head and tail components from thepresent statistics, under the assumption that the source has Gaussian-like extension with 1σ of
0.085◦ (after removing the effect of PSF), the loss and contamination rate from the θ2 cut are
estimated to be 23.5% and 2.7%, respectively.
Another concern is leakage from the θ2 cut region due to the limitation of angular resolution.
In the SkyPrism analysis (Vovk et al., 2018), we can calculate the effective area for a given spatial
source model, considering the energy dependence of PSF. We use the disk model with a radiusof 0.16◦ as a spatial model. The comparison in energy spectra is shown in Fig. 3.17. The spectra
derived with the SkyPrism method are in good agreement with the result derived in Sect. 3.3,
and the residual of each energy bin is lower than the statistical 1σ uncertainties.

Chapter 4

Spectral modeling and implication tothe γ-ray origin of SNR G106.3+2.7
The spatial coincidence of the MAGIC γ-ray emission and the 408 MHz radio continuum shown in
Fig. 3.3 (a) indicates that the VHE γ-ray emission is associated with the radio SNR G106.3+2.7.
On the other hand, the emission above 6 TeV is detected only in the tail region but not in thehead region, suggesting that the γ-ray origin is not the same for both regions, as well as, theemission up to 100 TeV observed in previous studies (Abdo et al., 2009b; Albert et al., 2020;
Amenomori et al., 2021b; Cao et al., 2021) comes mainly from the tail region.
Here, we attempt to model the γ-ray emission mechanism of the head and the tail region
individually, while the previous studies (e.g., Liu et al., 2020; Ge et al., 2021; Bao & Chen, 2021)
discussed the origin of γ rays using the spectrum up to 100 TeV of the whole region of this object.
The possible mechanisms of γ-ray radiation are electron bremsstrahlung, inverse Compton (IC)scattering, and hadronic emission, as described in Chapter 1. Considering the three cases where
each emission process dominates, we search for model parameters that reproduce the observed
data in each case and try to identify the origin of the head and tail emission.

4.1

Modeling

As mentioned above, we consider three cases: (i) bremsstrahlung-dominated, (ii) IC-dominatedmodel, and (iii) hadronic-dominated model. The electron synchrotron emitting in the radioto-X-ray band is also considered. The model emissions are computed with the naima framework (Zabalza, 2015).
The radio data for the head and tail regions are taken from Pineault & Joncas (2000). The
X-ray spectra for the head and tail regions are taken from results of the ”East” and ”West”
regions from Fujita et al. (2021), respectively, multiplying the intensity by the area of a circle
with a radius of 0.16◦ used in the MAGIC analysis. In this study, we do not consider the directcontributions from the compact Boomerang nebula, whose angular diameter is ∼ 0.05◦ , because
the γ-ray flux of the region is estimated to be ∼ 10% or less of the head region from the radio
and X-ray flux (Liu et al., 2020). At GeV range, Xin et al. (2019) and Liu et al. (2020) reported
66

CHAPTER 4. SPECTRAL MODELING AND IMPLICATION TO THE γ-RAY ORIGIN OF
SNR G106.3+2.7
67
the spectral points and upper limits assuming that the sources have a disk shape. They obtained
the radii of 0.20◦ and 0.25◦ for the disks. We scaled down their measurements by (0.16/0.20)2
for the head and (0.16/0.25)2 for the tail. We also assume that the emission above 10 TeV
measured with the air shower experiments (Abdo et al., 2009b; Albert et al., 2020; Amenomori
et al., 2021b; Cao et al., 2021) is only from the tail region, since the extracted spectra, shownin Fig. 3.8, suggest that the head contribution to the total flux above 10 TeV is less than 37.1%(2σ upper limit).
The chi-squared statistic for each model is calculated to confirm that model is statistically
acceptable. Since the results in the X-ray band (Fujita et al., 2021) and HAWC (Albert et al.,2020) have been given as a fitted power-law function, the flux and statistical uncertainty onlyat the normalization energy of the fit are considered in the calculation of chi-squared statistic.
Several experiments have given measurement results in the γ-ray band, each of which differs
beyond statistical error, and thus the influence of systematic error cannot be ignored. This
calculation also takes into account systematic as well as statistical errors for γ-ray data. As thesystematic errors for the Fermi-LAT and LHAASO results were not estimated in the previous
papers for this source, we estimate those of Fermi-LAT and LHAASO with the uncertainties ofthe effective area1 and absolute energy scale (Aharonian et al., 2021), respectively.
In this modeling, the distance to the SNR G106.3+2.7 from the Earth is assumed to be
0.8 kpc, as estimated by the radio observations (Kothes et al., 2001), rather than 3 kpc, as
estimated by the X-ray observations (Halpern et al., 2001b). We note that this difference inassumptions slightly affects the estimations of the SNR radius and the total energy of CRs by
a factor of (3/0.8 kpc) ∼ 4 and (3/0.8 kpc)2 ∼ 14, respectively, but not critically.

4.1.1

Bremsstrahlung-dominated model

The energy spectra of electrons are assumed to follow a power-law function with an exponential
cutoff. The target gas density of each region is estimated using the radio line data of HI and
12 CO (J = 1 − 0) (see Appendix A). As a result, we adopted n−3 for the
HI + nCO ∼ 100 cm
head region and nHI + nCO ∼ 200 cm−3 for the tail region. The following procedure obtained
the model parameters: the total amount of electrons is determined to reproduce the γ-ray datawith the given target gas density, and the magnetic field strength and electron cutoff energy are
determined such that the synchrotron reproduces the radio and X-ray data, respectively.
The top panels of Fig. 4.1 show the modeling results for the head (left) and tail (right).
Parameters for the modeling are summarized in Table 4.1. For both regions, models that can
reproduce X-ray data fails to reproduce γ-ray data. While the electron cutoff energy of ∼300 TeV is required to explain the γ-ray data with the bremsstrahlung emission, the synchrotron
spectrum computed with the cutoff parameter is in contradiction with the observed X-ray flux.
1

https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/Aeff_Systematics.html

CHAPTER 4. SPECTRAL MODELING AND IMPLICATION TO THE γ-RAY ORIGIN OF
SNR G106.3+2.7
68

4.1.2

IC-dominated model

As for the seed photon fields in the IC process, the cosmic microwave background, a galactic nearinfrared (NIR) radiation field, and a galactic far-infrared (FIR) radiation field are considered.
Using the model in the GALPROP package (Porter et al., 2008), the energy density of NIR
and FIR are estimated to be 0.1 eV cm−3 at T = 30 K and 0.3 eV cm−3 at T = 3000 K,
respectively. We do not consider the synchrotron self-Compton emission because its flux isabout two orders of magnitude lower than the IC emissions with the other seed photons. The
electron distribution and the magnetic-field strength are determined in the same manner as in
the bremsstrahlung-dominated model.
The modeling results is shown in the middle panels of Fig. 4.1 and the parameters aretabulated in Table 4.1. The broad-band spectrum of the head region can be explained well withthe leptonic model (χ2 /ndf = 5.0/7). In the case of the tail region, the leptonic model can
reproduce the observed data only in the radio, X-ray, Fermi -LAT and MAGIC band (χ2 /ndf
= 9.3/13), but fails when including air-shower experiments (χ2 /ndf = 103.1/31). To explain
the γ-ray emission above 10 TeV measured by air shower experiments, a high cutoff energy ofelectrons of ∼1200 TeV is required. However, the synchrotron spectrum produced with suchhigh cutoff energy is excluded by the observed X-ray flux. The χ2 /ndf for the model with thehigh cutoff energy is found to be ≫ 1 when considering the X-ray data.

4.1.3

Hadronic-dominated model

For the hadronic model, the γ-ray emission results from the decay of neutral pions produced
by inelastic pp-collisions. The energy spectra of protons are assumed to follow a power-law
function with an exponential cutoff. The target gas density of each region is set to 100 cm−3 forthe head region and 200 cm−3 for the head region, as in the bremsstrahlung-dominated model.
The proton spectrum (flux and energy cutoff) is determined to reproduce the γ-ray data, while
the electron spectrum is given such that the synchrotron radiation reproduces the radio and
X-ray data assuming a magnetic-field strength of 10 µG.
The modeling results show that the γ-ray spectra of both the head and the tail region canbe reproduced assuming a proton maximum energy of 35 TeV and 1 PeV, respectively (χ2 /ndf
= 5.3/7 and 39.9/31). The parent electron distribution should follow a power-law spectrum
different from the ones of protons (parameters shown in Table 4.1).

CHAPTER 4. SPECTRAL MODELING AND IMPLICATION TO THE γ-RAY ORIGIN OF
SNR G106.3+2.7
69

10 8
10 9
10 10
10 11
10 12
10 13
10 14

Bremsstrahlung-dominated model for tailBremsstrahlung (Ec = 2.0 TeV)
Synchrotron (Ec = 2.0 TeV)
Bremsstrahlung (Ec = 300 TeV)
Synchrotron (Ec = 300 TeV)

E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

Bremsstrahlung-dominated model for head

10 4 10 1 102 105 108 1011 1014
Energy [eV]

10 8
10 9
10 10
10 11
10 12
10 13
10 14

10 11
10 12
10 13
10 14

10 11

10 13
10 14

10 4 10 1 102 105 108 1011 1014
Energy [eV]

10 4 10 1 102 105 108 1011 1014
Energy [eV]

Synchrotron (Ec = 120 TeV)IC (Ec = 1200 TeV)

10 12
10 13
10 14

10 4 10 1 102 105 108 1011 1014
Energy [eV]

Hadronic-dominated model for tailE2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

10 12

Synchrotron (Ec = 300 TeV)

Synchrotron (Ec = 1200 TeV)

Hadronic-dominated model for head
10 11

Bremsstrahlung (Ec = 300 TeV)

IC-dominated model for tail10 4 10 1 102 105 108 1011 1014
Energy [eV]

HadronicIC
Synchrotron

Synchrotron (Ec = 1.5 TeV)IC (Ec = 120 TeV)

E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

IC-dominated model for head
IC
Synchrotron

Bremsstrahlung (Ec = 1.5 TeV)

10 11

HadronicIC
Synchrotron

10 12
10 13
10 14

10 4 10 1 102 105 108 1011 1014
Energy [eV]

Figure 4.1: Modeling of the SED of SNR G106.3+2.7. The left and right panels show the resultsof the head and the tail, respectively. The top, middle, and bottom panels show the results of
the Bremsstrahlung- and IC- and hadronic-dominated models, respectively. The white circles
show the radio flux of each region (Pineault & Joncas, 2000). The black bow-tie area showsa power-law fit and 1σ statistical errors measured by Suzaku (Fujita et al., 2021). The open
triangles show the Fermi -LAT measurements (Xin et al., 2019; Liu et al., 2020). The markers in
the TeV γ-ray band are the same as in Fig. 3.8, but those corresponding to the MAGIC data areshown in black here. The red, magenta, blue and green lines show the hadronic, Bremsstrahlung,IC, and synchrotron emission. The model parameters of each panel are summarized in Table 4.1.
In the top two and middle-right panels, the solid and dashed lines show the emissions with theelectron cutoff energy indicated in each legend.

head
tail
head
tail
head
tail

Bremsstrahlung
dominated

1.9
1.9
2.6
2.6
2.5
2.5

αe

Ec,e
We
B
[TeV]
[erg]
[µG]
2.0 (300)† 7.4 × 1044 (1.7 × 1045 )† 190
1.5 (300)† 5.9 × 1044 (1.4 × 1045 )† 230
360
1.4 × 10473†47
120 (1200)
1.6 × 10360
1.8 × 1046
1035
2.0 × 1046
10
1.7
1.7

αp

Ec,p
[TeV]
60
1000

Wp
[erg]
8.9 × 1045
8.2 × 1045

Ngas
[cm−3 ]
100
200
100
200

χ2 /ndf
25.3/7 (≫ 1)†
148.7/31 (≫ 1)†
5.0/7
103.1/31 (≫ 1)†
5.3/7
39.9/31† In the top-right panel of Fig. 4.1, the model curve using the value in the parentheses is shown with the dashed line.

Hadronic
dominated

IC-dominatedSourceModel

Table 4.1: Model parameters for reproducing the observed spectra. α and Ec are the power-law index and the cutoff energy of theparticle spectrum, respectively. W is the total energy of particles with energy above 1 GeV. The subscript e and p denote electronsand protons. B is the magnetic-field strength in unit of µG. Ngas is the gas density in unit of cm−3 .

CHAPTER 4. SPECTRAL MODELING AND IMPLICATION TO THE γ-RAY ORIGIN OF
SNR G106.3+2.7
70

CHAPTER 4. SPECTRAL MODELING AND IMPLICATION TO THE γ-RAY ORIGIN OF
SNR G106.3+2.7
71

4.2

Interpretation of the γ-ray origin

4.2.1

head region

The X-ray emission in the head region exhibits a softening of the spectral index with distancefrom the pulsar, suggesting the emission originates in electrons accelerated in and propagatedfrom the shock of the PWN (Ge et al., 2021). Our modeling result shows that X-ray and gammaray fluxes can be explained with IC-dominated model from the same electron population. It thus
implies that the γ-ray emission can originates in the PWN. Assuming the electron cutoff energy(Ec,e ) and the magnetic-field strength (B) used in the leptonic model for the head region, theelectron lifetime due to synchrotron losses is given by: ∼ 3.9 kyr (Ec,e /360 TeV)−1 (B/3 µG)−2 ,which is consistent with the age of the SNR estimated to be 3.9 kyr from the spectral breakin the radio spectrum of the PWN (Kothes et al., 2006) or 10 kyr from the pulsar spin-down
age (Halpern et al., 2001a). With the spin-down luminosity of Lsd = 2.2 × 1037 ergs−1 , thetotal electron energy from the PWN can be roughly approximated to be Lsd × tage = 2.8 ×
1048 erg (tage /4 kyr), which is sufficiently larger than the current model parameter (We =1.4 × 1047 erg).
The hadronic-dominated model also works for the head. The protons accelerated up to
60 TeV can explain the VHE γ-ray emission detected by MAGIC, given the presence of dense
HI clouds in the head region pointed out by Kothes et al. (2001). Although CO emission is
not prominent, the total intensity of HI and CO line emissions suggests the presence of gaseswith a total proton density of ∼ 100 cm−3 , which is sufficient for the pp emission, as derived
in Appendix A. Still electrons with a largely different spectral index are needed to explain theradio and X-ray emission. One of the simplest explanation would be that electrons are mainly
from the PWN, while the protons were accelerated in the shell. While SNRs can accelerate CRs
only for the limited periods (see Sect. 1.1.2), the shock wave of an SNR with an age of ∼ 4 kyr
is sufficient to produce protons up to 60 TeV (Ptuskin & Zirakashvili, 2003, 2005), which is
required by the modeling result.

4.2.2

tail region

The modeling described in the previous section suggests that it is difficult to explain the tailemission with the Bremsstrahlung- and IC-dominated models. On the other hand, the hadronicdominated model worked well. The γ-ray spectrum of the tail region can be reproduced assuming
a proton maximum energy of 1 PeV (χ2 /ndf = 39.9/31). Generally speaking, acceleration up
to 1 PeV can only be achieved at the early stages of the SNR evolution (∼ 0.1 kyr) (Ptuskin &Zirakashvili, 2003, 2005). However, as mentioned in Sect. 4.2.1, the age of this SNR has been
estimated to be 3.9 kyr from the spectral break in the radio spectrum of the PWN (Kothes et al.,2006) or 10 kyr from the pulsar spin-down age (Halpern et al., 2001a). This discrepancy in the
SNR age can be solved assuming a CR-escape scenario (e.g. Aharonian & Atoyan, 1996; Gabici& Aharonian, 2007). In this scenario, protons accelerated up to ∼ PeV energies at a young SNR
escape from acceleration regions and illuminate nearby clouds, which produce “delayed” γ-ray
emission. This scenario can also explain a proton index of 1.8, harder than 2.0 expected fromCHAPTER 4. SPECTRAL MODELING AND IMPLICATION TO THE γ-RAY ORIGIN OF
SNR G106.3+2.7
72electrons

SNR
tage ~ 4 kyr(Emax ~ 35 TeV)Radio & X-ray Synchrotron
protons

Observer

~ PeV

pp→π0
PeV protons escaped at tage ≲ 1 kyrare currently emitting γ rays at Cloud?

Cloud

π0-decay γ-ray

Figure 4.2: Schematic view of the interpretation of the non-thermal emissions from the tailregion. tage indicates the age of SNR (i.e., the time interval after the SN explosion).
the DSA theory, as described in Sect. 1.2.1. On the other hand, it requires high density clouds
spatially coinciding with the γ-ray morphology. Using the CGPS data of HI and 12 CO (J = 1−0)(see Appendix A), we confirmed a coincidence of the γ-ray emission with CO line emission in
the velocity range −6.41 to −3.94 km s−1 in the tail region, which was already pointed out
by Kothes et al. (2001) and Acciari et al. (2009). This supports the CR-escape scenario in the
tail region. Moreover, the scenario is consistent with the interpretation given in Albert et al.
(2020), Fujita et al. (2021), and Amenomori et al. (2021b). The authors estimated the diffusion√
length of CRs using the relation: ldiff = Dt, where D is the diffusion coefficient, and t is the
diffusion time. They then found, even assuming a small diffusion coefficient (D ∼ 1026 cm2 s−1
at GeV), that the diffusion length for CRs with an energy of O(100 TeV) in 5–10 kyr is larger(40–60 pc) than the size of the SNR (∼ 6 pc) and thus suggested the CRs are not confined in
the SNR.
Electrons may also escape in the same way as protons but be affected by radiative cooling,which is not considered in the modeling. However, the change in the spectral index due to thecooling effect is estimated to be at most 0.1–0.4 (Diesing & Caprioli, 2019), suggesting thatthe difference (∼ 0.7) between the proton and electron indices cannot be explained even byconsidering it. This fact implies that leptonic and hadronic emissions may happen at different
locations and thus under different physical conditions. For example, as shown in the schematic
view in Fig. 4.2, lepton emissions come from the SNR shell and hadron emissions come frominterstellar gas spatially separated from the SNR. This assumption can allow the unusual ratioof the total energy of CRs (Wp ≲ We ) because only the hadronic emission is affected by the
propagation effect (Gabici & Aharonian, 2007), and thus only Wp decreases.
The hard proton index found in the TeV band can also be explained with the SNR-cloud interactions (Inoue et al., 2012), as an alternative to the CR-escape scenario. The thermal instability
due to the shock-cloud interaction is expected to generate a clumpy ISM, where the magnetic field
may be amplified up to B = 1 mG. The penetration length of CRs in the clumpy ISM depends

CHAPTER 4. SPECTRAL MODELING AND IMPLICATION TO THE γ-RAY ORIGIN OF
SNR G106.3+2.7
73
on the energy and is given by lpd = 0.1η 1/2 (E/10 TeV)1/2 (B/100 µG)−1/2 (tage /1 kyr)1/2 pc,
where η is the gyro factor (δB 2 /B 2 ), E is the particle energy. Thus, it is acceptable that the
proton index is ≲ 0.5 harder than expected in DSA theory, as required in the modeling forthe tail region. However, this model cannot explain the maximum energy of ≳ PeV in SNRs
older than 1 kyr and the differences in the distribution of electrons and protons as mentioned
above (Diesing & Caprioli, 2019).

4.2.3

Remarks on the modeling

It should be noted that the data points of Milagro, HAWC, Tibet ASγ, and LHAASO, includedin the modeling of the tail spectrum, are from extraction regions, which include the head. Hence,
they are potentially contaminated if the head emits radiation > 10 TeV. Based on the present
data alone, we cannot strongly claim that the origins of the γ rays from the head and the tailare different. Thus, to accurately determine the emission mechanism, it is necessary to separate
the extraction regions at the head and tail also for spectral points above 10 TeV. In addition,
the integrated region of MAGIC-tail in this analysis may miss a fraction of the γ-ray emissions
observed by air shower experiments. Using the Gaussian extension at >6 TeV derived with
θ2 plot around tail, the survival rate is estimated to be 73.9–94.5% (1σ uncertainties). We
examined the effect on our model fit for the tail spectrum when using the scaled down fluxof air shower experiments by 26.1%. In the leptonic model, χ2 /ndf changed (from 103.1/31
to 96.3/31.1) slightly, indicating the model is still inconsistent with the observed data. In the
hadronic model, χ2 /ndf also changed (from 39.9/31 to 41.3/31), and the model still works. As
a result, the underestimation of the flux due to the angular cut does not affect our conclusion.
In conclusion, the spectral modeling indicates an evidence that the 100-TeV γ ray originatesfrom PeV proton accelerated by SNR G106.3+2.7. However, this conclusion is based on the qualitative interpretation by the CR-escape scenario, which has many uncertainties. To justify thisinterpretation, it would be necessary to pass two tests: whether the SNR can accelerate protons
up to ≥ 1 PeV in the past and whether the observed spectrum can be reproduced by considering CR diffusion. These questions will be approached using the time-evolution measurement of
the maximum energy in a single SNR, which will be attempted with SNR HB9 observations in
Chapter 5. We will then revisit this interpretation for SNR G106.3+2.7 in Chapter 6.

Chapter 5

γ-ray observations of SNR HB9
toward the measurement of the DSA
history
Although we have obtained an evidence that the SNR G106.3+2.7 examined in the previous
chapters is almost certainly a PeVatron, the interpretation by the CR-escape scenario remains
a qualitative discussion and requires a quantitative study. Also, a study of whether such SNRs
are universal would require a measurement of the time evolution of the maximum accelerationenergy within the above interpretation. In this chapter, we describe a method to measure
the time evolution of maximum energy in SNRs and attempt to establish the new method byadapting it to the SNR HB9 observations.

5.1

A measurement method of the maximum energy evolution
in an SNR

Systematic studies exploring the population of SNRs are expected to provide information on theevolution of the CR spectra of SNRs as a function of the SNR age (Ohira et al., 2011; Schure& Bell, 2013; Gaggero et al., 2018). Fig. 5.1 shows a correlation between γ-ray spectral cutoff(Ecut )and the age of SNRs (tage ). This correlation follows the power-law function:
Ecut = 1.3+1.1−0.63 TeV



tage
1 kyr

−0.81±0.24
,

(5.1)

which suggests that an acceleration up to PeV is not possible even by young SNRs with an ageof 100 yr (Suzuki et al., 2020, 2022). However, quantifying the evolution of DSA with this kindof method is challenging due to the large diversity of the surrounding environment of SNRs (e.g.,
Yasuda & Lee, 2019). This latter may affect the dominant radiative process for the gamma-ray
emissions and the estimation of the maximum particle energy.
Another way to study the evolution of DSA is a simultaneous observation of a single SNR,
specifically its shell part and nearby massive cloud regions. If a massive cloud exists in the74CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
75

Figure 5.1: Cutoff energy of γ-ray spectra of SNRs as a function of an estimated age (Suzuki
et al., 2022). The gray bow-tie area shows the fitting result with a simple power-law function.

SNR
8 kyr

0.4 kyr2 kyrtime evolution

Cloud
32 kyrProtons from the SNR illuminate the cloudand generate γ-ray via the π0-decay process.

8 kyr

2 kyr

SNR

CRs
Delayed bydiffusion time

γ-ray

Reflect the current age

Cloud“delayed”
γ-ray

Reflect past distribution

time evolution(based on Gabici+2007)

We can quantify a time evolution of CR distribution
by observing the two spectra simultaneously.

Figure 5.2: Schematic view of the measurement method for the time evolution of CR distribution
in a single SNR.

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
76
vicinity of the SNR, protons that escaped the SNR illuminate the cloud and generate γ-ray
emission via the π 0 -decay process (Aharonian & Atoyan, 1996; Gabici & Aharonian, 2007). The
delay of the timing of the γ-ray emission at the cloud region from that of the incident proton
escape depends on the propagation time and accordingly reflects the particle distribution in the
SNR at a specific epoch in the past. Fig. 5.2 shows (left) an example of time evolution in theenergy spectra at the shell and cloud regions and (right) a schematic view for this measurement.
A comparison of the spectra at the SNR shell and nearby clouds observed at roughly the sametime enables us to quantify the evolution of the DSA in the SNR.

5.1.1

Previous observations of delayed γ-ray emission

In observations of SNR W28, such “delayed” γ-rays have been detected from three cloud regionswith HESS and Fermi -LAT (Aharonian, 2008a; Hanabata et al., 2014). The observed γ-ray
spectra from these cloud regions, however, do not differ significantly from one another, andhence the maximum energy with the DSA at W28 has never been successfully measured as afunction of the SNR age. Since the angular distances from the SNR center to these individualclouds are almost the same in W28, the emissions from the clouds should originate in protons
accelerated at a similar epoch. The γ-ray emission from another cloud (named HESS J1801−233)
located to the north of W28 was suggested to originate in protons that have been re-accelerated
by a shock-cloud interaction (Cui et al., 2018), in which case the γ-ray spectrum does not reflect
the SNR age. Hanabata et al. (2014) reported detection of the fifth γ-ray emitting region, named
source W, located at the western part of W28. Although there is a possibility that the emission
is explained with protons escaping from the SNR, no cloud counterpart has been found. Thusfar, it is impossible to quantify the particle distribution at multiple epochs with the available
observations of W28.
One concern with the study using delayed γ-ray spectra is a large dependence on the diffusion
coefficient of CRs. The diffusion coefficient may be modified in the vicinity of an SNR due tothe effect of self-confinement and/or magnetic-field amplification caused by the generation of
turbulent plasma waves (Wentzel, 1974; Fujita et al., 2011; D’Angelo et al., 2018). In this
scenario, CRs traveling in the interstellar magnetic field generate plasma waves, which scatterthe CRs and prevent their diffusion. Furthermore, the CR pressure may increase the density
fluctuation of ISM, resulting in the growth of turbulence in the magnetic field. While delayed
γ-ray emissions are detected in the vicinity of several SNRs (e.g., Uchiyama et al. (2012)), thederived diffusion coefficients had large uncertainty, and so did the estimated amount of particle
acceleration at the SNR shell in the past. Fig. 5.3 shows the modeling result with the delayed
γ-ray emission for the observed spectrum around SNR W44. Their model spectra are almost
identical, even with diffusion coefficients that differ by one order of magnitude. This is mainly
because the SNRs discussed so far are relatively old. For older SNRs, clouds and SNR shells
often interact directly, and as a result, the observed data do not provide much information on
the diffusion coefficient (Fig. 5.3). Furthermore, even if clouds are significantly separated from
SNR shells, the delayed γ-ray spectra do not differ significantly from those of the current SNR
shells, and thus the data are not of much help for restricting the diffusion coefficients (Uchiyama

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
77

Figure 5.3: Modeling result to γ-ray emission from the cloud region around SNR W44 (Uchiyama
et al., 2012). The red and grey data points represent the energy spectrum at the cloud and theshell region, respectively. The blue lines show the model of the delayed γ-ray emission withthree cases of diffusion coefficient: D = [0.1, 1.0, 3.0] × 1028 cm2 s−1 (E/10 GeV)0.6 .
et al., 2012). This concern can be addressed by observing the spectral peak of the delayed γ-ray
spectrum at cloud regions. Fig. 5.4 shows an example of the delayed γ-ray spectrum and thephysical quantities associated with its features (i.e., spectral peak and index). The spectral
peak mainly depends on the diffusion coefficient as well as the time evolution of Emax and the
distance between the shell and cloud regions (dcl ). The energy spectrum at the shell part andthe two dimensional angular distance can give constraints on the evolution of Emax and dcl ,respectively. Thus, if its spectral peak is observable, the diffusion coefficient can be estimatedwith few uncertainties by fitting the delayed γ-ray spectrum.

5.1.2

Introduction of SNR HB9

Here, we focus on SNR HB9 (G160.9+2.6), which is relatively young (∼ 6.6 × 103 yr; Leahy& Tian (2007)) compared with other objects where delayed γ-rays have been observed. HB9
has two additional advantages for this type of study in the DSA evolution. Firstly, there aremolecular clouds in the vicinity of this SNR, but their locations do not overlap with the SNR in
the line of sight (Sezer et al., 2019), as shown in the left panel of Fig. 5.5. This allows to obtain

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
78
Before escaping (at Shell)
After escaping (at Cloud)Not reached yet wellspectral peakTime evolution of Emax
Diffusion coefficient (D0)Distance b/w SNR & Cloud

Already escaped

spectral index
Index at Shell
Index of diffusion (δ)→ Emax at tsedov

Figure 5.4: Example of the “delayed” γ-ray spectrum.
spatially-resolved spectra in the shell and cloud regions and to estimate the diffusion time using
the angular distance between the SNR and clouds. Secondly, HB9 has observable γ-ray emission
from the SNR shell (Araya, 2014), which is essential to estimate the current maximum energyof the accelerated particles at the SNR shock (the right panel in Fig. 5.5).
HB9 has a disk-like morphology with a radius of ∼ 1◦ in the radio band (e.g., Leahy & Roger,1991), and its spectral index is α = −0.47±0.06 in frequencies between 408 and 1420 MHz (Leahy& Tian, 2007), which is consistent with the typical value (α = 0.5) of SNRs (Reynolds et al.,2012). Its kinematic distance was estimated to be 0.8 ± 0.4 kpc from HI observations, yielding
the dynamical age estimate of 6.6 × 103 yr (Leahy & Tian, 2007). Sezer et al. (2019) found
an HI shell expanding toward the SNR with a velocity ranging between −10.5 and +1.8 km
s−1 and derived a kinematic distance to be 0.6 ± 0.3 kpc, which is consistent with the estimate
by Leahy & Tian (2007). In the γ-ray band, spatially extended emission was detected fromSNR HB9 along with its radio shell with the Fermi -LAT 5.5-year observations in the energy
band above 0.2 GeV (Araya, 2014). The γ-ray spectrum was characterized by a power law
with a photon index of 1.44 ± 0.25 accompanied by an exponential cutoff at 1.6 ± 0.6 GeV. The
spectrum is explained with emission from inverse Compton (IC) scattering of electrons with asimple power-law energy spectrum with a differential index of 2 and maximum electron energyof 500 GeV. Furthermore, Sezer et al. (2019) analyzed the Fermi -LAT 10-year data in an energy
range between 0.2 and 300 GeV and newly detected a point-like source near the SNR shell,
named PS J0506.5+4546. The spectrum of the point source can be characterized by a simple
power-law function with an index of 1.90 ± 0.19. The flux was (6.59 ± 3.47) × 10−10 cm−2 s−1 .
PS J0506.5+4546 was, however, suggested not to be related to the SNR shell, and its origin is
unclear (Sezer et al., 2019).
In this work, we study the γ-ray morphology of SNR HB9 and the spectra of the SNR shelland the nearby cloud regions, using 12-year observations with the Fermi -LAT, to quantify the(a) 12CO (J=1-0) line(c) CfA 1.2m CO (VLSR: −10.5 − +1.8 km s )

3
8

2
4
1

150.0

112.5

75.0

37.5

159

0.0

161
160
Galactic Longitude [degree]

159

48:00:00.0

X

0.065

0.036

0

-0.021

10:00.0

05:00.0

5:00:00.0

55:00.0

4:50:00.0-0.05(d) CfA 1.2m CO (B: +1.55 − +3.85 degree)

1.6

Figure 5.5: Skymaps in the vicinity of SNR HB9. (a) 12 CO (J = 1−0) line intensity map (Sezer
et-15.0al., 2019). Also shown with white contours is the radio continuum emission at 4850 MHz.
(b) GeV γ-ray map above 1 GeV (Araya,
2014). The green and magenta contours indicate the
1.2-10.0 continuum at 4850 MHz and X-ray emissions, respectively.
radio
-5.0

H

ll
she
ng
ndi
a
p
x
IE

0.8

K degree)162

VLSR [km s−1]

159

500

0.094

0.007645:00:00.0

600

0.12

47:00:00.0

Galactic Latitude [degree]

K km s−1

700

12

0.15

46:00:00.0

4

800

16
N

E

Declination (J2000)

900

0.18

Right Ascension (J2000)

K km s−1)(b) γ ray > 1 GeV−1K degree1

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
79

evolution of DSA. In addition, since we suspect that the γ-ray emission from PS J0506.5+4546
originatesin a molecular cloud located outside the radio shell, we search for spatial correlation
0.0
between the 12 CO (J = 1 − 0) line and0.4γ-ray emissions. In Section 5.2, we describe the Fermi 5.0
LAT
observations and show the results of our analysis. We present the modeling study tointerpretthe origin of the observed γ-ray
in Section 5.4.
10.0
0.0
162

161
160
Galactic Longitude [degree]

159

5.2

Fermi-LAT observations in the vicinity of HB9

5.2.1

Fermi -LAT

The Large Area Telescope (LAT) onboard Fermi satellite, launched in 2008, is capable of detecting γ rays in an energy range from ∼ 200 MeV to > 300 GeV (Atwood et al., 2009) and,as mentioned in Sect. 2.5 has better sensitivity at < 10 GeV than the MAGIC telescopes. Furthermore, the wide field of view of 2.4 sr can provide abundant data over the entire sky, which
is useful in investigating unknown γ-ray sources.
As shown in Fig. 5.6, the LAT detector consists of a tracker, a calorimeter, and an anticoincidence detector, with a total size of 1.8 m × 1.8 m × 0.72 m. The detection principle ofγ rays in the Fermi-LAT is similar to that of the IACT described in Chapter 2. In the LAT
detector, an incoming γ ray causes pair production at the tungsten sheet of the tracker part,
generating an electromagnetic shower in the calorimeter. The arrival direction of the gamma
ray can be estimated by capturing the electron track of the first electron pair in the tracker,
while its energy can be measured in the calorimeter. Each component is described below.
TrackerCHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
80
Tracker
TKR front section

ThermalBlanket12 ´ 3% X0

TKR back section

DAQ
Electronics

Anticoincidence
Detector

CAL

4 ´ 18% X0
2 ´ no W

8.6 X0Calorimeter

Figure 5.6: Left: Schematic view of the Fermi-LAT detector, based on Fig. 1 of Atwood et al.
(2009). Right: Schematic diagram of the tracker and calorimeter within the LAT detector,
adopted from (Ackermann et al., 2012b). X0 indicates the radiation length.
The tracker consists of a total of 18 layers of silicon strip detectors and the tungsten sheets,the later of which are installed only at the front 16 layers. Since only events where electron
pairs pass through three or more silicon strips are used for analysis, there are no tungsten
sheets in the last two layers. The reason there are no tungsten sheets in the last two layers
here is that only events in which electron pairs pass through three or more silicon sheets
can be used in the analysis. The thicknesses of the tungsten sheets are 0.03 cm for thefront 12 layers and 0.18 cm for the following 4 layers, which correspond to 3% and 18%of the γ-ray radiation length, respectively. The combination of the thicknesses has been
determined by balancing the angular resolution and the effective area. When an incident
γ ray causes pair production at the thin sheet, the multiple scattering effects of electrons
are suppressed compared to the case of thick sheets and thus the arrival direction canbe estimated with high accuracy. On the other hand, a thicker sheet provides a larger
effective area, which helps to detect high-energy γ rays.
Calorimeter
The calorimeter uses 96 (12 × 8) CsI(Tl) crystals, the each size of which is 2.7 cm × 2.0 cm
× 32.6 cm. Since we know which scintillator the shower particles pass through, we canreconstruct the electromagnetic shower and thus reduce the CR background events.
Anti-coincidence detector
Plastic scintillators cover the tracker and the calorimeter. They allow identifying thebackground events due to CRs passing through from the side of the detector.
Fig. 5.7 show the angular resolution of the LAT detector as a function of energy, the low
energy part of which roughly follow with the relation of E −0.8 (Ackermann et al., 2012b). It
should be noted that the angular resolution in the front 12 thin layers is about twice better thanthat of the thick 4 layers in the back.

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
81

Total
Front Layers
Back Layers

68% Containment angle [deg]

102

101

100
E 0.8
10 1
10 2

10 1

100

101
Energy [GeV]

102

103

Figure 5.7: Fermi-LAT angular resolution as a function of energy (https://www.slac.
stanford.edu/exp/glast/groups/canda/lat_Performance.htm). The black points represent
the resolution for all events, while the red (blue) points represent the events where first pair
production occurs in the front-thin (back-thick) layers of the tracker part.

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
82

5.2.2

Data selection and analysis

We analyze its 12-year data from 2008 August to 2020 August in the vicinity of SNR HB9. The
standard analysis software, Science Tools (version v11r5p31 ), is used. The “Source” selection
criteria and instrument responses (P8R2_SOURCE_V62 ) are chosen, considering a balance between
the precision and photon-count statistics. The zenith-angle threshold is set to 90◦ to suppress
the contamination of the background from the Earth rim. We employ the tool gtlike (in thebinned mode), using a standard maximum likelihood method (Mattox et al., 1996), for spatialand spectral analyses. We choose a square region of 15◦ ×15◦ with the center coinciding withthat of HB9 (Ra=75.25◦ , Dec=46.73◦ ) as the region of interest (ROI) for the (binned) maximum
likelihood analysis based on Poisson statistics. The pixel size is 0.1◦ .
The source spatial-distribution model includes all the sources in the fourth Fermi catalog (4FGL; Abdollahi et al., 2020) within the ROI and the two diffuse backgrounds, the Galactic (gll_iem_v7.fits) and extragalactic (iso_P8R3_SOURCE_V2_v1.txt) diffuse emissions. Regarding the emission of the SNR shell, Araya (2014) and Sezer et al. (2019) concluded that theradio template produced with the 4850 MHz radio continuum data from the Green Bank Telescope (Condon et al., 1994) is the best spatial model. Accordingly, we use the radio template
provided in the Science Tools as the spatial model. In the fitting of the maximum likelihood
analysis, all spectral parameters of HB9 SNR itself, 4FGL sources (Abdollahi et al., 2020) located within 5◦ from the center of HB9, and the two diffuse backgrounds are allowed to vary
freely. We do not use the data below 1 GeV in this analysis since the fitting results of delayed
γ-ray emission in this band suffer from the systematic uncertainty in the Galactic diffuse background model (Ackermann et al., 2012a). Note that the flux of the Galactic diffuse emission at
the molecular cloud region is larger in this energy range than that of the delayed γ-ray emission
from HB9, which will be estimated in Section 5.4.
The significance of a source is represented in this analysis by the Test Statistic (TS) definedas −2log(L0 /L), where L0 and L are the maximum likelihood values for the null hypothesis
and a model including additional sources, respectively (Mattox et al., 1996). The detection√
significance of the source can be approximated as TS when the number of counts is sufficiently
large.

5.3

Results

5.3.1

Morphology study

Figures 5.8 a and c show the background-subtracted γ-ray TS map created from the Fermi -LAT
12-year data above 1 GeV, where the background model consists of the Galactic and isotropic
extragalactic emissions and the contributions from the known Fermi sources (see the previous
subsection). The map is overlaid with cyan contours of the 12 CO (J = 1 − 0) line emission fromthe Dame survey data (Dame et al., 2001), which are integrated over a velocity range between
12

https://fermi.gsfc.nasa.gov/ssc/data/analysis/software
https://fermi.gsfc.nasa.gov/ssc/data/Cicerone/CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
83−10.4 and +2.6 km s−1 , and also green contours of 1420 MHz radio continuum emission obtained
from the CGPS survey with DRAO (Taylor et al., 2003). In Fig. 5.8 c, no significant emission
from the SNR shell was found, which is consistent with the fact that the SNR spectrum has acutoff below 10 GeV (Araya, 2014). Fig. 5.8 b is a similar γ-ray TS map to Fig. 5.8 a, but withthe contribution from the SNR shell is additionally subtracted, for which we use the 4850 MHz
radio-template model as a spatial model and assume a simple power-law spectrum. The γ-ray
excess in Fig. 5.8 b and c appears to be more extended than the point source J0506.5+4546
reported in the previous study and rather spatially coincident with the CO line emission.
We test whether the γ-ray (> 1 GeV) spatial distribution is correlated with CO line emissions
using the likelihood method. The 12 CO (J = 1 − 0) line image exhibits two distinctive regions(thus two molecular clouds), designated as R1 and R2 (Figure 5.8). For the test, we create a CO
template for each of R1 and R2, which is made from the 12 CO (J = 1−0) line image (Dame et al.,2001) integrated over a velocity range between −10.4 and +2.6 km s−1 and cut with a thresholdof > 4.5 K km s−1 (cyan thick-line contours in Figure 5.8). We assume the γ-ray spectra of thecloud regions and HB9 SNR shell to follow a simple power-law function of dN/dE ∝ E −Γ . Here,(i) the background model (the null hypothesis) consists of the radio template for the HB9 SNRshell as well as the Galactic and extragalactic diffuse emissions and 4FGL sources (Abdollahi
et al., 2020). To compare spatial models, we use the Akaike information criterion (AIC; Akaike,1974) defined as 2k − 2 log(L), where k is the number of the estimated parameters in the model.
Consequently, we find that the model with the smaller AIC value is favored and that the bettermodel gives a larger ∆AIC (= (AIC)0 − (AIC)m ), where (AIC)0 and (AIC)m are the AIC values
for the null hypothesis and a model including additional sources, respectively. When we apply
the CO template models to the two cloud regions (ii), ∆AIC is found to be 44.4, indicating asignificant correlation between the γ-ray and 12 CO (J = 1 − 0) line emissions (Table 5.1). Since
the γ-ray emission from R1 spatially coincides with PS J0506.5+4546 reported by Sezer et al.
(2019), we also apply a point source model to the same position as PS J0506.5+4546, instead of
the CO template model of R1 (iii). The resultant AIC is marginally (1.3σ level) improved fromthe case with the CO template. As a result, the γ-ray emission from R1 is consistent with PS
J0506.5+4546, while R2 is newly detected with a statistical significance of 6.1σ in this study.

0Right ascension [deg]43.0

78.0 77.0 76.0 75.0 74.0 73.0 72.0

43.0

45.0

46.0

44.0

R2

R1

Declination [deg]

44.0

45.0

46.0

47.0

5

TS

Right ascension [deg]1078.0 77.0 76.0 75.0 74.0 73.0 72.0

R2R1(b) >1 GeV (SNR subtracted)

Declination [deg]

15

43.0

44.0

45.0

46.0

47.0

20

Right ascension [deg]25

78.0 77.0 76.0 75.0 74.0 73.0 72.0

R2R1(c) >10 GeV

Figure 5.8: γ-ray TS maps in the vicinity of SNR HB9 observed with the Fermi -LAT. All maps are given with square bins of 0.05◦ , andGaussian smoothing with a kernel σ = 0.1◦ are applied. The energy ranges are (Panels a and b) 1 to 500 GeV and (c) 10 to 500 GeV. The
subtracted background consists of (a and c) the Galactic and extragalactic diffuse emissions and known γ-ray source contribution and (b)additionally the estimated γ-ray emission from the HB9 SNR shell (see text for detail). Green contours show the radio emission of SNR HB9
at 1420 MHz with DRAO (Taylor et al., 2003) and are linearly spaced in increments of 0.5 K from 5.5 K to 10.0 K. Cyan contours show the
12 CO (J=1-0) line intensity integrated over a velocity range between −10.4 and +2.6 km s−1 and are linearly spaced in increments of 1.0
K km s−1 from 4.5 to 10.5 K km s−1 . The two apparent CO-emission regions (R1 and R2) are indicated. The magenta cross in each panel
indicates the position of PS J0506.5+4546 (Sezer et al., 2019).

Declination [deg]47.0(a) >1 GeV

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
84

Table 5.1: Likelihood test results for spatial models and resultant spectral parameters of HB9 SNR shell and the two cloud regions, using theenergy range of 1–500 GeV. The following three models are considered: (i) the background model (null hypothesis); (ii) additionally including
γ-ray emission from R1 and R2 estimated by the CO template model (see text); and (iii) assuming a point-source model, where the positionof the source is the same as PS J0506.5+4546 (Sezer et al., 2019), for R1 instead of the CO template model, while √the CO template as in(ii) for R2. All sources of interest are fitted with a simple power-law function of dN/dE ∝ E −Γ . The “Flux” and “ TS” columns refer√ to
the integral flux in the energy band of 1–500 GeV in 10−10 cm−2 s−1 and the statistical significance (σ) of each source approximated as TS,
respectively.
R2
HB9 SNR itself √R1√√Model ∆AIC
Flux
Γ
TS
Flux
Γ
TS
Flux
Γ
TS
(i)
0
24.9 ± 2.1 2.54 ± 0.06 13.2
(ii)
44.4
2.0 ± 0.8 1.84 ± 0.18
4.5
4.8 ± 1.2 1.84 ± 0.14
6.1
24.3 ± 2.1 2.55 ± 0.10 12.8
(iii)
58.2
1.3 ± 0.5 1.77 ± 0.18
5.8
4.8 ± 1.2 1.84 ± 0.14
6.1
24.5 ± 2.1 2.54 ± 0.10 13.0

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
85

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
86

5.3.2

Spectral results

We extract Fermi -LAT energy spectra from the radio SNR shell region and two regions R1
and R2, individually, using the CO template model for the cloud regions. Figure 5.9 shows theresultant spectra for an energy range between 1 and 500 GeV. We find that the obtained Fermi LAT spectrum of the SNR shell is consistent with the γ-ray spectrum reported by Araya (2014)(Figure 5.9). The respective best-fit parameters are summarized in Table 5.1. Here, a potential
concern about the spectrum of R1 is that our assumption of diffuse γ-ray spatial distribution
may not be appropriate, given the spatial proximity with the already identified point source
PS J0506.5+4546 (Figure 5.8). However, the discrepancy of the determined spectral properties(flux and index) between the results of the two spatial-distribution models is smaller than the1σ uncertainty (Table 5.1). Therefore, we conclude that the difference in the results due tothe difference between the assumed spatial distributions is not significant. We note that thepower-law index of R1 is consistent with that of PS J0506.5+4546 reported by Sezer et al.
(2019).
As mentioned in Sect. 5.2.2, the uncertainty of the Galactic diffuse background model may
affect the spectral results. Here, we examine the spectral variations due to the backgroundmodel parameters (normalization and index). We perform the likelihood fitting in the following
four cases for the Galactic diffuse emission:· both normalization and index are free parameters,· only normalization is a free parameter,· only index is a free parameter,· both normalization and index are fixed.
Fig. 5.10 shows the energy spectra of the cloud regions obtained with those settings. The data
point at 1 GeV fluctuates by a 1.5σ level of statistical error, while the rest of the points are
not affected significantly. The resultant systematic variety of the spectral fitting results is at
most a 1σ level of statistical error. Thus, we conclude that the systematic uncertainty on the
background model is not critical in the modelling.

5.4

Spectral modelling and its implication

In this work, we explore the possibility that the R1 and R2 are attributed to the delayed γ-ray
emission from molecular clouds illuminated by the CRs accelerated in HB9. In the following
section, we calculate the delayed γ-ray spectra in the molecular cloud regions, R1 and R2, using
the method described in Section 5.4.1, and try to fit the observed spectra of the cloud regions
and the SNR shell simultaneously. According to Araya (2014), once the γ-ray spectrum of the
SNR shell is modeled with the hadronic emission, it would require enormous explosion energyof supernova or a dense density of the interstellar medium (ISM). Hence, we assume that theγ-ray emission in the SNR shell originates mainly from the leptonic processes. On the other
hand, we consider two cases for the origin of γ-ray emissions from the cloud regions: one where
the leptonic emission dominates (Leptonic-dominated model), and the other where the hadronic

E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY87

HB9 SNR shell [Araya 2014]
HB9 SNR shell (this work)

10 11
10 12
10 13

HB9 SNR Shell (this work)
R1 (this work)
R2 (this work)

10 11
10 12
10 13
109

1011
1010
Energy [eV]1012

Figure 5.9: Spectral energy distributions measured with the Fermi -LAT (data points) and fit
results in the γ-ray band. Top: Black and open green stars show the data of the SNR shellobtained in this work and Araya (2014), respectively. Bottom: Red and blue data represent
the spectra of cloud regions R1 and R2, respectively. Also shown as the shaded black region is
the fit result of the SNR shell.

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
88
0-decay from BGCR
norm = free & index = freenorm = fix & index = freenorm = free & index = fix
norm = fix & index = fix

E2dN/dE [erg cm 2s 1]

10 11

10 12

10 13

R2

0-decay from BGCR
norm = free & index = freenorm = fix & index = freenorm = free & index = fix
norm = fix & index = fix

10 11
E2dN/dE [erg cm 2s 1]

R1

10 12

10 13
109

1010
Energy [eV]

1011

1012

109

1010
Energy [eV]

1011

1012

Figure 5.10: Comparison in energy spectra of the cloud regions obtained with several background model parameters. Also shown the gray shaded region is the hadronic emission from thebackground CRs, which will be derived in Sect. 5.4.
emission dominates (Hadronic-dominated model) in the GeV band. We present the results ofthe modeling in Sect. 5.4.2 and discuss implications on the model parameters in Sect. 5.4.3.

5.4.1

Modeling of delayed γ-ray emission

The energy spectra of the γ-ray emissions from the molecular clouds around HB9 are calculated
following the method proposed by Gabici et al. (2009) and Ohira et al. (2011). They considered
the escape process of the CR protons, which are accelerated by an SNR, from a shock front
into interstellar space during the Sedov phase. Escaped CR protons emit high-energy γ-rays
via π 0 production through pp collision in molecular clouds in the vicinity of the HB9 SNR. In
addition, we also discuss the contribution from inverse Compton (IC) scattering and non-thermal
bremsstrahlung from CR electrons escaping from the SNR. We evaluate later in this section the
expected γ-ray fluxes from the R1 and R2 regions in conjunction with the energy spectra of theescaped CR protons and electrons.
The radio continuum contours (Figure 5.8) show a circular symmetric morphology of SNR
HB9, suggesting it has probably maintained a spherically symmetric structure throughout its
evolution. Therefore, we assume that the SNR shell itself and the escaped CR distribution are
spherically symmetric. Based on X-ray observations, the energy of the supernova explosion (ESN )was estimated to be (0.15–0.30)×1051 erg, and the relation between the explosion energy andthe hydrogen density of ISM (nISM ) was also estimated as (ESN /1051 erg)(nISM /1 cm−3 )−1 =
5 (Leahy & Tian, 2007). In this work, we adopt 0.3 × 1051 erg as the explosion energy, inwhich case the density is given as 0.06 cm−3 . If the initial velocity of the blast wave wasush,0 = 109 cm s−1 , HB9 entered the Sedov phase at the age tSedov = 3.6 × 102 yr when its radius
was 3.6 pc. Given that the observed radius of HB9 is ∼ 10 pc, which is larger than 3.6 pc, HB9
must be already in the Sedov phase. Leahy & Tian (2007) estimated the dynamical age of SNR
HB9 to be 6.6 × 103 yr, which we adopt as the age of the SNR in this discussion.
First of all, let us model the energy spectrum of CR protons. During the Sedov phase, themaximum energy of the CR protons accelerated at the SNR shock is determined by the timescale

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
89
of escape from the acceleration region, and decreases with time (Ptuskin & Zirakashvili, 2005;
Caprioli et al., 2009; Ohira et al., 2010). The temporal evolution of the maximum energy
depends on non-linear processes, such as amplifying the magnetic field at the shock, which isstill theoretically unclear. Thus, we adopt the phenomenological power-law dependence of thecutoff energy Eesc of the proton spectrum at the SNR shock on the age t of the SNR, as discussed
by Gabici et al. (2009) and Ohira et al. (2010):

Eesc (t) = Emax

−α

t

,

tSedov(5.2)

where Emax is the maximum energy of the CR protons at tSedov and is set to 3 PeV, andα is the power-law index that is determined so that the current maximum energy (Enow ) isequal to Eesc (t) at the current age tage = 6.6 × 103 yr (Leahy & Tian, 2007). Here, we adopt
Enow = 300 GeV, which is the same value as that of the electrons at the SNR shell at present (seealso Section 5.4.2). However, in general, the maximum energy of CRs in the SNR shell depends
on the particle species and is determined by the balance between acceleration and escape fromthe shell or energy loss due to radiative cooling. Araya (2014), assuming that radiative cooling
limits the maximum CR electron energy, estimated that the gyro factor would be a too large
value of ∼ 660 compared to the standard SNR value. In such a case, as shown in Ohira et al.
(2012), the escape process determines both the highest energies of CR electrons and protons.
Therefore, according to Ohira et al. (2012), we can assume that the current maximum energyof electrons is the same as that of the protons. The timescale for a particle with the energy
E, tesc (E), to escape into interstellar space from the supernova explosion, as a function of theparticle energy E is given by
tesc (E) = tSedov

E
Emax

−1/α
.

(5.3)

The distribution function fp,out per unit energy per unit volume of the CR protons escaping
into interstellar space can be obtained by solving the transport equation,∂fp,out(t, r, E) − DISM (E)∆fp,out (t, r, E) = qp,s (t, r, E),
∂t

(5.4)

where t is the age of the SNR, r is the distance from the SNR center, DISM is the diffusion
coefficient in the ISM, and qp,s is the injection rate of the CR protons from the SNR shock into
interstellar space per unit energy, unit volume, and unit time. We adopt the following form of
the diffusion coefficient DISM :
DISM (E) = D0

E
10 GeV

δ
,

(5.5)

where D0 is the diffusion coefficient of CRs at E = 10 GeV. We assume D0 = 3 × 1028 cm2 s−1
and the index δ = 1/3, the latter of which is consistent with the Galactic mean expected in the
CR propagation model (e.g., Blasi & Amato, 2012). Following Ohira et al. (2011), we assume

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
90
that CRs with the energy E are injected from the SNR shell Resc = ush,0 tSedov (tesc (E)/tSedov )2/5
at t = tesc (E). We also assume that the energy spectrum of the injected CRs is monochromaticwith the energy E and that particles start to escape at any given time (see equation (5.3)).
Taking account of these assumptions, the injection rate qs is given by
qp,s =

Nesc (E)δ(r − Resc (E))δ (t − tesc (E)) ,
4πr2

(5.6)

where r is the displacement from the center of the SNR, Nesc is the spectrum of all the CRs
that have escaped from the SNR up to the present time. In this work, we assume that the
total energy of the escaped CR protons is proportional to the explosion energy of the supernova
explosion ESN , namely
Z
Emax

ENesc (E)dE = ηESN ,

(5.7)

Enow

where η is the acceleration efficiency coefficient. Here, we assume that Nesc is a power-law
function of E with an index pesc . According to Ohira et al. (2010), it is shown that the index
pesc of CRs escaping from the SNR is steeper than the index pSNR of CRs confined in the SNR
shell. While the value of pesc is determined by the time evolution of the maximum energy of
CRs at the SNR shell and the CR production rate, we treat pesc as a parameter. By assuming
a power-law form of Nesc , the integral of the equation (5.7) yields,

Nesc (E) =


−1 


−pesc η(2−pesc )ESN × Emax 2−pesc − 1
E
×
2
Enow




Enow

Enow

ηESN
−2
ln(Emax /Enow ) E

pesc ̸= 2

.

(5.8)

pesc = 2

The propagation models (Obermeier et al., 2012; Yuan et al., 2017) expect the spectral index ofparticles escaped from Galactic CR origin to be ∼ 2.4, which is adopted as a fiducial parameterof pesc in this modeling. The effect of this parameter on the modeling will be discussed in
Section 5.4.3. The energy spectrum of the escaped CRs is obtained by combining the transport
equation (5.4) and equations (5.6) and (5.8). Specifically, the solution of equation (5.4) for
E > Eesc (t) is given by, according to Ohira et al. (2011)"!
!#(r − Resc )2(r + Resc )2
Nesc (E)fp,out (t, r, E) = 3/2
exp −− exp −
2
2
Rd,p
Rd,p
4π rResc Rd,p(5.9)

where Rd,p is the diffusion length of a CR proton with the energy E, defined as
Rd,p ≡p

4DISM (E) (t − tesc (E)).

(5.10)In order to obtain the energy spectrum of CR electrons fe,out , it is required to solve thefollowing transport equation with radiative cooling:∂fe,out∂(t, r, E) − DISM (E)∆fe,out (t, r, E) +(P (E)fe,out (t, r, E)) = qe,s (t, r, E),∂t∂E(5.11)

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
91
where P (E) is an energy loss rate of CR electrons, and qe,s is an injection rate of CR electrons,
which can be written asqe,s = Kep qp,s ,(5.12)
where Kep is the ratio of the total energy of CR electrons to CR protons. As a cooling process,we consider only the synchrotron radiation, which is the most dominant effect for electrons withenergies above O(GeV) (see Sect. 1.3).
The temporal evolution of the maximum energy of CR electrons in the SNR shell is different
from that of protons due to radiative cooling. After entering the Sedov phase, the maximum
energy of CR electrons is first determined by the balance between the radiative cooling andacceleration, and, after that, by the balance between acceleration and escape from the SNR
shell, similar to the proton (Ohira et al., 2012). When these two phases switch, CR electrons
start to escape from the SNR shell and be injected into the interstellar space. In this work, weparametrize this switching time te as follows:
te = ξe tSedov ,(5.13)

where ξe is the time in the unit of tSedov at which the electron starts to escape, and in thelimit of ξe → 1, the injection is identical to that of the proton (see equation (5.6)). After te ,since the time evolution of the maximum energy of electrons is expected same as that of the
proton (Ohira et al., 2012), we use equations (5.2) and (5.3). Originally, ξe is given based on
the assumed environment (Ohira et al., 2012), but we set ξe = 1 to consider the limit where theγ-ray flux of leptonic emissions is maximized in this modeling. The solution to equation (5.11)
can be derived using the method described in Appendix B.1.
Considering only the synchrotron radiation as the cooling process (i.e. P (E) = Qsyn E 2 ,see equation (B.7) of Appendix B.1 for details), the solution for Eesc (t) < E < Eesc (te ) can be
written as:"!
!#
Kep Nesc (Ec ) Ec2
1(r − Rc )2(r + Rc )2
exp −− exp −.
fe,out (t, r, E) = 3/2
2
2
Rd,e
Rd,e
4π rRc Rd,e E 2 1 − Qsyn tc Ec /α(5.14)
Once the energy spectrum of the CRs has been obtained, the flux of the γ-ray emission
is calculated based on the neutral pion decay process in molecular clouds for CR protons, andinverse Compton scattering and the relativistic bremsstrahlung for CR electrons. Here we assume
that the molecular cloud is “optically” thin for the CRs and that the cloud consists of a spatially
uniform gas with the total mass Mcl . For simplicity, we also assume that the molecular cloud isa sphere of radius Rcl . The spectrum of CRs in the molecular cloud can be written as follows:
Z dcl +Rcl
Ncl (E) =dcl −Rclπr  2
Rcl − (r − dcl )2 fs,out (r)dr,dcl(5.15)

where dcl is a distance to the center of the molecular cloud from the center of the SNR, and s
is species (e=electron, p=proton). Note that, for the parameters used in this work, a spatial
gradient of the CR distribution is not important, and only a few % difference occurs even if

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
92
Table 5.2: Fiducial parameters used to calculate the model spectra
SNR parameters
SN explosion energy
Initial shock velocity
Age of the SNR
Distance to the SNR
Acceleration efficiency†
Electron to proton flux ratio†
Maximum CR energy at t = tSedov
Current maximum energy of CRs
Magnetic field in the SNR
Particle index in the SNR
Particle index after escaping from SNR
ISM parameters
Number density
Diffusion coefficient at E = 10 GeV
Index of dependence on E of diffusion
Magnetic field in ISM
Cloud parameters
Distance to the cloud from SNR
Radius of the cloud
Average hydrogen number density†
Mass of the cloud†Magnetic field in cloud

Symbol
ESN
ush
tage
η
Kep
Emax
Enow
BSNR
pSNR
pesc
Symbol
nISM
D0
δ
BISM
Symbol
dcl
Rcl
nH
Mcl
Bcl

0.3 × 1051 erg
109 cm s−1
6.6 × 103 yr
0.8 kpc
0.1 (0.003)
0.02 (1)3 PeV300 GeV
8 µG
2.0
2.4
0.06 cm−3
3 × 1028 cm2 s−1
1/3
3 µG
R1
R2
17.8 pc
39.4 pc
3.6 pc
7.0 pc−3150 (1000) cm
200 (1000) cm−3
730 (4800) M⊙ 7100 (36000) M⊙
3 µG† The values in parentheses are used in the leptonic-dominated model.

3f
Ncl ∼ 4πRcls,out (dcl )/3 is used.
The energy spectrum of CRs in the SNR shell, Ns,shell , is calculated consistently with the
distribution function of escaped particles. By using the normalization equation (5.7) same as
equation (5.8), the spectra at E < Eesc of CR protons and electrons in the SNR shell can be
written as
−pSNRE
Np,shell (E) = Nesc (Enow ),(5.16)
Enowand
Ne,shell (E) = Kep Nesc (Enow )E
Enow

−pSNR,(5.17)

respectively, where pSNR is the index of CRs that have not yet escaped from the SNR shell.
We fix pSNR = 2 in this work because the index pSNR of HB9 has been well determined byobservations of radio continuum (see Sect. 5.1.2).
To calculate the spectra of non-thermal emissions, we use the radiative code naima (Zabalza,2015). As for the seed photon fields in the IC process, we assume the cosmic microwave background and Galactic far-infrared (FIR) radiation, the latter of which was not considered in themodeling of Araya (2014). The energy density of the FIR radiation is estimated to be 0.099 eV
cm−3 at T = 27 K, using the model in the GALPROP package (Porter et al., 2008).

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
93

5.4.2

Application to the cloud regions around HB9

In order to calculate the γ-ray flux at the cloud regions, R1 and R2, the distance between themolecular clouds and SNR is required. By fitting the CO intensity map integrated over thevelocity range between −10.4 and +2.6 km s−1 with a two-dimensional symmetric Gaussian, we
determine the center positions of clouds R1 and R2 to be (l, b) = (161.8 ± 0.1◦ , 2.8 ± 0.1◦ )and (162.6 ± 0.1◦ , 1.6 ± 0.1◦ ), respectively, in the galactic coordinates. Then, the projecteddistances between the centers of the SNR and clouds R1 and R2 are calculated to be 17.8 and39.4 pc, respectively, using the distance to the SNR from the Earth of 0.8 kpc (Sect. 5.1.2
and Table 5.2). We treat these projected distances as the actual distances in our discussion,
although they should be considered as the lower limit of the true distances in reality due to theuncertainty of their locations along the line of sight. We also obtain the radii of the molecular
clouds from the standard deviation of the fitted Gaussian, which are 3.6 pc and 7.0 pc for R1
and R2, respectively.
For the spectral modeling, the data points in the radio band for the HB9 SNR shell are
obtained from the literature (Dwarakanath et al., 1982; Reich et al., 2003; Leahy & Tian, 2007;
Roger et al., 1999; Gao et al., 2011), while the radio flux of the mean local background including
the cloud region at a frequency of 865 MHz (Reich et al., 2003) is used as an upper limit
for the cloud regions. In the X-ray band, only the thermal emission from the hot gas inside
the shell has been detected, while non-thermal emission has not been measured (e.g, Leahy &Aschenbach (1995); Sezer et al. (2019)). We adopt the 0.1–2.5 keV flux of the thermal emission
from HB9 (Leahy & Aschenbach, 1995) as an upper limit for the non-thermal emission in theenergy band.
We fit the spectra of R1 and R2 under two assumptions, the leptonic-dominated and hadronicdominated models, and simultaneously reproduce the spectrum of the SNR shell with the leptonic emission. The electron to proton flux ratio, Kep , is assumed to be 1 and 0.02 for the
leptonic-dominated and hadronic-dominated models, respectively, the latter of which is consistent with the ratio in the local CR abundance (e.g., Aguilar et al., 2016).
Figure 5.11 and 5.12 show the results of the leptonic-dominated and hadronic-dominatedmodel, respectively. The total electron energy and the magnetic field in the SNR shell are inagreement with those estimated by Araya (2014), while the electron maximum energy at the shell(corresponding to Enow ) is determined to be 300 GeV. Our obtained value of Enow is slightly
lower than that derived in the previous study (Araya, 2014) because the FIR radiation as a seedphoton in the IC process is newly taken into account in this work. In the leptonic-dominatedmodel, the γ-ray emissions via a bremsstrahlung process of relativistic electrons dominate in
the cloud regions, and thus the delayed γ-ray spectra for 1–500 GeV have a hard index of ∼1.3, which contradicts the observed one at the cloud regions (Figure 5.11). In the hadronicdominated model (Figure 5.12), the hadronic emission reproduces well the observed spectra
even though the assumed parameters in the calculation are typical ones for an SNR and ISM(Table 2). We note that the γ-ray flux from the background CRs is lower than the observed
data at higher energy than ∼ 3 GeV (see the middle and right panels of Figure 5.12). Thisis consistent with the fact that the systematic fluctuation due to the model parameters of theE2dN/dE [erg cm 2s 1]

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
94
10 10

HB9 shell

10 10

R1

10 10

10 11

10 11

10 11

10 12

10 12

10 12

10 13

10 13

10 13

10 14

10 14

10 14

10 15

10 15

10 4 100 104 108 1012
Energy [eV]
SynchrotronIC

10 4 100 104 108 1012
Energy [eV]Bremsstrahlung

10 15

R2

10 4 100 104 108 1012
Energy [eV]Total of -ray emission

Hadronic

Hadronic from BG

Figure 5.11: Broadband spectral energy distributions of the non-thermal emission from the HB9shell and the cloud regions with the leptonic-dominated model, using the parameters tabulated
in Table 5.2. The left, middle, and right panels show the results about the HB9 shell, R1 and R2,
respectively. The radio and X-ray data are taken from Dwarakanath et al. (1982); Reich et al.
(2003); Leahy & Tian (2007); Roger et al. (1999); Gao et al. (2011) and Leahy & Aschenbach(1995), respectively (see text for detail). The filled squares and open circles show the Fermi -LAT
data derived in Section 5.3.2 and Araya (2014), respectively. The lines represent each componentof the emission models: synchrotron (blue), electron bremsstrahlung (magenta), IC (green),
neutral pion decay (red), and the total γ-ray emissions (black). In order to demonstrate howmuch of the Galactic diffuse background at the cloud regions, which may affect the systematic
uncertainties in the Fermi analysis, we show the hadronic emission from the background CRs
with the energy spectrum JCR (E) = 2.2(E/GeV)−2.75 cm−2 s−1 GeV−1 sr−1 (e.g., Dermer, 1986)
as the shaded grey region.
Galactic diffuse emission is small at ≳ 3 GeV, as shown in Fig. 5.10. By calculating the energyof the proton corresponding to the peak of photon spectra in each cloud region with the fiducial
parameters, we estimate the epoch at which those protons escaped from the SNR shock by using
equation (5.3). We then obtain tref ≲ −5.8 × 101 yr for R1 and ≲ −6.3 × 102 yr for R2 datingback from now.

5.4.3

Accuracy of parameter determination

E2dN/dE [erg cm 2s 1]

We investigate in the following procedure how the hadronic-dominated model curve varies depending on the input parameters and summarize the results in Figure 5.13. First, we evaluate
10 10

HB9 shell

10 10

R1

10 10

10 11

10 11

10 11

10 12

10 12

10 12

10 13

10 13

10 13

10 14

10 14

10 14

10 15

10 4 100 104 108 1012
Energy [eV]
Synchrotron

10 15
IC

10 4 100 104 108 1012
Energy [eV]Bremsstrahlung

Hadronic

10 15
Total

R2

10 4 100 104 108 1012
Energy [eV]
Hadronic from BG

Figure 5.12: Same as Figure 5.11, but with the hadronic-dominated model.

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY95
the dependence of the model curve on D0 , which is the value of the diffusion coefficient at
E = 10 GeV (Panels a and b in Fig. 5.13). The Fermi -LAT spectra obtained in this work
are found to be well reproduced with the Galactic mean of D0 = 3 × 1028 cm2 s−1 . Orders ofmagnitude smaller D0 , in particular D0 = 3 × 1026 cm2 s−1 , however, clearly fail to reproduce
the observed spectra. In the previous studies for other SNRs (e.g., Fujita et al., 2009), the estimated values of D0 were ∼ 10 times smaller than the Galactic mean, which were explained inconjunction with self-confinement caused by the generation of turbulent plasma waves (Wentzel,
1974; Fujita et al., 2011; D’Angelo et al., 2018). Our diffusion coefficient value, which is close
to the Galactic mean, indicates that the excitation of such turbulent plasma waves at the distances to R1 and R2 is inefficient or that the wave damping has a significant effect. Second,
we find that the model with either δ = 1/3 or 1/2 is preferred over that with δ = 0 (Panels cand d in Fig. 5.13). Third, we also find that Emax above 10 TeV still explains the data points(Figures 5.13e and 5.13f). Given that there is a trend for a larger difference between the model
curves in the higher energy band, future observations in the TeV band will provide results moresensitive to determine the δ and Emax parameters.
In order to verify that the Galactic mean value of the diffusion coefficient (D0 ) is appropriate
to explain the γ-ray spectra, we also try to model the γ-ray spectra, fixing D0 at two orders ofmagnitude smaller D0 = 3 × 1026 cm2 s−1 . As can be seen in Figure 5.13(a, b), for smaller D0 ,the spectrum shifts to the higher energy side. Also, as shown in Figure 5.13(e, f), the spectrum

below O 1011 eV does not depend much on Emax . Therefore, the deviation of the spectrum
in the lower energy band due to the small D0 should be countered by δ dependence. Once weassume δ > 1.0 while other parameters are kept the same as in Table 5.2, the model curves
are roughly consistent with the observed γ-ray spectra (Figure 5.13(g, h)). However, δ > 1.0
is inconsistent with the CR propagation model (e.g., Blasi & Amato, 2012). Thus, the small
diffusion coefficient (i.e., D0 ∼ 1026 cm2 s−1 ) is unlikely.
We have assumed pesc = 2.4 in the modeling so far, but the discussion on the diffusion
coefficient may be affected by this parameter as well as δ. Here, we attempt to fit the observed
spectra by considering two cases where a flatter (pesc = 2.0) or steeper (pesc = 2.7) index isassumed. Note that these assumptions are not favored by the propagation model as mentioned
in Section 5.4.2. The results are shown in Figure 5.14 and 5.15. Although the model parameters
except for pesc and η are fixed at the values in Table 5.2, the observed spectra can be reproducedby giving a feasible value of η = 0.3 (0.05) for pesc = 2.0 (2.7) (Figure 5.14). Figure 5.15 shows
the dependency on D0 and δ for R2 as in Figure 5.13 (b) and (h). To explain the observed
spectrum with the diffusion coefficient two orders of magnitude smaller D0 = 3 × 1026 cm2 s−1than the Galactic mean, an unrealistic assumption that δ > 1.0 is required for both cases ofpesc = 2.0 and 2.7. This result is consistent with the case where pesc = 2.4, i.e., the diffusion
coefficient of the Galactic mean value is preferred regardless of the assumption on pesc .

10 12
10 13
10 14

10 11

(c) R1 -E2dN/dE [erg cm 2s 1]

=0
= 1/3
= 1/2

10 12
10 13
10 14
109
10 11

EMAX = 1×1013 eV
EMAX = 1×1014 eV
EMAX = 3×1015 eV

10 12
10 13
10 14

10 11
10 12
10 13
10 14

109

= 1/3
= 0.7
= 1.1
= 1.5

1010 1011 1012 1013 1014 1015
Energy [eV]

D0 = 3×1026
D0 = 3×1027
D0 = 3×1028
D0 = 3×1029

10 13
10 14
109
10 11

1010 1011 1012 1013 1014 1015
Energy [eV](d) R2 -

=0
= 1/3
= 1/2

10 12
10 13
10 14
109
10 11

1010 1011 1012 1013 1014 1015
Energy [eV](f) R2 - Emax

EMAX = 1×1013 eV
EMAX = 1×1014 eV
EMAX = 3×1015 eV

10 12
10 13
10 14

1010 1011 1012 1013 1014 1015
Energy [eV](g) R1 - (D0 = 3×1026)(b) R2 - D0

10 121010 1011 1012 1013 1014 1015
Energy [eV](e) R1 - Emax

109

10 11

1010 1011 1012 1013 1014 1015
Energy [eV]
E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

109

E2dN/dE [erg cm 2s 1]

D0 = 3×1026
D0 = 3×1027
D0 = 3×1028
D0 = 3×1029

E2dN/dE [erg cm 2s 1](a) R1 - D0

E2dN/dE [erg cm 2s 1]

10 11

E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
96

109
10 11

1010 1011 1012 1013 1014 1015
Energy [eV](h) R2 - (D0 = 3×1026)

10 12

= 1/3
= 0.7
= 1.1
= 1.5

10 13
10 14
109

1010 1011 1012 1013 1014 1015
Energy [eV]

Figure 5.13: Input parameter dependency of the model. In panels (a–f), the left (a, c, e, g) andright (b, d, f, h) panels show the results about R1 and R2, respectively. (a, b) Dependency on
D0 (diffusion coefficient). Dotted, dashed, solid, and dot-dashed curves indicate the models with
D0 of 3 × 1026 cm2 s−1 , 3 × 1027 cm2 s−1 , 3 × 1028 cm2 s−1 , and 3 × 1029 cm2 s−1 , respectively.
(c, d) Dependency on δ (degree of dependence of D0 on the energy). Dashed, solid, and dotdashed curves indicate the models with δ of 0, 1/3, and 1/2, respectively. (e, f ) Dependency on
Emax (maximum energy of the accelerated particles in the Sedov phase). Dashed, dot-dashed,and solid curves indicate the models with Emax of 1 × 1013 eV, 1 × 1014 eV, and 3 × 1015 eV,
respectively. (g, h) Dependency on δ with D0 of 3 × 1026 cm2 s−1 . Solid, dot-dashed, dashed,
and dotted curves indicate the models with δ of 1/3, 0.7, 1.1, and 1.5, respectively.

CHAPTER 5. γ-RAY OBSERVATIONS OF SNR HB9 TOWARD THE MEASUREMENT OF
THE DSA HISTORY
97

E2dN/dE [erg cm 2s 1]

10 10

HB9 shell

pesc = 2.0
pesc = 2.4
pesc = 2.7

10 11

10 10

R1

pesc = 2.0
pesc = 2.4
pesc = 2.7

10 11

10 10
10 11

10 12

10 12

10 12

10 13

10 13

10 13

10 14

10 14

10 14

10 15

10 4 100 104 108 1012
Energy [eV]

10 15

10 4 100 104 108 1012
Energy [eV]

R2

10 15

pesc = 2.0
pesc = 2.4
pesc = 2.7

10 4 100 104 108 1012
Energy [eV]10 11(a) pesc = 2.0

10 12
10 13
10 14
109

E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

Figure 5.14: Energy spectra with the hadronic-dominated model under different assumptions
on pesc . The dashed, solid, and dotted lines represent the model emissions with pesc of 2.0,
2.4, and 2.7, respectively. The model with pesc = 2.4 is the same as in Figure 5.12, while forthe models with pesc = 2.0 and 2.7, the acceleration efficiency is given to be η = 0.3 and 0.05,
respectively.

10 11
10 12

10 14

1010 1011 1012 1013 1014 1015
Energy [eV]

109

10 13
109

1010 1011 1012 1013 1014 1015
Energy [eV](d) pesc = 2.7 (D0 = 3×1026)
E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

10 121010 1011 1012 1013 1014 1015
Energy [eV]

D0 = 3×1026
D0 = 3×1027
D0 = 3×1028
D0 = 3×1029

10 13(c) pesc = 2.0 (D0 = 3×1026)
10 11(b) pesc = 2.7

10 11

= 1/3
= 0.7
= 1.1
= 1.5

10 12
10 13
109

1010 1011 1012 1013 1014 1015
Energy [eV]

Figure 5.15: Input parameter dependency of the model in the R2 region when pesc = 2.0(left) and 2.7 (right) are assumed. (a, b) Dependency on D0 (diffusion coefficient). (c, d)Dependency on δ fixing at D0 = 3 × 1026 cm2 s−1 .

Chapter 6

Discussion
With the delayed γ-ray emissions introduced in the previous chapter, we try to model the two
SNRs: SNR G106.3+2.7 (Sect. 6.1) and G335.2+0.1 (Sect. 6.2). The modeling results in bothcases will be discussed by comparing the results of HB9. In Sect. 6.3, we calculate the luminosityof CR escaped from SNRs using the parameters in the delayed γ-ray model and then compare
the measurement at Earth. This provides one test of whether the SNR paradigm is still valid.
Sect. 6.4 describes future prospects for γ-ray observations as well as neutrino observations.

6.1

Application of the delayed γ-ray model to SNR G106.3+2.7

The interpretation for the non-thermal emission from the tail region in SNR G106.3+2.7, as
mentioned in Sect. 4.2.2, proposes that the γ rays come from the molecular cloud region, while
the radio and X-ray emissions come from the SNR shell. We attempt to model the energy
spectrum of the tail region derived in Sect. 3.3.2 again using the delayed gamma-ray model
developed in Chapter 5, instead of the simple assumption used in Chapter 4. In this modeling,
we expect that the tail spectrum is explained by the sum of emissions from the two components,
the shell and cloud regions.
Fig. 6.1 shows the modeling result, and the given parameters are listed in Table 6.1. Combining the emissions from the two components mentioned above can explain all the observed
data. Despite the expectation that the synchrotron emission from the shell part explains theradio and X-ray fluxes (see Fig. 4.2), the new result indicates that the leptonic emissions fromthe shell part reproduce the radio and GeV γ-ray fluxes while the emissions from the cloud
explain the X-ray and TeV γ-ray data. Fig. 6.2 shows the sketch of the interpretation.
In addition, we find that the obtained model parameters are similar to the HB9 results except
for the diffusion coefficient as shown in Table 6.1. Fig. 6.3 shows the difference in the model
curve due to the assumption for (left) the diffusion coefficient and (right) the maximum energy
at the beginning of the Sedov phase. The model reproducing the γ-ray spectrum at the tailrequires a two orders of magnitude smaller diffusion coefficient (2 × 1026 cm2 s−1 ) than those of
the Galactic mean and HB9, as well as a cutoff energy of a few PeV.
Such the small diffusion coefficient can be explained with the effect of the CR-self confine-98

E2dN/dE [erg cm 2s 1]

CHAPTER 6. DISCUSSION

99

Shell + Cloud
Shell
Cloud

10 11
10 12
10 13
10 14

10 5

10 2

101

104 107
Energy [eV]

1010

1013

Figure 6.1: Modeling result to the tail spectrum of SNR G106.3+2.7 with the delayed γ-ray
emission (red line) as well as the shell emission (blue line). The black line indicates the sum of
the non-thermal emissions from the shell and cloud regions. The data points are the same as inFig. 4.1.

B = 4 µG

Shell

Synchrotron Radioep

Enow = 300 GeVtage = 3.9 kyrpeInverse Compton GeV γ-ray
Diffusion time~ 3.1 kyrObserverEmax = 3 PeV at tsedov

Synchrotron
coolingpp→π0p pπ0-decay TeV γ-ray

BISM = 3 µG
dcl = 15 pce
Bcl = 3 µG

Synchrotron X-ray

Cloud
Figure 6.2: Schematic view of the new interpretation for the tail emissions in SNR G106.3+2.7.
Some model parameters in Table 6.1 are also denoted.

100

D0 = 5×1025
D0 = 2×1026
D0 = 5×1026
D0 = 1×1027
D0 = 5×1027
D0 = 1×1028

10 11
10 12

E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

CHAPTER 6. DISCUSSION

10 13
10 14 9
10

1010

1011

1012
1013
Energy [eV]

1014

1015

10 11
10 12

Emax = 10 PeV
Emax = 5 PeV
Emax = 3 PeV
Emax = 1 PeV
Emax = 500 TeV
Emax = 100 TeV

10 13
10 14 9
10

1010

1011

1012
1013
Energy [eV]

1014

1015

Figure 6.3: D0 (left) and Emax (right) dependency of the delayed γ-ray spectrum at the SNR
G106.3+2.7-tail. The data points are the same as in 4.1. In the left panel, the values in thelegend indicate the diffusion constant at 10 GeV with units of cm2 s−1 .
ment (Wentzel, 1974; Fujita et al., 2011; D’Angelo et al., 2018). However, the difference in theobtained diffusion coefficients between SNR G106.3+2.7 and HB9 cannot be easily explainedeven by considering this effect. The diffusion coefficient is expected to depend on the strength of
the surrounding magnetic field (Gabici et al., 2009; Fujita et al., 2009) and its direction (Malkov
et al., 2013). In the latter case, the diffusion of CRs is considered to be asymmetric comparedto the former. If we obtain the delayed gamma-ray spectra from two or more clouds aroundan SNR, we can investigate the isotropy of the CR diffusion to resolve the two scenarios above.
Note that a large distance between the shell and the cloud may be assumed instead of the small
diffusion coefficient but is not preferred here because the case yields a too large flux at the shell
region compared to the cloud region.

6.2

Delayed γ-ray candidate: HESS J1626−490

HESS J1626−490 is one of the undefined TeV sources discovered by the HESS telescope. This
source is located close to SNR G335.2+0.1 but does not overlap with any radio and X-ray
sources in line of sight (Aharonian, 2008b). On contrary, Fermi-LAT observations revealed a
GeV emission spatially coinciding with SNR G335.2+0.1, and two other point sources near the
HESS emission (Abdollahi et al., 2020). Eger et al. (2011) have analyzed the HI and CO line
data as well as the X-ray data and then found the spatially coincidence between the HESS source
and the HI cloud associated with SNR G335.2+0.1, as shown in Fig. 6.4. The authors therefore
suggested that the TeV γ-ray emission originates in the protons escaped from SNR G335.2+0.1.
In this work we first attempt to model this unidentified TeV source with the delayed γ-ray
emission.
We briefly describe here the physical quantity of SNR G335.2+0.1 estimated by previous
studies. The angular diameter of this SNR is ∼ 0.33◦ (Whiteoak & Green, 1996), while thespectral index is found to be α = −0.46 (Clark et al., 1975). The distance from Earth isestimated to be about 1.8 kpc based on the HI observation (Eger et al., 2011). Therefore, thedynamical age is estimated to be ∼ 1.1 kyr, indicating this SNR is young. Although the location

CHAPTER 6. DISCUSSION

101

Table 6.1: Model parameters reproducing the observed spectra of MAGIC-tail in SNR
G106.3+2.7, HESS J1626−490, and HB9. The parameters of HB9 was derived in Chapter 5.
SNR parameters
SN explosion energy (1050 erg)
Initial shock velocity (109 cm s−1 )Age of the SNR, tage (yr)Distance to the SNR (kpc)
Acceleration efficiency
Electron to proton flux ratio
Maximum CR energy at tSedov (eV)Maximum CR energy at tage , Emax (eV)Temporal decay index of Ecut , αMagnetic field in the SNR (µG)
Particle index in the SNR
Particle index after escaping
ISM parameters
Number density (cm−3 )Diffusion coefficient at 10 GeV, D0 (cm2 s−1 )Energy slope of diffusion
Magnetic field in ISM (µG)
Cloud parameters
Distance to the Cloud from SNR (pc)
Radius of the Cloud (pc)
Average hydrogen number density (cm−3 )Reflected SNR age (yr)Magnetic field in Cloud (µG)

G106.3+2.7–Tail
1.0
1.0
3.9 × 103
0.8
0.1
0.02
3 × 1015
3 × 1011
3.15
4
2.0
2.0
G106.3+2.7–Tail
0.1
2 × 1026
1/3
3
G106.3+2.7–Tail
15.0
2.2
300
−3.1 × 103
3

G335.2+0.1
5.0
1.0
1.6 × 103
1.8
0.1
0.01
1 × 1015
5 × 1010
4.43
8
2.0
2.4
G335.2+0.1
1.0
1 × 1027
1/3
3
G335.2+0.1
13.6
7.5
200
−8.9 × 102
3

HB9–R23.0
1.0
6.6 × 103
0.8
0.1
0.02
3 × 1015
3 × 1011
3.16
8
2.0
2.4
HB9–R20.06
3 × 1028
1/3
3
HB9–R239.4
7.0
200−6.3 × 102
3

CHAPTER 6. DISCUSSION

102

1600
Radio continuumat 873 MHz

1500
TeV γ-ray
1300
1200
1100
1000

12CO (J=1-0)

900

HI brightness temperature [K]

1400

800
700

Figure 6.4: HI line intensity map in the vicinity of SNR G335.2+0.1 (McClure-Griffiths et al.,2005), which are integrated over a velocity range from −31 km s−1 to −18 km s−1 based on
Eger et al. (2011). The green, cyan, and magenta contours represent the emissions of radio continuum at 843 MHz (http://snrcat.physics.umanitoba.ca/downloadFITS.php?location=
FITS/radio/g335.2+00.1_8.73E08_MOST.fits), 12 CO (J = 1 − 0) line (Dame et al., 2001),
and TeV γ rays (Abdalla et al., 2018b), respectively. The black markers indicate the positionsof 4FGL sources (Abdollahi et al., 2020).

CHAPTER 6. DISCUSSION

Shell: G335.2+0.1

10 10

E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

10 10

103

10 11
10 12
10 13
10 14

10 7 10 4 10 1

Synchrotron

102

105

Energy [eV]
IC

108

1011

1014

Bremsstrahlung

Cloud: HESS J1626-490

10 11
10 12
10 13
10 14

10 7 10 4 10 1

Hadronic

102

105

Energy [eV]

Total

108

1011

1014

Hadronic from BG

Figure 6.5: Modeling result to the energy spectra of the shell (left) and cloud regions (right)
around SNR G335.2+0.1. The blue, green, magenta, and red lines represent the synchrotron,IC, electron bremsstrahlung, and π 0 -decay emissions, respectively. The black line indicates thesum of the non-thermal emissions from the shell and cloud regions. Also shown with the grey
shaded area is the hadronic emission from the background CRs. In the left panel, the data
points in the radio and GeV band are taken from the literatures (Clark et al., 1975; Whiteoak& Green, 1996; Green, 1988) and Abdollahi et al. (2020), respectively. In the right panel, theflux upper limit in the X-ray band is taken from Eger et al. (2011), while the TeV spectrum is
adopted from Aharonian (2008b).
of a pulsar PSR 1627−4845 is spatially consistent with SNR G335.2+0.1, its characteristic age(2.7 Myr) and estimated distance (6.8 kpc) significantly differ from those of the SNR mentioned
above. Therefore, this pulsar is not considered to be associated with the SNR (Kaspi et al.,1996).
Here, we model the observed spectra at the SNR G335.2+0.1 (shell) and HESS J1626-490(cloud), simultaneously, as with in Chapter 5. The radio data for the shell region are taken from
Clark et al. (1975), Whiteoak & Green (1996), and Green (1988). The Fermi-LAT spectrum of
4FGL J1626−4848c (Abdollahi et al., 2020) is adopted as the shell emission. We also adopt theupper limit on the X-ray flux between 0.7 and 10.0 keV derived by Eger et al. (2011) as well as the
TeV spectrum (Aharonian, 2008b), for the non-thermal emission from the HESS source (cloud)
region. The energy density of seed photon for IC is estimated to be 0.25 eV cm−3 at 40 K using
the model parameter in GALPROP (Porter et al., 2008). The cloud parameters (density, size,and location) for the HESS source region are adopted from the CO and HI observations (Egeret al., 2011).
Fig. 6.5 shows the modeling result to the energy spectra of G335.2+0.1 and HESS J1626−490.
The model parameters used in the modeling are tabulated in Tab. 6.1. We find the hadronic
emissions can reproduce both γ-ray spectra at the shell and cloud. Fig. 6.6 shows the model
variations with diffusion coefficient and cutoff energy. The diffusion coefficient to reproduce theobserved spectrum is found to be D ∼ 1027 cm2 s−1 , which is about one order smaller thanthe Galactic mean. As discussed in Sect. 6.1, the small diffusion constant can be explainedby the effect of the CR-self confinement (e.g., Wentzel, 1974). The maximum energy at tsedov

104
D0 = 1×1026 cm2/s
D0 = 1×1027 cm2/s
D0 = 1×1028 cm2/s
HESS J1626-490

10 10
10 11
10 12
10 13 810E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

CHAPTER 6. DISCUSSION

109 1010 1011 1012 1013 1014
Energy [eV]

Emax = 10 TeV
Emax = 100 TeV
Emax = 1 PeVHESS J1626-490

10 10
10 11
10 12
10 13 810

109 1010 1011 1012 1013 1014
Energy [eV]

Figure 6.6: D0 (left) and Emax (right) dependency of the delayed γ-ray model for HESS
J1626−490. The data points are taken from Aharonian (2008b).
is estimated to be ≳ 100 TeV. Although further observations with the Fermi-LAT or HESS
telescopes are needed to better identify the counterpart of the γ-ray emissions, such observations
would provide strong evidence that this SNR was a PeVatron in the past.

6.3

Discussion on the CR luminosity

The modeling results suggest that all three SNRs examined in this chapter are still allowed toaccelerate protons up to ≳ 1 PeV. Fig. 6.7 shows the time evolution of the maximum particle
energy for each SNR based on Table 6.1. The obtained temporal decay index of the maximum
energy (α = 3–4) is in good agreement with the phenomenological prediction (Ohira et al.,2011).
Here, we calculate intensities of the PeV CRs escaped from those SNRs and compare withthe one observed at Earth. With the CR energy spectrum (Dermer, 1986)
JCR = 2.2(E/GeV)−2.75 cm−2 s−1 GeV−1 sr−1 ,(6.1)

the CR luminosity at Earth can be approximated as follows:
LCR (E) ∼WCR
4π VGal 2∼
E JCR ,
tCR
c tCR

(6.2)

where WCR is the total energy of CRs, VGal is the volume of the Galaxy, and tCR is a diffusion
timescale for CRs. Assuming that the shape of Galaxy is a cylinder, VGal can be approximatedas follows

2 

R
H
2
68
3
,
(6.3)
VGal ∼ 2πR H = 2.09 × 10 cm
15 kpc
7 kpc
where H is the halo size of Galaxy and R is the radius of the Galaxy disk. The diffusion timescaleof CRs from the Galaxy can be estimated bytCR (E) ∼H2∼ 1.15 × 1016 s
D(E)



H
7 kpc

2 

E
1 GeV

−δ
,

(6.4)

CHAPTER 6. DISCUSSION

105

Maximum particle energy [eV]

1016

HB9
G106.3+2.7
G335.2+0.1

1015
1014
1013
1012
1011
1010
109

Phenomenological prediction [Ohira+2011]
Systematic study [Suzuki+2022]

102

103
SNR age [yr]104

Figure 6.7: Maximum energy of accelerated particles in SNRs as a function of the age. The data
points and solid lines represent the results obtained in this work. The circle and square markers
indicate the maximum energies at the shell and cloud region, respectively. The star markers
indicate the estimated time point (tSedov ) at entering the Sedov period. Although the modelin this work does not consider the younger part than tSedov , an extrapolation with tage basedon the theoretical expectation (e.g., Ohira et al., 2011) is shown for the sake of reference. The
shaded area enclosed by the red dotted lines represent the phenomenological prediction basedon Ohira et al. (2011), while the black dashed line represents the result derived by Suzuki et al.
(2022).

CHAPTER 6. DISCUSSION

106

Table 6.2: CR luminosity from SNRs computed with the model parametersof each SNR in Tab. 6.1.
Name
HB9
G106.3+2.7
G335.2+0.1

LSNR (10 TeV)
[erg s−1 ]2.88 × 1039
1.03 × 1039
2.33 × 1039

LSNR /LCR
at 10 TeV
1.08
0.39
0.87

LSNR (1 PeV)[erg s−1 ]4.56 × 1038
1.03 × 1039
3.69 × 1038

LSNR /LCR
at 1 PeV
0.54
1.22
0.44

where D is the diffusion coefficient and adopted the results of the Boron-to-Carbon ratio (Yuanet al., 2017). Substituting Eq. 6.3 and 6.4 into Equation 6.2 yields
40

LCR (E) ∼ 2.67 × 10

erg s

−1



E
1 GeV

−0.75+δ 

VGal
2.1 × 1068 cm3


,

(6.5)

where δ is the slope of the diffusion coefficient. In this calculation, we adopt the value (δ =0.500) estimated by Yuan et al. (2017). As a result, the CR luminosity is estimated to be
LCR = 2.67 × 1039 erg s−1 at 10 TeV and 8.45 × 1038 erg s−1 at 1 PeV.
On the other hand, the luminosity of CRs escaped from SNRs can calculated as follows.
First, the energy distribution of CRs escaped from an SNR is described by Eq. 5.8:

−1 


−pesc η(2−pesc )ESN × Emax 2−pesc − 1
E×pesc ̸= 2
2
Enow
Enow
Enow
Nesc (E) =
(5.8)

ηESN
−2

E
p
=
2
esc
ln(Emax /Enow )The total luminosity from SNRs can be approximated by multiplying the birth rate of SNe:LSNR ∼ Nesc (E) × E 2 × (SN rate).

(6.6)

As for the SN rate, 0.03 yr−1 (Tammann et al., 1994) is adopted.
Using the best-fit parameters shown in Table 6.1, we calculate the contribution to the CR
luminosity and summarize in Table 6.2. The calculated luminosity from SNRs (LSNR ) varies by
a factor of ∼ 3 but is roughly consistent with the measurement at Earth (LCR ).

6.4

Future prospect

6.4.1

Future γ-ray observations

As mentioned in Sect. 5.4.3, our work has suggested that further TeV γ-ray observations would
give a crucial constraint on the model parameters. In particular, the parameter Emax shouldbe estimated better to elucidate whether the SNR is a PeVatron. We simulate energy spectra at HB9-R2 expected to be obtained with the next generation γ-ray observatories, CTANorth (Bernlöhr et al., 2013) and LHAASO (Addazi et al., 2022). Since HB9-R2 is more extended (θsrc = 0.36◦ ) than PSFs of these instruments (∼ 0.05◦ for CTA-Northand ∼ 0.30◦
q22 )/θ 2 , to+ θsrcfor LHAASO), the sensitivities for a point source are scaled down by (θPSF
PSF

CHAPTER 6. DISCUSSION

E2dN/dE [erg cm 2s 1]

Model
Fermi-LAT data

CTA-North sensitivity (50 hr)
Expectation by CTA-North

10 11

10 11

10 12

10 12

10 13

10 13

10 14

(a) D0 = 3×1028, Emax = 3 PeV
109

E2dN/dE [erg cm 2s 1]

107

10 14

1010 1011 1012 1013 1014 1015
10 11

10 12

10 12

10 13

10 13

10 14

(c) D0 = 3×1027, Emax = 3 PeV
109

1010 1011 1012 1013 1014 1015
Energy [eV](b) D0 = 3×1028, Emax = 10 TeV
109

10 11

LHAASO sensitivity (1 yr)
Expectation by LHAASO

10 141010 1011 1012 1013 1014 1015(d) D0 = 3×1027, Emax = 10 TeV
1091010 1011 1012 1013 1014 1015
Energy [eV]

Figure 6.8: Energy spectra of HB9–R2 expected to be observed with the CTA (open circle) and
LHAASO experiments (open diamond). The red dashed line shows the 5σ-detection sensitivityof the 50-hr CTA and 1-yr LHAASO observations. Black solid line shows the model emissions
calculated using the parameters in Table 5.2 except for D0 and Emax , which are represented ineach panel. The black points represent the Fermi-LAT data dervied in Sect. 5.3.2.
account for the deterioration of the signal-to-background ratio in the extended source analysis.
Fig. 6.8 shows the expected energy spectra from a 50-hour observation by CTA-North and a
one-year observation by LHAASO. Although four sets of model parameters that have not been
ruled out in our study using the Fermi-LAT data are used for the simulation, the differences inthe expected spectra can be found.
It should be noted that the model spectra at the TeV band may be affected by the assumption
when electrons escape from the shock region. Under the assumption of Chapters 5 and 6, where
the electron and proton begin to escape simultaneously, the amount of escaped electrons is
maximized. Appendix B estimates the impact due to this assumption is at most a factor of ∼ 2.
Once Emax is determined better, the temporal decay index (α) in the time evolution can be
estimated with fewer uncertainties. The presently obtained value is α ∼ 3.1–4.4 for the three
SNRs, which is steeper than that (α = 0.81 ± 0.24) obtained in the previous study with thesystematic investigation of the gamma-ray spectra in many SNRs (Suzuki et al., 2022). As aresult, our result supports that these SNRs were a PeVatron. While the future observations for
HB9 in the northern sky are given here, future projects such as CTA-South (Bernlöhr et al.,2013), ALPACA (Asaba et al., 2018), and SWGO (Albert et al., 2019) are underway to observe γray sources with better sensitivities in the southern sky, such as SNR G335.2+0.1. Observations
with these future experiments would provide more precise results for the time evolution of DSA

CHAPTER 6. DISCUSSION

108

in SNRs, which may find the difference between the types of SNe (e.g., Kamijima & Ohira, 2021,2022).

6.4.2

Neutrino observations

Detection of neutrinos can strengthen the possibility of proton acceleration in SNRs. The
IceCube experiment has achieved the highest sensitivity in modern neutrino observations, andits upgrade project, IceCube-Gen2, is also underway to improve the sensitivity by a factor offive or more (Aartsen et al., 2021b). Here we estimate the neutrino flux at the tail region of
SNR G106.3+2.7 and compare it with the sensitivities of the IceCube and its future project.
Protons in SNRs generate π + , π − , and π 0 via p-p collision with a ratio of 1:1:1, as in Eq. 1.41.
As in Eq. 1.42, π 0 then decays to 2γ, while π ± decays to µ and νµ . Furthermore, µ decays to e±and νe . According to Ahlers & Murase (2014), the γ-ray flux and the neutrino fluxes for eachflavor produced in this process is expressed as
Eγ Jγ (Eγ ) ≃ exp(−d/λγγ ) ×

1X
Eν Jνα (Eν ),
3 ν

(6.7)

α

where Eγ (Eν ) is the γ-ray (neutrino) energy, Jγ (Jν ) is the differential flux for γ-rays (neutrinos),λγγ is the optical depth for the γγ absorption, and d is the distance from Earth. As for the effectof γγ absorption, we assume that it is negligible because d (∼ 1 kpc) for SNR G106.3+2.7 isenough small. Since the energy of neutrinos is about one quarter of pions, Eγ ≃ 2Eν . At Earth,after the propagation, the neutrino fluxes for each flavor are almost the same due to neutrino
oscillations (Cavasinni et al., 2006). Hence, Eq. 6.7 can be approximated as
Eγ Jγ (Eγ ) ≃ 2Eγ Jν (2Eγ ).

(6.8)

Using this relation between the γ-ray and neutrino flux, we estimate the energy spectrum ofmuon neutrino at the tail region of SNR G106.3+2.7 and show it in Fig. 6.9. The expected flux
is lower than the sensitivities of the neutrino observations and is also consistent with the upper
limit recently obtained by the IceCube experiment (Abbasi et al., 2022). Although 5σ detectionof the neutrino emission from this region would be challenging even with IceCube-Gen2, Fig. 6.9
suggests the experiment may find the excess with 2–3σ significance at the range of 10–100 TeV.

CHAPTER 6. DISCUSSION

109

Hadronic-dominated model for tailE2dN/dE [erg cm 2s 1]

IceCube (10 yr)

10 11

IceCube-Gen2 (10 yr)
Hadronic γ-ray
Neutrino

10 12
10 13
10 14109

1010

1011

1012 1013
Energy [eV]

1014

1015

Figure 6.9: Expected energy spectra of muon neutrinos (Fνµ + Fν µ ; red solid) as well as the
hadronic γ rays (red dotted) at the tail region of SNR G106.3+2.7. The blue solid (dashed) line
represents the 5σ differential sensitivities of 10-yr IceCube (IceCube-Gen2) observations at the
declination angle of 30◦ (Aartsen et al., 2021b).

Chapter 7

Summary
To elucidate the origin of CRs up to PeV, we investigated whether SNR PeVatrons (which accelerate particles up to PeV) exist in our Galaxy and whether SNRs can supply the total amount of
PeV CRs measured at Earth, using γ-ray observations. For this purpose, we attempted to validate the time evolution scenario, in which the maximum acceleration energy in SNRs depends
on the age and reaches PeV only for a limited period (tage < 1 kyr). Although theoretical works
have expected this scenario, there was no conclusive observational evidence. We measured thetime evolution of the maximum energy in a single SNR, for the first time, by simultaneously
observing γ-ray spectra from the SNR shell (reflecting the present particle distribution) andthe nearby clouds (reflecting the past particle distribution). This newly proposed method wasapplied to γ-ray observation data for the three SNRs, HB9, G106.3+2.7, and G335.2+0.1.
For the first application of this time-evolution measurement, we analyzed the GeV γ-ray
data in the vicinity of SNR HB9, using 12-yr Fermi-LAT observations. We newly detected theγ-ray emissions spatially correlated with the cloud regions with a statistical significance of 6.1σ
as well as the shell emission, the latter of which is consistent with previous studies. The γ-ray
spectra at the cloud regions can be fitted with a power-law function with an index of < 2.0,
which is the hardest spectrum ever observed at clouds around SNRs. The cloud spectra weresuccessfully reproduced with the model of the time evolution scenario, in which CRs escapedfrom the shell in the past and now illuminate the cloud and generate the delayed γ-ray emissions.
The results show that the SNR accelerated CRs up to higher energies (∼ 3 PeV) in the pastthan the maximum energy at present (300 GeV).
Expecting that the SNR PeVatron may still be hidden in the 100-TeV sources detectedby the LHAASO experiments, we focused on LHAASO J2226+6057/SNR G106.3+2.7 in thisthesis. Since this region is a PWN/SNR complex, the origin of 100-TeV γ-ray emission fromthis source was not yet identified due to the limited angular resolution in previous studies.
We observed it with the MAGIC telescopes, which achieve better angular resolution (0.07–
0.1◦ ) than the previous γ-ray observations for this source. We then found the extended TeV
emissions (> 6σ) spatially coinciding with the radio shell, which consists of the head and tailregions. In particular, we detected the emission above 6 TeV only in the tail region but not in
the head region, the latter of which contains the PWN. This fact suggests that the 100-TeV γ

110

CHAPTER 7. SUMMARY

111

rays originate in the SNR rather than the PWN. The spectral modeling in the two regions were
performed. Both leptonic and hadronic emissions can reproduce the observed spectrum of thehead region, while the leptonic model failed to reproduce the tail spectrum under the assumptionthat the observed emissions above 10 TeV come only from the tail region. We thus suggested
that the γ-ray spectrum up to 100 TeV in the tail region originates in proton acceleration up
to PeV in the SNR. Furthermore, we found that the time evolution scenario above explains theinconsistency between the proton cutoff energy of ∼ PeV and the SNR age of ∼ 4 kyr.
In addition, we performed the spectral modeling for HESS J1626−490/G335.2+0.1 with thetime evolution scenario. The energy spectrum of the unidentified TeV source, HESS J1626−490,
could be explained with the delayed γ-ray emissions by protons escaped from SNR G335.2+0.1.
The modeling results of the three SNRs yielded similar SNR parameters except for the
diffusion coefficient, indicating that the time evolution scenario is general for SNRs. Moreover,the proton acceleration up to ∼ PeV is suggested in the all three SNRs. The amount of theescaped protons from SNRs calculated with these model parameters is in agreement with the CR
flux observed at Earth within a factor of ∼ 3. The results give compelling evidence supportingthat Galactic SNRs are the origin of CRs up to PeV.
We also found this simultaneous model fitting to the shell and cloud regions allows for
estimating the diffusion coefficient around SNRs with unprecedented accuracy. The diffusion
coefficients around SNRs were theoretically expected to be smaller than the Galactic mean, but
the results for each SNR varied by about two orders of magnitude, indicating a significant environmental dependence. Further observations with the future project were expected to estimate
the above model parameters better. The luminosity frontier in TeV-PeV γ-ray observations will
provide a unified picture of CR propagation and acceleration in SNRs, leading to the discoveryof the CR origin.

Acknowledgements
First of all, I am deeply grateful to my supervisor, Prof. Hidetoshi Kubo, for his advice andcontinuous support. I would also like to express my gratitude to my research collaborators.
Dr. Takayuki Saito and Dr. Marcel Strzys analyzed the MAGIC data for SNR G106.3+2.7
with me. Advice and comments given by the MAGIC Collaboration members, in particular Dr.
Daniel Mazin and Dr. Marc Ribo, have been a great help in the MAGIC analysis. The MAGIC
Collaboration is also acknowledged for allowing the use of proprietary observational data. Some
MAGIC results presented in this thesis have not undergone the internal review procedure of
MAGIC publications. Dr. Hidetoshi Sano, Dr. Tsuyoshi Inoue, and Dr. Yasuo Fukui havehelped with the radio data analysis and its interpretation. Dr. Takaaki Tanaka gives insightfulcomments and suggestions for the Fermi-LAT data analysis. Dr. Wataru Ishizaki has provided
the model calculation code and discussed with me the application to the observation data of
SNR HB9. Special thanks to Dr. Toshihiro Fujii, Dr. Seiya Nozaki, and Kenta Terauchi for the
fruitful discussions. Finally, I appreciate the support and encouragement of all the people in the
Cosmic ray group at Kyoto University.

112Appendices113

Appendix A

Gas density in the vicinity of SNR
G106.3+2.7
We calculate the gas density in the two regions of SNR G106.3+2.7 with the following outline.
We use the data of HI line measured with the Dominion Radio Astronomy Observatory (DRAO)
Synthesis Telescope (Landecker et al., 2000) and 12 CO (J = 1 − 0) line measured with theFive College Radio Astronomy Observatory (FCRAO; Heyer et al., 1998) from the Canadian
Galactic Plane Survey (CGPS; Taylor et al., 2003) database. These observations were carriedout with the velocity resolution of 0.824 km s−1 at HI line and 0.98 km s−1 at CO line. The
R vmaxT (v)dv,
following relationship is used to calculate the column density: NH [cm−2 ] = X vmin
where v is the radial velocity, T (v) is the observed brightness temperature (K) and X is theconversion factor (Dickey & Lockman, 1990). HI-to-NHI and CO-to-NH2 are given by XHI =1.823 × 1018 (Dickey & Lockman, 1990) and XCO = 2.0 × 1020 (Bolatto et al., 2013). Fig. A.1
shows the radial profiles of HI and 12 CO (J = 1 − 0) line. There is a significant velocity
dependence of the column density, especially in the CO data, which is a concern because theuncertainty of the velocity range affects the calculation of the gas density. Here, we consider twocases on the velocity ranges that associates with SNR G106.3+2.7: (i) −7.23 to −5.59 km s−1
suggested by Kothes et al. (2001) and (ii) −6.41 to −3.94 km s−1 suggested by Acciari et al.
(2009); Albert et al. (2020). The clouds associated with the production of the observed γray emission are assumed to be a spherical region around the emission center with a radiusof 800 pc × tan(0.16◦ ) ∼ 2.2 pc estimated from the MAGIC data as shown in Table3.2. The
calculation results are summarized in Table A.1.
Table A.1: Gas densities of hydrogen atoms at the head and tail region. nHI and nCO areestimated with the HI line and 12 CO (J = 1 − 0) line data, respectively.
Velocity range [km s−1 ] −7.23 – −5.59 −6.41 – −3.94
nHI at head [cm−3 ]42
59
nCO at head [cm−3 ]73
66−3nHI at tail [cm ]3855
nCO at tail [cm−3 ]
137
191

114

4
3

150
100

2

50

1

0

14

12

10 8
6
4
Radial velocity [km s 1]

2

0

0

2.5

CO at head
CO at tail

2.0

(i)

115

5

(ii)

H2 Column density [×1020 cm 2]

200

(i)

5

(ii)

CO Brightness temperature [K]HI at head
HI at tail250

HI Column density [×1020 cm 2]

HI Brightness temperature [K]APPENDIX A. GAS DENSITY IN THE VICINITY OF SNR G106.3+2.7

4

1.5

3

1.0

2

0.5

1

0.0

14

12

10 8
6
4
Radial velocity [km s 1]

2

0

0

Figure A.1: The HI (left) and 12 CO (J = 1 − 0) (right) radial profile at the head and tailregion. In both panel, red and blue data represent the profile of the head and tail regions. The
green arrow labeled (i) indicates the velocity range pointed out by Kothes et al. (2001), while
the magenta arrow labeled (ii) indicates the range used in Acciari et al. (2009) and Albert et al.
(2020).
There is not big difference of the results between the integration velocity ranges. We use
100 cm−3 and 200 cm−3 as a gas density of head and tail regions for the modelling.

Appendix B

Electrons in the delayed γ-ray
emission model
B.1

Derivation of the spectrum of CR electrons

The transport equation for CR electrons is written by∂fe,out∂(t, r, E) − DISM (E)∆fe,out (t, r, E) +(P (E)fe,out (t, r, E)) = qe,s (t, r, E).
∂t
∂E

(B.1)

According to Atoyan et al. (1995), the Green function of equation (B.1) is
P (Et−t0 )|r − r0 |2
G(t, r, E; t0 , r0 ; E0 ) = 3/2
exp
−
2
3
rdif
π P (E)rdif!
δ (Et−t0 − E0 ) ,(B.2)

where Et−t0 is defined by the following implicit relation:t − t0 =Z Et−t0

Erdif =
andZ Et−t

0

∆u =
E

√

dE ′
,
P (E ′ )4∆u,

DISM (E ′ ) ′dE .
P (E ′ )

(B.3)
(B.4)

(B.5)

From linearity, the solution to equation (B.1) can be calculated as follows:Z
fe,out (r, t, E) =

Z
dE0

3

Z

d r0

dt0 G(t, r, E; t0 , r0 ; E0 )qe,s (r0 , t0 , E0 ) .

(B.6)

In this work, only synchrotron radiation is considered as the most important process inradiative cooling. Thus, the cooling function P (E) can be written as follows:
4
P (E) = σT c
3



E
me c2

2

116

2
BISM
≡ Qsyn E 2 ,
8π

(B.7)

APPENDIX B. ELECTRONS IN THE DELAYED γ-RAY EMISSION MODEL

117

where σT is the Thomson scattering cross section, me is the electron mass, and BISM is a valueof the magnetic field in the interstellar space. Evaluating equation (B.3) with using equation(B.7), we obtainE
Et−t0 =.
(B.8)
1 − Qsyn E(t − t0 )For equation (B.5), we can calculate using equation (5.5) and (B.7) as follows:"1−δ #2
rdif
E
DISM (E)
1−
.
=
∆u =
4
(1 − δ) Qsyn E
Et−t0

(B.9)

By using equations (5.2), (5.6), (5.12), (B.2), (B.8), and (B.9) and performing the integrationof equation (B.6), the solution to equation (B.1) can be obtained as follows:"!
!#
Kep Nesc (Ec ) Ec2
1(r − Rc )2(r + Rc )2
fe,out (t, r, E) = 3/2
× exp −
− exp −
,
2
2
Rd,e
Rd,e
4π rRc Rd,e E 2 1 − Qsyn tc Ec /αwhere tc = tc (t, E) is the time when the electron whose energy is E at the current time t is
injected and determined as a solution to an algebraic equation tc = tesc (Et−tc ), or equivalently
the below equation,
E= Eesc (tc ),(B.10)
1 − Qsyn E (t − tc )Ec (t, E) = Eesc (tc ), Rc (t, E) ≡ Resc (Ec ), and Rd,e (t, E) is the diffusion length of CR electronsv"u 1−δ #u 4DISM (E)E
Rd,e = t
1−.
(1 − δ) Qsyn E
Ec(B.11)

Here, we have assumed that there is only one solution to equation (B.10), however this is not
generally true. The condition for there to be only one solution can be written as follows:tSedov ≥ tc,crit ≡ tSedov
Qsyn tSedov Emax 1/(α−1)α(B.12)For the range of parameters we use in this work, this condition is satisfied, and then the assumption of only one solution is consequently justified.

B.2

Impact of electron parameter to the delayed γ-ray model

We demonstrate the impact of assumed electron parameters to the γ-ray emission model. According to Ohira et al. (2012), the estimate of ξe (= t/tsedov ) depends on the assumption of theenvironment of the magnetic field. The start time of escape of CR electrons is expressed by:te =9m2e c5/2 RS2
3/2
8ηg,free ηacc eEmax!−2/(α+2αB −1)(α+2α +2)/(α+2αB −1)

tsedov B,(B.13)

10 11(a) R1 - D0

10 12
10 13
10 14
109

e = 1.0
e = 3.22
e = 10.8
e = 16.6

E2dN/dE [erg cm 2s 1]

E2dN/dE [erg cm 2s 1]

APPENDIX B. ELECTRONS IN THE DELAYED γ-RAY EMISSION MODEL

1010 1011 1012 1013 1014 1015
Energy [eV]

10 11(b) R2 - D0

10 12

118
e = 1.0
e = 3.22
e = 10.8
e = 16.6

10 13
10 14
109

1010 1011 1012 1013 1014 1015
Energy [eV]

Figure B.1: ξe dependence of the delayed γ-ray model. Panels a and b show the energy spectraat the R1 and R2 cloud around SNR HB9. The data points are the same as in Fig. 5.13. The
difference of the line type indicates the difference in the assumption of ξe .
where
ηg,free ∼ 1: Gyrofactor during the free-expansion phase (tage < tsedov )ηacc ∼ 10: Numerical factor that depends on the shock compression ratio and the spatial
dependence of the diffusion coefficient
RS : SNR radius when tage = tsedov
α: Temporal decay index of the maximum energy of CR protons
αB : Temporal decay index of the magnetic field
In particular, αB depends on the environment of the magnetic field around SNRs, i.e., theassumption of the gyro-factor (ηg ) and the shock velocity (vsh ) in the SNR. Ohira et al. (2012)
assumed three magnetic field models and derived αB = α − 1, 9/10, and 3/5. Once we substitute
α = 3.16, tsedov = 357 yr, and RS = 3.65 pc for the SNR HB9 (Sect. 5.4), ξe is estimated tobe 3.22 for αB = α − 1, 10.8 for αB = 9/10, and 16.6 for αB = 3/5. Fig. B.1 shows the model
emissions computed with the parameters above. The TeV γ-ray fluxes computed with the larger
ξe are (up to a factor of ∼ 2) smaller than the case of ξe = 1 because of fewer escaped electrons.

Bibliography
M. G. Aartsen et al. (2013). ‘Measurement of the cosmic ray energy spectrum with IceTop-73’.
Phys. Rev. D 88(4):042004.
M. G. Aartsen et al. (2021a). ‘Detection of a particle shower at the Glashow resonance with
IceCube’. Nature 591(7849):220–224. [Erratum: Nature 592, E11 (2021)].
M. G. Aartsen et al. (2021b). ‘IceCube-Gen2: the window to the extreme Universe’. J. Phys.
G 48(6):060501.
R. Abbasi et al. (2022). ‘Searches for Neutrinos from LHAASO ultra-high-energy γ-ray sources
using the IceCube Neutrino Observatory’ .
H. Abdalla et al. (2018a). ‘H.E.S.S. observations of RX J1713.7−3946 with improved angular andspectral resolution: Evidence for gamma-ray emission extending beyond the X-ray emitting
shell’. Astron. Astrophys. 612:A6.
H. Abdalla et al. (2018b). ‘The H.E.S.S. Galactic plane survey’. Astron. Astrophys. 612:A1.
A. Abdo et al. (2007). ‘TeV Gamma-Ray Sources from a Survey of the Galactic Plane with
Milagro’. Astrophys. J. 664:L91–L94.
A. Abdo et al. (2009a). ‘Fermi LAT detection of pulsed gamma-rays from the Vela-like pulsars
PSR J1048-5832 and PSR J2229+6114’. Astrophys. J. 706:1331–1340.
A. A. Abdo et al. (2009b). ‘Milagro Observations of TeV Emission from Galactic Sourcesin the Fermi Bright Source List’. Astrophys. J. 700:L127–L131. [Erratum: Astrophys.
J.703,L185(2009)].
S. Abdollahi et al. (2020). ‘F ermi Large Area Telescope Fourth Source Catalog’. Astrophys. J.
Suppl. 247(1):33.
A. U. Abeysekara et al. (2017). ‘Observation of the Crab Nebula with the HAWC Gamma-Ray
Observatory’. Astrophys. J. 843(1):39.
V. A. Acciari et al. (2009). ‘Detection of Extended VHE Gamma Ray Emission from G106.3+2.7
with VERITAS’. Astrophys. J. 703:L6–L9.
F. Acero et al. (2010). ‘First detection of VHE gamma-rays from SN 1006 by H.E.S.S’. Astron.
Astrophys. 516:A62.
119

BIBLIOGRAPHY

120

F. Acero et al. (2017). ‘Prospects for Cherenkov Telescope Array Observations of the Young
Supernova Remnant RX J1713.7−3946’. Astrophys. J. 840(2):74.
M. Ackermann et al. (2012a). ‘Fermi-LAT Observations of the Diffuse Gamma-Ray Emission:
Implications for Cosmic Rays and the Interstellar Medium’. Astrophys. J. 750:3.
M. Ackermann et al. (2012b). ‘The Fermi Large Area Telescope On Orbit: Event Classification,
Instrument Response Functions, and Calibration’. Astrophys. J. Suppl. 203:4.
M. Ackermann et al. (2013). ‘Detection of the Characteristic Pion-Decay Signature in Supernova
Remnants’. Science 339:807.
A. Addazi et al. (2022). ‘The Large High Altitude Air Shower Observatory (LHAASO) Science
Book (2021 Edition)’. Chin. Phys. C 46:035001–035007.
M. Aguilar et al. (2015a). ‘Precision Measurement of the Helium Flux in Primary Cosmic Raysof Rigidities 1.9 GV to 3 TV with the Alpha Magnetic Spectrometer on the International
Space Station’. Phys. Rev. Lett. 115(21):211101.
M. Aguilar et al. (2015b). ‘Precision Measurement of the Proton Flux in Primary Cosmic Rays
from Rigidity 1 GV to 1.8 TV with the Alpha Magnetic Spectrometer on the International
Space Station’. Phys. Rev. Lett. 114:171103.
M. Aguilar et al. (2016). ‘Precision Measurement of the Boron to Carbon Flux Ratio in Cosmic
Rays from 1.9 GV to 2.6 TV with the Alpha Magnetic Spectrometer on the International
Space Station’. Phys. Rev. Lett. 117(23):231102.
M. Aguilar et al. (2018). ‘Observation of Complex Time Structures in the Cosmic-Ray Electron and Positron Fluxes with the Alpha Magnetic Spectrometer on the International Space
Station’. Phys. Rev. Lett. 121(5):051102.
F. Aharonian (2008a). ‘Discovery of very high energy gamma-ray emission coincident withmolecular clouds in the W28 (G6.4-0.1) field’. Astron. Astrophys. 481:401.
F. Aharonian (2008b). ‘HESS VHE Gamma-Ray Sources Without Identified Counterparts’.
Astron. Astrophys. 477:353–363.
F. Aharonian et al. (2021). ‘Calibration of the air shower energy scale of the water and air
Cherenkov techniques in the LHAASO experiment’. Phys. Rev. D 104(6):062007.
F. A. Aharonian & A. M. Atoyan (1996). ‘On the emissivity of π 0̂-̂decay gamma radiation in thevicinity of accelerators of galactic cosmic rays.’. Astronomy and Astrophysics 309:917–928.
M. Ahlers & K. Murase (2014). ‘Probing the Galactic Origin of the IceCube Excess withGamma-Rays’. Phys. Rev. D 90(2):023010.
H. S. Ahn et al. (2010). ‘Discrepant hardening observed in cosmic-ray elemental spectra’. Astrophys. J. Lett. 714:L89–L93.

BIBLIOGRAPHY

121

H. Akaike (1974). ‘A New Look at the Statistical Model Identification’. IEEE Transactions on
Automatic Control 19:716–723.
A. Albert et al. (2019). ‘Science Case for a Wide Field-of-View Very-High-Energy Gamma-Ray
Observatory in the Southern Hemisphere’ .
A. Albert et al. (2020). ‘HAWC J2227+610 and its association with G106.3+2.7, a new potential
Galactic PeVatron’. Astrophys. J. Lett. 896:L29.
J. Albert et al. (2007). ‘Unfolding of differential energy spectra in the MAGIC experiment’.
Nucl. Instrum. Meth. A 583:494–506.
J. Albert et al. (2008a). ‘FADC signal reconstruction for the MAGIC Telescope’. Nucl. Instrum.
Meth. A 594:407–419.
J. Albert et al. (2008b). ‘Implementation of the Random Forest Method for the Imaging Atmospheric Cherenkov Telescope MAGIC’. Nucl. Instrum. Meth. A 588:424–432.
J. Aleksic et al. (2012). ‘Performance of the MAGIC stereo system obtained with Crab Nebuladata’. Astropart. Phys. 35:435–448.
J. Aleksić et al. (2016). ‘The major upgrade of the MAGIC telescopes, Part I: The hardware
improvements and the commissioning of the system’. Astropart. Phys. 72:61–75.
J. Aleksić et al. (2016). ‘The major upgrade of the MAGIC telescopes, Part II: A performance
study using observations of the Crab Nebula’. Astropart. Phys. 72:76–94.
R. Alfaro et al. (2017). ‘All-particle cosmic ray energy spectrum measured by the HAWC
experiment from 10 to 500 TeV’. Phys. Rev. D 96(12):122001.
M. Amenomori et al. (2008). ‘The All-particle spectrum of primary cosmic rays in the wide
energy range from 10**14 eV to 10**17 eV observed with the Tibet-III air-shower array’.
Astrophys. J. 678:1165–1179.
M. Amenomori et al. (2021a). ‘First Detection of sub-PeV Diffuse Gamma Rays from the
Galactic Disk: Evidence for Ubiquitous Galactic Cosmic Rays beyond PeV Energies’. Phys.
Rev. Lett. 126(14):141101.
M. Amenomori et al. (2021b). ‘Potential PeVatron supernova remnant G106.3+2.7 seen in thehighest-energy gamma rays’. Nature Astron. 5(5):460–464.
D. Anderson, et al. (2016). ‘Precision timing calorimeter for high energy physics’. Nucl. Instrum.
Meth. A 824:670–673.
W. D. Apel et al. (2013). ‘KASCADE-Grande measurements of energy spectra for elementalgroups of cosmic rays’. Astropart. Phys. 47:54–66.
M. Araya (2014). ‘Fermi LAT Observation of Supernova Remnant HB9’. Mon. Not. Roy. Astron.
Soc. 444(1):860–865.

BIBLIOGRAPHY

122

T. Asaba et al. (2018). ‘The overview of the ALPACA Experiment’. PoS ICRC2017:827.
A. M. Atoyan, et al. (1995). ‘Electrons and positrons in the galactic cosmic rays’. Phys. Rev. D
52(6):3265–3275.
W. B. Atwood et al. (2009). ‘The Large Area Telescope on the Fermi Gamma-ray Space Telescope
Mission’. Astrophys. J. 697:1071–1102.
W. Baade & F. Zwicky (1934). ‘On Super-novae’. Proceedings of the National Academy of
Science 20(5):254–259.
A. Bamba, et al. (2003). ‘Fine structures of shock of SN 1006 with the Chandra observation’.
Astrophys. J. 589:827.
Y. Bao & Y. Chen (2021). ‘On the Hard Gamma-Ray Spectrum of the Potential PeVatron
Supernova Remnant G106.3 + 2.7’. Astrophys. J. 919(1):32.
A. R. Bell (1978). ‘The Acceleration of cosmic rays in shock fronts. I’. Mon. Not. Roy. Astron.
Soc. 182:147–156.
B. Bencheikh, et al. (1992). ‘A Simple light detector gain measurement technique’. Nucl.
Instrum. Meth. A 315:349–353.
C. Benn & S. Ellison (1998). ‘Brightness of the night sky over La Palma’. New Astronomy
Reviews 42(6):503–507.
K. Bernlöhr et al. (2013). ‘Monte Carlo design studies for the Cherenkov Telescope Array’.
Astropart. Phys. 43:171–188.
A. Biland, et al. (2007). ‘The Active Mirror Control of the MAGIC Telescope’. In 30th International Cosmic Ray Conference, vol. 3, pp. 1353–1356.
D. J. Bird et al. (1995). ‘Detection of a cosmic ray with measured energy well beyond the
expected spectral cutoff due to cosmic microwave radiation’. Astrophys. J. 441:144–150.
R. Blandford & J. Ostriker (1978). ‘Particle Acceleration by Astrophysical Shocks’. Astrophys.
J. Lett. 221:L29–L32.
P. Blasi & E. Amato (2012). ‘Diffusive propagation of cosmic rays from supernova remnants in
the Galaxy. II: anisotropy’. JCAP 01:011.
A. D. Bolatto, et al. (2013). ‘The CO-to-H2 Conversion Factor’. Ann. Rev. Astron. Astrophys.
51:207–268.
A. Bonardi, et al. (2014). ‘A new solution for mirror coating in γ-ray Cherenkov Astronomy’.
Exper. Astron. 38(1-2):1–9.
L. Breiman (2001). ‘Random Forests’. Machine Learning 45(1):5–32.

BIBLIOGRAPHY

123

T. Bretz, et al. (2009). ‘The drive system of the Major Atmospheric Gamma-ray Imaging
Cherenkov Telescope’. Astropart. Phys. 31:92–101.
Z. Cao, et al. (2021). ‘Ultrahigh-energy photons up to 1.4 petaelectronvolts from 12 γ-ray
Galactic sources’. Nature 594(7861):33–36.
Z. Cao et al. (2021). ‘Peta–electron volt gamma-ray emission from the Crab Nebula’. Science373(6553):425–430.
D. Caprioli, et al. (2009). ‘On the escape of particles from cosmic ray modified shocks’. Mon.
Not. Roy. Astron. Soc. 396:2065–2073.
M. Cardillo, et al. (2016). ‘Supernova remnant W44: a case of cosmic-ray reacceleration’. Astron.
Astrophys. 595:A58.
V. Cavasinni, et al. (2006). ‘TeV Neutrinos from SuperNova Remnants embedded in Giant
Molecular Clouds’. Astropart. Phys. 26:41–49.
D. J. H. Chung, et al. (1998). ‘Superheavy dark matter’. Phys. Rev. D 59:023501.
D. H. Clark, et al. (1975). ‘408 and 5000 MHz Observations of 28 New Galactic Supernova
Remnants’. Australian Journal of Physics Astrophysical Supplement 37:1–38.
J. J. Condon, et al. (1994). ‘A 4.85 GHz Sky Survey. III. Epoch 1986 and Combined (1986-1987)
Maps Covering 0degrees ¡ delta ¡+75degrees’. The Astronomical Journal 107:1829.
P. Cristofari (2021). ‘The Hunt for Pevatrons: The Case of Supernova Remnants’. Universe
7(9):324.
Y. Cui, et al. (2018). ‘Leaked GeV CRs from a Broken Shell: Explaining 9 Years of Fermi-LAT
Data of SNR W28’. Astrophys. J. 860(1):69.
T. M. Dame, et al. (2001). ‘The Milky Way in molecular clouds: A New complete co survey’.
Astrophys. J. 547:792–813.
M. D’Angelo, et al. (2018). ‘Diffuse gamma-ray emission from self-confined cosmic rays around
Galactic sources’. Mon. Not. Roy. Astron. Soc. 474(2):1944–1954.
A. Daum et al. (1997). ‘First results on the performance of the HEGRA IACT array’. Astropart.
Phys. 8:1–11.
M. de Naurois & D. Mazin (2015). ‘Ground-based detectors in very-high-energy gamma-ray
astronomy’. Comptes Rendus Physique 16:610–627.
C. D. Dermer (1986). ‘Secondary production of neutral pi-mesons and the diffuse galactic gamma
radiation’. Astronomy and Astrophysics 157(2):223–229.
J. M. Dickey & F. J. Lockman (1990). ‘HI in the galaxy’. Ann. Rev. Astron. Astrophys. 28:215–
261.

BIBLIOGRAPHY

124

R. Diesing & D. Caprioli (2019). ‘Spectrum of Electrons Accelerated in Supernova Remnants’.
Phys. Rev. Lett. 123(7):071101.
M. Doro et al. (2008). ‘The reflective surface of the MAGIC telescope’. Nucl. Instrum. Meth.
A 595:200–203.
B. Draine (2011). ‘Physics of the Interstellar and Intergalactic Medium’. Physics of the Interstellar and Intergalactic Medium by Bruce T. Draine. Princeton University Press, 2011.
ISBN: 978-0-691-12214-4 .
K. S. Dwarakanath, et al. (1982). ‘Observations of the supernova remnants HB 9 and IC 443 at
34.5 MHz’. Journal of Astrophysics and Astronomy 3:207–216.
P. Eger, et al. (2011). ‘A multi-wavelength study of the unidentified TeV gamma-ray source
HESS J1626-490’. Astronomy and Astrophysics 526:A82.
K. Fang, et al. (2022). ‘Evidence for PeV Proton Acceleration from Fermi-LAT Observations of
SNR G106.3+2.7’. Phys. Rev. Lett. 129(7):071101.
K. Fang & K. Murase (2021). ‘Multimessenger Implications of Sub-PeV Diffuse Galactic GammaRay Emission’. Astrophys. J. 919(2):93.
E. Fermi (1949). ‘On the Origin of the Cosmic Radiation’. Physical Review 75(8):1169–1174.
V. Fomin, et al. (1994). ‘New methods of atmospheric Cherenkov imaging for gamma-ray
astronomy. 1: The False source method’. Astropart. Phys. 2:137–150.
C. Fruck, et al. (2014). ‘A novel LIDAR-based Atmospheric Calibration Method for Improving
the Data Analysis of MAGIC’. In Proceedings, 33rd International Cosmic Ray Conference(ICRC2013): Rio de Janeiro, Brazil, July 2-9, 2013, p. 1054.
Y. Fujita, et al. (2021). ‘X-Ray Emission from the PeVatron-candidate Supernova Remnant
G106.3+2.7’. The Astrophysical Journal 912(2):133.
Y. Fujita, et al. (2009). ‘Molecular Clouds as a Probe of Cosmic-Ray Acceleration in a Supernova
Remnant’. Astrophys. J. Lett. 707:L179–L183.
Y. Fujita, et al. (2011). ‘Alfven Wave Amplification and Self-Containment of Cosmic-Rays
Escaping from a Supernova Remnant’. Mon. Not. Roy. Astron. Soc. 415:3434.
Y. Fukui et al. (2012). ‘A Detailed Study of the Molecular and Atomic Gas Toward the \gammaray SNR RX J1713.7-3946: Spatial TeV \gamma-ray and ISM Gas Correspondence’. Astrophys. J. 746:82.
Y. Fukui, et al. (2021). ‘Pursuing the origin of the gamma rays in RX J1713.7−3946 quantifying
the hadronic and leptonic components’ .
S. Funk (2015). ‘Ground- and Space-Based Gamma-Ray Astronomy’. Ann. Rev. Nucl. Part.
Sci. 65:245–277.

BIBLIOGRAPHY

125

S. Funk & J. A. Hinton (2009). ‘Monte-Carlo studies of the angular resolution of a future
Cherenkov gamma-ray telescope’. AIP Conf. Proc. 1085(1):878–881.
S. Gabici & F. Aharonian (2016). ‘Gamma-ray emission from young supernova remnants:
hadronic or leptonic?’. EPJ Web Conf. 121:04001.
S. Gabici & F. A. Aharonian (2007). ‘Searching for galactic cosmic ray pevatrons with multi-TeV
gamma rays and neutrinos’. Astrophys. J. Lett. 665:L131.
S. Gabici, et al. (2009). ‘Broad-band non-thermal emission from molecular clouds illuminatedby cosmic rays from nearby supernova remnants’. Monthly Notices of the Royal Astronomical
Society 396(3):1629–1639.
D. Gaggero, et al. (2018). ‘Time evolution of gamma rays from supernova remnants’. Mon. Not.
Roy. Astron. Soc. 475(4):5237–5245.
X. Y. Gao, et al. (2011). ‘A Sino-German λ6 cm polarization survey of the Galactic plane. V.
Large supernova remnants’. Astronomy and Astrophysics 529:A159.
M. Garcia-Munoz, et al. (1987). ‘Cosmic-Ray Propagation in the Galaxy and in the Heliosphere:
The Path Length Distribution at Low Energy’. The Astrophysical Journal Supplement Series
64:269.
C. Ge, et al. (2021). ‘Revealing a peculiar supernova remnant G106.3+2.7 as a petaelectronvolt
proton accelerator with X-ray observations’. The Innovation 2:100118.
Y. Génolini et al. (2017). ‘Indications for a high-rigidity break in the cosmic-ray diffusion
coefficient’. Phys. Rev. Lett. 119(24):241101.
Y. Génolini et al. (2019). ‘Cosmic-ray transport from AMS-02 boron to carbon ratio data:
Benchmark models and interpretation’. Phys. Rev. D 99(12):123028.
A. Giuliani et al. (2011). ‘Neutral pion emission from accelerated protons in the supernova
remnant W44’. Astrophys. J. Lett. 742:L30.
D. A. Green (1988). ‘A Revised Reference Catalogue of Galactic Supernova Remnants’. Astrophysics and Space Science 148(1):3–74.
K. Greisen (1956). Progress in Cosmic Ray Physics, Vol. 3. ed. J. G. Wilson. North-Holland,
Amsterdam.
J. Hahn et al. (2014). ‘Impact of aerosols and adverse atmospheric conditions on the data quality
for spectral analysis of the H.E.S.S. telescopes’. Astropart. Phys. 54:25–32.
J. Halpern, et al. (2001a). ‘Psr j2229+6114: discovery of an energetic young pulsar in the error
box of the egret source 3eg j2227+6122’. Astrophys. J. 552:L125.
J. Halpern, et al. (2001b). ‘A possible x-ray and radio counterpart of the high-energy gamma-ray
source 3eg j2227+6122’. Astrophys. J. 547:323.

BIBLIOGRAPHY

126

Y. Hanabata et al. (2014). ‘Detailed Investigation of the Gamma-Ray Emission in the Vicinityof SNR W28 with FERMI-LAT’. Astrophys. J. 786:145.
R. Hartman et al. (1999). ‘The Third EGRET catalog of high-energy gamma-ray sources’.
Astrophys. J. Suppl. 123:79.
T. Hassan et al. (2017). ‘Monte Carlo Performance Studies for the Site Selection of the Cherenkov
Telescope Array’. Astropart. Phys. 93:76–85.
S. Hayakawa (1956). ‘Supernova Origin of Cosmic Rays*)’. Progress of Theoretical Physics
15(2):111–121.
V. F. Hess (1912). ‘Über Beobachtungen der durchdringenden Strahlung bei sieben Freiballonfahrten’. Phys. Z. 13:1084–1091.
M. H. Heyer, et al. (1998). ‘The Five College Radio Astronomy Observatory CO Survey of theOuter Galaxy’. The Astrophysical Journal Supplement Series 115(2):241–258.
A. M. Hillas (1984). ‘The Origin of Ultra-High-Energy Cosmic Rays’. Annual Review of Astronomy and Astrophysics 22:425–444.
A. M. Hillas (1985). ‘Cerenkov Light Images of EAS Produced by Primary Gamma Rays andby Nuclei’. In 19th International Cosmic Ray Conference (ICRC19), Volume 3, vol. 3 ofInternational Cosmic Ray Conference, p. 445.
M. Holler et al. (2016). ‘Observations of the Crab Nebula with H.E.S.S. Phase II’. PoS
ICRC2015:847.
T. Inoue, et al. (2012). ‘Toward Understanding the Cosmic-ray Acceleration at Young Supernova
Remnants Interacting with Interstellar Clouds: Possible Applications to RX J1713.7-3946’.
Astrophys. J. 744:71.
P. S. Iroshnikov (1964). ‘Turbulence of a Conducting Fluid in a Strong Magnetic Field’. Soviet
Astronomy 7:566.
D. Ivanov (2020). ‘Energy Spectrum Measured by the Telescope Array’. PoS ICRC2019:298.
T. Jogler & S. Funk (2016). ‘Revealing W51c as a Cosmic ray Source Using Fermi-lat Data’.
Astrophys. J. 816(2):100.
G. Joncas & L. A. Higgs (1990). ‘The DRAO galactic-plane survey. II. Field at l= 105.’.
Astronomy and Astrophysics, Suppl. Ser. 82:113–144.
E. Kafexhiu, et al. (2014). ‘Parametrization of gamma-ray production cross-sections for pp
interactions in a broad proton energy range from the kinematic threshold to PeV energies’.
Phys. Rev. D 90(12):123014.
S. F. Kamijima & Y. Ohira (2021). ‘Escape of cosmic rays from perpendicular shocks in theinterstellar magnetic field’ .

BIBLIOGRAPHY

127

S. F. Kamijima & Y. Ohira (2022). ‘Escape of cosmic rays from perpendicular shocks in thecircumstellar magnetic field’ .
V. M. Kaspi, et al. (1996). ‘A Search for Radio Pulsars in Southern Supernova Remnants’. The
Astronomical Journal 111:2028.
N. Kawanaka & S.-H. Lee (2021). ‘Origin of Spectral Hardening of Secondary Cosmic-Ray
Nuclei’. Astrophys. J. 917(2):61.
T. W. Kephart & T. J. Weiler (1996). ‘Magnetic monopoles as the highest energy cosmic ray
primaries’. Astropart. Phys. 4:271–279.
A. N. Kolmogorov (1991). ‘The Local Structure of Turbulence in Incompressible Viscous Fluid
for Very Large Reynolds Numbers’. Proceedings of the Royal Society of London Series A
434(1890):9–13.
R. Kothes, et al. (2006). ‘The Boomerang PWN G106.6+2.9 and the Magnetic Field Structurein Pulsar Wind Nebulae’. Astrophys. J. 638(1):225–233.
R. Kothes, et al. (2001). ‘The Supernova Remnant G106.3+2.7 and Its Pulsar-Wind Nebula:
Relics of Triggered Star Formation in a Complex Environment’. The Astrophysical Journal560(1):236–243.
K. Koyama, et al. (1995). ‘Evidence for shock acceleration of high-energy electrons in the
supernova remnant SN1006’. Nature 378:255–258.
R. H. Kraichnan (1967). ‘Inertial Ranges in Two-Dimensional Turbulence’. Physics of Fluids
10(7):1417–1423.
T. L. Landecker et al. (2000). ‘The synthesis telescope at the dominion radio astrophysicalobservatory’. Astron. Astrophys. Suppl. Ser. 145:509.
D. A. Leahy & B. Aschenbach (1995). ‘ROSAT X-ray observations of the supernova remnant
HB 9.’. Astronomy and Astrophysics 293:853–858.
D. A. Leahy & R. S. Roger (1991). ‘Radio Emission from the Supernova Remnant G160.9+2.6(HB9)’. The Astronomical Journal 101:1033.
D. A. Leahy & W. Tian (2007). ‘Radio Spectrum and Distance of the SNR HB9’. Astron.
Astrophys. 461:1013.
T.-P. Li & Y.-Q. Ma (1983). ‘Analysis methods for results in gamma-ray astronomy’. Astrophys.
J. 272:317–324.
T. Linden & B. J. Buckman (2018). ‘Pulsar TeV Halos Explain the Diffuse TeV Excess Observedby Milagro’. Phys. Rev. Lett. 120(12):121101.
P. Lipari & S. Vernetto (2018). ‘Diffuse Galactic gamma ray flux at very high energy’. Phys.
Rev. D 98(4):043003.

BIBLIOGRAPHY

128

S. Liu, et al. (2020). ‘Hadronic vs leptonic models for γ-ray emission from VER J2227+608’.
Astrophys. J. 897(2):L34.
M. A. Malkov, et al. (2013). ‘Analytic Solution for Self-regulated Collective Escape of Cosmic
Rays from their Acceleration Sites’. Astrophys. J. 768:73.
D. Mattingly (2005). ‘Modern tests of Lorentz invariance’. Living Rev. Rel. 8:5.
J. R. Mattox et al. (1996). ‘The Likelihood Analysis of EGRET Data’. Astrophys. J. 461:396.
N. M. McClure-Griffiths, et al. (2005). ‘The Southern Galactic Plane Survey: HI observations
and analysis’. Astrophys. J. Suppl. 158:178–187.
R. Mirzoyan (1997). ‘On the Calibration Accuracy of Light Sensors in Atmospheric Cherenkov
Fluorescence and Neutrino Experiments’. In International Cosmic Ray Conference, vol. 7 ofInternational Cosmic Ray Conference, p. 265.
S. Nagataki (2004). ‘High-energy neutrinos produced by interactions of relativistic protons inshocked pulsar winds’. Astrophys. J. 600:883–904.
D. Nakajima (2013). ‘New Imaging camera for the MAGIC-I Telescope’. In 33rd International
Cosmic Ray Conference, p. 0787.
S. Nozaki et al. (2020). ‘Calibration and performance of the readout system based on switchedcapacitor arrays for the Large-Sized Telescope of the Cherenkov Telescope Array’. Proc. SPIE
Int. Soc. Opt. Eng. 11447:114470H.
A. Obermeier, et al. (2012). ‘The boron-to-carbon abundance ratio and Galactic propagation ofcosmic radiation’. Astrophys. J. 752:69.
Y. Ohira & K. Murase (2019). ‘Origin and Impacts of the First Cosmic Rays’. Phys. Rev. D100(6):061301.
Y. Ohira, et al. (2010). ‘Escape-limited Model of Cosmic-ray Acceleration Revisited’. Astron.
Astrophys. 513:A17.
Y. Ohira, et al. (2011). ‘Gamma-rays from molecular clouds illuminated by cosmic rays escaping
from interacting supernova remnants’. Monthly Notices of the Royal Astronomical Society410(3):1577–1582.
Y. Ohira, et al. (2012). ‘Escape of cosmic-ray electrons from supernova remnants’. Mon. Not.
Roy. Astron. Soc. 427:91.
M. Padovani, et al. (2020). ‘Impact of Low-Energy Cosmic Rays on Star Formation’. Space Sci.
Rev. 216(2):29.
J. S. Perkins & G. Maier (2009). ‘VERITAS Telescope 1 Relocation: Details and Improvements’.

BIBLIOGRAPHY

129

S. Pineault & G. Joncas (2000). ‘G106.3+2.7: A Supernova Remnant in a Late Stage of Evolution’. The Astronomical Journal 120(6):3218–3225.
T. A. Porter, et al. (2008). ‘Inverse Compton Origin of the Hard X-Ray and Soft Gamma-Ray
Emission from the Galactic Ridge’. Astrophys. J. 682:400–407.
V. S. Ptuskin & V. N. Zirakashvili (2003). ‘Limits on diffusive shock acceleration in supernova
remnants in the presence of cosmic-ray streaming instability and wave dissipation’. Astron.
Astrophys. 403:1–10.
V. S. Ptuskin & V. N. Zirakashvili (2005). ‘On the spectrum of high-energy cosmic rays producedby supernova remnants in the presence of strong cosmic-ray streaming instability and wave
dissipation’. Astron. Astrophys. 429:755–765.
W. Reich, et al. (2003). ‘35 cm observations of a sample of large supernova remnants’. Astronomy
and Astrophysics 408:961–969.
S. P. Reynolds, et al. (2012). ‘Magnetic fields in supernova remnants and pulsar-wind nebulae’.
Space Sci. Rev. 166:231–261.
S. Ritt, et al. (2010). ‘Application of the DRS chip for fast waveform digitizing’. Nucl. Instrum.
Meth. A 623:486–488.
R. S. Roger, et al. (1999). ‘The radio emission from the Galaxy at 22 MHz’. Astronomy andAstrophysics Supplement Series 137:7–19.
H. Sano et al. (2019). ‘Possible Evidence for Cosmic-Ray Acceleration in the Type Ia SNR
RCW 86: Spatial Correlation between TeV Gamma rays and Interstellar Atomic Protons’.
Astrophys. J. 876(1):37.
K. M. Schure & A. R. Bell (2013). ‘Cosmic ray acceleration in young supernova remnants’. Mon.
Not. Roy. Astron. Soc. 435:1174.
T. Schweizer, et al. (2002). ‘The optical calibration of the MAGIC telescope camera’. IEEE
Trans. Nucl. Sci. 49:2497–2503.
A. Sezer, et al. (2019). ‘Discovery of recombining plasma inside the extended gamma-ray supernova remnant HB9’. Mon. Not. Roy. Astron. Soc. 489(3):4300–4310.
Y. Shikaze et al. (2007). ‘Measurements of 0.2 to 20-GeV/n cosmic-ray proton and helium
spectra from 1997 through 2002 with the BESS spectrometer’. Astropart. Phys. 28:154–167.
J. Sitarek, et al. (2013). ‘Analysis techniques and performance of the Domino Ring Sampler
version 4 based readout for the MAGIC telescopes’. Nucl. Instrum. Meth. A 723:109–120.
H. Suzuki, et al. (2020). ‘Study on the escape timescale of high-energy particles from supernova remnants through thermal X-ray properties’. Publ. Astron. Soc. Jap. 72(5):Publications of the Astronomical Society of Japan, Volume 72, Issue 5, October 2020, 72,
https://doi.org/10.1093/pasj/psaa061.

BIBLIOGRAPHY

130

H. Suzuki, et al. (2022). ‘Observational Constraints on the Maximum Energies of Accelerated
Particles in Supernova Remnants: Low Maximum Energies and a Large Variety’. Astrophys.
J. 924(2):45.
G. A. Tammann, et al. (1994). ‘The Galactic Supernova Rate’. The Astrophysical Journal
Supplement Series 92:487.
T. Tanaka, et al. (2020). ‘Shock-Cloud Interaction in the Southwestern Rim of RX J1713.7−3946
Evidenced by Chandra X-ray Observations’. Astrophys. J. Lett. 900(1):L5.
T. Tanimori et al. (1994). ‘Observation of 7-TeV gamma-rays from the Crab using the large
zenith angle air Cherenkov imaging technique’. Astrophys. J. Lett. 429:L61–L64.
A. R. Taylor, et al. (2003). ‘The Canadian Galactic Plane Survey’. The Astronomical Journal
125(6):3145–3164.
S. Thoudam, et al. (2016). ‘Cosmic-ray energy spectrum and composition up to the ankle: thecase for a second Galactic component’. Astron. Astrophys. 595:A33.
M. Tluczykont, et al. (2014). ‘The HiSCORE concept for gamma-ray and cosmic-ray astrophysics
beyond 10 TeV’. Astropart. Phys. 56:42–53.
N. Tsuji, et al. (2021). ‘Systematic Study of Acceleration Efficiency in Young Supernova Remnants with Nonthermal X-ray Observations’. Astrophys. J. 907(2):117.
Y. Uchiyama, et al. (2010). ‘Gamma-ray Emission from Crushed Clouds in Supernova Remnants’. Astrophys. J. Lett. 723(1):L122.
Y. Uchiyama, et al. (2012). ‘Fermi-LAT Discovery of GeV Gamma-ray Emission from theVicinity of SNR W44’. Astrophys. J. Lett. 749:L35.
V. P. Utrobin & N. N. Chugai (2017). ‘Luminous Type IIP SN 2013ej with high-velocity 56 Niejecta’. Mon. Not. Roy. Astron. Soc. 472(4):5004–5010.
V. Verzi (2020). ‘Measurement of the energy spectrum of ultra-high energy cosmic rays using
the Pierre Auger Observatory’. PoS ICRC2019:450.
I. Vovk, et al. (2018). ‘Spatial likelihood analysis for MAGIC telescope data - From instrument
response modelling to spectral extraction’. Astron. Astrophys. 619:A7.
T. C. Weekes et al. (1989). ‘Observation of TeV gamma rays from the Crab nebula using theatmospheric Cerenkov imaging technique’. Astrophys. J. 342:379–395.
D. G. Wentzel (1974). ‘Cosmic-ray propagation in the galaxy: collective effects’. Ann. Rev.
Astron. Astrophys. 12:71–96.
J. B. Z. Whiteoak & A. J. Green (1996). ‘The MOST supernova remnant catalogue (MSC).’.
Astronomy and Astrophysics Supplement Series 118:329–380.

BIBLIOGRAPHY

131

Y. Xin, et al. (2019). ‘VER J2227+608: A Hadronic PeVatron Pulsar Wind Nebula?’. The
Astrophysical Journal 885(2):162.
S. Yanagita & K. Nomoto (1999). ‘Nucleosynthetic View of Supernovae Origin of the Galactic
Cosmic Rays’. Astrophysical Letters and Communications 38:461.
S. Yanagita, et al. (1990). ‘Supernovae Origin of Cosmic Rays’. In International Cosmic Ray
Conference, vol. 4 of International Cosmic Ray Conference, p. 44.
H. Yasuda & S.-H. Lee (2019). ‘Time Evolution of Broadband Nonthermal Emission fromSupernova Remnants in Different Circumstellar Environments’. Astrophys. J. 876(1):27.
Y. S. Yoon et al. (2017). ‘Proton and Helium Spectra from the CREAM-III Flight’. Astrophys.
J. 839(1):5.
Q. Yuan, et al. (2017).
95(8):083007.

‘Propagation of cosmic rays in the AMS-02 era’.

Phys. Rev. D

V. Zabalza (2015). ‘naima: a Python package for inference of relativistic particle energy distributions from observed nonthermal spectra’. Proc. of International Cosmic Ray Conference
2015 p. 922.
R. Zanin (2013). ‘MARS, the MAGIC analysis and reconstruction software’. In 33rd International Cosmic Ray Conference, p. 0773.

Every reasonable effort has been made to acknowledge the owners of copyright material. I would
be pleased to hear from any copyright owner who has been omitted or incorrectly acknowledged.