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

use crate::constants::physical::{C, HBAR};

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Wimp {
    pub mass_gev: f64,
    pub sigma_si_cm2: f64,
    pub sigma_sd_cm2: f64,
    pub sigma_v_cm3_s: f64,
}

impl Wimp {
    pub fn mass_kg(&self) -> f64 {
        self.mass_gev * 1.783e-27
    }

    pub fn relic_density(&self) -> f64 {
        let thermal = 3e-26;
        0.12 * thermal / self.sigma_v_cm3_s
    }

    pub fn is_thermal_relic(&self) -> bool {
        let ratio = self.sigma_v_cm3_s / 3e-26;
        ratio > 0.1 && ratio < 10.0
    }

    pub fn freeze_out_temperature_gev(&self) -> f64 {
        self.mass_gev / freeze_out_xf(self.mass_gev, self.sigma_v_cm3_s)
    }

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

    pub fn velocity_at_temperature_gev(&self, t_gev: f64) -> f64 {
        (3.0 * t_gev / self.mass_gev).sqrt() * C
    }

    pub fn de_broglie_wavelength(&self, v: f64) -> f64 {
        HBAR / (self.mass_kg() * v)
    }

    pub fn number_density_local(&self, rho_gev_cm3: f64) -> f64 {
        rho_gev_cm3 / self.mass_gev
    }

    pub fn scattering_rate(&self, rho_gev_cm3: f64, v_km_s: f64, n_target_per_kg: f64) -> f64 {
        let n_dm = rho_gev_cm3 * 1e6 / self.mass_gev;
        let v = v_km_s * 1e5;
        n_dm * v * self.sigma_si_cm2 * n_target_per_kg
    }
}

pub fn freeze_out_xf(mass_gev: f64, sigma_v: f64) -> f64 {
    Wimp::freeze_out_xf(mass_gev, sigma_v)
}

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

pub const UNITARITY_BOUND_MASS_GEV: f64 = 340e3;

pub const LEE_WEINBERG_BOUND_GEV: f64 = 2.0;

pub fn si_cross_section_neutrino_floor(mass_gev: f64) -> f64 {
    if mass_gev < 6.0 {
        1e-45
    } else if mass_gev < 100.0 {
        1e-49
    } else {
        1e-48
    }
}

pub fn neutralino_mass_estimate(m1: f64, m2: f64, mu: f64, tan_beta: f64) -> f64 {
    let sin_2beta = 2.0 * tan_beta / (1.0 + tan_beta * tan_beta);
    let mz = 91.2;
    let sw2 = 0.231;
    m1.min(m2).min(mu) - mz * mz * (1.0 - sw2) * sin_2beta / (2.0 * mu)
}