First direct detection constraint on mirror dark matter kinetic mixing using LUX 2013 data

TL;DR

This paper reports a direct detection search for mirror dark matter using LUX 2013 data, setting new constraints on the kinetic mixing parameter and excluding certain mirror electron temperature ranges.

Contribution

It provides the first direct detection constraints on mirror dark matter kinetic mixing using liquid xenon data, accounting for Earth shielding and modulation effects.

Findings

No evidence of mirror dark matter interactions was found.

Established new upper limits on the kinetic mixing parameter.

Excluded mirror electron temperatures above 0.3 keV at 90% confidence.

Abstract

We present the results of a direct detection search for mirror dark matter interactions, using data collected from the Large Underground Xenon experiment during 2013, with an exposure of 95 live-days $\times$ 118 kg. Here, the calculations of the mirror electron scattering rate in liquid xenon take into account the shielding effects from mirror dark matter captured within the Earth. Annual and diurnal modulation of the dark matter flux and atomic shell effects in xenon are also accounted for. Having found no evidence for an electron recoil signal induced by mirror dark matter interactions we place an upper limit on the kinetic mixing parameter over a range of local mirror electron temperatures between 0.1 and 0.6 keV. This limit shows significant improvement over the previous experimental constraint from orthopositronium decays and significantly reduces the allowed parameter space for the model. We exclude mirror electron temperatures above 0.3 keV at a 90% confidence level, for this model, and constrain the kinetic mixing below this temperature.

Tables

Table 1. Table 1 : Nuisance parameters used in the PLR test for a local mirror electron temperature 0.3 keV. The means and standard deviations of the Gaussian constraints are shown along with the value from the best fit to data.

Parameter	Constraint	Fit Value
Low-z-origin $\gamma$ counts	157 $\pm$ 78	160 $\pm$ 17
Other $\gamma$ counts	217 $\pm$ 108	179 $\pm$ 18
$\beta$ counts	65 $\pm$ 32	116 $\pm$ 17
127Xe counts	35 $\pm$ 18	41 $\pm$ 8
37Ar counts	10 $\pm$ 5	10 $\pm$ 7

Equations

$$\begin{array}{rl}\mathcal{L}=& {\mathcal{L}}_{SM}(e,u,d,\gamma ,W,Z,...)+\\ & {\mathcal{L}}_{SM}(e^{\prime },u^{\prime },d^{\prime },{\gamma }^{\prime },W^{\prime },Z^{\prime },...)+{\mathcal{L}}_{mix},\end{array}$$

$${\mathcal{L}}_{mix}=\frac{\epsilon }{2}{F}^{\mu \nu }{F}_{\mu \nu }^{\prime }+\lambda {\varphi }^{†}\varphi {\varphi }^{{}^{\prime }†}{\varphi }^{{}^{\prime }}.$$

$$T=\frac{\overline{m}{v}_{rot}^2}{2},$$

$$\frac{d\sigma }{d{E}_R}=\frac{\lambda }{E_R^2{v}^2},\lambda =\frac{2\pi {\epsilon }^2{\alpha }^2}{m_e}.$$

$$\begin{array}{rl}\frac{dR}{d{E}_R}=g_T{N}_T{n}_{e^{\prime }}^0\frac{\lambda }{v_c^0{E}_R^2}[1& +A_v\cos \omega (\mathrm{t}-{\mathrm{t}}_0)\\ & +A_{\theta }(\theta -\bar{\theta })].\end{array}$$

$$v_{cut}^4\approx \frac{16\pi }{m_e^2}{\alpha }^2\Sigma \mathrm{log\Lambda },$$

$$\Sigma (\psi )=\int n_{H{e}^{\prime }}\mathrm{dl}.$$

$$\frac{1}{v_c^0}=\frac{1}{N{v}_0\sqrt{\pi }}\int e^{-y^2/v_0^2}\mathrm{dcos}\psi ,$$

$$N={\int }_{|v|>v_{cut}}^{\infty }\frac{e^{-v^2/v_0^2}}{v_0^3{\pi }^{3/2}}{d}^{3}v.$$

$$n_{e^{\prime }}^0=N{n}_{e^{\prime }}^{far},$$

$$E=W(n_{e}+n_{\gamma})=W\bigg{(}\frac{S1_{c}}{g_{1}}+\frac{S2_{c}}{g_{2}}\bigg{)},$$

$$\varepsilon(90\%CL)=\varepsilon(0)\Bigg{(}\frac{nSig(90\%CL)}{nPDF(0)}\Bigg{)}^{\frac{1}{2}},$$

Full text

Introduction — The Standard Model (SM) is a gauge field theory with $SU(3)_{c}\bigotimes SU(2)\bigotimes U(1)$ gauge symmetry. It successfully describes known particles and their non-gravitational interactions, but does not contain a suitable dark matter candidate. One possibility for accommodating dark matter particles is that they exist in a hidden sector — a collection of particles and fields which do not interact via SM gauge boson forces, but do interact with SM particles gravitationally  Feng et al. (2008). Mirror dark matter is a special case where the hidden sector is exactly isomorphic to the SM Foot (2014), having the same gauge symmetry. Therefore, it contains mirror partners (denoted ′) of the SM particles with the same masses, lifetimes and self interactions. The full Lagrangian may be written as:

[TABLE]

where $\mathcal{L}_{SM}(e,...)$ and $\mathcal{L}_{SM}(e^{\prime},...)$ are the Langrangians for the SM and mirror sectors, respectively. The two sectors are related by a discrete $Z_{2}$ symmetry transformation, with the only allowed non-gravitational interactions given by:

[TABLE]

Here, the first term describes kinetic mixing of $U(1)_{Y}$ and mirror $U(1)_{Y}^{\prime}$, with field strength tensors $F_{\mu\nu},F_{\mu\nu}^{\prime}$ and kinetic mixing strength $\varepsilon$ Foot et al. (1991). The second term describes Higgs ($\phi$) – mirror Higgs ($\phi^{\prime}$) mixing, with strength determined by parameter $\lambda$. Kinetic mixing induces tiny ordinary electric charges, $\pm\varepsilon e$ for the mirror protons and electrons Holdom (1986). This allows very weak electromagnetic interactions between mirror and SM particles. The kinetic mixing parameter, $\varepsilon$, determines the strength of most mirror – SM particle couplings and is thus the target of experimental searches. The Higgs – mirror Higgs portal can be probed at colliders, through Higgs production and decays, but does not give observable signals in direct detection experiments Foot (2014).

Within the mirror dark matter model kinetic mixing is constrained theoretically to lie in the range; $10^{-11}\leq\varepsilon\leq 4\times 10^{-10}$ Foot (2014). In order for the mirror dark matter halo to be in equilibrium, heating from supernovae must balance energy loss from dissipative processes, giving the lower limit on $\varepsilon$ Foot and Volkas (2004). But if $\varepsilon$ is too high cosmic structure formation would be too heavily damped, giving the upper limit Foot and Vagnozzi (2016).

LUX Experiment — The Large Underground Xenon (LUX) experiment was a dual phase (liquid-gas) time projection chamber (TPC), containing a 250 kg active mass of liquid xenon. The main aim of LUX was to search for dark matter in the form of weakly interacting massive particles (WIMPs), placing limits on spin-independent WIMP-nucleon cross-sections for WIMP masses above 4 GeV Akerib et al. (2014, 2016a). Other studies include searches for spin-dependent WIMP-nucleon interactions Akerib et al. (2017a), electron recoil searches for solar axions and axionlike particles Akerib et al. (2017b) and sub GeV dark matter via the Bremsstrahlung and Migdal effects Akerib et al. (2019).

As described in Ref. Akerib et al. (2013), the LUX TPC was located in a low-radioactivity titanium cryostat, itself within a 6.1 m high, 7.6 m diameter water tank 1458 m underground at the Sanford Underground Research Facility, Lead, USA. Details of the detector calibration and performance are available in Ref. Akerib et al. (2018a). When a particle interacts in the liquid xenon, prompt scintillation photons (S1) and ionisation electrons are produced. The ionisation electrons are drifted upwards by a vertical electric field and extracted into the gas phase, where they produce an electroluminescence signal (S2). Photons from these signals are detected by two arrays of 61 photomultiplier tubes, above and below the active volume. The (x,y) position is obtained from the S2 light distribution in the top PMT array and the depth is found from the delay of the S2 relative to the S1 Akerib et al. (2018b), allowing for fiducialisation of the active volume.

The data used in this analysis was collected between 24th April and 1st September 2013, giving 118 kg $\times$ 95 live days total exposure. Four detector observables are used — $r,z,S1_{c},S2_{c}$, where $S1_{c}$ and $S2_{c}$ refer to amplitudes corrected to equalize the detector response throughout the active volume.

Signal Model — Mirror dark matter would exist as a multi-component plasma halo, assuming that the mirror electron temperature exceeds the binding energy of a mirror hydrogen atom and the cooling time exceeds the Hubble time Clarke and Foot (2016). This halo is predominantly composed of mirror electrons, $\rm e^{\prime}$, and mirror helium nuclei, $\rm He^{\prime}$. The $\rm He^{\prime}$ mass fraction is higher (and $\rm H^{\prime}$ lower) than for ordinary matter because freeze out happens earlier, due to a lower initial temperature in the mirror sector Foot (2014). Kinetic mixing allows electromagnetic interactions between mirror and SM particles, meaning that mirror electrons in the halo can scatter off Xe atomic electrons in the LUX detector.

For a dark matter halo in hydrostatic equilibrium, the local mirror electron temperature is given by  Foot and Volkas (2004):

[TABLE]

where $\overline{m}$ is the average mass of halo particles and $v_{rot}$ is the galactic rotational velocity. Arguments from early universe cosmology in the mirror model give a mirror helium mass fraction of 90% Ciarcelluti and Foot (2010) and assuming a completely ionized plasma $\overline{m}\approx 1.1$ GeV. Therefore, using $v_{rot}\approx 220$ kms*-1* and assuming the halo is in hydrostatic equilibrium, local mirror electron temperature $\sim 0.3$ keV is expected.

In such plasma dark matter models, it is important to consider capture of the dark matter by the Earth Foot (2019). Mirror dark matter is captured when it loses energy due to kinetic mixing interactions with normal matter. Once a significant amount has accumulated, further capture occurs due to mirror dark matter self interactions. Subsequently, mirror dark matter will thermalize with normal matter in the Earth to form an extended distribution, which can affect the incoming mirror dark matter via collisional shielding or deflection by a dark ionosphere. Interactions with the dark ionosphere are very difficult to model  Clarke and Foot (2016), but the collisional shielding, due to mirror particle interactions identical to the standard model version, can be accounted for. Here we follow the formalism presented in Ref. Foot (2019); Clarke and Foot (2016, 2017), first validating the calculations for NaI (as given in Foot (2019)) then performing the calculations for Xe.

The electron – mirror electron Coulomb scattering cross section for this process is given by Clarke and Foot (2016):

[TABLE]

Here $E_{R}$ is electron recoil energy, $v$ velocity of the incoming mirror electron, $m_{e}$ electron mass, $\varepsilon$ the kinetic mixing parameter and $\alpha$ the fine structure constant. The scattering rate, calculated by multiplying with the integral of the velocity distribution of the incoming mirror dark matter and Taylor expanding around the yearly average, is given by Foot (2019):

[TABLE]

Here $N_{T}$ is the number of target electrons, $n_{e^{\prime}}^{0}$ the number density of mirror electrons arriving at the detector and $v_{c}^{0}$ describes the modified velocity distribution at the detector due to shielding. The effective number of free electrons, $g_{T}$, is the number of electrons per target atom with atomic binding energy ($E_{b}$) less than recoil energy ($E_{R}$) — modelled as a step function for the atomic shells in xenon.

The $A_{v}\rm cos\omega(t-t_{0})$ term describes annual modulation resulting from the change of velocity of the Earth with respect to the dark matter halo. Here $\omega=2\pi/$year, $t_{0}=153$ days (2nd June) and modulation amplitude $A_{v}=0.7$ Foot (2019). The $A_{\theta}(\theta-\bar{\theta})$ term describes diurnal and annual modulation due to the rotation of the Earth and the variation of the Earth’s spin axis relative to the incoming dark matter wind. Here $\theta$ is the angle between the halo wind and the zenith at the detector location, $\bar{\theta}$ is the yearly average and the amplitude is $A_{\theta}=1$. The time variation of $\theta$ is examined in Clarke and Foot (2016). The mean modulation terms over the data taking period, accounting for the live time per day, are $A_{v}\rm\langle cos\omega(t-t_{0})\rangle=0.556$ and $A_{\theta}\langle\theta-\bar{\theta}\rangle=0.015$.

Equation 4 shows that $d\sigma/dE_{R}\propto 1/v^{2}$, so the collision length $\propto v^{2}$. This means that for sufficiently large incoming velocity, the effect of collisions becomes negligible (as scattering length exceeds the available distance). Therefore, above some cutoff velocity, $v_{cut}$, collisions do not need to be considered. Below this velocity collisions are important until mirror electron energy is reduced to $\sim$ 25 eV, after which energy loss to the captured mirror helium is no longer important. From energy loss considerations the cutoff velocity may be estimated as Foot (2019):

[TABLE]

where $\Lambda\sim T/E_{min}\approx$ 20, with minimum collisional energy loss $E_{min}$. Column density, $\Sigma$, is calculated by integrating the number density of captured mirror helium nuclei over the path of the incoming mirror dark matter particle:

[TABLE]

Here $\psi$ is the angle between the direction of the incoming mirror electron and the zenith at the detectors location and $l$ is the distance travelled.

The energy dependent term describing the velocity distribution is given by Foot (2019):

[TABLE]

where $v_{0}=\sqrt{2T/m_{e}}$ is the velocity dispersion. Dependence on recoil energy is through $y=MAX[v_{cut}(\psi),v_{min}(E_{R})]$, where $v_{min}(E_{R})=\sqrt{2E_{R}/m_{e}}$ is the minimum velocity needed to produce a recoil of energy $E_{R}$.

The dependence of $v_{c}^{0}$ on recoil energy is shown in Fig. 1. At low values of $E_{R}$ the average velocity exceeds the minimum, $|v|\gg v_{min}$, so most particles can produce recoils with energy $E_{R}$ and the integral becomes independent of $v_{min}$. For large $E_{R}$ the average particle velocity is lower than $v_{min}$, so the integral is suppressed, leading to a sharp rise in $v_{c}^{0}$.

The normalization, $N$, is given by:

[TABLE]

The number density of the high velocity component which arrives at the Earth is given by:

[TABLE]

where $n_{e^{\prime}}^{far}$ = 0.2 cm*-3* is the number density far from the Earth  Clarke and Foot (2017).

Both $v_{c}^{0}$ and $n_{e^{\prime}}^{0}$ depend on the mirror helium density at the Earth’s surface, $n_{He^{\prime}}(R_{E})$ (through column density), which is set to $n_{He^{\prime}}=5.8\times 10^{-11}\rm{cm}^{-3}$ Foot (2019). There is also dependence on electron recoil energy, $E_{R}$ (through $v_{min}$) and mirror electron temperature, $T$ (through $v_{0}$). Substituting Eq. 8 and Eq. 10 into Eq. 5 to calculate differential rate introduces dependence on the kinetic mixing parameter, $\varepsilon$ (through $\lambda$) and the target material (through $N_{T}$ and $g_{T}$). Calculation of the target independent parts $v_{c}^{0}$ and $n_{e^{\prime}}^{0}$ was validated by evaluating the differential rate for NaI. This was convolved with the expected detector resolution, assumed to be Gaussian with energy dependent width Bernabei et al. (2008), in order to reproduce Fig.4(a) from Ref.Foot (2019).

The differential rate of electron recoils in xenon could then be calculated using Eq. 5. If the shielding effects are not accounted for a Maxwellian velocity distribution is assumed for the mirror electrons, with the rate given by Eq.(6.4) of Ref. Clarke and Foot (2016). The differential energy spectra of electron recoils, calculated both with and without the shielding effects, are shown in Fig. 2 for a range of local mirror electron temperatures.

The low energy electron recoil response of the LUX detector was characterised using an internal tritium calibration, as described in Akerib et al. (2016b). The injection of tritiated methane into the gas circulation gave a large sample of electron recoils from beta decays in the energy range of interest, used to precisely measure light and charge yields in the detector. These yields show good agreement with the Noble Element Simulation Technique (NEST) package v2.0 Szydagis et al. (2018). Here we use NEST to model the distributions of the detector observables $r,z,S1_{c},S2_{c}$, taking into account the detector resolution and efficiency, for signal events simulated using the above energy spectra. The quantities $S1_{c}$ and $S2_{c}$ are measured in photons detected (phd), with the resulting distribution in $\rm{log}_{10}$ $S2_{c}$ vs. $S1_{c}$ is shown in Fig. 3(a), for mirror electron temperature $T=0.3$ keV and kinetic mixing $\varepsilon=10^{-10}$.

Background Model — Interactions of mirror dark matter particles within LUX would induce isolated low energy electron recoil events. Consequently, the signal being searched for competes with background events that arise from: Compton scattering of $\gamma$ rays from radioactive decay of isotopes in detector components, $\beta$ decay from 85mKr and Rn contaminants in the liquid xenon and X-rays following 127Xe electron capture where the coincident $\gamma$ ray escapes detection Akerib et al. (2015). Heavily down scattered decays from 238U chain, 232Th chain and 60Co generate additional $\gamma$ rays from the centre of a large copper block below the PMTs. The $\gamma$ rays can be modelled as two separate spatial distributions — one from below the bottom PMT array and one from the rest of the detector. Decays of 37Ar, by electron capture, within the fiducial volume are also included Akerib et al. (2016a). A fiducial radius of 18 cm is used to exclude low energy events from 210Pb on the detector walls. The full background model used in this analysis is shown in Fig. 3(b), with each component normalized to the expected value.

Data Analysis — A series of analysis cuts are applied to the data; events must also come from within a fiducial radius of 18 cm and z range of 8.5–48.6 cm above the bottom PMT array (drift time 305–38 $\mu$s). The S1 pulses in this analysis were required to have two PMTs in coincidence — at least two non adjacent PMTs must measure an integrated area exceeding 0.3 phd. This is imposed to prevent spontaneous photocathode emission from being misidentified as an S1 pulse, as discussed in Ref.  Akerib et al. (2018a). We also require $S1_{c}$ size 1–80 detected photons and the raw $S2$ size to exceed 165 detected photons. Corrected signal amplitudes $S1_{c}$, $S2_{c}$, account for non uniform temporal and spatial response throughout the detector, based on 83mKr calibrations. Position corrections mean that it is possible to have an S1 size below 2 phd, despite this two fold coincidence requirement. The data cuts leave 516 events in our region of interested, shown in Fig. 4 along with 90% signal contours. It should be noted that the signal model is not completely symmetric in $\rm{log_{10}}$ $S2_{c}$, so the contour containing 90% of the signal will not be exactly centred on the ER band. This is a threshold effect due to the exponential shape of the signal model and is more pronounced for the sharply peaked signal models with no shielding.

The energy deposited by an event is given by Aprile and Doke (2010):

[TABLE]

where $n_{e}$ and $n_{\gamma}$ are the number of electrons and photons produced, respectively and $W=(13.7\pm 0.2)$ eV is the work function for producing these quanta in liquid xenon. Gain factors $g_{1}=0.117\pm 0.003$ phd/photon and $g_{2}=12.1\pm 0.8$ phd/electron were determined from calibrations Akerib et al. (2017c).

Compatibility with the data is tested using a two sided profile likelihood ratio test with four physics observables; $S1_{c}$, $\rm{log_{10}}$ $S2_{c}$, $r$, $z$ Cowan et al. (2011). Simulated distributions of the signal model and background model were generated for each observable. The distribution of the test statistic, the ratio of the conditional maximum likelihood (with number of signal events fixed) to the global maximum likelihood, is found for a range of numbers of signal events. This is used to calculate the p-value for each number of signal events. The hypothesis test is then inverted to find the 90% confidence limit on the number of signal events observed in the data. Systematic uncertainties in the background rates are treated as nuisance parameters. As detailed in Ref. Akerib et al. (2015), an extensive screening campaign gave the radioactive content of detector components, which was further constrained using data. Internal backgrounds were estimated from direct measurements of LUX data and sampling the Xe during the run. These were used to project the background rates for the period of data taking and normalize the Monte Carlo spectra. Nuisance parameters had the estimated rate as the mean value with a Gaussian constraint from the uncertainty. The best fit model covers zero signal model contribution for all mirror electron temperatures. The input and fit value for each nuisance parameter is shown in Table 1, giving a total of $506\pm 32$ background events, compared to 516 events in the data. For $T=0.3$ keV, the background-only model gives KS test p-values of 0.27, 0.68, 0.71 and 0.60 for the projected distributions in $S1_{c}$, $\rm{log_{10}}$ $S2_{c}$, $r$ and $z$, respectively. For $T=0.3$ keV this results in a 90% confidence limit of 11 signal events, although it should be noted that the background events extend over a larger energy range than the signal.

The 90% confidence limit on kinetic mixing parameter is then calculated using:

[TABLE]

where $\varepsilon(0)$ is the arbitrary value of $\varepsilon$ used to generate the signal model, $nPDF(0)$ is the corresponding number of signal events and $nSig(90\%CL)$ is the 90% confidence limit on the number of signal events. The power of $1/2$ comes from the dependence of the rate on $\varepsilon^{2}$ in Eq. 4.

Results — We set a 90% confidence limit on the kinetic mixing parameter, $\varepsilon$, for the local mirror electron temperature range 0.1-0.9 keV, as shown in Fig. 5. The previous experimental constraint on $\varepsilon$ comes from invisible decays of orthopositronium in a vacuum Vigo et al. (2018).

If positronium – mirror positronium mixing were to occur, decay to missing photons would leave a missing energy signal. The upper limit placed on the branching fraction of orthopositronium to invisible states gives a 90% upper confidence limit on the kinetic mixing parameter of $\varepsilon\leq 3.1\times 10^{-7}$. The astrophysical constraint on kinetic mixing within the mirror dark matter theory; $10^{-11}\leq\varepsilon\leq 4\times 10^{-10}$ Foot (2014), is also shown.

In Ref. Aprile et al. (2015), the XENON100 collaboration examine the possibility of leptophilic dark matter models explaining the DAMA Bernabei et al. (2013) modulation signal. For each model the expected signal in xenon, given the DAMA modulation amplitude, is compared to XENON100 electron reocil data. This ruled out mirror dark matter as an explanation at a 3.6$\sigma$ confidence level, but there was no explicit search for mirror dark matter and no constraint was placed on the model itself.

Conclusion/Summary — We have presented the results of the first dedicated direct detection search for mirror dark matter. The effect of mirror dark matter capture by the Earth and subsequent shielding is included, for the first time, for a signal in Xe. A significant proportion of the parameter space allowed by the theory is excluded by this analysis. However the present theoretical treatment makes assumptions for the local mirror electron temperature (thermal equilibrium with nuclei in the halo) and density Clarke and Foot (2016, 2017). The effect of deflection by the captured dark ionosphere is not included and this could significantly alter the signal model. Furthermore, the extent of these shielding effects may have significant dependence on the detector elevation relative to sea level, if the captured distribution is assumed to be spherically symmetric.

Whilst there are possible caveats and extensions to this conceptually simple but phenomenologically complex mirror dark matter model, we have set limits based on the current model. This shows that it is possible to use direct detection experiments to probe low mass particles in a hidden sector.

Acknowledgements — The authors would like to thank Robert Foot for helpful correspondence.

This Letter was partially supported by the U.S. Department of Energy (DOE) under Award No. DE-AC02-05CH11231, DE-AC05-06OR23100, DE-AC52-07NA27344, DE-FG01-91ER40618, DE-FG02-08ER41549, DE-FG02-11ER41738, DE-FG02-91ER40674, DE-FG02-91ER40688, DE-FG02-95ER40917, DE-NA0000979, DE-SC0006605, DE-SC0010010, DE-SC0015535, and DE-SC0019066; the U.S. National Science Foundation under Grants No. PHY-0750671, PHY-0801536, PHY-1003660, PHY-1004661, PHY-1102470, PHY-1312561, PHY-1347449, PHY-1505868, and PHY-1636738; the Research Corporation Grant No. RA0350; the Center for Ultra-low Background Experiments in the Dakotas (CUBED); and the South Dakota School of Mines and Technology (SDSMT). Laboratório de Instrumentação e Física Experimental de Partículas (LIP)-Coimbra acknowledges funding from Fundação para a Ciência e a Tecnologia (FCT) through the Project-Grant PTDC/FIS-NUC/1525/2014. Imperial College and Brown University thank the UK Royal Society for travel funds under the International Exchange Scheme (IE120804). The UK groups acknowledge institutional support from Imperial College London, University College London and Edinburgh University, and from the Science & Technology Facilities Council for for PhD studentships R504737 (EL), M126369B (NM), P006795 (AN), T93036D (RT) and N50449X (UU). This work was partially enabled by the University College London (UCL) Cosmoparticle Initiative. The University of Edinburgh is a charitable body, registered in Scotland, with Registration No. SC005336.

This research was conducted using computational resources and services at the Center for Computation and Visualization, Brown University, and also the Yale Science Research Software Core.

We gratefully acknowledge the logistical and technical support and the access to laboratory infrastructure provided to us by SURF and its personnel at Lead, South Dakota. SURF was developed by the South Dakota Science and Technology Authority, with an important philanthropic donation from T. Denny Sanford. Its operation is funded through Fermi National Accelerator Laboratory by the Department of Energy, Office of High Energy Physics.