Solar Reflected Dark Matter (SRDM) event rates from differential rates

[plt109]

What?

Based on the StandardHaloModel class, the introduction of the SolarReflectedDMModel class in this PR allows the user to compute the SRDM differential rates from the differential flux generated by DaMaSCUS-SUN.

Why?

The SRDM differential rates computed using DaMaSCUS-SUN is already summed over all the xenon orbitals. The loss of orbital information means that one will not be able to take advantage of the orbital-dependent bonus quanta when performing Single Electron/S2only type analyses.

Instead of regenerating the spectra using DaMaSCUS-SUN which can be computationally costly, the differential flux as a function of dark matter speeds that is also an output of DaMaSCUS-SUN is injected into wimprates so that the differential event rates is handled by wimprates instead of DaMaSCUS-SUN.

This decoupling also allows the computation of the SRDM differential rates using different xenon ionization form factors.

This feature is already available on a forked version of wimprates, but it will be good to have it on the main repo because the xenon repo pulls from the main repo.

How?

The SRDM differential rates for DM-e interactions is computed using rate_srdm, which is similar to rate_dme that is used to compute the DM-e interaction differential rates for the usual halo dm.

Differential flux computed from DAMASCUS-Sun is equivalent to the following for the non-reflected halo dark matter dphi/dv = rho_mw/mw * v * f_observed

See eq 18 and surrounding text in Phys. Rev. D 105, 063020 (2022) (2102.12483v2) for details.

Testing?

The SRDM signal spectra can be generated using the following code snippet:t

minE = 10 # eV, but don't assign nu units yet
maxE = 600 # eV, but don't assign nu units yet
numE = 200
energies = np.linspace(minE, maxE, numE)

drs = {}
for n,l in wr.dme_shells:
    drs[wr.shell_str(n,l)] = wr.rate_srdm(
        energies*nu.eV, 
        n,l,
        mw=m_DM * nu.GeV/nu.c0**2,
        sigma_dme=s_cm2 * nu.cm**2, 
    )* (nu.kg * nu.eV * nu.day)

Screenshots (optional)

[JelleAalbers]

Thanks Pueh! I agree this would be a nice model to add to wimprates.

I took this opportunity to refactor the wimprates halo model:

• v_max is now a HaloModel method.
• The velocity part of the rate integral (inverse mean speed) is factored out and precomputed once, the first time the velocity distribution of a halo model is called. DM-electron scattering already used this trick, but for ordinary elastic recoils and the Migdal effect we did the integral every time. Factoring it out should greatly speed up these calculations. This changes wimprates' output at the 1e-3 level (probably a reasonable integration error), so I had to slightly adjust the tests.

I changed your SRDMModel into a HaloModel, SolarReflectedDMEModel. Of course, SRDM does not describe halo dark matter. It's "velocity distribution" just gives what you need to put into a rate calculation that assumes halo dark matter: differential flux / shm_number_density / v. This avoids duplicating the rate_dme calculation in rate_srdm, which also makes it easier to put in new form factors later. I think I get the same results (evaluated at just 10 points since I am lazy), but you probably want to check for yourself:

(Aside: adding up rates vs. recoil energy from different shells is potentially misleading, as events from different shells at the same recoil energy will yield quite a different number of produced electrons. But I've also made plots like this in older notes, so I won't complain 😅.)

In the new implementation you can speed up the computation by manually creating an SolarReflectedDMEModel and passing it as a halo_model to rate_dme, e.g.

import numericalunits as nu
import numpy as np
import wimprates as wr

e_rec = np.linspace(1, 600, 10) * nu.eV
mw = 9e-03 * nu.GeV/nu.c0**2
sigma_dme = 1e-38 * nu.cm**2

rates = dict()
halo_model = wr.SolarReflectedDMEModel(mw=mw, sigma_dme=sigma_dme)

for n, l in wr.dme_shells:
    rates[(n,l)] = wr.rate_dme(
        erec=e_rec,
        sigma_dme=sigma_dme,
        mw=mw,
        n=n,
        l=l,
        halo_model=halo_model)

Calling wr_rate_srdm in a loop is slower because you pay the price for precomputing the mean inverse speed for each shell. But wr_rate_srdm ensures that mw matches, the form factor is 1, etc, so it is more fool proof.

One remaining issue is that the user has to guess what (cross-section, mass) combinations we have fluxes for. If you pick an even slightly different cross-section you just get an error.

We should probably interpolate between the available fluxes. Something like the below should work if the speeds at which data is provided were the same for each point. Would that be easy to do? You'd have to provide a lot of points that are zero for some speeds of course. Maybe there is a better way to do the interpolation?

from scipy import interpolate

new_mass = 1.1e-3
new_xsec = 1e-37

keys = list(wr.halo.srdm_fluxes.keys())
# Keys are e.g. "mass9.490e-05_xsec1.100e-35"  Convert to two-tuples:
keys = [tuple([float(x[4:]) for x in k.split('_')]) for k in keys]
xys = np.asarray(keys)
values = np.asarray([
    df['Differential Flux'].values
    for df in wr.halo.srdm_fluxes.values()])
interpolator = interpolate.LinearNDInterpolator(np.log(xys), np.log(values))
result = np.exp(interpolator(np.log(np.array([new_mass, new_xsec]))))[0]

plt.plot(result)

What do you think?

(And I'm not sure why the tests are failing. Googling the error message gives something related to old nbconvert version that we shouldn't be using..)