PYTHIA 8 for CORSIKA 8

(1)

PYTHIA 8 for CORSIKA 8

Torbj¨ orn Sj¨ ostrand

Department of Astronomy and Theoretical Physics, Lund University

CORSIKA 8 Air-Shower Simulation and Development Workshop, Heidelberg, 12 – 15 July 2022

(2)

A new framework for hadronic collisions

Based on 2 articles by Marius Utheim & TS:

“A Framework for Hadronic Rescattering in pp Collisions”, Eur. Phys. J. C80 (2020) 907, arXiv:2005.05658

“Hadron Interactions for Arbitrary Energies and Species, with Applications to Cosmic Rays”,

Eur. Phys. J. C82 (2022) 21, arXiv:2108.03481

Models arbitrary hadron–hadron collisions at low energies.

Models arbitrary hadron-p/n collisions at any energy.

Initialization slow, ∼ 15 minutes,

? but thereafter works for any hadron–p/n at any energy, and

? initialization data can be saved, so only need to do once.

The Angantyr nuclear geometry part used to extend to hadron-nucleus at any energy.

Native C++ simplifies interfacing Pythia 8↔ Corsika 8.

So far limited comparisons with data.

Torbj¨orn Sj¨ostrand PYTHIA 8 for CORSIKA 8 slide 2/22

(3)

PYTHIA and the structure of an LHC pp collision

MPI MPI dˆσ0

·

· ·

··

Meson Baryon Antibaryon

·Heavy Flavour

Hard Interaction Resonance Decays MECs, Matching & Merging FSR

ISR*

QED Weak Showers Hard Onium

Multiparton Interactions Beam Remnants*

Strings

Ministrings / Clusters Colour Reconnections String Interactions Bose-Einstein & Fermi-Dirac Primary Hadrons

Secondary Hadrons Hadronic Reinteractions (*: incoming lines are crossed)

ISR*

Strings

ISR*

Strings

From Pythia 8.3 guide, arXiv:2203.11601, 315 pp

(4)

Space–time evolution

Pythia can now calculate production vertex of each particle, e.g. number of hadrons as a function of time for pp at 13 TeV:

time(fm/c) 1 10 10² 10³ 10⁴ 10⁵ 10⁶ 10⁷ 10⁸ 10⁹ 10¹⁰ 10¹¹ 10¹² 10¹³ 10¹⁴ 10¹⁵

hadn

0 20 40 60 80 100 120 140 160

180 Total number of hadrons

Primary hadrons Secondary hadrons Total number of final hadrons

S. Ferreres-Sol´e, TS, EPJC 78, 983

(5)

Hadronic rescattering

13 TeV nondiffractive pp events:

0 2 4 6 8 10

τ (fm) 0

5 10 15 20 25 30 35 40

dN/dτ

Invariant production time, rescattered or not rescattered not rescattered sum both rescattering off

Pythia now contains framework for hadronic rescattering:

1)Space–time motion and scattering opportunities 2)Cross section for low-energy hadron–hadron collisions

3)Final-state topology in such collisions

Already covered by other programs like UrQMD or SMASH, but then interfacing issues limits usefulness.

Rescattering not important enough to switch on for Corsika, but code extends Pythia modelling down to Ekin≈ 0.2 GeV.

(6)

Total cross sections (1)

Two examples of total and partial cross sections:

10^-1 10⁰ 10¹ 10² 10³ 10⁴ ECM (GeV)

0 20 40 60 80 100

σ (mb)

π⁺p cross sections as a function of energy total

LEHE NDLE HEel LEHE D/ELE HE

10^-1 10⁰ 10¹ 10² 10³ 10⁴ ECM (GeV)

0 20 40 60 80 100

σ (mb)

K⁺p cross sections as a function of energy total

LEHE NDLE HEel LEHE D/ELE HE

Combines low-energy (LE) and high-energy (HE) formalisms.

LE: process-specific with detailed structure, cf. UrQMD.

HE: smooth behaviour with pomeron + reggeon ansatz.

(7)

Total cross sections (2)

11/21 Motivation and background

The rescattering framework Results

Cross sections

Based on work by Pelaez, Rodas, Ruiz de Elvira et al.

(arXiv:1102.2183, arXiv:1907.13162, arXiv:1602.08404)

Marius Utheim Hadronic Rescattering in Pythia/Angantyr

When nothing else available, use Additive Quark Model, where nq,AQM= nd+ nu+ 0.6ns+ 0.2nc+ 0.07nb

σ^AB = n^A_q,AQM 3

n_q,AQM^B 3 σ^pp

(8)

Low-energy interaction types

Low-energy cross sections, roughly up to Ecm= 10 GeV, split as σ_tot= σel+σND+σ_SD(XB)+σ_{SD(AX )}+σDD+σCD+σexc+σann+σres+. . .

9/21 Motivation and background

The rescattering framework Results

Low-energy interactions

Elastic

Diffractive

Non-diffractive

Annihilation

Resonant

Marius Utheim Hadronic Rescattering in Pythia/Angantyr

Excitation≈ low-mass diffraction.

Also includes modelling of particle production in these types, where string fragmentation barely is applicable.

(9)

Multiparton interactions

Hadrons are composite

⇒ can contain many parton–parton interactions.

In part perturbatively calculable, but divergence for p⊥→ 0, impact-parameter description, colour reconnection, and more.

dσint

dp²_⊥ = Z Z

dx1dx2f1(x1, p_⊥²) f2(x2, p²_⊥) dˆσ dp²_⊥ dˆσ

dp²_⊥ ∝ α²_s(p_⊥²)

p_⊥⁴ → α²_s(p_⊥0² + p²_⊥) (p_⊥0² + p_⊥²)²

MPI machinery gradually switched on for Ecm > 10 GeV.

Related to multipomeron models, but matched to high-p⊥ jets.

(10)

Parton distribution functions

When Ecm > 10 GeV and we allow QCD 2→ 2 processes, we need to implement PDF sets for “all” hadrons, not only p:

0.0 0.2 0.4 0.6 0.8 1.0

x

0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

xv(x,Q

2 0)

Valence contents for some mesons at Q20= 0.26 qv

sKv qKv cDv qDv bBv qBv

0.0 0.2 0.4 0.6 0.8 1.0

x

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

xv(x,Q

2 0)

Valence contents for some baryons at Q₂₀= 0.26 dpv

upv sv⁺ uv⁺ cv^{+ +}^c uv^{+ +}^c bv^b⁺ uv^b⁺

Some generators omit hard processes, to simplify and speed up.

Meaningful at low energies, but more dubious at high energies, and e.g. for charm production.

(11)

Flexible beams and energies

Low energy: easy to switch between colliding hadrons and energies.

High energy: hard, since perturbative description (with PDFs) need initialization at desired energy to be used.

10⁰ 10¹ 10² 10³ 10⁴

ECM (GeV) 0

20 40 60 80 100

ncharged®

Rise of charged multiplicity with energy low-energy model

high-energy model interpolation

Smooth transition from low to high energies starting at 10 GeV.

High energies by interpolation in fixed grid of initialized energies for MPI.

Also fast switch between beam combinations

by array of such grids, which takes time, but can be reused.

(12)

Why not ANGANTYR?

Angantyr is a rather new module of Pythia,

intended for pp, pA and AA collisions at LHC and RHIC.

Provides a good description of overall event properties, even if some aspects are not (yet) fully described.

Seemingly ideal for applications such as collisions on atmospheric nuclei, but several limitations:

Only addresses interactions of p and n.

Only intended for large energies, say √

s > 100 GeV.

Initialization of nuclear geometry (and MPIs) very slow,

∼ 1 minute, and then only works for one fixed energy.

On the positive side, subsequent event generation reasonably fast, in particular relative to full hydrodynamics.

Also, involves no intra-nuclear cascades in hadron-nucleus collsion, only collisions related to incoming hadron.

(13)

ANGANTYR geometry (1)

A hadron passing through a nucleus can undergo a variable number nsub of hadron-nucleon subcollisions.

Use Angantyr nuclear geometry package, combined with total cross sections above, to find n_sub distributions:

0 2 4 6 8 10 12 14

nsub

10^-3 10^-2 10^-1 10⁰

Probability

Number of subcollisions for p N 10 GeV 100 GeV 1000 GeV 10000 GeV 100000 GeV

0 5 10 15 20 25 30 35 40

nsub

10^-3 10^-2 10^-1 10⁰

Probability

Number of subcollisions for p Pb 10 GeV 100 GeV 1000 GeV 10000 GeV 100000 GeV

Note an approximate geometrical series = straight line in log-plot

⇒ characterized by one number, hnsubi.

(14)

ANGANTYR geometry (2)

hnsubi depends on target A, hadron projectile h, and collision energy, but the latter two mainly in the combination σ_total^hp (s):

0 20 40 60 80 100 120 140

σtotal (mb) 1.0

1.5 2.0 2.5 3.0

nsub®

Average number of subcollisions for N target crude linear fit, low σ

crude linear fit, high σ all hadrons and energies

0 20 40 60 80 100 120 140

σtotal (mb) 1

2 3 4 5 6 7 8 9

nsub®

Average number of subcollisions for Pb target crude linear fit, low σ

crude linear fit, high σ all hadrons and energies

So need one parametrization ofhnsubi(σ^hp_total) for each A.

Constraint: expecthnsubi → 1 when σtotal → 0.

Currently¹⁴N,¹⁶O,⁴⁰Ar, ²⁰⁸Pb.

(Takes∼ 1 hour manual work to set up and fit for each A.)

(15)

Nuclear cross sections

In limit σ_total^hp → 0 expect σ_total^hA = A σ_total^hp (cf. ν beam).

For a “normal” σ_total^hp , nuclear geometry suggests σ^hA_total≈A^2/3σ_total^hp andhnsubi ≈A^1/3.

These two limits are consistent with figures on previous slide.

Assuming preserved total number of hp/hn subcollisions, whether free or bound nucleons, gives

σ_total^hA hnsubi = A σ^hp_total ⇒ σ_total^hA = A σ^hp_total

hnsubi .

For applications in code:

σ^hp(s) (total and partial) is reasonably time-consuming, so only do it once for a hadron of given energy,

hnsubi is quick and simple to evaluate,

so can afterwards easily get σ_total^hA (s) for different A.

(16)

ANGANTYR multiplicities

If several “independent” hp/hn subcollisions in an hA one, then naively dn^hA_charged

dy =hnsubidn_charged^hp dy

but net shift in A direction owing to momentum conservation.

10 5 0 5 10

y 0

2 4 6 8 10 12

dnch/dy

Charged rapidity distribution Angantyr Pythia scaled

10 5 0 5 10

y 0

2 4 6 8 10

dnch/dy

Angantyr vs. Pythia p N emulation Angantyr Pythia

Patched up by letting to hardest hadron in each subcollision go on to be the projectile in the next subcollision (+ some more).

(17)

Collision generation flow

Assume an incoming hadron h collides with a nucleus A.

1 np= Z , nn= A− Z .

2 Do big loop A times, unless a break before that.

3 After first time in the loop:

Keep going with probability 1− 1/hn^subi, else break.

Pick up the new product with largest longitudinal momentum, and assign it to represent the h.

Break if not enough hp/hn invariant mass for collision.

Set changed composition nondiffractive/other processes.

4 Pick target nucleon p or n in proportions np: nn. Subtract 1 from the chosen one.

5 Do the hp/hn collision, as an isolated event.

6 Copy new particles into common event and update history.

7 End of big loop.

8 Optionally do decays and/or compress the event record.

(18)

Atmospheric cascade evolution

incoming 10⁸GeV proton

0 250 500 750 1000 1250 1500 1750

X (g/cm²) 100

10¹ 102

10³ 10⁴

Probability

Atmospheric depth of interactions

Uniform p/n atmosphere Uniform nitrogen Exponential nitrogen Exponential nitrogen at 45

0 250 500 750 1000 1250 1500 1750

X (g cm²) 10²

10³ 10⁴ 10⁵ 10⁶

(1/nev)X 0dN

Number of muons at depth

0 250 500 750 1000 1250 1500 1750

X (g cm²) 10²

10³ 10⁴ 10⁵

(1/nev)X 0dN

Number of hadrons at depth

0 2 4 6 8 10 12 14

Ekin (GeV) 10⁴

10⁵ 10⁶ 10⁷

(1/nev)dn/dEkin

Kinetic energy spectrum of leptons at surface Uniform p/n atmosphere Uniform nitrogen Exponential nitrogen Exponential nitrogen at 45

(19)

The PythiaCR wrapper

In examples/main183.cc a hadronic cascade is traced through the atmosphere, but poor substitute for full Corsika tracking.

The tentative examples/main184.cc attempts to offer a wrapper where Corsika can do the tracking, but calls PythiaCR to

provide the hA collision cross section, perform an hA collision, or

perform an h decay.

Internally two Pythia instances

PythiaMain administrates an hA collision, and does an h decay,

PythiaColl does an hp/hn subcollision, and provides the hp/hn cross section.

(20)

PythiaCR methods

The public PythiaCR methods/references (currently) are PythiaCRconstructor,

init initializes all program elements, sigmaSetup calculates a hp cross section, sigmaColl calculates a hA cross section, based on the hp one above,

nextCollperforms an hA collision, nextDecay performs an h decay,

compressreduces the event record to final particles only, stat prints error statistics at the end of the run,

particleData(), rndm() references that can be used in the main program for particle data or random numbers.

(21)

The forward region

1000 2000 3000 4000 5000 6000

]-1 dN/dE [GeVine1/N

−10

10

−9

10

−8

10

−7

10

−6

10

−5

10

=13TeV photon s

LHCf

° φ=180

∆ > 10.94, η

Ldt=0.191nb-1

∫

Data QGSJET II-04 EPOS-LHC DPMJET 3.06 SIBYLL 2.3 PYTHIA 8.212

Energy [GeV]

1000 2000 3000 4000 5000 6000

MC/Data

0 1 2 3

4 1000 2000 3000 4000 5000 6000

]-1 dN/dE [GeVine1/N

−10

10

−9

10

−8

10

−7

10

−6

10

−5

10

=13TeV photon s

LHCf

° φ=20

∆

<8.99, η 8.81<

Ldt=0.191nb-1

∫

Energy [GeV]

1000 2000 3000 4000 5000 6000

MC/Data

0 1 2 3 4

Figure 4: Comparison of the photon spectra obtained from the experimental data and MC predictions. The top panels show the energy spectra, and the bottom panels show the ratio of MC predictions to the data. The hatched areas indicate the total uncertainties of experimental data including the statistical and the systematic uncertainties.

Acknowledgments

We thank the CERN staff and the ATLAS Collaboration for their essen- tial contributions to the successful operation of LHCf. This work was partly

290

supported by JSPS KAKENHI Grant Numbers JP26247037, JP23340076 and the joint research program of the Institute for Cosmic Ray Research (ICRR), University of Tokyo. This work was also supported by Istituto Nazionale di Fisica Nucleare (INFN) in Italy. Parts of this work were performed using the computer resource provided by ICRR (University of Tokyo), CERN and CNAF

295

(INFN).

References

[1] A. Aab et al. (Pierre Auger Collaboration), Nucl. Instrum. Methods Phys.

Res., Sect. A 798 (2015) 172.

14

Forward region important for cosmic-ray physics.

Current description poor, but two recent improvements:

Forbid popcorn mechanism for remnant diquarks only;

i.e. baryon then always produced at end of string, never meson

Allow harder fragmentation function for leading baryon, by separate a and b values.

Tuning for FPF by Max Fieg shows clear improvement.

(22)

Summary and outlook

Pythia now may offer a realistic alternative to current hadronic-interaction models.

Time-efficient access to perturbative activity for collisions of “any” incoming hadron at “any” energy, by rapid beam switching and energy interpolation of MPI parameters.

Not in any other generator (?), but does it matter?

Includes new low- and high-energy cross section descriptions, and new sets of PDFs for a wide range of hadrons.

Interface should allow easy calling from Corsika for both interactions and decays.

Limitations: no incoming nuclei, no nuclear target breakup, no hadronic γ interactions, . . . ,

. . . so work to be done to meet more requirements.

Feedback welcome. How interesting?