#include <wdm_decay_product.h>

Inheritance diagram for WdmDecayProductSpecies:

Public Member Functions
double	SigmaAt (double u_cut) const

void	SetParent (const DCDMSpecies *dcdm)

void	FillInjection (double a, double rho_dcdm, double J, double dJdlnq) const

void	RegisterBackgroundIndices (int &index_bg) override

void	RegisterIntegrationIndices (int &index_bi) override

void	SetBackgroundInitialConditions (const BackgroundICContext &ctx) override

void	ComputeBackground (double a, const double pvecback_B, double pvecback) override

double	Rho (const double *pvecback) const override

double	P (const double *pvecback) const override

double	PPrime (double a, double H, const double pvecback_B, const double pvecback) const override

double	RhoDotOverRho (const double *pvecback, double a_prime_over_a) const override

void	WriteBackgroundColumnTitles (BackgroundColumnWriter &w) const override

void	WriteBackgroundData (const double *pvecback, BackgroundColumnWriter &w) const override

double	GetOmega0 () const override

double	NeutrinoOmega0 () const override

double	NeffContribution (double) const override

bool	IsFreestreaming () const override

double	BackgroundAIni (double a_proposed, double) const override

void	CheckUltraRelativisticAtIc (const double *, double) const override

bool	IsUltraRelativisticAtIc (const double *, double) const override

void	PrintNeffInfo () const override

void	PrintMassInfo () const override

void	PerturbDerivs (const BaseSpecies::PerturbLayout &layout, double tau, const double y, double dy, const perturb_parameters_and_workspace &ppaw) const override

void	ApplyInitialConditions (const BaseSpecies::PerturbLayout &layout, double *y, const PerturbIcContext &ctx) override

StressEnergyContribution	StressEnergy (const BaseSpecies::PerturbLayout &layout, const perturb_vector pv, const double y, const double pvecback, const perturb_workspace ppw) const override

void	FillSources (const BaseSpecies::PerturbLayout &layout, const double y, const double dy, PerturbSourceContext &ctx) const override

void	WriteOutputColumns (PerturbColumnWriter &writer, const PerturbationsModule &mod, file_format fmt, TransferColumnSection section=TransferColumnSection::all) const override

void	RegisterTensorPerturbationIndices (BaseSpecies::PerturbLayout &, perturb_vector , const precision , int &, const perturb_workspace *, int) override

void	PerturbSynchronousToNewtonian (const BaseSpecies::PerturbLayout &, double *, const PerturbIcContext &) override

Public Member Functions inherited from NCDMBaseSpecies
std::unique_ptr< BaseSpecies::PerturbLayout >	CreatePerturbLayout () const override

bool	ClustersAsMatter () const override

bool	IsColdMatterSpecies () const override

double	TensorMasslessRelativisticRho (const double *pvecback) const override

void	WarnIfTooHeavyForHalofit (double m_ev_threshold) const override

void	SetBackgroundModule (const BackgroundModule *bgm) override

void	SetPerturbs (const perturbs *ppt) override

void	RegisterTransferSourceIndices (int &index_tp, const SourceRequestContext &ctx) override

void	PerturbTensorDerivs (const BaseSpecies::PerturbLayout &layout, double tau, const double y, double dy, const perturb_parameters_and_workspace &ppaw) const override

void	ContributeTensorGwSource (const BaseSpecies::PerturbLayout &layout, double a, const double y, perturb_workspace ppw) const override

void	MarkUsedInSources (const BaseSpecies::PerturbLayout &layout, const perturb_workspace ppw, int used_in_sources) const override

void	CopyPerturbationsAcrossSwitch (const BaseSpecies::PerturbLayout &old_layout, const BaseSpecies::PerturbLayout &new_layout, const double old_y, double new_y, const PerturbSwitchContext &ctx) const override

Public Member Functions inherited from BaseSpecies
virtual void	SetThermodynamicsModule (const ThermodynamicsModule *)

virtual std::optional< double >	GetParam (const std::string &) const

virtual void	BackgroundDerivs (double tau, const double y, double dy, const double *pvecback)

virtual void	FinalizeBackground (double a, double H, const double pvecback_B, double pvecback)

virtual double	FreestreamingRho (const double *pvecback) const

bool	IsPresent () const

virtual void	ProcessBackgroundTable (const double , int, int, const double )

virtual void	PerturbVectorDerivs (const PerturbLayout &, double, const double , double , const perturb_parameters_and_workspace &) const

virtual void	WriteTensorOutputColumnTitles (std::string &) const

virtual void	PrintVariables (PerturbColumnWriter &, double, const double , const PerturbationsModule &, const perturb_workspace ) const

virtual double	GetRadiationOmega0 () const

virtual double	DarkRadiationRhoToday (const double *) const

bool	ClustersAsMatterCached () const

virtual void	FinalizeMatterClassification ()

virtual bool	HasWarmMatter () const

Static Public Attributes
static constexpr double	kSparsityFloor = 1e-12

Protected Member Functions
double	GetDlnf0Dlnq (int, const double *) const override

Protected Member Functions inherited from NCDMBaseSpecies
virtual double	GetW0ForGwSource (int iq, const double *pvecback) const

virtual double	EvaluatePsdAnalytic (double q) const

virtual quadrature_method	DefaultQuadratureStrategy () const

virtual void	FillQuadratureParams (GBQuadParams &) const

void	BuildQuadratureAndMass (const NcdmSettings &settings)

Protected Member Functions inherited from BaseSpecies
void	ApplyFluidLikeNewtonianShift (double y, int idx_delta, int idx_theta, const double pvecback, const PerturbIcContext &ctx) const

Additional Inherited Members
Public Types inherited from BaseSpecies
enum class	EnergyType

Detailed Description

Warm decay products of cold dark matter (arXiv:2606.14849).

A DCDM parent decays χ → χ₁ + χ₂ into two identical daughters, each with mass m_d = ε·m_χ/2 and kick velocity v = √(1−ε²)·c. All daughters are born with the same physical momentum, so the comoving distribution f(q,τ) builds up at the moving cutoff q_cut(τ) = a·kQKick. f and g ≡ ∂f/∂ln q are integrated per momentum bin in the background ODE. Injection is a conservative cell-integrated Gaussian (erf) energy deposit in ln q: the energy-injection sum rule Σᵢ dqᵢ qᵢ² εᵢ Jᵢ · factor/a⁴ = aΓρ_dcdm holds exactly for every cutoff position, and the source is smooth in time (both evolvers integrate it).

Internal conventions: T_wdm = T_cmb, deg = 1, dimensionless kick momentum kQKick = p/T₀ = 10 (only v_kick is physical; the parent mass drops out). Dimensionless daughter mass M = kQKick·ε/v.

Perturbations (Task 4) evolve UNNORMALIZED multipoles ψ_ℓ = (Δf)_ℓ, which stay regular through f = 0. Synchronous gauge only; no tensor slots; no fluid approximation.

Member Function Documentation

◆ SigmaAt()

double WdmDecayProductSpecies::SigmaAt ( double u_cut ) const

Injection kernel width at u_cut: the analytic local quantile-bin width Δu_loc = [(F_hi−F_lo)/N] · H_fid(a) e^{+Γ t_fid(a)} / Γ (a = e^{u_cut}/kQKick), clamped to [kSigmaMin, kSigmaMax]. Smooth in u_cut; evaluated in log form so no overflow/inf is ever materialized (-ffast-math discipline). Floor = the background-table resolvability bound; cap = the placement-bias budget.

◆ SetParent()

void WdmDecayProductSpecies::SetParent ( const DCDMSpecies * dcdm )

inline

Wired by DCDM_WDM_Species; used by RhoDotOverRho / FillSources.

◆ FillInjection()

void WdmDecayProductSpecies::FillInjection	(	double	a,
		double	rho_dcdm,
		double *	J,
		double *	dJdlnq
	)		const

Injection source Jᵢ = dfᵢ/dτ (and optionally dJᵢ/dln q) for the current (a, ρ_dcdm). Conservative cell-integrated Gaussian (erf) energy deposit in u = ln q at u_cut = ln(a·kQKick): the kick's energy A = aΓρ_dcdm·a⁴/factor_ is split by GaussWeights over the cells (σ = SigmaAt(u_cut)), off-grid tails clamped into the edge bins, so Σ Wᵢ = 1 identically and the exact sum rule Σ dqᵢ qᵢ² εᵢ Jᵢ = A holds for every u_cut. dJ/dln q is the full derivative of J(u) = A·φ_σ(u−u_cut)/(q³ε): the cell-averaged raw-kernel derivative (normal-pdf differences at the cell edges, smooth, no edge special cases) MINUS the measure's own log-derivative (3 + q²/ε²)·J, which supplies the g-channel's ρ-weighted moment (Σ dq q² ε dJdlnq ≈ −3A non-relativistically, by integration by parts) that drives the daughter's fluid-limit metric-driven growth. dJdlnq may be nullptr. Background-thread only (uses mutable scratch).

◆ RegisterBackgroundIndices()

void WdmDecayProductSpecies::RegisterBackgroundIndices ( int & index_bg )

overridevirtual

Claim consecutive slots in pvecback. Called once during background_indices(). Implementation must assign index_bg_rho_ (and any other needed indices) and increment index_bg accordingly.

Implements BaseSpecies.

◆ RegisterIntegrationIndices()

void WdmDecayProductSpecies::RegisterIntegrationIndices ( int & index_bi )

overridevirtual

Claim slots in the ODE integration vector y. Called during background_indices(). Default: no integrated variables (analytic species do not need this).

Reimplemented from BaseSpecies.

◆ SetBackgroundInitialConditions()

void WdmDecayProductSpecies::SetBackgroundInitialConditions ( const BackgroundICContext & ctx )

overridevirtual

Set initial conditions for ODE-integrated background variables.

Parameters

ctx	Background initial-condition context (a_ini, rho_rad, and the integration vector to fill).

Reimplemented from BaseSpecies.

◆ ComputeBackground()

void WdmDecayProductSpecies::ComputeBackground	(	double	a,
		const double *	pvecback_B,
		double *	pvecback
	)

overridevirtual

Compute species background quantities at scale factor a. pvecback_B is the current ODE integration vector. Write density (and any other owned quantities) into pvecback.

Implements BaseSpecies.

◆ Rho()

double WdmDecayProductSpecies::Rho ( const double * pvecback ) const

inlineoverridevirtual

Energy density at current background state.

Implements BaseSpecies.

◆ P()

double WdmDecayProductSpecies::P ( const double * pvecback ) const

inlineoverridevirtual

Pressure at current background state.

Implements BaseSpecies.

◆ PPrime()

double WdmDecayProductSpecies::PPrime	(	double	a,
		double	H,
		const double *	pvecback_B,
		const double *	pvecback
	)		const

overridevirtual

Post-Friedmann. This species' conformal-time pressure derivative p' (a' / a = a*H). pvecback_B is the ODE state (IN); pvecback is fully populated through H. Default 0 (matter / Lambda).

Reimplemented from BaseSpecies.

◆ RhoDotOverRho()

double WdmDecayProductSpecies::RhoDotOverRho	(	const double *	pvecback,
		double	a_prime_over_a
	)		const

overridevirtual

Log conformal-time derivative of the background density, ρ̇/ρ (ρ̇ ≡ dρ̄/dτ). Single source of truth for the density-contrast gauge shift; mirrors the continuity equation used in BackgroundDerivs. Default: ρ̇/ρ = -3ℋ(Rho+P)/Rho. Species with extra source/sink terms (e.g. decay) override this.

Parameters

pvecback	Current background-quantity vector.
a_prime_over_a	ℋ = a'/a (base does not store H/a indices).

Reimplemented from BaseSpecies.

◆ WriteBackgroundColumnTitles()

void WdmDecayProductSpecies::WriteBackgroundColumnTitles ( BackgroundColumnWriter & ) const

overridevirtual

Write this species' background column titles into writer. Called by BackgroundModule::background_output_titles(). Default: no-op (species with no background output need not override).

Reimplemented from BaseSpecies.

◆ WriteBackgroundData()

void WdmDecayProductSpecies::WriteBackgroundData	(	const double *	,
		BackgroundColumnWriter &
	)		const

overridevirtual

Write this species' background data values into writer. pvecback is the full background vector at the current time step. Called by BackgroundModule::background_output_data(). Default: no-op.

Reimplemented from BaseSpecies.

◆ GetOmega0()

double WdmDecayProductSpecies::GetOmega0 ( ) const

inlineoverridevirtual

Returns the species' contribution to Omega0 (used for budget closure during construction). For composites, this defaults to summing children. For decay-product species starting at zero, override returning 0.

Reimplemented from NCDMBaseSpecies.

◆ NeutrinoOmega0()

double WdmDecayProductSpecies::NeutrinoOmega0 ( ) const

inlineoverridevirtual

This species' contribution to Omega0 of the neutrino/NCDM sector (fnu). Default 0; NCDM returns GetOmega0(); composites sum children.

Reimplemented from NCDMBaseSpecies.

◆ NeffContribution()

double WdmDecayProductSpecies::NeffContribution ( double ) const

inlineoverridevirtual

Contribution to N_eff at redshift z. Default 0; NCDM returns GetNeff(z); composites sum children.

Reimplemented from NCDMBaseSpecies.

◆ IsFreestreaming()

bool WdmDecayProductSpecies::IsFreestreaming ( ) const

inlineoverridevirtual

True for species that are free-streaming and massless at IC time (deep radiation domination). Used to accumulate rho_nu / fracnu in perturbation IC setup. Static at construction time. Default: false.

Reimplemented from BaseSpecies.

◆ BackgroundAIni()

double WdmDecayProductSpecies::BackgroundAIni	(	double	a_proposed,
		double	tol
	)		const

inlineoverridevirtual

NCDM may need to start integration earlier so it is relativistic at a_ini.

Reimplemented from NCDMBaseSpecies.

◆ CheckUltraRelativisticAtIc()

void WdmDecayProductSpecies::CheckUltraRelativisticAtIc	(	const double *	,
		double
	)		const

inlineoverridevirtual

Initial-time ultra-relativistic check (throws via class_test on failure). Default no-op; NCDM checks |w-1/3|; composites forward to children.

Reimplemented from NCDMBaseSpecies.

◆ IsUltraRelativisticAtIc()

bool WdmDecayProductSpecies::IsUltraRelativisticAtIc	(	const double *	,
		double
	)		const

inlineoverridevirtual

Initial-time ultra-relativistic predicate. Default true; NCDM checks |w-1/3| <= tol; composites AND over children.

Reimplemented from NCDMBaseSpecies.

◆ PrintNeffInfo()

void WdmDecayProductSpecies::PrintNeffInfo ( ) const

overridevirtual

Verbose N_eff line (background_verbose). Default no-op; NCDM prints; composites forward to children.

Reimplemented from NCDMBaseSpecies.

◆ PrintMassInfo()

void WdmDecayProductSpecies::PrintMassInfo ( ) const

overridevirtual

Verbose mass/omega line (background_verbose). Default no-op; NCDM prints; composites forward to children.

Reimplemented from NCDMBaseSpecies.

◆ PerturbDerivs()

void WdmDecayProductSpecies::PerturbDerivs	(	const BaseSpecies::PerturbLayout &	layout,
		double	tau,
		const double *	y,
		double *	dy,
		const perturb_parameters_and_workspace &	ppaw
	)		const

overridevirtual

Contribute to dy for the scalar perturbation ODE at conformal time tau. The PerturbScalarContext inside ppaw->ppw has pre-computed metric terms, cross-species state (delta_g, theta_g, theta_b, ...), and approximation flags. The layout is this species' PerturbLayout slot (ppw->pv->species_layouts[i]); composites receive their own layout and pass nested sub-layouts to children.

Implements BaseSpecies.

◆ ApplyInitialConditions()

void WdmDecayProductSpecies::ApplyInitialConditions	(	const BaseSpecies::PerturbLayout &	,
		double *	,
		const PerturbIcContext &
	)

overridevirtual

Write synchronous-gauge initial conditions for this species into y[]. y[] == ppw->pv->y. The module pre-computes ctx.delta_g_ic, ctx.delta_ur, etc. Guard each IC type: if (ctx.index_ic == ctx.p_mod->index_ic_ad_) { ... } Newtonian gauge transformation is handled by the module after this loop. Default: no-op.

Reimplemented from BaseSpecies.

◆ StressEnergy()

BaseSpecies::StressEnergyContribution WdmDecayProductSpecies::StressEnergy	(	const BaseSpecies::PerturbLayout &	layout,
		const perturb_vector *	pv,
		const double *	y,
		const double *	pvecback,
		const perturb_workspace *	ppw
	)		const

overridevirtual

All scalar stress-energy perturbations in one call.

Implements BaseSpecies.

◆ FillSources()

void WdmDecayProductSpecies::FillSources	(	const BaseSpecies::PerturbLayout &	,
		const double *	,
		const double *	,
		PerturbSourceContext &
	)		const

overridevirtual

Write per-k, per-tau source values for this species into the source table. Called in perturb_sources_member() after metric/temperature sources are set. Use ctx.p_mod->SetSourceValue(ctx.index_md, ctx.index_ic, ctx.p_mod->index_tp_XXX_, ctx.index_tau, ctx.index_k, value). All addressing (index_md, index_ic, index_k, index_tau) is in ctx. Default: no-op.

Reimplemented from BaseSpecies.

◆ WriteOutputColumns()

void WdmDecayProductSpecies::WriteOutputColumns	(	PerturbColumnWriter &	,
		const PerturbationsModule &	,
		file_format	,
		BaseSpecies::TransferColumnSection	= `TransferColumnSection::all`
	)		const

overridevirtual

Append this species' output columns (delta, theta) to writer. Called from perturb_output_titles (title mode) and perturb_output_data (data mode). In class_format, the module may call this separately for density and velocity sections to preserve the historical transfer-column ordering. Use writer.Add(title, mod.index_tp_XXX_, active) — writer looks up tk[tp_index].

Reimplemented from BaseSpecies.

◆ RegisterTensorPerturbationIndices()

void WdmDecayProductSpecies::RegisterTensorPerturbationIndices	(	BaseSpecies::PerturbLayout &	,
		perturb_vector *	,
		const precision *	,
		int &	,
		const perturb_workspace *	,
		int
	)

inlineoverridevirtual

Tensor modes: no slots registered — the daughter's tensor anisotropic stress is neglected (documented scope limit; it is empty until late times).

Reimplemented from NCDMBaseSpecies.

◆ PerturbSynchronousToNewtonian()

void WdmDecayProductSpecies::PerturbSynchronousToNewtonian	(	const BaseSpecies::PerturbLayout &	,
		double *	,
		const PerturbIcContext &
	)

inlineoverridevirtual

Synchronous gauge only (guarded at CreateAll); base NCDM transform would zero the injected ICs via GetDlnf0Dlnq() == 0.

Reimplemented from NCDMBaseSpecies.

◆ GetDlnf0Dlnq()

double WdmDecayProductSpecies::GetDlnf0Dlnq	(	int	,
		const double *
	)		const

inlineoverrideprotectedvirtual

Only consumed by tensor derivs / base gauge transform, both disabled.

Implements NCDMBaseSpecies.

Field Documentation

◆ kSparsityFloor

constexpr double WdmDecayProductSpecies::kSparsityFloor = 1e-12

staticconstexpr

Static-sparsity seed for the injection source. ndf15's numjac derives the Jacobian sparsity pattern from exact zeros and locks it permanently once it repeats (tools/evolver_ndf15.cpp: pattern from fabs(dFdy)!=0, locked after trust_sparse repeats, never re-derived). A moving compact source would lock a too-narrow pattern early and corrupt the grouped Jacobian when the footprint sweeps on. Flooring every weight (renormalized, Σw stays exactly 1) keeps every injection row structurally coupled to rho_dcdm at all times at a physically invisible level (~1e-12 of the instantaneous deposit). This trick only works because J is LINEAR in rho_dcdm (the floor is a fixed fraction of the same A that scales with rho_dcdm) — a nonlinear coupling would decouple the seeded pattern from the finite-difference Jacobian numjac actually probes, defeating the pattern-lock.

The documentation for this class was generated from the following files:

wdm_decay_product.h
wdm_decay_product.cpp

Public Member Functions

Static Public Attributes

Protected Member Functions

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ SigmaAt()

◆ SetParent()

◆ FillInjection()

◆ RegisterBackgroundIndices()

◆ RegisterIntegrationIndices()

◆ SetBackgroundInitialConditions()

◆ ComputeBackground()

◆ Rho()

◆ P()

◆ PPrime()

◆ RhoDotOverRho()

◆ WriteBackgroundColumnTitles()

◆ WriteBackgroundData()

◆ GetOmega0()

◆ NeutrinoOmega0()

◆ NeffContribution()

◆ IsFreestreaming()

◆ BackgroundAIni()

◆ CheckUltraRelativisticAtIc()

◆ IsUltraRelativisticAtIc()

◆ PrintNeffInfo()

◆ PrintMassInfo()

◆ PerturbDerivs()

◆ ApplyInitialConditions()

◆ StressEnergy()

◆ FillSources()

◆ WriteOutputColumns()

◆ RegisterTensorPerturbationIndices()

◆ PerturbSynchronousToNewtonian()

◆ GetDlnf0Dlnq()

Field Documentation

◆ kSparsityFloor