CERN Accelerating science

The modification of the No-Time-Counter scheme, used to simulate the interactions within the system's cells within the Direct Simulation Monte Carlo approach, for describing interactions in the MeV primordial plasma. First, we sample $N_{\text{sampled}}$ pairs to interact, Eq.~\eqref{eq:Nsampled}. For each pair, we compute its interaction weight and make an intermediate decision on whether it will interact using the criterion~\eqref{eq:interaction-acceptance}. Then, we sample the kinematics of the interacting particles, generate the final states resulting from the collision, and make the final decision of whether the interaction takes place from the Pauli principle~\eqref{eq:acceptance-pauli}. Finally, we update the local properties of the plasma: the EM plasma temperature and the number of EM particles, as well as neutrino flavor distributions by the oscillation probabilities, Eq.~\eqref{eq:oscillation-probabilities}.

The evolution of neutrinos and the EM plasma energy densities ratio under the scenario where the neutrino distribution shape is thermal (Eq.~\eqref{eq:distribution-thermal}), but has temperature $T_{\nu}$ different from the EM plasma $T_{\text{EM}}$. For the initial setup, we consider $T_{\nu}= 3.5\text{ MeV}$ and $T_{\text{EM}} = 3\text{ MeV}$. \textit{Top panel}: The energy densities ratio $\delta \rho_{\nu}$, given by Eq.~\eqref{eq:delta-rho}. The blue line shows the result of our DSMC approach, the green line denotes the prediction of the modified code from~\cite{Akita:2020szl,Akita:2024ork}, whereas the orange line is obtained using the method of integrated neutrino Boltzmann equations from~\cite{EscuderoAbenza:2020cmq}, which assumes that the shape of the neutrino distribution is perfectly thermal throughout the whole evolution. \textit{Bottom panel}: The analog of $\delta \rho_{\nu}$ but for the number densities of neutrinos and the EM particles, highlighting the deviation from the thermality of the neutrino spectrum throughout the evolution.

The behavior of the ratio~\eqref{eq:delta-rho} under the injection of 20 MeV neutrinos equally to all neutrino flavors at the temperature $T_{\text{EM}} = 3\text{ MeV}$. The total injected energy density is $\rho_{\nu,\text{inj}}/\rho_{\nu,\text{total}} = 5\%$. The blue line shows the prediction of the DSMC method, whereas the orange one corresponds to the integrated Boltzmann approach from~\cite{EscuderoAbenza:2020cmq}.

The temporal evolution of the plasma after the injection of neutrinos with energies $E_{\nu} = 70\text{ MeV}$ and the overall energy density $\rho_{\nu,\text{inj}}/\rho_{\nu,\text{total}} = 30\%$. The other parameters of the setup are similar to the one considered in Fig.~\ref{fig:injection-20-MeV}. \textit{Top panel}: the behavior of $\delta \rho_{\nu}$ with temperature, where we show the predictions of the DSMC (the blue curve) and the integrated approach from Ref.~\cite{EscuderoAbenza:2020cmq} (the orange curve). \textit{Bottom panel}: comparison of the shape of the neutrino energy distribution for the system from Fig.~\ref{fig:injection-70-MeV-30-perc} at the moment when $\delta \rho_{\nu} = 0$ during the equilibration, as obtained with the DSMC simulation (the blue curve) and assuming the equilibrium neutrino spectrum (the red curve).

The temporal evolution of the plasma after the injection of neutrinos with energies $E_{\nu} = 70\text{ MeV}$ and the overall energy density $\rho_{\nu,\text{inj}}/\rho_{\nu,\text{total}} = 30\%$. The other parameters of the setup are similar to the one considered in Fig.~\ref{fig:injection-20-MeV}. \textit{Top panel}: the behavior of $\delta \rho_{\nu}$ with temperature, where we show the predictions of the DSMC (the blue curve) and the integrated approach from Ref.~\cite{EscuderoAbenza:2020cmq} (the orange curve). \textit{Bottom panel}: comparison of the shape of the neutrino energy distribution for the system from Fig.~\ref{fig:injection-70-MeV-30-perc} at the moment when $\delta \rho_{\nu} = 0$ during the equilibration, as obtained with the DSMC simulation (the blue curve) and assuming the equilibrium neutrino spectrum (the red curve).

Comparison of the DSMC approach with the discretization codes for the setup of injection of 70 MeV neutrinos at $T_{\text{EM}} = 3\text{ MeV}$. Two configurations are considered: equal injection among the flavors (the top panels) and the injection solely into $\nu_{e}$ (the bottom panels). In both cases, the injected energy fraction is $\rho_{\nu,\text{inj}}/\rho_{\nu,\text{total}} = 5\%$. The left plots show the evolution of $\delta \rho_{\nu}$, given by Eq.~\eqref{eq:delta-rho}. In the plots, the blue lines are the DSMC predictions, the green lines denote the calculation by the discretization approach from~\cite{Akita:2020szl} (see also~\cite{Akita:2024ta}), the dashed blue line is the result obtained in~\cite{Boyarsky:2021yoh} (see text for discussions). The right plots are snapshots of the electron neutrino distribution spectrum at the temperature when $\delta \rho_{\nu} = 0$. In addition to the results from the DSMC and Ref.~\cite{Akita:2024ork}, we include the plot of the equilibrium neutrino distribution given by $E_{\nu}^{2}f_{\text{FD}}(E_{\nu},T_{\text{EM}})$ (solid red lines).

The same setup as in Fig.~\ref{fig:injection-20-MeV} but under an injection of 500 MeV neutrinos.

Impact of injection of heavy SM particles in the primordial plasma. \textit{Top panel}: the distribution of electron and muon neutrinos produced by decays of $K_{L}K_{S}$ pairs. When simulating their decay, we used the module of \texttt{SensCalc} tool~\cite{Ovchynnikov:2023cry}. For the chain of the decay products, we account for instant kinetic energy loss by charged particles. The continuous extension of the spectrum to $\simeq 200\text{ MeV}$ is caused by the direct decay of kaons into neutrinos. The increase at $E_{\nu} = 50\text{ MeV}$ follows from decays of secondary muons stopped in the plasma, whereas the sharp increase at $E_{\nu}\approx 34\text{ MeV}$ originates from decays of secondary pions. \textit{Bottom panel}: The evolution of the quantity $\delta \rho_{\nu}$ under the injection of $\mu^{+}\mu^{-}$ (the blue curve) and $K_{L}K_{S}$ (the green curve) in the primordial plasma at temperature $T_{\text{EM}} = 3\text{ MeV}$. The curves start at different temperatures $T_{\text{EM}}\neq 3\text{ MeV}$ because EM decays of these particles reheat the EM plasma.

The temporal evolution of the quantity $\delta \rho_{\nu}$ when varying numbers of neutrinos per cell $N_{\text{cell},\nu}$ and particles in the system $N$. \textit{Left panel}: fixing $N_{\text{cell},\nu} = 400$ and varying $N$. \textit{Right panel}: fixing $N = 3\cdot 10^{6}$ and varying $N_{\text{cell},\nu}$.

Evolution of the neutrino distribution function under the initial setup with equilibrium neutrinos and EM plasma at temperature $T = 3\text{ MeV}$ and non-equilibrium neutrinos with energies uniformly distributed in the range $300 \text{MeV} < E_{\nu} < 450 \text{MeV}$. Their total energy density is related to the total energy of equilibrium part as $\rho_{\nu_\alpha}^{\text{non-eq}}/\rho_{\nu_\alpha}^{\text{eq}} = 0.15$. The non-equilibrium part of the spectra rapidly loses its energy in the first steps of simulation, leading to the distortions of the spectra at high energies which are eventually equilibrated. \textit{Left plot}: snapshots of the binned neutrino distribution function as obtained at different iterations of the DSMC simulation. The iteration 0 corresponds to the initial setup, while the iteration 200 is the final state. \textit{Right plot}: the comparison of the neutrino distribution function between the initial (the blue line) and final (the green line) iterations. The final distribution approaches the analytic Fermi-Dirac neutrino distribution with the temperature equal to the temperature of the electromagnetic plasma $T_{\text{EM}}$, which we shown by the dashed green line.

The evolution of the neutrino distribution $dn_{\nu}/dE_{\nu}$ averaged over all flavors under the assumption of fully equilibrium initial conditions~\eqref{eq:rho-distribution-thermal} and~\eqref{eq:rho-ratio-thermal}. The ``Iter \#'' curves correspond to the number of the iteration. No significant changes are developed throughout the simulation. The minor changes are related to the quality of the sampler of the kinematics of the electrons via the Fermi-Dirac distribution. The dashed green line shows the analytic Fermi-Dirac distribution with the temperature equal to the temperature of the electromagnetic plasma $T_{\text{EM}}$.

Evolution of the neutrino distribution function $dn_{\nu}/dE_{\nu}$ averaged over all flavors under different initial setups, showing how DSMC drives it towards thermal equilibrium with the EM sector. \textit{Left panel}: with equilibrium neutrinos and EM plasma at temperature $T_{\text{EM}} = 3\text{ MeV}$ and non-equilibrium neutrinos with energies randomly distributed in the range $300 \text{ MeV} < E_{\nu} < 450 \text{ MeV}$. Their total energy density is related to the total energy of the equilibrium part as $\rho_{\nu_\alpha}^{\text{non-eq}}/\rho_{\nu}^{\text{eq}} = 0.15$. The non-equilibrium part of the spectra rapidly loses its energy in the first steps of simulation, leading to the distortions of the spectra at high energies, which are eventually equilibrated. The plot shows the snapshots of the binned neutrino distribution function as obtained at different iterations of the DSMC simulation. The iteration 0 corresponds to the initial setup, while the iteration 400 is the final state. For comparison, the long-dashed green line shows the Fermi-Dirac distribution $dn_{\nu}/dE_{\nu} = E_{\nu}^{2}f_{\text{FD}}(E_{\nu},T_{\text{EM,final}})$, being the thermal equilibrium of neutrinos with the EM plasma with the final temperature $T_{\text{EM,final}} \approx 3.15\text{ MeV}$. \textit{Right panel}: with equilibrium neutrinos having temperature $T_{\nu_{\alpha}} = 3.5\text{ MeV}$ and EM plasma at temperature $T_{\text{EM}} = 3\text{ MeV}$. The meaning of the lines is the same, while the number of iterations is 250.

The evolution of the ratio of the neutrino energy density to the EM energy density in DSMC simulation compared to the theoretical prediction from~\cite{EscuderoAbenza:2020cmq}, under an assumption that the shape of the neutrino distribution function is always thermal at each step of the simulation. The initial conditions for the setup are $T_{\nu_i} = 3.2$ MeV for every flavor and the temperature of the EM plasma is $T_{\text{EM}} = 3$ MeV. \textit{Left panel}: not including the expansion of the Universe. Due to the absence of expansion, the ratio approaches to the exact SM value. \textit{Right panel}: expansion included.