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

use crate::constants::cosmic::{
    HUBBLE_CONSTANT_SI, OMEGA_BARYON, OMEGA_DM, OMEGA_LAMBDA, OMEGA_MATTER,
};
use crate::constants::physical::{C, G, HBAR, K_B};
use std::f64::consts::PI;

pub const T_CMB_0: f64 = 2.7255;
pub const T_21CM_K: f64 = 0.068;
pub const A_10_S_INV: f64 = 2.85e-15;
pub const FREQ_21CM_MHZ: f64 = 1420.405;

pub fn brightness_temperature_21cm(x_hi: f64, delta_b: f64, t_spin_k: f64, z: f64) -> f64 {
    let t_cmb = T_CMB_0 * (1.0 + z);
    let factor = 27.0
        * x_hi
        * (1.0 + delta_b)
        * (OMEGA_BARYON * 0.7 * 0.7 / 0.023).sqrt()
        * ((1.0 + z) / 10.0 * 0.15 / (OMEGA_MATTER * 0.7 * 0.7)).sqrt();
    factor * (1.0 - t_cmb / t_spin_k)
}

pub fn spin_temperature(t_cmb_k: f64, t_gas_k: f64, t_lya_k: f64, x_c: f64, x_alpha: f64) -> f64 {
    let inv_t_s = (1.0 / t_cmb_k + x_c / t_gas_k + x_alpha / t_lya_k) / (1.0 + x_c + x_alpha);
    1.0 / inv_t_s
}

pub fn collisional_coupling(n_h_cm3: f64, t_gas_k: f64, z: f64) -> f64 {
    let kappa_hh = 3.1e-11 * t_gas_k.powf(0.357) * (-32.0 / t_gas_k).exp();
    let t_star = T_21CM_K;
    4.0 * t_star * n_h_cm3 * kappa_hh / (3.0 * A_10_S_INV * T_CMB_0 * (1.0 + z))
}

pub fn wouthuysen_field_coupling(j_alpha_cm2_s_hz_sr: f64, t_gas_k: f64) -> f64 {
    let s_alpha = 1.0 - 0.0631789 / t_gas_k + 0.115995 / (t_gas_k * t_gas_k);
    1.81e11 * s_alpha * j_alpha_cm2_s_hz_sr / (1.0 + 2.0 / t_gas_k)
}

pub fn gas_temperature_adiabatic(z: f64, z_dec: f64) -> f64 {
    T_CMB_0 * (1.0 + z_dec) * ((1.0 + z) / (1.0 + z_dec)).powi(2)
}

pub fn dm_heating_gas_temperature(t_gas_k: f64, epsilon_heat_ev_s_baryon: f64, dt_s: f64) -> f64 {
    let delta_t = epsilon_heat_ev_s_baryon * 1.602e-19 * dt_s / (1.5 * K_B);
    t_gas_k + delta_t
}

pub fn dm_baryon_scattering_heating_rate(
    sigma_0_cm2: f64,
    m_dm_gev: f64,
    m_b_gev: f64,
    rho_dm_gev_cm3: f64,
    t_dm_k: f64,
    t_gas_k: f64,
    v_rel_km_s: f64,
) -> f64 {
    let sigma_si = sigma_0_cm2 * 1e-4;
    let m_dm = m_dm_gev * 1.782_661_92e-27;
    let m_b = m_b_gev * 1.782_661_92e-27;
    let mu = m_dm * m_b / (m_dm + m_b);
    let n_dm = rho_dm_gev_cm3 * 1.782_661_92e-27 * 1e6 / m_dm;
    let v_th = ((K_B * t_gas_k / m_b + K_B * t_dm_k / m_dm).max(0.0)).sqrt();
    let v_rel = v_rel_km_s * 1e3;
    let v_eff = (v_rel * v_rel + v_th * v_th).sqrt();
    2.0 * mu / (m_dm + m_b) * n_dm * sigma_si * v_eff * K_B * (t_dm_k - t_gas_k)
}

pub fn power_spectrum_21cm(k_mpc: f64, delta_tb: f64, bias_hi: f64, p_matter: f64) -> f64 {
    let beam_suppression = (-0.5 * (k_mpc / 100.0).powi(2)).exp();
    delta_tb * delta_tb * bias_hi * bias_hi * p_matter * beam_suppression
}

pub fn optical_depth_21cm(n_hi_cm2: f64, t_spin_k: f64) -> f64 {
    5.9e-14 * n_hi_cm2 / t_spin_k
}

pub fn observed_frequency_mhz(z: f64) -> f64 {
    FREQ_21CM_MHZ / (1.0 + z)
}

pub fn dark_ages_redshift_range() -> (f64, f64) {
    (30.0, 200.0)
}

pub fn cosmic_dawn_redshift_range() -> (f64, f64) {
    (6.0, 30.0)
}

pub fn global_signal_absorption_trough_mhz(z_trough: f64) -> f64 {
    FREQ_21CM_MHZ / (1.0 + z_trough)
}

pub fn dm_21cm_excess_absorption(t_gas_standard_k: f64, t_gas_with_dm_k: f64, t_cmb_k: f64) -> f64 {
    let tb_standard = 27.0 * (1.0 - t_cmb_k / t_gas_standard_k);
    let tb_dm = 27.0 * (1.0 - t_cmb_k / t_gas_with_dm_k);
    tb_dm - tb_standard
}

pub fn hubble_rate_21cm(z: f64) -> f64 {
    HUBBLE_CONSTANT_SI
        * (OMEGA_DM * (1.0 + z).powi(3) + OMEGA_BARYON * (1.0 + z).powi(3) + OMEGA_LAMBDA).sqrt()
}

pub fn comoving_distance_21cm_mpc(z: f64, n_steps: usize) -> f64 {
    let dz = z / n_steps as f64;
    let mut dist = 0.0;
    for i in 0..n_steps {
        let zp = (i as f64 + 0.5) * dz;
        let h_z = hubble_rate_21cm(zp);
        dist += C / h_z * dz;
    }
    dist / 3.0857e22
}

pub fn dm_mass_in_21cm_volume_kg(z_min: f64, z_max: f64) -> f64 {
    let rho_crit = 3.0 * HUBBLE_CONSTANT_SI.powi(2) / (8.0 * PI * G);
    let rho_dm = rho_crit * OMEGA_DM;
    let z_mid = 0.5 * (z_min + z_max);
    let d_com = C * (z_max - z_min) / hubble_rate_21cm(z_mid);
    rho_dm * (1.0 + z_mid).powi(3) * d_com.powi(3)
}

pub fn de_broglie_smoothing_scale_mpc(m_dm_gev: f64, v_km_s: f64) -> f64 {
    let m_kg = m_dm_gev * 1.782_661_92e-27;
    let v = v_km_s * 1e3;
    HBAR / (m_kg * v) / 3.0857e22
}

pub fn matter_density_at_z(z: f64) -> f64 {
    let rho_crit = 3.0 * HUBBLE_CONSTANT_SI.powi(2) / (8.0 * PI * G);
    rho_crit * OMEGA_MATTER * (1.0 + z).powi(3)
}