Hadronic interaction models and IACTs
Michiko Ohishi
Institute for Cosmic Ray Research, University of Tokyo

In collaboration with
L. Arbeletche, V. de Souza (Univ. de São Paulo),
G. Maier (DESY), K. Bernlöhr(MPIK), A. Moralejo(IFAE),
J. Bregeon, L. Arrabito (LUPM, Univ. de Montpellier, CNRS/IN2P3),
T. Yoshikoshi (ICRR, Univ. of Tokyo)

Outline Hadronic showers and IACT observations
 Cosmic ray proton simulations in IACT field Difference of hadronic interaction models in shower
particles (focusing onto the effect on g-ray obs. ) p0 spectrum
 energy fraction consumed in electromagnetic (EM) components

 Effect on the observables in IACTs : CTA case
 Energy scale and shower rate of cosmic-ray proton Basic shower parameters and g-hadron separation MVAparameters In the viewpoint of model verification
 Difference in g-ray-like event rate Contribution from heavy nuclei
 gamma-ray sensitivity

 Summary1

Hadronic shower and IACT observations• IACT system detects Cherenkov photons fromsub-EM showers (primarily from p0) and muonscontained in a hadronic shower

Schematic diagram of a hadronic showerprotonp0• Energy spectra and angular distribution ofsecondary particles are different from model tomodel• Related studies in IACT field : Cherenkov photon density (Parsons+ 2011)
 Muon flux on the ground (Mitchell+ 2019) Nature of g-ray-like proton events (sub-EMshowers mimic gamma-ray showers)(Maier+ 2007, Sitarek+ 2017)

nmnm

ｈｈｈ2

Hadronic shower and IACT observations• IACT system detects Cherenkov photons fromsub-EM showers (primarily from p0) and muonscontained in a hadronic shower

Schematic diagram of a hadronic shower

protonp0• Energy spectra and angular distribution of
secondary particles are different from model tomodel• Related studies in IACT field : Cherenkov photon density (Parsons+ 2011)
 Muon flux on the ground (Mitchell+ 2019) Nature of g-ray-like proton events (sub-EMshowers mimic gamma-ray showers)(Maier+ 2007, Sitarek+ 2017)m-ringnmn

mshower*Simulated camera images from Mitchell+ (2019)ｈｈｈ3

Hadronic shower and IACT observations• IACT system detects Cherenkov photons from
sub-EM showers (primarily from p0) and muonscontained in a hadronic shower

Schematic diagram of a hadronic shower

protonp0

Sub-EM shower• Energy spectra and angular distribution of
secondary particles are different from model tomodel• Related studies in IACT field : Cherenkov photon density (Parsons+ 2011)
 Muon flux on the ground (Mitchell+ 2019) Nature of g-ray-like proton events (sub-EMshowers mimic gamma-ray showers)(Maier+ 2007, Sitarek+ 2017)
 Discrimination capability of model difference
depends on the array performance - this study is
focused on CTA, testing QGSJET-II-03 (currentlyused in CTA) and recent post-LHC models

In caseEp0 close to Ep

mng-like….
ｈｈｈ4

Current IACT systems and CTA (array scale)• Current IACT arrays (H.E.S.S., VERITAS, MAGIC): coverage of 0.03 km2• CTA :  4 km2 for South site (99 telescopes)
 0.6 km2 for North site (19 telescopes)→ Full containment of Cherenkov photons from g-ray and proton showers

Array configuration (South site), public at
https://www.cta-observatory.org/science/cta-performance/

Array scale and Cherenkov light-pool size5

Cosmic ray proton simulations in IACT field Current IACT systems

 Real cosmic-ray data (“OFF-source” data) are used as
background samples

Berge+ (2007)

 Real OFF-source data are also used in training of machine learning
for g-hadron separation and estimation of residual background in
the derivation of the sensitivity curve→ Usually protons are not simulated, as far as we are analyzing
gamma-ray sources …..
6

Cosmic ray proton simulations in IACT field Current IACT systems
 Real cosmic-ray data (“OFF-source” data) are used as backgroundsamples Real OFF-source data are also used in training of machine learning for ghadron separation and estimation of residual background CTA (and systems in design/construction phase) Monte Carlo (MC) simulation data are used for background estimation
 Usually cosmic-ray protons and electrons are simulated as backgrounds As for proton: currently interaction between cosmic-ray proton andnuclei in very-high-energy region is not perfectly understood several hadronic interaction models (QSGJET, EPOS, SIBYLL…) are in
use in VHE/UHE CR field Improvement of models with feedback from collider and CR
experiments is ongoing7

Proton simulation and g-ray sensitivity• g-ray sensitivity of an IACT system is mostly determined by• Significance of signal events to the background fluctuation (≥5s)• Signal-to-background ratio (≥5%)
Differential Sensitivity of CTA South arrayz = 20deg, 50h obs.

Ns≥5*1
Ng/NBG≥0.05

g-ray effective area

Ng≥10
Signal event
statistics

Significance to
BG fluctuation

Background rate

S/B ratio“Background”≈ CR proton+ electron
CTA Instrument Response Functions (IRFs), public at
https://www.cta-observatory.org/science/cta-performance/*1 Significance def. in
Li & Ma (1983), Eq. (17)8Difference of models in shower particles - p0 spectrum 10% ofprimaryenergyp0 spectra, Ep=1 TeV (mono)

primary •protonenergyAir shower simulation with CORSIKA toinvestigate difference of secondary
particles between different models• Used models:- QGSJET-II-03 in CORSIKA6.99(currently used in CTA)
- QGSJET-II-04, EPOS-LHC,SIBYLL2.3c in CORSIKA7.69- E<80 GeV: fixed low energy model
UrQMD (for all cases)• p0 spectrum- Spectrum at high energy end canaffect the rate of g-ray-like events- Harder spectrum tends to give
more g-ray-like BG events:
EPOS → SIBYLL→ QGSJET-II-03  QGSJET-II-04

9Difference of models in shower particles- Energy fraction in EM components -80%

• Energy fraction carried by g +e-+e+ (EM components) after the 3rd interaction(as for g-ray primary case, this fraction is close to 100%)• Similar pattern as p0 spectrum is seen; relation between model changes at 1 TeV• Energy fraction in EM which will be regarded as “g-ray-like” event depends on thesystem performance -- 80% was used in this study for CTA
Prob. of high EM fraction events VS true E
EEM/Eprimary distribution, Ep=1 TeV casecorrelate with g-ray-like event rate

10

Effect on the observables in IACTs : CTA simulation
Array configuration, South Site
N

Analysis tool: EventDisplay v500-rc04

Site

Paranal (Chile)

Array

4 LSTs, 25 MSTs, 70 SSTs
(configuration shown left)

Particle

Gamma, e-,proton: QGSJET-II-03 *1
proton :QGSJET-II-04
EPOS-LHC v3.4 /SIBYLL2.3c*2
Low Energy Model (E<80 GeV)
: fixed as UrQMD

Core range

2500 mViewcone

0 - 10 deg

Energy range

0.003 - 330 TeV (e-, gamma）
0.004 - 600 TeV (proton)Spec. index-2.0 *3*1 in CORSIKA 6.99, produced on GRID system in EU*2 in CORSIKA 7.69, produced on cluster in Japan*3 Reweighted in the analysis

11Energy scale and shower rate
• Difference in p0 production can lead to difference in E scale and CR proton rate
• 5% difference in reconstructed energy and 10% difference in CR proton rate
between models (before gamma-ray selection cuts)
Erec VS Etrue
CR proton shower rate (relative)(E is reconstructedassuming g-ray )

5%

10%

12

Difference in basic shower parameter distribution“Lateral size ofthe shower”

precut line

gamma

1 TeV<Erec<10 TeV“Longitudinal sizeof the shower”

Proton histograms
are normalized bynumber ofsimulated events• Most important shower characteristics for ghadron separation : WIDTH (lateral size of theshower)• MSCW : corrected and normalized WIDTH• Difference between models is seen at small
MSCW (g-ray-like region)“Height of showermaximum”

13

MVA parameters for g-hadron separation• Multivariate analysis (MVA) to introduce a single index of “gammaness” (or hadroness)- Boosted Decision Tree is used here, with precuts in basic shower parameters• EPOS and SIBYLL show worse separation, with more g-like events than QGS as expected
MVA parameter distribution
1.0 TeV <Erec< 5.6 TeV
BDT
trained for each model𝑸𝑪𝜸 / 𝑪𝒑 VS BDT cut value
QGSGoodseparation

Badseparation

Cg: cutacceptance for g

SIBYLL
EPOS

Histograms are normalized by
the area (num. of acceptedevents)

14

In the viewpoint of model verification with IACTs
MVA parameter distribution (1 TeV < Erec < 10 TeV)

p-likeContribution from e- is consideredg-like• Once we have real CR data,
we can test which model isthe closest to the reality bycomparing MC and real data :- Event rate
- Shower param. dist.
- g-hadron separation
parameter dist.
(relatively large factor 2difference)• Current IACT systems canalso contribute to model+100%
verification, though model
discrimination capability
depends on the system
performance (worse thanCTA).

identical trained BDT (QGSJETII-03) is used for all models

15Model verification: contribution from heavy nuclei?
MVA parameter distribution (1 TeV < Erec < 10 TeV) • Uncertainty in CR composition canaffect the model verification accuracy(EPOS)

g-likep-like• As far as treating g-ray-like events,
contribution from heavy nuclei isnegligibly small→ good verification measure• Helium and heavier nuclei do not
mimic g-rays because of their lateralsize and shower maximum heightHelium

He flux is assumed to be same level asproton (as an extremity)

gamma

He/p ratio

Height of shower max.

gamma

Lateral size of showerpHe

Histograms
normalized by theareapHe

16

Effect on the gamma-ray sensitivity
South site, LST+MST+SST array, z=20deg, average of North+South pointing

Signal event statistics

Contribution from e-

Low-energy interaction model
50h case+30%

Contribution from e-

Background rate(p+e-)

Differential sensitivity

+100%

17

Effect on the gamma-ray sensitivity…..
South site, LST+MST+SST array, z=20deg, average of North+South pointing

Signal event statistics

Contribution from e-Low-energy interaction model
50h case+30%

Contribution from e-

Background rate(p+e-)

Differential sensitivity+100%18

Summary Effect of difference in hadronic interaction models on gamma-ray sensitivity ofCTA south array (99tels, 4-LSTs + 25-MSTs + 70-SSTs) was estimated with MC
simulation data Tried models:- QGSJET-II-03 in CORSIKA6.99 (currently used for CTA IRF)
- QGSJET-II-04, EPOS-LHC, SIBYLL2.3c in CORSIKA7.69 (post-LHC models) 5% level difference in energy scale and 10% level difference in proton shower
rate were seen.
 As a preliminary result, difference in g-ray sensitivity between models wasestimated to be 30% level (with 10% statistical error from MC data); Relation
between models is consistent with p0 spectrum and EM fraction
 In the viewpoint of model verification, g-ray-like event rate is a relatively goodmeasure : almost free from uncertainty of cosmic-ray nuclei composition
 relatively large (factor 2) difference between models
 Current IACT systems can also contribute to model testing (discrimination
capability depends on the system performance )19

Backup slides

20Differential sensitivity – subsystems LST only

50h case

MST only

50h case

SST only

50h case21

CR spectra used in the background calculation
 CR proton𝑑𝑁
𝐸 −Γ= 𝐼0𝑑𝐸𝐸𝐶

2.62

𝐼0 = 9.8 × 10−6 cm−2 s −1 TeV −1 str −1 , 𝐸𝐶 = 1.0 TeV, Γ = CR electron(log 10 (𝐸/𝐸𝐶 )−𝜇)2𝑑𝑁
𝐸 −Γ
3𝐸 𝑑𝐸 = 𝐼0 𝐸× (1 + 𝑓 × (exp(exp(−))-1))
2𝜎2𝐶

𝐼0 = 2.385 × 10−9 cm−2 s −1 TeV −1 str −1 , 𝐸𝐶 = 1.0 TeV,
Γ = 0.43,
𝜇 = −0.101 , 𝜎 = 0.741, 𝑓 = 1.950

22

High EM fraction event prob. VS trueE
EM frac>70% prob.

EM frac>60% prob.

23

Differential sensitivity
South site, LST+MST+SST array, z=20deg, average of North+South pointing
Background rate(p+e-)

Differential sensitivity

50h case+30%+100%

24

CR charge measurement using Direct Cherenkov lightDCRelation on the
focal plane⚫ Detect Cherenkov photons which primaryparticle emits before the first interaction

宇宙線粒子

Number of Ch. photons per length:

Inelasticscattering𝑑𝑁𝑐
sin2 𝜃𝑐2= 2𝜋𝑍 𝛼 න𝑑𝜆𝑑𝑠𝜆2𝑒2α : 4𝜋𝜀 ℏ𝑐0dNc/ds  Z2

Direct Cherenkov → Charge (Z)Shower→ Energy ・Arrival direction・
core position ･Mass number (A)Schematic view of Direct
Cherenkov methods using
telescopes on the ground1cos𝜃𝑐 =𝑛𝛽⚫
→ We can estimate the primary charge
number using number of Cherenokov photons
detected on the ground⚫ Proposed by Kieda et al. (2001)⚫ H.E.S.S. measured cosmic ray iron spectrum
25
from 13 TeV to 200 TeV using this method(2007)
4