  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$^{\bigstar}$, 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.
