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
136
137

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

pub fn velocity_dependent_cross_section(
    v_km_s: f64,
    sigma_0_cm2_g: f64,
    v_ref_km_s: f64,
    power_index: f64,
) -> f64 {
    sigma_0_cm2_g * (v_km_s / v_ref_km_s).powf(power_index)
}

pub fn yukawa_cross_section(alpha_chi: f64, m_chi_gev: f64, m_phi_gev: f64, v_km_s: f64) -> f64 {
    let v_si = v_km_s * 1e3;
    let c = 3e8;
    let beta = v_si / c;
    let m_chi_kg = m_chi_gev * 1.782_661_92e-27;
    let mu = m_chi_kg / 2.0;
    let hbar = 1.054_571_817e-34;
    let a = m_phi_gev / (2.0 * alpha_chi * m_chi_gev);
    let b = m_chi_gev * beta / (2.0 * m_phi_gev);
    if b < 0.1 {
        let born_factor = a * a / (1.0 + a * a);
        let sigma = 8.0 * PI * alpha_chi * alpha_chi
            / (m_phi_gev * m_phi_gev * 1.782_661_92e-27 * 1.782_661_92e-27)
            * born_factor
            * hbar
            * hbar
            / (v_si * v_si * mu * mu + 1e-100);
        return sigma * 1e4;
    }
    let ln_term = (1.0 + b * b).ln();
    4.0 * PI * alpha_chi * alpha_chi
        / (m_chi_gev * m_chi_gev * 1.782_661_92e-27 * 1.782_661_92e-27 * beta.powi(4) * c.powi(4))
        * (ln_term - b * b / (1.0 + b * b))
        * hbar
        * hbar
        * 1e4
}

pub fn gravothermal_collapse_timescale(
    rho_core_kg_m3: f64,
    sigma_v_km_s: f64,
    r_core_kpc: f64,
    sigma_over_m_cm2_g: f64,
) -> f64 {
    let sigma_v = sigma_v_km_s * 1e3;
    let r_core = r_core_kpc * 3.0857e19;
    let kappa = sigma_over_m_cm2_g * 1e-1 * rho_core_kg_m3 * sigma_v;
    if kappa < 1e-50 {
        return f64::INFINITY;
    }
    150.0 * sigma_v / (G * rho_core_kg_m3 * r_core * kappa)
}

pub fn core_radius_sidm(
    rho_s_kg_m3: f64,
    r_s_kpc: f64,
    sigma_over_m_cm2_g: f64,
    age_gyr: f64,
) -> f64 {
    let t = age_gyr * 3.156e16;
    let sigma_v = (4.0 * PI * G * rho_s_kg_m3 * (r_s_kpc * 3.0857e19).powi(2) / 3.0).sqrt();
    let rate = sigma_over_m_cm2_g * 1e-1 * rho_s_kg_m3 * sigma_v;
    let n_scatter = rate * t;
    r_s_kpc * n_scatter.min(10.0).sqrt() * 0.5
}

pub fn isothermal_core_density(sigma_v_km_s: f64, r_core_kpc: f64) -> f64 {
    let sigma_v = sigma_v_km_s * 1e3;
    let r_core = r_core_kpc * 3.0857e19;
    9.0 * sigma_v * sigma_v / (4.0 * PI * G * r_core * r_core)
}

pub fn sidm_heat_conductivity(
    sigma_over_m_cm2_g: f64,
    rho_kg_m3: f64,
    sigma_v_km_s: f64,
    mean_free_path_kpc: f64,
) -> f64 {
    let sigma_v = sigma_v_km_s * 1e3;
    let mfp = mean_free_path_kpc * 3.0857e19;
    let l_eff = mfp.min(1.0 / (sigma_over_m_cm2_g * 1e-1 * rho_kg_m3 + 1e-100));
    rho_kg_m3 * sigma_v * l_eff * sigma_v * sigma_v / 3.0
}

pub fn mean_free_path_kpc(sigma_over_m_cm2_g: f64, rho_kg_m3: f64) -> f64 {
    let sigma_si = sigma_over_m_cm2_g * 1e-1;
    1.0 / (sigma_si * rho_kg_m3 + 1e-100) / 3.0857e19
}

pub fn scattering_rate_per_gyr(sigma_over_m_cm2_g: f64, rho_kg_m3: f64, sigma_v_km_s: f64) -> f64 {
    let sigma_si = sigma_over_m_cm2_g * 1e-1;
    let sigma_v = sigma_v_km_s * 1e3;
    sigma_si * rho_kg_m3 * sigma_v * 3.156e16
}

pub fn core_collapse_criterion(
    sigma_over_m_cm2_g: f64,
    rho_core_kg_m3: f64,
    sigma_v_km_s: f64,
    t_age_gyr: f64,
) -> bool {
    let rate = scattering_rate_per_gyr(sigma_over_m_cm2_g, rho_core_kg_m3, sigma_v_km_s);
    rate * t_age_gyr > 300.0
}

pub fn evaporation_rate(
    sigma_over_m_cm2_g: f64,
    m_sub_solar: f64,
    rho_host_kg_m3: f64,
    v_rel_km_s: f64,
    r_sub_kpc: f64,
) -> f64 {
    let sigma_si = sigma_over_m_cm2_g * 1e-1;
    let v_rel = v_rel_km_s * 1e3;
    let r_sub = r_sub_kpc * 3.0857e19;
    let m_sub = m_sub_solar * 1.989e30;
    let v_esc = (2.0 * G * m_sub / r_sub).sqrt();
    let frac = (-v_esc * v_esc / (2.0 * v_rel * v_rel)).exp();
    sigma_si * rho_host_kg_m3 * v_rel * PI * r_sub * r_sub * frac
}

pub fn bullet_cluster_offset_kpc(
    sigma_over_m_cm2_g: f64,
    sigma_dm_kg_m2: f64,
    v_collision_km_s: f64,
) -> f64 {
    let tau = sigma_over_m_cm2_g * 1e-1 * sigma_dm_kg_m2 / 1.989e30;
    let v = v_collision_km_s * 1e3;
    tau * v / 3.0857e19 * 1e-6
}

pub fn drag_deceleration(sigma_over_m_cm2_g: f64, rho_target_kg_m3: f64, v_rel_km_s: f64) -> f64 {
    let sigma_si = sigma_over_m_cm2_g * 1e-1;
    sigma_si * rho_target_kg_m3 * (v_rel_km_s * 1e3).powi(2)
}