hnl: Heavy Neutral Lepton Lifetime and Decay Widths

The hnl module computes the phenomenology of Heavy Neutral Leptons (HNLs) within the nuMSM (Neutrino Minimal Standard Model) framework. Given an HNL mass and a set of squared mixing angles [U²_e, U²_μ, U²_τ], the module calculates partial decay widths for every accessible final state — covering three-body leptonic decays, neutral and charged meson plus neutrino final states, and quark-level processes — then sums them to determine the total decay width, lifetime, and individual branching ratios. The implementation was written for the SHiP experiment sensitivity studies and includes full Majorana factors (factor-of-two enhancement from charge-conjugate channels) and three-loop QCD corrections for hadronic decay widths above 1 GeV.

from hnl import HNL

# Initialise with mass 1 GeV and mixing angles
b = HNL(1.0, [1e-8, 2e-8, 1e-9], debug=True)
# HNLbranchings instance initialized with couplings:
#      U2e   = 1e-08
#      U2mu  = 2e-08
#      U2tau = 1e-09
# and mass:
#      m = 1.0 GeV

lifetime = b.computeNLifetime()
print(f"Lifetime: {lifetime:.3e} s")              # e.g. 4.778e-05 s

br = b.findBranchingRatio("N -> pi+ mu-")
print(f"BR(N -> pi+ mu-): {br:.4f}")              # e.g. 0.1183

# Inspect all kinematically allowed channels
allowed = b.allowedChannels()
for ch, status in allowed.items():
    if status == "yes":
        print(ch, b.findBranchingRatio(ch))

Module-level helpers

`PDGname`

PDGname(particle) -> str

Translate a decay-product name used internally within the hnl module into the name recognised by ROOT’s TDatabasePDG. Specifically, the prime notation 1 (e.g. "eta1") is converted to "'" (e.g. "eta'"), and the generic neutrino name "nu" is resolved to "nu_e" to handle the inclusive three-neutrino channel.

particle

str

required

Internal particle name as used in HNLbranchings.decays channel strings, such as "eta1", "nu", "pi+", "mu-".

return

str

PDG database name for the particle, ready to be passed to ROOT.TDatabasePDG.Instance().GetParticle(name).

`mass`

mass(particle: str | None) -> float

Look up the mass of a particle from ROOT’s TDatabasePDG database. The name is first normalised through PDGname before the database query.

particle

str | None

required

Internal particle name string (see PDGname).

return

float

Particle mass in GeV as returned by TParticlePDG.Mass().

`lifetime`

lifetime(particle) -> float

Look up the lifetime of a particle from ROOT’s TDatabasePDG database.

particle

str

required

Internal particle name string (see PDGname).

return

float

Particle lifetime in seconds as returned by TParticlePDG.Lifetime().

`CKMmatrix`

Stores the CKM matrix elements from the PDG 2017 review, used when computing charged-current decay widths.

from hnl import CKMmatrix
ckm = CKMmatrix()
print(ckm.Vud)   # 0.9743
print(ckm.Vus)   # 0.2251

Attribute	Value
`Vud`	0.9743
`Vus`	0.2251
`Vub`	3.6 × 10⁻³
`Vcd`	0.2249
`Vcs`	0.9735
`Vcb`	4.11 × 10⁻²
`Vtd`	8.6 × 10⁻³
`Vts`	4.0 × 10⁻²
`Vtb`	0.999

`constants`

Stores physical constants used in HNL decay calculations.

Attribute	Description	Value
`GF`	Fermi constant	1.166379 × 10⁻⁵ GeV⁻²
`s2thetaw`	sin²θ_W (Weinberg angle)	0.23126
`heV`	ℏ in eV·s	6.58211928 × 10⁻¹⁶
`hGeV`	ℏ in GeV·s	6.58211928 × 10⁻²⁵
`decayConstant`	Meson decay constants f_h	dict keyed by meson name (GeV)

Meson decay constants (values in GeV):

Meson	f_h (GeV)
`pi+`, `pi0`	0.130
`eta`	1.2 × 0.130 = 0.156
`eta1` (η′)	0.152
`eta_c`	0.335
`rho+`, `rho0`	0.209
`omega`	0.195
`phi`	0.229
`D_s+`	0.249
`D*_s+`	0.315

`HNL`

HNL is the main public class. It inherits all decay width methods from HNLbranchings and adds the computeNLifetime convenience method.

Constructor

HNL(mass, couplings, debug=False)

mass

float

required

HNL mass in GeV.

couplings

list[float]

required

A list of three squared mixing angles [U²_e, U²_μ, U²_τ]. These are the squared moduli of the lepton-flavour mixing matrix elements.

debug

bool

default:"False"

When True, prints a summary of the configured mass and coupling values to stdout on initialisation.

`computeNLifetime`

computeNLifetime(system: str = "SI") -> float

Compute the HNL lifetime as ℏ / Γ_total, where Γ_total is the sum of all partial decay widths.

system

str

default:"\"SI\""

Unit system for the result. "SI" returns the lifetime in seconds. "FairShip" returns the lifetime in nanoseconds (multiplied by 10⁹).

return

float

HNL lifetime in seconds (SI) or nanoseconds (FairShip).

b = HNL(1.0, [1e-8, 2e-8, 1e-9])

tau_s  = b.computeNLifetime(system="SI")       # seconds
tau_ns = b.computeNLifetime(system="FairShip") # nanoseconds
print(f"τ = {tau_s:.3e} s = {tau_ns:.3e} ns")

`findBranchingRatio`

findBranchingRatio(decayString) -> float

Return the branching ratio for a specific decay channel. The method divides the partial width of the requested channel by the total HNL decay width. The special aggregate strings "N -> hadrons" and "N -> charged hadrons" are also accepted.

decayString

str

required

Decay channel identifier. Must be one of the strings listed in the Supported decay channels table below, or one of the aggregate strings "N -> hadrons" / "N -> charged hadrons". The process exits if an unrecognised string is supplied.

return

float

Branching ratio in [0, 1]. Returns 0 for kinematically forbidden channels.

`allowedChannels`

allowedChannels() -> dict[str, str]

Return a dictionary mapping every defined channel string to "yes" (kinematically allowed) or "no" (forbidden given the current HNL mass).

return

dict[str, str]

Dictionary with channel strings as keys and "yes" or "no" as values.

`HNLbranchings`

HNLbranchings is the base class that implements all partial width calculations. HNL inherits from it. You can use these methods directly if you need individual partial widths rather than branching ratios.

Internal width methods

Method	Description
`Width_3nu()`	Total width into three neutrinos (inclusive over all flavours)
`Width_nu_f_fbar(alpha, beta)`	Width into neutrino + fermion–antifermion pair via Z (and W for same-flavour)
`Width_l1_l2_nu2(alpha, beta)`	Width into two different-flavour charged leptons + neutrino via W
`Width_l_u_d(alpha, beta, gamma)`	Width into charged lepton + up-quark + down-quark via W
`Width_H0_nu(H, alpha)`	Width into neutral meson + neutrino
`Width_H_l(H, alpha)`	Width into charged meson + charged lepton
`Width_charged_leptons()`	Total leptonic width (charged leptons in final state)
`Width_neutral_mesons()`	Total width into neutral meson + neutrino
`Width_charged_mesons()`	Total width into charged meson + charged lepton
`Width_quarks_neutrino()`	Quark-level width with neutrino; zero below 1 GeV
`Width_quarks_lepton()`	Quark-level width with charged lepton; zero below 1 GeV
`NDecayWidth()`	Total HNL decay width
`QCD_correction()`	Three-loop QCD correction factor (only for `MN ≥ 1 GeV`)
`sqrt_lambda(a, b, c)`	Källén triangle function √λ(a,b,c); returns 0 in forbidden region
`integral(x1, x2, x3)`	Numerical Gaussian integral for W-mediated three-body phase space

Flavour indices throughout: 1 = e, 2 = μ, 3 = τ. For quarks in Width_nu_f_fbar: 4=u, 5=d, 6=s, 7=c, 8=b, 9=t.

Supported decay channels

The complete list of channel strings accepted by findBranchingRatio and returned by allowedChannels:

Channel string	Process type	Threshold (approx.)
`N -> nu nu nu`	3ν inclusive	0
`N -> e- e+ nu_e`	Z-mediated, same flavour	2 m_e
`N -> e- e+ nu_mu`	Z-mediated	2 m_e
`N -> e- e+ nu_tau`	Z-mediated	2 m_e
`N -> e- mu+ nu_mu`	W-mediated	m_e + m_μ
`N -> mu- e+ nu_e`	W-mediated	m_e + m_μ
`N -> pi0 nu_e`	Neutral meson	m_π⁰ ≈ 135 MeV
`N -> pi0 nu_mu`	Neutral meson	m_π⁰
`N -> pi0 nu_tau`	Neutral meson	m_π⁰
`N -> pi+ e-`	Charged meson	m_π⁺ + m_e
`N -> mu- mu+ nu_e`	Z-mediated	2 m_μ
`N -> mu- mu+ nu_mu`	Z-mediated, same flavour	2 m_μ
`N -> mu- mu+ nu_tau`	Z-mediated	2 m_μ
`N -> pi+ mu-`	Charged meson	m_π⁺ + m_μ
`N -> eta nu_e`	Neutral meson	m_η ≈ 548 MeV
`N -> eta nu_mu`	Neutral meson	m_η
`N -> eta nu_tau`	Neutral meson	m_η
`N -> rho0 nu_e`	Neutral vector meson	m_ρ⁰ ≈ 775 MeV
`N -> rho0 nu_mu`	Neutral vector meson	m_ρ⁰
`N -> rho0 nu_tau`	Neutral vector meson	m_ρ⁰
`N -> rho+ e-`	Charged vector meson	m_ρ⁺ + m_e
`N -> omega nu_e`	Neutral vector meson	m_ω ≈ 783 MeV
`N -> omega nu_mu`	Neutral vector meson	m_ω
`N -> omega nu_tau`	Neutral vector meson	m_ω
`N -> rho+ mu-`	Charged vector meson	m_ρ⁺ + m_μ
`N -> eta1 nu_e`	Neutral meson (η′)	m_η′ ≈ 958 MeV
`N -> eta1 nu_mu`	Neutral meson (η′)	m_η′
`N -> eta1 nu_tau`	Neutral meson (η′)	m_η′
`N -> phi nu_e`	Neutral vector meson	m_φ ≈ 1020 MeV
`N -> phi nu_mu`	Neutral vector meson	m_φ
`N -> phi nu_tau`	Neutral vector meson	m_φ
`N -> e- tau+ nu_tau`	W-mediated	m_e + m_τ
`N -> tau- e+ nu_e`	W-mediated	m_e + m_τ
`N -> mu- tau+ nu_tau`	W-mediated	m_μ + m_τ
`N -> tau- mu+ nu_mu`	W-mediated	m_μ + m_τ
`N -> D_s+ e-`	Charged pseudoscalar	m_Ds + m_e
`N -> D_s+ mu-`	Charged pseudoscalar	m_Ds + m_μ
`N -> D*_s+ e-`	Charged vector	m_D*s + m_e
`N -> D*_s+ mu-`	Charged vector	m_D*s + m_μ
`N -> eta_c nu_e`	Neutral charmonium	m_ηc ≈ 2984 MeV
`N -> eta_c nu_mu`	Neutral charmonium	m_ηc
`N -> eta_c nu_tau`	Neutral charmonium	m_ηc
`N -> hadrons`	Aggregate: all hadronic	—
`N -> charged hadrons`	Aggregate: charged hadrons only	—

The "eta1" channel corresponds to η′ (eta-prime). Internally, the prime symbol is stored as "1" in channel strings and converted by PDGname before PDG database queries.

Complete example

from hnl import HNL

# Create an HNL with mass 0.5 GeV, electron-dominated mixing
b = HNL(0.5, [1e-6, 0.0, 0.0], debug=True)

# Lifetime
tau = b.computeNLifetime()
print(f"Lifetime: {tau:.3e} s")

# Total decay width
gamma = b.NDecayWidth()
print(f"Total width: {gamma:.3e} GeV")

# Check which channels are open at 500 MeV
allowed = b.allowedChannels()
open_channels = [ch for ch, s in allowed.items() if s == "yes"]
print(f"Open channels at 500 MeV: {len(open_channels)}")

# Print all non-zero branching ratios
for ch in open_channels:
    br = b.findBranchingRatio(ch)
    if br > 0:
        print(f"  {ch:<30s}  BR = {br:.4f}")

# FairShip-unit lifetime for use in simulation
tau_ns = b.computeNLifetime(system="FairShip")
print(f"Lifetime in FairShip units: {tau_ns:.3e} ns")

Python Modules

Configuration

hnl: Heavy Neutral Lepton Lifetime and Decay Widths

Module-level helpers

`PDGname`

`mass`

`lifetime`

`CKMmatrix`

`constants`

`HNL`

Constructor

`computeNLifetime`

`findBranchingRatio`

`allowedChannels`

`HNLbranchings`

Internal width methods

Supported decay channels

Complete example

Build docs developers (and LLMs) love

Python Modules

Configuration

Documentation Index

​Module-level helpers

​PDGname

​mass

​lifetime

​CKMmatrix

​constants

​HNL

​Constructor

​computeNLifetime

​findBranchingRatio

​allowedChannels

​HNLbranchings

​Internal width methods

​Supported decay channels

​Complete example

Build docs developers (and LLMs) love

Module-level helpers

`PDGname`

`mass`

`lifetime`

`CKMmatrix`

`constants`

`HNL`

Constructor

`computeNLifetime`

`findBranchingRatio`

`allowedChannels`

`HNLbranchings`

Internal width methods

Supported decay channels

Complete example