Skip to content

add lepton asymmetry #39

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cgiovanetti merged 2 commits into from
Aug 28, 2026
Merged

add lepton asymmetry #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3435671
add lepton asymmetry
Aug 28, 2026
8db967e
fix docstring
Aug 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 10 additions & 1 deletion linx/abundances.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,7 @@ def __call__(
a_vec=None, t_vec=None,
eta_fac=jnp.asarray(1.), tau_n_fac = jnp.asarray(1.),
nuclear_rates_q=None, me = const.me,
xi_nu=jnp.asarray(0.),
Y_i=None, T_start=None, T_end=None, sampling_nTOp=150,
rtol=1e-6, atol=1e-9, solver=Kvaerno3(),
max_steps=4096,
Expand Down Expand Up @@ -158,6 +159,13 @@ def __call__(
distribution. If not specified, will be taken to be `q = 0`.
me : float, optional
Electron mass in MeV. Defaults to `const.me`.
xi_nu : float, optional
Electron-neutrino degeneracy parameter,
`xi_nu = mu_nu_e / T_nu_e`, i.e. the lepton asymmetry
chemical potential in units of the neutrino temperature.
Positive values favour n -> p and therefore lower Yp.
Defaults to 0. Note that this enters the weak rates only; it
does NOT backreact on `rho_nu_vec` / N_eff / H.
Y_i : tuple of float, optional
Initial abundances :math:`n_i/n_b` for species. Length must be equal to
`self.nuclear_net.max_i_species`. Must specify `T_start` and `T_end` if not `None`.
Expand Down Expand Up @@ -272,7 +280,8 @@ def __call__(

T_interval_nTOp, nTOp_frwrd, nTOp_bkwrd = self.weak_rates(
jnp.array([T_g_vec, T_nu_vec]),
T_start=T_start, T_end=T_end, sampling_nTOp=sampling_nTOp, me=me
T_start=T_start, T_end=T_end, sampling_nTOp=sampling_nTOp,
me=me, xi_nu=xi_nu
)

##################################
Expand Down
106 changes: 74 additions & 32 deletions linx/weak_rates.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,6 +134,7 @@ def __init__(self,
@eqx.filter_jit
def __call__(
self, T_vec_ref, T_start, T_end, sampling_nTOp, me=const.me,
xi_nu=jnp.asarray(0.),
):
"""
Evaluate n <-> p rates over range of EM temperatures.
Expand All @@ -150,6 +151,12 @@ def __call__(
Number of points between T_start and T_end to evaluate at.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter,
xi_nu = mu_nu_e / T_nu_e, i.e. the lepton asymmetry chemical
potential in units of the neutrino temperature. Positive
values favour n -> p and therefore lower Yp. Defaults to 0
(no lepton asymmetry).

Returns
-------
Expand All @@ -165,12 +172,14 @@ def __call__(
jnp.log10(T_start), jnp.log10(T_end), sampling_nTOp
)

nTOp_rates = self.nTOp_rates(T_interval, T_vec_ref, me)
nTOp_rates = self.nTOp_rates(T_interval, T_vec_ref, me, xi_nu)

return (T_interval, ) + nTOp_rates

@eqx.filter_vmap(in_axes=(None, 0, None, None))
def nTOp_rates(self, Tg, T_vec_ref, me=const.me):
@eqx.filter_vmap(in_axes=(None, 0, None, None, None))
def nTOp_rates(
self, Tg, T_vec_ref, me=const.me, xi_nu=jnp.asarray(0.)
):
"""
Dimensionless n <-> p rates, normalized to neutron decay width.

Expand All @@ -183,6 +192,9 @@ def nTOp_rates(self, Tg, T_vec_ref, me=const.me):
for computing the weak rates.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Returns
-------
Expand Down Expand Up @@ -233,8 +245,8 @@ def nTOp_rates(self, Tg, T_vec_ref, me=const.me):
pTOn_rate = 0.

y_CCR_vals = jnp.array([
self.dGamma_nTOp_dp(p_vals, x, xnu, me),
self.dGamma_pTOn_dp(p_vals, x, xnu, me)
self.dGamma_nTOp_dp(p_vals, x, xnu, me, xi_nu),
self.dGamma_pTOn_dp(p_vals, x, xnu, me, xi_nu)
])
CCR_rates = trapz(y_CCR_vals, p_vals)
nTOp_rate += CCR_rates[0] / lambda_0
Expand All @@ -243,15 +255,19 @@ def nTOp_rates(self, Tg, T_vec_ref, me=const.me):
if self.FM_corr:

y_FMCCR_vals = jnp.array([
self.ddelt_Gamma_nTOp_FM_dp(p_vals, x, xnu, me),
self.ddelt_Gamma_pTOn_FM_dp(p_vals, x, xnu, me)
self.ddelt_Gamma_nTOp_FM_dp(p_vals, x, xnu, me, xi_nu),
self.ddelt_Gamma_pTOn_FM_dp(p_vals, x, xnu, me, xi_nu)
])
FMCCR_rates = trapz(y_FMCCR_vals, p_vals)
nTOp_rate += FMCCR_rates[0] / lambda_0
pTOn_rate += FMCCR_rates[1] / lambda_0

if self.thermal_corr:

# NOTE: these tables are pre-tabulated assuming SBBN, i.e.
# xi_nu = 0, and so do not track a lepton asymmetry. The
# residual xi_nu dependence of these sub-percent
# corrections is neglected.
thermal_rates = jnp.array([
jnp.interp(
Tg,
Expand Down Expand Up @@ -543,7 +559,7 @@ def dlambda_dp_FM(self, pe, me=const.me):
* R_rad * self.Fermi(b, me)
)

def chi_Born(self, en, x, x_nu, sgnq, me=const.me):
def chi_Born(self, en, x, x_nu, sgnq, me=const.me, xi_nu=0.):
r"""
Integrand in momentum integral for Born weak rate.

Expand All @@ -559,18 +575,23 @@ def chi_Born(self, en, x, x_nu, sgnq, me=const.me):
Should have value +1 or -1, to switch between chi\_+ and chi\_-.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Notes
-----
See Pitrou+ 1801.08023 Eq. (79).

"""
Lepton asymmetry sign convention is mu = +xi_nu*T_nu for
nu_e (sgnq = +1, n -> p) and -xi_nu*T_nu for nubar_e
(sgnq = -1, p -> n).

# xi_nu = 0. # nu chemical potential set to zero for now.
"""

# sgnq = +1 corresponds to chi_plus.
e_nu = en - sgnq*(Q/me)
g_nu = expit(-x_nu*e_nu)
g_nu = expit(sgnq*xi_nu - x_nu*e_nu)
g_e = expit(-x*(-en))

return e_nu**2 * g_nu * g_e
Expand Down Expand Up @@ -607,7 +628,7 @@ def Fermi_sgn(self, sgnq, sgnE, b, me=const.me):
# n<->p Rates #
###############################

def dGamma_dp(self, p, x, x_nu, sgnq, me=const.me):
def dGamma_dp(self, p, x, x_nu, sgnq, me=const.me, xi_nu=0.):
"""
Integrand over momentum for n <-> p rate. including radiative
corrections.
Expand All @@ -627,6 +648,9 @@ def dGamma_dp(self, p, x, x_nu, sgnq, me=const.me):
+1 or -1, to select between n -> p or p -> n.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Returns
-------
Expand Down Expand Up @@ -656,16 +680,16 @@ def dGamma_dp(self, p, x, x_nu, sgnq, me=const.me):

return p**2 * (
(
self.chi_Born(en, x, x_nu, sgnq, me)
self.chi_Born(en, x, x_nu, sgnq, me, xi_nu)
* RC_term_plus
) + (
self.chi_Born(-en, x, x_nu, sgnq, me)
self.chi_Born(-en, x, x_nu, sgnq, me, xi_nu)
* RC_term_minus
)
)

@eqx.filter_vmap(in_axes=(None, 0, None, None, None))
def dGamma_nTOp_dp(self, p, x, xnu, me=const.me):
@eqx.filter_vmap(in_axes=(None, 0, None, None, None, None))
def dGamma_nTOp_dp(self, p, x, xnu, me=const.me, xi_nu=0.):
"""
Integrand over momentum for n -> p rate including radiative
corrections.
Expand All @@ -683,6 +707,9 @@ def dGamma_nTOp_dp(self, p, x, xnu, me=const.me):
electron mass.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Returns
-------
Expand All @@ -693,10 +720,10 @@ def dGamma_nTOp_dp(self, p, x, xnu, me=const.me):
See Pitrou+ 1801.08023 Eq. (101).
"""

return self.dGamma_dp(p, x, xnu, 1, me)
return self.dGamma_dp(p, x, xnu, 1, me, xi_nu)

@eqx.filter_vmap(in_axes=(None, 0, None, None, None))
def dGamma_pTOn_dp(self, p, x, xnu, me=const.me):
@eqx.filter_vmap(in_axes=(None, 0, None, None, None, None))
def dGamma_pTOn_dp(self, p, x, xnu, me=const.me, xi_nu=0.):
"""
Integrand over momentum for p -> n rate including radiative
corrections.
Expand All @@ -716,6 +743,9 @@ def dGamma_pTOn_dp(self, p, x, xnu, me=const.me):
+1 or -1, to select between n -> p or p -> n.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Returns
-------
Expand All @@ -726,15 +756,15 @@ def dGamma_pTOn_dp(self, p, x, xnu, me=const.me):
See Pitrou+ 1801.08023 Eq. (104).
"""

return self.dGamma_dp(p, x, xnu, -1, me)
return self.dGamma_dp(p, x, xnu, -1, me, xi_nu)


###########################
# Finite Mass Corrections #
###########################


def chi_FM(self, en, x, x_nu, sgnq, me=const.me):
def chi_FM(self, en, x, x_nu, sgnq, me=const.me, xi_nu=0.):
r"""
Integrand over momentum for finite mass correction to n <-> p rate.

Expand All @@ -753,6 +783,9 @@ def chi_FM(self, en, x, x_nu, sgnq, me=const.me):
+1 or -1 corresponding to chi\_+ or chi\_-.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Notes
-----
Expand Down Expand Up @@ -794,8 +827,8 @@ def chi_FM(self, en, x, x_nu, sgnq, me=const.me):
# Dimensionless energy of the neutrino.
en_nu = en - sgnq*Q / me

expit_neg = expit(-en_nu*x_nu)
expit_pos = expit(en_nu*x_nu)
expit_neg = expit(sgnq*xi_nu - en_nu*x_nu)
expit_pos = expit(en_nu*x_nu - sgnq*xi_nu)

res_e2p1 = (
2 * en_nu * expit_neg**2
Expand Down Expand Up @@ -855,7 +888,7 @@ def chi_FM(self, en, x, x_nu, sgnq, me=const.me):
)
return result

def ddelt_Gamma_FM_dp(self, p, x, znu, sgnq, me=const.me):
def ddelt_Gamma_FM_dp(self, p, x, znu, sgnq, me=const.me, xi_nu=0.):
"""
Integrand over momentum for finite mass corrections to the n <-> p
rate.
Expand All @@ -875,6 +908,9 @@ def ddelt_Gamma_FM_dp(self, p, x, znu, sgnq, me=const.me):
+1 or -1, to select between n -> p or p -> n.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Returns
-------
Expand Down Expand Up @@ -903,17 +939,17 @@ def ddelt_Gamma_FM_dp(self, p, x, znu, sgnq, me=const.me):

result = p**2 * (
(
self.chi_FM(en, x, znu, sgnq, me)
self.chi_FM(en, x, znu, sgnq, me, xi_nu)
* RC_term_plus
) + (
self.chi_FM(-en, x, znu, sgnq, me)
self.chi_FM(-en, x, znu, sgnq, me, xi_nu)
* RC_term_minus
)
)
return result

@eqx.filter_vmap(in_axes=(None, 0, None, None, None))
def ddelt_Gamma_nTOp_FM_dp(self, p, x, xnu, me=const.me):
@eqx.filter_vmap(in_axes=(None, 0, None, None, None, None))
def ddelt_Gamma_nTOp_FM_dp(self, p, x, xnu, me=const.me, xi_nu=0.):
"""
Integrand over momentum for finite mass corrections to the n -> p
rate.
Expand All @@ -931,6 +967,9 @@ def ddelt_Gamma_nTOp_FM_dp(self, p, x, xnu, me=const.me):
electron mass.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Returns
-------
Expand All @@ -941,10 +980,10 @@ def ddelt_Gamma_nTOp_FM_dp(self, p, x, xnu, me=const.me):
See Pitrou+ 1801.08023 Eq. (115).
"""

return self.ddelt_Gamma_FM_dp(p, x, xnu, 1, me)
return self.ddelt_Gamma_FM_dp(p, x, xnu, 1, me, xi_nu)

@eqx.filter_vmap(in_axes=(None, 0, None, None, None))
def ddelt_Gamma_pTOn_FM_dp(self, p, x, xnu, me=const.me):
@eqx.filter_vmap(in_axes=(None, 0, None, None, None, None))
def ddelt_Gamma_pTOn_FM_dp(self, p, x, xnu, me=const.me, xi_nu=0.):
"""
Integrand over momentum for finite mass corrections to the p -> n
rate.
Expand All @@ -962,6 +1001,9 @@ def ddelt_Gamma_pTOn_FM_dp(self, p, x, xnu, me=const.me):
electron mass.
me : float, optional
Electron mass in MeV. Defaults to const.me.
xi_nu : float, optional
Electron-neutrino degeneracy parameter, mu_nu_e / T_nu_e.
Defaults to 0 (no lepton asymmetry).

Returns
-------
Expand All @@ -972,4 +1014,4 @@ def ddelt_Gamma_pTOn_FM_dp(self, p, x, xnu, me=const.me):
See Pitrou+ 1801.08023 Eq. (115).
"""

return self.ddelt_Gamma_FM_dp(p, x, xnu, -1, me)
return self.ddelt_Gamma_FM_dp(p, x, xnu, -1, me, xi_nu)
Loading