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

use std::f64::consts::PI;

pub fn freeze_out_temperature(m_dm_gev: f64, sigma_v: f64) -> f64 {
    let xf = freeze_out_xf(m_dm_gev, sigma_v);
    m_dm_gev / xf
}

pub fn freeze_out_xf(m_dm_gev: f64, sigma_v: f64) -> f64 {
    let mpl_gev = 1.22e19;
    let g_star = effective_dof(m_dm_gev);
    let c0 = 0.038 * mpl_gev * m_dm_gev * sigma_v / g_star.sqrt();
    let log_c0 = c0.ln();
    let mut xf = log_c0;
    for _ in 0..20 {
        let new_xf = log_c0 - 0.5 * xf.ln();
        if (new_xf - xf).abs() < 0.001 {
            break;
        }
        xf = new_xf;
    }
    xf.max(1.0)
}

pub fn relic_density_omega_h2(m_dm_gev: f64, sigma_v: f64) -> f64 {
    let xf = freeze_out_xf(m_dm_gev, sigma_v);
    let mpl = 1.22e19;
    let g_star_s = effective_dof(m_dm_gev / xf);
    let numerator = 1.07e9 * xf;
    let denominator = mpl * g_star_s.sqrt() * sigma_v * 1e36;
    numerator / denominator
}

pub fn thermal_cross_section_for_omega(omega_h2: f64) -> f64 {
    3e-26 * 0.12 / omega_h2
}

pub fn effective_dof(t_gev: f64) -> f64 {
    if t_gev > 300.0 {
        106.75
    } else if t_gev > 100.0 {
        96.25
    } else if t_gev > 10.0 {
        86.25
    } else if t_gev > 1.0 {
        75.75
    } else if t_gev > 0.2 {
        61.75
    } else if t_gev > 0.1 {
        17.25
    } else if t_gev > 1e-3 {
        10.75
    } else if t_gev > 1e-4 {
        3.91
    } else {
        3.36
    }
}

pub fn entropy_dof(t_gev: f64) -> f64 {
    effective_dof(t_gev)
}

pub fn number_density_equilibrium(m_gev: f64, t_gev: f64, g_internal: f64) -> f64 {
    let x = m_gev / t_gev;
    if x > 30.0 {
        return 0.0;
    }
    if x < 0.1 {
        return g_internal * 1.2 * t_gev * t_gev * t_gev / (PI * PI) * 1e27;
    }
    g_internal * (m_gev * t_gev / (2.0 * PI)).powf(1.5) * (-x).exp() * 1e27
}

pub fn entropy_density(t_gev: f64) -> f64 {
    let gs = entropy_dof(t_gev);
    2.0 * PI * PI / 45.0 * gs * t_gev * t_gev * t_gev
}

pub fn yield_equilibrium(m_gev: f64, t_gev: f64, g_internal: f64) -> f64 {
    let n = number_density_equilibrium(m_gev, t_gev, g_internal);
    let s = entropy_density(t_gev);
    if s > 0.0 { n / s } else { 0.0 }
}

pub fn boltzmann_suppression(m_gev: f64, t_gev: f64) -> f64 {
    (-m_gev / t_gev).exp()
}

pub fn yield_today_from_freeze_out(m_gev: f64, sigma_v: f64) -> f64 {
    let xf = freeze_out_xf(m_gev, sigma_v);
    let s0 = entropy_density(1e-13);
    let mpl = 1.22e19;
    let gs = entropy_dof(m_gev / xf);
    xf / (0.264 * mpl * gs.sqrt() * sigma_v * s0)
}

pub fn coannihilation_effective_sigma_v(
    sigmas: &[(f64, f64, f64)],
    m_lightest: f64,
    t_gev: f64,
) -> f64 {
    let mut numerator = 0.0;
    let mut denominator = 0.0;
    for &(gi, mi, sigma_vi) in sigmas {
        let weight = gi * (-(mi - m_lightest) / t_gev).exp();
        numerator += weight * sigma_vi;
        denominator += weight;
    }
    if denominator > 0.0 {
        numerator / denominator
    } else {
        0.0
    }
}