The other day, I received a very interesting e-mail from Don Page, who saw one of my posts on Quora concerning neutrino-antineutrino annihilations. Don was wondering about the frequency of such events occurring in the relic cosmic neutrino background. More specifically, Don was wondering how often it may have happened that a neutrino-antineutrino pair produced more than one neutrino-antineutrino pairs?
Not sure if this was his intent, but Don sure pulled me down a fascinating rabbit hole. (No, I am not complaining. Such rabbit holes are fun.) In the end, we exchanged several e-mails, and I also consulted with our AI friends, notably the language models GPT4-o1 and Claude 3.5. As a result, I have been able to construct I think a reasonably coherent estimate.
I began with a few important assumptions:
- If we set the scale factor at $a=1$ at the present epoch (i.e., $a=1/(1+z)$ with the redshift $z=0$ at present), it was $a\approx 1/1091$ at the time of recombination;
- The present-day neutrino density (all species) is $n_\nu\approx 336~{\rm cm}{}^{-3}$ and it scales as $a^{-3}$;
- The neutrino background temperature is about $K_\nu\approx 1.95~{\rm K}$, hence the neutrino energy is $E_\nu\approx k_B T_\nu$, with $k_B= 8.6173303\times 10^{-5}~{\rm eV}/{\rm K}$, and with $T_\nu$ scaling as $a^{-1}$;
- We have (very roughly) $\Omega_{\Lambda}=0.7$, $\Omega_m=0.3$, with negligible $\Omega_\gamma$ and $\Omega_k$;
- The Hubble-parameter is about $H_0=70~{\rm km}/{\rm s}/{\rm Mpc}$.
Next, we have two expressions of importance, relating time-to-the-present $t$ and comoving distance $d$ to redshift $z$:
\begin{align}
t(z)&{}=\frac{1}{H_0}\int_0^{1/(1+z)}\frac{dx}{x\sqrt{\Omega_\Lambda+\Omega_kx^{-2}+\Omega_mx^{-3}+\Omega_\gamma x^{-4}}},\tag{1}\\
d(z)&{}=\frac{c}{H_0}\int_0^z\frac{dx}{\sqrt{\Omega_\Lambda+\Omega_k(1+x)^2+\Omega_m(1+x)^3+\Omega_\gamma(1+x)^4}}.\tag{2}
\end{align}
From the first of these expressions, we get
\begin{align}
dt &{}= \frac{dt}{d(1+z)^{-1}}\frac{d(1+z)^{-1}}{dz}dz \nonumber\\
&{}= -\frac{1}{H_0}\frac{(1+z)^{-1}dz}{\sqrt{\Omega_\Lambda+\Omega_k(1+z)^2+\Omega_m(1+z)^3+\Omega_\gamma(1+z)^4}},\tag{3}
\end{align}
which will allow us to change the integration variable from $t$ to $z$ when convenient.
The light cone volume at redshift $z$ is
\begin{align}
V=\tfrac{4}{3}\pi a^3d(z)^3=\tfrac{4}{3}\pi (1+z)^{-3}d(z)^3.\tag{4}
\end{align}
The $\nu\bar{\nu}$ cross-section is
\begin{align}
\sigma_\nu=G_F^2E_\nu^2,\tag{5}
\end{align}
with $G_F\approx 1.1663787\times 10^{-5}~{\rm GeV}^{-2}$ being the Fermi constant. As $E_\nu$ scales as $1+z$, $\sigma_\nu$ scales as $(1+z)^2$.
The corresponding event rate, assuming that the neutrinos are relativistic (moving at $\sim c$), is given by
\begin{align}
R=(1+z)^6n_\nu^2(1+z)^2\sigma_\nu c.\tag{6}
\end{align}
The ratio of events that may yield multiple $\nu\bar{\nu}$ pairs at redshift $z$ is
\begin{align}
f=G_F^2(1+z)^4E_\nu^4.\tag{7}
\end{align}
Putting it all together, the number of multiple pair events from the time of recombination $t_r$ corresponding to $z_0\approx 1090$ to the present $t=t_0$ at $z=0$ is
\begin{align}
N&{}=\int_{t_r}^{t_0}dt\,fRV\nonumber\\
&{}=\tfrac{4}{3}\pi n_\nu^2\,G_F^4E_\nu^6c^4H_0^{-4}\int_0^{z_0}\frac{(1+z)^{8}\displaystyle\left[\int_0^z\frac{dx}{\sqrt{\Omega_\Lambda+\Omega_k(1+x)^2+\Omega_m(1+x)^3+\Omega_\gamma(1+x)^4}}\right]^3\, dz}{\sqrt{\Omega_\Lambda+\Omega_k(1+z)^2+\Omega_m(1+z)^3+\Omega_\gamma(1+z)^4}}.\tag{8}
\end{align}
If we set $\Omega_k=\Omega_\gamma=0$ (good enough for our approximation), we get
\begin{align}
N&{}=\tfrac{4}{3}\pi n_\nu^2\,G_F^4E_\nu^6c^4H_0^{-4}\int_0^{z_0}\frac{(1+z)^{8}\displaystyle\left[\int_0^z\frac{dx}{\sqrt{\Omega_\Lambda+\Omega_m(1+x)^3}}\right]^3\, dz}{\sqrt{\Omega_\Lambda+\Omega_m(1+z)^3}}\\
&{}\approx 2.353\times 10^{-6}\int_0^{1090}\frac{(1+z)^{8}\displaystyle\left[\int_0^z\frac{dx}{\sqrt{\Omega_\Lambda+\Omega_m(1+x)^3}}\right]^3\, dz}{\sqrt{\Omega_\Lambda+\Omega_m(1+z)^3}}.\tag{9}
\end{align}
Using $\Omega_\Lambda=0.7$, $\Omega_m=0.3$, we get
\begin{align}
N\approx 1.1\times 10^{18}\tag{10}
\end{align}
multiple-pair events. However, if we change the upper integration limit from 1090 to 1, we get
\begin{align}
N_1\approx 0.000024,\tag{11}
\end{align}
so essentially no such events happened within our past light cone since $z=1$.
* * *
The calculation up to this point relied on several simplifying assumptions concerning the nature of neutrinos:
- We assumed that the neutrinos are ultrarelativistic, allowing us to use $v=c$ in the expression for $R$. If $E_\nu \lesssim m_\nu$, this assumption is not valid anymore and we need to use an appropriate value for $v$, namely $v=c\sqrt{1-[(m_\nu c^2)/(1+z)E_\nu]^2}$. This is further complicated by the fact that neutrino masses are presented in the form of a mass mixing matrix;
- The quantity $n_\nu$ is the total neutrino density. We assume all three species are present in equal numbers and that furthermore, they can annihilate, hence the use of $n_\nu^2$. Perhaps an expression like $n_i {\cal M}^{ij} \bar{n}_j$ would be more appropriate, where $n_i$ is the number density of the $i$-th neutrino species, $\bar{n}_i$ is the antineutrino number density, and ${\cal M}^{ij}=m^{ij}/|{\rm det}\,m|^{1/3}$ is the neutrino mass mixing matrix normalized to unit determinant;
- The production of two new pairs of neutrinos will be strongly suppressed if the collision center-of-mass energy becomes similar to, or less than, the rest mass-energy of the new pairs produced. This suppression could be implemented by scaling $f$, using a term such as $1/[\exp(4-E_\nu/m_\nu c^2)+1]$, though again things are made more complicated by the presence of the mass mixing matrix.
This yields a revised estimate in the form,
\begin{align}\hskip -0.4in
N=\tfrac{4}{3}\pi G_F^4E_\nu^6c^4H_0^{-4}\sum\limits_{i,j}n_i\frac{m^{ij}}{|{\rm det}\,m|^{1/3}}\bar{n}_j
\int_0^{z_0}
\frac{\sqrt{1-\left[\dfrac{m^{ij} c^2}{(1+z)E_\nu}\right]^2}}{e^{4-(1+z)E_\nu/m^{ij} c^2}+1}
\frac{(1+z)^{8}\displaystyle\left[\int_0^z\frac{dx}{\sqrt{\Omega_\Lambda+\Omega_m(1+x)^3}}\right]^3\, dz}{\sqrt{\Omega_\Lambda+\Omega_m(1+z)^3}}.\tag{12}
\end{align}
Then again, this part is a bit academic: we do not know the actual neutrino masses, and in any case, the overwhelming majority of events appear to have happened when the neutrino background was much hotter than today, so probably relativistic. The one refinement that I see that might be relevant is that the previous results (9)–(11) must be multiplied by $\tfrac{1}{3}$: that is because neutrino species do not mingle freely, and in fact if we were to use the identity matrix in place of the neutrino mass mixing matrix, this is the result that we would see.
* * *
This is the Maxima script I used to evaluate the numerical prefactor for $N$ and verifying its dimensions:
eV:1e-9*GeV$
J:1/(1.602176634e-19)*eV$
K:8.6173303e-5*eV$
hbar: 1.054571817e-34*J*s$
c: 299792458.0*m/s$
GF:1.1663787e-5*GeV^(-2)$
En:1.95*K$
sg:GF^2*En^2$
sg,GeV=1/ev(hbar*c/GeV,m=100*cm)$
H0:70*km/s/Mpc,km=1000*100*cm,Mpc=3.09e16*1e6*100*cm$
nn:336/cm^3$
4/3*float(%pi)*nn^2*GF^4*En^6*c^4/H0^4,GeV=1/ev(hbar*c/GeV),m=100*cm,eval;
* * *
This is the Maple script I used to evaluate the integral numerically:
N:=Int((1+z)^8*(Int(1/sqrt(OL+Om*(1+x)^3),x=0..z))^3/sqrt(OL+Om*(1+z)^3),z=0..z0):
2.353e-6*evalf(eval(N,[OL=0.7,Om=0.3,z0=1090]));
2.353e-6*evalf(eval(N,[OL=0.7,Om=0.3,z0=1]));
* * *
Finally, this Maxima script can be used to evaluate the normalized neutrino mass mixing matrix, with values from Wikipedia:
m1:matrix([1,0,0],[0,cos(t23),sin(t23)],[0,-sin(t23),cos(t23)]);
m2:matrix([cos(t13),0,sin(t13)*exp(-%i*tCP)],[0,1,0],[-sin(t13)*exp(%i*tCP),0,cos(t13)]);
m3:matrix([cos(t12),sin(t12),0],[-sin(t12),cos(t12),0],[0,0,1]);
trigsimp(determinant(m1.m2.m3));
expand(m1.m2.m3),t12=33.41*float(%pi)/180,t23=49.1*float(%pi)/180,
t13=8.54*float(%pi)/180,tCP=197*float(%pi)/180;
expand(determinant(m1.m2.m3)),t12=33.41*float(%pi)/180,t23=49.1*float(%pi)/180,
t13=8.54*float(%pi)/180,tCP=197*float(%pi)/180;
expand(m1.m2.m3),t12=33.41*float(%pi)/180,t23=49.1*float(%pi)/180,t13=8.54*float(%pi)/180,tCP=%pi;
