\documentclass{article}
\usepackage{amsmath} % for \text command
\usepackage{fullpage}
\usepackage{tikz}
\usepackage{hyperref}
\usepackage{amsfonts}
\newcommand*\diff{\mathrm{d}}
\usetikzlibrary{shapes}
\author{Anthony LaTorre}
\date{\today}
\title{Searching for Dark Matter with the Sudbury Neutrino Observatory}
\begin{document}
\maketitle
\section{Introduction}
\section{Estimating the Event rate in the SNO detector}
The event rate of self destructing dark matter events, $R$, in the SNO detector is given by first integrating over the detector.
\begin{equation}
R = \int_\mathrm{SNO} \mathrm{d}^3r \, R(r) 
\end{equation}
Next, we integrate over the earth where the dark matter annihilates:
\begin{equation}
R = \int_\mathrm{SNO} \mathrm{d}^3r \, \int_{r'} \mathrm{d}^3r' R(r') \mathrm{P}(\text{detect at r} | \text{DM scatters at r'})
\end{equation}
where we have assumed above that the dark matter annihilates immediately after
scattering. The event rate for scattering at a position $r'$ in the earth is:
\begin{equation}
R(r') = \Phi(r') \eta(r') \sigma(r')
\end{equation}
where $\Phi(r')$ is the flux at $r'$, $\eta(r')$ is the number density of
scatterers at $r'$, and $\sigma$ is the cross section for the dark matter to
scatter and annihilate. In general this will be a sum over the elemental
composition of the earth at $r'$, but for notational simplicity we will assume
a single cross section. We will also assume that the cross section is small
enough that the flux is essentially constant over the whole earth so that the
rate may be written as:
\begin{equation}
R(r') = \Phi \eta(r') \sigma(r').
\end{equation}
The rate may then be written as:
\begin{equation}
R = \Phi \int_\mathrm{SNO} \mathrm{d}^3r \, \int_{r'} \mathrm{d}^3r' \, \eta(r') \sigma(r') \mathrm{P}(\text{detect at r} | \text{DM scatters at r'})
\end{equation}
If we assume that the probability of detecting the dark matter is uniform
across the SNO detector we may write it as:
\begin{equation}
R = \Phi \int_\mathrm{SNO} \mathrm{d}^3r \, \int_{r'} \mathrm{d}^3r' \, \eta(r') \sigma(r') \mathrm{P}(\text{detect at SNO} | \text{DM scatters at r'})
\end{equation}
This assumption is pretty well motivated since for most values of the mediator
decay length the probability will be uniform across the detector. The only
value for which it might not be a good approximation is if the mediator decay
length is on the order of the detector radius in which case DM scattering in
the rock of the cavity walls might have a higher event rate at the edge of the
detector. Since the integrand no longer depends on $r$, we may write it as:
\begin{equation}
R = \Phi V_\text{fiducial} \int_{r'} \mathrm{d}^3r' \, \eta(r') \sigma(r') \mathrm{P}(\text{detect at SNO} | \text{DM scatters at r'})
\end{equation}
where $V_\text{fiducial}$ is the fiducial volume of the detector. The
probability that the mediator $V$ is emitted in a direction $\theta$ and
travels a distance $r$ in a spherical coordinate system centered on $r'$ may be
written as:
\begin{equation}
f(r,\theta) = \frac{\sin\theta}{4\pi}\frac{e^{-r/L_V}}{L_V}
\end{equation}
To transform this probability distribution to the coordinate system centered on the SNO detector we first transform it to a cartesian coordinate system:
\begin{align}
f(r,\theta) &= \frac{\sin\theta}{4\pi}\frac{e^{-r/L_V}}{L_V}\frac{1}{r^2sin\theta} \\
&= \frac{1}{4\pi r^2}\frac{e^{-r/L_V}}{L_V}
\end{align}
Then, the distribution is translated to the center of the detector, which
doesn't change the form since the radial coordinate $r$ is the same in both
coordinate systems. Finally, we switch back into spherical coordinates:
\begin{align}
f(r,\theta') &= \frac{1}{4\pi r^2}\frac{e^{-r/L_V}}{L_V}r^2\sin\theta' \\
&= \frac{1}{4\pi}\frac{e^{-r/L_V}}{L_V}\sin\theta'
\end{align}
where $\theta'$ is the polar angle in the SNO coordinate system. We can now write the rate as:
\begin{align}
R &= \Phi V_\text{fiducial} \int_r \mathrm{d}r \, \int_\theta \mathrm{d}\theta \, \int_\phi \mathrm{d}\phi \, \eta(r,\theta,\phi) \sigma(r,\theta,\phi) \frac{1}{4\pi}\frac{e^{-r/L_V}}{L_V}\sin\theta
\end{align}
We now assume that the number density of scatterers $\eta$ and the cross
section $\sigma$ are independent of the position in the earth. This is a good
approximation for certain values of $L_V$ since the integral will be dominated
by a single material. For example, if the mediator decay length $L_V$ is
approximately 1 meter, then the vast majority of the events in the detector
will be caused by dark matter scattering off of water. Similarly if the
mediator decay length is approximately 1 km then the majority of the events in
the detector will be caused by the dark matter scattering off of the norite
rock surrounding the detector. With this approximation, the rate may be
written:
\begin{align}
R &= \Phi V_\text{fiducial} \eta \sigma \int_r \mathrm{d}r \, \int_\theta \mathrm{d}\theta \, \int_\phi \mathrm{d}\phi \, \frac{1}{4\pi}\frac{e^{-r/L_V}}{L_V}\sin\theta \\
&= \Phi V_\text{fiducial} \eta \sigma \int_r \mathrm{d}r \, \frac{1}{2}\frac{e^{-r/L_V}}{L_V} \int_\theta \mathrm{d}\theta \, \sin\theta
\end{align}
The $\theta$ integral goes from $\theta_\text{min}$ to $\pi$:
\begin{align}
R &= \Phi V_\text{fiducial} \eta \sigma \int_r \mathrm{d}r \, \frac{1}{2}\frac{e^{-r/L_V}}{L_V} \int_{\theta_\text{min}}^\pi \mathrm{d}\theta \, \sin\theta
\end{align}
where $\theta_\text{min}$ is equal to:
\begin{equation}
\theta_\text{min} =%
\begin{cases}
0 & \text{if } r < \text{depth} \\
\pi - \arccos\left(\frac{\text{depth}^2 + r^2 - 2R\text{depth}}{2r(R-\text{depth})}\right) & \text{if } \text{depth} < r < 2R-\text{depth} \\
\end{cases}
\end{equation}
where $R$ is the radius of the earth and $\text{depth}$ is the distance from
the surface of the earth to the SNO detector.

\section{Cross Section}
In \cite{grossman2017} the differential scattering cross section for dark
matter off a nucleus is calculated as
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} = \frac{g_V^2 \epsilon^2 e^2}{4\pi v^2 (q^2 + m_V^2)^2} |F_D(q^2)|^2 Z^2 F^2(q),
\end{equation}
where $q$ is the momentum transferred, $g_V$ and $\epsilon$ are coupling
constants (FIXME: is this true?), $v$ is the velocity of the dark matter
particle, $m_V$ is the mass of the mediator, $F_D(q^2)$ is a form factor for
the dark matter to transition from a high angular momentum state to a lower
angular momentum state, $Z$ is the atomic number of the nucleus, and $F^2(q)$
is a nuclear form factor.

In the limit of low momentum transfer, the cross section is approximately
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} \simeq \frac{g_V^2 \epsilon^2 e^2}{4\pi v^2 m_V^4} |F_D(q^2)|^2 Z^2 F^2(q).
\end{equation}

For existing direct detection dark matter experiments, the relevant cross
section is (FIXME: is this true?)
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} \simeq \frac{g_V^2 \epsilon^2 e^2}{4\pi v^2 m_V^4} Z^2 F^2(q).
\end{equation}

A standard cross section can be defined as the total cross section in the zero
momentum limit\cite{pepin2016}
\begin{align}
\sigma_0 &= \int_0^{4\mu_T^2 v^2} \frac{\diff \sigma_\text{scatter}}{\diff q^2}\bigg\rvert_{q \rightarrow 0} \diff q^2 \\
&= \frac{\mu_T^2 g_V^2 \epsilon^2 e^2}{\pi m_V^4} Z^2,
\end{align}
where $\mu_T$ is the reduced mass of the WIMP and target nucleus.

Since different experiments use different detector targets, it is also useful
to define a standard cross section, $\sigma_p$ which is independent of the
nuclear target:
\begin{equation}
\sigma_p = \left(\frac{\mu_p}{\mu_T}\frac{1}{Z}\right)^2 \sigma_0.
\end{equation}

The direct detection cross section is then:
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} \simeq \frac{1}{4 \mu_p^2 v^2} \sigma_p Z^2 F^2(q).
\end{equation}
and the cross section for the dark matter to annihilate is:
\begin{equation}
\frac{\diff \sigma_\text{scatter}}{\diff q^2} \simeq \frac{1}{4 \mu_p^2 v^2} \sigma_p |F_D(q)|^2 Z^2 F^2(q).
\end{equation}

\subsection{Nuclear Form Factor}
The nuclear form factor, $F(q)$, characterizes the loss of coherence as the de
Broglie wavelength of the WIMP approaches the radius of the
nucleus\cite{caldwell2015}. The most commonly used form factor calculation used
in the direct detection community is that of Helm which is given by:
\begin{equation}
F(q) = 3\frac{j_1(q r_1)}{q r_1} e^{-\frac{(q s)^2}{2}},
\end{equation}
where $j_1$ is the spherical bessel function of the first order, $s$ is a
measure of the nuclear skin thickness, and $r_1$ is a measure of the nuclear
radius\cite{pepin2016}. The values used for these constants were
\begin{align}
s &= 0.9 \text{ fm} \\
a &= 0.52 \text{ fm} \\
c &= 1.23 A^\frac{1}{3} - 0.60 \text{ fm} \\
r_1 &= \sqrt{c^2 + \frac{7}{3}\pi^2 a^2 - 5 s^2}
\end{align}

\begin{figure}
\centering
\begin{tikzpicture}[scale=0.1]
% earth
\draw [thick,domain=120:150] plot[smooth] ({200*cos(\x)},{200*sin(\x)});
\begin{scope}[shift={(-100,100)},rotate=45]
% interaction
\node[star,star points=9,draw] at (20,20){};
% acrylic vessel
\draw [thick,domain={90+asin(75/600)}:{360+90-asin(75/600)}] plot[smooth] ({6*cos(\x)},{6*sin(\x)});
\draw [thick] ({6*cos(90+asin(75/600))},{6*sin(90+asin(75/600))}) -- ({6*cos(90+asin(75/600))},{6*sin(90+asin(75/600))+7.5}) --
      ({6*cos(90+asin(75/600))+2*0.75},{6*sin(90+asin(75/600))+7.5}) --
      ({6*cos(90+asin(75/600))+2*0.75},{6*sin(90+asin(75/600))});
% PSUP
\draw [domain=0:360] plot ({8.89*cos(\x)},{8.89*sin(\x)});
% cavity
\draw (-9.5,-10.5) --
      (-10.6,-10.5+5.6) --
      (-10.6,-10.5+14.93) --
      (-9.5,-10.5+30) --
      (9.5,-10.5+30) --
      (10.6,-10.5+14.93) --
      (10.6,-10.5+5.6) --
      (9.5,-10.5) --
      (-9.5,-10.5);
\draw[->,ultra thick] (-25,0) -- (25,0) node[right]{$x$};
\draw[->,ultra thick] (0,-25) -- (0,25) node[right]{$y$};
\end{scope}
\end{tikzpicture}
\end{figure}

\section{Event Reconstruction}
In order to reconstruct the physical parameters associated with an event we
compute a likelihood for that event given a proposed energy, position,
direction, and initial time. The likelihood may be written as:
\begin{equation}
\label{likelihood}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = P(\vec{q}, \vec{t} | E, \vec{x}, \vec{v}, t_0)
\end{equation}
where $E$, $\vec{x}$, $\vec{v}$ represent the initial particle's kinetic
energy, position, and direction respectively, $t_0$ represents the initial time
of the event, $\vec{q}$ is the charge seen by each PMT, and $\vec{t}$ is the
time recorded by each PMT.

In general the right hand side of Equation~\ref{likelihood} is not factorable
since for particle tracks which scatter there will be correlations between the
PMT hits. However, to make the problem analytically tractable, we assume that
the probability of each PMT being hit is approximately independent of the
others. With this assumption we can factor the right hand side of the
likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i P(\text{not hit} | E, \vec{x}, \vec{v}, t_0) \prod_j P(\text{hit}, q_j, t_j | E, \vec{x}, \vec{v}, t_0)
\end{equation}
where the first product is over all PMTs which weren't hit and the second
product is over all of the hit PMTs.

If we introduce the variable $n$ which represents the number of photoelectrons detected we can write the likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i P(n = 0 | E, \vec{x}, \vec{v}, t_0) \prod_j \sum_{n = 1}^{\infty} P(n, q_j, t_j | E, \vec{x}, \vec{v}, t_0)
\end{equation}

We can factor the right hand side of the likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i P(n = 0 | E, \vec{x}, \vec{v}, t_0) \prod_j \sum_{n = 1}^{\infty} P(q_j, t_j | n, E, \vec{x}, \vec{v}, t_0) P(n | E, \vec{x}, \vec{v}, t_0)
\end{equation}

If we now assume that the charge and time observed at a given PMT are
independent we can write the likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i P(n = 0 | E, \vec{x}, \vec{v}, t_0) \prod_j \sum_{n = 1}^{\infty} P(q_j | n, E, \vec{x}, \vec{v}, t_0) P(t_j | n, E, \vec{x}, \vec{v}, t_0) P(n | E, \vec{x}, \vec{v}, t_0)
\end{equation}

Since there are many photons produced in each event each of which has a small
probability to hit a given PMT, we will assume that the probability of
detecting $n$ photons at a given PMT is poisson distributed, i.e.
\begin{equation}
P(n | E, \vec{x}, \vec{v}, t_0) = e^{-\mu} \frac{\mu^n}{n!}
\end{equation}

We can therefore write the likelihood as:
\begin{equation}
\mathcal{L}(E, \vec{x}, \vec{v}, t_0) = \prod_i e^{-\mu_i} \prod_j \sum_{n = 1}^{\infty} P(q_j | n, E, \vec{x}, \vec{v}, t_0) P(t_j | n, E, \vec{x}, \vec{v}, t_0) e^{-\mu_j} \frac{\mu_j^n}{n!}
\end{equation}
where $\mu_i$ is the expected number of photoelectrons detected at the ith PMT
(given an initial particle's energy, position, and direction).

First, we'll calculate the expected number of photoelectrons for a single non-showering track which undergoes multiple scattering through small angles. In this case, we can calculate the expected number of photoelectrons as:
\begin{equation}
\mu_i = \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} P(\text{detected} | E, x, v)
\end{equation}
where $x$ is the position along the track and $\lambda$ is the wavelength of
the light.

If the particle undergoes many small angle Coulomb scatters, the net
angular displacement of the particle after a distance $x$ will be a Gaussian
distribution by the central limit theorem\cite{pdg2017}. The distribution of
the net angular displacement at a distance $x$ along the track is then given
by\footnote{This distribution will be correlated between different points along the track.}:
\begin{equation}
f(\theta,\phi) = \frac{\theta}{2\pi\theta_0^2}e^{-\frac{\theta^2}{2\theta_0^2}}
\end{equation}
where
\begin{equation}
\theta_0 = \frac{13.6 \text{ MeV}}{\beta c p}z\sqrt{\frac{x}{X_0}}\left[1 + 0.038\ln\left(\frac{x z^2}{X_0 \beta^2}\right)\right]
\end{equation}
where $p$, $\beta c$, and $z$ are the momentum, velocity, and charge of the
particle, and $X_0$ is the radiation length of the particle\cite{pdg2017}.

Now, we integrate over the angular displacement of the track around the original velocity:
\begin{align}
\mu_i &= \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \int_\theta \diff \theta \int_\phi \diff \phi P(\text{detected} | \theta, \phi, E, x, v) P(\theta, \phi | E, x, v) \\
\mu_i &= \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \int_\theta \diff \theta \int_\phi \diff \phi P(\text{detected} | \theta, \phi, E, x, v) f(\theta,\phi)
\end{align}
The probability of being detected can be factored into several different compontents:
\begin{align}
\mu_i &= \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \int_\theta \diff \theta \int_\phi \diff \phi P(\text{emitted towards PMT i} | \theta, \phi, E, x, v) f(\theta,\phi) P(\text{not scattered or absorbed} | \lambda, E, x, v) \epsilon(\eta) \mathrm{QE}(\lambda) \\
\label{eq:mui}
\mu_i &= \int_x \diff x \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} P(\text{not scattered or absorbed} | \lambda, E, x, v) \epsilon(\eta) \mathrm{QE}(\lambda) \int_\theta \diff \theta \int_\phi \diff \phi P(\text{emitted towards PMT i} | \theta, \phi, E, x, v) f(\theta,\phi)
\end{align}
where $\eta$ is the angle between the vector connecting the track position $x$
to the PMT position and the normal vector to the PMT, $\epsilon(\eta)$ is the
collection efficiency, and $\mathrm{QE}(\lambda)$ is the quantum efficiency of
the PMT.

The probability that a photon is emitted directly towards a PMT is given by a
delta function (we make the assumption here that the probability is uniform
across the face of the PMT):
\begin{equation}
P(\text{emitted towards PMT i} | \theta, \phi, E, x, v) = \delta\left(\frac{1}{n(\lambda)\beta} - \cos\theta'(\theta,\phi,x)\right) \frac{\Omega(x)}{4\pi}
\end{equation}
where $\theta'$ is the angle between the track and the PMT and $\Omega(x)$ is the solid angle subtended by the PMT.

In a coordinate system with the z axis aligned along the original particle velocity and with the PMT in the x-z plane, the angle $\theta'$ is defined by:
\begin{equation}
\cos\theta' = \sin\theta\cos\phi\sin\theta_1 + \cos\theta\cos\theta_1
\end{equation}
where $\theta_1$ is the angle between the PMT and the original particle velocity.

We can now solve the integral on the right hand side of Equation~\ref{eq:mui} as:
\begin{align}
P(\text{emitted towards PMT i}) &= \int_\theta \diff \theta \int_\phi \diff \phi \delta\left(\frac{1}{n(\lambda)\beta} - \cos\theta'(\theta,\phi,x)\right) \frac{\Omega(x)}{4\pi} \frac{\theta}{2\pi\theta_0^2}e^{-\frac{\theta^2}{2\theta_0^2}} \\
P(\text{emitted towards PMT i}) &= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2}\int_\theta \diff \theta \int_\phi \diff \phi \delta\left(\frac{1}{n(\lambda)\beta} - \cos\theta'(\theta,\phi,x)\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}} \\
P(\text{emitted towards PMT i}) &= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2}\int_\theta \diff \theta \int_\phi \diff \phi \delta\left(\frac{1}{n(\lambda)\beta} - \sin\theta\cos\phi\sin\theta_1 - \cos\theta\cos\theta_1\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}}
\end{align}

We now assume $\theta$ is small (which should be valid for small angle scatters), so that we can rewrite the delta function as:
\begin{align}
P(\text{emitted towards PMT i}) &= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2}\int_\theta \diff \theta \int_\phi \diff \phi \delta\left(\frac{1}{n(\lambda)\beta} - \theta\cos\phi\sin\theta_1 - \cos\theta_1\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}}
\end{align}

We can rewrite the delta function and solve the integral as:
\begin{align}
P(\text{emitted towards PMT i}) &= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2}\int_\theta \diff \theta \int_\phi \diff \phi \frac{1}{\left|\cos\phi\sin\theta_1\right|}\delta\left(\theta - \frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}} \\
&= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2} \frac{1}{\left|\sin\theta_1\right|} \int_\phi \diff \phi \frac{1}{\left|\cos\phi\right|} \int_\theta \diff \theta \delta\left(\theta - \frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}\right) \theta e^{-\frac{\theta^2}{2\theta_0^2}} \\
&= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2} \frac{1}{\left|\sin\theta_1\right|} \int_\phi \diff \phi \frac{1}{\left|\cos\phi\right|}\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}H\left(\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}\right)e^{-\frac{1}{2\theta_0^2}\left(\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\cos\phi\sin\theta_1}\right)^2} \\
&= \frac{\Omega(x)}{4\pi} \frac{1}{2\pi\theta_0^2} \frac{1}{\left|\sin\theta_1\right|}\sqrt{2\pi}\theta_0 e^{-\frac{1}{2\theta_0^2}\left(\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\sin\theta_1}\right)^2} \\
&= \frac{\Omega(x)}{4\pi} \frac{1}{\sqrt{2\pi}\theta_0} \frac{1}{\left|\sin\theta_1\right|} e^{-\frac{1}{2\theta_0^2}\left(\frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\sin\theta_1}\right)^2}.
\end{align}

To simplify this expression we can write
\begin{equation}
P(\text{emitted towards PMT i}) = \frac{\Omega(x)}{4\pi} \frac{1}{\sqrt{2\pi}\theta_0} \frac{1}{\left|\sin\theta_1\right|} e^{-\frac{\Delta^2(\lambda)}{2\theta_0^2}}
\end{equation}
where
\begin{equation}
\Delta(\lambda) = \frac{\frac{1}{n(\lambda)\beta}-\cos\theta_1}{\sin\theta_1}
\end{equation}

Plugging this back into Equation~\ref{eq:mui}
\begin{align}
\label{eq:mui-exact}
\mu_i &= \frac{1}{\sqrt{2\pi}\theta_0} \int_x \diff x \frac{\Omega(x)}{4\pi} \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} P(\text{not scattered or absorbed} | \lambda, E, x, v) \mathrm{QE}(\lambda) e^{-\frac{\Delta^2(\lambda)}{2\theta_0^2}}
\end{align}

Ideally we would just evaluate this double integral for each likelihood call,
however the double integral is too computationally expensive to perform for
every likelihood call (FIXME: is this true?). We therefore assume that the
second integral will be dominated by the Bessel function which has a
singularity when it's argument is zero, and rewrite Equation~\ref{eq:mui-exact}
as:
\begin{align}
\mu_i &= 2 \frac{1}{\sqrt{2\pi}\theta_0} \int_x \diff x \Omega(x) \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) P(\text{not scattered or absorbed} | \lambda_0, E, x, v) \mathrm{QE}(\lambda_0) e^{-\frac{\Delta^2(\lambda_0)}{4\theta_0^2}} \int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} K_0\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right)
\end{align}
where $\lambda_0$ is the wavelength at which $\Delta(\lambda) = 0$.

For small values of $\Delta$, the Bessel function may be approximated as:
\begin{equation}
K_0(x) \simeq -\log(x) + \log(2) - \gamma
\end{equation}

We may therefore approximate the expected charge as
\begin{multline}
\label{eq:mui-approx}
\mu_i = 2 \frac{1}{\sqrt{2\pi}\theta_0} \int_x \diff x \Omega(x) \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) P(\text{not scattered or absorbed} | \lambda_0, E, x, v) \mathrm{QE}(\lambda_0) e^{-\frac{\Delta^2(\lambda_0)}{4\theta_0^2}} \\
\int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right)
\end{multline}

The number of Cerenkov photons produced per unit length and per unit wavelength
is given by\cite{pdg2017}
\begin{equation}
\frac{\diff^2 N}{\diff x \diff \lambda} = \frac{2\pi\alpha z^2}{\lambda^2}\left(1 - \frac{1}{\beta^2 n^2(\lambda)}\right)
\end{equation}
where $\alpha$ is the fine-structure constant and $z$ is the charge of the
particle in units of the electron charge.

We can therefore write the second integral in Equation~\ref{eq:mui-approx} as
\begin{align}
\int_\lambda \diff \lambda \frac{\diff^2 N}{\diff x \diff \lambda} \left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right) &=
2\pi\alpha z^2 \int_\lambda \diff \lambda \frac{1}{\lambda^2}\left(1 - \frac{1}{\beta^2 n^2(\lambda)}\right) \left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right) \\
\label{eq:lambda}
&\simeq 2\pi\alpha z^2 \left(1 - \frac{1}{\beta^2 n^2(\lambda_0)}\right) \int_\lambda \diff \lambda \frac{1}{\lambda^2}\left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right)
\end{align}

Since the $\Delta$ function only depends on the wavelength through the index
which depends weakly on the wavelength, we can approximate the index of
refraction as:
\begin{equation}
n(\lambda) \simeq a + \frac{b}{\lambda^2}.
\end{equation}

The integral in Equation~\ref{eq:lambda} may then be solved analytically
\begin{multline}
\int_{\lambda_1}^{\lambda_2} \diff \lambda \frac{1}{\lambda^2}\left(-\log\left(\frac{\Delta^2(\lambda)}{4\theta_0^2}\right) + \log(2) - \gamma\right) =
\left[\log(4\theta_0^2) + \log(\sin^2\theta_1) + \log(2) - \gamma\right]\left(\frac{1}{\lambda_1}-\frac{1}{\lambda_2}\right) + \\
\left.\left(-4\sqrt{\frac{a}{b}}\arctan\left(\sqrt{\frac{a}{b}}\lambda\right) +
4\sqrt{\frac{1+a\Delta_0}{b\Delta_0}}\arctan\left(\sqrt{\frac{1+a\Delta_0}{b\Delta_0}}\lambda\right) -
\frac{1}{\lambda}\log\left[\left(\Delta_0+\frac{\lambda^2}{b+a\lambda^2}\right)^2\right]\right)\right|_{\lambda_1}^{\lambda_2}
\end{multline}
where $\lambda_1$ and $\lambda_2$ are chosen to cover the range where the
quantum efficiency is non-zero, typically between 300 nm and 600 nm.

For simplicity we will write this previous expression as $f(x)$
\begin{multline}
f(x) = \left[\log(4\theta_0^2) + \log(\sin^2\theta_1) + \log(2) - \gamma\right]\left(\frac{1}{\lambda_1}-\frac{1}{\lambda_2}\right) + \\
\left.\left(-4\sqrt{\frac{a}{b}}\arctan\left(\sqrt{\frac{a}{b}}\lambda\right) +
4\sqrt{\frac{1+a\Delta_0}{b\Delta_0}}\arctan\left(\sqrt{\frac{1+a\Delta_0}{b\Delta_0}}\lambda\right) -
\frac{1}{\lambda}\log\left[\left(\Delta_0+\frac{\lambda^2}{b+a\lambda^2}\right)^2\right]\right)\right|_{\lambda_1}^{\lambda_2}
\end{multline}

We can now write Equation~\ref{eq:mui-approx} as
\begin{equation}
\mu_i = 2 \frac{1}{\sqrt{2\pi}\theta_0} 2\pi\alpha z^2 \int_x \diff x \Omega(x) \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) P(\text{not scattered or absorbed} | \lambda_0, E, x, v) \mathrm{QE}(\lambda_0) e^{-\frac{\Delta^2(\lambda_0)}{4\theta_0^2}} \left(1 - \frac{1}{\beta^2 n^2(\lambda_0)}\right) f(x)
\end{equation}

The probability that a photon travels to the PMT without being scattered or absorbed can be calculated as follows
\begin{align}
P(\text{not scattered or absorbed} | \lambda, x) &=
P(\text{not scattered} | \lambda, x) P(\text{not absorbed} | \lambda, x) \\
&= \int_l^\infty\frac{1}{s(\lambda)}e^{-\frac{x}{s(\lambda)}}\int_l^\infty\frac{1}{a(\lambda)}e^{-\frac{x}{a(\lambda)}} \\
&= e^{-\frac{l}{s(\lambda) + a(\lambda)}}
\end{align}
where $l$ is the distance to the PMT from the position $x$, $s(\lambda)$ is the
scattering length, and $a(\lambda)$ is the absorption length.

We can therefore write the expected charge as
\begin{equation}
\mu_i = 2 \frac{1}{\sqrt{2\pi}\theta_0} 2\pi\alpha z^2 \int_x \diff x \Omega(x) \frac{1}{\left|\sin\theta_1\right|} \epsilon(\eta) e^{-\frac{l(x)}{s(\lambda) + a(\lambda)}} \mathrm{QE}(\lambda_0) e^{-\frac{\Delta^2(\lambda_0)}{4\theta_0^2}} \left(1 - \frac{1}{\beta^2 n^2(\lambda_0)}\right) f(x)
\end{equation}

The last integral is calculated numerically when the likelihood is evaluated.

We now return to the likelihood and calculate the probability of observing a
given time. In principle, this depends on the number of photons hitting a PMT
since the PMT hit will only register the \emph{first} photoelectron which
crosses threshold. However, since this is expected to be a small effect, we
assume that the probability of observing a given time is independent of the
number of photons which hit the PMT, i.e.
\begin{equation}
P(t_j | n, E, \vec{x}, \vec{v}, t_0) \simeq P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0)
\end{equation}

We first condition on the \emph{true} time at which the photon hits the PMT
\begin{equation}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0) = \int_{t_j'} \diff t P(t_j | t_j') P(t_j' | n \geq 1, E, \vec{x}, \vec{v}, t_0)
\end{equation}
where we used the fact that the probability of a measured time only depends on the true PMT hit time.

Now, we integrate over the track
\begin{equation}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0) = \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x P(t_j', x | n \geq 1, E, \vec{x}, \vec{v}, t_0)
\end{equation}
where $x$ here stands for the event that a photon emitted at a distance $x$ along the track makes it to the PMT.

We now use Bayes theorem to rewrite the last probability
\begin{align}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0) &= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x P(t_j', x | n \geq 1, E, \vec{x}, \vec{v}, t_0) \\
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x P(t_j' | x, n \geq 1, E, \vec{x}, \vec{v}, t_0) P(x | n \geq 1, E, \vec{x}, \vec{v}, t_0) \\
\end{align}
The first term in the integral is just a delta function (up to slight differences due to dispersion) since we are assuming direct light
\begin{align}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0)
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) P(x | n \geq 1, E, \vec{x}, \vec{v}, t_0) \\
\end{align}

We now use Bayes theorem to rewrite the last term
\begin{align}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0)
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \frac{P(n \geq 1 | x, E, \vec{x}, \vec{v}, t_0) P(x | E, \vec{x}, \vec{v}, t_0)}{P(n \geq 1 | E, \vec{x}, \vec{v}, t_0)} \\
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \frac{P(x | E, \vec{x}, \vec{v}, t_0)}{P(n \geq 1 | E, \vec{x}, \vec{v}, t_0)} \\
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \frac{P(x | E, \vec{x}, \vec{v}, t_0)}{1 - e^{-\mu_j}} \\
&= \int_{t_j'} \diff t_j' P(t_j | t_j') \int_x \diff x \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \frac{\mu_j(x)}{1 - e^{-\mu_j}} \\
&= \frac{1}{1 - e^{-\mu_j}} \int_x \diff x \mu_j(x) \int_{t_j'} \diff t_j' P(t_j | t_j') \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \\
\end{align}

We assume the transit time spread is equal to a gaussian (we ignore the pre and late pulsing)
\begin{align}
P(t_j | n \geq 1, E, \vec{x}, \vec{v}, t_0)
&= \frac{1}{1 - e^{-\mu_j}} \int_x \diff x \mu_j(x) \int_{t_j'} \diff t_j' \frac{1}{\sqrt{2\pi}\sigma_t} e^{-\frac{(t_j-t_j')^2}{2\sigma_t^2}} \delta\left(\frac{l(x)n(\lambda_0)}{c}-t_j'\right) \\
&= \frac{1}{1 - e^{-\mu_j}} \frac{1}{\sqrt{2\pi}\sigma_t} \int_x \diff x \mu_j(x) e^{-\frac{(t_j-t_0(x))^2}{2\sigma_t^2}}
\end{align}
where in the last expression we define
\begin{equation}
t_0(x) \equiv \frac{l(x)n(\lambda_0)}{c}
\end{equation}

\begin{thebibliography}{9}
\bibitem{grossman2017}
    Grossman, et al. \textit{Self-Destructing Dark Matter}. \href{https://arxiv.org/abs/1712.00455}{{\tt arXiv:1712.00455}}. Dec 2017.
\bibitem{pepin2016}
    M. Pepin. \textit{Low-Mass Dark Matter Search Results and Radiogenic Backgrounds for the Cryogenic Dark Matter Search}. \url{http://hdl.handle.net/11299/185144}. Dec 2016.
\bibitem{caldwell2015}
    T. Caldwell. \textit{Searching for Dark Matter with Single Phase Liquid Argon}. \url{https://repository.upenn.edu/dissertations/1632}. 2015.
\bibitem{pdg2017}
    C. Patrignani et al. (Particle Data Group), Chin. Phys. C, 40, 100001 (2016) and 2017 update.
\end{thebibliography}
\end{document}