diff options
| -rw-r--r-- | .gitignore | 5 | ||||
| -rw-r--r-- | Makefile | 11 | ||||
| -rwxr-xr-x | calculate_limits.py | 419 | ||||
| -rw-r--r-- | sddm.tex | 500 |
4 files changed, 935 insertions, 0 deletions
@@ -17,3 +17,8 @@ test-find-peaks *.xzdab *.tar.gz *.hdf5 +calculate_limits +*.aux +*.log +*.pdf +*.out @@ -7,3 +7,14 @@ default: all cd src && "$(MAKE)" $@ cd utils && "$(MAKE)" $@ +all: sddm.pdf + +calculate_limits: calculate_limits.c + +sddm.pdf: sddm.tex + pdflatex $^ + +clean: + rm calculate_limits + +.PHONY: clean sddm.pdf diff --git a/calculate_limits.py b/calculate_limits.py new file mode 100755 index 0000000..f297a69 --- /dev/null +++ b/calculate_limits.py @@ -0,0 +1,419 @@ +#!/usr/bin/env python3 +import numpy as np +from scipy.integrate import quad +from scipy.stats import expon +from scipy.special import spherical_jn, erf +from numpy import pi +from functools import lru_cache + +# on retina screens, the default plots are way too small +# by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1 +# Qt5 will scale everything using the dpi in ~/.Xresources +import matplotlib +matplotlib.use("Qt5Agg") + +# speed of light (exact) +SPEED_OF_LIGHT = 299792458 # m/s + +# depth of the SNO detector (m) +# currently just converted 6800 feet -> meters +h = 2072.0 + +# radius of earth (m) +# from the "Earth Fact Sheet" at https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html +# we use the volumetric mean radius here +R = 6.371e6 + +# Approximate dark matter velocity in m/s. The true distribution is expected to +# be a Maxwell Boltzmann distribution which is modulated annually by the +# earth's rotation around the sun, but we just assume a single constant +# velocity here. From Lewin and Smith Appendix B +DM_VELOCITY = 230e3 + +# Approximate earth velocity in the galactic rest frame (?) +# from Lewin and Smith Equation 3.6 +EARTH_VELOCITY = 244e3 + +# Approximate dark matter density in GeV/m^3. From Tom Caldwell's thesis. +# from Lewin and Smith page 91 +DM_DENSITY = 0.4e6 + +# Number density of scatterers in the Earth. +# +# FIXME: Currently just set to the number density of atoms in water. Need to +# update this for rock, and in fact this will change near the detector since +# there is water outside the AV. +DENSITY_WATER = 1e3 # In kg/m^3 + +# mean density of norite rock which is the rock surrounding SNO +# probably conservative +NORITE_DENSITY = 3e3 # In kg/m^3 + +# atomic masses for various elements +# from https://www.angelo.edu/faculty/kboudrea/periodic/structure_mass.htm +element_mass = { + 'H':1.0, + 'C':12.011, + 'O':15.9994, + 'Na':22.98977, + 'Mg':24.305, + 'Al':26.98154, + 'Si':28.0855, + 'K':39.0983, + 'Ca':40.08, + 'Mn':54.9380, + 'Fe':55.847, + 'Ti':47.90 +} + +# composition and mean density of norite rock which is the rock surrounding SNO +# the composition is from Table 3.2 in the SNOLAB User's handbook +water = {'composition': + {'H':20, + 'O':80}, + 'density':1e3} + +# composition and mean density of the mantle +# from www.knowledgedoor.com/2/elements_handbook/element_abundances_in_the_earth_s_mantle.html +# density from hyperphysics.phy-astr.gsu.edu/hbase/Geophys/earthstruct.html +mantle = {'composition': + {'O':44.33, + 'Mg':22.17, + 'Si':21.22, + 'Fe':6.3}, + 'density':4.4e3} + +# composition and mean density of norite rock which is the rock surrounding SNO +# the composition is from Table 3.2 in the SNOLAB User's handbook +norite = {'composition': + {'H':0.15, + 'C':0.04, + 'O':46.0, + 'Na':2.2, + 'Mg':3.3, + 'Al':9.0, + 'Si':26.2, + 'K': 1.2, + 'Ca':5.2, + 'Mn':0.1, + 'Fe':6.2, + 'Ti':0.5}, + 'density':3e3} + +# Fiducial volume (m) +FIDUCIAL_RADIUS = 5 + +# Fiducial volume (m^3) +FIDUCIAL_VOLUME = 4*pi*FIDUCIAL_RADIUS**3/3 + +# proton mass from the PDG (2018) +PROTON_MASS = 0.938 # GeV + +# proton mass from the PDG (2018) +ATOMIC_MASS_UNIT = 0.931 # GeV + +# mass of Xenon in atomic units +XENON_MASS = 131.293 + +# mass of Neon in atomic units +NEON_MASS = 20.18 + +# mass of argon in atomic units +ARGON_MASS = 39.948 + +# mass of germanium in atomic units +GERMANIUM_MASS = 72.64 + +# mass of tungsten in atomic units +TUNGSTEN_MASS = 183.84 + +# mass of oxygen in atomic units +OXYGEN_MASS = 15.999 + +# mass of silicon in atomic units +SILICON_MASS = 28.0855 + +# mass of iron in atomic units +IRON_MASS = 55.845 + +# mass of magnesium in atomic units +MAGNESIUM_MASS = 24.305 + +# galactic escape velocity (m/s) +# from Tom Caldwell's thesis page 25 +ESCAPE_VELOCITY = 244e3 + +# conversion constant from PDG +HBARC = 197.326978812e-3 # GeV fm + +# Avogadros number (kg^-1) +N0 = 6.02214085774e26 + +def get_probability(r, l): + """ + Returns the probability of a dark photon decaying in the SNO detector from + a dark matter decay distributed uniformly in the Earth. Assumes that the + depth of SNO is much larger than the dimensions of the SNO detector. + """ + if r <= h: + theta_min = 0 + elif r <= 2*R - h: + theta_min = pi - np.arccos((h**2 + r**2 - 2*R*h)/(2*r*(R-h))) + else: + return 0 + return (1 + np.cos(theta_min))*np.exp(-r/l)/(2*l) + +@lru_cache() +def get_probability2(l): + p, err = quad(get_probability, 0, min(2*R-h,10*l), args=(l,), epsabs=0, epsrel=1e-7, limit=1000) + + return p + +def get_event_rate(m, cs0, l, A): + """ + Returns the event rate of leptons produced from self-destructing dark + matter in the SNO detector for a given dark matter mass m, a cross section + cs, and a mediator decay length l. + """ + # For now we assume the event rate is constant throughout the earth, so we + # are implicitly assuming that the cross section is pretty small. + flux = DM_VELOCITY*DM_DENSITY/m + + #p, err = quad(get_probability, 0, min(2*R-h,10*l), args=(l,), epsabs=0, epsrel=1e-7, limit=1000) + p = get_probability2(l) + + cs = get_cross_section(cs0, m, A) + + # FIXME: factor of 2 because the DM particle decays into two mediators? + return p*cs*(N0/A)*flux*FIDUCIAL_VOLUME + +def get_event_rate_sno(m, cs0, l, composition): + """ + Returns the event rate of leptons produced from self-destructing dark + matter in the SNO detector for a given dark matter mass m, a cross section + cs, and a mediator decay length l. + """ + rate = 0.0 + + for element, mass_fraction in composition['composition'].items(): + rate += mass_fraction/100*get_event_rate(m,cs0,l,element_mass[element]) + + return rate*composition['density'] + +def get_nuclear_form_factor(A, e): + """ + Returns the nuclear form factor for a WIMP-nucleus interaction. + + From Tom Caldwell's thesis page 24. + + Also used Mark Pepin's thesis page 50 + """ + # mass of xenon nucleus + mn = A*ATOMIC_MASS_UNIT + + # calculate approximate size of radius + s = 0.9 # fm + a = 0.52 # fm + c = 1.23*A**(1/3) - 0.60 # fm + r1 = np.sqrt(c**2 + (7/3)*pi**2*a**2 - 5*s**2) + + q = np.sqrt(2*mn*e) + + if q*r1/HBARC < 1e-10: + return 1.0 + + # Helm form factor + # from Mark Pepin's thesis page 50 + f = 3*spherical_jn(1,q*r1/HBARC)*np.exp(-(q*s/HBARC)**2/2)/(q*r1/HBARC) + + return f + +def get_cross_section(cs0, m, A): + """ + Returns the WIMP cross section from the target-independent WIMP-nucleon + cross section cs0 at zero momentum transfer. + + From Tom Caldwell's thesis page 21. + """ + # mass of xenon nucleus + mn = A*ATOMIC_MASS_UNIT + + # reduced mass of the nucleus and the WIMP + mr = (m*mn)/(m + mn) + + # reduced mass of the proton and the WIMP + mp = (m*PROTON_MASS)/(m + PROTON_MASS) + + return cs0*A**2*mr**2/mp**2 + +def get_differential_event_rate_xenon(e, m, cs0, A): + """ + Returns the event rate of WIMP scattering in the Xenon 100T detector. + """ + # mass of nucleus + mn = A*ATOMIC_MASS_UNIT + + # reduced mass of the nucleus and the WIMP + mr = (m*mn)/(m + mn) + + cs = get_cross_section(cs0, m, A) + + v0 = DM_VELOCITY + + vesc = ESCAPE_VELOCITY + + # earth's velocity through the galaxy + ve = EARTH_VELOCITY + + # minimum wimp velocity needed to produce a recoil of energy e + vmin = np.sqrt(mn*e/2)*(mn+m)*SPEED_OF_LIGHT/(mn*m) + + f = get_nuclear_form_factor(A, e) + + x = vmin/v0 + y = ve/v0 + z = vesc/v0 + + # Equation 3.49 in Mark Pepin's thesis + k0 = (pi*v0**2)**(3/2) + + # Equation 3.49 in Mark Pepin's thesis + k1 = k0*(erf(z)-(2/np.sqrt(pi))*z*np.exp(-z**2)) + + # From Mark Pepin's CDMS thesis page 59 + if x <= z - y: + I = (k0/k1)*(DM_DENSITY/(2*ve))*(erf(x+y) - erf(x-y) - (4/np.sqrt(pi))*y*np.exp(-z**2)) + elif x <= y + z: + I = (k0/k1)*(DM_DENSITY/(2*ve))*(erf(z) - erf(x-y) - (2/np.sqrt(pi))*(y+z-x)*np.exp(-z**2)) + else: + return 0 + + return (N0/A)*(mn/(2*m*mr**2))*cs*f**2*I*SPEED_OF_LIGHT**2 + +def get_event_rate_xenon(m, cs, A, threshold): + """ + Returns the event rate of WIMP scattering in a dark matter detector using + an element with atomic mass A. Rate is in Hz kg^-1. + """ + # mass of nucleus + mn = A*ATOMIC_MASS_UNIT + + # reduced mass of the nucleus and the WIMP + mr = (m*mn)/(m + mn) + + vesc = ESCAPE_VELOCITY + + # earth's velocity through the galaxy + ve = EARTH_VELOCITY + + # max recoil they looked for was 40 keV + emax = 2*mr**2*((ve+vesc)/SPEED_OF_LIGHT)**2/mn + + rate, err = quad(get_differential_event_rate_xenon, threshold, emax, args=(m,cs,A), epsabs=0, epsrel=1e-2, limit=1000) + + return rate + +if __name__ == '__main__': + import matplotlib.pyplot as plt + + ls = np.logspace(-1,8,1000) + + cs = 1e-50 # cm^2 + + # FIXME: should use water density for L < 1 m, silicon density for L ~ 1 + # km, and iron for L >> 1 km + rate = np.array([get_event_rate_sno(1, cs*1e-4, l, water) for l in ls]) + + colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] + + plt.figure() + plt.subplot(111) + plt.plot(ls, rate/np.max(rate),color=colors[0]) + plt.xlabel("Mediator Decay Length (m)") + plt.ylabel("Event Rate (arbitrary units)") + plt.axvline(x=FIDUCIAL_RADIUS, color=colors[1], ls='--', label="SNO radius") + plt.axvline(x=h, ls='--', color=colors[2], label="Depth of SNO") + plt.axvline(x=R, ls='--', color=colors[3], label="Earth radius") + plt.gca().set_xscale('log') + plt.gca().set_xlim((ls[0], ls[-1])) + plt.title("Event Rate of Self-Destructing Dark Matter in SNO") + plt.legend() + + # threshold is ~5 keVr + xenon_100t_threshold = 5e-6 + + # threshold is ~1 keVr + cdms_threshold = 1e-6 + + # threshold is ~100 eV + # FIXME: is this correct? + cresst_threshold = 1e-7 + + ms = np.logspace(-2,3,200) + cs0s = np.logspace(-50,-40,200) + + mm, cs0cs0 = np.meshgrid(ms, cs0s) + + rate1 = np.empty(mm.shape) + rate2 = np.empty(mm.shape) + rate3 = np.empty(mm.shape) + rate4 = np.empty(mm.shape) + rate5 = np.empty(mm.shape) + rate6 = np.empty(mm.shape) + for i in range(mm.shape[0]): + print("\r%i/%i" % (i+1,mm.shape[0]),end='') + for j in range(mm.shape[1]): + rate1[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, XENON_MASS, xenon_100t_threshold) + rate2[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1.0, water) + rate3[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, GERMANIUM_MASS, cdms_threshold) + rate4[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, TUNGSTEN_MASS, cresst_threshold) + rate5[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1e3, norite) + rate6[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1e6, mantle) + print() + + # Fiducial volume of the Xenon1T detector is 1042 +/- 12 kg + # from arxiv:1705.06655 + xenon_100t_fiducial_volume = 1042 # kg + + # Livetime of the Xenon1T results + # from arxiv:1705.06655 + xenon_100t_livetime = 34.2 # days + + plt.figure() + plt.subplot(111) + plt.gca().set_xscale('log') + plt.gca().set_yscale('log') + CS1 = plt.contour(mm,cs0cs0,rate1*3600*24*xenon_100t_fiducial_volume*xenon_100t_livetime,[10.0], colors=[colors[0]]) + CS2 = plt.contour(mm,cs0cs0,rate2*3600*24*668.8,[10.0],colors=[colors[1]]) + CS3 = plt.contour(mm,cs0cs0,rate3*3600*24*70.10,[10.0],colors=[colors[2]]) + # FIXME: I used 2.39 kg day because that's what CRESST-3 reports in their paper + # but! I only use Tungsten here so do I need to multiply by the mass fraction of tungsten? + CS4 = plt.contour(mm,cs0cs0,rate4*3600*24*2.39,[10.0],colors=[colors[3]]) + CS5 = plt.contour(mm,cs0cs0,rate5*3600*24*668.8,[10.0],colors=[colors[1]], linestyles=['dashed']) + CS6 = plt.contour(mm,cs0cs0,rate6*3600*24*668.8,[10.0],colors=[colors[1]], linestyles=['dotted']) + plt.clabel(CS1, inline=1, fmt="XENON1T", fontsize=10, use_clabeltext=True) + plt.clabel(CS2, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1 m)", fontsize=10) + plt.clabel(CS3, inline=1, fmt="CDMSLite", fontsize=10) + plt.clabel(CS4, inline=1, fmt="CRESST-3", fontsize=10) + plt.clabel(CS5, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1 km)", fontsize=10) + plt.clabel(CS6, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1000 km)", fontsize=10) + plt.xlabel(r"$m_\chi$ (GeV)") + plt.ylabel(r"WIMP-nucleon scattering cross section ($\mathrm{cm}^2$)") + plt.title("Self-Destructing Dark Matter Limits") + + x = np.linspace(0,300e-6,1000) + + # reproducing Figure 2.1 in Tom Caldwell's thesis + plt.figure() + plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(XENON_MASS, x)**2,x)), label="Xe") + plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(NEON_MASS, x)**2,x)), label="Ne") + plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(ARGON_MASS, x)**2,x)), label="Ar") + plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(GERMANIUM_MASS, x)**2,x)), label="Ge") + plt.xlabel("Recoil Energy (keV)") + plt.ylabel(r"$F^2(E)$") + plt.gca().set_yscale('log') + plt.legend() + plt.gca().set_ylim((1e-5,1)) + plt.gca().set_xlim((0,300)) + plt.title("Nuclear Form Factors for various elements") + plt.show() diff --git a/sddm.tex b/sddm.tex new file mode 100644 index 0000000..c7da970 --- /dev/null +++ b/sddm.tex @@ -0,0 +1,500 @@ +\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} |