For masses around 50 meV then, this gives a lifetime of $3.5 ×10^8$ seconds or around 11 years. This is more than enough time for a relic neutrino to decay since the universe was 1 second old.\\

In general, the model-independent contraint on neutrino decays is [42]

\begin{equation}
\frac{\tau_\nu}{m_\nu} \geq  10^5 \  \mathrm{s\cdot eV}^{-1}.
\end{equation}

If neutrinos did decay into photons, or into some exotic matter that we can't detect, then the number of C$ν$B neutrinos remaining to be detected would be

\begin{equation}
\end{equation}

where $N_0$ is the number after decoupling. This $λ_ν$ is given by [43]

\begin{equation}
\lambda_\nu = \int \frac{dt}{\tau_\nu} = \int_{0}^{z_{dec}} \frac{dz}{(1+z)H(z)\gamma(z)\tau_\nu^0},
\end{equation}

where $τ_ν^0$ is the proper lifetime of neutrinos, and $H(z)$ is the Hubble constant at some redshift $z$. This integral may be difficult to perform, but if we make the approximation that most neutrinos have decayed at later times (as $z$ goes to 0), we can insert a Dirac delta function and obtain

\begin{equation}
\lambda_\nu\approx \frac{1}{H_0\tau_\nu^0} = \frac{t_0}{\tau_\nu^0}.
\end{equation}

So, the approximate number of relic neutrinos left today would be given by

\begin{equation}
\end{equation}

If $τ_ν^0$ is much smaller than the age of the universe, then we would not be able to detect the C$ν$B at all - a very lamentable thought. If, as we hope, relic neutrinos are observed, this will place an extremely strong constraint on neutrino decays, as their lifetime would have to be at least larger than the age of the universe.\\

If neutrinos all decayed into the lightest species, this could either amplify or diminish the detection rate. If neutrinos are normally ordered, then all neutrinos would decay into $ν_1$ neutrinos, and instead of summing over all species, we would simply have thrice the amount of $ν_1$ neutrinos. Similarly for the inverted hierarchy, except with $ν_3$ neutrinos. As we will see in the capture rate chapter, the ordering will result in either doubling the capture rate, or diminishing it almost completely, if neutrinos do indeed decay into other neutrinos.\\

The final consideration made in this subsection is what the consequences of a non-standard thermal history would be. Recall from \eqref{n0} that the number density is dependent on the effective neutrino temperature cubed - a very sensitive quantity. If some ``dark radiation" (some species of matter that behaved as radiation but which we have not yet detected) was present, and it too annihilated between the decoupling of neutrinos and photons (similarly to the electrons and positrons), then the ratio calculated in \eqref{Trelation} would be altered. Owing to the cubic nature of the relationship, a doubling of the neutrino effective temperature would result in the number density today being 8 times larger. Similarly, a colder temperature would diminish the number density on earth. This method of entropy injection from the annihilation of some species would actually result in a colder C$ν$B, and thus a smaller number density.

Recall that

\begin{equation}
\frac{T_{\gamma,0}}{T_{\nu,0}}\propto\left(\frac{g_*(T_{dec})}{g_*(T_0)}\right)^{1/3}=\left(\frac{11/2+\Delta g_*}{2}\right)^{1/3},
\end{equation}

so any dark radiation would increase this ratio. Since $T_γ,0$, the present-day temperature of the CMB, is a known observable, this would mean the neutrino effective temperature would decrease, and in turn so would the number density. This $Δg_*$ is quite well constrained, as $N_eff$ would of course be affected if there was any other radiative species present before recombination. \\

This chapter has aided us in achieving a more accurate prediction of what the cosmic neutrino background may look like around earth today, and which effects would be most noticeable. While annual modulation is one of the most natural characteristics of a C$ν$B, its effect is too small to be measurable in the near future. Distortions to the number density from non-instantaneous decoupling and gravitational clustering are competing effects, with the dominant one decided by the mass of the neutrino: with neutrino masses above approximately 10 meV, clustering effects overshadow the spectral distortion, while relativistic neutrinos would only exhibit the effects of spectral distortion. Other non-standard physics may explain any deviations from these predictions, but as we shall discuss later, it may be difficult to distinguish between models.

\newpage
\section{Types of Neutrino Detectors}

For the purposes of this text, supernovae are the main astrophysical source we will consider. The reason for this is that supernovae are short events, and thus data from both electromagnetic observations and neutrinos (and perhaps in future gravitational waves) can be compared to gain insights. Though much remains to be learned about supernovae, and data extracted is often model-dependent, the addition of precision neutrino astronomy will assist in making great strides in the field.\\
The main source of information regarding these events is SN1987A, a supernova that occurred in the Large Magellanic Cloud near the Milky Way [44]. From this, we learned that supernova neutrinos have energies in the range of tens of MeV - many orders of magnitude larger than the relic neutrinos, as seen in \eqref{lowenergy}.\\
Owing to this energy difference, the detectors used for supernova neutrinos and those for the C$ν$B utilise different capture mechanisms. This chapter does not pretend to act as a review of all detection methods, but simply attempts to elucidate some of the most popular methods used.

\subsection{Astrophysical Neutrino Detectors}

Owing to the small masses seen in \eqref{massupperlim}, supernova neutrinos are extremely relativistic, with speeds close to that of light. Their high energy allows for certain interactions that relic neutrinos are not energetic enough to achieve.\\

\underline{Cherenkov Radiation}\\

Cherenkov radiation is a phenomenon wherein charged particles travelling faster than light in a medium (i.e $v>1/n$ when $c=1$) emit radiation in a characteristic conal shape, analogous to the sonic boom achieved by objects moving faster than the speed of sound in a medium. All decelerating charged particles emit radiation, but the conal structure is unique to Cherenkov radiation. One way that Cherenkov radiation can be used to detect neutrinos is if a neutrino transfers energy to a charged particle in a detector.\\
Consider some charged particle moving with a speed $v$, emitting radiation at some angle $θ$ from its axis of propagation, symmetrically around the axis.

\begin{figure}[!h]
\centering
\includegraphics[width=0.3\linewidth]{Cherenkov}
\caption{As the particle travels $vt$, the emitted light travels $\frac{t}{n}$. Drawn using [6].}
\end{figure}

In some time $t$ (short enough such that its velocity does not change considerably in this time),the particle would have moved $vt$ while the radiation would have moved $t/n$. Thus, we find

\begin{equation}
\cos\theta=\frac{t/n}{vt}=\frac{1}{nv},
\end{equation}

and since $nv>1$ for all times of emission, this angle is always acute, and thus as the particle travels we see the characteristic cone of light. For example, for a very relativistic particle $v≊1$ in water with $n≊1.33$, the cone will make an angle of $≊41.25^∘$.\\
In order for this to be energetically possible, we need the kinetic energy, given by

\begin{equation}
\end{equation}

to obey $v>1/n$. So the condition for Cherenkov radiation is that the charged particle have energy larger than

\begin{equation}
\end{equation}

As an illustrative example, for an electron with mass 0.511\ MeV in water with refractive index 1.333, we find $E_K^thr=264.06 keV$.\\
Remember that this is the energy that the electron (or any other charged particle) must have, not the neutrino we wish to detect. Neutrinos can scatter off electrons via

\begin{equation}
\nu_\alpha + e^- \ \ \rightarrow \ \ \nu_\alpha + e^-.
\end{equation}

This technique is extremely useful as we gain information not only about the neutrino's energy, but also about its direction of incidence as charged particles (in this case, the electrons) paths can be tracked and the Cherenkov cones seen. The interaction cross section can be calculated using the electroweak theory above. For electron neutrinos, there will be both charged and neutral current interactions, but for muon and tau neutrinos there will of course only be the neutral current. Thus, this method favours the detection of electron neutrinos. Though the threshhold energy does not translate to exactly the required energy of the neutrino, in elastic scattering with electrons, it is pretty close. If we consider an initially stationary electron hit by a neutrino which afterwards comes to rest, we would have

\begin{equation}
\end{equation}

or, even more simply, $E_K,e=E_K,ν$. Thus, in this limiting case, the kinetic energy of the detected electron is the exact same as that of the incoming neutrino. Of course, if the neutrino does not come to rest, or if the electron has some initial velocity, there can be slight deviations. However, we can take the threshhold energy for Cherenkov radiation as a good estimate for the kinetic energy needed in the incident neutrino. Thus, for even the lightest charged particle (the electron) in any realistic material, relic neutrinos would not be able to trigger Cherenkov radiation, while supernova neutrinos can.\\
The other way that a neutrino can create a Cherenkov cone is via weak interactions with nuclei in the medium,

\begin{equation}
\nu+X(A,Z) \ \ \rightarrow \ \ Y(A,Z+1) + e^-,
\end{equation}

and similarly for the production of positrons by the incidence of an anti-neutrino. In this case, the incident (anti-)neutrino must have enough energy to overcome the binding energy of the $X$ nucleus as well as provide enough kinetic energy to the (positron) electron in order for it to be above $E_K^thr$. Once again, taking water as an example, the anti-neutrino could interact with the hydrogen nucleus (proton) or the neutrino with the oxygen nucleus, turning them into a neutron or fluorine nucleus, respectively. The former needs $≈1.29$\ MeV while the latter needs around 16\ MeV. For supernova neutrinos these values are certainly attainable, even if the latter is only for some, but for cosmic neutrinos it simply is not energetically possible. Finally, note that only the ejected electron might create Cherenkov light, as the energy threshhold for nuclei, as seen in \eqref{thresh} is around $m_nuc$, while our supernova neutrinos are of the order of tens of MeV.\\

Another variation of these detectors are ``Long-String" Cherenkov detectors, named for the long strings of photomultiplier tubes (PMTs) placed in the medium and used to detect the Cherenkov radiation. The advantage of this variation is the increased time sensitivity, owing to the array of PMTs, but the drawback is the loss of individual interaction information - only statistical information is gathered [45].\\

\underline{Scintillation Detectors}\\

As we have seen, anti-neutrinos can interact with nuclei via the weak interaction. One of these interactions is inverse beta decay,

\begin{equation}
\bar{\nu}_\alpha + p \ \rightarrow \ n + e^+.
\end{equation}

The energy threshold for this interaction is $m_n+m_e-m_p≊1.8  MeV$, and so while essentially impossible for cosmic neutrinos, this reaction is well within the capabilities of supernova neutrinos. Given a detector full of free electrons, the positrons then go on to undergo electron-positron annihilation,

\begin{equation}
e^- + e^+ \ \rightarrow \ 2\gamma_e
\end{equation}

and, after a short delay, the neutrons undergo neutron capture:

\begin{equation}
p + n \ \rightarrow \ d + \gamma_d.
\end{equation}

At the energies of supernova neutrinos, the electron-positron annihilation will only produce photons, and not other massive bosons. These $γ_e$ photons will have energies around $m_e≊0.511 MeV$, while the energy released by the fusion of $p$ and $n$ means that the $γ_d$ photons will have energies around $m_p+m_n-m_d≊2.22 MeV$.\\

Detectors that use the detection of these photons to observe neutrinos are known as scintillation detectors, and consist of large volumes of the interacting material, usually in liquid form, surrounded by photomultiplier tubes in order to enhance and capture the photonic signal [45]. The clear signal of a detection is the characteristic chain of events: a pair of $≊0.511 MeV$ photons, followed after a few hundred microseconds by the $≊2.22 MeV$ signal. This technique was the first to successfully detect neutrinos back in 1956.\\

Because the photons are emitted isotropically, they cannot directly give information about the direction of the incoming anti-neutrino. However, with very small segmentation, information from both $γ_e$ and $γ_d$ may give some idea of the incident direction [45].\\

There are many other variants of detectors, but they all function using similar methods: either produced photons are measured or other particles produced. For example, liquid argon time-projection chambers use the interaction

\begin{equation}
\nu_e +\ ^{40}\mathrm{Ar} \ \rightarrow \ ^{40}\mathrm{K}^* + e^- \ \rightarrow \ ^{40}\mathrm{K} + e^- + \gamma_K,
\end{equation}

and detect the $γ_K$ arising from the de-excitation of the krypton produced during the charged current interaction.  In the same experimental setup, anti-neutrinos can be detected using

\begin{equation}
\bar{\nu}_e +\ ^{40}\mathrm{Ar} \ \rightarrow \ ^{40}\mathrm{Cl}^* + e^- \ \rightarrow \ ^{40}\mathrm{Cl} + e^- + \gamma_{Cl},
\end{equation}

and the characteristic de-excitation photons detected once again. This same design is used in so-called heavy nuclei detectors, which use substances such as iron and lead.

\subsection{Mechanisms for Detecting Relic Neutrinos}

\underline{Elastic Scattering}\\

As we have seen, the factor of $G_F=1.1664×10^-5  GeV^-2$ which dictates weak interactions (in the 4-Fermi approximation, which is appropriate for the low energies at which we are working) is a small quantity: that is, terms including higher orders of $G_F$ are smaller. Thus, it seems to be in our best interest to analyse interactions that are of the order $G_F$.

The first mechanism we can consider is ``neutrino optics" - if the relic neutrinos are travelling through a medium with inter-atomic spacing smaller than their de Broglie wavelength, this medium can be characterised by a refractive index $n$, with $n=p'/p$, $p$ being the neutrino's original momentum and $p'$ its altered form. This $n$ is found to be $n=1+𝒪(G_F)$ [46], and so at first this seems promising. However, it was shown that the force of this refraction would be proportional to $∇̅n_ν$, and in the case of an experiment on earth, the neutrino density is negligibly different from constant, and so this force vanishes [47]. \\

The second option available at this order of $G_F$ is known as the Stodolsky effect. Consider a relic neutrino scattering off an electron. At low energy, we can use our 4-Fermi approximate Lagrangian, and from \eqref{4FermiLag} the Hamiltonian (density) is clearly just the negative of this. Only considering electrons and neutrinos, we have

\begin{equation}
\mathcal{H} = \frac{G_F}{\sqrt{2}}\bar{e}\gamma^\mu(1-\gamma^5)e\bar{\nu_e}\gamma^\mu(1-\gamma^5)\nu_e.
\end{equation}

In theory, there will be an induced energy difference $ΔE$ owing to the spin up and down electrons [46]. In the non-relativistic limit, the electron current yields a factor of $σ̅_̅e̅·v̅_̅e̅$, while the neutrino current is found to be proportional to $n_ν-n_ν̅$. Thus, not only is the interaction dependent on the alignment of the electrons' spins, but a lepton asymmetry in the neutrino sector is needed. Thus, for the standard model wherein the lepton asymmetry is negligible, this mode too seems unrealistic. \\

At the next order of $G_F^2$, perhaps our luck will be better. In this case, we consider elastic scattering off of some large target whose nuclei have a mass number $A$, following [48]. As we discussed in the section above on annual modulation, the earth can be thought of as moving through the sea of neutrinos that is the C$ν$B. If the earth was stationary relative to the C$ν$B, the average momentum exchange between relic neutrinos and our target would be 0, as they would all cancel out. The motion of earth induces a dipole distortion in the momentum of the cosmic neutrinos. If the relic neutrino loses all of its momentum in the interaction, then the momentum exchange in earth's frame is

\begin{equation}
\braket{p}_{\mathrm{rel}} = v_{\mathrm{earth}}p_\nu \ \ \ \ ; \ \ \ \ \braket{p}_{\mathrm{nonrel}} = v_{\mathrm{earth}}m_\nu,
\end{equation}

for massless and massive species, respectively. For the massive case, it is assumed that the neutrino that interacts with the target has clustered in the Milky Way.

This imparting of momentum onto the target would cause a small acceleration:

\begin{equation}
a=\frac{F}{m_T}=\frac{\braket{p}}{m_T\Delta t}.
\end{equation}

This $1/Δt$ is the collision rate. Let us consider a target of mass 1 gram, such that $m_T=1$ and $N_T=N^A/A$, with the numerator being Avogadro's constant. Since the collision rate is $1/Δt = ϕ_νN_T σ_νN$ with $ϕ_ν=n_νv_ν$ the neutrino flux, we have

\begin{equation}
a = \frac{N^Av_{\mathrm{earth}}E_\nu n_\nu v_\nu \sigma_{\nu N}}{A},
\end{equation}

with $E_ν=p_ν$ for the relativistic neutrinos and $m_ν$ for the non-relativistic case. The acceleration of this 1g mass target will be the sum of the effects from the relativistic and non-relativistic neutrino species, as well as from both neutrinos and anti-neutrinos. Recall also from \eqref{4Fermiapprox} that the cross section is given approximately by

\begin{equation}
\sigma_{\mathrm{rel}} \approx G_F^2 p_\nu^2 \ \ \ \ ; \ \ \ \ \sigma_{\mathrm{nonrel}} \approx G_F^2 m_\nu^2.
\end{equation}

At leading order in the instantaneous decoupling limit, the number density for each species (adding both neutrinos and anti-neutrinos) is $≈112cm^-3$.

For the sake of numerical calculation, let us consider the minimal inverted hierarchy case, wherein we have 1 relativistic species and 2 massive species with mass around 50 meV and velocity approximately 0.01, using \eqref{lowenergy}. Then we have

\begin{equation}
a = 2a_{\mathrm{nonrel}} + a_{\mathrm{rel}} = \frac{2N^AG_F^2v_{\mathrm{earth}}n_0}{A} \left[2\times 0.01\times m_\nu^3+\braket{p_\nu}^3\right],
\end{equation}

where $⟨p_ν|$⟩ is the average momentum, calculated in (<ref>). We will take $A$ to be around 56, to correspond with using iron as a target. Using $N^A=6.022\times10^{23}$ [49] and $v_{\mathrm{earth}}\approx400$ km.s$^{-1} = 1.33\times10^{-3}$,[This is a maximum, using the fact that the sun is moving at around 369 km.s$^{-1}$ relative to the C$\nu$B [31] and the earth is moving with a velocity of around 30 km.s$^{-1}$ around the sun [32], $v_{\mathrm{earth}}$ can range between $\approx$ 330 and 400 km.s$^{-1}$.] this leads to a numerical acceleration of $a = 6.815 \times 10^{-52}$ cm.s$^{-2}$: our target will move roughly $1\times10^{-16}$ cm over the age of the universe. According to [48], a more precise calculation leads to $a \approx 10^{-27}f_c\rho_\nu$ cm.s$^{-2}$, with $f_c$ the amplification owing to clustering (around 1.12 for our 50 meV neutrinos) and $\rho_\nu$ the energy density of neutrinos in grams per cm$^3$.

Clearly, elastic scattering is not our best bet for detecting relic neutrinos.

High Energy Cosmic Rays

Cosmic “rays" are in fact particles such as protons and nuclei that come from distant sources and have extremely high energies. Owing to the presence of the CMB, there is an upper limit on the possible energies that these cosmic rays can have. As an illustrative example, let's consider very high energy protons.

As these protons are travelling towards us, they are moving through a sea of relic photons and there is a high chance they will undergo some kind of interaction. For example,

\begin{equation}
p \ + \ \gamma_{\mathrm{CMB}} \ \rightarrow \ p \ + \ \pi^0,
\end{equation}

where $\pi^0$ is the neutral pion. Let us consider 4-momentum conservation. Then we have

\begin{equation}
(p_i+p_\gamma)^2 = (p_f+p_\pi)^2 \ \ \rightarrow -m_p^2 + 2p_i\cdot p_\gamma = -(m_p+m_\pi)^2,
\end{equation}

where in the second part we are interested in the threshold, wherein the outgoing particles are approximately at rest. For simplicity, let us consider that the proton and CMB photon were moving towards each other along the z-axis. Then, since the proton is extremely relativistic, we can take their 4-momenta to be

\begin{equation}
p_i = (E_i,0,0,E_i) \ \ \ \ ; \ \ \ \ p_\gamma=(E_\gamma,0,0,-E_\gamma).
\end{equation}

Plugging this in, we have

\begin{equation}
-m_p^2 - 4E_iE_\gamma = -(m_p+m_\pi)^2 \ \ \rightarrow \ \ E_i=\frac{(m_p+m_\pi)^2-m_p^2}{4E_\gamma}.
\end{equation}

Using that the mass of the proton is 938.2721 MeV [50], the mass of the neutral pion is 134.9768 MeV [51], and the energy of the CMB photon is (on average) $E_\gamma=T_\gamma=T_{\mathrm{CMB}}=0.2348$ meV [27], we have the cutoff energy of the cosmic ray as around $E_i=2.89\times 10^{20}$ eV. More precise calculations give the GZK (Greisen-Zatsepin-Kuzmin) cutoff of $E_{\mathrm{GZK}}\approxeq 5 \times 10^{19}$ eV [52].

Experimentally, however, cosmic rays with energies higher than this cutoff (so called “Oh-My-God" particles) have been observed, which seems extremely unlikely. One explanation is that the cosmic ray must have come from a very nearby source, such that it has a higher chance of not yet having interacted with a relic photon. Unfortunately, there are no known nearby sources of such high energy rays. One explanation involves the C$\nu$B: of course, as we know, neutrinos only interact weakly, and thus would not interact with the CMB. So an extremely high energy neutrino, coming from a distant source, would be able to pass through the fog of relic photons without its energy being diminished. There is, of course, a chance that this cosmic ray neutrino interacts with a relic neutrino from the C$\nu$B. This interaction is a neutral current, and thus would create a Z-boson which could then decay into many possible combinations of particles. The energy threshold, calculating as before, would be

\begin{equation}
E_i = \frac{m_Z^2}{2m_\nu} = 4.158 \times 10^{21}  \left(\frac{1\mathrm{eV}}{m_\nu}\right) \mathrm{eV},
\end{equation}

where the relic neutrino was taken to be at rest, and the mass of the Z-boson much larger than the neutrino. Theoretically, this cosmic ray neutrino could interact with a relic neutrino very close to earth, emitting particles with energies higher than the GZK cutoff and allowing them to arrive on earth without interacting with a CMB photon, such that it appears as if they have come from a nearby source.

Another way in which the C$\nu$B might be seen is in the power spectrum of observed cosmic rays [53]. If a cosmic ray (say, for example, a proton) were to interact with a relic neutrino (continuing our example, via inverse beta decay), this should cause a change to the cosmic ray spectrum, depending on the mass of the neutrino. Performing a similar calculation to the GZK cutoff, we find that the threshold energy for a proton in this case would be around $E_i\approxeq 1.695\times 10^{15} \left(\frac{1\mathrm{eV}}{m_\nu}\right)$ eV. For example, a neutrino mass of $m_\nu = 5$ meV could explain the “knee" seen in the spectrum of energies at around $10^{17.5}$ eV [54] . Unfortunately, both of these methods are not only explicable by a relic neutrino population, and many other theories exist to explain both the kinks in the cosmic ray energy spectrum and the existence of rays with energies higher than the GZK cutoff [54]. Thus, perhaps this is not the most fool-proof method for detecting the C$\nu$B.

Capture by Beta-Decaying Nuclei

The most promising mode of detection seems to stem from the use of inverse beta decay, known as neutrino capture by beta-decaying nuclei (NCB), proposed by Weinberg [55]. Recall that our relic neutrinos' energies are extremely small, with maximal values around the neutrinos' masses. Thus, an important feature of NCB is that the nuclei used must readily undergo beta decay, such that there is no energy threshold on the reaction.

The difference between regular beta decay, which the nuclei will also go through, and beta decay owing to neutrino capture is a difference in the possible energies of the emitted electron. As we shall see, the gap between the beta decay endpoint (the maximum possible energy for an emitted electron during a regular beta decay event) and the energy of an electron as the result of neutrino capture is $E_{\mathrm{C\nu B}}+m_{\mathrm{lightest}}$, with $m_{\mathrm{lightest}}$ being the mass of the lightest neutrino species.

Since regular beta decay is the main source of background noise, it is preferable to choose an isotope for the experiment wherein the ratio between the C$\nu$B neutrino capture rate and the inverse of the half life of the isotope is maximised. Some candidates include tritium, $^3$H and Rhenium-187, $^{187}$Re [54]. Rhenium-187 is an interesting candidate, as its half life is massive - over 40 billion years [56]. This is useful as it will minimise background noise in the neutrino capture experiment. Tritium's half life, on the other hand, is much smaller at around 12.32 years [57]. On the other side of the coin, the cross section of neutrino capture on tritium (as we will calculate explicitly later) is found to be of the order $10^{-45}$ cm$^2$, while for rhenium-187 it is much lower at around $10^{-52}$ cm$^2$ [54]. So even though striking a balance is key, as we shall see in the next chapter, a cross section 7 orders smaller will make the capture rate - an already very small quantity - 7 orders smaller. For this reason, albeit creating a large amount of background noise for our experiment, it seems that NCB on tritium is our best bet for detecting the C$\nu$B, as such an interaction would not be explicable in any other way, and as we shall observe in the next chapter, the characteristic peaks in the emitted electron's spectrum will be a smoking gun indication.

For astrophysical neutrinos - particularly from supernovae - the detectors of choice seem to utilise the Cherenkov effect, as they are able to offer not only energetic measurements, but also track the angle of incidence, to some degree. While scintillation detectors are adept at the former, owing to the isotropic nature of the emission of de-excitation photons, the incident direction of the neutrino is more difficult to detect.

For relic neutrinos, as we saw, capturing them on tritium seems to be the best option for the foreseeable future. While interactions with cosmic ray ultra-high energy neutrinos may explain the existence of “Oh-My-God" particles, and distortions to the cosmic ray energy spectrum may prove the existence of the C$\nu$B, it is simply another indirect measurement - we are still not able to detect the energies, masses, and fermionic type of these elusive ancient particles. Elastic scattering - while indeed being a direct probe - is not a viable candidate owing to the miniscule reaction rates and the expected lack of lepton asymmetry. Should a lepton asymmetry be found in the universe, perhaps the Stodolsky effect will hold more promise. Other isotopes could indeed be used for capturing cosmic neutrinos, but factoring in half life and reaction rate among other variables, tritium seems to be the best bet.

§ THE CAPTURE RATE AND ENERGY SPECTRA

As we have seen in the previous chapter, the best chance we have of detecting the extremely elusive cosmic neutrinos lies in inverse beta decay. The PTOLEMY experiment, projected to create a detector made up of 100 grams of tritium ($^3\mathrm{H}$) [58], is the leading hope for the C$\nu$B's imminent detection. The reasons for the choice of tritium are that firstly, and most importantly, the isotope needs to readily undergo inverse beta decay without any energy threshold, as the relic neutrinos' energies are extremely small. Secondly, the isotope must have a fairly large half-life, so that the experimental setup does not decay too quickly - tritium's half life is $t_{1/2}\approx 12.32$ years. Finally, a low Q-value (the energy released by the reaction) is preferred. This is because the energy difference between this inverse beta decay and regular beta decay is very small (as we will see kinematically), and so working at smaller energies will make it easier to observe. The reaction utilised is given by

\begin{equation}
\nu_i\ +\ ^3\mathrm{H} \ \rightarrow \ ^3\mathrm{He}\ +\ e^- \label{truereaction}
\end{equation}

and in calculating the scattering amplitude for this reaction, we will be able to obtain the expected capture rate - that is, the rate at which we should detect these relic neutrinos.

§.§ Cross Section of Neutrino-Tritium Interaction

The full theory needed to calculate the scattering amplitude of the reaction in (<ref>) is the electroweak theory described in chapter 3. However, since this reaction is occurring at energies much smaller than the weak boson masses, we can safely utilise 4-Fermi theory, the low-energy effective field theory at tree level.

To begin, following the lead of [43], we consider the process

\begin{equation}
\nu_i + n \ \rightarrow \ p + e^- \label{reaction}.
\end{equation}

From our theory, we then get the matrix element

\begin{equation}
i\mathcal{M}_i = -i\frac{G_F}{\sqrt{2}}V_{ud}U^*_{ei}\left[\bar{u}_e\gamma^\alpha(1-\gamma^5)u_{\nu_i}\right]\left[\bar{u}_p\gamma_{\alpha}(f-g\gamma^5)u_n\right],
\end{equation}

where the factors of $f$ and $g$ have been added post-hoc to account for the fact that the neutron and proton are not elementary particles but have internal structure, consisting of quarks and gluons. These are known as the nucleonic form factors, with $f$ being for the proton and $g$ the neutron. These form factors are actually momentum dependent, but owing to the extremely small momenta of the relic neutrinos, we can take the values of $f$ and $g$ to be those in the limit as the transfer momentum goes to 0.

Squaring this, we have

\begin{equation}
\end{equation}

\begin{equation}
\mathcal{N}_1^{\alpha \beta} = {\rm tr}\left[\gamma^{\alpha}(1-\gamma^5)u_{\nu}\bar{u}_{\nu}\gamma^{\beta}(1-\gamma^5)u_e\bar{u}_e \right],
\end{equation}

\begin{equation}
\mathcal{N}_2^{\gamma \delta} = {\rm tr}\left[\gamma^{\gamma}(f-g\gamma^5)u_{n}\bar{u}_{n}\gamma^{\delta}(f-g\gamma^5)u_p\bar{u}_p \right].
\end{equation}

We need to sum over all the possible spins that the proton, neutron and electron may have (as well as averaging over the neutron spins), so we need to use the well-known completeness relations that are the outer products of the spinors:

\begin{equation}
\sum_{s_j=\pm 1/2}u_j\bar{u}_j=(\slashed{p}_j+m_j),\label{outerprod}
\end{equation}

but we don't want to sum over the neutrino spins yet. Thus, we need to derive their outer product without summation using spinor formalism.

Two-component spinors (or bi-spinors) have 2 indices each ranging from 1 to 2. Spinor indices are raised and lowered using the Levi-Civita tensor in 2D[We use the convention $\epsilon^{12}=1$ so $\epsilon_{12}=-1$.]. Next, we define the “soldering forms", which can be used to change vectors into two-component spinors, using the Pauli matrices $\bar{\sigma}$:

\begin{equation}
\sigma^\mu_{\alpha\dot{\beta}} = (1,\bar{\sigma}) \ \ \ \ ;  \ \ \ \ \sigma_{\mu\alpha\dot{\beta}} = (1,-\bar{\sigma}),
\end{equation}

\begin{equation}
\tilde{\sigma}^{\mu\dot{\alpha}\beta} = (1,-\bar{\sigma}) \ \ \ \ ;  \ \ \ \ \tilde{\sigma}_\mu^{\dot{\alpha}\beta} = (1,\bar{\sigma}).
\end{equation}

In the above, $\mu$ is our vector index ranging from 0 to 3, while $\alpha$ and $\beta$ are our spinor components, ranging from 1 to 2. Keep in mind that the “1" is actually a 2D identity matrix.

We also know that our neutrinos satisfy the Dirac equation, so writing their 4-component spinors (and their Dirac conjugates $\bar{u}_s =u_s^\dagger\gamma^0)$ as

\begin{equation}
u_s(\bar{p}) = \begin{pmatrix} x_{s,\alpha}(\bar{p}) \\ y_s^{\dagger\dot{\alpha}}(\bar{p}) \end{pmatrix} \ \ \ \ ;  \ \ \ \ \bar{u}_s^T(\bar{p}) = \begin{pmatrix} y_s^\alpha(\bar{p}) \\ x_{s,\dot{\alpha}}^\dagger(\bar{p}) \end{pmatrix},
\end{equation}

with $s$ being their spin, and selecting our basis (by setting a phase) to be

\begin{equation}
(\chi_{\frac{1}{2}})_\alpha = (\eta^\dagger_{\frac{1}{2}})^{\dot{\alpha}} = \begin{pmatrix}
1 \\ 0
\end{pmatrix} \ \ \ \ ;  \ \ \ \ (\chi_{-\frac{1}{2}})_\alpha = (\eta^\dagger_{-\frac{1}{2}})^{\dot{\alpha}} = \begin{pmatrix}
0 \\ 1 \label{basisspinors}
\end{pmatrix},
\end{equation}

we can solve the Dirac equation for our spinors $u_{s,\alpha}(\bar{p})$. We can also solve for the antiparticle's spinors $v_{s,\alpha}(\bar{p})$, but these are not necessary for our ends. Plugging $u_{s,\alpha}(\bar{p})$ into the Dirac equation and solving similarly to in [17], we find

\begin{equation}
x_{s,\alpha}(\bar{p}) = \tensor{\left(\sqrt{p\cdot\sigma}\right)}{_\alpha^\beta}(\chi_s)_\beta \ \ \ \ ;  \ \ \ \ \ y_s^{\dagger\dot{\alpha}}(\bar{p}) = \tensor{\left(\sqrt{p\cdot\tilde{\sigma}}\right)}{^{\dot{\alpha}}_{\dot{\beta}}}(\eta^\dagger_s)^{\dot{\beta}},
\end{equation}

\begin{equation}
x^\dagger_{s,\dot{\alpha}}(\bar{p}) = (\chi_s^\dagger)_{\dot{\beta}}\tensor{(\sqrt{p\cdot\sigma})}{^{\dot{\beta}}_{\dot{\alpha}}} \ \ \ \ ;  \ \ \ \ y_s^\alpha(\bar{p}) = (\eta_s)^\beta \tensor{(\sqrt{p\cdot\tilde{\sigma}})}{_\beta^\alpha}.
\end{equation}

These give the most general solutions for our spinor and antispinor. Let us now consider what is known as the helicity basis - this is when we choose the spin's direction to align with the momentum's: $\hat{s} = \hat{p}$. This is named as such because the helicity, defined as $\hat{h} = \frac{\bar{\sigma} \cdot \bar{p}}{|\bar{p}|}$ (which is a conserved quantity) satisfies the eigenvalue equation [59]

\begin{equation}
\hat{h}u_s = \pm2su_s
\end{equation}

for the case of a spin-$\frac{1}{2}$ particle. Thus, we see that $\bar{p}\cdot\bar{\sigma} = \pm2s|\bar{p}|$. For the basis spinors in (<ref>), we have that

\begin{equation}
\bar{p}\cdot\bar{\sigma} \chi_s = 2s|\bar{p}| \chi_s \ \ \ \ ;  \ \ \ \ \bar{p}\cdot\bar{\sigma} \eta_s = 2s|\bar{p}| \eta_s.
\end{equation}

Plugging this in, and using that $p\cdot\sigma = p^\mu \sigma_{\mu\alpha\dot{\beta}} = E\mathds{1} - \bar{p}\cdot\bar{\sigma}$ and similarly $p\cdot\tilde{\sigma} = E\mathds{1} + \bar{p}\cdot\bar{\sigma}$, we have

\begin{equation}
x_{s,\alpha}(\bar{p}) = \sqrt{E-2s|\bar{p}|}(\chi_s)_\alpha \ \ \ \ ;  \ \ \ \ y_s^{\dagger\dot{\alpha}}(\bar{p}) = \sqrt{E+2s|\bar{p}|}(\eta^\dagger_s)^{\dot{\alpha}},
\end{equation}

\begin{equation}
x^{\dagger\dot{\alpha}}_s(\bar{p}) = \sqrt{E-2s|\bar{p}|}(\chi_s^\dagger)^{\dot{\alpha}}  \ \ \ \ ;  \ \ \ \ y_{s,\alpha}(\bar{p}) = \sqrt{E+2s|\bar{p}|}(\eta_s)_\alpha.
\end{equation}

Thus, we have our spinor and its conjugate:

\begin{equation}
u_s = \begin{pmatrix}
\sqrt{E-2s|\bar{p}|}(\chi_s)_\alpha \\ \sqrt{E+2s|\bar{p}|}(\eta^\dagger_s)^{\dot{\alpha}}
\end{pmatrix} \ \ ;  \ \
\bar{u}_s = \begin{pmatrix}
\sqrt{E+2s|\bar{p}|}\epsilon^{\alpha\beta}(\eta_s)_\beta & \sqrt{E+2s|\bar{p}|}\epsilon_{\dot{\alpha}\dot{\beta}}(\chi^\dagger_s)^{\dot{\beta}}
\end{pmatrix}. \label{refme}
\end{equation}

The outer product of the spinors for each neutrino species is then

\begin{equation}
u_{\nu_i}\bar{u}_{\nu_i} = \begin{pmatrix}
\sqrt{E^2-4s^2|\bar{p}|^2}(\chi_s)_\alpha\epsilon^{\beta\gamma}(\eta_s)_\gamma &
(E-2s|\bar{p}|)(\chi_s)_\alpha\epsilon_{\dot{\beta}\dot{\gamma}}(\chi^\dagger_s)^{\dot{\gamma}} \\
(E+2s|\bar{p}|)(\eta^\dagger_s)^{\dot{\alpha}}\epsilon^{\beta\gamma}(\eta_s)_\gamma &
\sqrt{E^2-4s^2|\bar{p}|^2}(\eta^\dagger_s)^{\dot{\alpha}}\epsilon_{\dot{\beta}\dot{\gamma}}(\chi^\dagger_s)^{\dot{\gamma}}
\end{pmatrix}.
\end{equation}

Luckily, this is equivalent to a far neater expression [59]:

\begin{equation}
u_{\nu_i}\bar{u}_{\nu_i}=\frac{1}{2}\bigl(\slashed{p}_{\nu_i}+m_{\nu_i}\bigl)\bigl(1+2s_{\nu}\gamma^5\slashed{S}_{\nu_i}\bigl), \label{neutouter}
\end{equation}

where we have specified $s=s_\nu$, and where the spin vector $S_{\nu_i}$ is given by

\begin{equation}
S_{\nu_i}^{\mu}=\left(\frac{|\bar{p}_{\nu}|}{m_{\nu_i}}, \frac{E_{\nu}}{m_{\nu_i}}\frac{\bar{p}_{\nu}}{|\bar{p}_\nu|}\right),
\end{equation}

again specifying now $E=E_\nu$ , $\bar{p} = \bar{p}_\nu$, etc. Note that summing (<ref>) over spins reduces it to (<ref>), as must be.

Plugging in (<ref>) and the other outer products, we have

\begin{equation}
\frac{1}{2}\sum_{s_n,s_p,s_e=\pm 1/2}|\mathcal{M}_i|^2=\frac{G_F^2}{4}|V_{ud}|^2|U_{ei}|^2\tilde{\mathcal{N}}_1^{\alpha\beta}\tilde{\mathcal{N}}_{2\alpha\beta},
\end{equation}

\begin{equation}
\tilde{\mathcal{N}}_1^{\alpha\beta}=\frac{1}{2}\mathrm{tr}\left[\gamma^\alpha\bigl(1-\gamma^5\bigl)\bigl(\slashed{p}_{\nu_i}+m_{\nu_i}\bigl)\bigl(1+2s_\nu\gamma^5\slashed{S}_{\nu_i}\bigl)\gamma^\beta\bigl(1-\gamma^5\bigl)\bigl(\slashed{p}_e+m_e\bigl)\right],
\end{equation}

\begin{equation}
\tilde{\mathcal{N}}_2^{\gamma\delta}=\mathrm{tr}\left[\gamma^\gamma\bigl(f-g\gamma^5\bigl)\bigl(\slashed{p}_n+m_n\bigl)\gamma^\delta\bigl(f-g\gamma^5\bigl)\bigl(\slashed{p}_p+m_p\bigl)\right].
\end{equation}

Multiplying these two expressions, and taking the traces, we have

\begin{multline}
\tilde{\mathcal{N}}_1^{\alpha\beta}\tilde{\mathcal{N}}_{2\alpha\beta}=32\left\{\left(g+f\right)^2\left[\left(p_e\cdot p_p\right)\left(p_{\nu_i} \cdot p_n\right)\right]+\left(g-f\right)^2\left[\left(p_e\cdot p_n\right)\left(p_{\nu_i}\cdot p_p\right)\right]+\left(g^2-f^2\right)m_nm_p\left(p_e\cdot p_{\nu_i}\right)\right\} \\
-64s_\nu m_{\nu_i} \left\{\left(g+f\right)^2\left[\left(p_e\cdot p_p\right)\left(S_{\nu_i} \cdot p_n\right)\right]+\left(g-f\right)^2\left[\left(p_e\cdot p_n\right)\left(S_{\nu_i}\cdot p_p\right)\right]+\left(g^2-f^2\right)m_nm_p\left(p_e\cdot S_{\nu_i}\right)\right\}.
\end{multline}

Before continuing, we select a particular frame in which to work. We will work in the neutron's rest frame, such that

\begin{equation}
p_n^\mu = (m_n, 0),\ \ \ \ p_{\nu}^\mu=(E_{\nu}, \bar{p}_{\nu}),\ \ \ \ p_p^{\mu}=(E_p, \bar{p}_p),\ \ \ \ p_e^{\mu}=(E_e, \bar{p}_e).
\end{equation}

In this case we have

\begin{align}
\tilde{\mathcal{N}}_1^{\alpha\beta}\tilde{\mathcal{N}}_{2\alpha\beta}&=32m_nE_pE_eE_{\nu_i}\biggl\{\left(g+f\right)^2\left(1-\frac{\bar{p}_e\cdot\bar{p}_p}{E_eE_{p}}\right)+(g-f)^2\left(1-\frac{\bar{p}_{\nu}\cdot \bar{p}_p}{E_{\nu_i}E_p} \right) \nonumber \\
&\ \ \ \ +(g^2-f^2)\frac{m_p}{E_p}\left(1-\frac{\bar{p}_e\cdot\bar{p}_{\nu}}{E_eE_{\nu_i}}\right) \biggl\} -64s_{\nu}m_nE_pE_eE_{\nu_i}\biggl\{v_{\nu_i}(g+f)^2\left(1-\frac{\bar{p}_e\cdot\bar{p}_p}{E_eE_p}\right) \nonumber \\
&\ \ \ \ +(g-f)^2\left(v_{\nu_i}-\frac{\bar{p}_{\nu}\cdot\bar{p}_p}{|\bar{p}_{\nu}|E_p}\right) +(g^2-f^2)\frac{m_p}{E_p}\left(v_{\nu_i}-\frac{\bar{p}_{\nu}\cdot\bar{p}_e}{|\bar{p}_{\nu}|E_e} \right) \biggl\} \nonumber \\
&=32m_nE_pE_eE_{\nu_i}\biggl\{2\left(g^2+f^2 \right)\left(1-2s_{\nu}v_{\nu_i}\right)
+\left(g^2-f^2 \right)\frac{m_p}{E_p}\left(1-2s_{\nu}v_{\nu_i}\right) \nonumber \\
&\ \ \ \ +\left(f^2-g^2\right)\frac{m_p}{E_p}\left(v_{\nu_i}-2s_{\nu}\right)v_e\cos\theta_{e\nu}
-\left(g+f\right)^2(1-2s_{\nu}v_{\nu_i})v_ev_p\cos\theta_{ep} \nonumber \\
&\ \ \ \ -\left(g-f\right)^2(v_{\nu_i}-2s_{\nu})v_p\cos\theta_{\nu p} \biggl\},
\label{NN}
\end{align}

where $v_j=|\bar{p}_j|/E_j$ is the velocity of each particle and $\cos\theta_{jk}=\bar{p}_j\cdot\bar{p}_k/(|\bar{p}_j||\bar{p}_k|)$ is the angle between the $j$th and $k$th particles. Next, we can use momentum conservation ($\bar{p}_\nu=\bar{p}_p+\bar{p}_e$) to rewrite all these angles in terms of one: the angle between the incoming neutrino and outgoing electron, $\cos\theta_{e\nu}$:

\begin{equation}
\cos\theta_{ep}=-\frac{|\bar{p}_e|}{|\bar{p}_p|}+\frac{|\bar{p}_{\nu}|}{|\bar{p}_p|}\cos\theta_{e\nu},\ \ \ \
\cos\theta_{\nu p}=\frac{|\bar{p}_{\nu}|}{|\bar{p}_p|}-\frac{|\bar{p}_e|}{|\bar{p}_p|}\cos\theta_{e\nu}.
\end{equation}

Using the above equations then, our squared matrix element becomes

\begin{align}
&\frac{1}{2}\sum_{s_n,s_p,s_e=\pm 1/2}|\mathcal{M}_i|^2= \nonumber \\
&\ \ \ \ 8G_F^2|V_{ud}|^2|U_{ei}|^2m_nE_pE_eE_{\nu_i} \biggl\{2\left(g^2+f^2 \right)\left(1-2s_{\nu}v_{\nu_i}\right)
+\left(g^2-f^2 \right)\frac{m_p}{E_p}\left(1-2s_{\nu}v_{\nu_i}\right) \nonumber \\
&\ \ \ \ +\left(f^2-g^2\right)\frac{m_p}{E_p}\left(v_{\nu_i}-2s_{\nu}\right)v_e\cos\theta_{e\nu}
+\left(g+f\right)^2(1-2s_{\nu}v_{\nu_i})v_e\left(\frac{|\bar{p}_e|}{E_p}-\frac{|\bar{p}_{\nu}|}{E_p}\cos\theta_{e\nu} \right) \nonumber \\
&\ \ \ \  -\left(g-f\right)^2(v_{\nu_i}-2s_{\nu})\biggl(\frac{|\bar{p}_{\nu}|}{E_p}-\frac{|\bar{p}_e|}{E_p}\cos\theta_{e\nu} \biggl) \biggl\}.
\label{M2}
\end{align}

Finally we reach our cross section. Just before though, note that the Mandelstam variable $t$ is given by

\begin{equation}
t = (p_e-p_{\nu_i})^2=(E_e-E_{\nu_i})^2-|\bar{p_e}-\bar{p}_{\nu_i}|^2=(E_e-E_{\nu_i})^2-|\bar{p}_e|^2-|\bar{p}_{\nu_i}^2|+2|\bar{p}_e||\bar{p}_{\nu_i}|\cos\theta_{e\nu_i}
\end{equation}

and also that the centre of mass neutrino momentum $\bar{p}_{\nu_i}^{\mathrm{com}}$ is related to the neutrino momentum in the neutron's rest frame $\bar{p}_{\nu_i}$ by

\begin{equation}
m_n^2|\bar{p}_{\nu_i}|^2 = s|\bar{p}_{\nu_i}^{\mathrm{com}}|^2,
\end{equation}

where $s$ is the Mandelstam variable $s=(p_n+p_{\nu_i})^2$.

Thus, using that

\begin{equation}
\frac{d\sigma_i}{dt}=\frac{1}{64 \pi s}\frac{1}{|\bar{p}_{\nu}^{\rm com}|^2}\ \frac{1}{2}\sum_{s_n,s_e,s_p=\pm 1/2}|\mathcal{M}_i|^2,
\end{equation}

we have

\begin{equation}
\frac{d\sigma_i}{d\cos\theta_{e\nu_i}} = \frac{d\sigma_i}{dt} \frac{dt}{d\cos\theta_{e\nu_i}} = \frac{1}{32\pi}\frac{1}{m_n^2}\frac{|\bar{p}_e|}{|\bar{p}_{\nu_i}|}\  \frac{1}{2}\sum_{s_n,s_e,s_p=\pm 1/2}|\mathcal{M}_i|^2.
\end{equation}

So, plugging in our matrix element, we have

\begin{align}
\frac{d \sigma_i}{d \cos\theta_{e\nu_i}}&=\frac{G_F^2}{4\pi}|V_{ud}|^2|U_{ei}|^2\frac{E_pE_e|\bar{p}_e|}{m_nv_{\nu_i}} \biggl\{2\left(g^2+f^2 \right)\left(1-2s_{\nu}v_{\nu_i}\right)
+\left(g^2-f^2 \right)\frac{m_p}{E_p}\left(1-2s_{\nu}v_{\nu_i}\right) \nonumber \\
&\ \ \ \ +\left(f^2-g^2\right)\frac{m_p}{E_p}\left(v_{\nu_i}-2s_{\nu}\right)v_e\cos\theta_{e\nu}
+\left(g+f\right)^2(1-2s_{\nu}v_{\nu_i})v_e\left(\frac{|\bar{p}_e|}{E_p}-\frac{|\bar{p}_{\nu}|}{E_p}\cos\theta_{e\nu} \right) \nonumber \\
&\ \ \ \  -\left(g-f\right)^2(v_{\nu_i}-2s_{\nu})\biggl(\frac{|\bar{p}_{\nu}|}{E_p}-\frac{|\bar{p}_e|}{E_p}\cos\theta_{e\nu} \biggl) \biggl\}.
\end{align}

and finally integrating over $\cos\theta_{e\nu_i}$ and taking $v_{\nu_i}$ over, we have

\begin{align}
\sigma_iv_{\nu_i}&=\frac{G_F^2}{2\pi}|V_{ud}|^2|U_{ei}|^2\frac{E_pE_e|\bar{p}_e|}{m_n}\biggl\{2\left(g^2+f^2 \right)\left(1-2s_{\nu}v_{\nu_i}\right)
+\left(g^2-f^2 \right)\frac{m_p}{E_p}\left(1-2s_{\nu}v_{\nu_i}\right) \nonumber \\
&\ \ \ \  +\left(g+f\right)^2(1-2s_{\nu}v_{\nu_i})v_e\frac{|\bar{p}_e|}{E_p}  -\left(g-f\right)^2(v_{\nu_i}-2s_{\nu})\frac{|\bar{p}_{\nu}|}{E_p} \biggl\}.
\end{align}

Note that actually, when integrating over $\cos\theta_{e\nu_i}$ an approximation was made, since

\begin{align}
E_p &\simeq m_p\left(1+\frac{|\bar{p}_p|}{2m_p^2}\right)\simeq m_p\left(1+ \frac{|\bar{p}_e|^2}{2m_p^2}-\frac{|\bar{p}_e||\bar{p}_{\nu}|}{m_p^2}\cos\theta_{e\nu_i} \right).\label{protonE}
\end{align}

When integrating, this last term actually falls away. However, this change is only around $10^{-18}$ of the proton's energy, and so we assume that $E_p$ is left unchanged. In fact, to a very high accuracy we can take $E_p$ to be $m_p$, as the other terms in  (<ref>) are at least $10^{5}$ times smaller.

For consistency, in this regard we are able to discard the last 2 terms, proportional to $\frac{|\bar{p}_e|}{E_p}$ and $\frac{|\bar{p}_\nu|}{E_p}$, as they are around $10^{-5}$ and $10^{-13}$ times the leading order terms. Thus, we are left with

\begin{align}
\sigma_iv_{\nu_i}&=\frac{G_F^2}{2\pi}|V_{ud}|^2|U_{ei}|^2\frac{m_pE_e|\bar{p}_e|}{m_n}\left(3g^2+f^2 \right)\left(1-2s_{\nu}v_{\nu_i}\right).
\end{align}

To obtain the full expression for $\sigma_iv_{\nu_i}$ of the reaction in (<ref>), we need to now shoehorn some added factors in. Firstly, the most obvious change is that the neutron and proton masses become the $^3$H and $^3$He masses respectively, and the nucleonic form factors are replaced by transition probabilities, with $f^2$ becoming the Fermi transition probability $\braket{f_F}^2$ and $3g^2$ becoming $\frac{g_A^2}{g_V^2}\braket{g_{GT}}^2$ [60], with $\braket{g_{GT}}$ the Gamow-Teller transition probability and $g_V$ and $g_A$ the vector and axial coupling constants, respectively.

Next, the Fermi function $F(Z,E_e)$ is needed to account for the Coulombic attraction between the outgoing electron and helium nucleus, and is given by

\begin{equation}
F(Z,E_e)=\frac{2\pi\alpha ZE_e/|\bar{p}_e|}{1-e^{-2\pi\alpha ZE_e/|\bar{p}_e|}},
\end{equation}

where $Z$ is the atomic number of the nucleus - in this case $Z=2$ - and $α$ is the fine structure constant. Putting this all in, we have

\begin{align}
\sigma_iv_{\nu_i}&=\frac{G_F^2}{2\pi}|V_{ud}|^2|U_{ei}|^2\frac{m_\mathrm{^3He}E_e|\bar{p}_e|}{m_\mathrm{^3H}}F(2,E_e)\left(\braket{f_F}^2 + \frac{g_A^2}{g_V^2}\braket{g_{GT}}^2 \right)\left(1-2s_{\nu}v_{\nu_i}\right).\label{crosssection}
\end{align}

The term including $v_ν_i$ is important, and is the next-to leading order term. Even for the most massive possible neutrinos, this is still of the order $10^-2$, and so for the accuracy to which we are working, it is not negligible. For massless neutrinos, it is even a leading order term. The minute terms can be neglected as other miniscule effects have been as well, such as the differences incurred by using the full electroweak theory, as well as those brought about by the fact that atomic tritium is used - not free nuclei.

\subsection{The Capture Rate of Relic Neutrinos on Tritium}

Considering some lone tritium atom in a sea of neutrinos, the usual rate for this atom to capture a neutrino would be $Γ=⟨nσv|$⟩, with $n$ the number density of the surrounding neutrinos, $\sigma$ the interaction cross section and $v$ the neutrino's velocity. In actuality, this is just an approximation, with the real capture rate being

\begin{equation}
\Gamma_i = \sum_{s_\nu=\pm 1/2} \int \frac{d^3p_\nu}{(2\pi)^3}f_{\nu_i}(\bar{p}_\nu)\sigma_{\nu_i}(\bar{p}_\nu,s_\nu)v_{\nu_i}.
\end{equation}

Considering an experimental setup that involves a mass $M_T$ of atomic tritium with atomic mass $m^{\mathrm{atom}}_{^3\mathrm{H}}$, there are now $N_T=\frac{M_T}{m^{\mathrm{atom}}_{^3\mathrm{H}}}$ interaction points. Atop this, there are 3 species in the mass basis which need to be summed over. Thus, the full capture rate of cosmic neutrinos on this setup of tritium is

\begin{equation}
\Gamma_{\mathrm{C\nu B}} = N_T\sum_{i=1}^{3}\sum_{s_\nu=\pm 1/2} \int \frac{d^3p_\nu}{(2\pi)^3}f_{\nu_i}(\bar{p}_\nu)\sigma_{\nu_i}(\bar{p}_\nu,s_\nu)v_{\nu_i}.
\end{equation}

The factor of $\sigma_{\nu_i}(\bar{p}_\nu,s_\nu)v_{\nu_i}$ was just calculated in (<ref>) above. As discussed in chapters 4 and 5, besides for non-standard physics, the distribution function can be written as

\begin{equation}
f_{\nu_i} = f_0(1+\delta f_{\nu_i}^d +\delta f_{\nu_i}^c)
\end{equation}

for both left- and right-handed neutrinos if they are Majorana fermions, and for just left-handed neutrinos if they are Dirac particles. As a reminder, $f_0$ is the massless Fermi-Dirac distribution function, with $\delta f_{\nu_i}$ being the change to this distribution function owing to non-instantaneous decoupling and gravitational clustering.

Before continuing, we need to consider the kinematics of the situation to see how the ejected electron's energy and momentum are related to the incoming neutrino's momentum. To begin, we start by analysing normal beta decay by a tritium nucleus : $^3\mathrm{H} \ \rightarrow \ ^3\mathrm{He} + e^- + \bar{\nu}_i$. In this case, the kinetic energy of the emitted electron, $K_e=E_e-m_e$, is given by

\begin{equation}
K_e=\frac{(m_{\mathrm{^3H}}-m_e)^2-m_{\nu_i}^2-m_{\mathrm{^3He}}^2-2E_{\nu_i}E_{\mathrm{^3He}}+2|\bar{p}_{\nu}||\bar{p}_{\mathrm{^3He}}|\cos \theta_{\mathrm{^3He}\nu }}{2m_{\mathrm{^3H}}}.
\end{equation}

Recall that during beta decay, the emitted electron can take on many possible values, and has a characteristic spectrum - this was one of the first signs that a neutrino was needed in the standard model. The electron's possible energy does have a cutoff, however: the “endpoint", the maximal possible kinetic energy that the electron can have, is achieved when the electron is emitted anti-parallel to the helium and anti-neutrino, and when the anti-neutrino emitted is in the lightest mass eigenstate. In this case then, $\cos \theta_{\mathrm{^3He}\nu } = 1$, $m_{\nu_i}=m_{\mathrm{lightest}}$, and the maximisation of $E_e$ with respect to $|\bar{p}_\nu|$ leads to a minimisation of $E_\nu+E_{^3\mathrm{He}}$ with respect to $|\bar{p}_\nu|$, which gives

\begin{equation}
\frac{|\bar{p}_\nu|}{|\bar{p}_{^3\mathrm{He}}|} = \frac{m_{\nu_i}}{m_{^3\mathrm{He}}}.
\end{equation}

Plugging these in, we have

\begin{equation}
K_{\mathrm{end}} = \frac{(m_{\mathrm{^3H}}-m_e)^2-(m_{\rm lightest}+m_{\mathrm{^3He}})^2}{2m_{\mathrm{^3H}}},\label{Kend}
\end{equation}

independent of everything except masses. We can define $K_{\mathrm{end}}^0$, the maximum possible kinetic energy for the emitted electron if the lightest species of neutrinos is massless:

\begin{equation}
K_{\mathrm{end}}^0 = \frac{(m_{\mathrm{^3H}}-m_e)^2-m_{\mathrm{^3He}}^2}{2m_{\mathrm{^3H}}}. \label{Kend0}
\end{equation}

Next, let us move on to the reaction at hand: neutrino capture. In this case, the kinetic energy of the emitted electron is slightly altered, and is given by

\begin{align}
K_e^{\rm{C\nu B}}&=\frac{(E_{\nu_i}+m_{\mathrm{^3H}}-m_e)^2-|\bar{p}_{\nu}|^2+2|\bar{p}_{\nu}||\bar{p}_e|\cos \theta_{e\nu}-m_{\mathrm{^3He}}^2}{2(E_{\nu_i}+m_{\mathrm{^3H}})} \nonumber \\
&\simeq \frac{(E_{\nu_i}+m_{\mathrm{^3H}}-m_e)^2-m_{\mathrm{^3He}}^2}{2(E_{\nu_i}+m_{\mathrm{^3H}})},
\end{align}

where in the second line we have neglected terms proportional to $|\bar{p}_{\nu}|^2$ and $2|\bar{p}_{\nu}||\bar{p}_e|$, as these terms are much smaller than the mass scale of tritium and helium. By similar logic, keeping only the leading order term in $E_{\nu_i}$, we have that

\begin{equation}
K_e^{\rm{C\nu B}} \approxeq K_{\mathrm{end}}^0 + E_{\nu_i}.
\end{equation}

Thus, for the $E_e$ appearing in our capture rate, we have $E_e \approxeq m_e + K_{\mathrm{end}}^0 + E_{\nu_i}$, with $|\bar{p}_e| = \sqrt{E_e^2-m_e^2}$.

In order to simplify our calculation, we can make the approximation

\begin{equation}
E_e \approxeq m_e + K_{\mathrm{end}}^0 + m_{\nu_i},
\end{equation}

so that we can remove both $E_e$ and $|\bar{p}_e|$ from the integral over $p_{\nu}$. This approximation is certainly justified: for massive neutrinos, their masses are much larger than the characteristic momentum seen in (<ref>), and for massless neutrinos, this selfsame momentum in (<ref>) is obviously much smaller than the energy scale of $m_e$.

Going back to our capture rate then, we now have

\begin{multline}
\Gamma_{\mathrm{C\nu B}} = N_T\sum_{i=1}^{3}\sum_{s_\nu=\pm 1/2}\frac{G_F^2}{2\pi}|V_{ud}|^2|U_{ei}|^2\frac{m_\mathrm{^3He}E_e|\bar{p}_e|}{m_\mathrm{^3H}}F(2,E_e)\left(\braket{f_F}^2 + \frac{g_A^2}{g_V^2}\braket{g_{GT}}^2 \right) \times \\
\int \frac{d^3p_\nu}{(2\pi)^3}f_0(1+\delta f_{\nu_i}^d +\delta f_{\nu_i}^c)\left(1-2s_{\nu}v_{\nu_i}\right),
\end{multline}

or rather

\begin{multline}
\Gamma_{\mathrm{C\nu B}} = N_T\sum_{i=1}^{3}\sum_{s_\nu=\pm 1/2}\frac{G_F^2}{2\pi}|V_{ud}|^2|U_{ei}|^2\frac{m_\mathrm{^3He}E_e|\bar{p}_e|}{m_\mathrm{^3H}}F(2,E_e)\left(\braket{f_F}^2 + \frac{g_A^2}{g_V^2}\braket{g_{GT}}^2 \right)\times \\
(n_0+\delta n_{\nu_i} -2s_\nu n_0 \braket{v_{\nu_i}} - 2s_\nu\braket{\delta v_{\nu_i}}),\label{capturerate}
\end{multline}

with $n_0$ given by (<ref>), $\delta n_{\nu_i} = \delta n_{\nu_i}^d + \delta n_{\nu_i}^c$ given by (<ref>) and whose numerical values were seen in chapter 5, $\braket{v_{\nu_i}}$ the normalised expectation value of the neutrino's momentum using the massless Fermi-Dirac distribution, and

\begin{equation}
\braket{\delta v_{\nu_i}} = \int\frac{d^3p_\nu}{(2\pi)^3} f_0(\bar{p}_\nu)(\delta f_{\nu_i}^d + \delta f_{\nu_i}^c)v_{\nu_i}
\end{equation}

the unnormalised expectation of the change to the velocity owing to the clustering and spectral distortive effects.

Some interesting quantities to seperate are

\begin{align}
\delta \Gamma_i^d&=N_T\frac{G_F^2}{2\pi}|V_{ud}|^2|U_{ei}|^2\frac{m_\mathrm{^3He}}{m_\mathrm{^3H}}\left(\langle f_F \rangle^2+\frac{g_A^2}{g_V^2}\langle g_{GT} \rangle^2 \right)F(2,E_e)E_e|\bar{p}_e|
\sum_{s_{\nu}=\pm\frac{1}{2}}( \delta n_{\nu_i}^d  -2s_\nu \langle \delta v_{\nu_i}^d \rangle),  \nonumber \\
\delta \Gamma_i^c&=N_T\frac{G_F^2}{2\pi}|V_{ud}|^2|U_{ei}|^2\frac{m_\mathrm{^3He}}{m_\mathrm{^3H}} \left(\langle f_F \rangle^2+\frac{g_A^2}{g_V^2}\langle g_{GT} \rangle^2 \right)F(2,E_e)E_e|\bar{p}_e|
\sum_{s_{\nu}=\pm\frac{1}{2}}( \delta n_{\nu_i}^c  -2s_\nu \langle \delta v_{\nu_i}^c \rangle),
\end{align}

the changes to the capture rate for each species $i$ owing to the effects of spectral distortion due to non-instantaneous decoupling and gravitational clustering.

\subsection{Numerical Results}

As we have seen, particularly in this chapter and the discussions in chapter 5, there are many factors that play a role in the expected capture rate of relic neutrinos. The two main factors, making the largest difference to the situation are the choice of hierarchy (normal or inverted) and the type of fermion (Dirac or Majorana). Thus, this is how our results will be divided, with other smaller effects being considered within each category. In this section, we consider only standard physics, with the added effects of new physics being covered in a later subsection.\\

The numerical values calculated were done so using both Python and Matlab, and the values arrived at were in agreement to the accuracy stated below. Below are the values used for the numerical calculations that were not already stated.

\begin{center}
\begin{tabular}{||c|c||}
$m_{^3\mathrm{H}}$ & $2.808921\times 10^9$ eV [61]\\
$m^{\mathrm{atom}}_{^3\mathrm{H}}$ & $2.809432\times 10^9$ eV [61]\\
$m_{^3\mathrm{He}}$ & $2.808391\times 10^9$ eV [61]\\
$m_{e}$ & $5.1099891\times 10^5$ eV [62]\\
$\braket{f_F}$ & 0.9998 [60]\\
$\braket{g_{GT}}$ & $\sqrt{3}\times (0.9511\pm0.0013)$ [60]\\
$g_A$ & 1.2723 [60]\\
$g_V$ & 1 [60]\\
$V_{ud}$ & 0.974 [51]\\
$|U_{ej}|^2$ & [0.681 , 0.297 , 0.0222]  [22]\\
\end{tabular}
\end{center}

Using these values, we can calculate a characteristic cross section for this reaction following the lead in [43]:

\begin{equation}
\bar{\sigma} = \frac{G_F^2}{2\pi}|V_{ud}|^2\frac{m_\mathrm{^3He}E_e|\bar{p}_e|}{m_\mathrm{^3H}}F(2,E_e)\left(\braket{f_F}^2 + \frac{g_A^2}{g_V^2}\braket{g_{GT}}^2 \right) \approxeq 3.835 \times 10^{-45} \mathrm{cm}^2. \label{charcrosssec}
\end{equation}

For our numerical results we considered a tritium target weighing 100 grams, in line with the planned project PTOLEMY. From the mass differences seen in \eqref{massdiff} and the constraint in \eqref{massupperlim}, the possible masses for neutrinos are limited but not fixed. For the sake of this calculation, for both the normal and inverted hierarchies, the minimal allowed masses were used: that is, the lightest species was taken to be massless.

This was done for two reasons. Firstly, it allows us to analyse relic neutrinos that are still relativistic today, and secondly, it is the only way that the effects of non-instantaneous decoupling could be observed behind the shroud of gravitational clustering, as we saw in chapter 5.

Finally, before stating the results, for the sake of the calculation the $⟨δv_ν_i|$⟩ terms were neglected. Working to next-to leading order accuracy, both the velocity and effects altering the distribution function (non-instantaneous decoupling and gravitational clustering) were included, but their product would be an even smaller correction, akin to annual modulation, helicity flipping and anistropy, already neglected above.

For the normal ordering case, with masses given by (<ref>), for the case where neutrinos are Majorana particles, we have the total capture rate for each massive species $\Gamma_i^M$, its deviation originating from the spectral distortion from non-instantaneous decoupling $\delta \Gamma_i^{Md}$, and that from the gravitational clustering effects $\delta \Gamma_i^{Mc}$, given by (considering $100$ g of tritium)

\begin{align}
\Gamma_1^M\simeq 5.48\ {\rm yr^{-1}},  \ \ \ \  \Gamma_2^M\simeq 2.40\ {\rm yr^{-1}},  \ \ \ \   \Gamma_3^M\simeq 0.200 \ {\rm yr^{-1}},
\label{ValueM}
\end{align}

\begin{align}
\delta \Gamma_1^{Md}\simeq 0.061\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_2^{Md}\simeq 0.024\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_3^{Md}\simeq 1.6 \times 10^{-3} \ {\rm yr^{-1}},
\label{DGMd}
\end{align}

\begin{align}
\delta \Gamma_1^{Mc}\simeq 0\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_2^{Mc}\simeq 0.013\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_3^{Mc}\simeq 0.021 \ {\rm yr^{-1}},
\label{DGMc}
\end{align}

and therefore $\Gamma_{C\nu B}^M\approxeq 8.08{\rm yr^{-1}}$.

Similarly for the Dirac case, we have

\begin{align}
\Gamma_1^D\simeq 5.48\ {\rm yr^{-1}},  \ \ \ \ \Gamma_2^D\simeq 1.27\ {\rm yr^{-1}},  \ \ \ \  \Gamma_3^D\simeq 0.101 \ {\rm yr^{-1}},
\end{align}

\begin{align}
\delta \Gamma_1^{Dd}\simeq 0.061\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_2^{Dd}\simeq 0.012\ {\rm yr^{-1}},  \ \ \ \   \delta \Gamma_3^{Dd}\simeq 8.0 \times 10^{-4} \ {\rm yr^{-1}},
\end{align}

\begin{align}
\delta \Gamma_1^{Dc}\simeq 0\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_2^{Dc}\simeq 6.3 \times 10^{-3}\ {\rm yr^{-1}},  \ \ \ \   \delta \Gamma_3^{Dc}\simeq 0.011\ {\rm yr^{-1}},
\end{align}

and so $\Gamma_{C\nu B}^D\approxeq 6.851{\rm yr^{-1}}$.

It is interesting to note that the ratios between the capture rates for Dirac and Majorana neutrinos are

\begin{align}
\Gamma_1^M/\Gamma_1^D=1, \ \ \ \ \Gamma_2^M/\Gamma_2^D\simeq 1.89, \ \ \ \ \Gamma_3^M/\Gamma_3^D\simeq 1.98,
\end{align}

as expected. For relativistic neutrinos, owing to the factor of $1-2s_\nu v_\nu$, even in the Majorana case where right-helical states contain some left-chiral component, and therefore can interact, right handed neutrinos would not contribute as $1-2s_\nu v_\nu=0$ when the spin is $\frac{1}{2}$ and the velocity is 1.

For non-relativistic (massive) neutrinos, this factor of approximately 2 (altered slightly by clustering and velocity effects) is of course owing to the fact that for Dirac neutrinos only left-helical neutrinos can interact with the tritium, whereas for Majorana neutrinos both left- and right-helical neutrinos can be captured.

It is interesting to note that in the Majorana case, the capture rate's mass dependence only comes in in the clustering effect - the only other term that the neutrino mass affects is the velocity term, which cancels out in the Majorana case as the left- and right-helical components are summed over.

For the inverted hierarchy case, with masses given by (<ref>), once again we have the total capture rate for each massive species $\Gamma_i^M$, its deviation originating from the spectral distortion from non-instantaneous decoupling $\delta \Gamma_i^{Md}$, and that from the gravitational clustering effects $\delta \Gamma_i^{Mc}$, given by (considering $100$ g of tritium)

\begin{align}
\Gamma_1^M\simeq 6.13\ {\rm yr^{-1}},  \ \ \ \  \Gamma_2^M\simeq 2.67\ {\rm yr^{-1}},  \ \ \ \   \Gamma_3^M\simeq 0.178 \ {\rm yr^{-1}},
\end{align}

\begin{align}
\delta \Gamma_1^{Md}\simeq 0.061\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_2^{Md}\simeq 0.024\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_3^{Md}\simeq 1.6 \times 10^{-3} \ {\rm yr^{-1}},
\end{align}

\begin{align}
\delta \Gamma_1^{Mc}\simeq 0.65\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_2^{Mc}\simeq 0.28\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_3^{Mc}\simeq 0 \ {\rm yr^{-1}},
\end{align}

and so $\Gamma_{C\nu B}^M\approxeq 8.978{\rm yr^{-1}}$. Note that we take the same values of $\delta n_1^c$ for $m_1=49.3\ {\rm meV}$ and $\delta n_2^c$ for $m_2=50\ {\rm meV}$ [63].

Next, for the case of Dirac neutrinos, we have

\begin{align}
\Gamma_1^D\simeq 3.10\ {\rm yr^{-1}},  \ \ \ \ \Gamma_2^D\simeq 1.35\ {\rm yr^{-1}},  \ \ \ \  \Gamma_3^D\simeq 0.178 \ {\rm yr^{-1}},
\end{align}

\begin{align}
\delta \Gamma_1^{Dd}\simeq 0.031\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_2^{Dd}\simeq 0.012\ {\rm yr^{-1}},  \ \ \ \   \delta \Gamma_3^{Dd}\simeq 1.6 \times 10^{-3}\ {\rm yr^{-1}},
\end{align}

\begin{align}
\delta \Gamma_1^{Dc}\simeq 0.33\ {\rm yr^{-1}},  \ \ \ \  \delta \Gamma_2^{Dc}\simeq 0.14\ {\rm yr^{-1}},  \ \ \ \   \delta \Gamma_3^{Dc}\simeq 0\ {\rm yr^{-1}}.
\end{align}

Once again, we see the ratios between the capture rates of each species for the Dirac and Majorana cases are

\begin{align}
\Gamma_1^M/\Gamma_1^D \simeq 1.98, \ \ \ \ \Gamma_2^M/\Gamma_2^D\simeq 1.98, \ \ \ \ \Gamma_3^M/\Gamma_3^D = 1,
\end{align}

again displaying unity for relativistic species and approximately 2 for the massive species.

As stated in chapter 5, we could also consider helicity flipping. In the case of complete helicity flipping, which in the Dirac case entails half of the left-helical population to become right-helical, at leading order this does not make a difference as helicities are summed over. However, the next-to leading order term which involves the velocity now cancels out. In the relativistic case, this makes no difference as the $\sum_{s_\nu=\pm 1/2}(1-2s_\nu v_{\nu_i})$ factor is equal to 2 either way. In the non-relativistic case however, the small correction of $v_{\nu_i}$ is lost. In a more realistic case, where there is only partial helicity flipping, and only clustered neutrinos' helicities being flipped, only a fraction of this will be lost, changing the capture rate at the order $\approx10^{-4} - 10^{-3} \mathrm{yr}^{-1}$.

As stated already, the mass dependence of this capture rate is tied to its velocity dependence, which cancels out in the Majorana case. We can see in figure <ref> how the capture rate is affected in the Dirac case for the lightest species, with the Majorana capture rate as a benchmark.

The relationship between the capture rate and the mass of the lightest neutrino, with Dirac cases shown in blue and Majorana in red. The left panel depicts the normal hierarchy case (that is, the lightest neutrino is $m_1$) and the right panel is for the inverted hierarchy case (where $m_3$ is the smallest mass). Taken from [63].

§.§ Beta Decay Rate's Energy Dependence

The main source of background noise in our chosen experiment is the regular beta decay of the tritium nuclei. In order to see how the spectrum of this noise will look, we consider the three body decay process of tritium into helium-3, an electron and an anti-neutrino:

\begin{equation}
^3\mathrm{H}(M,0)\rightarrow\ ^3\mathrm{He}(E',\bar{p}') + e^-(E_e,\bar{p}_e) + \bar{\nu}_e(E_\nu,\bar{p}_\nu).
\end{equation}

We wish to analyse the dependence of the decay rate $\Gamma$ on the electrons energy: $\frac{d\Gamma}{dE_e}$. To do so, we begin with the differential decay rate [17]:

\begin{equation}
d\Gamma = \frac{1}{2E_i}|\mathcal{M}|^2d\Pi_{LIPS},
\end{equation}

where $E_i=M$, $|\mathcal{M}|^2$ is the squared matrix element, and $d\Pi_{LIPS}$ is the Lorentz invariant phase space. Expanding this yields

\begin{equation}
d\Gamma = \frac{1}{2^9\pi^5M}\int\frac{d^3p'd^3p_ed^3p_\nu}{E'E_eE_\nu}|\mathcal{M}|^2\delta^4(p_H-p_{He}-p_e-p_\nu).
\end{equation}

We can then use the 3 dimensional delta function to kill the $d^3p'$ integral. This will leave $E'=\sqrt{(\bar{p}_e+\bar{p}_\nu)^2+M'^2}=\sqrt{\bar{p}_e^2+\bar{p}_\nu^2+2|\bar{p}_e||\bar{p}_\nu|\cos\theta+M'^2}$, where $M'$ is the mass of the helium-3.

The remaining delta function is used to remove one angular integral. We have that

\begin{equation}
\delta(M-E'-E_e-E_\nu)=\delta(M-E_e-E_\nu-(\bar{p}_e^2+\bar{p}_\nu^2+2|\bar{p}_e||\bar{p}_\nu|\cos\theta+M'^2)^{1/2}).
\end{equation}

We then use the property of delta functions that $\delta(f(x))=\frac{\delta(x-x*)}{f'(x*)}$, for every root $x*$ of $f(x)$ to obtain

\begin{equation}
\delta(M-E'-E_e-E_\nu)=\delta(\cos\theta-\frac{(M-E_e-E_\nu)^2-\bar{p}_e^2-\bar{p}_\nu^2-M'^2}{2|\bar{p}_e||\bar{p}_\nu|})\frac{M-E_e-E_\nu}{|\bar{p}_e||\bar{p}_\nu|}.
\end{equation}

From here we will use the notation that $p_i=\bar{p}_i$, unless otherwise specified. We have from the above energy delta function that $E'=M-E_e-E_\nu$. Plugging this in and doing the remaining 3 trivial angular integrals yields

\begin{equation}
d\Gamma = \frac{2\cdot2\pi\cdot4\pi}{2^9\pi^5M}\int\frac{dp_edp_\nu}{E_eE_\nu(M-E_e-E_\nu)}\frac{p_e^2p_\nu^2(M-E_e-E_\nu)}{p_ep_\nu}|\mathcal{M}|^2.
\end{equation}

We now change our variable of integration from momentum to energy, and since $E_i^2=p_i^2+m_i^2$, $p_idp_i=E_idE_i$ and we have that

\begin{equation}
d\Gamma = \frac{1}{(2\pi)^38M}\int dE_edE_\nu|\mathcal{M}|^2.
\end{equation}

Next we consider the form of $|\mathcal{M}|^2$. Using symmetry, the general form, up to momenta squared, is given by

\begin{equation}
\end{equation}

where in the above we used 4-momenta and used bars to clearly denote 3-momenta. We also used the fact that $E'=M-E_e-E_\nu$.

We can use 4-momentum conservation, $p_H-p_{He}=p_e+p_\nu$ to rewrite $\bar{p}_e\cdot\bar{p}_\nu$. We have that

\begin{equation}
(p_H-p_{He})^2=(p_e+p_\nu)^2 \rightarrow M^2+M'^2-2M(M-E_e-E_\nu) = m_e^2+m_\nu^2+2E_eE_\nu-2\bar{p}_e\cdot\bar{p}_\nu,
\end{equation}

and so

\begin{equation}
\bar{p}_e\cdot\bar{p}_\nu = \frac{1}{2}(M^2-M'^2+m_e^2+m_\nu^2-2ME_e+2E_\nu E_e-2ME_\nu),
\end{equation}

and thus we will insert into the $d\Gamma$ integral that

\begin{equation}
|\mathcal{M}|^2 = A+\frac{B}{2}(M'^2-M^2-m_e^2-m_\nu^2+2ME_e+2ME_\nu)+CM(M-E_e-E_\nu).
\end{equation}

Now recall that we want an expression for $\frac{d\Gamma}{dE_e}$. So we have now that

\begin{equation}
\frac{d\Gamma}{dE_e}= \frac{1}{(2\pi)^38M}\int_{E_\nu^{min}}^{E_\nu^{max}}dE_\nu\left[A+\frac{B}{2}(M'^2-M^2-m_e^2-m_\nu^2+2ME_e+2ME_\nu)+CM(M-E_e-E_\nu)\right].
\end{equation}

This integration is trivial and we obtain

\begin{multline}
\frac{d\Gamma}{dE_e} = \frac{1}{(2\pi)^38M}[(A+\frac{B}{2}(M'^2-M^2-m_e^2-m_\nu^2+2ME_e)+CM(M-E_e))(E_\nu^{max}-E_\nu^{min})\\
\end{multline}

where we have used that $(E_\nu^{max})^2-(E_\nu^{min})^2=(E_\nu^{max}-E_\nu^{min})(E_\nu^{max}+E_\nu^{min})$.

Now, following the notation used in [64] and [65], we define

\begin{equation}
\end{equation}

where $p_i$ is the 4-momentum of a decay product. Naming $E_1=E_\nu$, $E_2=E'$ and $E_3=E_e$, specifically, we have

\begin{equation}
\end{equation}

Similarly, we have $m_{23}^2=M^2+m_\nu^2-2ME_\nu$, and so the maximisation/minimisation of $m_{23}^2$ corresponds with that of $E_\nu$, as the masses are constant.
Expanding $m_{23}^2=(p'+p_e)^2$, we have that it is maximised/minimised when $\bar{p}'$ is parallel/anti-parallel to $\bar{p}_e$.

This leads to the maximisation/minimisation conditions:

\begin{equation}
\end{equation}

where $E^*$ is the energy of each particle in the rest frame of $m_{12}$. For example,

\begin{equation}
(m_{12}-p_2)^2=m_{12}^2-2m_{12}E^{*}_2+m_2^2=m_1^2\ \ \ \rightarrow \ \ \ E^*_2=E'^*=\frac{m_{12}^2-m_\nu^2+M'^2}{2m_{12}}.
\end{equation}

Similarly $E_e^*=\frac{M^2-m_{12}^2-m_e^2}{2m_{12}}$.

Adding and subtracting $m_{23}^{2max}$ and $m_{23}^{2min}$ gives 2 factorisations:

\begin{equation}
\end{equation}

\begin{equation}
\end{equation}

After plugging in the expressions for each $E^*$, and defining $E_e^{max}=\frac{1}{2M}(M^2+m_e^2-(m_\nu+M')^2)$ as the maximum electron energy, we are left with

\begin{equation}
E_\nu^{max}-E_\nu^{min}=\frac{2Mp_e}{m_{12}^2}\sqrt{(E_e^{max}-E_e)(E_e^{max}-E_e+\frac{2m_\nu M'}{M})},
\end{equation}

\begin{equation}
\end{equation}

In the 4-Fermi approximation, we saw previously that $|\mathcal{M}|^2$ is proportional to $E_eE_\nu$ with a constant of order $G_F^2$. Thus, we can set $A=C=0$ and we get

\begin{multline}
\frac{d\Gamma}{dE_e} = \frac{B}{(2\pi)^3}\frac{p_e}{4m_{12}^2}\sqrt{y(y+\frac{2m_\nu M'}{M})}(\frac{M^2}{m_{12}^2}(M-E_e)(y+\frac{m_\nu}{M}(M'+m_\nu))\\
\end{multline}

wherein $y=E_e^{max}-E_e$. We then perform some algebra on the last bracket above:

\begin{equation}
M^2-M'^2+m_e^2+m_\nu^2-2ME_e=-M^2-m_e^2+(m_\nu+M')^2-2m_\nu^2-2M'm_\nu+2ME_e = 2M(-E_e^{max}+E_e-\frac{m_\nu^2+M'm_\nu}{M}),
\end{equation}

and, plugging this back in and using the definition of $y$, we add the 2 terms and we obtain

\begin{equation}
\frac{d\Gamma}{dE_e} = \frac{B}{(2\pi)^3}\frac{Mp_e}{4m_{12}^4}\sqrt{y(y+\frac{2m_\nu M'}{M})}(ME_e-m_e^2)(y+\frac{m_\nu}{M}(M'+m_\nu)).
\end{equation}

Taking out some factors and cancelling, we arrive at the form

\begin{equation}
\frac{d\Gamma}{dE_e} = B'E_ep_e\frac{1-\frac{m_e^2}{E_eM}}{(1-\frac{2E_e}{M}+\frac{m_e^2}{M^2})^2}\sqrt{y(y+\frac{2m_\nu M'}{M})}(y+\frac{m_\nu}{M}(M'+m_\nu))
\end{equation}

and, in defining $H(E_e,m_\nu)$, we have $\frac{d\Gamma}{dE_e} = B'E_ep_eH(E_e,m_\nu)$.

Now, we know that $B'$ must be of order $G_F^2$ as it is a Fermi interaction, we know it must have a factor $|V_{ud}|^2$ and $|U_{ej}|^2$ for each species in the mass basis and it must be multiplied by $N_T$ for the number of beta-decaying sites. Adding the small correction factors of Fermi and Gamow-Teller transition probabilities owing to the tritium and helium-3, a factor of $\frac{M'}{M}$ from the kinematics, as well as a Fermi function to account for the Coulombic attraction of the outgoing electron to the nucleus, we are left with

\begin{equation}
\frac{d\Gamma}{dE_e} = \sum_{j=1}^{3}|U_{ej}|^2\frac{\bar{\sigma}}{\pi^2}N_TH(E_e,m_\nu),\label{betadec}
\end{equation}

where $\bar{\sigma}$ is defined as in (<ref>) as

\begin{equation}
\bar{\sigma} = \frac{G_F^2}{2\pi}|V_{ud}|^2F(Z,E_e)\frac{M'}{M}E_ep_e (\braket{f_F}^2+(\frac{g_A}{g_V})^2\braket{g_{GT}}^2),
\end{equation}

which is obtained when computing the cross section of the 4-Fermi interaction. Of course, we must remember that this expression is not perfectly relativistic, as even though $H(E_e,m_\nu)$ is, this $\bar{\sigma}$ comes from an effective Lagrangian.

§.§ The Expected Energy Spectra

From the previous subsections, we can now model both neutrino capture and beta decay and see graphically how they may look in an experiment like PTOLEMY. In this experiment, recall that it is the emitted electron that is observed, and its energy is what signals a C$\nu$B interaction. As we will see below, one of the main challenges for observing the C$\nu$B is the distinction of the signal from the background noise that is the constant supply of electrons from usual beta decay.

Recall from the kinematical discussions above that the maximum possible energy that an electron can have as a result of beta decay - the “beta decay endpoint" - is

\begin{equation}
E_{\mathrm{end}}\approxeq K_{\mathrm{end}}^0+m_e-m_{\mathrm{lightest}}.\label{Eend}
\end{equation}

A more accurate value is given by $m_e+K_{\mathrm{end}}$, taken from (<ref>), but (<ref>) is correct to linear order and since $m_{\mathrm{lightest}}\ll m_e$, this approximation is valid.

On the other hand, when an electron is emitted owing to the tritium capturing a relic neutrino, we have

\begin{equation}
E^{\mathrm{C\nu B}}\approxeq K_{\mathrm{end}}^0+m_e +E_{\nu_i},
\end{equation}

which is clearly always larger that (<ref>). Again, this was approximated to linear order as above. What is important is that the difference in energy between the two is

\begin{equation}
\Delta E_e = E^{\mathrm{C\nu B}} - E_{\mathrm{end}} \approxeq E_{\nu_i} + m_{\mathrm{lightest}}.
\end{equation}

Thus, the “smoking gun" characteristic signal of detections of the C$\nu$B would be an emitted electron energy spectrum with peaks a distance $\Delta E_e$ away from the beta decay endpoint.

From (<ref>) and its derivation, we know the spectrum of usual beta decay, and from (<ref>) we know the spectrum of electrons emitted owing to C$\nu$B interactions. Unfortunately, in reality, the energy resolution of the detector is not perfect, and can only measure the electron's energy to a certain accuracy. In order to take this experimental fact into account, instead of plotting the theoretical, exact spectra, we must rather plot Gaussian-smeared versions of them, with the FWHM (full width at half max - the spread of the Gaussian at the half-way point of its height) $\Delta$, where $\Delta$ is the energy resolution of the detector. These Gaussian-smeared spectra are given by [43]

\begin{align}
\frac{d\tilde{\Gamma}_i}{dE_e} &=\frac{1}{\sqrt{2\pi}\sigma}\int^{\infty}_{-\infty}dE_e'\  \Gamma_i(E_e')\ \delta[E_e'-(E_{\rm end}+E_{\nu_i}+m_{\rm lightest})]\exp\left[-\frac{(E_e'-E_e)^2}{2\sigma^2} \right],
\label{tildeGi}
\\
\frac{d\tilde{\Gamma}_{\beta}}{dE_e}&=\frac{1}{\sqrt{2\pi}\sigma}\int^{\infty}_{-\infty}dE_e' \ \frac{d \Gamma_{\beta}}{dE_e}(E_e')\ \exp\left[-\frac{(E_e'-E_e)^2}{2\sigma^2} \right]\label{tildebeta},
\end{align}

with $\Gamma_i(E_e')$ from (<ref>) and $\frac{d \Gamma_{\beta}}{dE_e}(E_e')$ from (<ref>). Here $\sigma$ is the standard deviation of the Gaussian (not a cross section), and can be related to $\Delta$.

To see how, consider a normalised continuous Gaussian distribution with standard deviation $\sigma$, given by $y=\frac{1}{\sqrt{2\pi}\sigma}e^{-x^2/2\sigma^2}$, with $y$ its height and $x$ half its width. Inverting this relation yields $x=\sigma\sqrt{-\ln(2\pi\sigma^2y^2)}$. At the maximum (when $x=0$), $y=\frac{1}{\sqrt{2\pi}\sigma}$ so at half this maximum $y=\frac{1}{2\sqrt{2\pi}\sigma}$. At this height, the full width $\Delta$ is given by twice $x$, so

\begin{equation}
\Delta = 2x|_{\mathrm{HM}} = 2\sigma\sqrt{-\ln(\frac{2\pi\sigma^2}{8\pi\sigma^2})} = \sigma\sqrt{8\ln(2)}.
\end{equation}

The results of these smeared spectra (<ref>) and (<ref>) are seen in figures <ref> and <ref> - this is what we would actually expect to observe.

The expected spectra for both the normal and inverted hierarchies for the case wherein neutrinos are Dirac particles, and the energy resolution of the detector is 20 meV - a realistic but perhaps optimistic figure. The solid lines show the normal ordering case, with orange being the C$\nu$B spectrum and blue being the background noise from the tritium beta decay. Similarly the dashed lines are the inverted hierarchy case, with red the C$\nu$B spectrum and green the background noise from the tritium beta decay. The vertical axis describes the capture rate for each energy while the horizontal axis denotes the emitted electron's kinetic energy $K_e$, less the kinetic energy of the endpoint of beta decay for the case of the lightest neutrino being massless $K_0$, which is denoted $K_{\mathrm{end}}^0$ in the body of this work - see (<ref>). The lightest mass neutrino was chosen in both cases to be massless, and the values shown are for the case of a 100g tritium experiment. Plotted using Python.

The expected spectra for both the normal and inverted hierarchies for the case wherein neutrinos are Dirac particles, and the energy resolution of the detector is 0.4 meV - the maximum necessary to fully resolve all C$\nu$B signals. The solid lines show the normal ordering case, with orange being the C$\nu$B spectrum and blue being the background noise from the tritium beta decay. Similarly the dashed lines are the inverted hierarchy case, with red the C$\nu$B spectrum and green the background noise from the tritium beta decay. The vertical axis describes the capture rate for each energy while the horizontal axis denotes the emitted electron's kinetic energy $K_e$, less the kinetic energy of the endpoint of beta decay for the case of the lightest neutrino being massless $K_0$, which is denoted $K_{\mathrm{end}}^0$ in the body of this work - see (<ref>). The lightest mass neutrino was chosen in both cases to be massless, and the values shown are for the case of a 100g tritium experiment. Plotted using Python.

These figures are plotted for two specific values of the energy resolution - 20 meV, a realistic  but perhaps optimistic figure, and 0.4 meV - the threshold resolution necessary to fully resolve all C$\nu$B signals, as we shall see.  The beta decay background near the end point is larger in the normal hierarchy case than that in the inverted hierarchy case since the lightest neutrinos contribute to the endpoint, proportionally to $|U_{ei}|^2$, and $|U_{ei}|^2$ is obviously larger in the normal ordering case, for when $i=1$.

In order to see how small a resolution is needed, we need to consider the capture rate within the energy range $\Delta$ around a C$\nu$B signal's $E^{\mathrm{C\nu B}}\approxeq K_{\mathrm{end}}^0+m_e +E_{\nu_i}$ for both the signal and the background noise:

\begin{align}
\tilde{\Gamma}_i(\Delta)=\int^{E^{{\rm C\nu B},i}_e+\Delta/2}_{E^{{\rm C\nu B},i}_e-\Delta/2}dE_e\frac{d\tilde{\Gamma}_i}{dE_e}(E_e),
\label{OGi} \\
\tilde{\Gamma}_{\beta, i}(\Delta)=\int^{\langle E^{{\rm C\nu B},i}_e \rangle+\Delta/2}_{\langle E^{{\rm C\nu B},i}_e \rangle-\Delta/2}dE_e\frac{d\tilde{\Gamma}_\beta}{dE_e}(E_e).
\label{OGbeta}
\end{align}

This $\langle E^{{\rm C\nu B},i}_e \rangle$ is defined as the average energy for a C$\nu$B signal electron, so $E_{\nu_i}=\sqrt{\braket{p_0}^2+m_{\nu_i}^2}$ with $\braket{p_0}$ taken from \eqref{lowenergy}.

From this, we can define a signal-to-noise ratio:

\begin{equation}
r^i_{\mathrm{SN}}(\Delta) = \frac{\tilde{\Gamma}_i(\Delta)}{\tilde{\Gamma}_{\beta, i}(\Delta)}.
\end{equation}

This ratio is how to tell quantitatively if a C$ν$B signal will be drowned out by the normal beta decay. An $r^i_SN$ value of 1 would mean that there are as many signals in that energy range from beta decay as there are from C$ν$B interactions - this is the tipping point of observation. An $r^i_SN$ that is $≫1$ would mean extreme certainty of observation, and $≪1$ would mean any C$ν$B signal would be completely masked behind the veil of tritium beta decay.

In order to confirm the effects of non-instantaneous decoupling and gravitational clustering, we need to measure the C$ν$B signals to within 1\% precision. In order for this to ever be possible, we would need an $r^i_SN$ of at least 100. As it turns out, the energy resolution needed to achieve $r^i_SN=1$ or $r^i_SN=100$ is almost identical, as $r^i_SN$ is an exponentially rising function of $Δ$.

In Tables \ref{tb:DeltaNH} and \ref{tb:DeltaIH}, we show the energy resolution required to distinguish between the C$ν$B signals and the beta decay background in both the Dirac and Majorana cases for $r^i_SN=1$ and $r^i_SN=100$, considering both the normal and inverted mass orderings.
In both the Dirac and Majorana cases, we find almost the same required $Δ$ - this makes sense, as the position of the C$ν$B peaks are dependent on the neutrino masses, not their fermionic type. In the normal ordering case, the required $Δ$ values for $ν_1$ with $m_ν_1=0 meV$, $ν_2$ with $m_ν_2=8.6 meV$ and $ν_3$ with $m_ν_1=50 meV$ are $0.46 meV, 4.2 meV$ and $19 meV$ respectively.

In the inverted hierarchy case, owing to the mass and mass difference scales, it is difficult to distinguish between the signals for $ν_1$ and $ν_2$ due to their nearly degenerate masses. Because of this, we estimate the required $Δ$ to distinguish the combined signal for $ν_1$ and $ν_2$ from the $β$-decay background, as seen in table \ref{tb:DeltaIH}. The required $Δ$ values for $ν_1$ and $ν_2$ with $m_ν_1=49.3 meV$ and $m_ν_2=50 meV$, and $ν_3$ with $m_ν_3=0 meV$ are $23 meV$ and $0.46 meV$ respectively.

In both cases, if the lightest neutrino is taken to be massive rather, the required energy resolution to detect the lightest neutrino becomes larger and the difficulty in doing so lessened. Thus, the required energy resolution calculated herein is the absolute minimum needed.

\begin{table}[h]
\begin{center}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
NH case& $r^M_{\mathrm{SN}}=1$ & $r^M_{\mathrm{SN}}=100$ & $r^D_{\mathrm{SN}}=1$ & $r^D_{\mathrm{SN}}=100$ \\
\hline
$\nu_1\ (m_{\nu_1}=0\ {\rm meV})$ & $\Delta =0.83\ {\rm meV}$ &  $\Delta =0.46\ {\rm meV}$ & $\Delta =0.83\ {\rm meV}$ & $\Delta =0.46\ {\rm meV}$   \\
$\nu_2\ (m_{\nu_2}=8.6\ {\rm meV})$ & $\Delta =5.3\ {\rm meV}$ &  $\Delta =4.2\ {\rm meV}$ & $\Delta =5.1\ {\rm meV}$ & $\Delta =4.1\ {\rm meV}$   \\
$\nu_3\ (m_{\nu_3}=50\ {\rm meV})$ & $\Delta =21\ {\rm meV}$ &  $\Delta =19\ {\rm meV}$ & $\Delta =21\ {\rm meV}$ & $\Delta =18\ {\rm meV}$   \\
\hline
\end{tabular}
\caption{Table showing the energy resolutions needed to achieve certain signal-to-noise ratios $r_{\mathrm{SN}}$ for both the cases wherein neutrinos are Dirac and Majorana particles, for each massive species in the normal hierarchy. Note that, as expected, the massless species requires the smallest resolution of under 0.5 meV, while the most massive only requires around 20 meV. Taken from [63].}
\label{tb:DeltaNH}
\end{center}
\end{table}

\begin{table}[h]
\begin{center}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
IH case& $r^M_{\mathrm{SN}}=1$ & $r^M_{\mathrm{SN}}=100$ & $r^D_{\mathrm{SN}}=1$ & $r^D_{\mathrm{SN}}=100$ \\
\hline
$\nu_1+\nu_2 $ & $\Delta =29\ {\rm meV}$ &  $\Delta =23\ {\rm meV}$ & $\Delta =28\ {\rm meV}$ & $\Delta =22\ {\rm meV}$   \\
\hspace{-0.1cm}$(m_{\nu_1, \nu_2}=49.3, 50\ {\rm meV})$ \hspace{-0.3cm} & &   &  &    \\
$\nu_3\ (m_{\nu_3}=0\ {\rm meV})$ & $\Delta =0.83\ {\rm meV}$ &  $\Delta =0.46\ {\rm meV}$ & $\Delta =0.83\ {\rm meV}$ & $\Delta =0.46\ {\rm meV}$   \\
\hline
\end{tabular}
\caption{Table showing the energy resolutions needed to achieve certain signal-to-noise ratios $r_{\mathrm{SN}}$ for both the cases wherein neutrinos are Dirac and Majorana particles, for each massive species in the inverted hierarchy. Note that, as expected, the massless species requires the smallest resolution of under 0.5 meV, while the most massive only requires around 20 meV. Taken from [63].}
\label{tb:DeltaIH}
\end{center}
\end{table}

Finally we can consider the energy resolution needed to distinguish between the signals from $ν_1$ and $ν_2$ neutrinos in the inverted hierarchy case. For this purpose, we can define

\begin{align}
r^{12}_{\mathrm{SN}}(\Delta)&=\frac{\int^{E^{{\rm C\nu B},1}_e+\Delta/2}_{E^{{\rm C\nu B},1}_e-\Delta/2}dE_e\frac{d\tilde{\Gamma}_1}{dE_e}(E_e)}{\int^{E^{{\rm C\nu B},1}_e+\Delta/2}_{E^{{\rm C\nu B},1}_e-\Delta/2}dE_e\frac{d\tilde{\Gamma}_2}{dE_e}(E_e)}, \label{r12}
\end{align}

and the inverse would be $r^21_SN(Δ)$. We find numerically that for $r^12_SN(Δ)=100$ we would need $Δ≊0.50$ meV and for $r^21_SN(Δ)=100$ we would need $Δ≊0.42$ meV [63]. Thus, to both distinguish between $ν_1$ and $ν_2$ in the inverted hierarchy, as well as distinguishing between $ν_3$ and the beta decay background, we would need an energy resolution of about $Δ≈0.4$ meV.

\subsection{Beyond the Standard Model}

We have now established what standard physics tells us about the expected capture rate and energy spectra in a capture-by-tritium experiment. Any deviations from this outside of the usual error margins would be a clear sign of new physics.

As discussed in chapter 5, if a large deviation was detected, the leading theories to explain it include sterile neutrinos, lepton asymmetry, neutrino decay and a change to the C$ν$B's temperature.\\

The number density for sterile neutrinos would be given by \eqref{sterilen}, accounting for the remainder of the observed $N_eff$. From the latest calculation in [28], the theoretical value for $N_eff$ is 3.044 when accounting for only 3 flavours, while the observed value from [7] is maximally 3.33, to within a 95\% confidence level. It is important to note that this maximal experimental value depends on the fixing of other parameters in the data, and as such can be increased. The reason why the observed effective number of species might not be closer to 4 can be seen in \eqref{Neff}: if we consider a fourth neutrino species as well, it would be altered to

\begin{equation}
N_{\mathrm{eff}}=\left(\frac{1.40102}{1.39797}\right)^4 \left(4+\frac{\delta\rho_{\nu_e}}{\rho_{\nu_0}}+\frac{\delta\rho_{\nu_\mu}}{\rho_{\nu_0}}+\frac{\delta\rho_{\nu_\tau}}{\rho_{\nu_0}}+\frac{\delta\rho_{\nu_\chi}}{\rho_{\nu_0}}\right),
\end{equation}

where $χ$ is some random symbol for the fourth flavour of neutrino. Thus, if $δρ_ν_χ/ρ_ν_0$ is both negative and fairly large relative to 1, it is possible to achieve an $N_eff$ of 3.33.

Defining $ΔN_eff = N_eff^exp-N_eff^theo ≈0.284$, we can begin to see what the capture rate for these sterile neutrinos might be. Recall that, in order to explain certain excesses in laboratory experiments, the parameters chosen for these sterile neutrinos are a mass between 0.3 eV and 3 eV  (based on $Δm_41^2 ≈(0.1 - 10) eV^2$) and $|U_e4|^2≈|U_μ4|^2≈0.03$. For the sake of calculation, we can consider a 1 eV neutrino species.

According to [30], the clustering effect $δn^c_ν_s$ would actually be the new leading order term, being approximately equal to 10. Putting this all together in \eqref{sterilen} yields

\begin{equation}
n_4 \approxeq \Delta N_{\mathrm{eff}}(1+\delta n^c_{\nu_4})n_0 \approxeq 3.124n_0.
\end{equation}

For a mass this large, the velocity and the effect of non-instantaneous decoupling would be negligible, and to leading order, the capture rate would be approximately

\begin{equation}
\Gamma_4^M = N_T\frac{G_F^2}{\pi}|V_{ud}|^2|U_{e4}|^2\frac{m_\mathrm{^3He}E_e|\bar{p}_e|}{m_\mathrm{^3H}}F(2,E_e)\left(\braket{f_F}^2 + \frac{g_A^2}{g_V^2}\braket{g_{GT}}^2 \right)n_4 \approxeq  3.124\frac{|U_{e4}|^2}{|U_{e1}|^2}\Gamma_1^M \approxeq 0.754 \mathrm{yr}^{-1},\label{sterilerate}
\end{equation}

with the rate in the Dirac case being roughly half of this. While this capture rate appears small, it is larger than those of some active neutrinos seen in the section of numerical results. On top of this, however, is the fact that sterile neutrinos (should they exist) would be far easier to detect, as their large mass would cause their peak in the electron energy spectrum to be much further away from both the beta decay endpoint and those of other active species, circumventing the need for extremely good energy resolution. Atop this, including a fourth massive species in one's model allows for parameters to be varied so as to fit the data. In the model with only one added species, the data from Planck can be reinterpreted and $N_eff$ can be close to 4 [37], making \eqref{sterilerate} even larger, by a factor of around 3.

For the case of a keV-scale sterile neutrino to account for dark matter, the detection prospects are far less likely. Though the characteristic peak would be unmistakable, the capture rate would simply be far too small. As stated in chapter 5 already, in this case, $|U_e4|^2$ would be of the order $10^-11$, constrained to be maximally $<10^-9$ [43]. Even though, owing to their much larger mass, they would cluster extremely effectively, this would only increase their number density at earth by around a factor of 800, $n_4≊800n_0$ and so in the best case scenario, the capture rate would be

\begin{equation}
\Gamma_4^M \approxeq 800 \frac{|U_{e4}|^2}{|U_{e1}|^2}\Gamma_1^M \approxeq 6.4\times 10^{-6} \mathrm{yr}^{-1},
\end{equation}

corresponding to roughly, on average, one detection every 156 thousand years. Clearly this experiment would not be the best way to detect these keV neutrinos, should they exist.\\

Next we can turn to considering the effect a lepton asymmetry would have on the capture rate. If there was some primordial asymmetry between neutrinos and anti-neutrinos, a chemical potential would arise as an attempt to balance it out, entering into the distribution function in the usual way. As we showed in \eqref{asym}, we saw that including this effect leads to additional terms in the number density:

\begin{equation}
n_{\nu} = T_\nu^3\left(\frac{3\zeta(3)}{4\pi^2} +\frac{\xi}{12} +\frac{\ln(2)}{2\pi^2}\xi^2\right),
\end{equation}

with $ξ$ related to the chemical potential by $ξ=μ/T_ν$. Thus, the enhancement factor owing to the chemical potential is given by

\begin{equation}
\frac{n_\nu}{n_0} \approxeq 1 + \frac{\pi^2\xi}{9\zeta(3)} + \frac{2\ln(2)\xi^2}{3\zeta(3)}.
\end{equation}

This factor is actually always there, but in the standard model, $ξ$ would be negligibly small. If we consider other scenarios however, we can see that the chemical potential can indeed affect the capture rate.

For the Dirac case, where only left-helical neutrinos are detected, this value of $ξ$ could be either positive or negative: that is, there could be more neutrinos or anti-neutrinos in the universe today. Thus, for certain values of $ξ$, this factor might actually diminish the capture rate of our relic neutrinos. For example, for some arbitrary value of $ξ=±0.1$, we have

\begin{equation}
n_\nu = n_0\left(1\pm\frac{\pi^2}{90\zeta(3)} + \frac{2\ln(2)}{300\zeta(3)}\right) \approxeq n_0(1\pm 0.0912 + 0.0038),
\end{equation}

changing the capture rate by almost 10\% in either direction. For Majorana neutrinos however, both the left-helical neutrinos and the right-helical neutrinos are detected, and these two components - taking the place of ``neutrinos" and ``anti-neutrinos" - will have opposite signs for their chemical potentials. Thus, when calculating the capture rate and summing over spins, the linear term will cancel out and we are left with only the quadratic contribution. Thus, in this case, the capture rate is always enhanced by a lepton asymmetry. Considering again some arbitrary value of $ξ=±0.1$, we would get $n_ν=n_0(1+0.0038)$, increasing the capture rate around 0.4\%.

Observations have set limits on what $ξ$ can actually take on for its values. The current bound is $-0.091<ξ<0.051$ [66], meaning that this can change the capture rate maximally in the Dirac case by either diminishing it by 7.98\% or increasing it by 4.75\%, and in the Majorana case increasing it by approximately 0.32\%.\\

Next up we may consider the case wherein all neutrinos decay into the lightest species. Of course, it is possible to consider some scenario wherein only some of the more massive species have decayed, but for the sake of the calculation, let us consider the case where in the past almost 14 billion years, all relic neutrinos have decayed into the lightest species. In this case then, the lightest eigenstate of neutrino has 3 times as many neutrinos in it, with the others being unpopulated. At leading order, ignoring all other higher-order terms and effects, this will affect the capture rate purely via the elements of the PMNS matrix. Instead of summing over all 3 species, we simply have one term with the PMNS element corresponding to the lightest eigenstate, but multiplied by 3. That is, $∑_i=1^3|U_ei|^2  → 3|U_eℓ|^2$, with $ℓ$ the lightest mass eigenstate. Owing to unitarity, $∑_i=1^3|U_ei|^2=1$ and so the change to the capture rate would be a factor of $3|U_eℓ|^2$. As we shall see, in any scenario this leads to a smaller chance of detection.

Conversely, the presence of multiple signals with different masses will put a new upper constraint on the lifetime of neutrinos - the present age of the universe.\\

Let us consider first the normal ordering case. In this scenario, at first glance, the capture rate seems to be significantly improved :

\begin{equation}
\Gamma^{M,D}_{C\nu B} = \Gamma^{M,D}_{C\nu B,0}\times 3|U_{e1}|^2 \approxeq 2.043 \  \Gamma^{M,D}_{C\nu B,0},
\end{equation}

essentially doubling the total capture rate. However, recall from the previous section that the lightest neutrino species is also the most difficult to observe owing to the poor energy resolution of cosmic neutrino detectors and the experimental noise of regular beta decay. Thus, even though the capture rate would be doubled, unless the energy resolution is improved significantly or unless the lightest species' mass is actually around 20 meV, none of these captured neutrinos would actually be observed.\\

In the inverted ordering case, where the lightest mass is $m_3$, the diminishment of the capture rate is far more apparent. In this case, we have

\begin{equation}
\Gamma^{M,D}_{C\nu B} = \Gamma^{M,D}_{C\nu B,0}\times 3|U_{e3}|^2 \approxeq 0.0666 \  \Gamma^{M,D}_{C\nu B,0},
\end{equation}

a diminishment of the capture rate by over 93\%. Atop this, the added effect of the lightest species' being difficult to detect owing to the noise of beta decay and poor energy resolution leaves us with an essentially invisible signal.\\

Finally, we may consider a change to the C$ν$B's temperature and how this might affect the capture rate. As we have seen many times now, the number density is related to the effective temperature by $n∝T_ν^3$. Any changes to this temperature would then affect the capture rate. In fact, we have already seen this in the section on the distortive effects of non-instantaneous decoupling: as the temperature was changed from the leading order calculation of $T_ν= 1.9454$ K to the more accurate $T_ν= 1.9496$ K, the number density (and thus, the capture rate) was increased by around 0.65\%, seen in \eqref{newnumdens}.

A more striking change may come about if there was some ``dark radiation" (an exotic, yet undetected species of matter that was relativistic at the epoch considered) which annihilated between the times of neutrino and photon decoupling. In this case, the presence of the dark radiation would change the effective degrees of freedom and, as we saw at the end of chapter 5, would change the ratio between the neutrino effective temperature and the CMB's temperature today:

\begin{equation}
T_{\nu,0}\propto T_{\gamma,0}\left(\frac{g_*(T_0)}{g_*(T_{dec})}\right)^{1/3}=T_{\gamma,0}\left(\frac{2}{11/2+\Delta g_*}\right)^{1/3}.
\end{equation}

Thus, the presence of any more species would decrease the temperature, and in turn decrease the number density and capture rate by a factor of

\begin{equation}
n_\nu = n_0\left(\frac{T^{\mathrm{new}}_\nu}{T_\nu^{\mathrm{std}}}\right)^3 = \frac{22}{22+\Delta g_*}n_0.
\end{equation}

For some numerical examples, the inclusion of another fermionic spin-$1/2$ particle (with its antiparticle) would reduce the capture rate by about 13.7\% \footnote{This hypothetical particle does not interact strongly, and so does not have colour degrees of freedom.}, while the inclusion of 3 of these would reduce the capture rate by approximately 32.3\%.\\

One final point must be made here. The calculations leading up to this section are what we \itshape expect \normalfont to be detected. In this section, then, we provided some explanations, should the capture rate deviate from our predictions. Unfortunately though, should some deviation be found, it seems extremely difficult, for the most part, to differentiate between these models. As we saw in our numerical calculations, the capture rate varies considerably based on factors like whether neutrinos are Dirac or Majorana particles, as well as their mass. Other standard factors like gravitational clustering are only approximate, and rely on numerical simulations, whose error margins can certainly be argued to be large. On top of this, the elements of the PMNS matrix have current error margins of around 10\% at the 3$σ$ level. The question is, how then, given this miasma of factors and parameters, are we expected to eke out information and conclusions from data collected?

While the detection of the C$ν$B on tritium may not instantly answer all of our questions about the Dirac versus Majorana debacle, the questions of sterile neutrinos and lepton asymmetry and other such mysteries, it will certainly provide many more answers than we currently have available. For example, the observed peaks in the energy spectrum of the emitted electrons will certainly provide tight bounds on the neutrino masses. The capture rates - applying Occam's razor - should provide a sign as to whether they are Majorana or Dirac neutrinos, owing to the factor of 2 difference for non-relativistic species. The presence of sterile neutrinos would have the smoking gun of a characteristic peak in the electrons' spectra, and neutrino decay would be seen by the presence of only a single peak. In conclusion, while competing theories of non-standard physics may arise to assuage the tension if the capture rates are not as expected, varying parameters to fix their models to the data, this first step into observing the C$ν$B will already provide some interesting, if not concrete, information.

\newpage
\section{Lensing as a Multipurpose Tool}

Gravitational lensing is a natural corollary of general relativity, and has been experimentally confirmed since 1919 [67]. Fortunately, it is more than just a consequence - it is a useful tool we can use to our advantage in our quest to learn about neutrinos and our universe. We begin by covering the basics of lensing, then showing how we can use it to massively increase the neutrino influx from a supernova event. After that we show how the C$ν$B can actually be used to observe the evolution of the ``lens", and in the next chapter we will show how lensed supernova neutrinos might help us uncover their masses.

\subsection{A Brief Introduction to Lensing}

To begin our study of gravitational lensing, let us consider a spacetime with a single object of mass $M$ surrounded by vacuum. We can describe this spacetime using the Schwarzschild metric or a Minkowski metric perturbed by a Newtonian potential:

\begin{equation}
\end{equation}

wherein $ϕ=-GM/r$. Recall that this is also the weak field limit of the Schwarzschild metric. Now, since light travels on null geodesics, we can obtain the ``effective" speed of light for when $ds^2=0$. In that case we have

\begin{equation}
(1-2\phi)d\bar{x}^2=(1+2\phi)dt^2 \ \ \ \rightarrow \ \ \ \left|\frac{d\bar{x}}{dt}\right|=\sqrt{\frac{1+2\phi}{1-2\phi}}.
\end{equation}

To linear order then, $|dx̅/dt|≊1+2ϕ$. Then, since the refractive index $n$ is defined as the ratio of the speed of light in vacuum to the actual speed,

\begin{equation}
\end{equation}

Since the potential is always $<0$, we can write $n=1+2|ϕ|$.\\
The same effects seen in usual geometric optics are observed owing to the presence of gravity. Similarly to how light is slowed and its path bent when entering a glass prism, analogously the presence of our mass bends and slows light in its gravitational field. The first of these connections is the deflection angle, $α$, which is related to the gradient of the refractive index [68]:

\begin{equation}
\bar{\alpha}=-\int\bar{\nabla}_{\perp}nd\ell \ = \ 2\int\bar{\nabla}_{\perp}\phi d\ell.
\end{equation}

So, consider our point mass $M$ placed at the origin, and our light will move in the $x-z$ plane, as seen in figure \ref{fig:lensing}. Then the $r$ in our potential is $√(x^2+z^2)$, and we have

\begin{equation}
\bar{\alpha}=2\int\bar{\nabla}_{\perp}\phi d\ell = 2\int\bar{\nabla}_x\frac{-GM}{\sqrt{x^2+z^2}} dz.
\end{equation}

Most of the bending will occur near the point of closest approach to the mass, and so this is the most important region of the path. At the point of closest approach, the direction perpindicular to the path is $x$. So, performing the differentiation and taking the source to be infinitely far away and us (the observer) to be as well in the other direction, we have

\begin{equation}
\bar{\alpha}=2\int_{-\infty}^{\infty}\frac{GM\bar{x}}{(x^2+z^2)^{3/2}}dz = \frac{2GM\bar{x}}{x^2}\frac{z}{\sqrt{x^2+z^2}}\bigg|^\infty_{-\infty}.
\end{equation}

So the magnitude of the angle of deflection will be

\begin{equation}
\alpha=\frac{4GM}{x}\lim_{z\rightarrow\infty}\frac{z}{\sqrt{x^2+z^2}} = \frac{4GM}{x}.\label{alpha}
\end{equation}

This $x$ is measured at the point of closest approach, and so is known as the impact parameter $b$. Thus, the deflection angle comes out to be $α=4GM/b=2r_s/b$, where $r_s$ is the Schwarzschild radius of the lensing object.\\

\begin{figure}
\centering
\includegraphics[width=\linewidth]{Lensing}
\caption{Lensing by a point mass, drawn using [6].}
\label{fig:lensing}
\end{figure}

In the above, we assumed that the source, lens and observer were collinear, but this may not always be the case. Depicted in figure \ref{fig:lensing}, we can relate the angle of the \itshape apparent \normalfont source $θ$ to that of the actual source $β$ by simple geometry:

\begin{equation}
\theta D_S = \beta D_S+\alpha D_{LS},\label{lensinit}
\end{equation}

where $D_S$ is the distance to the source, $D_LS$ the distance from the lens to the source, and $α$ our deflection angle as before. In the above, all angles are taken to be small enough such that the $sinx≊x$ approximation can be used. Thus, for the case of a single point mass, we have $α$ from \eqref{alpha} and

\begin{equation}
\theta-\beta=\frac{D_{LS}}{D_S}\frac{4GM}{b}.
\end{equation}

$b$, the impact parameter, is the distance of closest approach and geometrically (again using the small angle approximation) $b≊D_Lθ$, with $D_L$ the distance to the lens. An important distance is the Einstein radius: the distance from the point mass that the apparent source appears, when the source, lens and observer form a syzygy. In this case, $β=0$ and we have

\begin{equation}
\theta=\frac{D_{LS}}{D_S}\frac{4GM}{D_L\theta} \ \ \rightarrow \ \ \theta=\sqrt{\frac{4GMD_{LS}}{D_SD_L}}.
\end{equation}

The Einstein ``radius", though actually an angle, is commonly symbolised by $θ_E$ and is given by (under the small angle approximation)

\begin{equation}
\theta_E=\theta-0=\sqrt{\frac{4GMD_{LS}}{D_SD_L}}=\sqrt{\frac{4GMy}{D_L(1+y)}},\label{einsteinring}
\end{equation}

where $y=D_LS/D_L$. This angle, multiplied by $D_L$, gives the radius of the ``ring" formed by the source when lensed. This is because in the case of a syzygy, the light from the source is lensed towards the observer symmetrically around the lens, forming a circular image. Lensed light from sources that do not lie exactly on the observer-lens line of sight do not form rings, and instead form arcs or arclets [68]. The Einstein ring is the most extreme case of a lensing event.\\

One extremely useful consequence of lensing is that the source is ``magnified" - that is, more light (or any particles) reaches us, as even though the surface brightness is conserved [69], the apparent solid angle is increased. The magnification $μ$, when the lens is circularly symmetric (which is the case for a point mass), is defined as

\begin{equation}
\mu = \frac{\theta}{\beta}\frac{d\theta}{d\beta}\label{mag}.
\end{equation}

To this end, for our point mass lens, we can use \eqref{alpha}, \eqref{lensinit}, \eqref{einsteinring} and $b≊D_Lθ$ to see that

\begin{equation}
\beta = \theta-\frac{\theta_E^2}{\theta}.\label{lens}
\end{equation}

Solving this yields $θ=1/2(β±√(β^2+4θ_E^2))$. For a point mass, any source is always lensed twice, creating two images or arc(lets)\footnote{In the case of an Einstein ring, the two arcs are semicircles and thus appear as one ring.}(for more complex lenses, this is not the case). These two images occur at what are commonly denoted $θ_±$, and unless $θ_±=θ_E$, one image is always outside the ring and one inside. Thus, using \eqref{mag} and \eqref{lens}, the full magnification is given by the sum of the magnification of the two images:

\begin{equation}
\mu = \mu_++\mu_- = \left(1-\left(\frac{\theta_E}{\theta_+}\right)^4\right)^{-1} + \left(1-\left(\frac{\theta_E}{\theta_-}\right)^4\right)^{-1} = \frac{\frac{\beta^2}{\theta_E^2}+2}{\frac{\beta}{\theta_E}\sqrt{\frac{\beta^2}{\theta_E^2}+4}}.\label{mag2}
\end{equation}

%This increased intensity of incoming particles is owing purely to the increase in solid angle we view of the source, and the source's surface brightness is unaltered. To see this, consider a lens of mass $M$ which is no longer a point, but is circularly symmetric along the axis of the source-observer line, as seen in fig ??. If we consider the case without the lensing object, then at the point of the object the source's radiation would have some impact parameter $b_0$, while in the presence of the lens it becomes $b$.
For a point mass this takes on a maximum value as $β→0$, in which case $μ→θ_E/β$. Physically, $β$ cannot be smaller than the radius of the source over the distance to it, so $μ_max=θ_ED_S/R_S$.\\

Owing to the presence of the lens, the light from the source must take a different, longer path to reach us, resulting in a change to the travel time known as the Shapiro time delay. Since the light from the source can travel along two different paths (or more, for more complex mass distributions) there is another time delay between the two detections, $Δt = t_–t_+$. The light from the innermost apparent source arrives last, as it has the longest path to travel, while that coming from the outside of the ring arrives first. Thus, $Δt$ is always positive. These time delays will be calculated for supernova neutrinos below.\\

Before continuing, for lensing objects such as galaxies, the point-mass model is perhaps a little too far from reality. For more complex lensing objects that have some structure, there are other models available. We will discuss one here: the singular isothermal sphere model.

In this model, galaxies are modelled with their constituents behaving like particles in an ideal gas - an approximation that is quite fair, as most astrophysical objects do not regularly interact. Thus, the equation of state of these constituents is of course the usual $P=NT/V$, where the Boltzmann constant $k_B=1$ and we can rewrite this as

\begin{equation}
P=\frac{\rho T}{m},\label{idealgas}
\end{equation}

with $ρ$ and $m$ the mass density and mass of the constituents, respectively\footnote{Of course in reality not every object in a galaxy has the same mass, and this mass is usually taken to be an average stellar mass.}. In thermal equilibirium, the temperature is related to the velocity dispersion by equating the thermal and kinetic energies:

\begin{equation}
\frac{3}{2}T = \frac{1}{2}m\sigma_v^2.\label{energies}
\end{equation}

The word ``isothermal" refers to the fact that the temperature is taken to be constant throughout the galaxy, and so in turn is the velocity dispersion. Of course, a more general model could be used for a radially-dependent temperature, but since most of the galaxy consists of interstellar gas which is approximately isothermal, it is not a terrible approximation.

Next, from hydrostatic equlibirium (wherein the force of thermal pressure balances that of gravity) we have that

\begin{equation}
\frac{1}{\rho}\frac{dP}{dr} = -\frac{GM}{r^2}.
\end{equation}

So, since $M=4π/3ρr^3$, we have that $dP/dr =-4πGρ^2r/3$ and so $P=-2πGρ^2r^2/3$. And so, using this along with \eqref{idealgas} and \eqref{energies}, we have

\begin{equation}
\rho = \frac{\sigma_v^2}{2\pi Gr^2}.
\end{equation}

In more general lenses that are not point masses, the surface mass density is important, as this essentially describes the ``lens" nature of the object, and the surface mass density distribution is all that is needed to derive the characteristics of a lens. It is given by

\begin{equation}
\Sigma(r) = \int \rho(r,z)dz,
\end{equation}

where $r$ is the radial distance from the centre of the lens and $z$ the distance along the observer-lens axis. So, in the case of a singular isothermal sphere (SIS), we have

\begin{equation}
\Sigma(r) = \int \frac{\sigma_v^2}{2\pi Gr^2} dz = \int_0^\pi \frac{\sigma_v^2}{2\pi Gr} d\theta = \frac{\sigma_v^2}{2Gr},
\end{equation}

where we used that $dz≊rdθ$. To obtain the deflection angle, we need to sum over all deflections from all the masses in the plane. This gives a two dimensional integral:

\begin{equation}
\alpha(r) = 4G \int d^2r' \frac{(\bar{r}-\bar{r}')\Sigma(r')}{|\bar{r}-\bar{r}'|^2}.
\end{equation}

To solve this, we put $r'$ into radial coordinates, such that $d^2r' = r'dr'dθ'$ and then use circular symmetry to integrate over the angular coordinate. Next, we use a change of variables from $r-r'$ to $r'$, similar to integrals in electrostatics. Then, we are left with

\begin{equation}
\alpha(r) = \frac{8\pi G}{r} \int_0^r \Sigma(r')r'dr'.
\end{equation}

This can be written as $α(r) = 4GM(r)/r$, which is remniscient of \eqref{alpha}, except instead of $M$ being a known value, it is now given by $M(r) = 2π∫_0^r Σ(r')r'dr'.$ Thus, for our SIS, we have

\begin{equation}
\alpha(r) = \frac{4GM(r)}{r} = \frac{4G}{r}\cdot2\pi\int_0^r \frac{\sigma_v^2}{2Gr'}r'dr' = 4\pi \sigma_v^2.
\end{equation}

This deflection angle can then be inserted into \eqref{lensinit} to obtain the Einstein radius, and that in turn can be used with the lens equation to calculate the magnification and time delays.

As already stated, there are many models used to describe lenses. As long as the surface distribution $Σ(r)$ is known, all other quantities can be derived from there.

\subsection{Lensing of Massive Particles}

While the previous subsection served as our springboard, for the case of massive particles we require the full Schwarzschild metric to describe our spacetime.

The Schwarzschild metric is given by

\begin{equation}
ds^2=-d\tau^2=-(1-\frac{2M}{r})dt^2+(1-\frac{2M}{r})^{-1}dr^2+r^2d\Omega^2, \label{SC}
\end{equation}

where $dΩ^2=r^2dθ^2+r^2sin^2θdϕ^2$, and for this calculation we have taken both $c$ and $G$ to be unity.\\
Let us consider some particle moving in this spacetime around our central mass. We can always align our coordinate system such that the particle's orbit is at $θ=π/2$, so that $dθ=0$ and $sin^2θ=1$.

Dividing \eqref{SC} by $dτ^2$ gives us

\begin{equation}
g_{\mu\nu}\frac{dx^\mu}{d\tau}\frac{dx^\nu}{d\tau} = -(1-\frac{2M}{r})\left(\frac{dt}{d\tau}\right)^2+ (1-\frac{2M}{r})^{-1}\left(\frac{dr}{d\tau}\right)^2-r^2\left(\frac{d\phi}{d\tau}\right)^2.\label{r}
\end{equation}

Now, requiring that the line element is invariant, we can use calculus of variations:

\begin{equation}
\delta s = 0 \ \ \rightarrow \ \ \delta\int ds = 0 \ \ \rightarrow \ \ \delta\int\sqrt{ds^2} = 0 \ \ \rightarrow \delta\int\sqrt{g_{\mu\nu}\frac{dx^\mu}{d\tau}\frac{dx^\nu}{d\tau}}d\tau = 0.
\end{equation}

Thus, in order for this to be invariant, the integrand must obey the Euler-Lagrange equations. However, if $ℒ$ is conserved, so is $ℒ^1/2$. Thus, instead of considering the above, our integrand can simply be taken to be

\begin{equation}
\mathcal{L} = -(1-\frac{2M}{r})\left(\frac{dt}{d\tau}\right)^2+ (1-\frac{2M}{r})^{-1}\left(\frac{dr}{d\tau}\right)^2-r^2\left(\frac{d\phi}{d\tau}\right)^2.
\end{equation}

Applying the E-L equations, $d/dτ∂ℒ/∂ẋ=∂ℒ/∂x$ for the cases $x=t$ and $ϕ$ yield

\begin{equation}
\frac{d}{d\tau}\left[\left(1-\frac{2M}{r}\right)\frac{dt}{d\tau}\right]=0 \ \ \rightarrow \ \ \left(1-\frac{2M}{r}\right)\frac{dt}{d\tau}=E,
\end{equation}
\begin{equation}
\frac{d}{d\tau}\left[r^2\frac{d\phi}{d\tau}\right] = 0 \ \ \rightarrow \ \ r^2\frac{d\phi}{d\tau}=J,
\end{equation}

where $E$ and $J$ are constants of integration and represent conserved quantities (the former being energy and the latter angular momentum per unit mass). Rearranging \eqref{r}, we have that

\begin{equation}
\left(\frac{dr}{d\tau}\right)^2 =\left(1-\frac{2M}{r}\right)^2\left(\frac{dt}{d\tau}\right)^2 -\left(1-\frac{2M}{r}\right)-r^2\left(1-\frac{2M}{r}\right)\left(\frac{d\phi}{d\tau}\right)^2,
\end{equation}

\begin{equation}
\rightarrow \ \ \left(\frac{dr}{d\tau}\right)^2 + \left(1-\frac{2M}{r}\right)\left(1+\frac{J^2}{r^2}\right)=E^2.
\end{equation}

Our next act of manipulation involves $u=1/r$ and $λ=E^2-1/J^2$. Since $dr/dτ=-1/u^2du/dτ$ and $du/dτ=du/dϕdϕ/dτ=du/dϕJu^2$, we can rewrite our previous result as

\begin{equation}
\left(\frac{du}{d\phi}\right)^2=\lambda+\frac{2Mu}{J^2}-u^2+2Mu^3.
\end{equation}

Now, we can write this as

\begin{equation}
\left(\frac{du}{d\phi}\right)^2=2M(u-u_0)(u-u_1)(u-u_2),\label{uphi}
\end{equation}

where $u_i$ are the roots of $λ+2Mu/J^2-u^2+2Mu^3$. Now, though this is the fully general relativistic equation, in the Newtonian limit, the $u^3$ term would disappear, and so we can obtain the Newtonian roots of $λ+2Mu/J^2-u^2$ as

\begin{equation}
u_{1,2}=\frac{M}{b^2v_0^2}\left(1\pm\sqrt{1+\frac{\lambda b^4v_0^4}{M^2}}\right),
\end{equation}

where we have used that the angular momentum per unit mass $J$ at the point of closest approach to the lensing object is equal to the impact parameter $b$ times the velocity of the particle at infinity, $v_0$.

Now these $u_1$ and $u_2$ can be used in our relativistic case approximately, and we can rewrite $u_0$ in terms of them as well. We have that

\begin{equation}
2Mu^3-u^2+\frac{2M}{J^2}u+\lambda = 2M(u-u_0)(u-u_1)(u-u_2),
\end{equation}

and expanding the right hand side and equating coefficients of powers of $u$ by uniqueness of power series we get from the coefficient of $u^2$ that

\begin{equation}
\end{equation}

Going back to \eqref{uphi}, we have that

\begin{equation}
\phi-\phi_0=2\int_{0}^{u_1}\frac{du}{\sqrt{2M(u-u_0)(u-u_1)(u-u_2)}},
\end{equation}

where the factor of 2 in front of the integral comes from the fact that we have the path from the source to the point of closest approach, and that from the point of closest approach to the observer. Using our relationship between the $u_i$, and taking the initial angle to be $-π$ (that is, the source is approximately infinitely far away) we have

\begin{equation}
\phi=2\int_{0}^{u_1}\frac{du}{\sqrt{(u_1-u)(u-u_2)(1-2M(u+u_1+u_2))}}-\pi.
\end{equation}

We now make the change of variables $u=u_2+(u_1-u_2)sin^2x$, such that at $u=0$ we have $x=x_*=arcsin√(u_2/u_2-u_1)$ and at $u=u_1$ we have $x=π/2$. Using that $du=2(u_1-u_2)sinxcosxdx$, we have

\begin{equation}
\phi + \pi = 4\int_{x_*}^{\frac{\pi}{2}}\frac{dx}{\sqrt{1-2M(u_1+u_2+u_1\sin^2x+u_2\cos^2x)}}.
\end{equation}

Now, $2Mu_i$ goes as $M^2/b^2v_0^2$, which in all cases of interest is a small quantity. Thus, we will take $M/bv_0$ to be our small quantity and perform perturbation theory. \\
For example, in the solar system case, considering a solar-mass star and a particle with relativistic speed at an impact parameter of characteristic distance around 1 AU, this quantity is of the order $10^-8$. For a super massive black hole at the centre of our galaxy, using a galactic characteristic distance of around 8 kpc, this quantity is as small as $10^-11$.

So, using the binomial approximation this integral becomes

\begin{equation}
\phi+ \pi = 4\int_{x_*}^{\frac{\pi}{2}}dx(1+M(u_1+u_2+u_1\sin^2x+u_2\cos^2x) +\frac{3M^2}{2}(u_1+u_2+u_1\sin^2x+u_2\cos^2x)^2+...).
\end{equation}

Performing this integration, and keeping terms only linear in $M/bv_0$, we have

\begin{equation}
\phi+ \pi \approxeq 4\left[\frac{\pi}{2}-x_*+M(u_1+u_2)\left(\frac{\pi}{2}-x_*\right) +\frac{M}{2}u_1\left(\frac{\pi}{2}-x_*+\sin x_*\cos x_*\right) +\frac{M}{2}u_2\left(\frac{\pi}{2}-x_*-\sin x_*\cos x_*\right)\right].
\end{equation}

Now, we know that $sinx_*=√(u_2/u_2-u_1)$. Plugging in our expressions for $u_i$, and once again using the binomial expansion owing to our small parameter, we have

\begin{equation}
\sin x_*=\sqrt{\frac{\sqrt{1+\frac{\lambda b^4v_0^4}{M^2}}-1}{2\sqrt{1+\frac{\lambda b^4v_0^4}{M^2}}}} \approxeq \frac{1}{\sqrt{2}}\left(1-\frac{M}{2\sqrt{\lambda}b^2v_0^2}\right).
\end{equation}

Next we use some trigonometric identities to find $x_*$. First, we see that

\begin{equation}
\sin\left(\frac{x_*}{2}-\frac{\pi}{8}\right)\cos\left(\frac{x_*}{2}+\frac{\pi}{8}\right) = \frac{1}{2}\sin x_* -\frac{1}{2}\sin\frac{\pi}{4} = -\frac{M}{4\sqrt{2\lambda}b^2v_0^2}
\end{equation}

and so if we take $x_*$ close to $π/4$, using the small angle approximation for the sine function we obtain $x_*≊π/4-M/2√(λ)b^2v_0^2$. Using this, as well as $cosx_*≊1/√(2)(1+M/2√(λ)b^2v_0^2)$, $λ= 1/b^2$ (since $J=bv_0$ and $E^2-1=v_0^2$), and that $u_1+u_2=2M/b^2v_0^2$, we have to linear order in $M$

\begin{equation}
\phi + \pi \approxeq 2\pi - \pi + \frac{2M}{bv_0^2}+\frac{2M}{b} + ...,
\end{equation}

\begin{equation}
\rightarrow \ \ \phi = \frac{2M}{bv_0^2}(1+v_0^2) +\mathcal{O}\left(\left(\frac{M}{bv_0}\right)^2\right). \label{massdeflec}
\end{equation}

This deflection angle for massive particles reduces exactly to that of massless particles when $v_0^2=1$, equation \eqref{alpha}. Similarly, in the Newtonian limit when $v_0^2≪1$, we obtain the Newtonian deflection angle: $ϕ=2M/bv_0^2$.

Plugging in \eqref{massdeflec} as our $α$, we can get the Einstein radius for massive particles as

\begin{equation}
\theta_E = \sqrt{\frac{2M(1+v_0^2)D_{LS}}{v_0^2D_SD_L}} = \sqrt{\frac{2My(1+v_0^2)}{v_0^2D_L(1+y)}},
\end{equation}

which can be used to calculate the magnification in \eqref{mag}.

Though this derivation was done using a point mass, the result for the SIS model has a similar effect from the mass [70], with a factor of $1+v_0^2/v_0^2$ entering the deflection angle.

\subsection{Lensing of Supernova Neutrinos}

Though we have the cosmic neutrino background, there is no constant source of high energy neutrinos that we may study. On top of this, owing to their very high velocities, they will not cluster gravitationally to assist us with detection. One intense source of these high energy neutrinos seems to be supernovae, but supernovae that are close enough to be of use are rare, and as we shall see, those that are lensed are even less frequent. Adding to our troubles are that these intense bursts have a pulse time of only around 10 seconds [71].\\

Despite these hinderences, with some luck, a supernova event lensed on its way to earth can provide some interesting information, and has distinct characteristics. Deflection by the sun would cause too small an angular deviation for detection, as current neutrino detectors' angular resolution is not very small. For example, a neutrino skimming our sun would have a deflection angle of $α=4GM/R≊8.48 μrad$, while detectors like Super Kamiokande have angular resolution of the order of degrees [71].

Similarly, lensing by very far away (extragalactic) sources (such as other galaxies) would cause the flux to be diminished too greatly to be compensated by the amplification owing to magnification. The most suitable lens seems to be the black hole at the centre of the Milky Way, as it is extremely massive and relatively close.\\

Though any lensing would provide a wealth of information, for our purposes here we will focus on the case of a supernova occurring on the opposite side of our galaxy. The two main characteristic effects of lensing are the amplification of the flux and a time dispersion of their arrival. The usual characteristic behaviour of a supernova can be described as a time-dependent luminosity, modelled as [72]

\begin{equation}
\end{equation}

where $E_B$ is the energy released by the supernova, done so over the time period $t_P$. The factor of 6 arises to account for the 3 flavours of neutrinos and their anti-particles. The exponential decay following the initial pulse has an e-fold time equal to that of the initial pulse, and using $t_P≈3s$, this model is consistent with observations. Thus, the effects of magnification and time dispersion will be clearly seen if the detected luminosity pattern differs significantly from our model.\\

As shown above in \eqref{mag}, the magnification goes as

\begin{equation}
\mu = \frac{\frac{\beta^2}{\theta_E^2}+2}{\frac{\beta}{\theta_E}\sqrt{\frac{\beta^2}{\theta_E^2}+4}},
\end{equation}

going as $θ_E/β$ at very small $β$ values, and with its maximal value $μ_max=θ_ED_S/R_S$. From \eqref{einsteinring}, this gives

\begin{equation}
\mu_{max}=\sqrt{\frac{4GMyD_S^2}{R_S^2D_L(1+y)}} =\frac{1}{R_S}\sqrt{4GMD_Ly(1+y)}.
\end{equation}

As a numerical estimate, for a usual source around 20 km in diameter, and using that the supermassive black hole at the centre of our galaxy is around $3.61×10^6$ solar masses and a distance approximately 8 kpc from us [71], we get $μ_max=2.3×10^11√(y(1+y))$. Thus, varying $y$ from very small values (around 0.01) to 2 (approximately the edge of the Milky Way) only changes $μ_max$ by an order of magnitude. From this, we can see that a perfect supernova-black hole-earth syzygy would result in an absolutely massive amplification of the number of neutrinos reaching earth.\\

Next, we need to consider the time delays of these neutrinos. For this we once again use the Schwarzschild metric. Of course, this spacetime implies that the only object in the spacetime is the black hole, but for the case of neutrinos this is a valid approximation, as the Milky Way is extremely ``transparent" to neutrinos: the mean free path of neutrinos is much larger than the size of the galaxy.\\
If the mean free path is given by
\begin{equation}
\ell_{mfp}=\frac{1}{\sigma n},
\end{equation}

with $σ$ of the order $σ̅≈10^-45cm^2$ as seen in \eqref{charcrosssec} and $n$ for the interstellar medium at its densest in molecular clouds is around $10^6cm^-3$ [73], we have an $ℓ_mfp$ of around $10^20$pc, while the Milky Way's diameter is only around 30kpc.\\

So, using our Schwarzschild metric, we can consider the geodesic equations for $t,r$ and $ϕ$. Recall that $g_00=-(1-2GM/r)$ and $g_11=(1-2GM/r)^-1$. Writing them suggestively, we have:

\begin{equation}
\frac{d^2t}{d\tau^2}+\frac{\partial_rg_{00}}{g_{00}}\frac{dr}{d\tau}\frac{dt}{d\tau}=0,
\end{equation}
\begin{equation}
\frac{d^2\phi}{d\tau^2}+\frac{2}{r}\frac{d\phi}{d\tau}\frac{dr}{d\tau}=0,
\end{equation}
\begin{equation}
\frac{d^2r}{d\tau^2}+\frac{\partial_rg_{11}}{2g_{11}}\left(\frac{dr}{d\tau}\right)^2 -\frac{r}{g_{11}}\left(\frac{d\phi}{d\tau}\right)^2 -\frac{\partial_rg_{00}}{2g_{11}}\left(\frac{dt}{d\tau}\right)^2 =0.\label{rgeo}
\end{equation}

These first two equations can be written as complete derivatives, leading to conserved quantities. Dividing the $t$ equation by $dt/dτ$, and using chain rule that $d/drdr/dτ=d/dτ$, we have:

\begin{equation}
\frac{d}{d\tau}\left[\ln\left(\frac{dt}{d\tau}g_{00}\right)\right]=0 \ \ \rightarrow \ \ \frac{dt}{d\tau} = \frac{C}{g_{00}}.\label{tgeo}
\end{equation}

We can use $C$ to define our parameter $τ$, and we can choose $C=1$. Performing similar manipulation, for the $ϕ$ geodesic equation we have

\begin{equation}
\frac{d}{d\tau}\left[\ln\left(\frac{d\phi}{d\tau}r^2\right)\right]=0 \ \ \rightarrow \ \ \frac{d\phi}{d\tau}=\frac{J}{r^2},
\end{equation}

where $J$ is another conserved quantity, and is once again the angular momentum per unit mass. Plugging these in to \eqref{rgeo}, we have

\begin{equation}
\frac{d^2r}{d\tau^2}+\frac{\partial_rg_{11}}{2g_{11}}\left(\frac{dr}{d\tau}\right)^2 -\frac{J^2}{g_{11}r^3}-\frac{\partial_rg_{00}}{2g_{11}g_{00}^2}=0.
\end{equation}

Multiplying by $2g_11dr/dτ$ and again using chain rule, we have that

\begin{equation}
\frac{d}{d\tau}\left[g_{11}\left(\frac{dr}{d\tau}\right)^2+\frac{J^2}{r^2}+\frac{1}{g_{00}}\right] = 0 \ \ \rightarrow \ \ \frac{dr}{d\tau}=\sqrt{-\frac{E}{g_{11}}-\frac{J^2}{r^2g_{11}}-\frac{1}{g_{00}g_{11}}},\label{Eeqn}
\end{equation}

where $E$ is another constant of integration, and is the conserved quantity related to the energy of the particle. To see this, we consider the line element of our metric (recalling that $θ=π/2$ and that $dθ=0$):
\begin{equation}
ds^2=g_{00}dt^2+g_{11}dr^2+r^2d\phi^2 \ \ \rightarrow \ \ ds^2=d\tau^2\left(\frac{1}{g_{00}}-E-\frac{J^2}{r^2}-\frac{1}{g_{00}}+\frac{J^2}{r^2}\right)=-Ed\tau^2.
\end{equation}

Thus, for massless particles with $ds^2=0$, $E=0$ and for massive particles we have $E>0$ as they are timelike.

Using \eqref{tgeo} and \eqref{Eeqn}, we have

\begin{equation}
\frac{dt}{dr}=\sqrt{\frac{g_{11}}{g_{00}^2\left(-E-\frac{J^2}{r^2}-\frac{1}{g_{00}}\right)}}.\label{tr}
\end{equation}

If we consider \eqref{tr} at infinity, $g_00$ and $g_11$ go to -1 and 1 respectively and we have

\begin{equation}
\left(\frac{dr}{dt}\right)^2\big|_\infty=-E+1,
\end{equation}

and so $E$ is related to the particle's velocity at infinity $v_0$ by $E=1-v_0^2$.\\
Now, to calculate the time delay we need to integrate \eqref{tr}. For our case, we are considering extremely relativistic neutrinos, and so $E≊0$. Similarly, for our case of interest, $GM/r$ is a small quantity. Numerically, for our central black hole we have $GM≈1.73×10^-7$pc [71], and $r$ on the order of kpc, so this is a very small quantity. Thus, it is reasonable to consider the time delay to linear order in $GM/r$. In this case, $g_00=-1+2GM/r$ and $g_11≊1+2GM/r$. Then we have

\begin{equation}
\frac{dt}{dr}\approxeq-\left(1+\frac{3GM}{r}\right)\left(1+\frac{2GM}{r}-\frac{J^2}{r^2}\right)^{-1/2}.
\end{equation}

Next, if the particle's closest approach is at $r_0$, at $r_0$, $dr/dt=0$. So, at $r_0$ we have that $J^2=r_0^2(1+2GM/r_0)$. But $J$ is constant, so this is true for all $r$. Plugging this in, we have

\begin{equation}
\frac{dt}{dr}\approxeq -\left(1+\frac{3GM}{r}\right)\left(1-\frac{r_0^2}{r^2} +2GM(\frac{1}{r}-\frac{r_0}{r^2})\right)^{-1/2}
\end{equation}

\begin{equation}
=-\frac{1}{\sqrt{1-\frac{r_0^2}{r^2}}}\left(1+\frac{3GM}{r}\right)\left(1+\frac{2GM}{r+r_0}\right)^{-1/2} = -\frac{1}{\sqrt{1-\frac{r_0^2}{r^2}}}\left(1+\frac{GM(2r+3r_0)}{r(r+r_0)}\right),
\end{equation}

and so, integrating from the source $r_i$ to the point of closest approach $r_0$, and then from $r_0$ to the observer $r_f$, we have

\begin{equation}
t=\int_{r_0}^{r_i}\frac{dr}{\sqrt{1-\frac{r_0^2}{r^2}}}\left(1+\frac{GM(2r+3r_0)}{r(r+r_0)}\right) + \int_{r_0}^{r_f}\frac{dr}{\sqrt{1-\frac{r_0^2}{r^2}}}\left(1+\frac{GM(2r+3r_0)}{r(r+r_0)}\right).
\end{equation}

Luckily, this integral can be performed analytically. The first term in each integral gives what is known as the ``geometric" contribution: i.e, the time taken for the particle to trace out its path. Though called geometric, it of course still comes to be owing to gravity, as without the lensing object, the path would simply be a straight line. Integrating this term gives

\begin{equation}
\end{equation}

the standard Pythagorean distance. The second term gives the ``gravitational" contribution (despite both terms owing to gravity) and is the gravitational time dilation owing to the lensing mass. Integrating gives

\begin{equation}
t_{grav}=GM\left[\sqrt{\frac{r_i-r_0}{r_i+r_0}}+\sqrt{\frac{r_f-r_0}{r_f+r_0}} + 2\ln\left(\frac{r_i+\sqrt{r_i^2-r_0^2}}{r_0}\right) + 2\ln\left(\frac{r_f+\sqrt{r_f^2-r_0^2}}{r_0}\right)\right].
\end{equation}

Defining $η=βD_S$, we have $r_i=√(D_LS^2+η^2)=D_L√(y^2+β^2(1+y)^2)$ and $r_f=D_L$, where as before $y=D_LS/D_L$. $r_0$ takes on two values, for the two images, $r_0=D_L|θ_±|$. Thus, $Δt$, the time dispersion of the signals, is given by $t_–t_+$, which is $t=t_geo+t_grav$ with $θ_±$ substituted in from solving \eqref{lens}:

\begin{equation}
\Delta t \approxeq 2GM\left[2\frac{\beta}{\theta_E}\sqrt{1+\frac{\beta^2}{4\theta^2_E}} +\ln\left(\frac{1+\frac{\beta^2}{2\theta^2_E}-\frac{\beta}{\theta_E}\sqrt{1+\frac{\beta^2}{4\theta^2_E}}}{1+\frac{\beta^2}{2\theta^2_E}+\frac{\beta}{\theta_E}\sqrt{1+\frac{\beta^2}{4\theta^2_E}}}\right)\right].
\end{equation}