  New neutrino interactions beyond the Standard Model (BSM) have been of much
interest in not only particle physics but also cosmology and astroparticle
physics. We numerically investigate the time delay distribution of
astrophysical neutrinos that interact with the cosmic neutrino background.
Using the Monte Carlo method, we develop a framework that enables us to
simulate the time-dependent energy spectra of high-energy neutrinos that
experience even multiple scatterings en route and to handle the sharp increase
in the cross section at the resonance energy. As an example, we focus on the
case of secret neutrino interactions with a scalar mediator. While we find the
excellent agreement between analytical and simulation results for small optical
depths, our simulations enable us to study even optically-thick cases that are
not described by the simplest analytic estimates. Our simulations are used to
understand effects of cosmological redshifts, neutrino spectra and flavors. The
developments will be useful for probing BSM neutrino interactions with not only
current neutrino detectors such as IceCube and Super-Kamiokande but also future
neutrino detectors such as IceCube-Gen2 and Hyper-Kamiokande.
