darkmatter 0.0.2

Dark matter simulation engine — gravitational fields, particle dynamics, halo stability, and cosmological constants
Documentation 
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

use crate::constants::physical::C;
use std::f64::consts::PI;

pub fn reduced_mass(m1: f64, m2: f64) -> f64 {
    m1 * m2 / (m1 + m2)
}

pub fn max_recoil_energy(m_dm: f64, m_nucleus: f64, v: f64) -> f64 {
    let mu = reduced_mass(m_dm, m_nucleus);
    2.0 * mu * mu * v * v / m_nucleus
}

pub fn recoil_spectrum_rate(
    rho_local: f64,
    m_dm: f64,
    m_nucleus: f64,
    sigma_0: f64,
    v_0: f64,
    e_r: f64,
) -> f64 {
    let mu = reduced_mass(m_dm, m_nucleus);
    let v_min = (m_nucleus * e_r / (2.0 * mu * mu)).sqrt();
    if v_min > v_0 * 3.0 {
        return 0.0;
    }
    let n_dm = rho_local / m_dm;
    let form_factor = helm_form_factor_sq(e_r, m_nucleus);
    let a_sq = m_nucleus / (1.66054e-27);
    let prefactor = n_dm * sigma_0 * a_sq * a_sq / (2.0 * mu * mu);
    prefactor * form_factor * velocity_integral(v_min, v_0)
}

pub fn helm_form_factor_sq(e_r: f64, m_nucleus: f64) -> f64 {
    let a_nucleon = m_nucleus / 1.66054e-27;
    let hbar_c = 197.327e-15;
    let r_n = 1.2e-15 * a_nucleon.powf(1.0 / 3.0);
    let s = 0.9e-15;
    let r_eff = (r_n * r_n - 5.0 * s * s).sqrt();
    let q = (2.0 * m_nucleus * e_r).sqrt();
    let qr = q * r_eff / (hbar_c * 1e15 / C);
    if qr < 1e-10 {
        return 1.0;
    }
    let j1 = qr.sin() / (qr * qr) - qr.cos() / qr;
    let ff = 3.0 * j1 / qr * (-q * q * s * s / (2.0 * hbar_c * hbar_c * 1e30 / (C * C))).exp();
    ff * ff
}

fn velocity_integral(v_min: f64, v_0: f64) -> f64 {
    let v_esc = 544e3;
    if v_min >= v_esc {
        return 0.0;
    }
    let k = PI.sqrt() * v_0;
    let erf_term =
        crate::utils::math::erf_approx(v_min / v_0) - crate::utils::math::erf_approx(v_esc / v_0);
    let exp_term = (-v_esc * v_esc / (v_0 * v_0)).exp();
    k * (erf_term.abs() - 2.0 * (v_esc - v_min) / (PI.sqrt() * v_0) * exp_term)
}

pub fn annual_modulation_amplitude(v_earth_orbital: f64, v_sun: f64, v_0: f64) -> f64 {
    v_earth_orbital * v_sun / (v_0 * v_0)
}

pub fn spin_independent_cross_section(
    sigma_p: f64,
    a: f64,
    m_dm: f64,
    m_proton: f64,
    m_nucleus: f64,
) -> f64 {
    let mu_p = reduced_mass(m_dm, m_proton);
    let mu_a = reduced_mass(m_dm, m_nucleus);
    sigma_p * (mu_a / mu_p).powi(2) * a * a
}

pub fn spin_dependent_cross_section(
    sigma_p: f64,
    j_nucleus: f64,
    s_p: f64,
    s_n: f64,
    a_p: f64,
    a_n: f64,
) -> f64 {
    let factor = (j_nucleus + 1.0) / j_nucleus;
    sigma_p * factor * (s_p * a_p + s_n * a_n).powi(2)
}

pub fn event_rate_per_kg_day(
    rho_local_gev_cm3: f64,
    m_dm_gev: f64,
    sigma_cm2: f64,
    a_target: f64,
    v_0_km_s: f64,
) -> f64 {
    let rho = rho_local_gev_cm3 * 1e6;
    let n_dm = rho / m_dm_gev;
    let v_0 = v_0_km_s * 1e5;
    let n_target = 6.022e23 / a_target * 1000.0;
    n_dm * n_target * sigma_cm2 * v_0 * 86400.0
}

pub fn exclusion_limit_sigma(
    n_events: f64,
    exposure_kg_day: f64,
    rho_local_gev_cm3: f64,
    m_dm_gev: f64,
    a_target: f64,
    v_0_km_s: f64,
) -> f64 {
    let rho = rho_local_gev_cm3 * 1e6;
    let n_dm = rho / m_dm_gev;
    let v_0 = v_0_km_s * 1e5;
    let n_target = 6.022e23 / a_target * 1000.0;
    n_events / (exposure_kg_day * n_dm * n_target * v_0 * 86400.0)
}