# Sibyll★

Felix Riehn  Anatoli Fedynitch  Ralph Engel  Instituto Galego de Física de
Altas Enerxías (IGFAE), Universidad de Santiago de Compostela, 15782 Santiago
de Compostela, Spain  Institute of Physics, Academia Sinica, Taipei City,
11529, Taiwan  Karlsruhe Institute of Technology, Institute for Astroparticle
Physics, 76021 Karlsruhe, Germany  Karlsruhe Institute of Technology,
Institute of Experimental Particle Physics, 76021 Karlsruhe, Germany

###### Abstract

In the last decade, an increasing number of datasets have revealed a
consistent discrepancy between the number of muons measured in ultra-high-
energy extensive air showers (EAS) and the numbers predicted by simulations.
This gap persists despite incorporating Large Hadron Collider (LHC) data into
the tuning of current hadronic interaction models, leading to the phenomenon
often termed the “muon puzzle”. To gain a deeper understanding of the
potential origins of this muon puzzle, we have developed Sibyll★, a series of
phenomenologically modified versions of Sibyll 2.3d. In these models, we have
increased muon production by altering $\rho^{0}$, baryon–antibaryon pair, or
kaon production in hadronic multiparticle production processes. These variants
remain within bounds from provided by accelerator measurements, including
those from the LHC and fixed-target experiments, notably NA49 and NA61,
showing a level of consistency comparable to Sibyll 2.3d. Our findings show
that these modifications can increase the muon count in EAS by up to 35%,
while minimally affecting the depth of shower maximum ($X_{\rm max}$) and
other shower variables. Additionally, we assess the impact of these
modifications on various observables, including inclusive muon and neutrino
fluxes and the multiplicities of muon bundles in deep underground and
water/ice Cherenkov detectors. We aim for at least one of these model variants
to offer a more accurate representation of EAS data at the highest energies,
thereby enhancing the quality of Monte Carlo predictions used in training
neural networks. This improvement is crucial for achieving more reliable data
analyses and interpretations.

###### keywords:

cosmic rays , extensive air showers , muon production , hadronic interactions
, muon puzzle , atmospheric leptons

††journal: Astroparticle Physics

##  1 Introduction

The phenomenon known as the “muon puzzle” in extensive air showers (EAS) –
referring to a significant discrepancy between the number of muons observed in
EAS with respect to those predicted by standard simulation models – has led to
a critical reassessment of our understanding of hadronic interactions at
ultra-high energies [1, 2]. To address this issue, both conventional [3, 4, 5,
6, 7] and exotic [8, 9, 10, 11, 12, 13] modifications of multiparticle
production, some including extensions beyond the Standard Model of particle
physics, have been proposed.

While the muon puzzle is of great scientific interest for learning more about
hadronic multiparticle production, it also gives rise to large systematic
uncertainties in the analysis and interpretation of EAS data. This is
especially critical for training machine learning algorithms [14, 15, 16, 17]
to reconstruct, for example, the depth of shower maximum $X_{\rm max}$ and the
mass composition of cosmic rays from the ground-level signals of EAS surface
arrays. Having a better understanding of the origin of the muon puzzle and
being able to simulate air showers with a more realistic muon number are
highly desirable for reducing the systematic uncertainties and, in some cases,
for the first time, obtaining consistent measurements using different
observables.

The typically highly indirect way of deriving information on muon numbers or
densities with EAS detectors that are not built for measuring muons separately
on the first hand, like the Pierre Auger Observatory [18, 19], Telescope Array
[20] or IceTop [21], makes a straightforward experimental test of the various
proposals to address the muon puzzle very difficult. Detailed simulations of
realistic air showers incorporating these proposals, covering different
primary particles and a wide range in energy, would have to be compared with
existing data sets. However, most of the proposed modifications have either
not been worked out to such a level of detail or have not been implemented in
a hadronic model capable of generating realistic air showers.

In this work, we introduce Sibyll★ (in the following also abbreviated as S★),
a set of phenomenological modifications to the Sibyll 2.3d [22] hadronic
interaction model. These modifications explore various conventional mechanisms
with the objective of increasing muon production in EAS. The modified versions
of Sibyll can be directly utilized in realistic EAS simulations, ensuring that
energy, momentum, and quantum numbers are conserved on an event-by-event
basis. One important advantage of these modified versions is that,
intentionally, a minimum number of changes has been applied to keep the
overall model predictions as much as possible identical to that of Sibyll
2.3d. This will allow a direct comparison with already existing simulations
and data analyses/interpretations made with Sibyll 2.3d. Furthermore, making a
direct connection between a change of EAS predictions and the underlying
modifications to the hadronic interaction model will be possible.

This article is organized as follows. In Sec. 2 the algorithms of making
modifications to the particle output of Sibyll 2.3d are introduced and the
tuning of the relevant parameters is discussed. It is shown that the
predictions of the modified variants of the interaction model are compatible
with a selected set of data of accelerator experiments that have direct
sensitivity to the modifications. In many cases the agreement of the model
predictions with the data even improves. In total, four variants of the Sibyll
interaction model are presented. These model variants are then used in Secs. 3
and 4 to study the model predictions for air showers and inclusive fluxes,
respectively, and to compare them with Sibyll 2.3d as reference. It is argued
that, at the level of comparison possible without including a detector
simulation, a model variant like, for example, Sibyll★(mix) should properly
reproduce the muon data of the Pierre Auger Observatory. The impact of the
modifications on the predicted inclusive muon and neutrino fluxes is discussed
with emphasis of possible observations with IceCube. Finally, a summary of the
results is given in Sec. 5.

##  2 Sibyll★

###  2.1 Ad-hoc Event Modification

We construct the custom models by modifying events generated with Sibyll 2.3d.
Once the initial event generation is complete, we let all hadronic resonances
with lifetimes shorter than that of K${}^{0}_{\rm s}$ decay, except for
$\pi^{0}$. In the next step we replace particles according to the selection
probabilities described below.

To replace a particle, for example a $\pi^{0}$ by a $\rho^{0}$, in a given
generated event without violating energy-momentum conservation one has to
change the momenta of (at least) two particles that together reach a
sufficiently large invariant mass. Therefore, we iterate through the list of
generated particles (event stack) to identify appropriate pairs or triples of
pions, while for each pion considering only the five nearest neighboring pions
in rapidity. If the sampling criterion for replacing a particle with another
one, or replacing both particles, is fulfilled and there is sufficient
invariant mass, we replace these pions with a pair of new particles conserving
total momentum, invariant mass, and charge. The final momenta are calculated
from the invariant mass, the masses of the new particles, and a small
transverse momentum component drawn from an exponential distribution in
transverse mass. While energy, momentum, and charge are conserved by
construction, the (iso)spin conservation is not explicitly maintained111by
favoring $\rho^{0}$ over the charged $\rho$ states it is in fact explicitly
broken.

The acceptance rate for particle exchanges in our model depends on both the
total center-of-mass (c.m.) energy, $\sqrt{s}$, of the interaction and the
fractional longitudinal momentum of the considered secondary particle $x_{\rm
F}\equiv p_{z}/p_{z,\mathrm{max}}$ (with momenta $p$ in the c.m. frame). The
exchange probability is parameterized by

$P_{i}~{}=~{}P_{i,0}\cdot|x_{\rm F}|^{\epsilon_{i}}\cdot f(\sqrt{s},E_{\rm
thr})\ .$ (1)

Whether the modifications are applied mainly to forward or central particle
production depends on the chosen value for the exponent $\epsilon_{i}$ in the
$x_{\rm F}$-dependence. If $\epsilon_{i}=0$, all particles receive equal
weight, effectively preserving the shape of the original distribution in
longitudinal phase space. In the case of $\epsilon_{i}=1$, the modification is
proportional to $x_{\rm F}$ and the forward part of the $x_{\rm F}$-spectrum
undergoes significant enhancement (with maximal enhancement at $x_{\rm F}\to
1$). And so on for larger values of $\epsilon_{i}$.

The energy dependence of $f(\sqrt{s},E_{\rm thr})$ follows a logarithmic form.
It is specifically chosen to ensure that the modification probability is
precisely zero below a set threshold energy, $E_{\rm thr}$, and reaches its
nominal value of $P_{0,i}$ at lab energies of $10^{19}\,$eV ($1.37\times
10^{5}\,$GeV in the c.m. frame). This parameterization of the energy
dependence ensures a gradual change in particle production. At low energies,
where fixed-target experiments effectively constrain particle production
across the entire phase space, no event modification is permitted. At LHC
energies, only the central region is well-constrained by measurements, leaving
room for modifications in particle production for the forward phase space. The
constraints from the LHCf experiment [23, 24, 25, 26] only apply to forward
neutral particle production, whereas FASER [27] and SND [28] or future
experiments at the FPF [29] may provide stricter constraints on forward
charged particle production by observing muons and neutrinos [30]. Finally, at
the UHECR energy scale, where laboratory experiments offer no constraints,
particle production undergoes substantial modification.

As an alternative, we implement a more significant increase in the exchange
rate at energies above $13\,$TeV in the c.m. frame, expressed as $P_{i}\to
P_{i,0}+P_{i,\mathrm{HE}}\cdot f_{\rm HE}(\sqrt{s}=13\,\mathrm{TeV})$. This
mode simulates a rapid transition to new physics beyond the c.m. energy of the
LHC.

We employ this algorithm for all interactions in an EAS if one can assume the
enhancement should be universal by its nature. Only in the case of $\rho^{0}$
enhancement, directly related to the leading particle effect of pions, we
restrict it solely to pion-air interactions.

###  2.2 Tuning Enhancement Parameters with Accelerator Measurements

Figure 1: Upper panel: fraction of projectile energy carried by $\rho^{0}$
(left), anti protons (center) and charged kaons (right) in $\pi^{-}$C
collisions [31, 32]. The values of the energy fraction in kaons were obtained
by integrating the published $x_{\rm F}$-spectra. Lower panels: multiplicities
of $\pi^{+}$, anti protons and charged kaons in pp collisions [33, 34, 35].
Lines are Sibyll 2.3d, Sibyll 2.1 and different variants of Sibyll★.

Leveraging the algorithm described above, we construct distinct variants of
Sibyll 2.3d with the explicit goal of increasing muon production in EAS.
Within an EAS, the bulk of muons arise from meson decay, which in turn are
predominantly generated through the interactions of other hadrons (the
hadronic cascade). The total number of muons in a shower, therefore, depends
critically on the energy retained within the hadronic cascade at each
interaction step. More precisely, increasing the fraction of energy ($f_{\rm
had}$) carried by hadrons that reinteract in the shower (which, to a very good
approximation, encompasses all hadrons except neutral pions and $\eta$
resonances) leads to a corresponding increase in the total muon count within
the shower. The important point here is that the hadronic energy is changed
across the whole hadronic cascade, from the highest energy interactions down
to the low energy interactions that produce the pions that will then decay
into muons. This way small changes in each interaction accumulate to produce a
large effect in the muons at the ground. The opposite approach, dramatically
changing the highest energy interaction and leaving the low energy
interactions unchanged, also serves to increase the number of muons at the
ground but would affect the fluctuations in the number of muons as well.
These, however, were observed to be described well by current models [36].

We select three specific modifications recognized for their efficacy in
increasing the hadronic energy and enhancing muon production: $\rho^{0}$
production, baryon–antibaryon pair production, and kaon production enhancement
[3, 4, 5, 6, 12]. We denote these variants as S★($\rho^{0}$), S★($\bar{p}$),
and S★(K±,0).

In the $\rho^{0}$ variant, $\pi^{0}$ are directly substituted with $\rho^{0}$.
For the baryon pair and kaon pair variants, charge-neutral combinations of two
or three pions are replaced with $p\bar{p}$ or $n\bar{n}$ pairs, and K+K- or
K0 $\bar{\mathrm{K}}^{0}$ pairs respectively. We adjust the parameters for
each variant to ensure a sufficiently good description of laboratory
measurements. Representative examples include the energy fractions in
$\rho^{0}$, antiprotons, and kaons measured by the NA61 experiment [32, 31]
(upper panels in Fig. 1), as well as the collection of multiplicity
measurements of pions, antiprotons, and kaons in proton-proton interactions
(lower panels in Fig. 1 [33, 34, 35]). Sibyll 2.3d does not always describe
the lab. data well. For example in comparison with the energy fractions in
$\rho^{0}$, $\bar{p}$ and $K^{\pm}$ measured in the NA61 experiment (c.m.
energy around $17\,$GeV), the predictions by Sibyll 2.3d are consistently
lower than the observation (see upper panels in Fig. 1). To allow the improved
description of these data we set the threshold energy $E_{\rm thr}$ in Eq. (1)
to $5\,$GeV in the c.m. frame.

Figure 2: Change in the hadronic energy fraction between Sibyll 2.3d and the
enhanced variants. In the left panel pp interactions are shown. In the panel
on the right $\pi$C interactions are shown.

Within the phase space extending beyond the scope of laboratory measurements,
we set the parameters to achieve an approximately equal hadronic energy
fraction ($f_{\rm had}$) among the different variants at primary energies
around $10^{19}\,$eV. We target a value of $\approx 0.82$, an increase of
$\approx 10\%$ over Sibyll 2.3d (see Fig. 2). $\rho^{0}$ production has a
strong impact on muon production in EAS by redirecting the energy from the
electromagnetic channel into hadrons through the substitution of $\pi^{0}$ by
$\rho^{0}$ in particular for forward $\pi^{0}$. However, the data only support
an enhancement of forward $\rho$ production for meson projectiles [37, 38, 39]
and not for protons. Consequently, we apply the modifications in the
$\rho^{0}$ variant exclusively for the interaction of pion projectiles. In
contrast, the modifications introduced for baryon pair and kaon production
encompass all projectile types. The availability of proton-proton measurements
at significantly higher energies (LHC energy scale) compared to measurements
with pion projectiles constrains the baryon pair and kaon variants to much
higher energies than the $\rho^{0}$ variant. This requires a more rapid
increase of the modifications in these two variants at energies beyond the LHC
to attain the targeted hadronic energy fraction at ultra high energies.

We construct a fourth variant, S★(mix), that combines both $\rho^{0}$ and
baryon pair production enhancements. In this scenario, a more moderate
increase in $\rho$ production and no rapid increase in baryon production at
high energies are needed to achieve similar effects on the showers. Table 1
presents a detailed overview of the parameters employed in the different
variants.

Table 1: Parameters in different variants of Sibyll★. $\rho$-mix and $\bar{p}$-mix refer to the subcomponents of the mixed variant.  Label | $P_{i,0}$ | forward weight $\epsilon_{i}$ | projectiles | $P_{i,\mathrm{HE}}$
---|---|---|---|---
$\rho^{0}$ | 0.9 | 0.3 | mesons | -
$\bar{p}$ | 0.5 | 0.7 | all | 0.25
$K^{\pm,0}$ | 0.5 | 0.8 | all | 0.3
$\rho$-mix | 0.95 | 0.4 | mesons | -
$\bar{p}$-mix | 0.5 | 0.7 | all | -

The above algorithm introduces only gradual changes both to the overall
production rate (and energy) of the particles but also in the shape of the
respective production spectra. The effect on the spectra at $14\,$TeV is
demonstrated in Fig.3 for the case of pp interactions. In the figure the
spectra of the energy fraction in the lab. frame ($x_{\rm L}$) for pions,
kaons and protons are shown. The spectra are weighted by $x_{\rm L}$. Both
baryon and kaon enhancement lead to substantial changes in the production
spectra at large $x_{\rm L}$ in the forward region which leads to the overall
increase in the fraction of energy carried by hadrons (see Fig.2). As the
forward phasespace is outside the coverage of the LHC experiments the observed
production rates (Fig. 1) are not affected.

Figure 3: Spectra of the energy fraction in the lab. frame for pions, kaons
and protons in pp interactions at $14\,$TeV for Sibyll 2.3d and different
variants of S★.

##  3 Predictions for Extensive Air Showers

In the following, we explore the impact of the presented Sibyll variants on
muon production predictions in EAS. First, we review the muon production in
air showers and later discuss the impact on specific observables. All
predictions were calculated using CORSIKA v7.7420 [40] with FLUKA 2021.2.9
[41] as the low-energy interaction model.

###  3.1 Expectation from the Heitler-Matthews Model

As indicated before, muons in EAS are produced by the decay of mesons in the
hadronic cascade. The larger the hadronic cascade (the more mesons) the more
muons will be produced in a shower. The size of the hadronic cascade is
determined by the energy of the primary particle and the fraction of energy
that is emitted into photons in each interaction222A small fraction of mesons
and subsequently muons also arises from the electromagnetic cascade via
photon-nucleus interactions.. The more energy is retained in the hadronic
cascade, the greater the meson abundance and subsequent muon production. For a
proton primary in a shower with energy $E_{0}$, the average number of muons,
$\langle N_{\mu}^{p}(E_{0})\rangle$, is given by

$\langle
N^{p}_{\mu}(E_{0})\rangle=\left(\frac{E_{0}}{\varepsilon}\right)^{\beta}.$ (2)

Here, $\beta$ depends on the energy transfer from the hadronic to the
electromagnetic shower component, and $\varepsilon$ represents the critical
energy – the point where mesons decay rather than interact, effectively
terminating the cascade. In the Heitler-Matthews model [42], assuming
identical particle yields and only pion production in each hadronic
interaction, $\beta$ relates to the ratio of the logarithms of the
multiplicity of neutral pions (decaying into two photons) to all pions. This
simplified case assumes the cascade ends around $130\,$GeV, where pion decay
becomes more likely than interaction. However, real EAS involve various hadron
species with differing critical energies, influencing both $\beta$ and
$\varepsilon$. In simulations with only pions, $\beta$ is around 0.88, while
full Monte Carlo simulations yield values closer to 0.9 [43].

For nuclear primaries, the Heitler-Matthews model utilizes the superposition
principle. This assumes an EAS initiated by a nucleus with $A$ nucleons is
equivalent to $A$ individual proton-initiated showers with energy $E_{0}/A$.
Consequently, the average number of muons for nuclear primaries, $\langle
N_{\mu}^{A}(E_{0})\rangle$, can be written as

$\langle N^{A}_{\mu}(E_{0})\rangle=A\cdot\langle
N^{p}_{\mu}(E_{0}/A)\rangle=A^{1-\beta}\cdot\langle N^{p}_{\mu}(E_{0})\rangle$
(3)

This dependence implies a shrinking separation between different mass
primaries as $\beta$ and the overall muon count increase.

The various Sibyll variants aim to enhance muon production through distinct
mechanisms, primarily by increasing energy retention within the hadronic
cascade (effectively raising $\beta$):

  * 1.

$\rho^{0}$ production: Directly substitutes $\pi^{0}$ with $\rho^{0}$
(decaying into charged pions), altering the charged/neutral pion ratio but not
the critical energy or overall balance of hadron species.

  * 2.

Baryon pair enhancement: Due to baryon number conservation, a baryon cascade
terminates only when the baryons become non-relativistic (around $1\,$GeV).
Creating more baryons within the hadronic cascade lowers the critical energy
and significantly enhances the number of low-energy muons produced in the
final stage of the shower.

  * 3.

Kaon enhancement: Increases production of charged and neutral kaons, leading
to a higher critical energy since kaons have a shorted lifetime than pions
($\varepsilon\sim 1\,$TeV). This increases $\beta$ and muon count, but with
slightly higher-energy muons.

Analyzing the specific predictions of these enhanced models in terms of muon
production and their impacts on EAS observables will be the focus of the next
sections.

###  3.2 Muons at Production

Figure 4:  Left panel: Energy spectrum of muons produced at the shower
maximum. Right panel: Distribution of muon production as a function of depth
(relative to the depth of maximum production). Both figures show the average
over 100 proton showers with a zenith angle of $67^{\circ}$ and a primary
energy of $10\,$EeV.

The zenith angle of the primary, the distance of the observer from the shower
axis (core distance), the observation height, and even atmospheric conditions
(temperature, season) can slightly influence the number of muons detectable at
a specific point on the surface. However, given the energy, transverse
momentum, and depth distribution of muons emerging from the hadronic cascade,
local muon yields can be reliably calculated [44, 45]. Therefore, analyzing
the distributions of muons at their production point is crucial for
understanding the origin of the muon puzzle.

The left panel of Fig. 4 presents the energy spectrum of muons produced around
the shower maximum. The bottom of the figure shows the ratios with respect to
Sibyll 2.3d, revealing the differences between the variants. While $\rho^{0}$
and kaon enhancements increase high-energy muons, baryon enhancement is more
effective at lower energies. All S★ variants exhibit a similar change in the
profile of muon production (shown on the right in Fig. 4). The position of the
production maximum ($X_{\mu}^{\rm max}$) remains largely unaffected (less than
$5\,$g$/$cm2 change, not shown here).

###  3.3 Energy Dependence

Figure 5:  Top Panel: Energy dependence of the surface muon number as a
function of cosmic ray energy. The expectation for nuclei can be obtained
through the superposition model in Eq. (3). The trend mirrors that of the
hadronic energy fraction (Fig. 2). Bottom Panels: Muon energy spectra at three
representative primary cosmic ray energies. Despite the distinct primary
energies, the qualitative features of the muon spectra remain remarkably
consistent. In the high-energy limit $E_{\mu}/E_{p}\rightarrow 1$, dominated
by the first proton-air interaction, the different S★ models produce nearly
identical results.

Figure 5 illustrates the energy dependence of the muon multiplicity, following
the trend indicated by the hadronic energy fraction (compare with Fig. 2). The
flat ratios of the energy spectra in the bottom panels indicate that the
relative increase of the muon number compared to the Sibyll 2.3d reference is
independent of energy thresholds (up to a few tens of GeV). However, at higher
energy thresholds, such as those encountered in underground muon studies or
muon bundle analyses, the S★ variants may begin to demonstrate more
significant differences.

###  3.4 Inclined Showers at the Pierre Auger Observatory

Figure 6: Comparison of the predicted values for $\langle X_{\rm max}\rangle$
and $\langle N_{\mu}\rangle$ from various Sibyll★ variants to the measurements
obtained from the Pierre Auger Observatory [36].

Measurements of the average number of muons in inclined showers at the Pierre
Auger Observatory have unveiled a deficit in the muon count predicted by air
shower simulations [19, 36]. This deficit becomes particularly apparent when
the depth of shower maximum ($\langle X_{\rm max}\rangle$) is also considered.
Since both $\langle X_{\rm max}\rangle$ and $\langle\ln N_{\mu}\rangle$ are
proportional to $\langle\ln{A}\rangle$, the model predictions would fall along
a line in the $\langle X_{\rm max}\rangle$-$\langle\ln N_{\mu}\rangle$ plane
for any given composition scenario. The slope of this line, as dictated by Eq.
(3), is proportional to the exponent $\beta$. The muon deficit in simulations
is revealed by the fact that the Auger data point is well off the predicted
lines for any model. Assuming a composition inferred from $\langle X_{\rm
max}\rangle$ measurements, the deficit for Sibyll 2.3d is around 26% [46].

In Figure 6 the predicted values for $\langle X_{\rm max}\rangle$ and $\langle
N_{\mu}\rangle$ from the S★ variants and the measurements obtained from the
Pierre Auger Observatory are shown. Both the $\rho^{0}$ (red line) and mixed
variants (magenta line) demonstrate a sufficient increase in muon number to
describe the Auger data. The baryon pair (blue line) and kaon enhanced
variants (green line) also lead to a substantial increase in muon count, but
not to the extent required to fully describe the data. It is worth noting that
the slope $\beta$ increases for all the S★ variants. While it is around 0.93
for Sibyll 2.3d, it rises to a value closer to 0.95 for the S★ variants.

Figure 7: $N_{\mu}$ and $X_{\rm max}$ for proton showers at $67^{\circ}$
across Sibyll variants. For $N_{\mu}$ (left panel) there is a substantial
increase in the average, up to 35%, for the mixed and $\rho$ variant. The gray
line represents the required increase in muon count to align with the data for
a mixed composition seen by the Pierre Auger Observatory. The variation in the
average of the shower maximum (panel at the top on the right) between Sibyll
2.3d and its variants is less than $7\,$g$/$cm2. The fluctuations (lower
panels) decrease by up to 20% for $N_{\mu}$ (left) but are not affected for
$X_{\rm max}$ (right).

In Figure 7 the changes in the average (top) and the fluctuations (bottom) in
$N_{\mu}$ (left) and $X_{\rm max}$ (right) across S★ variants are shown for
proton showers at a zenith angle of 67∘. The predicted value of $\langle
X_{\rm max}\rangle$, shown in the right panel on the top, is almost unaltered
due to the changes in muon production. At $10\,$EeV the differences between
all S★ variants and Sibyll 2.3d remain below $7\,$g$/$cm2. In the top left
panel of Fig. 7 the relative change in the predicted number of muons is shown.
It is important to note that due to the change in $\beta$ to $\beta+\delta$
(with $\delta>0$) and the dependence of the muon count on $A^{1-\beta}$ in the
superposition model (see Eq. (3)), a $c$-fold increase in muon count for
proton showers translates to only a $c\,A^{-\delta}$-fold increase for a
shower of mass number $A$. Since the primary composition in the Auger data is
not pure protons ($A=1$), the necessary increase in muon count to describe the
data (denoted by the gray line in Fig. 7) is approximately 40% for pure proton
showers.

Finally, it is worth emphasizing that the shower-to-shower fluctuations in
$X_{\rm max}$ (shown in the lower panel on the right) remain unaffected by the
enhancements implemented in the S★ variants. Fluctuations in the number of
muons (shown in the panel on the lower left) on the other hand decrease by
approximately 20% in case of S★(mix) and around 10% for S★($\rho^{0}$). For
the pure baryon and kaon enhanced variants the fluctuations decrease only
slightly.

###  3.5 Vertical Showers at Sea Level: Lateral Distribution

Most cosmic ray experiments primarily detect vertical cosmic rays. Ground-
based detectors achieve mass composition sensitivity by separating the muonic
and electromagnetic signal components [47]. As mentioned previously, the
relative contributions of these components reaching the ground depend on the
location of the experiment. As a generic case study, we examine the impact of
the S★ variants on the muon component for vertical showers reaching a sea-
level experiment in the US standard atmosphere.

Figure 8:  Ratio between the average number of muons in iron and proton
showers (left panel) and the slope $\beta$ (right panel), both as functions of
the distance from the shower axis. The inset numbers in the left panel show
the total muon ratio between iron and proton showers (first number) and the
muon density ratio at $1000\,$m (second number in parentheses).

According to the superposition model (Eq. (3)), the ratio between the total
number of muons in a shower of a primary with $A$ nucleons and a proton shower
is $A^{1-\beta}$. Since enhancing muon production via hadronic particle
production requires an increase in the hadronic energy (leading to larger
$\beta$), this inevitably reduces the mass resolution that is the difference
between proton and iron induced showers. For Sibyll 2.3d with $\beta$ around
0.93, the separation between proton and iron in the number of muons is close
to 30%. For the S★ variants, it can be as low as 17% (see legend in Fig. 8).
Note that this connection between muon production and mass resolution does not
apply to modifications affecting only the first high-energy interactions since
such changes would not influence the development of the hadronic cascade
($\beta$), only its initial conditions. However, this scenario is constrained
or even excluded by: 1) Current model predictions matching the fluctuations in
the number of muons observed at the Pierre Auger Observatory [36]; 2) The muon
deficit being observed across several orders of magnitude in energy [2]. That
is to say estimates of the mass resolution based on the total number of muons
with models that produce less muons than are seen in the data are likely an
overestimation.

The situation may not be as dramatic for specific experiments, however. As the
different variants alter the muon energy spectrum, the lateral distribution of
muons reaching the ground will change slightly. Notably, $\beta$ or the mass
separation will vary with the distance from the core. This effect is shown in
Fig. 8. At around $1000\,$m, the separation for Sibyll 2.3d is already around
50%, and for the variants, it increases to around 30% (with a correspondingly
lower $\beta$). The reason for this difference as a function of lateral
distance is that different phases of the shower development dominate at
different distances [22, 44, 45, 48, 49].

Figure 9: Lateral distribution of muons at the ground. Proton primaries are
shown in the left panel, iron primaries are shown on the right. The ratios in
the lower panels are with respect to proton primaries for Sibyll 2.3d. Between
the variants of S★ and Sibyll 2.3d there is only a slight variation in the
shape of the LDF. The largest difference in the shape is between Sibyll 2.3d
and EPOS-LHC.

In Fig. 9 the lateral distribution of muons and the ratio of the distributions
with regard to Sibyll 2.3d are shown for proton and iron primaries. The
largest difference in shape (slope) are between Sibyll 2.3d and EPOS-LHC [50].
The variants of S★ and QGSJet-II 04 [51] are all remarkably similar in shape
to Sibyll 2.3d. The largest difference between the models is by far in the
absolute number of muons in the relevant lateral distance range of about $300$
to $3000$ m.

###  3.6 Vertical Showers at the South Pole

The lowest-energy measurement of the muon deficit comes from IceTop, the
surface detector of the IceCube neutrino observatory, measuring muon densities
from $3\,$PeV to $100\,$PeV [21, 16].

Figure 10:  Ratio of muon counts for vertical proton showers at the South Pole
for different S★ variants relative to Sibyll 2.3d. Full symbols are muon
densities at $600\,$m lateral distance (surface muons). Empty symbols are
total muon counts with energies above $0.5\,$TeV (in-ice muons). The left
panel corresponds to $1\,$PeV primary protons, and the right panel to
$10\,$PeV primary protons.

Figure 10 presents the ratio of muon counts for the S★ variants relative to
Sibyll 2.3d for vertical proton showers at the South Pole. Full symbols
represent muon densities at $600\,$m lateral distance (surface muons), while
empty symbols indicate total muon counts with energies above $0.5\,$TeV (in-
ice muons).

For S★ variants that align with Auger data at $10\,$EeV (the $\rho^{0}$ and
mixed variants), surface muons at $1\,$PeV increase by 10-15%, while in-ice
muons remain unchanged. The situation differs at $10\,$PeV. Here, the
$\rho^{0}$ variant exhibits a roughly 15% increase in both surface and in-ice
muons, while the mixed variant shows a 20% increase in surface muons but only
a 5% increase in in-ice muons.

The analysis by the IceCube collaboration found that Sibyll 2.1 is consistent
with both surface and in-ice muon data, assuming the mass composition from the
Global Spline Fit (GSF) cosmic ray flux model [21, 16, 52]. However, Sibyll
2.3d predicts a roughly 20% increase in surface muons for proton primaries
compared to Sibyll 2.1. While this increase is smaller for a mixed
composition, it is unlikely that Sibyll 2.3d describes the surface data. For
in-ice muons, Sibyll 2.1 and Sibyll 2.3d predict similar numbers. The S★
variants, with a 30-40% increase in surface muons relative to Sibyll 2.1, are
even less likely to describe the surface and in-ice data consistently. Note
however that the effect of snow and ice on the muons at the South Pole is
complex so our definition of in-ice and surface muons only approximately
corresponds to the quantities measured by IceCube.

##  4 Predictions for Inclusive Fluxes and Underground Muons

In this section, we study what implications the S★ variants have for inclusive
atmospheric muon and neutrino flux calculations, as well as muon bundles
observed in deep underground and Cherenkov detectors in water or ice.
Calculations were performed using the MCEq code [53, 54] with DPMJET-III-19.3
[55] as the low-energy interaction model, potentially leading to slight
variations compared to CORSIKA calculations. The Global Spline Fit (GSF)
primary flux model [52] is used throughout this work to emphasize differences
induced by S★ variants.

Sibyll 2.3d is extensively used for the modeling of high-energy atmospheric
neutrino fluxes in IceCube [56]. Except for an inconsistency in the seasonal
variations of the neutrino flux [57], no significant deviations have been
observed. Since inclusive flux predictions depend equally on the cosmic ray
flux model, interpreting deviations from measured spectra cannot be
unequivocally attributed to the hadronic model alone. Stronger model
discrimination might be achieved by combining multiple observables, including
muon and neutrino spectra, angular distributions, deep underground muon rates,
muon bundle multiplicities, and seasonal variations, potentially shedding
light on the origin of the muon excess.

###  4.1 Inclusive Fluxes

Figure 11: Inclusive muon and neutrino fluxes calculated for a vertical zenith
angle. These are conventional fluxes predominantly originating from decays of
charged pions and kaons. Notably, S★ modifications have a negligible impact,
except for the S★(K±,0) model’s enhancement of neutrino fluxes due to
increased kaon decays. Calculations at inclined zenith angles yield
qualitatively similar results with minor differences. The data-driven
daemonflux model [58], derived from accelerator, cosmic-ray, and surface muon
data, is included for comparison, with error estimates extrapolated to higher
energies.  Figure 12: Calculated zenith distributions of tracks (mostly muon
neutrino charged-current events) and electron neutrino events contained within
the instrumented volume of IceCube, using effective areas from public data
releases [59, 60] scaled to 10 years of livetime. No specific zenith
dependence is observed between S★ versions, while daemonflux predicts more
events from the horizon due to increased neutrinos from muon decay. As
expected from Fig. 11, the S★(K±,0) model predicts 15-20% higher neutrino
rates.

The prediction for the inclusive spectra of muons, muon neutrinos and electron
neutrinos (and their respective anti-particles) are shown in Fig. 11. Apart
from the S★(K±,0) model, the predictions from the S★ variants demonstrate no
significant difference compared to Sibyll 2.3d. To some extent this result is
expected, since inclusive fluxes are dominated by the decay of (charged) pions
and kaons produced in the first interaction and less so by particles produced
in the sub-cascades as it is the case in air showers. Therefore S★($\rho^{0}$)
and S★($\bar{p}$), which only lead to small changes in the production of pions
in individual interactions at these energies, show no effect. In S★(K±,0), the
enhanced strangeness production in proton-air interactions leads to more kaon
parents decaying into neutrinos, but with little effect on muon fluxes since
these, at the energies shown here, predominantly originate from pions. The
gray line and band in Fig. 11 is the prediction by a data-driven model, the
so-called daemonflux. This model is fit to various muon measurements performed
at the surface, taking into account systematic uncertainties and constraints
from lab. measurements [58]. In the energy range from $10\,$GeV to $1\,$TeV it
represents the “world data” on inclusive fluxes to well within 10%. The muon
fluxes in all Sibyll models (and most other interaction models [61]) are about
30% lower than the daemonflux model that is than the experimentally observed
flux of muons. A preliminary KM3NeT Collaboration study observes a similar 30%
deviation between the underwater reconstructed muon rate and Sibyll 2.3d [62].
A too low muon prediction is also supported by comparisons with underground
laboratory data [63]. However, daemonflux only incorporates surface muon data
up to a few TeV and predictions in the hundreds TeV range are affected by
large extrapolation errors. Whether the discrepancy observed here in the
inclusive flux of muons at low energies is related to the muon excess observed
in air showers at ultra-high energies is not evident. While both require an
increased production of pions to achieve consistency between data and
predictions, they are sensitive to different of parts of the atmospheric
cascades. For one, as shown in Fig. 11, the S★ models, which do solve the muon
puzzle at ultra-high energies, are not the answer to the atmospheric muon
puzzle.

For both neutrino flavors, all Sibyll variants and the daemonflux model agree
well within errors. A similar observation holds for the zenith distribution of
predicted neutrino events in IceCube (Fig. 12). No significant deviation from
the zenith distribution predicted by Sibyll 2.3d is observed, except for a
15-20% higher neutrino rate in the S★(K±,0) case. Some of the IceCube
$\nu_{\mu}$-focused analyses [64, 65] have hinted at a higher atmospheric
neutrino flux normalization, however no specific measurement has been
performed. The main difference from the daemonflux model is observed at the
horizon, where the muon production distance is far enough from the surface for
decay into neutrinos even at TeV energies to occur, so the higher muon flux
estimate propagates into neutrinos.

###  4.2 Underground Muons

Figure 13: Vertical equivalent muon fluxes in water, obtained using the MUTE
code, which combines MCEq and PROPOSAL codes for underground muon flux
calculations [63].  Figure 14: Muon bundle multiplicity distribution in water
for different slant depths. The zenith angle at the surface has been
accordingly adjusted as $\theta=\arccos(1{\rm km.w.e.}/X)$. The distribution
of muon bundle multiplicities (left panel) strongly depends on the primary
cosmic ray flux model. When scaled to a reference model (here Sibyll 2.3d, see
right panels), differences between the S★ emerge, which become more pronounced
in the rate of very high multiplicity bundles with $\sim 100$ muons.

Muons observed in laboratories deep underground, in water, or in ice exhibit
significant energy loss. The overburden acts as a high-pass filter,
selectively absorbing muons with lower energies ($<1\,$TeV). This makes them
valuable probes for the high-energy portion of the surface spectrum,
potentially distinguishing between different mechanisms that enhance muon
numbers in S★ variants. Figure 13 presents calculated vertical equivalent muon
fluxes in water, obtained using the MUTE code, which combines MCEq and
PROPOSAL [66] codes for underground muon flux calculations [63]. The vertical
equivalent intensity refers to the rate of single muons above a low energy
threshold. Due to median muon energies being in the hundreds of GeV range, the
threshold has minimal impact. While S★ models exhibit more variation here than
at the surface, the changes are still too subtle to be discerned given
uncertainties in interaction and cosmic ray flux models. As expected from
surface flux comparisons (Fig. 11), the data-driven daemonflux model predicts
a 20-30% higher flux, resembling the recent preliminary results from KM3NeT’s
ORCA and ARCA detectors [62].

Figure 14 shows distributions of average muon bundle multiplicities,
representing the muonic cores of air showers that reach specific depths
underwater or in rock. Muon bundle multiplicity distributions strongly depend
on the primary cosmic ray spectrum and mass composition, particularly around
the energy range of the CR knee [67]. However, with better composition
constraints from other measurements, ratios between low and high-multiplicity
events at fixed zenith angles (or slant depths) could serve as a proxy to
discriminate between S★ models, complementing surface muon multiplicities. At
large depths and/or inclined zenith angles the differences between S★ and the
default Sibyll 2.3d reach up to 50% or more at near-horizontal incidence
angles.

While muon bundles can be studied in ice [68], water-based detectors like
KM3NeT [69] and Baikal-GVD [70] might offer superior multiplicity separation
due to a larger fraction of direct Cherenkov light. Upcoming, technologically
advanced, large-volume neutrino telescopes, such as P-ONE [71] and TRIDENT
[72], are anticipated to significantly enhance bundle multiplicity
measurements. Furthermore, large scale radio arrays in the future may provide
the opportunity to directly constrain the flux of atmospheric muons at PeV
energies [73].

##  5 Summary

The new Sibyll★ models are a set of phenomenological modifications to the
well-known Sibyll 2.3d hadronic interaction model. They aim to increase the
number of muons in extensive air showers (EAS) to address the ”muon puzzle.”
These modifications increase the fraction of energy retained in the hadronic
cascade, leading to more meson production and, consequently, more muons. At
the same time, predictions for other air shower observables are intentionally
not changed, or changes kept to a minimum.

The models are exploring different conventional mechanisms for increasing muon
production: enhanced $\rho^{0}$ production, increased baryon–antibaryon pair
production, and increased kaon production. A mixed variant is integrating both
$\rho^{0}$ and baryon pair production enhancements as a baseline approach for
comparison to data. The different variants have been adjusted to describe the
data on particle production that are available from accelerator experiments,
while also achieving a desired increase in the hadronic energy fraction at
ultra-high energies.

We observe that the S★ variants can increase the number of muons in EAS by up
to 35%, depending on the variant and shower configuration. This occurs with a
very small impact on the average depth of shower maximum $\langle X_{\rm
max}\rangle$, typically less than $7\,$g$/$cm2. The $\rho^{0}$ and mixed
variants can provide a sufficient increase in muon number to match
measurements from the Pierre Auger Observatory, but are unlikely to improve
compatibility with IceCube muon data.

Notably, modifications that address the muon deficit imply a change in the
slope parameter $\beta$ that describes the energy dependence of the muon
number. This is understood within the Heitler-Matthews model and confirmed by
our numerical studies. A direct consequence is a reduced sensitivity of the
muon number of showers as mass sensitive observable. However, a good mass
separation capability can be recovered at larger lateral distances from the
shower core.

Regarding inclusive muon and neutrino fluxes, the S★ modifications have little
effect on inclusive muon and neutrino energy spectra compared to Sibyll 2.3d.
The only exception is the S★(K${}^{\pm,0})$ variant, which increases neutrino
fluxes due to the enhanced kaon production and subsequent decays. The
predicted zenith distributions of tracks and cascades in IceCube show no
specific dependence on S★ versions. Therefore, given that the IceCube
observations are consistent with the prediction by Sibyll 2.3d, the
predictions by all the S★ variants will similarly be consistent.

While the changes in the surface flux are subtle, the S★ variants exhibit
larger variations in the underground muon fluxes, especially at high depths
and inclined zenith angles. S★ models predict different distributions of muon
bundle multiplicities in water or ice, offering a potential discriminator
between them. The next generation large-volume neutrino telescopes could
significantly improve measurements of these multiplicities.

All-in-all, the S★ models provide a promising framework for addressing the
muon puzzle in EAS. It allows a more consistent description of EAS
observables, for high energy cosmic rays and neutrino observatories. By
representing accurately the correlations between observables in the same event
and for different events, in a coherent way conserving energy-momentum and
most quantum-numbers, it allows the derivation of full MC sets to train
machine learning algorithms. The different variants of S★ induce similar
changes for the number of muons and the shower maximum but can be
distinguished by other measurements.

The code is available as part of the chromo package333https://github.com/impy-
project/chromo [74]. A standalone version is available upon request from one
of the authors (FR).

#### Acknowledgements

The authors acknowledge many fruitful discussions with colleagues of the
Pierre Auger and IceCube Collaborations. We are grateful to Jenni Adams, Sofia
Andringa and Ruben Conceição for comments on the manuscript of this article.
FR and RE are supported in part by the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No.
101065027 and by BMBF grant No. 05A2023VK4, respectively. The authors jointly
acknowledge the valuable computing resources provided by the Academia Sinica
Grid Computing Center (ASGC), supported by the Institute of Physics at
Academia Sinica. AF is supported by the Academia Sinica Grand Challenge Seed
Program Grant AS-GCS-113-M04.

#### Dedication

The authors dedicate this work to their late colleague Thomas K. Gaisser, with
whom they have enjoyed many years of fruitful collaboration on developing the
hadronic interaction model Sibyll and on calculating predictions for
astroparticle physics applications.

#### Declaration of generative AI and AI-assisted technologies in the writing
process

Generative AI has been used to improve text flow and grammar of parts of the
existing text.

## References

  * [1] J. Albrecht, et al., The Muon Puzzle in cosmic-ray induced air showers and its connection to the Large Hadron Collider, Astrophys. Space Sci. 367 (3) (2022) 27.  arXiv:2105.06148, doi:10.1007/s10509-022-04054-5.
  * [2] J. C. Arteaga-Velzquez, An update on the combined analysis of muon data from nine air-shower experiments at cosmic-ray energies above 1 PeV, PoS ICRC (2023) 466.
  * [3] S. Ostapchenko, QGSJET-II: physics, recent improvements, and results for air showers, EPJ Web Conf. 52 (2013) 02001.  doi:10.1051/epjconf/20125202001.
  * [4] H.-J. Drescher, Remnant break-up and muon production in cosmic ray air showers, Phys. Rev. D 77 (2007) 056003.
  * [5] T. Pierog, K. Werner, Muon production in extended air shower simulations, Phys. Rev. Lett. 101 (2008) 171101.
  * [6] P. K. F. Grieder, The effect of $n\bar{n}$-production on particle spectra in vertically incident and inclined showers derived from simulations, Proc. of 13th Int. Cosmic Ray Conf., Denver 4 (1973) 2467.
  * [7] S. Baur, H. Dembinski, M. Perlin, T. Pierog, R. Ulrich, K. Werner, Core-corona effect in hadron collisions and muon production in air showers, Phys. Rev. D 107 (9) (2023) 094031.  arXiv:1902.09265, doi:10.1103/PhysRevD.107.094031.
  * [8] G. R. Farrar, J. D. Allen, A new physical phenomenon in ultra-high energy collisions, EPJ Web Conf. 53 (2013) 07007.  arXiv:1307.2322, doi:10.1051/epjconf/20135307007.
  * [9] J. Alvarez-Muñiz, L. Cazon, R. Conceição, J. D. de Deus, C. Pajares, M. Pimenta, Muon production and string percolation effects in cosmic rays at the highest energies (2012).  arXiv:1209.6474.
  * [10] M. Rybczyński, Z. Włodarczyk, Puzzle of muons in extensive air showers, Int. J. Mod. Phys. D 28 (08) (2019) 1950097.  arXiv:1902.09988, doi:10.1142/S0218271819500974.
  * [11] L. A. Anchordoqui, H. Goldberg, T. J. Weiler, Strange fireball as an explanation of the muon excess in Auger data, Phys. Rev. D 95 (6) (2017) 063005\.  arXiv:1612.07328, doi:10.1103/PhysRevD.95.063005.
  * [12] L. A. Anchordoqui, C. G. Canal, F. Kling, S. J. Sciutto, J. F. Soriano, An explanation of the muon puzzle of ultrahigh-energy cosmic rays and the role of the Forward Physics Facility for model improvement, JHEAp 34 (2022) 19–32.  arXiv:2202.03095, doi:10.1016/j.jheap.2022.03.004.
  * [13] J. Manshanden, G. Sigl, M. V. Garzelli, Modeling strangeness enhancements to resolve the muon excess in cosmic ray extensive air shower data, JCAP 02 (2023) 017.  arXiv:2208.04266, doi:10.1088/1475-7516/2023/02/017.
  * [14] A. Aab, et al., Deep-learning based reconstruction of the shower maximum $X_{max}$ using the water-Cherenkov detectors of the Pierre Auger Observatory, JINST 16 (07) (2021) P07019.  arXiv:2101.02946, doi:10.1088/1748-0221/16/07/P07019.
  * [15] A. Aab, et al., Extraction of the muon signals recorded with the surface detector of the Pierre Auger Observatory using recurrent neural networks, JINST 16 (07) (2021) P07016.  arXiv:2103.11983, doi:10.1088/1748-0221/16/07/P07016.
  * [16] S. Verpoest, Multiplicity of TeV muons in extensive air showers detected with IceTop and IceCube, PoS ICRC2023 (2023) 207.  arXiv:2307.14689, doi:10.22323/1.444.0207.
  * [17] R. Conceição, B. S. González, A. Guillén, M. Pimenta, B. Tomé, Muon identification in a compact single-layered water Cherenkov detector and gamma/hadron discrimination using machine learning techniques, Eur. Phys. J. C 81 (6) (2021) 542.  arXiv:2101.10109, doi:10.1140/epjc/s10052-021-09312-4.
  * [18] A. Aab, et al., Reconstruction of inclined air showers detected with the Pierre Auger Observatory, JCAP 1408 (08) (2014) 019.  arXiv:1407.3214, doi:10.1088/1475-7516/2014/08/019.
  * [19] A. Aab, et al., Muons in air showers at the Pierre Auger Observatory: Mean number in highly inclined events, Phys. Rev. D 91 (3) (2015) 032003, [Erratum: Phys. Rev.D91,no.5,059901(2015)].  arXiv:1408.1421, doi:10.1103/PhysRevD.91.059901,10.1103/PhysRevD.91.032003.
  * [20] R. U. Abbasi, et al., Study of muons from ultrahigh energy cosmic ray air showers measured with the Telescope Array experiment, Phys. Rev. D 98 (2) (2018) 022002.  arXiv:1804.03877, doi:10.1103/PhysRevD.98.022002.
  * [21] R. Abbasi, et al., Density of GeV muons in air showers measured with IceTop, Phys. Rev. D 106 (3) (2022) 032010.  arXiv:2201.12635, doi:10.1103/PhysRevD.106.032010.
  * [22] F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, T. Stanev, Hadronic interaction model SIBYLL 2.3d and extensive air showers, Phys. Rev. D 102 (2020) 063002.  doi:10.1103/PhysRevD.102.063002.
URL https://link.aps.org/doi/10.1103/PhysRevD.102.063002

  * [23] O. Adriani, et al., The LHCf detector at the CERN Large Hadron Collider, JINST 3 (2008) S08006.  doi:10.1088/1748-0221/3/08/S08006.
  * [24] A. Tiberio, The LHCf experiment at the Large Hadron Collider: status and prospects, PoS ICRC2023 (2023) 444.  doi:10.22323/1.444.0444.
  * [25] G. Piparo, Measurement of the very forward $\pi^{0}$ and $\eta$ meson productions in p-p collisions at $\sqrt{s}$=13 TeV with the LHCf detector, PoS ICRC2023 (2023) 447.  doi:10.22323/1.444.0447.
  * [26] O. Adriani, et al., Measurement of energy flow, cross section and average inelasticity of forward neutrons produced in $\sqrt{s}$ = 13 TeV proton-proton collisions with the LHCf Arm2 detector, JHEP 07 (2020) 016.  arXiv:2003.02192, doi:10.1007/JHEP07(2020)016.
  * [27] H. Abreu, et al., First Direct Observation of Collider Neutrinos with FASER at the LHC, Phys. Rev. Lett. 131 (3) (2023) 031801.  arXiv:2303.14185, doi:10.1103/PhysRevLett.131.031801.
  * [28] R. Albanese, et al., Observation of Collider Muon Neutrinos with the SND@LHC Experiment, Phys. Rev. Lett. 131 (3) (2023) 031802.  arXiv:2305.09383, doi:10.1103/PhysRevLett.131.031802.
  * [29] J. L. Feng, et al., The Forward Physics Facility at the High-Luminosity LHC, J. Phys. G 50 (3) (2023) 030501.  arXiv:2203.05090, doi:10.1088/1361-6471/ac865e.
  * [30] F. Kling, L. J. Nevay, Forward neutrino fluxes at the LHC, Phys. Rev. D 104 (11) (2021) 113008.  arXiv:2105.08270, doi:10.1103/PhysRevD.104.113008.
  * [31] A. Aduszkiewicz, et al., Measurement of meson resonance production in $\pi^{-}+$ C interactions at SPS energies, Eur. Phys. J. C77 (9) (2017) 626.  arXiv:1705.08206, doi:10.1140/epjc/s10052-017-5184-z.
  * [32] H. Adhikary, et al., Measurement of hadron production in $\pi$–C interactions at 158 and 350 GeV/c with NA61/SHINE at the CERN SPS, Phys. Rev. D 107 (6) (2023) 062004.  arXiv:2209.10561, doi:10.1103/PhysRevD.107.062004.
  * [33] A. M. Sirunyan, et al., Measurement of charged pion, kaon, and proton production in proton-proton collisions at sqrt(s) = 13 TeV (2017).  arXiv:1706.10194.
  * [34] S. Chatrchyan, et al., Study of the inclusive production of charged pions, kaons, and protons in $pp$ collisions at $\sqrt{s}=0.9$, 2.76, and 7 TeV, Eur.Phys.J. C72 (2012) 2164.  arXiv:1207.4724, doi:10.1140/epjc/s10052-012-2164-1.
  * [35] E. Albini, P. Capiluppi, G. Giacomelli, A. M. Rossi, Mean charged hadron multiplicities in high-energy collisions, Nuovo Cim. A32 (1976) 101.  doi:10.1007/BF02730044.
  * [36] A. Aab, et al., Measurement of the Fluctuations in the Number of Muons in Extensive Air Showers with the Pierre Auger Observatory, Phys. Rev. Lett. 126 (15) (2021) 152002.  arXiv:2102.07797, doi:10.1103/PhysRevLett.126.152002.
  * [37] M. Adamus, et al., Inclusive production of pi0 mesons in pi+ p, k+ p, and p p interactions at 250-gev/c, Z. Phys. C35 (1987) 7.
  * [38] N. M. Agababyan, et al., Inclusive production of vector mesons in pi+ p interactions at 250-gev/c, Z.Phys. C46 (1990) 387–395.  doi:10.1007/BF01621026.
  * [39] M. Aguilar-Benitez, W. W. W. Allison, A. A. Batalov, E. Castelli, P. Cecchia, et al., Inclusive particle production in 400-gev/c p p interactions, Z.Phys. C50 (1991) 405–426.  doi:10.1007/BF01551452.
  * [40] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, T. Thouw, CORSIKA: a Monte Carlo code to simulate extensive air showers, Wissenschaftliche Berichte, Forschungszentrum Karlsruhe FZKA 6019 (1998).
  * [41] T. T. Böhlen, F. Cerutti, M. P. W. Chin, A. Fassò, A. Ferrari, P. G. Ortega, A. Mairani, P. R. Sala, G. Smirnov, V. Vlachoudis, The FLUKA Code: Developments and Challenges for High Energy and Medical Applications, Nuclear Data Sheets 120 (2014) 211–214.  doi:10.1016/j.nds.2014.07.049.
  * [42] J. Matthews, A Heitler model of extensive air showers, Astropart. Phys. 22 (2005) 387–397.
  * [43] R. Engel, D. Heck, T. Pierog, Extensive air showers and hadronic interactions at high energy, Ann. Rev. Nucl. Part. Sci. 61 (2011) 467–489.  doi:10.1146/annurev.nucl.012809.104544.
  * [44] L. Cazon, R. Conceição, M. Pimenta, E. Santos, A model for the transport of muons in extensive air showers, Astropart. Phys. 36 (2012) 211–223.  arXiv:1201.5294, doi:10.1016/j.astropartphys.2012.05.017.
  * [45] L. Cazon, R. Conceição, F. Riehn, Universality of the muon component of extensive air showers, JCAP 03 (2023) 022.  arXiv:2210.13407, doi:10.1088/1475-7516/2023/03/022.
  * [46] J. Bellido, Depth of maximum of air-shower profiles at the Pierre Auger Observatory: Measurements above $10^{17.2}$ eV and Composition Implications, PoS ICRC2017 (2018) 506, [,40(2017)].  doi:10.22323/1.301.0506.
  * [47] A. Aab, et al., The Pierre Auger Observatory Upgrade - Preliminary Design Report (2016).  arXiv:1604.03637.
  * [48] S. Müller, R. Engel, T. Pierog, M. Roth, Impact of muon detection thresholds on the separability of primary cosmic rays, Astroparticle Physics 97 (2018) 174–185.  doi:https://doi.org/10.1016/j.astropartphys.2017.11.005.
URL https://www.sciencedirect.com/science/article/pii/S0927650517301305

  * [49] I. C. Mariş, et al., Influence of low energy hadronic interactions on air- shower simulations, Nucl. Phys. Proc. Suppl. 196 (2009) 86–89.
  * [50] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, K. Werner, EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider, Phys. Rev. C92 (3) (2015) 034906.  arXiv:1306.0121, doi:10.1103/PhysRevC.92.034906.
  * [51] S. Ostapchenko, Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model, Phys. Rev. D83 (2011) 014018.  arXiv:1010.1869, doi:10.1103/PhysRevD.83.014018.
  * [52] H. P. Dembinski, R. Engel, A. Fedynitch, T. Gaisser, F. Riehn, T. Stanev, Data-driven model of the cosmic-ray flux and mass composition from 10 GeV to $10^{11}$ GeV, PoS ICRC2017 (2018) 533.  arXiv:1711.11432, doi:10.22323/1.301.0533.
  * [53] A. Fedynitch, R. Engel, T. K. Gaisser, F. Riehn, T. Stanev, Calculation of conventional and prompt lepton fluxes at very high energy, EPJ Web Conf. 99 (2015) 08001.  arXiv:1503.00544, doi:10.1051/epjconf/20159908001.
  * [54] A. Fedynitch, F. Riehn, R. Engel, T. K. Gaisser, T. Stanev, Hadronic interaction model SIBYLL $2.3\mathrm{c}$ and inclusive lepton fluxes, Phys. Rev. D 100 (2019) 103018.  doi:10.1103/PhysRevD.100.103018.
URL https://link.aps.org/doi/10.1103/PhysRevD.100.103018

  * [55] A. Fedynitch, Cascade equations and hadronic interactions at very high energy, Ph.D. thesis, KIT (Nov. 2015).
  * [56] M. G. Aartsen, et al., Searches for Sterile Neutrinos with the IceCube Detector, Phys. Rev. Lett. 117 (7) (2016) 071801.  arXiv:1605.01990, doi:10.1103/PhysRevLett.117.071801.
  * [57] R. Abbasi, et al., Observation of seasonal variations of the flux of high-energy atmospheric neutrinos with IceCube, Eur. Phys. J. C 83 (9) (2023) 777.  arXiv:2303.04682, doi:10.1140/epjc/s10052-023-11679-5.
  * [58] J. P. Yañez, A. Fedynitch, Data-driven muon-calibrated neutrino flux, Phys. Rev. D 107 (12) (2023) 123037.  arXiv:2303.00022, doi:10.1103/PhysRevD.107.123037.
  * [59] M. G. Aartsen, et al., Evidence for Astrophysical Muon Neutrinos from the Northern Sky with IceCube, Phys. Rev. Lett. 115 (8) (2015) 081102.  arXiv:1507.04005, doi:10.1103/PhysRevLett.115.081102.
  * [60] M. G. Aartsen, et al., Atmospheric and astrophysical neutrinos above 1 TeV interacting in IceCube, Phys. Rev. D 91 (2) (2015) 022001.  arXiv:1410.1749, doi:10.1103/PhysRevD.91.022001.
  * [61] A. Fedynitch, F. Riehn, R. Engel, T. K. Gaisser, T. Stanev, The hadronic interaction model Sibyll-2.3c and inclusive lepton fluxes (2018).  arXiv:1806.04140.
  * [62] A. Romanov, P. Kalaczyński, Comparison of the atmospheric muon flux measured by the KM3NeT detectors with the CORSIKA simulation using the Global Spline Fit model, PoS ICRC2023 (2023) 338.  doi:10.22323/1.444.0338.
  * [63] A. Fedynitch, W. Woodley, M.-C. Piro, On the Accuracy of Underground Muon Intensity Calculations, Astrophys. J. 928 (1) (2022) 27.  arXiv:2109.11559, doi:10.3847/1538-4357/ac5027.
  * [64] M. G. Aartsen, et al., Observation and Characterization of a Cosmic Muon Neutrino Flux from the Northern Hemisphere using six years of IceCube data, Astrophys. J. 833 (1) (2016) 3.  arXiv:1607.08006, doi:10.3847/0004-637X/833/1/3.
  * [65] M. G. Aartsen, et al., Searching for eV-scale sterile neutrinos with eight years of atmospheric neutrinos at the IceCube Neutrino Telescope, Phys. Rev. D 102 (5) (2020) 052009.  arXiv:2005.12943, doi:10.1103/PhysRevD.102.052009.
  * [66] M. Dunsch, J. Soedingrekso, A. Sandrock, M. Meier, T. Menne, W. Rhode, Recent Improvements for the Lepton Propagator PROPOSAL, Comput. Phys. Commun. 242 (2019) 132–144.  arXiv:1809.07740, doi:10.1016/j.cpc.2019.03.021.
  * [67] M. Ambrosio, et al., Measurement of the residual energy of muons in the Gran Sasso underground laboratories, Astropart. Phys. 19 (2003) 313–328.  arXiv:hep-ex/0207043, doi:10.1016/S0927-6505(02)00217-7.
  * [68] M. G. Aartsen, et al., Characterization of the Atmospheric Muon Flux in IceCube, Astropart. Phys. 78 (2016) 1–27.  arXiv:1506.07981, doi:10.1016/j.astropartphys.2016.01.006.
  * [69] S. Adrian-Martinez, et al., Letter of intent for KM3NeT 2.0, J. Phys. G 43 (8) (2016) 084001.  arXiv:1601.07459, doi:10.1088/0954-3899/43/8/084001.
  * [70] A. V. Avrorin, et al., Deep-Underwater Cherenkov Detector in Lake Baikal, J. Exp. Theor. Phys. 134 (4) (2022) 399–416.  doi:10.1134/S1063776122040148.
  * [71] M. Agostini, et al., The Pacific Ocean Neutrino Experiment, Nature Astron. 4 (10) (2020) 913–915.  arXiv:2005.09493, doi:10.1038/s41550-020-1182-4.
  * [72] Z. P. Ye, F. Hu, W. Tian, Q. C. Chang, Y. L. Chang, Z. S. Cheng, J. Gao, T. Ge, G. H. Gong, J. Guo, X. X. Guo, X. G. He, J. T. Huang, K. Jiang, P. K. Jiang, Y. P. Jing, H. L. Li, J. L. Li, L. Li, W. L. Li, Z. Li, N. Y. Liao, Q. Lin, J. Lin, F. Liu, J. L. Liu, X. H. Liu, P. Miao, C. Mo, I. Morton-Blake, T. Peng, Z. Y. Sun, J. N. Tang, Z. B. Tang, C. H. Tao, X. L. Tian, M. X. Wang, Y. Wang, Y. Wang, H. D. Wei, Z. Y. Wei, W. H. Wu, S. S. Xian, D. Xiang, D. L. Xu, Q. Xue, J. H. Yang, J. M. Yang, W. B. Yu, C. Zeng, F. Y. D. Zhang, T. Zhang, X. T. Zhang, Y. Y. Zhang, W. Zhi, Y. S. Zhong, M. Zhou, X. H. Zhu, G. J. Zhuang, A multi-cubic-kilometre neutrino telescope in the western Pacific Ocean, Nature Astronomy (Oct. 2023).  doi:10.1038/s41550-023-02087-6.
  * [73] L. Pyras, C. Glaser, S. Hallmann, A. Nelles, Atmospheric muons at PeV energies in radio neutrino detectors, JCAP 10 (2023) 043.  arXiv:2307.04736, doi:10.1088/1475-7516/2023/10/043.
  * [74] H. Dembinski, A. Fedynitch, A. Prosekin, Chromo: An event generator frontend for particle and astroparticle physics, PoS ICRC2023 (2023) 189.  doi:10.22323/1.444.0189.


# Sharp inequalities involving multiplicative chaos sums

Grigori A. Karagulyan  Faculty of Mathematics and Mechanics, Yerevan State
University, Alex Manoogian, 1, 0025, Yerevan, Armenia<EMAIL_ADDRESS>Institute of Mathematics of NAS of RA, Marshal Baghramian ave., 24/5, Yerevan,
0019, Armenia<EMAIL_ADDRESS>
###### Abstract.

The present note is an essential addition to author’s recent paper [Kar4],
concerning general multiplicative systems of random variables. Using some
lemmas and the methodology of [Kar4], we obtain a general extreme inequality,
with corollaries involving Rademacher chaos sums and those analogues for
multiplicative systems. In particular we prove that a system of functions
generated by bounded products of a multiplicative system is a convergence
system.

###### Key words and phrases:

multiplicative system, martingale difference, Khintchine’s inequality, Azuma-
Hoeffding inequality, lacunary subsystem, Rademacher chaos

###### 2010 Mathematics Subject Classification:

42C05, 42C10, 42A55, 60G42, 60G48

The work was supported by the Science Committee of RA, in the frames of the
research project 21AG‐1A045

##  1\. Introduction

A sequence of bounded random variables $\phi_{n}$, $n=1,2,\ldots$ (finite or
infinite) is said to be multiplicative if the equality

(1.1) ${\textbf{E}}\left[\prod_{j\in A}\phi_{j}\right]=0$

holds for all nonempty subsets $A\subset\mathbb{N}$ of positive integers. The
multiplicative systems were introduced by Alexits (see [Alex], chap 3.2), who
also proved in [Alex1] that the uniformly bounded multiplicative systems are
convergence systems, i.e. the condition $\sum_{k}a_{k}^{2}<\infty$ implies
almost sure convergence of series $\sum_{k}a_{k}\phi_{k}$ (see [KoRe, Kom,
Gap2] for further extensions of this result). The multiplicative systems were
considered also in the context of central limit theorems and law of iterated
logarithm ([Mor, Ber, Fuk1, Fuk2]). Well-known examples of multiplicative
sequences are mean zero independent random variables and more general, the
martingale-differences. The trigonometric sequences $\\{\cos
2n_{k}\pi(x+x_{k})\\}$ on $(0,1)$ with positive integers $n_{k}$, satisfying
$n_{k+1}\geq 2n_{k}$, are known to be non-martingale examples of
multiplicative systems (see [Zyg1], chap. 5). A wide class of multiplicative
systems was recently introduced in [Rub]. Namely, under certain symmetry
conditions on $f\in L^{\infty}(\mathbb{R}/\mathbb{Z})$, the functions $f(nx)$,
$n=1,2,\ldots$ form a multiplicative system on $(0,1)$.

Recall the definition of multiplicative error given in [Kar4]. Denote
$\mathbb{Z}_{n}=\\{1,2,\ldots,n\\}$ and let
$\mathbb{Z}_{\infty}=\\{1,2,\ldots\\}$ be the set of all positive integers. We
consider systems of random variables
$\phi=\\{\phi_{k}:\,k\in\mathbb{Z}_{n}\\}$ (here $n$ can be either finite or
infinite) satisfying

(1.2) $A_{k}\leq\phi_{k}\leq B_{k},\text{ where }A_{k}<0<B_{k}.$

Given $d\leq n\leq\infty$ the $d$-multiplicative error of the system $\phi$ is
the quantity

(1.3) $\mu_{d}(\phi)=\sum_{A\subset\mathbb{Z}_{n},\,\\#(A)\leq
d}\quad\prod_{j\in A}(C_{j})^{-1}\left|{\textbf{E}}\left[\prod_{j\in
A}\phi_{j}\right]\right|,$

where $C_{k}=\min\\{-A_{k},B_{k}\\}$. We say $\phi$ is $d$-multiplicative if
$\mu_{d}(\phi)=0$ or equivalently (1.1) holds for any finite nonempty
$A\subset\mathbb{N}$ of cardinality $\\#(A)\leq d$. If $d=\infty$ this gives
the classical definition of multiplicative systems. It was proved in [Kar4]
(see Lemma 2.1 below) that if the system $\phi$ is defined on $[0,1]$ and
$\mu=\mu_{d}(\phi)<\infty$, then the functions of $\phi$ can be extended up to
$[0,1+\mu]$ to a $d$-multiplicative system, obeying same condition (1.2). We
call a system of random variables $d$-independent if any collection of its $d$
elements is independent. It is clear that $d$-independent mean zero random
variables are $d$-multiplicative.

The main result of [Kar4] reads as follows.

######  Theorem A ([Kar4]).

Let $\Phi:{\mathbb{R}}\to{\mathbb{R}}^{+}$ be a convex function and
$\phi=\\{\phi_{k}:\,k=1,2,\ldots,n\\}$ be a system of random variables
satisfying (1.2). Then for any integer $d\leq n$ and a choice of real
coefficients $a_{1},\ldots,a_{n}$ it holds the inequality

(1.4)
$\textbf{E}\left[\Phi\left(\sum_{k=1}^{n}a_{k}\phi_{k}\right)\right]\leq(1+\mu_{d}(\phi))\textbf{E}\left[\Phi\left(\sum_{k=1}^{n}a_{k}\xi_{k}\right)\right],$

where $\xi_{k}$, $k=1,2,\ldots,n$, are $\\{A_{k},B_{k}\\}$\- valued
$d$-independent mean zero random variables.

Applying this inequality for $d$-multiplicative systems, i.e.
$\mu_{d}(\phi)=0$, with the parameters $A_{k}=-1$ and $B_{k}=1$, we
immediately obtain the following extremal property of Rademacher random
variables.

######  Theorem B ([Kar4]).

Let $\Phi:{\mathbb{R}}\to{\mathbb{R}}^{+}$ be a convex function. If
$\phi=\\{\phi_{k}:\,k=1,2,\ldots,n\\}$, is a multiplicative system of random
variables, satisfying $\|\phi_{k}\|_{\infty}\leq 1$, then for any choice of
coefficients $a_{1},\ldots,a_{n}$ we have

(1.5)
$\textbf{E}\left[\Phi\left(\sum_{k=1}^{n}a_{k}\phi_{k}\right)\right]\leq{\textbf{E}}\left[\Phi\left(\sum_{k=1}^{n}a_{k}r_{k}\right)\right],$

where $r_{k}$, $k=1,2,\ldots,n$, are Rademacher independent random variables.

Applying inequalities (1.4) and (1.5), paper [Kar4] provides extensions of
some probability inequalities for general multiplicative systems. In
particular it was proved the following generalization of Azuma-Hoeffding’s
martingale inequality [Azu, Hoe].

######  Theorem C ([Kar4]).

If a system of random variables $\phi=\\{\phi_{k}:\,k=1,2,\ldots,n\\}$
satisfies (1.2), then it holds the inequality

(1.6)
$\left|\left\\{\sum_{k=1}^{n}\phi_{k}>\lambda\right\\}\right|\leq(1+\mu_{n}(\phi)\exp\left(-\frac{2\lambda^{2}}{\sum_{k=1}^{n}(B_{k}-A_{k})^{2}}\right),\quad\lambda>0.$

The present note is an addition to paper [Kar4]. Essentially using the lemmas
and the methodology of [Kar4] it provides an extension of inequality (1.4) in
more general settings, with corollary inequalities involving Rademacher chaos
and its analogue for multiplicative systems. We say a function
$G:{\mathbb{R}}^{n}\to{\mathbb{R}}^{+}$ is convex, if it is convex with
respect to each variable separately. Hence the main result of the present note
is the following.

######  Theorem 1.1.

Let a system of random variables $\phi=\\{\phi_{k}:\,k=1,2,\ldots,n\\}$
satisfy (1.2), $d\leq n$ and $G:{\mathbb{R}}^{n}\to{\mathbb{R}}^{+}$ be a
convex function. Then it holds the inequality

(1.7)
$\textbf{E}\left[G\left(\phi_{1},\ldots,\phi_{n}\right)\right]\leq(1+\mu_{d}(\phi))\textbf{E}\left[G\left(\xi_{1},\ldots,\xi_{n}\right)\right],$

where $\xi_{k}$, $k=1,2,\ldots,n$, are $\\{A_{k},B_{k}\\}$\- valued
$d$-independent mean zero random variables.

######  Corollary 1.1.

Let $G:{\mathbb{R}}^{n}\to{\mathbb{R}}^{+}$ be a convex function. If
$\phi=\\{\phi_{k}:\,k=1,2,\ldots,n\\}$, is a system of random variables,
satisfying $\|\phi_{k}\|_{\infty}\leq 1$, (i.e. $B_{k}=-A_{k}=1$ in (1.2)),
then we have

(1.8)
$\textbf{E}\left[G\left(\phi_{1},\ldots,\phi_{n}\right)\right]\leq(1+\mu_{d}(\phi))\textbf{E}\left[G\left(r_{1},\ldots,r_{n}\right)\right],$

where $r_{k}$, $k=1,2,\ldots,n$, are Rademacher $d$-independent random
variables.

######  Remark 1.1.

We remark that if $l<n$, then the sequences of $\\{A_{k},B_{k}\\}$\- valued
mean zero $d$-independent random variables can be probabilistically different,
whereas when $d=n$ such a sequence is uniquely determined. So if $d=n$ in
Corollary 1.1, then $r_{k}$, $k=1,2,\ldots,n$, become ordinary Rademacher
independent random variables. In fact, the case $d=n$ is the most interesting
case in Theorem 1.1 and Corollary 1.1.

######  Remark 1.2.

We will mainly be interested in convex functions $G$ obtained as follows. Let
$X$ be a Banach space and $\mathcal{A}\subset 2^{\mathbb{Z}_{n}}$ be a family
of nonempty subsets of $\mathbb{Z}_{n}=\\{1,2,\ldots,n\\}$. To each
$A\in\mathcal{A}$ we assign a coefficient $x_{A}\in F$ and consider the
function

(1.9) $P(\textbf{t})=\sum_{A\in\mathcal{A}}x_{A}\prod_{k\in
A}t_{k}:\mathbb{R}^{n}\to X,$

with respect to variable $\textbf{t}=(t_{1},\ldots,t_{n})\in\mathbb{R}^{n}$.
We will call such functions quasi-polynomials. Given sequence of quasi
polynomials $P_{j}$, $j=1,2,\ldots,m$, and a one variable convex function
$\Phi:\mathbb{R}^{+}\to\mathbb{R}^{+}$. It is an easy task to check that

(1.10) $G(t_{1},\ldots,t_{n})=\Phi\left(\max_{1\leq j\leq
m}\|P_{j}(t_{1},\ldots,t_{n})\|_{X}\right)$

is a convex function on $\mathbb{R}^{n}$.

Let $\phi_{n}$, $n=1,2,\ldots$, be a multiplicative system. For a given
nonempty $A\subset\mathbb{N}$ denote

(1.11) $\phi_{A}(x)=\prod_{k\in A}\phi_{k}(x).$

If $\\{\phi_{n}\\}$ coincides with the Rademacher system, $\phi_{n}=r_{n}$,
then $\\{\phi_{A}\\}$ becomes the classical Walsh system that will be denoted
by $w_{A}(x)$. The subsystem $\\{w_{A}:\,A\subset\mathbb{N},\,\\#(A)=d\\}$ of
Walsh functions, where only $d$-products of Rademacher functions are involved,
is called Rademacher chaos of order $d$. The Rademacher chaos of first order
($d=1$) coincides with the Rademacher functions. The Rademacher chaos systems
are very well investigated in different points of view. It is of interest the
comparisons of the different $L^{p}$-norms of finite sums

(1.12) $S(x)=\sum_{A\subset\mathbb{N},\,\\#(A)=d}b_{A}w_{A}(x),$

where $b_{A}$ is finitely supported sequence of real numbers. There are many
other interesting properties of Rademacher chaos systems that doesn’t have the
complete system of Walsh functions (see [Ast], chap. 6). The following
generalization of Khintchine [Khi] inequality for the Rademacher chaos proved
independently by Bonami [Bon] and Kiener [Kie] is a crucial tool for the
analysis of Boolean functions (see [Don]).

######  Theorem D (Bonami-Kiener).

For any integers $d\geq 2$ and a Rademacher chaos sum (1.12) it holds the
bound

(1.13) $\left\|S\right\|_{p}\leq(p-1)^{d/2}\left\|S\right\|_{2},\quad p>2.$

It was also proved in [Bon] that the constant growth in (1.13) as $p\to\infty$
is optimal in the sense that for any $p>2$ there exists a sum (1.12) such that
$\|S\|_{p}\geq C_{d}p^{d/2}\|S\|_{2}$ with a constant $C_{d}>0$, depending
only on $d$. Using a standard argument one can deduce from inequalities (1.13)
the bound

(1.14) $\left\|S\right\|_{p}\leq C_{p,q,d}\left\|S\right\|_{q},$

where $1\leq q<p<\infty$. It is of interest to find the optimal constants in
(1.14). The case of the classical Rademacher functions (when $d=1$) is well
investigated. The analogue of inequality (1.13) for the simple Rademacher sums
is Khintchine’s classical inequality and there were found the values of all
sharp constants in (1.14) if either $p$ or $q$ is equal $2$ (see [Ste, Sza,
Haa, You]). To the best of our knowledge the sharp values of these constants
if $d\geq 2$ are not known for any combination of the parameters $p>q$. Some
particular estimates of constants (1.14) one can find in papers [Jan, Don,
Lar, Ivan]. The following inequalities are immediately obtained from (1.7).

######  Corollary 1.2.

Let $\phi_{n}$, $n=1,2,\ldots$, be a system of random variables with
$\|\phi_{n}\|_{\infty}\leq 1$, $2\leq d\leq\infty$, and
$\Phi:\mathbb{R}_{+}\to\mathbb{R}_{+}$ be a convex function. If $A_{n}$ is a
sequence of different nonempty finite subsets of $\mathbb{N}$ with
$\\#(A_{n})\leq d$, then for any finitely supported sequence of vectors
$b_{n}$ from a Banach space $X$ there hold the bounds

(1.15)
$\displaystyle\textbf{E}\left[\Phi\left(\left\|\sum_{k}b_{k}\phi_{A_{k}}\right\|_{X}\right)\right]\leq(1+\mu_{d}(\phi))\textbf{E}\left[\Phi\left(\left\|\sum_{k}b_{k}w_{A_{k}}\right\|_{X}\right)\right],$
(1.16)
$\displaystyle\textbf{E}\left[\Phi\left(\max_{n}\left\|\sum_{k}b_{k}\phi_{A_{k}}\right\|_{X}\right)\right]\leq(1+\mu_{d}(\phi))\textbf{E}\left[\Phi\left(\max_{n}\left\|\sum_{k}b_{k}w_{A_{k}}\right\|_{X}\right)\right].$

Letting $\Phi(t)=t^{p}$ in (1.15) and $X$ to be the Banach space of real line
$\mathbb{R}$, we will get inequality (1.13) for more general sums

(1.17) $S(x)=\sum_{A\subset\mathbb{N},\,\\#(A)=d}b_{A}\phi_{A}(x),$

where $\phi_{A}$ are the product random variables (1.11) generated by a
multiplicative systems.

######  Theorem 1.2.

$\\{\phi_{n}:\,n=1,2,\ldots\\}$ with $\|\phi_{n}\|_{\infty}\leq 1$, then for
any integers $d\geq 2$ and any sum (1.17) it holds the bound

(1.18) $\left\|S\right\|_{p}\leq(p-1)^{d/2}\left\|S\right\|_{2},\quad p>2.$

Recall that a system of functions $\\{\phi_{n}\\}$ is said to be a convergence
system if any series

(1.19) $\sum_{n}a_{n}\phi_{n}(x)$

with coefficients satisfying $\sum_{n}a_{n}^{2}<\infty$ converges a.e., if
a.e. convergence holds after any rearrangement of the terms of (1.19), then we
say $\phi_{n}(x)$ is an unconditional convergence system. As it was proved by
Stechkin [Ste], if a system of functions $\phi_{n}\in L^{2}$ satisfies a
Khintchine type inequality

(1.20) $\left\|\sum_{k=1}^{n}a_{k}\phi_{k}\right\|_{p}\leq
c\left(\sum_{k=1}^{n}a_{k}^{2}\right)^{1/2}$

for some $p>2$ and $c>0$, then $\phi_{n}(x)$ is an unconditional convergence
system. Using this result and the bound (1.18) we conclude the following.

######  Corollary 1.3.

Let $\phi=\\{\phi_{n}$, $n=1,2,\ldots\\}$, be an arbitrary system of random
variables with $\|\phi_{n}\|_{\infty}\leq M$ and $\mu_{d}(\phi)<\infty$,
$d\geq 1$. Then for any sequence $A_{n}\subset\mathbb{N}$, $n=1,2,\ldots$, of
different nonempty sets with $\\#(A_{n})\leq d$ the random variables
$\\{\phi_{A_{n}}\\}$ form an unconditional convergence system.

######  Remark 1.3.

An interesting case in Corollary 1.3 is when $\\{\phi_{k}\\}$ is a bounded
martingale difference system, that is
$\textbf{E}\left(\phi_{n}|\phi_{1},\ldots,\phi_{n-1}\right)=0,\quad
n=2,3,\ldots.$

######  Remark 1.4.

Inequalities (1.15) and (1.16) show certain extremal properties of the
Rademacher chaos or the Walsh functions. Some extremal properties of simple
Rademacher functions one can find in [Ast], chap. 7 and 9. Those are mostly
concerned to independent random variables and no one consider Rademacher chaos
or maximal function of partial sums as we have in (1.15) or (1.16).

##  2\. Proof of Theorem 1.1

A real-valued function $f$ defined on $[a,b)$ is said to be a step function if
it can be written as a finite linear combination of indicator functions of
intervals $[\alpha,\beta)\subset[a,b)$. The following lemmas were proved in
[Kar4].

######  Lemma 2.1 ([Kar4], Lemma 2.2).

Let $\phi=\\{\phi_{k}:\,k=1,2,\ldots,n\\}$ be a system of measurable functions
on $[0,1)$ satisfying (1.2) and ${\mathfrak{M}}$ be an arbitrary family of
subsets of ${\mathbb{Z}}_{n}$. Then the functions $\phi_{k}$ can be extended
up to the interval $[0,1+\mu)$ such that

(2.1) $\displaystyle A_{k}\leq\phi_{k}(x)\leq B_{k}\text{ if }x\in[0,1+\mu),$
(2.2) $\displaystyle\int_{0}^{1+\mu}\prod_{j\in A}\phi_{j}=0\text{ for all
}A\in{\mathfrak{M}},$

where

(2.3) $\mu_{\mathfrak{M}}(\phi)=\sum_{A\in\mathfrak{M}}\quad\prod_{j\in
A}(C_{j})^{-1}\left|{\textbf{E}}\left[\prod_{j\in A}\phi_{j}\right]\right|,$

is the multiplicative error corresponding to the family $\mathfrak{M}$.
Moreover, each $\phi_{k}$ is a step function on $[1,1+\mu)$.

######  Lemma 2.2 ([Kar4], Lemma 2.3).

Let measurable functions $\phi_{k}$, $k=1,2,\ldots,n$, defined on $[0,1)$,
satisfy (1.2). Then for any $\delta>0$ one can find step functions $f_{k}$,
$k=1,2,\ldots,n$, on $[0,1)$ with $A_{k}\leq f_{k}\leq B_{k}$ such that

(2.4) $|\\{|\phi_{k}-f_{k}|>\delta\\}|<\delta,\quad k=1,2,\ldots,n,$

and

(2.5) $\int_{0}^{1}f_{n_{1}}f_{n_{2}}\ldots f_{n_{\nu}}=0\text{ as
}\\{n_{1},n_{2},\ldots,n_{\nu}\\}\in{\mathfrak{M}},$

where ${\mathfrak{M}}$ is the family of finite nonempty sets
$A\subset\mathbb{N}$, satisfying (1.1).

######  Lemma 2.3 ([Kar4], Lemma 2.4).

If $g_{k}$, $k=1,2,\ldots,n$, is a $d$-multiplicative system of nonzero random
variables such that each $g_{k}$ takes two values, then $g_{k}$ are
$d$-independent.

The following lemma is a version of Lemma 2.5 of [Kar4] with the same proof.
For the sake of completeness its complete proof will be provided.

######  Lemma 2.4.

Let $f_{1},f_{2},\ldots,f_{n}$ be real-valued step functions on $[0,1)$
satisfying $A_{k}\leq f_{k}(x)\leq B_{k}$ and let
$G:{\mathbb{R}}^{n}\to{\mathbb{R}}^{+}$ be a convex function. Then there are
$\\{A_{k},B_{k}\\}$-valued step functions $\xi_{1},\xi_{2},\ldots,\xi_{n}$ on
$[0,1)$ such that

(2.6)
$\int_{0}^{1}\xi_{n_{1}}\xi_{n_{2}}\ldots\xi_{n_{\nu}}=\int_{0}^{1}f_{n_{1}}f_{n_{2}}\ldots
f_{n_{\nu}},$

for any choice of
$\\{n_{1},n_{2},\ldots,n_{\nu}\\}\subset{\mathbb{Z}}_{n}=\\{1,2,\ldots,n\\}$,
and it holds the inequality

(2.7)
$\int_{0}^{1}G\left(f_{1},\ldots,f_{n}\right)\leq\int_{0}^{1}G\left(\xi_{1},\ldots,\xi_{n}\right).$

###### Proof.

Let $\Delta_{j}$, $j=1,2,\ldots,m$ be the constancy intervals of functions
$f_{k}$. Let $\Delta=[\alpha,\beta)$ be one of those intervals. Observe that
the point

(2.8) $c=\frac{B_{1}\alpha-
A_{1}\beta}{B_{1}-A_{1}}+\frac{1}{B_{1}-A_{1}}\cdot\int_{\alpha}^{\beta}f_{1}(t)dt$

is in the closure of the interval $\Delta$. Then we define $\xi_{1}$ on
$\Delta$ as

(2.9)
$\xi_{1}(x)=B_{1}\cdot{\textbf{1}}_{[\alpha,c)}(x)+A_{1}\cdot{\textbf{1}}_{[c,\beta)}(x),\quad
x\in\Delta.$

Applying this to each $\Delta_{j}$, we will get a function $\xi_{1}$ defined
on entire $[a,b)$ and one can check

(2.10) $\int_{\Delta_{j}}\xi_{1}(t)dt=\int_{\Delta_{j}}f_{1}(t)dt,\quad
j=1,2,\ldots,m.$

Since each $f_{k}$ is constant on the intervals $\Delta_{j}$,
$j=1,2,\ldots,m$, from (2.10) we conclude that

(2.11) $\int_{0}^{1}\xi_{1}f_{n_{2}}\ldots
f_{n_{\nu}}=\int_{0}^{1}f_{1}f_{n_{2}}\ldots f_{n_{\nu}},$

for any collection of integers $1<n_{2}<\ldots<n_{\nu}\leq n$. We may also
claim that

(2.12)
$\int_{0}^{1}G\left(f_{1},f_{2},\ldots,f_{n}\right)\leq\int_{0}^{1}G\left(\xi_{1},f_{2},\ldots,f_{n}\right).$

Indeed, let us fix an interval $\Delta_{j}$ and suppose that $f_{k}(t)=c_{k}$
on $\Delta_{j}$. Applying (2.10) and the Jessen inequality, we get

(2.13)
$\displaystyle\int_{{\Delta_{j}}}G\left(f_{1}(t),f_{2}(t),\ldots,f_{n}(t)\right)dt$
$\displaystyle=G\left(c_{1},c_{2},\ldots c_{n}\right)|{\Delta_{j}}|$ (2.14)
$\displaystyle=G\left(\frac{1}{|{\Delta_{j}}|}\int_{{\Delta_{j}}}\xi_{1}\,,c_{2},\ldots,c_{n}\right)|{\Delta_{j}}|$
(2.15)
$\displaystyle\leq\int_{{\Delta_{j}}}G\left(\xi_{1}(t),f_{2}(t),\ldots,f_{n}(t)\right)dt,$

then the summation over $j$ implies (2.12). Applying the same procedure to the
new system $\xi_{1},f_{2},\ldots,f_{n}$ we can similarly replace $f_{2}$ by
$\xi_{2}$. Continuing this procedure we will replace all functions $f_{k}$ to
$\xi_{k}$ ensuring the conditions of lemma. ∎

###### Proof of Theorem 1.1.

At first we suppose that $\\{\phi_{k}\\}$ is a $d$-multiplicative system.
Without loss of generality we can suppose that $\phi_{k}$ are defined on
$[0,1)$. Applying Lemma 2.2, we find an $d$-multiplicative system of step
functions $f_{k}$ on $[0,1)$ satisfying

(2.16) $|\\{x\in[0,1):\,|f_{k}(x)-\phi_{k}(x)|>\delta\\}|<\delta.$

Then we apply Lemma 2.4 and get $A_{k},B_{k}$-valued step functions $\xi_{k}$
defined on $[0,1)$ and satisfying (2.6) and (2.7) ($a=0,b=1$). Since
$\\{f_{k}\\}$ is $d$-multiplicative, in view of (2.6) so we will have for
$\\{\xi_{k}\\}$ and then by Lemma 2.3 $\\{\xi_{k}\\}$ forms
$A_{k},B_{k}$-valued $d$-independent random variables. Observe that

(2.17)
$\int_{0}^{1}G\left(\phi_{1},\ldots,\phi_{n}\right)\leq\varepsilon(\delta)+\int_{0}^{1}G\left(f_{1},\ldots,f_{n}\right),$

where by (2.16) we have

(2.18) $\varepsilon=\varepsilon(\delta)\to 0\text{ as }\delta\to 0.$

Here we also use the continuity of the convex function $G$. Thus, from (2.7)
we obtain

(2.19) $\displaystyle\int_{0}^{1}G\left(\phi_{1},\ldots,\phi_{n}\right)$
$\displaystyle\leq\varepsilon+\int_{0}^{1}G\left(f_{1},\ldots,f_{n}\right)\leq\varepsilon+\int_{0}^{1}G\left(\xi_{1},\ldots,\xi_{n}\right).$

Therefore, taking into account (2.18), one can easily get (1.7). Indeed, the
system $\\{\xi_{k}\\}$ in fact depends on $\delta$. So let
$\\{\xi_{k}^{(m)}\\}$ be the system corresponding to $\delta=1/m$. We can
suppose that there is a partition $0=x_{0}^{(m)}\leq x_{1}^{(m)}\leq\ldots\leq
x_{2^{n}}^{(m)}=1$ such that each function $\xi_{k}^{(m)}$ is equal to a
constant (either $A_{k}$ or $B_{k}$) on each interval
$[x_{j-1}^{(m)},x_{j}^{(m)})$, $1\leq j\leq 2^{n}$, and this constant is
independent of $m$. Then we find a sequence $m_{k}$ such that each sequence
$x_{j}^{(m_{k})}$, $k=1,2,\ldots$, is convergent. One can check that this
generates a limit sequence of $d$-multiplicative $\\{A_{k},B_{k}\\}$-valued
random variables $\xi_{k}$, $k=1,2,\ldots,n$, satisfying (1.7).

Now suppose that $\\{\phi_{n}\\}$ is an arbitrary system. Applying Lemma 2.1
with $\mathfrak{M}=\\{A\subset\mathbb{Z}_{n}:\,\\#(A)\leq d\\}$, we may extend
$\\{\phi_{n}\\}$ up to $[0,1+\mu]$, satisfying the conditions of the lemma.
Consider the $d$-multiplicative system $\psi_{n}(x)=g_{n}((1+\mu)x)$ on
$[0,1)$. We have

(2.20) $\displaystyle\int_{0}^{1}G\left(\phi_{1},\ldots,\phi_{n}\right)$
$\displaystyle\leq\int_{0}^{1+\mu}G\left(\phi_{1},\ldots,\phi_{n}\right)$
(2.21)
$\displaystyle=(1+\mu)\int_{0}^{1}G\left(\psi_{1},\ldots,\psi_{n}\right)\leq(1+\mu)\int_{0}^{1}G\left(\xi_{1},\ldots,\xi_{n}\right),$

where $\xi_{n}$ is a $d$-independent system generated by Lemma 2.3 applied on
$\psi_{n}$. This completes the proof of theorem.

∎

##  3\. Trigonometric subsystems

The Walsh functions have a canonical numeration, $w_{n}$, $n=0,1,2,\ldots$.
That is $w_{0}\equiv 1$ and for $n\geq 1$ with a dyadic decomposition

(3.1) $n=2^{k_{1}}+2^{k_{2}}+\cdots+2^{k_{d}},\quad 0\leq
k_{1}<k_{2}<\cdots<k_{d}.$

we have

$w_{n}(x)=r_{k_{1}+1}(x)r_{k_{2}+1}(x)\ldots r_{k_{d}+1}(x).$

We call $d$ a power of the number $n$ and write $\rho(n)=d$. Let
$L_{\Phi}(0,1)$ denote the Orlicz space of functions on $(0,1)$, corresponding
to a Young function $\Phi:\mathbb{R}_{+}\to\mathbb{R}_{+}$ that is a convex
function satisfying

$\lim_{t\rightarrow
0+}\frac{\Phi(t)}{t}=\lim_{t\rightarrow\infty}\frac{t}{\Phi(t)}=0.$

The Luxemburg norm of a function $f\in L_{\Phi}(0,1)$ is defined by

$\|f\|_{\Phi}=\inf\left\\{\lambda:\,\lambda>0,\,\int_{0}^{1}\Phi\left(\frac{|f|}{\lambda}\right)\leq
1\right\\}<\infty.$

######  Corollary 3.1.

If $\Phi:\mathbb{R}_{+}\to\mathbb{R}_{+}$ is a Young function, then for any
finitely supported sequence of coefficients $b_{n}$, $n=1,2,\ldots,$ we have
the bound

(3.2) $\displaystyle\left\|\sum_{n:\,\rho(n)\leq d}b_{n}\cos 2\pi
n(x+x_{n})\right\|_{\Phi}\leq 2^{d-1}\left\|\sum_{n:\,\rho(n)\leq
d}b_{n}w_{n}(x)\right\|_{\Phi},$ (3.3) $\displaystyle\left\|\max_{m\geq
1}\left|\sum_{n\leq m:\,\rho(n)\leq d}b_{n}\cos 2\pi
n(x+x_{n})\right|\right\|_{\Phi}\leq 2^{d-1}\left\|\max_{m\geq
1}\left|\sum_{n\leq m:\,\rho(n)\leq d}b_{n}w_{n}(x)\right|\right\|_{\Phi}$

###### Proof.

Using the trigonometric sums formulae and decomposition (3.1), the function
$\cos 2\pi n(x+x_{n})$ may be written as a sum of $2^{d-1}$ products of the
form

(3.4)
$t_{1}\big{(}2^{k_{1}}(x+\alpha_{1})\big{)}t_{2}\big{(}2^{k_{2}}(x+\alpha_{2})\big{)}\cdots
t_{s}\big{(}2^{k_{d}}(x+\alpha_{d})\big{)},$

where each $t_{j}(x)$ is either $\sin 2\pi x$ or $\cos 2\pi x$. On the other
hand for any choice of trigonometric functions $t_{j}(x)$ and numbers
$\alpha_{j}$ the system $t_{j}(2^{j}(x+\alpha_{j}))$ is multiplicative and
obeys the conditions of Corollary 1.2. Thus, inequalities (1.15) and (1.16)
imply (3.2) and (3.3) respectively. ∎

Then using the Bonami-Kiener inequality (1.13), from (3.2) we immediately
obtain the following

######  Corollary 3.2.

For any finitely supported sequence $b_{n}$, $n=1,2,\ldots,$ we have the bound

(3.5) $\left\|\sum_{n:\,\rho(n)=d}b_{n}t_{n}(x)\right\|_{p}\leq
2^{d-1}(p-1)^{d/2}\left\|\sum_{n:\,\rho(n)=d}b_{n}t_{n}(x)\right\|_{2},\quad
p>2.$

Note that inequality (3.5) with another constant $c(d)$ instead of $2^{d-1}$
was proved in [LoRo] (see Theorem 5.13). The method of proof in [LoRo] is the
Riesz products technique, quite a different from ours. An examination of its
proof may provide a value $c(d)\geq 16^{d}$ for the mentioned constant. We
have a question concerning the possible optimal constant in (3.2).

######  Question 1.

Letting $\Phi(t)=t^{p}$, $1<p<\infty$, prove inequality (3.2) with the right
side constant independent of $d$.

One can check that in the case

(3.6) $\int_{0}^{\infty}\frac{dt}{\Phi(t)}=\infty,$

the statement of this problem doesn’t hold, since for the trigonometric and
Walsh Dirichlet polynomials

(3.7) $D_{2^{n}}^{T}(x)=\sum_{k=1}^{2^{n}}\cos 2\pi kx,\quad
D_{2^{n}}^{W}(x)=\sum_{k=1}^{2^{n}}w_{k}(x)$

we have

$\frac{\|D_{2^{n}}^{T}\|_{\Phi}}{\|D_{2^{n}}^{W}\|_{\Phi}}\to\infty\text{ as
}n\to\infty.$

Polynomials (3.3) tell nothing us about the optimal constant in (3.3). So in
this case we suggest the following.

######  Question 2.

Prove inequality (3.3) with the right side constant independent of $d$ for any
Young function $\Phi(t)$ .

## References