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

use crate::cosmology::relic::{effective_dof, entropy_dof};
use std::f64::consts::PI;

pub fn comoving_number_density_eq(x: f64, g_dm: f64) -> f64 {
    let factor = g_dm * x.powf(1.5) * (-x).exp() / (2.0 * PI).powf(1.5);
    factor.max(0.0)
}

pub fn collision_term_annihilation(y: f64, y_eq: f64, sigma_v_cm3_s: f64, s_entropy: f64) -> f64 {
    -sigma_v_cm3_s * s_entropy * (y * y - y_eq * y_eq)
}

pub fn boltzmann_rhs(x: f64, y: f64, m_dm_gev: f64, sigma_v_cm3_s: f64, g_dm: f64) -> f64 {
    let t = m_dm_gev / x;
    let g_star_s = entropy_dof(t);
    let mpl = 1.22e19;
    let s = 2.0 * PI * PI / 45.0 * g_star_s * t * t * t;
    let h = (PI * PI * effective_dof(t) / 90.0).sqrt() * t * t / mpl;
    let y_eq = comoving_number_density_eq(x, g_dm) / (s + 1e-100);
    let lambda = sigma_v_cm3_s * s / (h * x + 1e-100);
    -lambda * (y * y - y_eq * y_eq) / x
}

pub fn solve_boltzmann_freeze_out(
    m_dm_gev: f64,
    sigma_v_cm3_s: f64,
    g_dm: f64,
    x_start: f64,
    x_end: f64,
    n_steps: usize,
) -> Vec<(f64, f64)> {
    let dx = (x_end - x_start) / n_steps as f64;
    let t0 = m_dm_gev / x_start;
    let g_star_s0 = entropy_dof(t0);
    let s0 = 2.0 * PI * PI / 45.0 * g_star_s0 * t0 * t0 * t0;
    let y_eq0 = comoving_number_density_eq(x_start, g_dm) / (s0 + 1e-100);
    let mut y = y_eq0;
    let mut x = x_start;
    let mut result = Vec::with_capacity(n_steps + 1);
    result.push((x, y));
    for _ in 0..n_steps {
        let k1 = dx * boltzmann_rhs(x, y, m_dm_gev, sigma_v_cm3_s, g_dm);
        let k2 = dx * boltzmann_rhs(x + 0.5 * dx, y + 0.5 * k1, m_dm_gev, sigma_v_cm3_s, g_dm);
        let k3 = dx * boltzmann_rhs(x + 0.5 * dx, y + 0.5 * k2, m_dm_gev, sigma_v_cm3_s, g_dm);
        let k4 = dx * boltzmann_rhs(x + dx, y + k3, m_dm_gev, sigma_v_cm3_s, g_dm);
        y += (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
        y = y.max(1e-30);
        x += dx;
        result.push((x, y));
    }
    result
}

pub fn relic_density_from_solution(y_final: f64, m_dm_gev: f64, x_final: f64) -> f64 {
    let t_final = m_dm_gev / x_final;
    let g_star_s = entropy_dof(t_final);
    let s_now = 2891.2;
    let s_final = 2.0 * PI * PI / 45.0 * g_star_s * t_final * t_final * t_final;
    let n_now = y_final * s_now;
    let rho_crit_gev_cm3 = 1.054e-5;
    let h2 = 0.49;
    n_now * m_dm_gev / (rho_crit_gev_cm3 * h2) * s_now / (s_final + 1e-100)
}

pub fn freeze_out_x_from_solution(solution: &[(f64, f64)], g_dm: f64, m_dm_gev: f64) -> f64 {
    for (x, y) in solution.iter().skip(1) {
        let t = m_dm_gev / x;
        let g_star_s = entropy_dof(t);
        let s = 2.0 * PI * PI / 45.0 * g_star_s * t * t * t;
        let y_eq = comoving_number_density_eq(*x, g_dm) / (s + 1e-100);
        if *y > 2.5 * y_eq && y_eq > 0.0 {
            return *x;
        }
    }
    solution.last().map_or(25.0, |(x, _)| *x)
}

pub fn coannihilation_effective_sigma_v(
    sigma_11: f64,
    sigma_12: f64,
    sigma_22: f64,
    g1: f64,
    g2: f64,
    delta_m_over_m: f64,
    x: f64,
) -> f64 {
    let r = g2 / g1 * (1.0 + delta_m_over_m).powf(1.5) * (-x * delta_m_over_m).exp();
    let g_eff = 1.0 + r;
    (sigma_11 + 2.0 * r * sigma_12 + r * r * sigma_22) / (g_eff * g_eff)
}

pub fn freeze_in_yield(coupling_sq: f64, m_mediator_gev: f64, m_dm_gev: f64, mpl_gev: f64) -> f64 {
    let g_star = effective_dof(m_mediator_gev);
    let phase_space = (1.0 - (m_dm_gev / m_mediator_gev).powi(2)).max(0.01).sqrt();
    135.0 * coupling_sq * phase_space
        / (1.74 * 8.0 * PI * PI * PI * g_star.powf(1.5) * m_mediator_gev)
        * mpl_gev
}

pub fn freeze_in_relic_density(coupling_sq: f64, m_mediator_gev: f64, m_dm_gev: f64) -> f64 {
    let mpl = 1.22e19;
    let y_fi = freeze_in_yield(coupling_sq, m_mediator_gev, m_dm_gev, mpl);
    let s_now = 2891.2;
    let rho_crit = 1.054e-5;
    y_fi * m_dm_gev * s_now / rho_crit
}

pub fn boltzmann_rhs_decay(
    x: f64,
    y: f64,
    m_parent_gev: f64,
    gamma_decay: f64,
    m_dm_gev: f64,
    g_parent: f64,
) -> f64 {
    let t = m_parent_gev / x;
    let g_star = effective_dof(t);
    let mpl = 1.22e19;
    let h = (PI * PI * g_star / 90.0).sqrt() * t * t / mpl;
    let n_eq = g_parent * (m_parent_gev * t / (2.0 * PI)).powf(1.5) * (-m_parent_gev / t).exp();
    let depletion = (-y * h * x / (n_eq + 1e-300)).exp().min(1.0);
    gamma_decay * n_eq * m_dm_gev * depletion / (h * x * m_parent_gev * 1e-100_f64.max(h))
}

pub fn entropy_injection_ratio(t_reheat_gev: f64, t_decay_gev: f64) -> f64 {
    if t_decay_gev >= t_reheat_gev {
        return 1.0;
    }
    let g_rh = entropy_dof(t_reheat_gev);
    let g_d = entropy_dof(t_decay_gev);
    (g_rh / g_d).powf(1.0 / 3.0) * t_reheat_gev / t_decay_gev
}