Source-linked AI summary
micrOMEGAs5.0 : freeze-in
Geneviève Bélanger, Fawzi Boudjema, Andreas Goudelis, Alexander Pukhov, Bryan Zaldivar
TL;DR
Calculating freeze-in dark matter abundances in general extensions of the Standard Model is complicated by high-temperature loss of equilibrium. This paper develops formalism for different freeze-in scenarios in micrOMEGAs and finds that proper treatment of bath-particle statistics can shift predicted relic abundance by a factor of two.
Problem
Calculating freeze-in abundances in general Standard Model extensions is complicated by high-temperature loss of equilibrium.
Source: §6 · p. 28Method
The paper develops formalism for different freeze-in types and implements it in micrOMEGAs, accounting for bath-particle phase-space distributions beyond the Maxwell-Boltzmann approximation.
Source: §1 · p. 4, §6 · p. 29Results
Proper bath-particle statistics can change the predicted relic abundance by a factor of two, especially for annihilation processes initiated by bosons.
Source: §1 · p. 4Takeaways & Limitations
The Maxwell-Boltzmann approximation is not always adequate for predicting freeze-in abundances.
Source: §6 · p. 29The routine assumes all odd feebly interacting particles decay into the lightest one, and can overestimate the dark matter yield when such particles decay predominantly into bath particles.
Source: §3.2 · p. 19
Abstract
from arXiv · showhide
We present a major upgrade of the micrOMEGAs dark matter code to compute the abundance of feebly interacting dark matter candidates through the freeze-in mechanism in generic extensions of the Standard Model of particle physics. We develop the necessary formalism in order to solve the freeze-in Boltzmann equations while making as few simplifying assumptions as possible concerning the phase-space distributions of the particles involved in the dark matter production process. We further show that this formalism allows us to treat different freeze-in scenarios and discuss the way it is implemented in the code. We find that, depending on the New Physics scenario under consideration, the effect of a proper treatment of statistics on the predicted dark matter abundance can range from a few percent up to a factor of two, or more. We moreover illustrate the underlying physics, as well as the various novel functionalities of micrOMEGAs, by presenting several example results obtained for different dark matter models.
1 Introduction
The paper extends micrOMEGAs to calculate freeze-in abundances for FIMP dark matter and treats bath-particle statistics beyond the Maxwell–Boltzmann approximation. It presents the freeze-in scenarios and functionalities supported by the code.
1 Introduction
- Feebly coupled dark matter, with coupling strength typically O(10−10 −10−12), can evade standard WIMP-search bounds yet reproduce the measured density through scattering or decay without reaching thermal equilibrium [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33].Unlike freeze-out, the predicted relic density scales proportionally with the thermally averaged self-annihilation cross-section; production is often concentrated near the DM or mediator mass, but can be sensitive to high temperatures [23] [34] [35] [36].Source: §1 · p. 3
- Although FIMP Boltzmann equations had been solved in models, no public computational tool performed the calculation, motivating this extension of micrOMEGAs [38] [39] [40] [41] [43] [44].Source: §1 · p. 3
- Freeze-in production can proceed through mediator decays, SM annihilation into DM, or a feebly coupled DM–NLDSP pair; when the NLDSP or mediator is also feebly interacting, coupled Boltzmann equations may be needed.A WIMP NLDSP’s freeze-out abundance can instead determine the DM abundance.Source: §1 · p. 4
- micrOMEGAs adds FIMP freeze-in abundance calculations for simple or complex Standard Model extensions, including multiple mediators or dark-sector particles and FIMPs that later decay into WIMP dark matter.The code starts from a very small initial FIMP number density and tracks production by bath interactions or decays of mediator or dark-sector particles, whether or not they are in thermal equilibrium with the SM.Source: §1 · p. 4
- Up to a factor two, proper Fermi–Dirac and Bose–Einstein bath statistics can shift the predicted relic density, especially for boson-initiated annihilation.The paper illustrates micrOMEGAs results in a Z′-mediated Dirac-fermion model, the minimal singlet scalar model, and a model with an additional dark-sector particle.Source: §1 · pp. 4–5
- The code lets users specify FIMPs, reheating temperature, and included production processes, supports up to two DM candidates, and retains freeze-out routines; FIMP detection signatures are left for future work.Freeze-in is also available for the lightest or next-to-lightest dark-sector particle via freeze-out.Source: §1 · p. 5
2 Relic density of FIMPs
FIMPs can be produced from decays or annihilations of other particles without reaching thermal equilibrium, and their abundance freezes in when production subsides. Their annihilations into bath particles are negligible, and stronger visible-sector coupling implies a larger FIMP density.
2 Relic density of FIMPs
- FIMPs can be produced by decays or annihilations of other particles while remaining out of thermal equilibrium, typically until those particles are Boltzmann suppressed.This picture assumes an initially very small FIMP number density.Source: §2 · p. 5
- Negligible FIMP annihilations into bath particles allow the abundance to freeze in at a constant value that can match CMB measurements.Source: §2 · p. 5
- Freeze-in production can involve intermediate, potentially long-lived particles, and stronger FIMP coupling to the visible sector yields a larger dark matter density.Source: §2 · p. 5
2.1 Notations and useful formulas
This section establishes the phase-space and cosmological notation used to formulate freeze-in production. It summarizes particle distributions, thermodynamic densities, and the expansion and entropy relations connecting time and temperature.
2.1 Notations and useful formulas
- Particle number and energy densities are obtained from the phase-space distribution, with internal degrees of freedom, momentum, and energy specifying the relevant quantities.The density expressions are introduced for a species i with distribution f_i.Source: §2.1 · p. 6
- In kinetic equilibrium, particles follow Fermi–Dirac or Bose–Einstein distributions; the notation η_i incorporates the chemical potential, while η ≈ 0 for FIMPs.The text sets η_s = 1 for bosons and η_s = −1 for fermions.Source: §2.1 · p. 6
- For T ≳ 100 eV in the radiation-dominated Universe, energy and entropy densities are expressed using the effective degrees of freedom g_eff and h_eff.Source: §2.1 · p. 6
- The Hubble rate uses the Planck mass, with a rough estimate H(T) ≈ 10^-18(T/GeV)^2 GeV; entropy conservation relates time and temperature.Entropy conservation is stated as ds/dt = −3Hs.Source: §2.1 · p. 6
2.2 General freeze-in Boltzmann equation
The section formulates the general Boltzmann evolution and integrated collision term, then applies freeze-in assumptions to calculate the dark matter yield and distinguish production scenarios.
2.2 General freeze-in Boltzmann equation
- The general Boltzmann equation tracks a species’ number density through integrated reaction collision terms, built from phase-space distributions, squared matrix elements, and combinatorial factors.For 2 → 2 reactions, the factor is 1/2 for identical incoming particles and 1 otherwise.Source: §2.2 · p. 7
- Freeze-in assumes feeble visible-sector couplings and negligible initial dark matter abundance, so only dark matter production processes are retained.If annihilations can be ignored, a nonzero initial abundance instead adds to the freeze-in contribution; the analysis assumes negligible initial abundance throughout.Source: §2.2 · p. 7
- For exactly stable dark matter, the approximation is subject to a condition involving the thermally averaged annihilation cross section, which micrOMEGAs does not check.The text notes that the negligible-abundance assumption may not hold for general dark matter coupling strengths.Source: §2.2 · p. 8
- The final yield is obtained by integrating production from the reheating temperature TR to today; the abundance depends on TR when it is near the mediator or dark matter mass, or production is high-temperature dominated.The present micrOMEGAs version omits channels producing three or four dark matter particles, such as mediator-pair production followed by mediator decays.Source: §2.2 · p. 8
- The section distinguishes freeze-in through 1 → 2 decays of bath-equilibrated or non-equilibrated particles, including late decays, and 2 → 2 bath-particle annihilations; the lightest odd-sector particle is assumed to be a FIMP.Annihilation final states may be two dark matter particles or one dark matter particle and one bath particle.Source: §2.2 · p. 8
2.3 1 →2 (2 →1) processes
Section 2.3 develops collision terms for freeze-in production through decays and inverse annihilations, retaining the phase-space distributions and quantum statistics of the particles involved. It treats mediators in and out of chemical and kinetic equilibrium with the bath.
2.3.1 Mathematical detour
- The medium-modified decay rate incorporates final-state statistical factors (1 ∓ fa)(1 ∓ fb); in the dilute limit it recovers the usual in-vacuum decay rate.The derivation uses the Lorentz invariance of the phase-space element while retaining distribution functions in the comoving frame.Source: §2.3.1 · p. 9
2.3.2 Back to the collision term
- The integrated Y → a,b collision term includes the parent distribution and is expressed using a generalized ˜K1 function; its inverse a,b → Y term follows from the same expression.When all ηi = 0, ˜K1 reduces to the usual K1(mY/T).Source: §2.3.2 · pp. 10–11
2.3.3 The case of a mediator in thermal equilibrium with the SM bath
- The approximate equilibrium-mediator yield is accurate to about 1% for mY > 3 GeV, and its reheating-temperature dependence is below 1% when TR > 2mY.The approximation is given for a bosonic mediator.Source: §2.3.3 · pp. 11–12
- For an equilibrated mediator, Bose-Einstein statistics increase DM abundance by 3.5%, whereas Fermi-Dirac statistics decrease it by 2.5% relative to Maxwell-Boltzmann statistics.For a fermion mediator, a small mass splitting can make final-state statistics effects much larger.Source: §2.3.3 · p. 12
2.3.4 The case of a mediator not in chemical equilibrium with the SM bath
- When mediator couplings to bath particles are sufficiently suppressed, the mediator never reaches chemical equilibrium and decays at T ≪ mY; assuming kinetic equilibrium or kinetic decoupling changes DM abundance only at the percent level.The comparison finds kinetic equilibrium with the bath largely irrelevant for the resulting abundance in this regime.Source: §2.3.4 · p. 13
Mediator in kinetic equilibrium with the SM bath
- For a mediator kept in kinetic equilibrium with the bath, its out-of-equilibrium population requires solving for its yield and chemical-potential parameter self-consistently.The treatment defines a bath-modified average width and boost factor to rewrite the mediator evolution equation.Source: Mediator in kinetic equilibrium with… · pp. 13–14
Kinetically decoupled mediator
- For a kinetically decoupled mediator, the unintegrated Boltzmann equation tracks its phase-space distribution without energy-redistribution terms, and the numerical method solves coupled evolution equations on a momentum grid.The equations are solved using the method of characteristics, then interpolated to calculate the DM number density.Source: Kinetically decoupled mediator · pp. 15–16
- Fermi-Dirac statistics lower the final DM abundance by about 25% for fermion-coupled mediators, while boson-coupled mediators raise it by up to 80%; kinetic treatments agree at percent level.The fermion result holds whether the mediator decays rapidly into DM or its decay is delayed.Source: §2.3.4 · p. 13, Kinetically decoupled mediator · p. 16
2.4 2 →2 processes
The section formulates 2 →2 freeze-in production, distinguishing reactions with two FIMPs from those producing one FIMP and describing how s-channel resonances are treated to match mediator decays. For single-FIMP production, micrOMEGAs uses an empirical correction factor in place of a costly full integration.
2.4 2 →2 processes
- When both final particles are FIMPs, their statistical factors can be approximated as unity, yielding an expression similar to standard thermal freeze-out.The expression can be manipulated and rewritten using the ˜K1 function.Source: §2.4 · p. 16
- Producing one FIMP alongside a bath particle requires the full 12-dimensional collision-term integration, so micrOMEGAs instead inserts an empirical correction factor in the integrand.The bath-particle distribution is set to zero in this treatment.Source: §2.4 · p. 16
- 2% precision: the empirical correction factor reproduces the full numerical result for the section 4.3 model, assuming X thermalises rapidly.The comparison is between the corrected approximation and the full numerical calculation.Source: §2.4 · p. 16
2.4.1 The special case of an s-channel resonance
- For mediator-dominated production, micrOMEGAs temperature-dependently rescales the mediator’s total width so each 2 →2 process matches the mediator-decay result.The correction accounts for the mediator chemical potential and replaces its vacuum branching ratio into DM with the thermal-bath branching ratio; an overall cross-section factor would also alter off-resonance contributions. The described treatment covers DM pair production or final states containing only FIMPs, not χX final states with X in the bath.Source: §2.4.1 · p. 17
- For a mediator with a small decay width, the branching-ratio ansatz accounts for its survival from production temperature Tc until decay temperature Td.The survival probability has an integral representation; in practice, micrOMEGAs detects s-channel resonances, tabulates mediator functions assuming kinetic equilibrium with the SM bath, computes an effective width, and substitutes it into the matrix element.Source: §2.4.1 · pp. 17–18
3 Functions of micrOMEGAs
Section 3 describes micrOMEGAs helper functions and routines for freeze-in abundance calculations, along with user responsibilities and limitations that affect their use. It covers both the tools available and assumptions made by the code.
3.1 Constants and auxiliary functions
- Helpers expose radiation energy and entropy degrees of freedom, the radiation-era Hubble rate (valid for T ≳100 eV), the entropy derivative, and statistics-aware decay and 1→2 kernels.gEff(T) includes only SM particles; Stat2 returns the S function for mediator decay at fixed momentum, while K1to2 returns the statistics-aware ˜K1 function.Source: §3.1 · pp. 18–19
3.2 Freeze-in routines
- Users must declare feebly interacting particles; micrOMEGAs does not check thermal equilibrium, and a non-empty feeble-particle list affects freeze-out calculations.The freeze-out routines exclude listed odd particles, potentially affecting coannihilation; in specified bath/feeble-DM cases, the reported abundance is rescaled by the ratio of the lightest-particle masses. The code also does not check whether odd-particle decay rates are much smaller than H(TFO).Source: §3.2 · pp. 19, 21
- darkOmegaFiDecay computes abundance from a named particle’s decays, with a switch for whether that particle is in kinetic equilibrium with the SM.The KE=1 option is faster. The routine assumes all odd FIMPs decay to the DM, so it can overestimate yield when final-state FIMPs decay mainly to bath particles; it omits mediator decays to DM pairs and supports only 2-body final states.Source: §3.2 · pp. 19–20
- darkOmegaFi22 computes abundance from a specified 2→2 process, and can directly integrate the collision term with vegas=1 for precision checks.That integration is slow; for an s-channel resonance, production temperatures from darkOmegaFiDecay and darkOmegaFi22 correspond to mediator decay and mediator creation, respectively.Source: §3.2 · p. 20
- darkOmegaFi sums 2→2 processes with bath particles initially and at least one odd feeble particle finally, substituting decay contributions where needed to avoid on-shell poles.For two dark sectors, it computes the total abundance, with contributions separable using sort2FiDm. printChannelsFi can write channel contributions to a file.Source: §3.2 · pp. 20–21
4 Sample models and numerical results
The sample models illustrate micrOMEGAs5.0 freeze-in results and show that exact bath-particle statistics can affect predicted relic abundances, with the impact depending on the production scenario.
4.1 Vector portal
- The vector portal has off-shell, on-shell, and EFT regimes: production proceeds through annihilation, resonant mediator production, and a point-like interaction, respectively.The regimes are defined by mZ′ < 2mχ < TR, 2mχ < mZ′ < TR, and 2mχ < TR < mZ′; the EFT description assumes production starts at TR.Source: §4.1 · p. 22
- In the vector portal, exact statistics generally lower the relic density by about 4% in the EFT regime and 25% in other regimes, with a 4% increase in one on-shell case.The exception occurs when the Z′ couples mostly to SM fermions: suppressed mediator production is offset by a smaller decay width.Source: §4.1 · p. 23, Figure 3 · p. 24
4.2 Scalar portal
- In the singlet scalar model, exact bath statistics raise the abundance by nearly a factor of two when the Higgs mediator cannot be produced on shell.For light singlets produced mainly by Higgs decay, statistical effects are instead about 3.5%; the relic-density constraint can be satisfied for sub-GeV or weak-scale dark matter.Source: §4.2 · p. 24, Figure 4 · p. 25
4.3 Extended dark sector
- In the extended dark sector, heavy-lepton decays dominate at large mass splitting, while scattering takes over as the splitting decreases.Statistics matter most when decays dominate at moderate splitting; Fermi-Dirac distributions can significantly lower the relic density, while the observed abundance is obtained in the scattering-dominated regime for the chosen parameters.Source: §4.3 · p. 26, Figure 6 · p. 27
5 Installation and sample output
The section gives installation and compilation guidance for micrOMEGAs 5.0 and demonstrates a sample freeze-in relic-density calculation with its contributing channels.
5 Installation and sample output
- Users download and unpack the code, enter the micromegas_5.0 directory, select an included model directory, and compile the main code.Source: §5 · p. 27
- The sample code computes relic density through decay or scattering processes and is run with the sample data file ./main data1.par.Source: §5 · pp. 27–28
- 1.220E+00 is the sample freeze-in Omega h^2, with 6.680E-01 attributed to B, b -> ~x1,~x1.The sample output also lists the other channels contributing to the relic density.Source: §5 · p. 28
- The sample output reports an h total width of 3.067960E-03 and a branching fraction of 7.722552E-01 for h -> b,B.Source: §5 · p. 28
6 Discussion and conclusions
The paper presents micrOMEGAs methods for computing freeze-in dark matter abundances across different production scenarios, accounting for complications from high-temperature production, multiple feebly coupled particles, initial conditions, and quantum statistics. It finds that Maxwell-Boltzmann approximations can substantially alter predicted abundances and identifies limitations of the current code and areas for future development.
6 Discussion and conclusions
- Freeze-in may require tracking production back to TR, contributions from all particles in the spectrum, simultaneous equations for multiple FIMPs, and a specified initial dark matter abundance.Unlike freeze-out, freeze-in can be dominated by high temperatures, and bath-particle statistics can substantially affect the abundance.Source: §6 · p. 29
- Up to a factor 2, the dark matter abundance can be mispredicted by using Maxwell-Boltzmann bath distributions; infrared-dominated production through mediator decays peaks around mY /3.The paper presents a new micrOMEGAs version and formalism for different types of freeze-in.Source: §6 · p. 29
- For production from a heavier state that freezes in itself, the predicted abundance does not depend crucially on its kinetic equilibrium, and micrOMEGAs provides two methods to estimate this result’s validity.The code also modifies the mediator width in 2 →2 scattering when a resonance produced at Tc ≫Td decays at Td.Source: §6 · p. 29
- The paper uses “freeze-in” for both infrared- and ultraviolet-dominated scenarios, although the latter depend on the reheating temperature TR and may conflict with the mechanism’s original infrared-dominated spirit [32] [33].The authors state that the two cases do not present substantial differences from a thermodynamical standpoint.Source: §6 · p. 29
- The code does not check whether dark matter couplings satisfy the freeze-in assumption; the intermediate regime between freeze-in and freeze-out requires a more general Boltzmann equation and is planned for a future upgrade.The paper hopes the code will support study of freeze-in beyond the simplest scenarios and feebly coupled particle searches.Source: §6 · p. 30
A Treatment of t-channel propagators
At freeze-in energies, t-channel poles can destabilize 2→2 cross-section calculations; micrOMEGAs handles them with a small-angle approximation, though the method has stated limitations.
A Treatment of t-channel propagators
- At energies much larger than particle masses, t-channel poles can cause numerical instabilities in freeze-in 2→2 cross sections.For δ13 below 10^-10, double-precision integration is insufficient.Source: A Treatment of t-channel propagators · p. 30
- For t-channel diagrams, micrOMEGAs calculates the cross section for cos θ13 < 1 − 10^-10, fits small-angle values to determine A and C, then integrates the approximation symbolically.This method is intended to produce stable results.Source: A Treatment of t-channel propagators · p. 31
- The method can miss the leading contribution when the smallest-mass t-channel propagator sets δ13 but a heavier propagator dominates.A second failure occurs for l = 2 and δ13 < 10^-30, when a pole negligible near cos θ13 ≈ 1 − 10^-10 can still noticeably affect the total cross section; both issues were slated for a future release.Source: A Treatment of t-channel propagators · p. 31