# Rogue wave patterns associated with Adler-Moser polynomials in the nonlinear
Schrödinger equation

Bo Yang1, Jianke Yang2 1 School of Mathematics and Statistics, Ningbo
University, Ningbo 315211, China
2 Department of Mathematics and Statistics, University of Vermont, Burlington,
VT 05401, U.S.A

###### Abstract

We report new rogue wave patterns in the nonlinear Schrödinger equation. These
patterns include heart-shaped structures, fan-shaped sectors, and many others,
that are formed by individual Peregrine waves. They appear when multiple
internal parameters in the rogue wave solutions get large. Analytically, we
show that these new patterns are described asymptotically by root structures
of Adler-Moser polynomials through a dilation. Since Adler-Moser polynomials
are generalizations of the Yablonskii-Vorob’ev polynomial hierarchy and
contain free complex parameters, these new rogue patterns associated with
Adler-Moser polynomials are much more diverse than previous rogue patterns
associated with the Yablonskii-Vorob’ev polynomial hierarchy. We also compare
analytical predictions of these patterns to true solutions and demonstrate
good agreement between them.

##  I Introduction

Rogue waves are unusually large, mysterious, and suddenly appearing surface
water waves that can be dangerous even to large ships Peli_book2009 . Their
counterparts in optics and many other physical fields have also been reported
Solli2007 ; Wabnitz2017 ; RomeroRos2022 . Due to their unexpected nature and
potential damage, rogue waves have been heavily studied in the physical and
mathematical communities in recent years. Physically, a lot of laboratory
experiments on rogue waves have been performed in diverse fields such as
optical fibers Kibler2010 ; Xu2019 ; Frisquet2016 ; Baronio2018 , water tanks
Chabchoub2011 ; Chabchoub2012a ; Chabchoub2012b ; Chabchoub2013 , superfluid
helium Helium2008 , plasma Bailung2011 ; Plasma2016 , and Bose-Einstein
condensates Engles2023 . Mathematically, rogue wave studies have been greatly
facilitated by the fact that, many integrable equations that govern diverse
physical processes, such as the nonlinear Schrödinger (NLS) equation for wave
packet evolution in deep water, optical fibers, plasma, and Bose-Einstein
condensates Benney ; Ablowitzbook ; Hasegawabook ; Saito1984 ; Gross ;
Pitaevskii ; Kevrekidis2022 , and the Manakov system for light transmission in
randomly birefringent optical fibers Bergano , admit explicit solutions that
exhibit rogue-wave characteristics. The first such solution was reported by
Peregrine in the NLS equation Peregrine1983 . This Peregrine solution starts
from a slightly perturbed constant-amplitude background wave. Then, it
develops a localized peak that is three times the height of the background.
Afterwards, this peak decays and merges to the background again. This is a
rogue wave since its transient peak is much higher than its original wave
amplitude, and this peak appears and disappears unexpectedly. It turns out
that the Peregrine solution is just the simplest (fundamental) rogue wave in
the NLS equation. More intricate NLS rogue waves that could reach even higher
transient peak amplitudes or multiple peaks were later discovered AAS2009 ;
DGKM2010 ; AKA2011 ; KAAN2011 ; GLML2012 ; OhtaJY2012 ; Miller2019 . In
addition, rogue waves were also derived for many other integrable systems such
as the Manakov system BDCW2012 ; ManakovDark ; LingGuoZhaoCNLS2014 ;
Chen_Shihua2014 ; Chen_Shihua2015 ; ZhaoGuoLingCNLS2016  and so on. These
analytical expressions of rogue waves shedded much light on the intricate
rogue wave dynamics and guided their observations in laboratory experiments
Kibler2010 ; Xu2019 ; Frisquet2016 ; Baronio2018 ; Chabchoub2011 ;
Chabchoub2012a ; Chabchoub2012b ; Chabchoub2013 ; Helium2008 ; Bailung2011 ;
Engles2023 .

Pattern formation in rogue waves is an important question, since such
information allows prediction of later rogue wave shapes from earlier rogue
wave forms. One of the simplest rogue patterns is a rogue triplet, which
comprises three fundamental rogue waves forming a triangle in the space-time
plane. Such rogue triplets have been reported theoretically in many integrable
equations DGKM2010 ; AKA2011 ; OhtaJY2012 ; Miller2019 ; Chen_Shihua2014 ;
Chen_Shihua2015  and observed experimentally in both water tanks and randomly
birefringent optical fibers Chabchoub2013 ; Baronio2018 . Beyond rogue
triplets, more sophisticated rogue patterns comprising Peregrine waves forming
shapes such as pentagons, heptagons, and rings have also been reported for the
NLS equation KAAN2011 ; HeFokas ; KAAN2013 ; YangNLS2021 . In addition, it was
revealed in YangNLS2021  that those sophisticated patterns would reliably
arise when one of the internal parameters in NLS rogue wave solutions gets
large, and their shapes are predicted asymptotically by root structures of the
Yablonskii-Vorob’ev polynomial hierarchy through a dilation and rotation (each
nonzero root of the hierarchy predicts the spatial-temporal location of a
Peregrine wave). In YangNLS2021b , it was further revealed that those NLS
rogue patterns associated with root structures of the Yablonskii-Vorob’ev
polynomial hierarchy are universal and would arise in many other integrable
systems, such as the derivative NLS equations, the Boussinesq equation, and
the Manakov system.

Given the importance of the NLS equation in diverse physical disciplines, an
important question is whether this equation admits other types of rogue
patterns. If so, what special polynomials can be used to predict those
patterns?

In this paper, we show that the NLS equation does admit many new rogue
patterns. These patterns arise when _multiple_ internal parameters in rogue
wave solutions get large, and their shapes are predicted asymptotically by
root structures of Adler-Moser polynomials through a simple dilation. Adler-
Moser polynomials are generalizations of the Yablonskii-Vorob’ev polynomial
hierarchy, and their root structures are much more diverse than those of the
Yablonskii-Vorob’ev polynomial hierarchy. As a consequence, we report many new
rogue patterns in the NLS equation, such as heart-shaped structures, fan-
shaped sectors, and others. Our analytical predictions for these new patterns
based on root structures of Adler-Moser polynomials are compared to true rogue
solutions, and excellent agreement is demonstrated.

##  II Preliminaries

The NLS equation is

$\textrm{i}u_{t}+\frac{1}{2}u_{xx}+|u|^{2}u=0.$ (1)

This equation governs nonlinear wave packet evolution in numerous physical
systems such as deep water, optical fibers, plasma, and two-component Bose-
Einstein condensates Benney ; Ablowitzbook ; Hasegawabook ; Saito1984 ;
Kevrekidis2022 .

###  II.1 General bilinear rogue wave solutions

Rogue wave solutions in this NLS equation have been derived by various methods
before DGKM2010 ; GLML2012 ; OhtaJY2012 ; Miller2019 . The most explicit forms
of these rogue wave solutions are the ones that were derived by the bilinear
method in OhtaJY2012  and then further simplified in YangNLS2021 . These
simplified bilinear rogue wave solutions of the $N$-th order are

$u_{N}(x,t)=\frac{\sigma_{1}}{\sigma_{0}}e^{\textrm{i}t},$ (2)
$\sigma_{n}=\det_{\begin{subarray}{l}1\leq i,j\leq
N\end{subarray}}\left(\begin{array}[]{c}\phi_{2i-1,2j-1}^{(n)}\end{array}\right),$
(3) $\phi_{i,j}^{(n)}=\sum_{\nu=0}^{\min(i,j)}\frac{1}{4^{\nu}}\hskip
1.70709ptS_{i-\nu}(\textbf{\emph{x}}^{+}(n)+\nu\textbf{\emph{s}})\hskip
1.70709ptS_{j-\nu}(\textbf{\emph{x}}^{-}(n)+\nu\textbf{\emph{s}}),$ (4)
$x_{1}^{\pm}=x\pm\textrm{i}t\pm n,\quad x_{2k}^{\pm}=0,$ (5)
$x_{2k+1}^{+}=\frac{x+2^{2k}(\textrm{i}t)}{(2k+1)!}+a_{2k+1},\quad
x_{2k+1}^{-}=(x_{2k+1}^{+})^{*},$ (6)

where $S_{k}(\mbox{\boldmath$x$})$ with
$\emph{{x}}=\left(x_{1},x_{2},\ldots\right)$ are Schur polynomials defined by
the generating function

$\sum_{k=0}^{\infty}S_{k}(\mbox{\boldmath$x$})\epsilon^{k}=\exp\left(\sum_{k=1}^{\infty}x_{k}\epsilon^{k}\right),$
(7)

$\textbf{\emph{s}}=(0,s_{2},0,s_{4},\cdots)$ are coefficients from the
expansion

$\sum_{j=1}^{\infty}s_{j}\lambda^{j}=\ln\left[\frac{2}{\lambda}\tanh\left(\frac{\lambda}{2}\right)\right],$
(8)

the asterisk * represents complex conjugation, and
$a_{3},a_{5},\cdots,a_{2N-1}$ are free irreducible complex parameters which
control the shape of this rogue wave solution.

When $N=1$, the above solution is $u_{1}(x,t)=\hat{u}_{1}(x,t)\hskip
1.13791pte^{\textrm{i}t}$, where

$\hat{u}_{1}(x,t)=1-\frac{4(1+2\textrm{i}t)}{1+4x^{2}+4t^{2}}.$ (9)

This is the fundamental rogue wave in the NLS equation that was discovered by
Peregrine in Peregrine1983 .

###  II.2 Adler-Moser polynomials and their root structures

Adler-Moser polynomials were proposed by Adler and Moser Adler_Moser1978 , who
expressed rational solutions of the Korteweg-de Vries equation in terms of
those polynomials. In a different context of point vortex dynamics, it was
discovered unexpectedly that the zeros of these polynomials also form
stationary vortex configurations when the vortices have the same strength but
positive or negative orientations, and the numbers of those positive and
negative vortices are consecutive triangular numbers Aref2007FDR ;
Clarkson2009 .

Adler-Moser polynomials $\Theta_{N}(z)$ can be written as a determinant
Clarkson2009

$\Theta_{N}(z)=c_{N}\left|\begin{array}[]{cccc}\theta_{1}(z)&\theta_{0}(z)&\cdots&\theta_{2-N}(z)\\\
\theta_{3}(z)&\theta_{2}(z)&\cdots&\theta_{4-N}(z)\\\
\vdots&\vdots&\vdots&\vdots\\\
\theta_{2N-1}(z)&\theta_{2N-2}(z)&\cdots&\theta_{N}(z)\end{array}\right|,$
(10)

where $\theta_{k}(z)$ are Schur polynomials defined by

$\sum_{k=0}^{\infty}\theta_{k}(z)\epsilon^{k}=\exp\left(z\epsilon+\sum_{j=1}^{\infty}\kappa_{j}\epsilon^{2j+1}\right),$
(11)

$\theta_{k}(z)\equiv 0$ if $k<0$, $c_{N}=\prod_{j=1}^{N}(2j-1)!!$, and
$\kappa_{j}\hskip 1.13791pt(j\geq 1)$ are arbitrary complex constants. Note
that our $\kappa_{j}$ constant is slightly different from that in Clarkson2009
by a factor of $-1/(2j+1)$, and this different parameter definition will be
more convenient for our purpose. The determinant in (10) is a Wronskian since
we can see from Eq. (11) that $\theta^{\prime}_{k}(z)=\theta_{k-1}(z)$, where
the prime denotes differentiation. In addition, these $\Theta_{N}(z)$
polynomials are monic with degree $N(N+1)/2$, which can be seen by noticing
that the highest $z$ term of $\theta_{k}(z)$ is $z^{k}/k!$, and the
determinant in (10) with $\theta_{k}(z)$ replaced by its highest $z$ term can
be explicitly calculated as $z^{N(N+1)/2}$ OhtaJY2012 . Adler-Moser
polynomials reduce to the Yablonskii-Vorob’ev polynomial hierarchy when all
$\kappa_{j}$ constants are set as zero except for one of them YangNLS2021 .
Thus, we can view Adler-Moser polynomials as generalizations of the
Yablonskii-Vorob’ev polynomial hierarchy.

The first few Adler-Moser polynomials are

$\displaystyle\Theta_{1}(z)=z,$
$\displaystyle\Theta_{2}(z)=z^{3}-3\kappa_{1},$
$\displaystyle\Theta_{3}(z)=z^{6}-15\kappa_{1}z^{3}+45\kappa_{2}z-45\kappa_{1}^{2},$
$\displaystyle\Theta_{4}(z)=z^{10}-45\kappa_{1}z^{7}+315\kappa_{2}z^{5}-1575\kappa_{3}z^{3}$
$\displaystyle+4725\kappa_{1}\kappa_{2}z^{2}-4725\kappa_{1}^{3}z-4725\kappa_{2}^{2}+4725\kappa_{1}\kappa_{3}.\mbox{\hskip
28.45274pt}$

Root structures of Adler-Moser polynomials are important to us, since we will
link them to rogue wave patterns in the later text. Due to the free complex
parameters $\\{\kappa_{j}\\}$ in them, their root structures will be
understandably very diverse — much more diverse than root structures of
Yablonskii-Vorob’ev hierarchy polynomials. Indeed, when setting all
$\\{\kappa_{j}\\}$ as zero except for one of them, we get root structures of
Yablonskii-Vorob’ev hierarchy polynomials which are in the shape of triangles,
pentagons, heptagons, and so on. When we continuously change those
$\\{\kappa_{j}\\}$ values, we will get root structures which smoothly deform
from one type of Yablonskii-Vorob’ev root structure to another, such as from a
triangle to a pentagon. In this process, uncountably infinite new root shapes
will be generated. These roots are generically simple roots. Indeed, if a root
happens to be a multiple root, it will split into simple roots when the
complex parameters $\\{\kappa_{j}\\}$ are slightly perturbed. For this reason,
we will focus on the case when all roots of $\Theta_{N}(z)$ are simple in this
article. In this case, $\Theta_{N}(z)$ will have $N(N+1)/2$ roots.

Of the uncountably infinite root structures of Adler-Moser polynomials, we
illustrate only three of them for brevity. These three samples are for
$\Theta_{5}(z;\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4})$, with three sets
of $(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4})$ values as

$(\textrm{i},\textrm{i},\textrm{i},\textrm{i}),\hskip
2.84544pt(5\textrm{i}/3,\textrm{i},5\textrm{i}/7,-5\textrm{i}/9),\hskip
2.84544pt(1,1,1,1).$ (12)

Their root structures are displayed in Fig. 1(a, b, c), respectively. In these
panels, every root is a simple root. The (a) panel shows a heart-shaped
structure, (b) shows a fan-shaped circular sector, and (c) shows a two-arc
structure combined with a triangle.

Figure 1: Three sample root structures of $\Theta_{5}(z)$ for parameter values
$(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4})$ given in Eq. (12). In all
panels, $-7\leq\mbox{Re}(z),\mbox{Im}(z)\leq 7$.

##  III Analytical predictions for rogue patterns with multiple large internal
parameters

Rogue wave solutions $u_{N}(x,t)$ in Eq. (2) contain $N-1$ free internal
complex parameters $a_{3},a_{5},\cdots,a_{2N-1}$. If only one of those
parameters is large, then the resulting rogue pattern is predicted by root
structures of Yablonskii-Vorob’ev hierarchy polynomials, see YangNLS2021 . In
this section, we consider patterns of these rogue solutions when _multiple_ of
these internal parameters are large.

Specifically, suppose parameters $a_{3},a_{5},\cdots,a_{2N-1}$ in $u_{N}(x,t)$
are of the following form

$a_{2j+1}=\kappa_{j}\hskip 1.13791ptA^{2j+1},\quad 1\leq j\leq N-1,$ (13)

where $A\gg 1$ is a large positive constant, and
$(\kappa_{1},\kappa_{2},\dots,\kappa_{N-1})$ are $O(1)$ complex constants not
being all zero. Suppose also that roots of the Adler-Moser polynomial
$\Theta_{N}(z)$ with parameters $\\{\kappa_{j}\\}$ are all simple. Then, our
analytical prediction on the pattern of this rogue wave solution $u_{N}(x,t)$
is given by the following theorem.

Theorem 1. _If all roots of $\Theta_{N}(z)$ are simple, then the $N$-th order
rogue wave $u_{N}(x,t)$ in Eq. (2) with its internal large parameters
$(a_{3},a_{5},\cdots,a_{2N-1})$ as given by Eq. (13) would asymptotically
split into $N(N+1)/2$ fundamental (Peregrine) rogue waves of the form
$\hat{u}_{1}(x-\hat{x}_{0},t-\hat{t}_{0})\hskip 1.42271pte^{\textrm{i}t}$,
where $\hat{u}_{1}(x,t)$ is given in Eq. (9), and positions
$(\hat{x}_{0},\hat{t}_{0})$ of these Peregrine waves are given by_

$\hat{x}_{0}+\textrm{i}\hskip 1.42271pt\hat{t}_{0}=z_{0}A,$ (14)

_with $z_{0}$ being every one of the $N(N+1)/2$ simple roots of
$\Theta_{N}(z)$. The error of this Peregrine wave approximation is
$O(A^{-1})$. Expressed mathematically, when
$(x-\hat{x}_{0})^{2}+(t-\hat{t}_{0})^{2}=O(1)$, we have the following solution
asymptotics_

$u_{N}(x,t;a_{3},a_{5},\cdots,a_{2N-1})=\hat{u}_{1}(x-\hat{x}_{0},t-\hat{t}_{0})\hskip
1.42271pte^{\textrm{i}t}+O\left(A^{-1}\right).$

_When $(x,t)$ is not in the neighborhood of any of these Peregrine waves,
$u_{N}(x,t)$ would asymptotically approach the constant-amplitude background
$e^{\textrm{i}t}$ as $A\to+\infty$._

This theorem indicates that the rogue pattern is asymptotically a simple
dilation of the root structure of the underlying Adler-Moser polynomial by a
factor of $A$, with each root predicting the location of a Peregrine wave in
the $(x,t)$ plane according to Eq. (14). Thus, this theorem establishes a
direct connection between rogue patterns and root structures of Adler-Moser
polynomials.

One may notice that in the present case of multiple large parameters, the
rogue pattern is a simple dilation of the root structure of an Adler-Moser
polynomial, while in the previous case of a single large parameter as studied
in YangNLS2021 , the rogue pattern was a dilation _and rotation_ of the root
structure of a Yablonskii-Vorob’ev hierarchy polynomial. The reason our
current rogue pattern does not involve rotation to the root structure is that,
the Adler-Moser polynomial contains free complex constants $\\{\kappa_{j}\\}$,
which automatically put its root structure in proper orientation to match the
rogue pattern. Comparatively, a Yablonskii-Vorob’ev hierarchy polynomial does
not contain such free complex constants, and thus the orientation of its root
structure is fixed. In this case, in order for its root structure to match the
orientation of the rogue wave, a proper rotation is needed.

##  IV Numerical confirmation

Now, we numerically verify Theorem 1 by comparing its predictions with true
rogue-wave solutions. This comparison will be done only for fifth-order rogue
waves $u_{5}(x,t)$ for brevity. Such fifth-order solutions have internal
complex parameters $(a_{3},a_{5},a_{7},a_{9})$.

We will do this comparison on three examples. Internal parameter values in
these three examples are of the form (13) with $A=5$, which is large as
desired, and their $(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4})$ values are
given in Eq. (12). These $\kappa_{j}$ values are used since root structures of
Adler-Moser polynomials $\Theta_{5}(z)$ for these values have been displayed
in Fig. 1. For these three sets of internal parameters, true rogue wave
solutions are plotted in the upper three panels of Fig. 2, respectively. It is
seen that each panel comprises 15 lumps (Peregrine waves) in the $(x,t)$
plane. In the first panel, these 15 Peregrine waves form a heart-shaped
structure, with another mini-heart in its interior. In the second panel, these
15 Peregrine waves form a fan-shaped structure. In the third panel, these 15
Peregrine waves form two vertically-oriented arcs plus a smaller triangle on
their right side.

Our analytical predictions $|u_{5}^{(p)}(x,t)|$ for these rogue waves from
Theorem 1 can be assembled into a simple formula,

$\left|u_{5}^{(p)}(x,t)\right|=1+\sum_{j=1}^{15}\left(\left|\hat{u}_{1}(x-\hat{x}_{0}^{(j)},t-\hat{t}_{0}^{(j)})\right|-1\right),$
(15)

where $\hat{u}_{1}(x,t)$ is the Peregrine wave given in (9), and their
positions $(\hat{x}_{0}^{(j)},\hat{t}_{0}^{(j)})$ are given by (14) with
$z_{0}$ being every one of the $N(N+1)/2=15$ simple roots of the Adler-Moser
polynomial $\Theta_{5}(z)$. These predicted solutions for the same
$(a_{3},a_{5},a_{7},a_{9})$ values as in the true solutions are plotted in the
lower three panels of Fig. 2. When compared to the root structures of the
Adler-Moser polynomial $\Theta_{5}(z)$ in Fig. 1, our predicted rogue patterns
in these lower panels are obviously a simple dilation of those root
structures, by a factor of $A=5$, with each root replaced by a Peregrine wave,
as Theorem 1 says.

When comparing the true rogue solutions in the upper row to their analytical
predictions in the lower row, we can clearly see that they agree with each
other very well. In fact, one can hardly notice the difference between them,
which is an indication that our prediction in Theorem 1 is highly accurate.

Quantitatively, we have also measured the error of our analytical predictions
versus the $A$ value, similar to what we did in Fig. 5 of Ref. YangNLS2021 .
That error analysis confirmed that the error does decay in proportion to
$A^{-1}$, as Theorem 1 predicts. Thus, Theorem 1 is fully confirmed
numerically. Details of this quantitative comparison are omitted here for
brevity.

Figure 2: Comparison between true rogue solutions $|u(x,t)|$ (upper row) and
their analytical predictions (lower row) for $N=5$ and $A=5$. From left to
right columns:
$(\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4})=(\textrm{i},\textrm{i},\textrm{i},\textrm{i}),\
(5\textrm{i}/3,\textrm{i},5\textrm{i}/7,-5\textrm{i}/9),\ (1,1,1,1)$. In all
panels, $-30\leq x,t\leq 30$.

##  V Proof of Theorem 1

In this section, we prove the analytical predictions on NLS rogue patterns in
Theorems 1. The main idea of our proof resembles that in Ref. YangNLS2021  for
a single large internal parameter case.

To derive the large-parameter asymptotics of the rogue wave solution
$u_{N}(x,t)$ in Eq. (2), we need asymptotic expressions for the determinant
$\sigma_{n}$ in Eq. (3). For this purpose, we first use determinant identities
and the Laplace expansion to rewrite $\sigma_{n}$ as OhtaJY2012

$\displaystyle\sigma_{n}=\sum_{0\leq\nu_{1}<\nu_{2}<\cdots<\nu_{N}\leq
2N-1}\det_{1\leq i,j\leq
N}\left[\frac{1}{2^{\nu_{j}}}S_{2i-1-\nu_{j}}(\textbf{\emph{x}}^{+}(n)+\nu_{j}\textbf{\emph{s}})\right]$
$\displaystyle\times\det_{1\leq i,j\leq
N}\left[\frac{1}{2^{\nu_{j}}}S_{2i-1-\nu_{j}}(\textbf{\emph{x}}^{-}(n)+\nu_{j}\textbf{\emph{s}})\right],$
(16)

where $S_{k}\equiv 0$ if $k<0$.

When internal parameters $(a_{3},a_{5},\cdots,a_{2N-1})$ are of the form (13)
with $A\gg 1$, and $x,t=O(A)$ or smaller, we have

$\displaystyle
S_{k}(\textbf{\emph{x}}^{+}(n)+\nu\textbf{\emph{s}})=S_{k}\left(x_{1}^{+},\nu
s_{2},x_{3}^{+},\nu s_{4},\cdots\right)$
$\displaystyle=S_{k}\left(x_{1}^{+},0,\kappa_{1}A^{3},0,\kappa_{2}A^{5},\cdots\right)\left[1+O(A^{-2})\right]$
$\displaystyle=S_{k}(\hat{\textbf{v}})\left[1+O(A^{-2})\right],$ (17)

where
$\hat{\textbf{v}}=\left(x+\textrm{i}t+n,0,\kappa_{1}A^{3},0,\kappa_{2}A^{5},0,\cdots\right)$.
From the definition (7) of Schur polynomials, one can see that the polynomial
$S_{k}(\hat{\textbf{v}})$ is related to $\theta_{k}(z)$ in (11) as

$S_{k}(\hat{\textbf{v}})=A^{k}\theta_{k}(\hat{z}),$ (18)

where $\hat{z}\equiv A^{-1}(x+\textrm{i}t+n)$.

The dominant contribution in the Laplace expansion (V) of $\sigma_{n}$ comes
from two index choices, $\nu=(0,1,\cdots,N-1)$, and $\nu=(0,1,\cdots,N-2,N)$.

With the first index choice, in view of Eqs. (17)-(18), the determinant
involving $\textbf{\emph{x}}^{+}(n)$ inside the summation of (V) is
asymptotically

$\alpha\hskip
1.70709ptA^{\frac{N(N+1)}{2}}\Theta_{N}(\hat{z})\left[1+O\left(A^{-2}\right)\right],$
(19)

where $\alpha=2^{-N(N-1)/2}c_{N}^{-1}$. Let us define
$(\hat{x}_{0},\hat{t}_{0})$ by Eq. (14), i.e.,
$z_{0}=A^{-1}(\hat{x}_{0}+\textrm{i}\hat{t}_{0})$, where $z_{0}$ is a simple
root of the Adler-Moser polynomial $\Theta_{N}(z)$. Then, when $(x,t)$ is in
the $O(1)$ neighborhood of $(\hat{x}_{0},\hat{t}_{0})$, we expand
$\Theta_{N}(\hat{z})$ around $\hat{z}=z_{0}$. Recalling $\Theta_{N}(z_{0})=0$,
we get

$\Theta_{N}(\hat{z})=A^{-1}\left[(x-\hat{x}_{0})+\textrm{i}(t-\hat{t}_{0})+n\right]\Theta^{\prime}_{N}(z_{0})\left[1+O\left(A^{-1}\right)\right].$

Inserting this equation into (19), the determinant involving
$\textbf{\emph{x}}^{+}(n)$ inside the summation of (V) becomes

$\left[(x-\hat{x}_{0})+\textrm{i}(t-\hat{t}_{0})+n\right]\hskip
1.70709pt\alpha\hskip
1.70709ptA^{\frac{N(N+1)}{2}}\Theta_{N}^{\prime}(z_{0})\left[1+O\left(A^{-1}\right)\right].$

Similarly, the determinant involving $\textbf{\emph{x}}^{-}(n)$ inside this
summation becomes

$\left[(x-\hat{x}_{0})-\textrm{i}(t-\hat{t}_{0})-n\right]\hskip
1.70709pt\alpha\hskip
1.70709ptA^{\frac{N(N+1)}{2}}\Theta_{N}^{\prime}(z_{0}^{*})\left[1+O\left(A^{-1}\right)\right].$

Next, we consider the contribution from the second index choice of
$\nu=(0,1,\cdots,N-2,N)$. For this index choice, the determinant involving
$\textbf{\emph{x}}^{+}(n)$ inside the summation of (V) becomes

$\frac{1}{2}\alpha\hskip
1.70709ptA^{\frac{N(N+1)-2}{2}}\Theta_{N}^{\prime}(\hat{z})\left[1+O\left(A^{-2}\right)\right].$

When $(x,t)$ is in the $O(1)$ neighborhood of $(\hat{x}_{0},\hat{t}_{0})$, the
above term is asymptotically equal to

$\frac{1}{2}\alpha\hskip
1.70709ptA^{\frac{N(N+1)-2}{2}}\Theta_{N}^{\prime}(z_{0})\left[1+O\left(A^{-1}\right)\right].$

Similarly, the determinant involving $\textbf{\emph{x}}^{-}(n)$ inside the
summation of (V) becomes

$\frac{1}{2}\alpha\hskip
1.70709ptA^{\frac{N(N+1)-2}{2}}\Theta_{N}^{\prime}(z_{0}^{*})\left[1+O\left(A^{-1}\right)\right].$

Summarizing the above two dominant contributions in the Laplace expansion (V),
we find that

$\displaystyle\sigma_{n}(x,t)=\alpha^{2}\hskip
1.70709pt\left|\Theta_{N}^{\prime}(z_{0})\right|^{2}A^{N(N+1)-2}$ (20)
$\displaystyle\times\left[\left(x-\hat{x}_{0}\right)^{2}+\left(t-\hat{t}_{0}\right)^{2}-2\textrm{i}n\left(t-\hat{t}_{0}\right)-n^{2}+\frac{1}{4}\right]$
$\displaystyle\times\left[1+O\left(A^{-1}\right)\right].$

Since the root $z_{0}$ has been assumed simple,
$\Theta_{N}^{\prime}(z_{0})\neq 0$. Thus, the above leading-order asymptotics
for $\sigma_{n}(x,t)$ does not vanish. Therefore, when $A$ is large and
$(x,t)$ in the O(1) neighborhood of $\left(\hat{x}_{0},\hat{t}_{0}\right)$, we
get from (20) that

$\displaystyle
u_{N}(x,t)=\frac{\sigma_{1}}{\sigma_{0}}e^{\textrm{i}t}=e^{\textrm{i}t}\left(1-\frac{4[1+2\textrm{i}(t-\hat{t}_{0})]}{1+4(x-\hat{x}_{0})^{2}+4(t-\hat{t}_{0})^{2}}\right)$
$\displaystyle\hskip 39.83368pt+O\left(A^{-1}\right),$

which is a Peregrine wave $\hat{u}_{1}(x-\hat{x}_{0},t-\hat{t}_{0})\hskip
1.42271pte^{\textrm{i}t}$, and the error of this Peregrine prediction is
$O\left(A^{-1}\right)$. Theorem 1 is then proved.

##  VI Conclusions and discussions

In this paper, we have reported many new rogue patterns in the NLS equation
which are predicted by root structures of new special polynomials.
Specifically, we have shown that when multiple internal parameters in the
rogue wave solutions are large, many new rogue patterns would arise, including
heart-shaped structures, fan-shaped structures, and others. Analytically,
these rogue patterns are determined by root structures of Adler-Moser
polynomials. If all roots of the Adler-Moser polynomial are simple, then the
rogue pattern is simply a dilation of the Adler-Moser-polynomial’s root
structure, with each root replaced by a Peregrine wave. Since Adler-Moser
polynomials contain free complex parameters, their root structures would be
very diverse. As a result, NLS rogue waves could assume much more varied
spatial-temporal patterns beyond those reported earlier.

Since NLS rogue waves have been observed in diverse physical systems
Kibler2010 ; Xu2019 ; Chabchoub2011 ; Chabchoub2012a ; Chabchoub2012b ;
Chabchoub2013 ; Helium2008 ; Bailung2011 ; Engles2023 , these new NLS rogue
patterns open up more varieties of rogue dynamics which could be verified in
experiments too.

The previous NLS rogue patterns associated with root structures of the
Yablonskii-Vorob’ev polynomial hierarchy were later found to be universal and
would appear in many other integrable systems YangNLS2021 ; YangNLS2021b . The
present rogue patterns associated with Adler-Moser polynomials are expected to
be universal as well, and they should arise in other integrable systems too
when multiple internal parameters in rogue waves of those integrable equations
are large. This prospect will be pursued in the near future.

Our rogue-pattern predictions in this article were made under the assumption
that all roots of the Adler-Moser polynomial are simple. A very interesting
question is what will happen if some roots of the Adler-Moser polynomial are
not simple. This question is still open that merits further studies.

## Acknowledgment

The work of B.Y. was supported in part by the National Natural Science
Foundation of China (GrantNo.12201326), and the work of J.Y. was supported in
part by the National Science Foundation (U.S.) under award number DMS-1910282.

## References

## References

  * (1) C. Kharif, E. Pelinovsky, and A. Slunyaev, _Rogue Waves in the Ocean_ (Springer, Berlin, 2009).
  * (2) D.R. Solli, C. Ropers, P. Koonath, and B. Jalali, Nature 450, 1054 (2007).
  * (3) S. Wabnitz (ed.), _Nonlinear Guided Wave Optics: A Testbed for Extreme Waves_ (IOP Publishing, Bristol, 2017).
  * (4) A. Romero-Ros, G. C. Katsimiga, S. I. Mistakidis, B. Prinari, G. Biondini, P. Schmelcher, and P. G. Kevrekidis, Phys. Rev. A 105, 053306 (2022).
  * (5) B. Kibler, J. Fatome, C. Finot, G. Millot, F. Dias, G. Genty, N. Akhmediev and J. M. Dudley, Nature Phys. 6, 790 (2010).
  * (6) G. Xu, K. Hammani, A. Chabchoub, J. M. Dudley, B. Kibler, and C. Finot, Phys. Rev. E 99, 012207 (2019).
  * (7) B. Frisquet, B. Kibler, P. Morin, F. Baronio, M. Conforti, G. Millot, and S. Wabnitz, Sci. Rep. 6, 20785 (2016).
  * (8) F. Baronio, B. Frisquet, S. Chen, G. Millot, S.Wabnitz, and B. Kibler, Phys. Rev. A 97, 013852 (2018).
  * (9) A. Chabchoub, N.P. Hoffmann, and N. Akhmediev, Phys. Rev. Lett. 106, 204502 (2011).
  * (10) A. Chabchoub, N. Hoffmann, M. Onorato, and N. Akhmediev, Phys. Rev. X 2, 011015 (2012).
  * (11) A. Chabchoub, N. Hoffmann, M. Onorato, A. Slunyaev, A. Sergeeva, E. Pelinovsky and N. Akhmediev, Phys. Rev. E 86, 056601 (2012).
  * (12) A. Chabchoub and N. Akhmediev, Phys. Lett. A 377, 2590 (2013).
  * (13) A.N. Ganshin, V.B. Efimov, G.V. Kolmakov, L.P. Mezhov-Deglin and P.V.E. McClintock, Phys. Rev. Lett. 101, 065303 (2008).
  * (14) H. Bailung, S.K. Sharma and Y. Nakamura, Phys. Rev. Lett. 107, 255005 (2011).
  * (15) Y.Y. Tsai, J.Y. Tsai and L. I, Nat. Phys. 12, 573 (2016).
  * (16) A. Romero-Ros, G.C. Katsimiga, S.I. Mistakidis, S. Mossman, G. Biondini, P. Schmelcher, P. Engels, and P.G. Kevrekidis, arXiv:2304.05951 [nlin.PS] (2023).
  * (17) D.J. Benney and A.C. Newell, J. Math. Phys. 46, 133 (1967).
  * (18) M.J. Ablowitz and H. Segur, _Solitons and the Inverse Scattering Transform_ (SIAM, Philadelphia, 1981).
  * (19) A. Hasegawa and Y. Kodama, _Solitons in Optical Communications_ (Clarendon Press, Oxford, 1995).
  * (20) M. Saito, S. Watanabe and H. Tanaca, J. Phys. Soc. Japan 53, 2304 (1984).
  * (21) E.P. Gross, Il Nuovo Cimento. 20, 454 (1961).
  * (22) L.P. Pitaevskii, Sov. Phys. JETP. 13, 451 (1961).
  * (23) A. Romero-Ros, G. C. Katsimiga, S. I. Mistakidis, B. Prinari, G. Biondini, P. Schmelcher, and P. G. Kevrekidis, Phys. Rev. A 105, 053306 (2022).
  * (24) S.G. Evangelides, L.F. Mollenauer, J.P. Gordon, and N.S. Bergano, J. Lightwave Technol. 10, 28 (1992).
  * (25) D.H. Peregrine, J. Aust. Math. Soc. B 25, 16 (1983).
  * (26) N. Akhmediev, A. Ankiewicz and J.M. Soto-Crespo, Phys. Rev. E 80, 026601 (2009).
  * (27) P. Dubard, P. Gaillard, C. Klein and V.B. Matveev, Eur. Phys. J. Spec. Top. 185, 247 (2010).
  * (28) A. Ankiewicz, D. Kedziora and N. Akhmediev, Phys. Lett. A 375, 2782 (2011).
  * (29) D.J. Kedziora, A. Ankiewicz and N. Akhmediev, Phys. Rev. E 84, 056611 (2011).
  * (30) B.L. Guo, L.M. Ling and Q.P. Liu, Phys. Rev. E 85, 026607 (2012).
  * (31) Y. Ohta and J. Yang, Proc. R. Soc. Lond. A 468, 1716 (2012).
  * (32) D. Bilman and P.D. Miller, Commun. Pure Appl. Math. 72, 1722 (2019).
  * (33) F. Baronio, A. Degasperis, M. Conforti and S. Wabnitz, Phys. Rev. Lett. 109, 044102 (2012).
  * (34) F. Baronio, M. Conforti, A. Degasperis, S. Lombardo, M. Onorato and S. Wabnitz, Phys. Rev. Lett. 113, 034101 (2014).
  * (35) L. Ling, B. Guo and L. Zhao, Phys. Rev. E 89, 041201(R) (2014).
  * (36) S. Chen, J. M. Soto-Crespo, and P. Grelu, Opt. Express 22, 27632 (2014).
  * (37) S. Chen and D. Mihalache, J. Phys. A 48, 215202 (2015).
  * (38) L. Zhao, B. Guo and L. Ling, J. Math. Phys. 57, 043508 (2016).
  * (39) J.S. He, H.R. Zhang, L.H. Wang, K. Porsezian and A.S. Fokas, Phys. Rev. E 87, 052914 (2013).
  * (40) D.J. Kedziora, A. Ankiewicz and N. Akhmediev, Phys. Rev. E 88, 013207 (2013).
  * (41) B. Yang and J. Yang, Physica D 419, 132850 (2021).
  * (42) B. Yang and J. Yang, Physica D 425, 132958 (2021).
  * (43) M. Adler and J. Moser, Commun. Math. Phys. 61, 1 (1978).
  * (44) H. Aref, Fluid Dynam. Res. 39, 5 (2007).
  * (45) P.A. Clarkson, Stud. Appl. Math. 123, 37 (2009).


# Probing the cosmic ray energy spectrum at $10^{12}$–$10^{16}$ eV with two
HiSPARC scintillators

K. van Dam<EMAIL_ADDRESS>B. van Eijk  J.J.M. Steijger  Nikhef National
Institute for Subatomic Physics, Science Park 105, 1098 XG Amsterdam, The
Netherlands  Faculty of Science and Technology, University of Twente, 7500 AE
Enschede, The Netherlands

###### Abstract

The high school project on astrophysics research with cosmics (HiSPARC)
employs a large number of small detection stations that sample the footprint
of extensive cosmic ray air showers. The majority of these stations has two
0.5 ${\mathrm{m}}^{2}$ scintillation detectors. A new method is presented
which enables probing the cosmic ray flux with a single two-scintillator
station in five energy decades at $10^{12}$, $10^{13}$, $10^{14}$, $10^{15}$
and $10^{16}$ eV. The method is based on the energy dependence of the
distribution of the number of particles passing through a single detector. A
relatively short data taking period of approximately one month is sufficient
to probe this energy range. The flux values agree well with measurements by
other experiments. For the first time, the cosmic ray flux at $10^{12}$ and
$10^{13}$ eV is derived at sea level.

###### keywords:

HiSPARC , Cosmic rays , Scintillation detector , Cosmic ray energy spectrum

††journal: Astroparticle Physics

##  1 Introduction

Cosmic rays are charged particles from space with energies ranging from GeV to
hundreds of EeV. The majority of cosmic rays are protons (90%) and helium
nuclei (9%) with a remainder of heavier elements [1]. The cosmic ray flux
drops rapidly with increasing energy. The flux as function of energy ($E$) in
the range from GeV to $10^{14}$ eV is approximately given by [2]:

$I(E)\approx 1.8\cdot 10^{4}\ (E/1\ \text{GeV})^{-\alpha}\frac{\text{\\# of
nucleons}}{\text{m}^{2}\ \text{s sr GeV}}$ (1)

with $\alpha\approx 2.7$. At an energy of $\sim$$10^{15.5}$ eV the spectrum
steepens. This region is known as ’the knee’. At $\sim$$10^{18.5}$ eV the
slope of the spectrum becomes less steep, this is called ’the ankle’. Cosmic
rays with energies below $\sim$$10^{14}$ eV can be measured using high
altitude balloons or spacecraft. For higher energies this approach becomes
unfeasible because of the low cosmic ray flux, the limited collection area and
lifetime of space based instruments. Above $\sim$$10^{14}$ eV cosmic rays are
studied using the Earth’s atmosphere. When a cosmic ray enters the atmosphere
it will most likely interact with a nitrogen or oxygen nucleus. In this
interaction new particles may be produced which, in turn, will interact with
other atmospheric particles. The particle generation continues until there is
insufficient energy left for the creation of additional particles. This
phenomenon is known as the development of an extensive air shower (EAS). The
particle footprint of an EAS at the ground is sampled using detector arrays.
The direction and energy of the cosmic ray can be reconstructed from particle
arrival times and multiplicities. Also atmospheric radiation generated by the
EAS, such as atmospheric Cherenkov and fluorescence, and radio waves, can be
used to reconstruct the shower size and direction [1]. Cosmic rays with
energies below $\sim$$10^{14}$ eV do not result in large air shower
footprints. The particle densities are too small to reconstruct the direction
and energy of the primary cosmic ray.

If an even lower energy cosmic ray (e.g. $10^{11}$ eV) hits the Earth’s
atmosphere, the air shower dies out before reaching the ground. The majority
of shower particles are absorbed at high altitude. Only energetic remnants, of
which mainly muons (and neutrinos), are able to reach ground level. The flux
of perpendicularly incident muons at sea level is $\sim$70
${\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{sr}}^{-1}$ above
1 GeV [2]. If an energetic muon decays before reaching the ground, the decay
electron (positron) emits gamma radiation due to Bremsstrahlung. Energetic
gamma rays will interact with matter via pair production. The two processes
(Bremsstrahlung and pair production) generate an electromagnetic shower.
Figure 1 shows an example of such a shower in a $10^{11}$ eV proton CORSIKA
[3] simulation. The upper right figure shows the longitudinal profile, i.e.
the number of particles per type as function of atmospheric depth. The
original shower dies out and, at a depth $\sim$900
$\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$, a second, electromagnetic (e.m.),
shower appears. The bottom right plot shows the energy per particle type as
function of the atmospheric depth. In one of the first interactions, $\sim$40
GeV of the 100 GeV proton is forwarded to a single muon. At a height of 1.6 km
(855 $\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$) this muon decays and a large
fraction of its energy is transferred to the electron (positron) which
initiates an e.m. shower. The left figure shows the footprint created by the
shower at the ground. The particle density is very small. The position of the
shower core is shifted with respect to the original direction of the cosmic
ray proton due to absorption of the shower particles other than the muon, and
neutrino(s) escaping detection. Muon decay is the dominant source for
electrons at sea level. The total flux of perpendicular incident electrons and
positrons above 10 MeV is $\sim$30
${\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{sr}}^{-1}$ [2].

Figure 1: Example of a muon decay induced shower in a $10^{11}$ eV proton
CORSIKA [3] simulation. The upper right figure shows the longitudinal profile,
i.e. the number of particles per particle type as function of atmospheric
depth. The original shower dies out and, at a depth $\sim$900
$\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$, a second (very large) shower emerges.
The bottom right plot shows the energy per particle type as function of the
atmospheric depth. In one of the first interactions, $\sim$40 GeV of the
available 100 GeV is transferred to a single muon. At a height of 1.6 km (855
$\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$) this muon decays and a large part of
its energy is passed on to an electron (positron) which ignites an e.m.
shower. The left figure shows the e.m. shower footprint at the ground. The
position of the shower core is slightly shifted from the direction of the
primary proton as particles are absorbed and neutrino(s) remain undetected.

The high school project on astrophysics research with cosmics (HiSPARC [4])
employs a large number of small EAS detection stations. The stations are
predominantly hosted by high schools. Each station consists of two or four
scintillators connected to an electronics unit that digitizes the analogue PMT
signals. A station uses a GPS receiver for accurate timing and position
information. The majority of high schools employs a two-detector station. As
neighboring schools often also have HiSPARC equipment, multiple stations are
used to sample the same EAS. Data sets are combined to reconstruct the air
shower.

A new method is presented to enable the determination of the cosmic ray flux
at energies below $10^{14}$ eV using two scintillators at sea level. The
method does not require the full reconstruction of individual EASs but is
based on measuring the energy dependence of the distribution of the number of
particles passing through a detector.

First, the experimental setup is introduced after which the new method is
detailed. Finally, the results are compared with measurements by other
experiments.

##  2 HiSPARC station

A HiSPARC detector consists of a rectangular scintillator (100 cm $\times$ 50
cm $\times$ 2 cm) glued to a light-guide which is attached to a PMT. The
detector is made light tight with a thin aluminum foil and pond liner. The
assembly is placed inside a roof box and mounted on top of a roof (figure 2).
When a charged particle traverses the scintillator, it creates a light pulse
which is converted into an electric pulse by the PMT. This pulse is sampled
and digitized at 400 MHz. The PMT is calibrated such that single charged
particles generate pulses with an amplitude of $\sim$150 mV. A 30 mV threshold
rejects PMT noise and the response to low energy gamma rays. The detector
response is simulated using GEANT4 [5] and verified against experimental data.
The detection efficiency for minimum ionizing particles (MIP) is $\sim$100%.
The detector has a much lower sensitivity for gamma rays in EASs ($\sim$10%)
[4].

To discriminate against single particle sources, HiSPARC applies the
’coincidence method’. By using two (or four) detectors and selecting
coincident PMT pulses (i.e. with a time difference smaller than 1.5
$\mathrm{\SIUnitSymbolMicro s}$) the majority of random coincidences are
rejected. If two or more PMT pulses exceed the threshold within the trigger
window, the pulses are stored. Random coincidences occur due to gamma rays
from decays of nearby environmental radionuclides [6] and muons generated by
low energy cosmic rays of which the other components of the air shower are
absorbed in the atmosphere, etc. A rare event like the one in figure 1 has a
chance to result in a trigger.

The distance between the detectors in a two-detector station varies per setup.
In this paper a detector separation of 4.95 meter is used. In a four-detector
station two electronics units are combined in a master-slave configuration.
Similar trigger conditions apply as for a two-detector station. In a four-
detector station the direction can be triangulated with a resolution of
$\sim$6$\mathrm{\SIUnitSymbolDegree}$. An estimate of the shower energy can
also be made. An extensive description is presented in [4].

Figure 2: A two-detector HiSPARC station with two scintillators placed in roof
boxes mounted on top of a roof. The detector separation is typically
$\sim$5$-$10 m. In between the two detectors a GPS antenna is placed for
accurate timing and location information. Photo courtesy of Arne de Laat.

##  3 Pulse integral distribution (PID)

The value of the PMT pulse integral, i.e. the area under the pulse, is
directly proportional to the number of scintillation photons that reach the
PMT, and thus to the number of particles traversing the detector (MIPs). The
energy loss of a MIP follows Landau’s theory [7]. The Landau distribution has
a peak at the most probable energy loss with a pronounced tail towards higher
energies. If multiple particles simultaneously traverse a detector, the energy
loss is described by a sum of Landau distributions. The shape of the pulse
integral distribution is then determined by the number of MIPs that traverse
the detector within the trigger time window ($1.5$ $\mathrm{\SIUnitSymbolMicro
s}$). The energy of the primary cosmic ray is directly proportional to the
number of EAS particles in the footprint at the ground.

###  3.1 Analyzing the PID

In a first approach all cosmic rays are assumed to stem from the zenith. The
particle densities can be approximated by the lateral density profile from
AGASA [8]:

$\displaystyle S(r)\propto$
$\displaystyle\left(\frac{r}{r_{\textup{M}}}\right)^{-1.2}\left(1+\frac{r}{r_{\textup{M}}}\right)^{-(\eta-1.2)}$
$\displaystyle\times\left[1+\left(\frac{r}{1000}\right)^{2}\right]^{-0.6}$ (2)

with distance to the shower core $r$ in meter, $\eta=3.84$ and Molière radius
$r_{\textup{M}}=91.6$ m. The energy associated with this lateral density
profile is given by [9]:

$\textup{E}=a\cdot 10^{17}\cdot(S(600))^{b}\quad[\textup{eV}]$ (3)

with $a=2.03$ and $b=1.0$. The dashed lines in figure 3 show EAS particle
densities as function of distance to the shower core for primary energies of
$10^{13}$ (blue), $10^{14}$ (red), $10^{15}$ (green) and $10^{16}$ eV
(orange). Obviously, as the EAS energy increases, the size of the footprint
increases as well. Next, the particle density per unit area is rounded to an
integer. This is illustrated by the solid lines in figure 3. As an example a
$10^{13}$ eV EAS footprint with a regular pattern is shown in figure 4. Within
each ring the particle density is constant. Each dot represents a single
particle.

If a detector with a size of 1 ${\mathrm{m}}^{2}$ is randomly projected at the
density map, the probability to detect $x$ number of particles scales with the
surface area of each ring. For $x=1$ the probability is equal to the surface
area of the outermost ring divided by the surface area of the entire disc.
This is illustrated in figure 5 which displays the probability density to find
a number of particles as a function of EAS energy. Each probability density
represents a single energy PID. When many particles simultaneously traverse
the detector, a large pulse integral is generated. If only one particle
traverses the scintillator, the value of the pulse integral follows a single
Landau distribution.

In order to predict the PIDs more accurately, a Monte Carlo simulation is
constructed that includes multiple EAS zenith angles and accurate EAS
evolution and detector description.

Figure 3: The dashed lines show the lateral density profiles for $10^{13}$
(blue), $10^{14}$ (red), $10^{15}$ (green) and $10^{16}$ eV (orange) showers
described by eq. 3.1. The solid lines show the multiplicity rounded to the
nearest integer. In the calculation only EASs from the zenith are considered.
Figure 4: Example of a $10^{13}$ eV perpendicular incident EAS footprint. Each
dot represents a single particle. Within each ring the particle density is
taken constant. When a detector with a size of $1\times 1$ ${\mathrm{m}}^{2}$
is randomly projected on the density map, the probability to detect $x$ number
of particles scales with the surface area of each ring. For $x=1$ the
probability is equal to the surface area of the outermost ring divided by the
surface area of the entire disc. Figure 5: The probability to detect a number
of particles in the footprint at a fixed energy. For low energy EASs the
probability to detect a single MIP is relatively large compared to higher
multiplicities. With increasing energy the slope of the distributions changes
significantly.

###  3.2 Detector separation

The small particle densities at low energy EASs ($\lesssim 10^{13}$ eV)
generate a steeply falling PID, i.e. there is a relatively high probability to
detect just 1 or 2 particles; higher multiplicities at this energy become
rare. Increasing the distance between the detectors then translates into an
energy cut. Since the size and shape of the PID directly depends on the number
of particles that ’simultaneously’ traverse the detector, it provides an
indirect measure for the energy of the primary cosmic ray. Figure 3 shows that
for $10^{13}$ eV showers (blue curve) the probability to detect particles
beyond 20 m from the shower core becomes extremely small. Thus, if two HiSPARC
detectors are placed 40 m apart, the EAS detection probability at this energy
vanishes. Since it is impossible to select single energy EASs, a PID obtained
by a station can be thought of as a sum of single energy PIDs. Introducing an
energy cut by increasing the detector separation will affect the shape of the
station’s PID as single low energy PIDs are excluded from the collection.

Figure 6 shows the measured distribution of arrival time differences for
particles in the two detectors. The times at which the PMT pulses exceed the
threshold are taken as arrival times. The plateau of random coincidences
(green striped region) extends to small time differences (blue horizontal
line). The peak is caused by particles belonging to EASs. For time differences
smaller than 300 ns the random coincidences are obviously indistinguishable
from EASs. However, the PID distribution obtained by selecting random
coincidences can be subtracted (proportional to the number of background
events - blue crossed region) from the PID for time differences smaller than
300 ns.

Figure 6: Distribution of arrival time differences between particles that
satisfy the trigger condition in a two-detector HiSPARC station (red). The
deviation of the peak value from zero indicates a small timing offset between
the two detectors. The plateau of events (green striped region) is
extrapolated towards small time differences (blue horizontal line) and is the
result of random coincidences (PMT noise, etc.). For time differences smaller
than 300 ns (blue crossed region) the random coincidences are
indistinguishable from air shower events. Their contribution is however small.

A small experiment has been carried out in which the PID is measured for
different detector separations (5, 15 and 35 m). The result is shown in figure
7. The number of single MIP events decreases with increasing detector
separation. The number of large pulse integrals stays approximately constant
since these events are caused by high energy showers.

Figure 7: PID for different detector separations. By increasing the distance
between the detectors, the station will not trigger on low energy EAS
footprints due to their small particle density (Fig. 5). The number of single
MIP detections decreases proportionally to the increasing detector separation.
The number of large pulse integrals (multiple simultaneous particles)
traversing the detector stays approximately constant since these events are
caused by high energy showers which are still detected at large detector
separations.

###  3.3 Systematic uncertainties

The relation between the shape of the PID and the energy of the primary cosmic
ray can then be used to probe the cosmic ray flux at a certain energy range.
To quantify systematic differences between detectors, the PIDs of detector
combinations in a four-detector station (diamond shaped configuration) [4]
have been compared. Each detector in this station can be paired with three
other detectors. The four sides of the diamond are 10 m. The long diagonal
combination is discarded. This results in ten PID distributions. Figure 8
shows the ten distributions of the four detectors. All combinations with
’detector one’ are displayed in blue, all combinations with detector two are
displayed in red, etc. The PIDs obtained with the same detector are very
similar. The difference between detectors is larger and becomes more visible
at large pulse integrals (see orange and green lines). This indicates that
there are instrumental differences caused by variations in the number of
scintillation photons reaching the PMT, large pulse PMT response, etc. After
closer inspection, systematic differences between the PIDs up to $\sim$35 Vns
($\sim$12 MIPs) are small.

Figure 8: PIDs of ten detector combinations in a four-detector (diamond shaped
configuration) station. Each color represents one detector. Per detector two
or three PIDs are displayed. The detector separation for each combination is
10 m. Combinations with the same detector (e.g. red lines) are very similar.
The difference between the detectors is larger, especially at large pulse
integrals (see orange and green lines).

##  4 EAS Simulations

A large sample of EASs for primary cosmic rays with an energy of $10^{12}$,
$10^{13}$, $10^{14}$, $10^{15}$ and $10^{16}$ eV is generated with CORSIKA.
Only proton initiated showers are considered. The showers have been generated
with zenith angles ranging from 0$\mathrm{\SIUnitSymbolDegree}$ to
60$\mathrm{\SIUnitSymbolDegree}$ in steps of
3.75$\mathrm{\SIUnitSymbolDegree}$. For high energy hadronic interactions
QGSJET-II [10] was selected. Hadronic interactions below 80 GeV are simulated
using GHEISHA [11] and electromagnetic interactions are incorporated with EGS4
[12]. The full particle shower was simulated (no ’thinning’ [13]). The
location of the EAS shower cores were randomly chosen within a circle with a
radius of 150 m. The two-detector station is at the center of the circle. The
arrival direction (azimuth and zenith) was chosen isotropically (random points
on the surface of a unit sphere). The response of the scintillator and light-
guide to particles traversing the detector was simulated with GEANT4 [5]. A
parameterized PMT response was used [4]. Finally, the HiSPARC trigger
conditions were applied.

Figure 9 shows the normalized, simulated PIDs for the five different energies.
At higher energies the relative abundance of large pulses increases
considerably. By fitting a linear combination of the simulated distributions
to the experimentally observed spectrum, the cosmic ray flux at a fixed energy
interval can be estimated.

Figure 9: Simulated PIDs at five different energies. As expected from figure 5
the PID depends on the energy of the EASs. Contrary to the model, a full
detector simulation is carried out with isotropically selected EASs and
uniformly chosen core distances. Only proton induced EASs are considered. The
difference becomes evident at large pulse integrals.

###  4.1 Effective surface area and solid angle

EASs at large zenith angles result in lower particle densities at ground level
because of increased absorption due to the larger path length through the
atmosphere. Moreover, the particle density observed in the detection station
strongly depends on the distance to the shower core. These two effects need to
be quantified in order to calculate the cosmic ray flux. For this an effective
surface area and effective solid angle are introduced. The fraction of events
that results in a trigger is defined as the ’EAS detection efficiency’
($\varepsilon$). The left plot in figure 10 shows $\varepsilon$ for $10^{15}$
eV air showers as function of the zenith angle and distance to the shower
core. The maximum efficiency occurs at small $r$ and small zenith angle
$\theta$.

Integrating over the solid angle yields the effective solid angle as a
function of core distance (upper right plot):

$\Omega(r)=\int_{0}^{2\pi}\int_{0}^{\pi/2}\varepsilon(\theta,r)\sin\theta
d\theta d\phi$ (4)

The effective solid angle can accurately be parameterized using the
exponentially modified Gaussian distribution:

$\displaystyle\Omega(r;\alpha,\mu,\sigma,\lambda)=$
$\displaystyle\quad\alpha\exp\left[\frac{\lambda}{2}(2\mu+\lambda\sigma^{2}-2r)\right]$
$\displaystyle\times\mathrm{erfc}\left(\frac{\mu+\lambda\sigma^{2}-r}{\sqrt{2}\sigma}\right)$
(5)

with $\alpha$ a scaling parameter, $\mu$ and $\sigma$ the mean and standard
deviation of the Gaussian part of the distribution and $\lambda$ the rate of
the exponential part. The complementary error function, $\mathrm{erfc}(x)$, is
given by:

$\textup{erfc}(x)=\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-y^{2}}dy$ (6)

Integrating the detection efficiency over the surface area (polar coordinates
$r$ and $\zeta$) yields the effective surface area as a function of zenith
angle (bottom right plot in figure 10):

$\textup{A}(\theta)=\int_{0}^{2\pi}\int_{0}^{\infty}\varepsilon(\theta,r)rdrd\zeta$
(7)

The effective surface area can be parameterized using the following formula:

$\textup{A}(\theta)=a\exp\left(-b\left(\frac{1}{\cos\theta}-1\right)\right)$
(8)

with $a$ and $b$ fit parameters. Since, in order to obtain the flux, the
number of events needs to be divided by both the effective solid angle and the
effective surface area; there is no need to evaluate them separately. Instead,
the two are combined:

$\displaystyle\textup{A}\Omega$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\infty}\Omega(r)rdrd\zeta$ (9)
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{\pi/2}A(\theta)\sin\theta d\theta
d\phi$ (10)

Table 1 shows the fit parameters that describe the effective solid angle and
surface area for the simulated energies (for $10^{15}$ eV proton showers
$\textup{A}\Omega$ is $8.58\cdot 10^{3}$
${\mathrm{m}}^{2}\text{\,}\mathrm{sr}$). The $\textup{A}\Omega$ values can be
calculated from eqs. 9 or 10, or by direct summation of the two-dimensional
histogram in figure 10. The last column in table 1 lists the averages and
standard deviations of these three methods.

Figure 11 shows the five $\textup{A}\Omega$ values. The interpolation is
defined by the following equation:

$\log_{10}(\textup{A}\Omega)=ax^{2}+bx+c$ (11)

with $x=\log_{10}(E)$, $a=-0.239$, $b=8.31$, $c=-66.8$. Extrapolation leads to
an estimate at $10^{11}$ eV of $\textup{A}\Omega=4.2\cdot 10^{-5}$ (red
circle). A precise direct estimate at this energy from simulation is difficult
because it is computationally expensive to collect a sufficiently large data
set; the number of EASs that satisfy the selection criteria is extremely
small. The number of triggers at an EAS energy of $10^{11}$ eV is
approximately 400 times smaller than at $10^{12}$ eV. This is consistent with
the estimate obtained by extrapolation of eq. 11.

Figure 10: Left: The detection efficiency ($\varepsilon$) for $10^{15}$ eV proton induced showers as function of core distance ($r$) and zenith angle ($\theta$). Top right: the detection efficiency integrated over the zenith angle yields the effective solid angle as function of core distance ($\Omega(r)$, blue points) which can accurately be described using the parametrisation in eq. 4.1 (black line). Bottom right: the detection efficiency integrated over the core distance yields the effective surface area as function of zenith angle ($\textup{A}(\theta)$, blue points) can be described by the parametrisation in eq. 8 (black line). Table 1: The parameters in eq. 4.1 (effective solid angle $\Omega(r)$) and eq. 8 (effective surface area $\textup{A}(\theta)$) are listed as a function of energy. The last column gives the value of the combination $\textup{A}\Omega$ (eq. 9 or 10). Energy [eV] | $\Omega(r)$ |  |  |  | $\textup{A}(\theta)$ |  |  $\textup{A}\Omega$ [${\mathrm{m}}^{2}\text{\,}\mathrm{sr}$]
---|---|---|---|---|---|---|---
| $\alpha$ | $\mu$ | $\sigma$ | $\lambda$ | $a$ | $b$ |
$10^{12}$ | $1.73\cdot 10^{-5}$ | $-2.3$ | 1.0 | $8.00\cdot 10^{-2}$ | $4.86\cdot 10^{-2}$ | 9.66 | $(2.75\pm 0.09)\cdot 10^{-2}$
$10^{13}$ | $8.71\cdot 10^{-3}$ | $-3.1$ | 3.6 | $1.22\cdot 10^{-1}$ | 7.85 | 6.89 | $5.50\pm 0.14$
$10^{14}$ | $5.44\cdot 10^{-1}$ | $-10.4$ | 8.1 | $8.78\cdot 10^{-2}$ | $6.30\cdot 10^{2}$ | 7.90 | $(4.23\pm 0.21)\cdot 10^{2}$
$10^{15}$ | 4.40 | $-14.2$ | 24.5 | $7.25\cdot 10^{-2}$ | $9.21\cdot 10^{3}$ | 4.76 | $(8.58\pm 0.28)\cdot 10^{3}$
$10^{16}$ | 6.92 | $-50.7$ | 63.3 | $2.71\cdot 10^{-2}$ | $6.33\cdot 10^{4}$ | 3.62 | $(7.33\pm 0.13)\cdot 10^{4}$
Figure 11: The blue curve (eq. 11) shows the interpolation of the
$\textup{A}\Omega$ values (blue dots) listed in in table 1. $\textup{A}\Omega$
at $10^{11}$ eV (red circle, $\textup{A}\Omega=4.2\cdot 10^{-5}$
${\mathrm{m}}^{2}\text{\,}\mathrm{sr}$) is obtained by extrapolation.

##  5 Fitting single energy PIDs to experimental data

The data set used for figure 6 was used to obtain a PID. The contribution to
the PID from random coincidences was subtracted (see section 3.2). A Bayesian
method was used to fit a combination of simulated single energy PIDs in the
range $10^{12}$$-$$10^{16}$ eV to the experimental data. Only pulse integrals
smaller than 35 Vns are considered to limit systematic differences between the
detectors (figure 8). Signals are affected by gamma rays, and Cherenkov light
generated in the light-guide [4]. Especially at lower multiplicities this
contribution becomes apparent. Pulse integrals below 6 Vns (2 MIPs) are
therefore discarded as well.

As shown in figure 11, the probability that a $10^{11}$ eV shower results in a
footprint that triggers the station is negligible. The slope of PIDs with
energies $10^{17}$ eV and beyond becomes rather similar to that of $10^{16}$
eV in the pulse integral range between 6 to 35 Vns whereas the flux at those
energies rapidly decreases. The experimentally observed PID will therefore be
restricted to the sum of five single energy contributions
($10^{12}$$-$$10^{16}$ eV):

$\displaystyle z(n_{i},\vec{x})=$ $\displaystyle\
x_{1}\cdot\textup{PID}_{12}(n_{i})+x_{2}\cdot\textup{PID}_{13}(n_{i})$
$\displaystyle+x_{3}\cdot\textup{PID}_{14}(n_{i})+x_{4}\cdot\textup{PID}_{15}(n_{i})$
$\displaystyle+x_{5}\cdot\textup{PID}_{16}(n_{i})$ (12)

Here $\vec{x}\equiv(x_{1},x_{2},x_{3},x_{4},x_{5})$, $n_{i}$ is the bin value,
$z(n_{i})$ is the expected number of events in each bin and each
$\textup{PID}(n_{i})$ is a single energy model (fig. 9). The uncertainty in
the number of entries in the bins of the experimentally observed PID is
described by a Poisson distribution. The probability to obtain a number of
counts $y_{i}$ in bin $n_{i}$ given $\vec{x}$ is given by:

$p(y_{i},n_{i}|\vec{x})=P(y_{i},z(n_{i},\vec{x}))$ (13)

with the Poisson distribution:

$P(y,z)=\frac{z^{y}}{y!}e^{-z}$ (14)

Bayes’ rule can be exploited to define a probability to obtain the fit
parameters $\vec{x}$.

$p(\vec{x}|N,Y)=\frac{p(N,Y|\vec{x})p(\vec{x})}{p(N,Y)}$ (15)

Here $N$ and $Y$ are the combined collection of elements $n_{i}$ and $y_{i}$.
The denominator $p(N,Y)$ is a normalization constant. The $p(\vec{x})$ is
known as the prior probability distribution. The prior states that the
parameters $x_{1}$, $x_{2}$, $x_{3}$, $x_{4}$ and $x_{5}$ cannot become
negative. The function $p(\vec{x}|N,Y)$ is the posterior probability
distribution. The $p(N,Y|\vec{x})$ function is the likelihood (the product of
conditional probabilities):

$p(N,Y|\vec{x})=\prod_{i=1}^{M}p(y_{i},n_{i}|\vec{x})=\mathcal{L}$ (16)

Since the slope differences between the PIDs occur at higher multiplicities, a
gradually larger weight is assigned with increasing multiplicity. This leads
to a weighted likelihood:

$\mathcal{\widehat{L}}=\prod_{i=1}^{M}p(y_{i},n_{i}|\vec{x})^{w(n_{i})}$ (17)

with $w(n_{i})$ the weight function:

$w(n_{i})=1.161^{n_{i}}$ (18)

which corrects for the number of entries in the lower ($1.7\cdot 10^{4}$ at 6
Vns) and higher bins (231 at 35 Vns). By taking the logarithm of this weighted
likelihood, the product converts into:

$\displaystyle\log(\mathcal{\widehat{L}})$
$\displaystyle=\sum_{i=1}^{M}w(n_{i})\cdot\log\left(\frac{{z_{i}}^{y_{i}}}{y_{i}!}e^{-z_{i}}\right)$
(19)
$\displaystyle=\sum_{i=1}^{M}w(n_{i})\cdot(y_{i}\log(z_{i})-z_{i}-\log(y_{i}!))$
(20) $\displaystyle=\sum_{i=1}^{M}w(n_{i})\cdot(y_{i}\log(z_{i})-z_{i}-C)$
(21)

with $z_{i}=z(n_{i},\vec{x})$ and $C$ is a constant.

Instead of directly maximizing the posterior probability distribution, a range
of parameters ($\vec{x}$) are explored. This is done using a Markov Chain
Monte Carlo (MCMC) algorithm [14]. Figure 12 shows (lower dimensional) subsets
of the sampled (five dimensional) posterior probability distribution. The
upper (on the diagonal) Gaussian shaped histograms show the sample selection
in the dimension of the fit parameters. The median of the histogram (red
lines) is taken as the best fit. The uncertainty is shown by the standard
deviation (green lines). The best fit values are $(3.01\pm 0.35)\cdot 10^{4}$,
$(2.32\pm 0.05)\cdot 10^{5}$, $(1.09\pm 0.03)\cdot 10^{5}$, $(5.63\pm
0.20)\cdot 10^{4}$, $(2.63\pm 0.16)\cdot 10^{4}$ for $10^{12}$ to $10^{16}$ eV
resp. The other subplots display the relation between two fit parameters.
There is some interdependence between the fit parameters of neighboring
energies. This is especially pronounced at fit parameters $x_{4}$ and $x_{5}$
(fourth sub-plot in bottom row) due to the relatively small slope differences
between the $10^{15}$ and $10^{16}$ eV PIDs (figure 9). The resulting fit to
the experimental data is shown in figure 13.

In addition to the statistical uncertainty in the experimental data, the
single energy models also have an intrinsic statistical uncertainty. Creating
larger data sets for the single energy PIDs is currently limited by the
generation of simulated EASs which is computationally expensive. The
uncertainty has been estimated by resampling the single energy PIDs. The
number of events in each bin was randomly resampled following Poisson
statistics. A combination of these new PIDs was fitted to the experimental
data as well. This procedure was carried out multiple times. The mean and
standard deviation of the best fit parameters are listed in table 2.

Table 2: Best fit for the number of events with their uncertainties. Energy [eV] | Number of events
---|---
$10^{12}$ | $(4.09\pm 1.76)\cdot 10^{4}$
$10^{13}$ | $(2.04\pm 0.23)\cdot 10^{5}$
$10^{14}$ | $(1.25\pm 0.20)\cdot 10^{5}$
$10^{15}$ | $(5.65\pm 1.32)\cdot 10^{4}$
$10^{16}$ | $(2.43\pm 1.04)\cdot 10^{4}$
Figure 12: Subsets of the sampled posterior probability distributions. The
upper Gaussian shaped histograms show the sample selection in the dimension of
the fit parameters ($\vec{x}$). The median of the histogram (red lines) is
taken as the best fit. The uncertainty in the fit is shown by the standard
deviation (green lines). The other subplots display the relation between two
fit parameters. There is some interdependence between the fit parameters of
neighboring energies. This is more pronounced for $x_{4}$ and $x_{5}$ (fourth
sub-plot in bottom row). See also figure 9. Figure 13: A linear combination of
five simulated PIDs (colored lines, see figure 9) fitted to the experimentally
obtained PID (thick grey line). The black line shows the best fit. The bottom
plot gives the ratio of the experimental data and the fit.

##  6 Cosmic ray flux

The estimated number of events per energy decade derived from the fit can be
used in combination with the effective surface area and solid angle to obtain
the cosmic ray flux as a function of energy. The fluxes are calculated using:

$F=\frac{N_{\textup{events}}}{\textup{A}\Omega\cdot t\cdot\Delta E}$ (22)

with $t$ the duration of the experiment ($t=29$ days $=2.506\cdot 10^{6}$ s)
and $\Delta E$ the width of the energy bin (e.g. $10^{15.5}-10^{14.5}\approx
2.85\cdot 10^{15}$ eV). The flux values are listed in table 3. Figure 14 shows
the cosmic ray energy spectrum (circles) from several experiments [15, 16, and
references therein] together with the flux values derived from HiSPARC data
(red dots). The uncertainties are smaller than the dot size. The grey line
represents the function in equation 1. The HiSPARC values agree well with the
other measurements. Note that the reference flux values below $10^{14}$ eV are
measured using spacecraft (e.g. Proton satellite [17]).

Table 3: Cosmic ray flux values obtained using a two-detector HiSPARC station. Energy [eV] | Flux [${\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{sr}}^{-1}\text{\,}{\mathrm{GeV}}^{-1}$]
---|---
$10^{12}$ | $(2.08\pm 0.93)\cdot 10^{-4}$
$10^{13}$ | $(5.19\pm 0.60)\cdot 10^{-7}$
$10^{14}$ | $(4.14\pm 0.70)\cdot 10^{-10}$
$10^{15}$ | $(9.22\pm 2.20)\cdot 10^{-13}$
$10^{16}$ | $(4.64\pm 2.01)\cdot 10^{-15}$
Figure 14: Cosmic ray fluxes at $10^{12}$$-$$10^{16}$ eV from HiSPARC data
(red dots) compared to results from other experiments [15, 16, and references
therein] (circles). The uncertainties are smaller than the dot size. The
fluxes obtained by HiSPARC agree well with the other measurements. The grey
line represents eq. 1.

##  7 Discussion and conclusion

A HiSPARC station (at sea-level) with only two 0.5 ${\mathrm{m}}^{2}$
scintillator detectors was used to probe the cosmic ray energy spectrum in the
range $10^{12}$$-$$10^{16}$ eV (five energy decades). One month of data shows
to be sufficient to obtain flux values that are in good agreement with the
results from dedicated (space-based) experiments.

The presented analysis is very different from those applied by other ground-
based cosmic ray observatories. The method does not require the reconstruction
of the shower core and/or shower size of individual EASs but relies on
deriving the energy dependent particle multiplicity distributions in a single
scintillator detector. PMT pulse integrals are used as a proxy to estimate the
number of particles simultaneously traversing the detector. In comparing
simulated PIDs at $10^{15}$ and $10^{16}$ eV (and beyond $10^{16}$ eV) with
experimental data, the simulation starts to suffer from lack of statistics;
showers become very large exceeding the presently available CPU power and data
storage facilities. Re-doing the analysis with simulations in which ’thinning’
is applied has not been explored but may offer an alternative. On the other
hand, as becomes apparent from figure 9, the difference in slope of the PIDs
towards higher particle multiplicities decreases with increasing energy.
Moreover, at these large multiplicities, the dynamical range of the readout
electronics and non-linearities in detector response require careful
evaluation. By extending the distance between the detectors, the contribution
from small, low energy ($10^{12}$$-$$10^{13}$ eV) EASs can be reduced,
improving the sensitivity of the analysis towards higher energy EASs.

Since the HiSPARC network consists of more than 120 stations in which the
separation between the scintillators varies between 5 and 17 m while for
several stations more than 10 years of data are stored, statistics can easily
be increased and systematic uncertainties can be investigated in detail.

Finally, the cosmic ray flux below $10^{12}$ eV could possibly be derived by
increasing the altitude of the HiSPARC station, compensating for the shift of
the shower maximum towards higher altitudes. The two-detector station in
Windhoek (Namibia) is at an altitude of almost 2 km…

## Acknowledgements

We thank Hogeschool lecturer and high school teacher Kees van der Velden for
carrying out the measurements at the roof of his high school presented in
figure 7. Also the support of the contributing HiSPARC schools, their teachers
and especially their students is gratefully acknowledged. Without their help,
it would have been impossible to keep the stations ’up and running’.

## References

  * [1] T. K. Gaisser, R. Engel, E. Resconi, Cosmic rays and particle physics, Cambridge University Press, 2016.
  * [2] J. Beringer, J. Arguin, R. Barnett, K. Copic, O. Dahl, D. Groom, C. Lin, J. Lys, H. Murayama, C. Wohl, et al., Review of particle physics, Physical Review D-Particles, Fields, Gravitation and Cosmology 86 (1) (2012) 010001.
  * [3] D. Heck, J. Knapp, J. N. Capdevielle, et al., CORSIKA: a Monte Carlo code to simulate extensive air showers., Forschungszentrum Karlsruhe GmbH, Karlsruhe (Germany), 1998.
  * [4] K. van Dam, B. van Eijk, D. Fokkema, et al., The HiSPARC experiment, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 959 (2020) 163577.
  * [5] J. Allison, K. Amako, J. Apostolakis, et al., Recent developments in Geant4, Nucl. Instrum. Methods Phys. Res. A 835 (2016) 186–225.
  * [6] K. van Dam, Increased radioactivity during precipitation measured by the HiSPARC experiment, Physica Scripta 95 (2020) 074011.
  * [7] L. Landau, On the energy loss of fast particles by ionization, J. Phys.(USSR) 8 (1944) 201–205.
  * [8] M. Takeda, N. Sakaki, K. Honda, M. Chikawa, M. Fukushima, N. Hayashida, N. Inoue, K. Kadota, F. Kakimoto, K. Kamata, et al., Energy determination in the akeno giant air shower array experiment, Astroparticle Physics 19 (2003) 447–462.
  * [9] H. Dai, K. Kasahara, Y. Matsubara, M. Nagano, M. Teshima, On the energy estimation of ultra-high-energy cosmic rays observed with the surface detector array, Journal of Physics G: Nuclear Physics 14 (1988) 793.
  * [10] S. Ostapchenko, QGSJET-II: towards reliable description of very high energy hadronic interactions, in: Nuclear Physics B - Proceedings Supplements, Vol. 151, 2006, pp. 143–146.
  * [11] H. Fesefeldt, GHEISHA The Simulation of Hadronic Showers, Tech. rep., RWTH Aachen, PITHA-85/02 (1985).
  * [12] W. R. Nelson, D. W. O. Rogers, H. Hirayama, The EGS4 code system, Tech. rep., Stanford Linear Accelerator Center, Stanford, California (1985).
  * [13] A. M. Hillas, , in: Proceedings International Cosmic Ray Conference, Vol. 1, 1981, p. 193.
  * [14] D. Foreman-Mackey, D. W. Hogg, D. Lang, J. Goodman, emcee: the mcmc hammer, Publications of the Astronomical Society of the Pacific 125 (2013) 306.
  * [15] Cosmic Ray Spectrum website by Dr. William Hanlon, https://www.physics.utah.edu/~whanlon/spectrum.html.
  * [16] J. W. Cronin, T. K. Gaisser, S. P. Swordy, Cosmic rays at the energy frontier, Scientific American 276 (1997) 44–49.
  * [17] N. Grigorov, Y. V. Gubin, I. Rapoport, I. Savenko, B. Yakovlev, V. Akimov, V. Nesterov, Energy spectrum of primary cosmic rays in the 10 11-10 15 ev energy range according to the data of proton-4 measurements, in: 12th International conference on cosmic rays, 16-25th August 1971, Hobart, Tasmania, Australia, 1971.