SIMPonium bound states of complex scalar dark matter: Relic density and astrophysical signatures

We study a strongly interacting complex scalar dark matter candidate $(χ)$, subject to an attractive potential mediated by a vector boson $A^μ$. Such interactions allow $χ$ to form bound states: SIMPonium. In this work, we systematically investigate the bound state dynamics of our dark sector, including the formation and decay of SIMPonium and it’s influence on the thermal history. Our analysis shows in absence of any bound state formation $χ$ freezes out at $x\approx 16$ and the presence of SIMPonium in the thermal bath slightly modifies the freeze out behaviour of the free $χ$ particles, which freezes out at $x \approx 20$. While the bound state itself remains in chemical equilibrium for a longer duration and freezes out at a significantly later time, $x \approx 250$. We compute the indirect energy spectra arising from free dark matter annihilation and SIMPonium decay, where the resulting Standard Model particles subsequently produce both final state radiation and radiative decay spectra. We find that the total differential photon flux from dark matter with mass $m_χ= 150,\mathrm{MeV}$ lies in the range $E_γ^2 \frac{dϕ}{dΩ, dE_γ} \in [10^{-27},,10^{-17}] ,\mathrm{MeV,cm^{-2},s^{-1},sr^{-1}}$ , for photon energies in the interval $E_γ\in [10^{-3},,10^{2}] ,\mathrm{MeV}$. The predicted signal is therefore exceedingly feeble and remains well below the sensitivity of current experimental facilities.

💡 Research Summary

The authors investigate a strongly interacting dark sector consisting of a complex scalar particle χ that is stabilized by a Z₂ symmetry and a dark U(1) gauge symmetry. The χ particles interact via a mass‑less dark photon Aμ, which mediates an attractive Coulomb‑like potential V(r)=−α/r (α=gχ²/4π). This long‑range attraction allows χ and its antiparticle χ̄ to form bound states, dubbed “SIMPonium” (Bₙ), analogous to hydrogenic atoms. The study systematically incorporates the formation, de‑excitation, ionization, and decay of these bound states into the thermal history of the dark sector.

Model and Wavefunctions
The Lagrangian includes the kinetic term for Aμ, the covariant derivative Dμ=∂μ+igχAμ acting on χ, a quartic self‑interaction λχ|χ|⁴ that drives the 4→2 number‑changing process characteristic of the SIMP paradigm, and a Higgs portal term λχH|χ|²|H|² that provides the only connection to Standard Model (SM) fields. The authors solve the non‑relativistic Schrödinger equation for the Coulomb potential, obtaining analytic expressions for the discrete bound‑state wavefunctions (n=1–3) and the continuum scattering states. These wavefunctions are used to compute overlap integrals that determine the radiative bound‑state formation cross section σ_BSF (χχ̄→Bₙ+A) and the inverse ionization cross section.

Boltzmann System
The thermal evolution is described by coupled Boltzmann equations for the yields Yχ=nχ/s and YB=nB/s, where s is the entropy density. The equations contain: (i) the 4→2 SIMP annihilation χχχχ→χχ, (ii) the 3→2 number‑changing process χχχ→χχ, (iii) bound‑state formation and dissociation, (iv) de‑excitation of excited bound states, (v) bound‑state decay into dark photons (which are invisible) and into SM particles via the Higgs portal, and (vi) ionization processes. Thermal averages ⟨σv⟩, ⟨σv²⟩, ⟨σv³⟩ are computed analytically in the appendices and numerically integrated.

Freeze‑out Results
In the absence of bound states, χ freezes out at x≡mχ/T≈16. When SIMPonium is included, the additional radiative formation channel slightly delays the freeze‑out of free χ to x≈20, because χχ̄→Bₙ+A provides extra depletion. The bound states themselves remain in chemical equilibrium much longer, decoupling only at x≈250. However, their final yield is suppressed (YB/Yχ≈10⁻³), so the relic density is essentially set by the free χ abundance. The authors verify that for a wide parameter range (mχ up to 1 GeV, gχ≈0.5–0.9, λχ≈0.5–0.7) the observed ΩDMh²≈0.12 can be reproduced.

Indirect Detection Prospects
Two observable channels are considered: (1) annihilation of free χ into SM particles (χχ̄→SM SM) and (2) decay of SIMPonium into SM particles via the Higgs portal. The latter proceeds through Bₙ→SM SM, producing primary SM states that subsequently undergo final‑state radiation (FSR) and radiative decays (e.g., π⁰→γγ). The authors compute the full photon spectrum dΦ/dEγ for three astrophysical environments—galaxy clusters, the Galactic Center, and dwarf spheroidal galaxies—by folding the spectra with the appropriate J‑factors (for annihilation) and D‑factors (for decay). For a benchmark point mχ=150 MeV, gχ≈0.7, the resulting differential photon flux lies in the range

Eγ² dΦ/dΩ dEγ ≈ 10⁻²⁷ – 10⁻¹⁷ MeV cm⁻² s⁻¹ sr⁻¹

for photon energies from 10⁻³ MeV to 10² MeV. This signal is many orders of magnitude below the sensitivities of current γ‑ray and X‑ray telescopes (Fermi‑LAT, INTEGRAL, COMPTEL, etc.), rendering the model effectively invisible to present indirect‑detection experiments.

Conclusions and Outlook
The paper demonstrates that bound‑state effects in a SIMP framework are quantitatively important for the freeze‑out dynamics but do not dramatically alter the relic density. Moreover, the indirect‑detection signatures are extremely suppressed because the dark photon is massless and kinetically decoupled from the SM, leaving only the Higgs portal as a feeble bridge. The authors suggest that introducing a small dark‑photon mass or enhancing kinetic mixing could amplify observable signals, a direction worth exploring in future work. Overall, the study provides a comprehensive treatment of bound‑state formation, evolution, and phenomenology in a minimal SIMP model, highlighting both its theoretical consistency and the challenges it poses for experimental verification.