darkmatter 0.0.1

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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

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

pub fn sommerfeld_coulomb(alpha: f64, v_over_c: f64) -> f64 {
    let eps = PI * alpha / v_over_c;
    if eps < 1e-3 {
        return 1.0 + eps;
    }
    eps / (1.0 - (-2.0 * eps).exp())
}

pub fn sommerfeld_yukawa_s_wave(alpha: f64, m_dm_gev: f64, m_phi_gev: f64, v_over_c: f64) -> f64 {
    let epsilon_v = v_over_c / alpha;
    let epsilon_phi = m_phi_gev / (alpha * m_dm_gev);
    if epsilon_phi < 1e-10 {
        return sommerfeld_coulomb(alpha, v_over_c);
    }
    let a_val = epsilon_v / epsilon_phi;
    let c_val = 6.0 * epsilon_phi / (PI * PI);
    let sinh_a = if a_val > 500.0 {
        return 1.0;
    } else {
        (PI * a_val).sinh()
    };
    let cosh_a = (PI * a_val).cosh();
    let cos_c = (2.0 * PI * (c_val - 0.25 * a_val * a_val).abs().sqrt()).cos();
    PI * sinh_a / (cosh_a - cos_c)
}

pub fn sommerfeld_yukawa_resonant(
    alpha: f64,
    m_dm_gev: f64,
    m_phi_gev: f64,
    v_over_c: f64,
    n: u32,
) -> f64 {
    let epsilon_v = v_over_c / alpha;
    let epsilon_phi = m_phi_gev / (alpha * m_dm_gev);
    let epsilon_res = 6.0 / (PI * PI * n as f64 * n as f64);
    let delta = epsilon_phi - epsilon_res;
    let numerator = epsilon_v.powi(2) + delta * delta;
    if numerator < 1e-100 {
        return 1.0 / (epsilon_v * epsilon_v + 1e-100);
    }
    epsilon_res / (epsilon_v * epsilon_v + delta * delta)
}

pub fn sommerfeld_p_wave(alpha: f64, v_over_c: f64) -> f64 {
    let eps = PI * alpha / v_over_c;
    if eps < 1e-3 {
        return 1.0 + 2.0 * eps * eps / 3.0;
    }
    let s = sommerfeld_coulomb(alpha, v_over_c);
    s * (1.0 + eps * eps / 4.0)
}

pub fn unitarity_bound(m_dm_gev: f64, v_over_c: f64) -> f64 {
    let m_kg = m_dm_gev * 1.782_661_92e-27;
    let p = m_kg * v_over_c * 3e8;
    4.0 * PI * HBAR * HBAR / (p * p)
}

pub fn effective_cross_section(sigma_0_cm2: f64, alpha: f64, v_km_s: f64) -> f64 {
    let v_over_c = v_km_s * 1e3 / 3e8;
    let s = sommerfeld_coulomb(alpha, v_over_c);
    let sigma_max = unitarity_bound_cm2(alpha, v_km_s);
    let sigma_enhanced = sigma_0_cm2 * s;
    sigma_enhanced.min(sigma_max)
}

fn unitarity_bound_cm2(alpha: f64, v_km_s: f64) -> f64 {
    let v = v_km_s * 1e3;
    let m_dm = alpha * 3e8 * HBAR / (v * 1e-15);
    let ub = unitarity_bound(m_dm / 1.782_661_92e-27, v / 3e8);
    ub * 1e4
}

pub fn thermal_average_sommerfeld_coulomb(alpha: f64, x: f64) -> f64 {
    let n_points = 200;
    let v_max = 10.0 / x.sqrt();
    let dv = v_max / n_points as f64;
    let mut integral = 0.0;
    for i in 0..n_points {
        let v = (i as f64 + 0.5) * dv;
        if v < 1e-15 {
            continue;
        }
        let s = sommerfeld_coulomb(alpha, v);
        let weight = (x / PI).sqrt() * 4.0 * PI * v * v * (-x * v * v).exp();
        integral += s * weight * dv;
    }
    integral.max(1.0)
}

pub fn freeze_out_boost(alpha: f64, m_dm_gev: f64, t_fo_gev: f64) -> f64 {
    let x_fo = m_dm_gev / t_fo_gev;
    thermal_average_sommerfeld_coulomb(alpha, x_fo)
}

pub fn velocity_weighted_sigma_v(sigma_0_cm3_s: f64, alpha: f64, v_dispersion_km_s: f64) -> f64 {
    let n_points = 100;
    let v_max = 4.0 * v_dispersion_km_s;
    let dv = v_max / n_points as f64;
    let sigma = v_dispersion_km_s * 1e3 / 3e8;
    let mut integral = 0.0;
    let mut norm = 0.0;
    for i in 0..n_points {
        let v_km_s = (i as f64 + 0.5) * dv;
        let v_oc = v_km_s * 1e3 / 3e8;
        if v_oc < 1e-15 {
            continue;
        }
        let s = sommerfeld_coulomb(alpha, v_oc);
        let mb = 4.0 * PI * v_oc * v_oc * (-v_oc * v_oc / (2.0 * sigma * sigma)).exp();
        integral += s * mb * dv;
        norm += mb * dv;
    }
    if norm < 1e-100 {
        return sigma_0_cm3_s;
    }
    sigma_0_cm3_s * integral / norm
}

pub fn boost_factor_dwarf(alpha: f64, v_dispersion_km_s: f64) -> f64 {
    velocity_weighted_sigma_v(1.0, alpha, v_dispersion_km_s)
}

pub fn boost_factor_cluster(alpha: f64, v_dispersion_km_s: f64) -> f64 {
    velocity_weighted_sigma_v(1.0, alpha, v_dispersion_km_s)
}

pub fn saturation_velocity(alpha: f64, m_dm_gev: f64, m_phi_gev: f64) -> f64 {
    let epsilon_phi = m_phi_gev / (alpha * m_dm_gev);
    alpha * epsilon_phi.sqrt() * 3e8 / 1e3
}