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
156
157
158
159
160
161
162
163
164
165
166
167

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

pub const ALPHA_EM: f64 = 1.0 / 137.035_999;

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct DarkPhoton {
    pub mass_gev: f64,
    pub kinetic_mixing: f64,
}

impl DarkPhoton {
    pub fn mass_kg(&self) -> f64 {
        self.mass_gev * 1.782_661_92e-27
    }

    pub fn decay_width_to_electrons_gev(&self) -> f64 {
        let me_gev = 0.000_511;
        if self.mass_gev < 2.0 * me_gev {
            return 0.0;
        }
        let beta = (1.0 - 4.0 * me_gev * me_gev / (self.mass_gev * self.mass_gev)).sqrt();
        self.kinetic_mixing * self.kinetic_mixing * ALPHA_EM * self.mass_gev / 3.0
            * (1.0 + 2.0 * me_gev * me_gev / (self.mass_gev * self.mass_gev))
            * beta
    }

    pub fn lifetime_seconds(&self) -> f64 {
        let width = self.decay_width_to_electrons_gev();
        if width <= 0.0 {
            return f64::INFINITY;
        }
        HBAR / (width * 1.782_661_92e-27 * C * C)
    }

    pub fn decay_length_meters(&self, gamma_lorentz: f64) -> f64 {
        let tau = self.lifetime_seconds();
        if tau.is_infinite() {
            return f64::INFINITY;
        }
        C * tau * gamma_lorentz
    }

    pub fn production_cross_section_cm2(&self, s_gev2: f64) -> f64 {
        let eps2 = self.kinetic_mixing * self.kinetic_mixing;
        let m2 = self.mass_gev * self.mass_gev;
        4.0 * PI * ALPHA_EM * ALPHA_EM * eps2 * s_gev2
            / ((s_gev2 - m2).powi(2)
                + (self.decay_width_to_electrons_gev() * self.mass_gev).powi(2))
            * 0.389_379e-27
    }
}

pub fn kinetic_mixing_from_loop(coupling: f64, mass_heavy_gev: f64) -> f64 {
    coupling * ALPHA_EM / (6.0 * PI) * (mass_heavy_gev / 1.0).ln()
}

pub fn dark_fine_structure_constant(g_dark: f64) -> f64 {
    g_dark * g_dark / (4.0 * PI)
}

pub fn portal_coupling_higgs(lambda_hp: f64, v_higgs_gev: f64, v_dark_gev: f64) -> f64 {
    lambda_hp * v_higgs_gev * v_dark_gev
}

pub fn invisible_branching_ratio(width_invisible_gev: f64, width_visible_gev: f64) -> f64 {
    width_invisible_gev / (width_invisible_gev + width_visible_gev)
}

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct ChameleonField {
    pub coupling_matter: f64,
    pub lambda_scale_ev: f64,
    pub n_power: f64,
}

impl ChameleonField {
    pub fn effective_mass_ev(&self, rho_kg_m3: f64) -> f64 {
        let rho_ev4 = rho_kg_m3 / 1.782_661_92e-36 * 1e9;
        let lambda4 = self.lambda_scale_ev.powi(4);
        let n = self.n_power;
        ((n + 1.0)
            * lambda4
            * (self.coupling_matter * rho_ev4 / lambda4).powf((n + 2.0) / (n + 1.0)))
        .powf(0.5 / (n + 2.0))
    }

    pub fn compton_wavelength_meters(&self, rho_kg_m3: f64) -> f64 {
        let m_ev = self.effective_mass_ev(rho_kg_m3);
        if m_ev <= 0.0 {
            return f64::INFINITY;
        }
        HBAR / (m_ev * 1.782_661_92e-36 * C)
    }

    pub fn thin_shell_factor(&self, rho_inside: f64, rho_outside: f64, radius: f64) -> f64 {
        let phi_inside = self.field_value_ev(rho_inside);
        let phi_outside = self.field_value_ev(rho_outside);
        let delta_phi = (phi_outside - phi_inside).abs();
        let m_pl_ev = 1.221e28;
        delta_phi / (self.coupling_matter * m_pl_ev * radius)
    }

    fn field_value_ev(&self, rho_kg_m3: f64) -> f64 {
        let rho_ev4 = rho_kg_m3 / 1.782_661_92e-36 * 1e9;
        let lambda4 = self.lambda_scale_ev.powi(4);
        let n = self.n_power;
        lambda4.powf(1.0 / (n + 1.0)) / (self.coupling_matter * rho_ev4).powf(1.0 / (n + 1.0))
    }
}

#[derive(Clone, Copy, Debug, PartialEq)]
pub struct SymmetronField {
    pub mu_gev: f64,
    pub lambda_coupling: f64,
    pub m_scale_gev: f64,
}

impl SymmetronField {
    pub fn critical_density_kg_m3(&self) -> f64 {
        self.mu_gev * self.mu_gev * self.m_scale_gev * self.m_scale_gev * 1.782_661_92e-27
            / (HBAR.powi(3) * C.powi(3))
            * 1e-18
    }

    pub fn symmetry_restored(&self, rho_kg_m3: f64) -> bool {
        rho_kg_m3 > self.critical_density_kg_m3()
    }

    pub fn vev_gev(&self, rho_kg_m3: f64) -> f64 {
        let rho_crit = self.critical_density_kg_m3();
        if rho_kg_m3 >= rho_crit {
            return 0.0;
        }
        (self.mu_gev * self.mu_gev / self.lambda_coupling * (1.0 - rho_kg_m3 / rho_crit)).sqrt()
    }

    pub fn effective_mass_gev(&self, rho_kg_m3: f64) -> f64 {
        let rho_crit = self.critical_density_kg_m3();
        if rho_kg_m3 >= rho_crit {
            return (self.mu_gev * self.mu_gev * (rho_kg_m3 / rho_crit - 1.0)).sqrt();
        }
        (2.0 * self.mu_gev * self.mu_gev * (1.0 - rho_kg_m3 / rho_crit)).sqrt()
    }
}

pub fn yukawa_potential(g: f64, r: f64, m_mediator_kg: f64) -> f64 {
    let inv_range = m_mediator_kg * C / HBAR;
    -g * g / (4.0 * PI * r) * (-inv_range * r).exp()
}

pub fn mediator_range_meters(mass_gev: f64) -> f64 {
    HBAR / (mass_gev * 1.782_661_92e-27 * C)
}

pub fn sommerfeld_enhancement(alpha: f64, v_over_c: f64) -> f64 {
    let x = PI * alpha / v_over_c;
    x / (1.0 - (-2.0 * x).exp())
}

pub fn dark_sector_temperature(xi: f64, t_photon_gev: f64) -> f64 {
    xi * t_photon_gev
}

pub fn delta_neff_from_dark_radiation(g_star_dark: f64, xi: f64) -> f64 {
    g_star_dark * (xi).powi(4) * 4.0 / 7.0
}