1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
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}
|