# The Greisen Function and its Ability to Describe Air-Shower Profiles

Maximilian Stadelmaier  Institute of Physics of the Czech Academy of Sciences,
Prague Vladimír Novotný  Institute of Physics of the Czech Academy of
Sciences, Prague Institute of Particle and Nuclear Physics, Faculty of
Mathematics and Physics, Charles University, Prague Jakub Vícha  Institute of
Physics of the Czech Academy of Sciences, Prague

(March 2023)

###### Abstract

Ultrahigh-energy cosmic rays are almost exclusively detected through extensive
air showers, which they initiate upon interaction with the atmosphere. The
longitudinal development of these air showers can be directly observed using
fluorescence detector telescopes, such as those employed at the Pierre Auger
Observatory or the Telescope Array. In this article, we discuss the properties
of the Greisen function, which was initially derived as an approximate
solution to the electromagnetic cascade equations, and its ability to describe
the longitudinal shower profiles. We demonstrate that the Greisen function can
be used to describe longitudinal air-shower profiles, even for hadronic air
showers. Furthermore we discuss the possibility to discriminate between
hadrons and photons from the shape of air-shower profiles using the Greisen
function.

##  I Introduction

Extensive air showers are created by cosmic rays upon interaction with the
atmosphere [1, 2]. They can be detected at the ground using surface detector
arrays, or directly observed at night using fluorescence detector telescopes.
To reconstruct the shower development and shower observables, a profile
function needs to be fitted to the detector data. Gaisser and Hillas proposed
an empiric function [3] to describe the longitudinal development of proton air
showers as an alternative to the Constant Intensity Cut method [4, 5], which
is used in surface detector experiments to take into account the atmospheric
attenuation of particles in air showers from different zenith angles. It was
shown in [6] that the Gaisser–Hillas (GH) function can be used to approximate
a system of particles being created and absorbed in an extended
Heitler–Matthews model [7], and can be adjusted to very closely match the
Greisen function, which is an approximation for the solutions to the
electromagnetic cascade equations [8]. Both the Pierre Auger Observatory and
the Telescope Array use the GH function to describe their fluorescence
detector data [9, 10].

In this article, we will discuss the Greisen function and its properties, such
as its connection to the shower age. We will demonstrate the usability of the
Greisen function as an alternative to the GH function to fit longitudinal
shower profiles and present its performance to reconstruct the depth of the
shower maximum as well as the primary energy using Monte Carlo (MC)
simulations of air showers.

##  II The Greisen Function

The average longitudinal development of electromagnetic air showers can be
very well described analytically [11]. This description holds in good
approximation also for hadronic showers, initiated by ionized nuclei, which
make up the upper end of the cosmic-ray energy spectrum [12]. The solutions to
the cascade equations derived by Rossi and Greisen under Approximation
A111Approximation A is the high-energy approximation to the cascade equations,
in which only bremsstrahlung and pair production are considered as relevant
processes. to describe extensive air showers were used to motivate important
properties in the context of air-shower physics, such as the shower age

$\displaystyle s=\frac{3t}{t+2\ln(E_{0}/E_{\text{cut}})}$ (1)

that describes the development of an electromagnetic shower, initiated by a
primary particle of the energy $E_{0}$, after $t$ radiation lengths,
considering only the particles above an energy of $E_{\text{cut}}$. Note that
$s=1$ at $t=\ln(E_{0}/E_{\text{cut}})$; this value is usually assigned with
the shower maximum. For electromagnetic showers a reasonable choice for
$E_{\text{cut}}$ is close to ${\approx}87\,\text{MeV}$, which is the energy
above which electromagnetic particles on average lose more energy in radiative
shower processes than to scattering and ionization. Furthermore, it was
demonstrated that the relative rate of change222We use the notation
$\lambda_{1}$, even though we do not mean to suggest that this quantity is to
be understood as a (wave) length, to adhere to historic convention.
$\lambda_{1}$ in the number $N$ of particles as a function of the surpassed
radiation lengths $t$,

$\displaystyle\lambda_{1}=\frac{1}{N(t)}\frac{\partial N(t)}{\partial t},$ (2)

is similar for all electromagnetic showers at high primary energies. Greisen
introduced the approximation [8]

$\displaystyle\lambda_{1}\simeq\frac{1}{2}\left(s-1-3\ln s\right),$ (3)

which is in good agreement with the exact solution. Originally, the parameter
$s$ describes the spectra $n$ of electromagnetic particles in a shower, that
is approximately given by $n_{\upgamma}\sim n_{\text{e}^{\pm}}\sim E^{-(s+1)}$
for particles at energies $E\ll E_{0}$, but from Eq. 3 there exists a relation
between $\lambda_{1}$ and $s$.

It is straightforward to combine Eqs. 1, 2 and 3 and to solve the resulting
expression for $N(t)$ by integration. This yields

$\displaystyle N(t)=N_{0}\exp{\left[t\left(1-\tfrac{3}{2}\ln
s\right)\right]},$ (4)

with a constant $N_{0}$. The maximum number $N_{\text{max}}$ of particles
above the energy of $E_{\text{cut}}=98\,\text{MeV}$ in a cascade initiated by
a particle of energy $E_{0}$ was derived in [13] under Approximation
B333Approximation B of the cascade equations augments Approximation A by a
term concerning Coulomb scattering. and found to be

$\displaystyle
N_{\text{max}}=\frac{0.31}{\sqrt{\ln(E_{0}/E_{\text{cut}})}}\frac{E_{0}}{E_{\text{cut}}}.$
(5)

Eq. 4 has its maximum at $t_{\text{max}}=\ln(E_{0}/E_{\text{cut}})$, which
evaluates to $N_{0}\,E_{0}/E_{\text{cut}}$. Thus, solving for $N_{0}$ using
Eq. 5 yields the Greisen function, which reads as

$\displaystyle
N(t)=\frac{0.31}{\sqrt{\beta}}\exp{\left[t\left(1-\tfrac{3}{2}\ln
s\right)\right]},$ (6)

using the short notation $\beta=\ln(E_{0}/E_{\text{cut}})=t_{\text{max}}$. The
Greisen function, introduced in [14], is thus an approximate solution to the
electromagnetic cascade equations, given in [11], combining aspects of both
Approximation A and Approximation B. There was no strict derivation given by
Kenneth Greisen himself, but _a-posteriori_ derivations (such as the one
presented here) were provided in [15] and [16].

In the following, we have to overcome two major shortcomings of the Greisen
function as written in Eq. 6. Firstly, the Greisen function in its classical
form cannot accurately describe the point of the first interaction of a cosmic
ray with the atmosphere, since by construction the cascade is initiated always
at $t=0$. Secondly, the scale of the Greisen function is only accurate for
average electromagnetic showers. We will introduce a parameter $\epsilon$ to
account for this issue and demonstrate that the Greisen function generalized
this way is able to describe the longitudinal profile of hadronic showers and
the corresponding shower-to-shower fluctuations.

##  III The Modified Greisen Function

The classical Greisen function, which is given in Eq. 6, assumes that a shower
starts at $t=0$. Furthermore, the Greisen function is technically only able to
describe electromagnetic showers. In this section, we introduce minor
modifications to the function to describe the longitudinal development of both
hadronic and electromagnetic air showers.

Firstly, we introduce a non-zero point of the first interaction at a slanted
atmospheric depth $X_{1}$, that will be described by $t_{1}=X_{1}/X_{0}$ using
the electromagnetic radiation length444Equivalently, we use $t=X/X_{0}$ and
$t_{\text{max}}=X_{\text{max}}/X_{0}$. $X_{0}\simeq
37\,\text{g}\,\text{cm}^{-2}$. Thus, the shower age $s$ will be given as

$\displaystyle s=\frac{3t^{\prime}}{t^{\prime}+2\beta}\,\Theta(t^{\prime}),$
(7)

with $t^{\prime}=t-t_{1}$ and the Heaviside function $\Theta$. To maintain the
property of the shower age, which is supposed to be 1 at the maximum of the
shower, and to keep number of radiation lengths required to reach the maximum
of the shower the same as before, we redefine $\beta$ in accordance with the
previous modification as

$\displaystyle\beta=\ln(E_{0}/E_{\text{cut}})=t_{\text{max}}-t_{1}.$ (8)

Finally, we introduce the factor $\epsilon$, which is defined in units of
energy deposit per step length, and which can be interpreted as the effective
energy loss per particle and step length at the shower maximum555Here we
ignore the factor of 0.31 from Eq. 5.. Thus, the modified Greisen profile
reads as

$\displaystyle
N(t)\equiv\frac{\mathrm{d}E}{\mathrm{d}X}(t)=\frac{\epsilon}{\sqrt{\beta}}\exp\left[{(t-t_{1})\left(1-\tfrac{3}{2}\ln
s\right)}\right],$ (9)

with $N(t)=0$ for $t\leq t_{1}$. For the sake of simplicity, here and in the
following, in the text we abbreviate the energy deposit
$\mathrm{d}E/\mathrm{d}X$ with the symbol $N$, analogously to the number of
particles.

##  IV Calibration of the Greisen Profile

If Eq. 9 is used to describe the longitudinal profiles of (hadronic) showers
in terms of deposited energy rather than a number of particles, it is
necessary to examine viable (effective) numerical values of the energy
$E_{\text{cut}}$ as well as of $\epsilon$.

We investigate the behaviour of $\epsilon$ and $E_{\text{cut}}$ using the MC
values of $X_{\text{max}}$, $X_{1}$, $E_{0}$, and the maximum energy deposit
$N(t_{\text{max}})$ of the longitudinal profiles of simulated air showers. The
simulations were produced using the Sibyll 2.3d [17], Epos-LHC [18], and the
Qgsjet II-04 [19] models of hadronic interactions with different primary
particles at different primary energies. All simulations were produced in the
Conex event generator [20] at version v7.60. We produced 1000 simulated
showers at primary energies of $10^{18.5}\,\text{eV}$, $10^{19}\,\text{eV}$,
and $10^{19.5}\,\text{eV}$ with gamma-ray, proton, and iron nuclei primary
particles, each.

Shower-to-shower fluctuations of (hadronic) air showers will severely affect
the maximum number of particles produced in a shower as well as the absolute
depth of the shower maximum. These fluctuations, however, can be accurately
reproduced by the behaviour of the Greisen function. Rewriting Eq. 8 in terms
of the slanted atmospheric depth of the shower maximum $X_{\text{max}}$ and
the average radiation length $X_{0}$,

$\displaystyle\beta=\frac{X_{\text{max}}-X_{1}}{X_{0}},$ (10)

we find that

$\displaystyle
E_{\text{cut}}=E_{0}\,\mathrm{e}^{-\beta}=E_{0}\exp{\left[-\frac{X_{\text{max}}-X_{1}}{X_{0}}\right]}.$
(11)

Note that the explicit dependence of $E_{\text{cut}}$ on $E_{0}$ is cancelled
by the dependence of the average $X_{\text{max}}$ on the logarithm of the
primary energy in Eq. 8. We choose an effective value for the radiation length
of $X_{0}=40\,\text{g}\,\text{cm}^{-2}$ as a compromise for the different
considered primary particles.

In a similar manner, using the numerical maximum of a given shower profile
according to Eq. 9, the parameter $\epsilon$ can be identified as

$\displaystyle\epsilon=N(t_{\text{max}})\,\sqrt{\beta}\,\mathrm{e}^{-\beta}.$
(12)

Figure 1: The behaviour of $\epsilon$ and $E_{\text{cut}}$ calculated from the
Monte-Carlo values of $X_{1}$, $X_{\text{max}}$, $N(t_{\text{max}})$, and
$E_{0}$ from simulated air showers produced with the Sibyll2.3d model of
hadronic interactions. $E_{\text{cut}}$ and $\epsilon$ were obtained using Eq.
11 and Eq. 12, respectively. The markers show data points from individual
shower simulations, the line shows a log-linear regression.

As can be seen in Fig. 1, we observe a strong correlation for the MC values of
$\epsilon$ and $E_{\text{cut}}$ (and thus $\beta$), assuming Greisen-like
longitudinal profiles. Furthermore, we show in Fig. 8 and Fig. 9 that the
distributions of $\epsilon$ and $E_{\text{cut}}$ as well as their relation are
approximately the same even for all three hadronic interaction models. This
behaviour indicates that shower-to-shower fluctuations of hadronic showers are
not at all random, but follow certain regularities. For example, a deeper than
average value of $X_{\text{max}}$ corresponds to smaller than average value of
$N(t_{\text{max}})$ (and vice-versa), if all other parameters are fixed.
Furthermore, we expect larger values of $\beta$ for photon-induced showers
than for protons or iron nuclei. Between different primaries there is a
gradual transition in the shape of the shower from very hadronic (iron-like)
to proton-like and lastly electromagnetic showers, where the behaviour of
$\epsilon$ and $E_{\text{cut}}$ can be described by a power law,

$\displaystyle\frac{\epsilon}{\text{PeV}/\text{g}\,\text{cm}^{-2}}\simeq\left(\frac{E_{\text{cut}}}{10^{16.8}\,\text{eV}}\right)^{0.97},$
(13)

with residuals on average within 0.2%. Note that if expressed in terms of
$X_{\text{max}}$ and $X_{1}$ (cf. Eq. 11), $E_{\text{cut}}$ is not an explicit
parameter of the Greisen function. Eq. 13 thus expresses the universal
relation between the maximum energy deposit and the extent of the shower in
terms of radiation lengths after the first interaction.

Even though $\epsilon$ appears as a pre-factor in the modified Greisen
function, the numerical values of $\epsilon$ from a best fit are independent
of the integrated profile (cf. Eqs. 13 and 11 and Fig. 9) and thus of the
energy of the primary particle. Because $E_{\text{cut}}$ does not depend on
the primary energy, $\epsilon$ solely depends on the shape of the shower. And
thus, as can be seen in Fig. 1, depends on the amount of hadronization
occuring during the shower development, which relates to the primary mass. The
integrated profile (and thus the calorimetric energy deposit) is mainly
governed by the value of $\beta$ (cf. Eq. 17). For smaller values of
$\epsilon$ (e.g. photon-like showers), the shower takes longer to reach its
maximum in terms of radiation lengths. In case of showers with a significant
amount of hadronization, where multiple cascades are effectively in
superposition (i.e. iron-like showers), we expect larger numerical values for
$\epsilon$, corresponding to an earlier shower maximum. From superposition and
the Heitler–Matthews model666The model stimates the difference of the averages
in $X_{\text{max}}$ for proton and iron showers to be
$\langle{X}_{\text{max}}\rangle_{\text{p}}-\langle{X}_{\text{max}}\rangle_{\text{Fe}}\simeq
150\,\text{g}\,\text{cm}^{-2}$. [7] one would expect
$\epsilon_{\text{Fe}}/\epsilon_{\text{p}}\simeq 50$ (cf. Eqs. 9 and 12);
however, the average ratio obtained from simulations is much smaller. From the
ratio $\epsilon_{\text{Fe}}/\epsilon_{\text{p}}\simeq 10$, we estimate a
difference in $\beta$ of about 2.4 radiation lengths between the average
proton and iron shower to reach the shower maximum. This is in accordance with
simulations, which imply777The effect of the depth of the first interaction is
neglected.
$\langle{X}_{\text{max}}\rangle_{\text{p}}-\langle{X}_{\text{max}}\rangle_{\text{Fe}}\simeq
100\,\text{g}\,\text{cm}^{-2}$.

##  V Fitting Simulated Data

Figure 2: Example shower profiles with corresponding best fits. The showers
were initiated by a (_top to bottom_) gamma-ray, proton, and iron primary
particle with a primary energy of $10^{19}\,\text{eV}$ each, using the
Sibyll2.3d model of hadronic interactions. The best fit values corresponding
to both the Greisen function (_red_) and the Gaisser-Hillas function (_orange
dashed_) are given in each panel along with the MC values of $X_{\text{max}}$.
Additionally, the best fit values of $X_{\text{max}}$ are depicted by vertical
lines. The profile tails (_gray data points_) were disregarded for the fits.
Below each panel, the relative deviation $\delta=(f-d)/d$ of the function
value $f$ and the simulated profile data $d$ is shown for both functions in
the respective color. Figure 3: $\chi^{2}$-distributions for the Greisen and
the GH functions fitted to simulated Sibyll2.3d air showers with primary
energies of $10^{19}\,\text{eV}$. The distributions for the Greisen (GH)
functions are shown as a full (dashed) line, for each of the three primary
particles. The respective mean of each distribution is indicated by a vertical
line.

To test the ability of the Greisen function to describe longitudinal profiles
of air showers, we examine and fit simulated shower profiles. The simulation
library contains the same configuration of primary energies, particles, and
hadronic interaction models as mentioned before. We compare the results
against results from fitting the same showers to the very commonly used GH
function, which in terms of the slanted depth $X$ reads as

$\displaystyle\begin{split}N(X)=N_{\text{max}}&\left(\frac{X-X_{1}}{X_{\text{max}}-X_{1}}\right)^{\tfrac{X_{\text{max}}-X_{1}}{\Lambda}}\\\
&\hskip
50.00008pt\times\exp{\left[\frac{X_{\text{max}}-X}{\Lambda}\right]},\end{split}$
(14)

with $N(X)=0$ for $X\leq X_{1}$, the depth of the shower maximum
$X_{\text{max}}$, the maximum value $N_{\text{max}}$ of the function, and a
characteristic length $\Lambda$, which is related but not equal to the
electromagnetic interaction length. Approximately one expects $\Lambda\simeq
3\,X_{0}/2$ (cf. Eq. 16 and [15, 16]).

To describe the shape of the shower profiles, we use $N_{\text{max}}$,
$X_{\text{max}}$, $X_{1}$, and $\Lambda$ as free parameters for the GH
function, and $\epsilon$, $X_{\text{max}}$, $X_{1}$, and $X_{0}$ for the
Greisen function. We estimate the uncertainty of the individual MC data points
to $0.3\,\text{PeV}/(\text{g}\,\text{cm}^{-2})$. The tails of the profiles
($\mathrm{d}E/\mathrm{d}X\leq 0.8\,\text{PeV}/(\text{g}\,\text{cm}^{-2})$) are
not used in any of the fits. The constant threshold for the tails as well as
the uncertainty for the data points were estimated so that
$\chi^{2}/\text{ndf}$ is $\sim 1$ for a mixture of proton and iron showers
with primary energies of $10^{19}\,\text{eV}$. Examples of fitted profiles of
the photon, proton, and iron nucleus induced showers are given in Fig. 2.

The $\chi^{2}$-distributions for the Greisen and GH functions fitted to
simulated data from Sibyll2.3d showers at a primary energy of
$10^{19}\,\text{eV}$ are depicted in Fig. 3. Additional
$\chi^{2}$-distributions from showers simulated with different hadronic
interaction models and primary energies are given in Fig. 10 (note that the
estimated uncertainty of the individual data points, as well as the limit for
fitting the tails was not adjusted for the different primary energies). We
find that on average for all energies and hadronic interaction models, the
average values of $\chi^{2}/\text{ndf}$ are smaller for the Greisen function
fit, than for the fit using a GH function, thus implying a better match of the
function to the profile data. For iron showers the difference of
$\chi^{2}/\text{ndf}$ is the largest, with values being approximately $20\%$
smaller for the Greisen function. The difference for proton showers is
approximately $10\%$, while for photon showers the average
$\chi^{2}/\text{ndf}$ differs by less than $1\%$.

Furthermore, we investigate the ability of the Greisen function to recover the
depth of the shower maximum as well as the calorimetric energy deposit of the
shower. The distributions of the residuals
$X_{\text{max}}^{\text{rec}}-X_{\text{max}}^{\text{MC}}$ as a function of the
Conex MC values of $X_{\text{max}}$ are depicted in Fig. 4. We observe that
$X_{\text{max}}$ can be recovered very accurately from the simulated profile
data using both functions, with an average precision of about
$4\,\text{g}\,\text{cm}^{-2}$ for both fit functions. The performance of the
Greisen function of “finding” the right depth of the shower maximum is thus
approximately equal as of the GH function.

The calorimetric energy deposit of the shower can be obtained by integrating
the fitted profile function from $X_{1}$ up to $\infty$. To obtain the
calorimetric energy deposit from the best-fit of both functions, we integrate
numerically888The result for the calorimetric energy deposit as obtained from
the fitted and numerically integrated Greisen function is approximately the
same as using the formula given in Eq. 17. from the respective best-fit value
of $X_{1}$ to $2000\,\text{g}\,\text{cm}^{-2}\lesssim\infty$. The relative
residuals of the recovered calorimetric energy deposit $E_{\text{cal}}$ with
respect to the simulated calorimetric energy deposit
$E^{\text{MC}}_{\text{cal}}$ are depicted in Fig. 5 as a function of the Conex
MC values of $X_{\text{max}}$. The accuracy and precision of the recovered
calorimetric energy to estimate the primary energy is the same for the Greisen
and the GH function for all primary particles (and all hadronic interaction
models). The performance of the Greisen function to estimate the primary
energy of the particle initiating the shower is thus equal to the performance
of the GH function.

Additionally, we present two-dimensional distributions of $X_{1}$ and $X_{0}$
obtained from the Greisen function fit, as well as $X_{1}$ and $\Lambda$ from
the GH function fit in Fig. 11. When fitting simulated data, both the Greisen
and the GH function appear to have the same “pathology” to produce mostly
negative best-fit values for $X_{1}$. Even though the best-fit values for
$X_{1}$ obtained from the Greisen function are less negative than from the GH,
this shows that the values obtained for $X_{1}$ from fitted data cannot be
treated as the point of the first interaction. Besides the absolute scale, the
distributions for the best-fit values of $X_{1}$ and $X_{0}$ ($\Lambda$)
depicted in Fig. 11 appear to be very similar for both the Greisen and the GH
function.

Figure 4: Distributions of the residuals of the reconstructed values of
$X_{\text{max}}$ as a function of the Conex MC values of $X_{\text{max}}$
using the Greisen function (_left_) and the GH function (_right_). The showers
were simulated with primary energies of $10^{18.5}\,\text{eV}$,
$10^{19}\,\text{eV}$, and $10^{19.5}\,\text{eV}$, using the Sibyll2.3d model
of hadronic interactions. The overall mean and standard deviation of the
distribution is given in the upper right corner. The distributions of
residuals for the individual primary particles are colored accordingly.

Figure 5: Distributions of the difference of the recovered calorimetric energy
$E_{\text{cal}}$ and the simulated calorimetric energy deposit
$E^{\text{MC}}_{\text{cal}}$ as a function of the Conex MC values of
$X_{\text{max}}$ using the Greisen function (_left_) and the GH function
(_right_). The showers were simulated with primary energies of
$10^{18.5}\,\text{eV}$, $10^{19}\,\text{eV}$, and $10^{19.5}\,\text{eV}$,
using the Sibyll2.3d model of hadronic interactions. The overall mean and
standard deviation of the distribution is given in the upper right corner. The
distributions of residuals for the individual primary particles are colored
accordingly.

##  VI Mass-composition sensitivity of the Greisen function

In terms of the performance to obtain $X_{\text{max}}$ from simulated data,
the Greisen function is not second to the GH function (cf. Fig. 4). Studying
the depth of the shower maximum, the same sensitivity to the mass composition
as expected from the GH function can thus be achieved using the Greisen
function.

Given the fact that $\epsilon$ is independent of the primary energy of the
shower, and that the MC distributions of $\epsilon$ (cf. Fig. 9 (_left_)) are
dependent on the type of primary particle, it is tempting to examine the
primary-mass sensitivity of the best-fit values of $\epsilon$. In Fig. 6 we
show the two-dimensional distributions of the best-fit results of
$X_{\text{max}}$ and $\epsilon$. To remove the direct dependence on the
primary energy, instead of the true $X_{\text{max}}$, here we use

$\displaystyle
X_{\text{max}}^{19}:=X_{\text{max}}-58\,\text{g}\,\text{cm}^{-2}\,\lg\left(E_{\text{cal}}/10^{19}\,\text{eV}\right),$
(15)

assuming a constant decadal elongation rate of approximately
$58\,\text{g}\,\text{cm}^{-2}$, and using $E_{\text{cal}}$ as obtained from
the numerical integration of the best-fit function.

From Fig. 6 it is obvious that the separation of the distributions of
individual primary particles increases when the best-fit values of $\epsilon$
are considered alongside with $X_{\text{max}}$. Numerically, the means of the
distributions of $X_{\text{max}}^{19}$ obtained from a fit to proton and iron
shower profiles are approximately 1.5 average standard deviations apart; on
the diagonal line, which combines the information of $\epsilon$ and
$X_{\text{max}}$, the means of the proton and iron distributions are separated
by almost 1.9 average standard deviations999For the estimation, see Eq. 21..
In case of photon-hadron separation, the distance improves from $1.4$ to
$1.6$. Thus, using Conex simulations, we see a clear improvement in terms of
the separation of primary particles when employing the combination of
$\epsilon$ and $X_{\text{max}}$ from the Greisen function over
$X_{\text{max}}$ only. Additional indicators for photon-like showers are the
obtained values for $\chi^{2}$ and $X_{0}$ (cf. Figs. 3 and 11). The Greisen
function thus might be useful when trying to identify ultrahigh-energy photons
in fluorescence detector data.

Figure 6: Two-dimensional distributions of the best-fit values of
$X_{\text{max}}^{19}$ and $\epsilon$ using the Greisen function to fit
simulated longitudinal profiles of showers with primary energies of
$10^{18.5}\,\text{eV}$, $10^{19}\,\text{eV}$, and $10^{19.5}\,\text{eV}$. The
curved lines show the estimated $1\sigma$ extent of the respective
distributions. All showers were simulated using the Sibyll2.3d model of
hadronic interactions.

To compare against the behaviour of the GH function, in Fig. 12 we show the
two-dimensional distributions of $X_{\text{max}}^{19}$ and
$N_{\text{max}}^{19}=N_{\text{max}}/(E_{\text{cal}}/10^{19}\,\text{eV})$, as
obtained from the GH function fitted to simulated air-shower data. As can be
seen from Fig. 12, $N_{\text{max}}^{19}$ does not yield additional information
about the mass of the primary particle on its own, while $\epsilon$ does.

##  VII Fitting the Greisen function with fixed shape

To boost the performance of the fitting procedure given only poor data, the GH
function can be used with constraints to fix the shape of the function to an
expected shape of the profile data [21]. These constraints are realized by
reparametrizing the GH function in terms of “$L$” and “$R$” and then
constraining these new parameters, which define the width and skewness of the
function, respectively (see [21] for details). However, the average values of
$L$ and $R$ depend on the primary particle of the shower [22, 23]. Using the
Greisen function, the shape of the profile can be determined simply by fixing
or constraining the parameter $\epsilon$. Moreover, one can easily choose
whether the function should resemble an average gamma-ray, proton, or iron
shower, depending on the corresponding value of $\epsilon$ (cf. Fig. 9).

To demonstrate, we fix
$\epsilon=10^{-6.2}\,\text{PeV}/(\text{g}\,\text{cm}^{-2})$ as a compromise
between iron-like and proton-like shower profiles and fit the simulated data
with only three free parameters, namely $X_{1}$, $X_{\text{max}}$, and
$X_{0}$. As can be seen from Fig. 7, the $X_{\text{max}}$ bias is minimal for
proton showers, using the Greisen function with an hadron-like fixed shape,
but increases for gamma-ray-induced showers and heavy nuclei. The estimated
calorimetric energy deposit, however, is almost unaffected (the bias for
hadronic showers changes by $\approx 0.5\%$), when using only three free
parameters. Lastly, the best-fit values of $X_{\text{max}}$ and $X_{1}$ become
highly correlated for the Greisen function with fixed $\epsilon$, as it is
expected from Eqs. 10 and 12.

Figure 7: The $X_{\text{max}}$-bias (_top_) and $E_{\text{cal}}$-bias
(_middle_) as a function of the Conex MC value of $X_{\text{max}}$, as well as
the correlation of $X_{\text{max}}$ and $X_{1}$ using a fixed-shape Greisen
function with only three free parameters to fit simulated air-shower data. The
showers were simulated using the Sibyll2.3d model of hadronic interactions
with primary energies of $10^{18.5}\,\text{eV}$, $10^{19}\,\text{eV}$, and
$10^{19.5}\,\text{eV}$.

##  VIII Discussion and Summary

In this article, we present the Greisen function in its original form and
discuss its relation to the shower age parameter $s$. Furthermore, we present
a way to derive the function from literature. We show that with slight
modifications the Greisen function can be rewritten to match individual
simulated air-shower profiles, even from hadronic primaries. We confirm this
statement using simulated air showers from different primary particles at
different energies, and using different hadronic interaction models. Contrary
to popular belief, the Greisen function matches the simulated air-shower
profiles even somewhat better than the most commonly used Gaisser–Hillas
profile function. In contrary to the Gaisser–Hillas function, which was
introduced as an alternative to the Constant Intensity Cut method, the Greisen
function was derived to describe the full longitudinal profiles of air
showers.

We analyse the performance of the Greisen function to recover air-shower
observables from simulated profile data assuming an ideal detector. In this
analysis, we show that the Greisen function yields approximately the same
performance as the Gaisser–Hillas function to determine the calorimetric
energy deposit and the depths of the shower maxima from simulated showers at
different primary energies using different hadronic interaction models.

We identified the shape-parameter $\epsilon$ of the Greisen function, which is
primary-mass sensitive and can help distinguishing different types of primary
particles, additionally to the slanted depth of the shower maximum
$X_{\text{max}}$. Lastly, we demonstrate that fixing $\epsilon$ can elegantly
fix the shape of the Greisen function and thus shower profiles can be fitted
even with only three free parameters. In this case, while the recovered values
of $X_{\text{max}}$ are slightly biased, the accuracy and precision of the
obtained calorimetric energy deposit are unaffected.

We conclude that the Greisen function proves itself useful to describe air-
shower profiles and bears additional potential for photon-hadron separation as
well as mass-composition studies, and could thus be used in the search for
light particles in air-shower data.

##  IX Acknowledgements

The authors would like to thank Alexey Yushkov, Eva Santos, Armando di Matteo,
and Darko Veberič for fruitful comments and discussion as well as the referee
for valuable comments and critique. This work was partially supported by the
Ministry of Education Youth and Sports of the Czech Republic and by the
European Union under the grant FZU researchers, technical and administrative
staff mobility, registration number CZ.02.2.69/0.0/0.0/18_053/0016627. This
work was supported by the Czech Science Foundation – Grant No. 21-02226M.

##  X Code Availability

The analysis code for this article is available upon request under:

    gitlab.com/stadelmaier/greisen_fct_fit

## References

  * Rossi [1933] B. Rossi, Über die eigenschaften der durchdringenden korpuskularstrahlung im meeresniveau, Z. Phys. 82, 151 (1933).
  * Auger  _et al._ [1939] P. Auger, P. Ehrenfest, R. Maze, J. Daudin, and A. F. Robley, Extensive cosmic ray showers, Rev. Mod. Phys. 11, 288 (1939).
  * Gaisser and Hillas [1977] T. K. Gaisser and A. M. Hillas, Reliability of the Method of Constant Intensity Cuts for Reconstructing the Average Development of Vertical Showers, Proc. Int. Cosmic Ray Conf. 8, 353 (1977).
  * Hersil  _et al._ [1961] J. Hersil, I. Escobar, D. Scott, G. Clark, and S. Olbert, Observations of extensive air showers near the maximum of their longitudinal development, Phys. Rev. Lett. 6, 22 (1961).
  * Alvarez-Muniz  _et al._ [2002] J. Alvarez-Muniz, R. Engel, T. K. Gaisser, J. A. Ortiz, and T. Stanev, Atmospheric shower fluctuations and the constant intensity cut method, Phys. Rev. D 66, 123004 (2002), arXiv:astro-ph/0209117 .
  * Montanus [2012] J. M. C. Montanus, Intermediate models for longitudinal profiles of cosmic showers, Astropart. Phys. 35, 651 (2012), arXiv:1106.1073 [astro-ph.HE] .
  * Matthews [2005] J. Matthews, A Heitler model of extensive air showers, Astropart. Phys. 22, 387 (2005).
  * Greisen [1960] K. Greisen, Cosmic ray showers, Ann. Rev. Nucl. Sci. 10, 63 (1960).
  * Aab  _et al._ [2014] A. Aab _et al._ (Pierre Auger), Depth of Maximum of Air-Shower Profiles at the Pierre Auger Observatory: Measurements at Energies above $10^{17.8}$ eV, Phys. Rev. D 90, 122005 (2014), arXiv:1409.4809 [astro-ph.HE] .
  * Abbasi  _et al._ [2018] R. U. Abbasi _et al._ (Telescope Array), Depth of Ultra High Energy Cosmic Ray Induced Air Shower Maxima Measured by the Telescope Array Black Rock and Long Ridge FADC Fluorescence Detectors and Surface Array in Hybrid Mode, Astrophys. J. 858, 76 (2018), arXiv:1801.09784 [astro-ph.HE] .
  * Rossi and Greisen [1941] B. Rossi and K. Greisen, Cosmic-ray theory, Rev. Mod. Phys. 13, 240 (1941).
  * Kampert and Unger [2012] K.-H. Kampert and M. Unger, Measurements of the Cosmic Ray Composition with Air Shower Experiments, Astropart. Phys. 35, 660 (2012), arXiv:1201.0018 [astro-ph.HE] .
  * Tamm and Belenky [1946] I. Tamm and S. Belenky, The energy spectrum of cascade electrons, Phys. Rev. 70, 660 (1946).
  * Greisen [1956] K. Greisen, The extensive air showers, Progress in Cosmic Ray Physics 3, 3 (1956).
  * Lipari [2009] P. Lipari, The Concepts of ’Age’ and ’Universality’ in Cosmic Ray Showers, Phys. Rev. D 79, 063001 (2009), arXiv:0809.0190 [astro-ph] .
  * Stadelmaier [2022] M. Stadelmaier,  _On air-Shower universality and the mass composition of ultra-high-energy cosmic rays_ , Ph.D. thesis, KIT, Karlsruhe, KIT, Karlsruhe, Dept. Phys. (2022).
  * Riehn  _et al._ [2018] F. Riehn, H. P. Dembinski, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev, The hadronic interaction model SIBYLL 2.3c and Feynman scaling, PoS ICRC2017, 301 (2018), arXiv:1709.07227 [hep-ph] .
  * Pierog  _et al._ [2015] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider, Phys. Rev. C 92, 034906 (2015), arXiv:1306.0121 [hep-ph] .
  * Ostapchenko [2011] S. Ostapchenko, Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model, Phys. Rev. D 83, 014018 (2011), arXiv:1010.1869 [hep-ph] .
  * Bergmann  _et al._ [2007] T. Bergmann, R. Engel, D. Heck, N. N. Kalmykov, S. Ostapchenko, T. Pierog, T. Thouw, and K. Werner, One-dimensional Hybrid Approach to Extensive Air Shower Simulation, Astropart. Phys. 26, 420 (2007), arXiv:astro-ph/0606564 .
  * Aab  _et al._ [2019] A. Aab _et al._ (Pierre Auger), The Pierre Auger Observatory: Contributions to the 36th International Cosmic Ray Conference (ICRC 2019): Madison, Wisconsin, USA, July 24- August 1, 2019 (2019), arXiv:1909.09073 [astro-ph.HE] .
  * Andringa  _et al._ [2011] S. Andringa, R. Conceicao, and M. Pimenta, Mass composition and cross-section from the shape of cosmic ray shower longitudinal profiles, Astropart. Phys. 34, 360 (2011).
  * Andringa  _et al._ [2012] S. Andringa, R. Conceicao, F. Diogo, and M. Pimenta, Sensitivity to primary composition and hadronic models from average shape of high energy cosmic ray shower profiles (2012), arXiv:1209.6011 [hep-ph] .

##  Appendix A Additional Expressions

### Alternative Form of the Greisen Function

The Greisen function can be rewritten as

$\displaystyle
N(t)=\frac{\epsilon}{\sqrt{\beta}}\,\mathrm{e}^{t^{\prime}\left(1-\tfrac{3}{2}\ln
s\right)}=\frac{\epsilon}{\sqrt{\beta}}\left(\frac{3t^{\prime}}{t^{\prime}+2\beta}\right)^{-\tfrac{3}{2}t^{\prime}}\,\mathrm{e}^{t^{\prime}},$
(16)

with $t^{\prime}=t-t_{1}$ and $\beta=t_{\text{max}}-t_{1}$.

### Integral of the Greisen Function

Numerically, we find that the integrated Greisen profile, which yields the
calorimetric energy deposit and is thus a good estimator for the primary
energy, can be approximated within ${\sim}0.5\%$ by

$\displaystyle
E_{\text{cal}}=\int\limits_{X_{1}}^{\infty}\frac{\epsilon}{\sqrt{\beta}}\,\exp\left[(t-t_{1})\left(1-\tfrac{3}{2}\ln
s\right)\right]\mathrm{d}X\simeq 3.1\,\epsilon\,\mathrm{e}^{\beta}\,X_{0},$
(17)

where identified the pre-factor of ${\sim}3.1$ numerically as suitable for
different values of $X_{0}$ and $X_{1}$. Thus, we find that

$\displaystyle\epsilon\simeq\frac{E_{\text{cut}}}{3.1\,X_{0}},$ (18)

which is in good agreement with Eq. 13.

### Proton-Iron Separation Using the Greisen Function

The mean and standard devations of the distributions of $X_{\text{max}}^{19}$
as well as of $\lg(\epsilon/(\text{PeV}/(\text{g}\,\text{cm}^{-2})))$ depicted
in Fig. 6 are given in Table 1. In units of the mean standard deviation of the
distributions, the average values of $X_{\text{max}}^{19}$ of protons and iron
nuclei induced showers are separated by

$\displaystyle\xi=\frac{|\mu_{\text{p}}(X_{\text{max}}^{19})-\mu_{\text{Fe}}(X_{\text{max}}^{19})|}{\sqrt{\sigma_{\text{p}}(X_{\text{max}}^{19})^{2}+\sigma_{\text{Fe}}(X_{\text{max}}^{19})^{2}}}\simeq
1.52.$ (19)

For $\hat{\epsilon}=\lg(\epsilon/(\text{PeV}/(\text{g}\,\text{cm}^{-2})))$,
this distance evaluates to

$\displaystyle\zeta=\frac{|\mu_{\text{p}}(\hat{\epsilon})-\mu_{\text{Fe}}(\hat{\epsilon})|}{\sqrt{\sigma_{\text{p}}(\hat{\epsilon})^{2}+\sigma_{\text{Fe}}(\hat{\epsilon})^{2}}}\simeq
1.18.$ (20)

Assuming no covariance between $\hat{\epsilon}$ and $X_{\text{max}}^{19}$, we
find

$\displaystyle\sqrt{\xi^{2}+\zeta^{2}}=1.93.$ (21)

Note that Eq. 21 is a generous estimate, since it is clear from Fig. 6 that
the covariance between $\hat{\epsilon}$ and $X_{\text{max}}^{19}$ is not zero.

Table 1: Mean $\mu$ and standard deviations $\sigma$ of the distributions of $X_{\text{max}}^{19}$ and $\epsilon$ as depicted in Fig. 6. ($\mu$, $\sigma$) | $X_{\text{max}}^{19}/(\text{g}\,\text{cm}^{-2})$ | $\lg(\epsilon/(\text{PeV}/(\text{g}\,\text{cm}^{-2})))$
---|---|---
$\upgamma$ | (984.4, 107.1) | (-11.312, 4.33)
p | (812.0, 62.2) | (-7.217, 2.57)
Fe | (711.0, 22.5) | (-4.124, 0.49)

##  Appendix B Additional Figures

We show the effect of shower-to-shower fluctuations for showers simulated with
different hadronic interaction models and at different primary energies. Fig.
8 depicts the universal relation between $\epsilon$ and $E_{\text{cut}}$ even
for the Epos-LHC and QgsjetII-04 hadronic interaction models, using different
primary particles and different primary energies. In Fig. 9 we present the
individual distributions of $\epsilon$ and $E_{\text{cut}}$ for all three
hadronic interaction models. Note that the histograms of the distributions for
photon showers in Fig. 9 are slightly truncated, as for individual showers
$\epsilon$ ($E_{\text{cut}}$) reaches values down to
${\approx}1\,\text{keV}/\text{g}\,\text{cm}^{-2}$ (${\approx}1\,\text{MeV}$)
(cf. Figs. 1 and 8).

Figure 8: Relation of $\epsilon$ and $E_{\text{cut}}$ for simulated iron
(blue), proton (red), and gamma-ray (orange) showers using the Epos-lhc
(_left_) and QgsjetII-04 (_right_) model of hadronic interactions.
$E_{\text{cut}}$ and $\epsilon$ were obtained using Eq. 11 and Eq. 12,
respectively. In each panel we show results from 9000 showers using gamma-
rays, protons, and iron nuclei primary particles with simulated energies of
$10^{18.5}\,\text{eV}$, $10^{19.0}\,\text{eV}$, and $10^{19.5}\,\text{eV}$
(1000 for each configuration).

Figure 9: Distributions of $\epsilon$ (_left_) and $E_{\text{cut}}$ (_right_)
calculated according to Eqs. 11 and 12 for showers simulated using the (_top
to bottom_) Sibyll2.3d, Epos-LHC, and the QgsjetII-04 models of hadronic
interactions. In each panel we show results from 9000 showers using gammay-
rays, protons, and iron nuclei as primary particles with simulated energies of
$10^{18.5}\,\text{eV}$, $10^{19.0}\,\text{eV}$, and $10^{19.5}\,\text{eV}$
(1000 for each configuration).

Figure 10: $\chi^{2}$-distributions for the Greisen and the GH functions
fitted to simulated Epos-LHC air showers with primary energies of
$10^{19.5}\,\text{eV}$ (_left_), and fitted to simulated QgsjetII-04 air
showers with primary energies of $10^{18.5}\,\text{eV}$ (_right_). The
distributions for the Greisen (GH) functions are shown as a full (dashed)
line, for each of the three primary particles. The respective mean of each
distribution is indicated by a vertical line.

Figure 11: Two-dimensional distribution of the best-fit values of $X_{1}$ and
$X_{0}$ using the Greisen function to fit simulated longitudinal profiles
(_left_), and $X_{1}$ and $\Lambda$ using the GH function (_right_). The
showers were simulated with primary energies of $10^{18.5}\,\text{eV}$,
$10^{19}\,\text{eV}$, and $10^{19.5}\,\text{eV}$, using the Sibyll2.3d model
of hadronic interactions. The curved lines show the estimated $1\sigma$ extent
of the respective distributions.  Figure 12: Two-dimensional distribution of
the best-fit values of $X_{\text{max}}^{19}$ and $N_{\text{max}}^{19}$ using
the GH function to fit simulated longitudinal profiles of showers with primary
energies of $10^{18.5}\,\text{eV}$, $10^{19}\,\text{eV}$, and
$10^{19.5}\,\text{eV}$. All showers were simulated using the Sibyll2.3d model
of hadronic interactions. The curved lines show the estimated $1\sigma$ extent
of the respective distributions.