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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

use crate::constants::physical::{C, G};
use std::f64::consts::PI;

#[derive(Clone, Debug, PartialEq)]
pub struct SpikeProfile {
    pub m_bh_kg: f64,
    pub rho_0: f64,
    pub r_0: f64,
    pub gamma_sp: f64,
}

pub fn spike_density(r: f64, m_bh_kg: f64, rho_0: f64, r_0: f64, gamma_sp: f64) -> f64 {
    let r_sp = spike_radius(m_bh_kg, rho_0, r_0, gamma_sp);
    if r >= r_sp || r <= 0.0 {
        return 0.0;
    }
    let r_isco = 6.0 * G * m_bh_kg / (C * C);
    if r < r_isco {
        return 0.0;
    }
    rho_0 * (r_0 / r).powf(gamma_sp) * (1.0 - r_isco / r).powi(3)
}

pub fn spike_radius(m_bh_kg: f64, rho_0: f64, r_0: f64, gamma_sp: f64) -> f64 {
    let r_h = G * m_bh_kg / (200.0 * 1e3 * 1e3);
    let factor = m_bh_kg / (4.0 * PI * rho_0 * r_0.powf(gamma_sp));
    factor.powf(1.0 / (3.0 - gamma_sp)).min(r_h)
}

pub fn spike_slope_from_nfw(gamma_init: f64) -> f64 {
    (9.0 - 2.0 * gamma_init) / (4.0 - gamma_init)
}

pub fn spike_slope_from_initial(gamma_init: f64) -> f64 {
    spike_slope_from_nfw(gamma_init)
}

pub fn annihilation_plateau_density(m_dm_gev: f64, sigma_v_cm3_s: f64, t_bh_gyr: f64) -> f64 {
    let m_dm_kg = m_dm_gev * 1.782_661_92e-27;
    let t_s = t_bh_gyr * 3.156e16;
    let sigma_v_si = sigma_v_cm3_s * 1e-6;
    m_dm_kg / (sigma_v_si * t_s)
}

pub fn spike_density_with_annihilation(
    r: f64,
    m_bh_kg: f64,
    rho_0: f64,
    r_0: f64,
    gamma_sp: f64,
    rho_plateau: f64,
) -> f64 {
    let rho_sp = spike_density(r, m_bh_kg, rho_0, r_0, gamma_sp);
    if rho_sp > rho_plateau {
        rho_plateau
    } else {
        rho_sp
    }
}

pub fn annihilation_luminosity_spike(
    profile: &SpikeProfile,
    sigma_v_cm3_s: f64,
    m_dm_gev: f64,
    t_bh_gyr: f64,
    n_steps: usize,
) -> f64 {
    let rho_plat = annihilation_plateau_density(m_dm_gev, sigma_v_cm3_s, t_bh_gyr);
    let r_isco = 6.0 * G * profile.m_bh_kg / (C * C);
    let r_sp = spike_radius(
        profile.m_bh_kg,
        profile.rho_0,
        profile.r_0,
        profile.gamma_sp,
    );
    let sigma_v_si = sigma_v_cm3_s * 1e-6;
    let m_dm_kg = m_dm_gev * 1.782_661_92e-27;
    let log_r_min = r_isco.ln();
    let log_r_max = r_sp.ln();
    let d_log_r = (log_r_max - log_r_min) / n_steps as f64;
    let mut luminosity = 0.0;
    for i in 0..n_steps {
        let log_r = log_r_min + (i as f64 + 0.5) * d_log_r;
        let r = log_r.exp();
        let rho = spike_density_with_annihilation(
            r,
            profile.m_bh_kg,
            profile.rho_0,
            profile.r_0,
            profile.gamma_sp,
            rho_plat,
        );
        luminosity += rho * rho * r * r * r * d_log_r;
    }
    luminosity * 4.0 * PI * sigma_v_si / (2.0 * m_dm_kg * m_dm_kg)
}

pub fn adiabatic_growth_density(
    r_f: f64,
    rho_initial: &dyn Fn(f64) -> f64,
    m_bh_kg: f64,
    m_halo_fn: &dyn Fn(f64) -> f64,
    n_iter: usize,
) -> f64 {
    let mut r_i = r_f;
    for _ in 0..n_iter {
        let m_f = G * m_bh_kg + m_halo_fn(r_f) * G;
        let m_i = m_halo_fn(r_i) * G;
        if m_i < 1e-100 {
            break;
        }
        r_i = r_f * m_f / m_i;
    }
    let rho_i = rho_initial(r_i);
    rho_i * (r_i / r_f).powi(2) * m_halo_fn(r_i) / (m_halo_fn(r_f) + m_bh_kg)
}

pub fn dynamical_relaxation_spike_slope(gamma_sp: f64) -> f64 {
    1.5_f64.min(gamma_sp)
}

pub fn spike_enclosed_mass(
    m_bh_kg: f64,
    rho_0: f64,
    r_0: f64,
    gamma_sp: f64,
    r_max: f64,
    n_steps: usize,
) -> f64 {
    let r_isco = 6.0 * G * m_bh_kg / (C * C);
    let log_r_min = r_isco.ln();
    let log_r_max = r_max.ln();
    let d_log_r = (log_r_max - log_r_min) / n_steps as f64;
    let mut mass = 0.0;
    for i in 0..n_steps {
        let log_r = log_r_min + (i as f64 + 0.5) * d_log_r;
        let r = log_r.exp();
        let rho = spike_density(r, m_bh_kg, rho_0, r_0, gamma_sp);
        mass += rho * r * r * r * d_log_r;
    }
    mass * 4.0 * PI
}

pub fn isco_radius(m_bh_kg: f64) -> f64 {
    6.0 * G * m_bh_kg / (C * C)
}

pub fn schwarzschild_radius(m_bh_kg: f64) -> f64 {
    2.0 * G * m_bh_kg / (C * C)
}

pub fn gravitational_capture_radius(m_bh_kg: f64, v_dm_km_s: f64) -> f64 {
    let v = v_dm_km_s * 1e3;
    2.0 * G * m_bh_kg / (v * v)
}