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

use crate::constants::cosmic::HUBBLE_CONSTANT_SI;
use crate::constants::physical::{G, HBAR};
use std::f64::consts::PI;

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct WarmDarkMatter {
    pub mass_kev: f64,
    pub temperature_ratio: f64,
}

impl WarmDarkMatter {
    pub fn mass_kg(&self) -> f64 {
        self.mass_kev * 1e3 * 1.782_661_92e-36
    }

    pub fn free_streaming_scale_mpc(&self) -> f64 {
        0.2 * (1.0 / self.mass_kev) * (self.temperature_ratio / 1.0).powf(1.0 / 3.0)
    }

    pub fn free_streaming_mass_solar(&self) -> f64 {
        let lambda_fs = self.free_streaming_scale_mpc() * 3.0857e22;
        let rho_crit = 3.0 * HUBBLE_CONSTANT_SI * HUBBLE_CONSTANT_SI / (8.0 * PI * G);
        let omega_dm = 0.265;
        4.0 / 3.0 * PI * (lambda_fs / 2.0).powi(3) * rho_crit * omega_dm / 1.989e30
    }

    pub fn thermal_velocity_today_km_s(&self) -> f64 {
        0.012 * (1.0 / self.mass_kev) * self.temperature_ratio
    }

    pub fn phase_space_density_bound(&self, rho_gev_cm3: f64, sigma_km_s: f64) -> f64 {
        let rho_si = rho_gev_cm3 * 1.782_661_92e-27 * 1e6;
        let sigma_si = sigma_km_s * 1e3;
        let m = self.mass_kg();
        rho_si / (m * sigma_si * sigma_si * sigma_si)
    }
}

pub fn transfer_function_wdm(k_mpc: f64, mass_kev: f64, omega_dm: f64, h: f64) -> f64 {
    let alpha =
        0.049 * (mass_kev / 1.0).powf(-1.11) * (omega_dm / 0.25).powf(0.11) * (h / 0.7).powf(1.22);
    let ratio = alpha * k_mpc;
    (1.0 + (ratio).powf(2.0 * 1.12)).powf(-5.0 / (1.12 * 2.0))
}

pub fn half_mode_mass_solar(mass_kev: f64) -> f64 {
    2.7e8 * (mass_kev / 1.0).powf(-3.33)
}

pub fn suppression_power_spectrum(k_mpc: f64, mass_kev: f64, omega_dm: f64, h: f64) -> f64 {
    let t = transfer_function_wdm(k_mpc, mass_kev, omega_dm, h);
    t * t
}

pub const LYMAN_ALPHA_LOWER_BOUND_KEV: f64 = 3.3;

pub fn tremaine_gunn_bound_ev(sigma_km_s: f64, r_kpc: f64) -> f64 {
    let sigma_si = sigma_km_s * 1e3;
    let r_si = r_kpc * 3.0857e19;
    let bound = (9.0 * PI * HBAR * HBAR * HBAR
        / (G.sqrt() * sigma_si * r_si * r_si * 1.782_661_92e-36))
        .powf(0.25);
    bound * 1e36 / 1.782_661_92
}

pub fn wdm_halo_mass_function_suppression(m_solar: f64, half_mode_mass: f64) -> f64 {
    let ratio = half_mode_mass / m_solar;
    (1.0 + 2.7 * ratio).powf(-0.99)
}

pub fn sterile_neutrino_mixing_angle(mass_kev: f64, omega_dm: f64) -> f64 {
    2.3e-5 * (omega_dm / 0.12) * (1.0 / mass_kev).powf(1.8)
}

pub fn xray_flux_from_decay(
    mass_kev: f64,
    sin2_2theta: f64,
    distance_mpc: f64,
    dm_mass_solar: f64,
) -> f64 {
    let gamma = 1.38e-29 * sin2_2theta * (mass_kev / 1.0).powi(5);
    let n_particles = dm_mass_solar * 1.989e30 / (mass_kev * 1e3 * 1.782_661_92e-36);
    let d = distance_mpc * 3.0857e22;
    gamma * n_particles / (4.0 * PI * d * d)
}

pub fn dodelson_widrow_mass_kev(omega_dm: f64) -> f64 {
    1.8 * (omega_dm / 0.12).powf(1.0 / 1.8)
}

pub fn resonant_production_mass_kev(lepton_asymmetry: f64, omega_dm: f64) -> f64 {
    0.3 * (omega_dm / 0.12).powf(0.5) * (1e-3 / lepton_asymmetry).powf(0.5)
}