Sibyll^(bigstar)

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.