Unraveling the Unexpected Seasonal Variation of Multi-Muon Events at NO$\nu$A

The seasonal variation of single muons is a well-understood phenomenon, mainly driven by a positive correlation with temperature fluctuations in the atmospheric profiles. However, the rate of multi-muon events recorded by various experiments has shown an intriguing opposite seasonal modulation that remains unexplained by previous research. For the first time, we quantitatively describe the phase and amplitudes of the seasonal variation for cosmic multi-muon events detected by the NO$\nu$A Near Detector. We can further explain the amplitude dependence for multi-muon events across various multiplicities. For our analysis, we use the general-purpose Monte Carlo code FLUKA-CERN 4.4.1, which provides a more realistic description of the detector, atmospheric profiles, and muon propagation underground. Finally, we compare our results with those obtained from the latest CORSIKA version 7.8010, utilizing the most up-to-date high-energy hadronic interaction models EPOS LHC-R and QGSJETIII-01. Our findings provide a fresh perspective on seasonal muon flux modulation and offer key constraints for cosmic-ray interaction models and underground detector studies.