© The Authors 2025

1. Introduction

According to the Lambda cold dark matter (ΛCDM) paradigm, the Universe is composed of two species of matter: baryonic matter and cold dark matter. The latter’s potential interactions with the fundamental constituents of the standard model of particle physics are below the sensitivity of modern experiments, in particular its couplings with photons. Additionally, the ΛCDM includes a dark energy component responsible for the accelerated expansion of the Universe, which is parameterized by the cosmological constant Λ. Various cosmological probes have provided tight constraints on the dark matter, baryonic matter, and dark energy content of the Universe, for example, through the analysis of the cosmic microwave background (CMB) (Planck Collaboration XXIV 2016; Planck Collaboration VI 2020; Hinshaw et al. 2013), Type Ia supernovae (Abbott et al. 2019; DES Collaboration 2024; Scolnic et al. 2022; Brout et al. 2022), cluster abundance (Ghirardini et al. 2024, G24 hereafter, Bocquet et al. 2024, Costanzi et al. 2021, Lesci et al. 2022), weak lensing shear (Asgari et al. 2021; Abbott et al. 2022; Amon et al. 2022; Miyatake et al. 2023; Dalal et al. 2023), and galaxy clustering (Zhao et al. 2022; Adame et al. 2025).

Despite the constraints on the density of dark energy and dark matter obtained by state-of-the-art cosmological studies, their nature remains elusive. Consequently, potential extensions of ΛCDM with viable dark matter candidates have been explored using observations of the large-scale structure, such as axion-like particles (ALPs), and self-interacting dark matter (Randall et al. 2008; Conlon et al. 2017). The ALPs are promising potential dark matter candidates. Depending on their properties, ALPs can also partially explain dark energy. ALPs are scalar fields with typically tiny masses, compared to known elementary particles, that move in a periodic potential. Historically, the so-called quantum chromodynamics (QCD) axion was proposed by Peccei & Quinn (1977a,b) to solve a fine-tuning problem of QCD, known as the strong charge parity (CP) problem. Although predicted by QCD, the strong CP problem refers to unobserved violations of simultaneous charge conjugation (C) and parity (P). Particles with the same underlying physics naturally occur in string theory due to the compactification of large extra dimensions and are called axions or ALPs (Svrcek & Witten 2006). The typical particle masses can span the range from 10⁻³³ eV to 10⁻⁵ eV. In principle, multiple ALPs with a large spectrum of mass might exist in string theory, called the axiverse (Arvanitaki et al. 2010a). Upper bounds on the number of ALP species have been found by superradiance effects of stellar and supermassive black holes (Stott & Marsh 2018) and supernovae (Gendler et al. 2024). ALPs with masses m_a ≳ 10⁻²² eV are candidates to compose a large fraction of dark matter, while ALPs with masses m_a ≲ 10⁻³³ eV represent candidates for dark energy.

In the dark energy scenario (m_a ∼ 10⁻³³ eV), the axion field exhibits slow-roll behavior due to the small mass of ALPs. As it is still rolling today, the axion field has not yet started to oscillate around the minimum of its potential. Thus, it behaves similarly to a fluid with negative pressure, which is similar to dark energy (Hložek et al. 2015, 2018; Passaglia & Hu 2022). For masses m_a ≲ 10⁻³³ eV, the axion field is frozen in, and only its potential energy contributes to the vacuum energy of the Universe. This scenario is indistinguishable from the ΛCDM dark energy with the equation of state w = −1.

In the case of small masses around 10⁻²⁴ eV, the bosonic nature and negligible self-couplings of ALPs enable them to form Bose-Einstein condensates (BECs) on scales determined by their thermal de Broglie wavelength. The simulations by Schive et al. (2014a) have demonstrated the formation of core-like structures in the center of halos, whose size could be used to obtain constraints on the ALP mass from dwarf galaxies. Over the past decades, a number of simulations of ALPs as dark matter candidates (also called wave-like dark matter) have been performed with increasing resolution and while assuming a cosmology with a dark matter sector fully comprised by ALPs (e.g., Mocz et al. 2018; Nori & Baldi 2018; May & Springel 2021, 2023). With an increasing precision of ALP mass constraints, mixed ALP-cold dark matter simulations became relevant since they can be used to constrain the fraction of ALPs in the Universe (Schwabe et al. 2020; Laguë et al. 2024; Dome et al. 2025). In this axion mass regime, the thermal de Broglie wavelength extends to the cosmological scales, i.e., several kilo- or megaparsecons, making the axion condensates observable in studies of the matter power spectrum. A BEC on these scales forms large core-like structures, smooths out the dark matter distribution, and suppresses the formation of small-scale halos. ALPs have been proposed as a potential solution to the small-scale crisis in cosmology, particularly for addressing discrepancies in structure formation at galactic and subgalactic scales. However, the existence of the small-scale crisis itself has been questioned, with alternative explanations suggesting that baryonic physics or observational limitations may account for the perceived discrepancies (Nori & Baldi 2018; Bullock & Boylan-Kolchin 2017; De Laurentis & Salucci 2022). We refer to these cosmologically interesting low-mass axions with masses m_a ≲ 10⁻²⁴ eV as ultralight axions (ULAs). In the literature, the resulting dark matter is referred to as fuzzy dark matter due to its smoothing effect (e.g., Hu et al. 2000; Hui et al. 2017).

Constraining the coupling constant between ALPs and photons is also feasible, but it requires direct X-ray observations of the intracluster medium (ICM), which have the potential to detect a possible interaction between axions and photons. Such an interaction is only possible in magnetic fields since the axion coupling term is proportional to the scalar product of electric and magnetic fields. Measurements of spectral features in galaxy clusters have been used to constrain the properties of dark matter candidates such as ALPs and decaying dark matter (Bulbul et al. 2014; Conlon et al. 2017; Reynolds et al. 2020). Unlike the studies performed with the probes of structure growth, these analyses can determine bounds on the axion-photon coupling constant and mixing angles, but they cannot constrain the relic density of dark matter candidates.

The viable mass ranges for ULAs and abundances constituting dark matter and dark energy can be constrained through observations of the large-scale structure as well as the CMB and subgalactic scales. In the dark matter regime, where the ULAs make up all of the dark matter by assumption, strong bounds on the ULA parameter space have been reported in the literature. Lyman-alpha forest observations place a lower-mass bound of 2 × 10⁻²⁰ eV, at a 95% confidence level (Rogers & Peiris 2021). The observations of ultrafaint dwarf galaxies tighten this bound to 3 × 10⁻¹⁹ eV, at a 99% confidence level (Dalal & Kravtsov 2022), while superradiance by supermassive black holes can be used to exclude a window of masses around 7 × 10⁻²⁰ eV (Arvanitaki et al. 2015; Cardoso et al. 2018; Stott & Marsh 2018; Stott 2020; Hoof et al. 2025).

The possibility that ULAs may comprise only a fraction of the total abundance of dark matter nevertheless remains. For the first time, CMB data from the Planck satellite (Planck Collaboration XIII 2016) provided strong bounds in the mass regime 10⁻³² eV to 10⁻²⁵ eV by performing a binned analysis in different ULA mass bins. They found varying upper bounds on the ULA abundance, reaching down to Ω_a ≲ 0.01 in the intermediate mass range of 10⁻³⁰ − 10⁻²⁷ eV (Hložek et al. 2015, 2018). The same mass regime has also been constrained by Laguë et al. (2022) and Rogers et al. (2023) using Planck 2018 (Planck Collaboration VI 2020) and Baryon Oscillation Spectroscopic Survey (BOSS; (Alam et al. 2017)) data, improving existing constraints. Kobayashi et al. (2017) used the Lyman-α forest to place bounds on the fraction of ULAs in the higher-mass regime of 10⁻²³ eV to 10⁻²⁰ eV, and Winch et al. (2024) closed the gap by using UV-bright galaxy abundance to put an upper bound on the ULA fraction in the mass regime 10⁻²⁶ eV to 10⁻²³ eV. They find that ULAs cannot make up more than 22% of all dark matter in this mass range.

In this work, we use for the first time the observed abundance of galaxy clusters to constrain the parameter space of ULAs. Since the formation of halos is influenced by bosonic condensation, we expect cluster counts to be sensitive to ULAs (Diehl & Weller 2021). We particularly focus on the mass regime 10⁻³² eV to 10⁻²⁴ eV. Clusters of galaxies trace the highest peaks in the low-redshift matter density field through their number counts and therefore serve as an important late-time probe of the structure formation in the Universe. Measurements of their number count as a function of mass and redshift bin per unit volume, the so-called cluster halo mass function (HMF), is a sensitive probe of cosmology. The cluster mass function can constrain the ULA mass and abundance parameter space since the growth of structure and bosonic condensation influence the formation of halos.

To this end, we used the clusters of galaxies detected by the extended ROentgen Survey with an Imaging Telescope Array (eROSITA) on board the Spektr-RG (SRG) mission (Sunyaev et al. 2021) during its first All-Sky Survey in the western Galactic hemisphere (eRASS1). eROSITA, launched in 2019, scans the sky in the soft X-ray band with its highest sensitivity in the 0.2 to 2.3 keV energy band (Predehl et al. 2021). eRASS1 provides different cluster catalogs depending on purity and completeness requirements. The primary eRASS1 catalog comprises 12 247 galaxy clusters, with a sample purity of 86% and a coverage of 12 791 deg² (Bulbul et al. 2024; Kluge et al. 2024). In this work, we use the cosmology sample containing 5 259 securely detected clusters in the redshift range of 0.1 ≤ z ≤ 0.8, which have a sample purity of 96%. The eRASS1 cosmological sample has provided competitive results for standard cosmologies (see G24 for the details of the analysis) as well as modified gravity (Artis et al. 2024) and the growth of structures (Artis et al. 2025). The eRASS1 cosmology sample is the largest ICM-selected galaxy cluster sample ever used to infer cosmology from number counts. The covered cluster mass range extends down into the galaxy group regime, opening the possibility of probing cosmological models sensitive to below-cluster scales. The presence of ULAs has the strongest imprint on the number density of low-mass halos, making the eRASS1 cosmology sample a promising tool to constrain the properties of ULAs. In contrast to the main cosmological inference G24, in this work, to constrain the ULA parameter space, we used the Boltzmann solver axionCAMB instead of CAMB to compute a matter power spectrum that is sensitive to two ULA parameters, particularly the mass of the ULA, m_a, and the energy density fraction of ULAs, Ω_a, in units of the critical density of the Universe.

This paper is organized as follows. A more detailed overview of the ULAs and their effect on cosmology is given in Sect. 2, including their influence on structure formation. The underlying cosmological sample of the eRASS1 survey and the optical and weak lensing data used for the analysis are summarized in Sect. 3. Section 4 reviews the eRASS1 cluster cosmological pipeline for Bayesian inference in light of the ULA parameters. The results are presented and discussed in Sects. 5 and 6. For better readability, we use natural units c ≡ ℏ ≡ 1 throughout this work.

2. Ultralight axion cosmology

In this section, we first provide a review of the motivation for ULAs as dark matter candidates. We then explain the implications of their existence on the hierarchical structure formation and cosmology.

2.1. Physical origin of axions

Axions are hypothetical pseudoscalar particles with typically low mass (between 10⁻³³ eV and 10⁻⁵ eV) and weak couplings to the Standard Model particles (i.e., quadratic coupling terms can be neglected). Pseudoscalar hereby refers to the property that the particle behaves like a scalar under transformations, except parity transformations, under which it switches signs. Although the axion was originally introduced in QCD to solve the strong CP problem explicitly, the term has been extended to a number of (pseudo) Nambu-Goldstone bosons. Nambu-Goldstone bosons are particles that occur as massless degrees of freedom after breaking a spontaneous symmetry. Axions are pseudo-Nambu-Goldstone bosons as they acquire a very small mass via further symmetry breaking. These particles appear in various theories of high energy physics such as string theory (e.g., Svrcek & Witten 2006; Conlon 2006; Arvanitaki et al. 2010b) or supersymmetric extensions of the Standard Model (e.g., Kim & Carosi 2010; Marsh 2016, for review). To distinguish such general axions from the QCD axion, they are often referred to as ALPs or ultralight axions in the cosmologically relevant low-mass regime. The common properties of these particles are well described in Section 2.2 in Marsh (2016).

We heuristically illustrate axion physics using the QCD axion introduced by Peccei & Quinn (1977a,b) as a generic example. We refer the reader to Marsh (2016) for further details. A non-Abelian SU(3) gauge theory with a gluon gauge field $G_{\mu \nu }^{\mathrm{a}}$ naturally exhibits a topological term with some coupling θ in its Lagrangian density, which is of the form

$$\begin{array}{c}\hfill \mathcal{L}\supset \frac{\theta }{32{\pi }^2}{G}_{\mu \nu }^a{\stackrel{\sim }{G}}^{\mathrm{a}\mu \nu },\end{array}$$(1)

where ${\stackrel{\sim }{G}}^{\mathrm{a}\mu \nu }={ϵ}^{\mu \nu \alpha \beta }{G}_{\alpha \beta }^{\mathrm{a}}/2$ is the dual of the gluon gauge field tensor, μ and ν refer to spacetime coordinates and a is some color index. This term is CP-violating and cannot be transformed away. If QCD is CP-violating, it leads to a non-vanishing neutron electric dipole moment d_n ≈ θ (3.6 × 10⁻¹⁶e cm) (Crewther et al. 1979). However, experiments place strong constraints on the neutron electric dipole moment (Abel et al. 2020), which results in a fine-tuned value of θ ≲ 10⁻¹⁰. In other words, QCD appears to be non-CP-violating, although it should. This is known as the strong CP problem.

Peccei & Quinn (1977a,b) solved this problem by generalizing the coupling constant θ to a scalar field that naturally cancels any contribution to the CP-violating term. To this end, they introduced a new global U_PQ(1) symmetry, known as the Peccei-Quinn symmetry. The spontaneous breaking of this symmetry at a high energy scale f_a leads to a vacuum manifold that is isomorphic to the one-sphere S¹. The emerging Nambu-Goldstone boson (which will become the axion later) is an angular degree of freedom in an effectively flat potential. Thus, at this point, the new particle is massless. At some lower energy scale, Λ_a (not to be confused with the cosmological constant Λ), non-perturbative QCD effects (e.g., instantons, topological effects) break the flatness of the vacuum manifold, the massless Nambu-Goldstone boson lives in. Due to the angular nature of the Nambu-Goldstone boson, only a discrete angular shift symmetry is preserved. Ultimately, the Nambu-Goldstone boson acquires a mass and is canonically called an axion. Simplifying further, the field rolls into the circular minimum of a sombrero-like potential during the first spontaneous symmetry breaking. The Nambu-Goldstone boson is the angular degree of freedom that moves freely along the minimum of the potential. The angular shift symmetry of this potential is then broken by nonperturbative effects, leading to local minima and maxima along the former circular minimum of the potential. Particles moving in a potential with minima acquire a mass; thus, the Nambu-Goldstone boson becomes a massive axion particle.

The axion mass is related to the symmetry-breaking scales f_a and Λ_a. A typical but not unique effective axion potential can take the form

$$\begin{array}{c}\hfill V(\varphi )={\mathrm{\Lambda }}_{\mathrm{a}}^4[1-cos\left(\frac{\varphi }{f_{\mathrm{a}}}\right)].\end{array}$$(2)

Considering small oscillations around one of the local minima gives the effective mass of the axion by approximating the local environment of the minimum with a quadratic function:

$$\begin{array}{cc}\hfill V(\varphi \approx {\varphi }_{\min })& \approx \frac{1}{2}{\left(\frac{{\mathrm{\Lambda }}_{\mathrm{a}}^2}{f_{\mathrm{a}}}\right)}^2{\varphi }^2+{\mathrm{\Lambda }}_{\mathrm{a}}^4\mathcal{O}\left(\frac{{\varphi }^4}{f_{\mathrm{a}}^4}\right)\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill & =\frac{1}{2}{m}_{\mathrm{a}}^2{\varphi }^2+\left(\text{self-interactions}\right).\hfill \end{array}$$(4)

The mass of the axion can be read off as m_a = Λ_a²/f_a. Typical values of f_a lie between 10⁹ GeV and 10¹⁷ GeV in the QCD case and, in the general case, above 10¹⁰ GeV with typical values around the grand unified theory (GUT) scale at 10¹⁶ GeV (Marsh 2016). In general, self-interactions are thus suppressed by m_a²/f_a² ≲ 10⁻⁹⁰ for the ULA mass regime considered in this work, and we make the assumption that self-interactions are negligible in our analysis.

The underlying mechanisms of the QCD axion also apply to axions in string theory models. However, in string theory, the periodicity of the axion potential originates from the compactification of extra dimensions, i.e., extra dimensions are wrapped up, which makes them unobservable by macroscopic observers. The need for extra dimensions in string theory, which have not been observed yet, leads to a natural production mechanism of axion particles with a potentially wide range of possible axion masses.

2.2. Cosmological implications of ultralight axions

This section reviews the implications of an ultralight axion fluid on cosmology. For in-depth discussion and derivations, we refer the reader to the derivations presented in Marsh (2016), Hložek et al. (2015), and Hložek et al. (2018). The ULA field enters the Einstein-Hilbert action of general relativity as a massive scalar field in curved spacetime. On a Friedmann-Lemaître-Robertson-Walker metric background, the Klein-Gordon equation (equation of motion) of the ULA is

$$\begin{array}{c}\hfill \ddot{\varphi }+3H\dot{\varphi }+m_{\mathrm{a}}^2{\varphi }^2=0.\end{array}$$(5)

Here, $H(t)=\dot{a}(t)/a(t)$ is the Hubble parameter, ϕ is the ULA field, and dots refer to derivatives in time. This equation describes a harmonic oscillator in ϕ with a damping term proportional to the Hubble parameter H. For this illustration, we neglect any source terms from other cosmological fluid components or self-interactions, which would enter the equation on the right-hand side and replace the zero. For H ≪ m_a, i.e., if the ULA mass is orders of magnitude larger than the Hubble parameter, the damping term is negligible, and the ULA field oscillates freely. If the ULA’s mass is sufficiently large, this describes the matter-like late-time behavior of the ULA. For the opposite case, H ≫ m_a, the field is overdamped and frozen in, showing the same behavior as dark energy. The transition between the two regimes defines a scale a_osc which can be approximated by H(a_osc)≈m_a. In this regime, the ULA is slow-rolling. Thus, a single ULA species may transition from a dark energy-like behavior to a dark matter-like behavior. As the universe evolves in the opposite direction from matter-dominated to dark energy-dominated era, no single ULA species can describe both dark matter and dark energy. We therefore distinguish two distinct phenomenological regimes for ULAs: the dark matter regime, where they behave as dark matter, and the dark energy regime, where they show the same observed effect of dark energy. We review the two regimes in the following subsections.

2.3. Dark matter regime of ultralight axions

Introducing perturbations to the ULA field and choosing a gauge yields an equation for the evolution of perturbations in terms of the ULA overdensity δ_a = δρ_a/ρ_a. In the simple case of a ULA-dominated universe, it takes the form

$$\begin{array}{c}\hfill {\ddot{\delta }}_{\mathrm{a}}+2H{\dot{\delta }}_{\mathrm{a}}+(\frac{k^2{c}_{\mathrm{s},\mathrm{a}}^2}{a^2}-4\pi G{\rho }_{\mathrm{a}}){\delta }_{\mathrm{a}}=0,\end{array}$$(6)

where k is the wave number, a the scale factor, ρ_a describes the background density field of the ULA, δρ_a its perturbation. Other potential fluid components would add source terms to the equation. The effective sound speed of ULAs (c_s, a) becomes

$$\begin{array}{c}\hfill c_{\mathrm{s},\mathrm{a}}^2=\frac{k^2}{k^2+4m_{\mathrm{a}}^2{a}^2}.\end{array}$$(7)

Equation (6) defines a scale-dependent Jeans scale k_J for ULAs. As described in Marsh (2016), large-scale perturbations (k < k_J) grow and behave similarly to CDM. Small-scale perturbations (k > k_J) oscillate, making them distinct from CDM perturbations on the same scales. As described in Marsh et al. (2012), modes entering the horizon while the ULA sound speed in Eq. (7) is near the speed of light will suppress structure formation, while modes entering the horizon while the ULA sound speed is approaching zero will cluster like ordinary cold dark matter. Effectively, this translates to suppressing the matter power spectrum on small physical scales (large values of the wave number k).

The properties of ULAs and their cosmological effects can be described with two quantities: The particle mass of the ULA m_a, and the relative energy density of ULAs with respect to the critical density of the universe, Ω_ax = ρ_ax/ρ_crit, also known as the relic density of ULAs. Here, the critical density of the universe is defined as ρ_crit = 3H₀²/(8πG). The ΛCDM cosmology is recovered when Ω_a → 0. Furthermore, the ULAs’ suppression imprint on the matter power spectrum becomes negligible for high ULA masses m_a ≫ 10⁻²² eV, effectively recovering a ΛCDM cosmology in the context of this analysis. Fig. 1 demonstrates the suppression effect of ULA dark matter with different logarithmic particle masses log₁₀(m_a[eV]) and relic densities Ω_a on the matter power spectrum. Suppression of the power can be seen on small scales (large k), whereby the effect is stronger with higher ULA fractions or lower ULA masses.

Fig. 1. Effect of ULAs on the matter power spectrum recovered assuming the Planck Collaboration VI (2020) cosmological values (ΩΛ ≈ 0.7, Ωm + Ωa ≈ 0.3) at z = 0.1. On the left, the change in the matter power spectrum with varying ULA mass is shown when the relative ULA density is fixed to Ωa/(Ωcdm + Ωa) = 0.5. The panel on the right displays the dependence of the power spectrum on varying ULA abundance at a fixed ULA mass of ma = 10−25 eV. The concordance ΛCDM model is shown in a black curve on both panels. Small-scale suppression becomes stronger with decreasing ULA mass and increasing ULA abundance. For higher ULA masses of ma ≫ 10−22 eV, the matter power spectrum becomes indistinguishable from the one in a ΛCDM cosmology, as is the case for Ωa → 0.

In this work, we explore whether a fraction of dark matter is composed of ULAs of an unknown mass in the dark matter regime. We still need to postulate a cosmological constant, dark energy, in this context. As we assume a flat universe, the dark energy density is given by Ω_Λ = 1 − Ω_m − Ω_a, while we fit for the total matter density Ω_m and the ULA density Ω_a. The results are provided in Sect. 5.

2.4. Dark energy regime of ultralight axions

If the ULA is so light that its mass becomes comparable with or smaller than the Hubble parameter H₀, the oscillation scale a_osc approaches or even exceeds the current scale factor of the universe a(z = 0) = 1, indicating that the ULA field has not started oscillating yet. Indeed, this kind of field behaves as dark-energy-like, which can be heuristically illustrated with the equation of motion. Neglecting the m_a²ϕ² term in Eq. (5) and assuming that H(t) varies slowly enough in time exhibits exponential modes for the ULA field ϕ. This behavior is associated with an equation of state w_ax ≈ −1 via the density and pressure of the ULA field. In terms of the de Broglie wavelength associated to the ULA BEC, in the extreme case of H₀ ≈ m_a which occurs around ULA masses of m_a ≈ 10⁻³³ eV the de Broglie wavelength is comparable to the size of the universe leading to an overall energy shift, typically associated to a dark energy-like component. However, even at higher ULA masses, as long as the de Broglie wavelength of the ULA BEC exceeds its Jeans scale, ULAs cannot form halos and are thus still referred to as dark energy-like ULAs. Furthermore, if the ULA BEC de Broglie wavelength is comparable to the size of the universe, this affects the Hubble rate and therefore also the distance measurements. This results in a redshift-dependent transition in the Hubble rate, with the nature of the transition dependent on the ULA mass. Wherever the de Broglie wavelength exceeds the size of the universe, ULAs should be modeled as a dark energy component in the Hubble rate. In the opposite case of a de Broglie wavelength smaller than the size of the universe, ULAs should be modeled as a dark matter component in the Hubble rate. The transition between the two regimes is determined by the oscillation redshift. Fig. 2 illustrates the dependence of the oscillation redshift on the ULA mass. We emphasize here that the regime change in the Hubble rate is not related to the typically referenced regime change between dark energy-like ULAs’ and dark matter-like ULAs’ contribution to halo formation, which is determined by the respective Jeans scale.

Fig. 2. Demonstration of the ULA mass and redshift dependent regime change in the Hubble rate. In the pink regime, ULAs contribute to the Hubble rate e in a manner similar to dark matter; in the blue region, ULAs contribute to the Hubble rate in a dark energy-like fashion. The orange shaded region shows the eRASS1 cluster redshift range, while the gray shaded region shows the redshift range of lensed galaxies from the DES, HSC, and KiDS data, described in Sect. 3.2. Above a ULA mass of ma ∼ 10−31.5 eV, ULAs should be modeled as dark matter in the Hubble rate, resulting in correct distance measurements during the parameter inference. For ULAs with masses below ma ∼ 10−31.5 eV, the distance measurements are affected by the regime change via the Hubble rate, where ULAs are treated as dark energy.

3. Multiwavelength data

3.1. eROSITA cluster catalog and X-ray observable

This work uses the cosmology sample compiled from the first eROSITA All-Sky Survey (Merloni et al. 2024; Bulbul et al. 2024; Kluge et al. 2024). This sample is constructed by employing an extent likelihood cut of ℒ_ext > 6 and a redshift range of 0.1 ≤ z ≤ 0.8. The publicly available DESI Legacy Survey Data Release 10 (LS DR10) data are used for the confirmation, identification, and photometric redshift measurements of the X-ray-selected cluster candidates detected in the eROSITA All-Sky Survey’s western Galactic hemisphere (Kluge et al. 2024). The common footprint between eRASS1 and the LS DR10-South footprint is 12 791 deg², and contains 5,259 securely confirmed clusters (Kluge et al. 2024). Its purity level is estimated to be 96%.

The eROSITA X-ray data are reprocessed with the eSASS software as described in Merloni et al. (2024) by applying a correction for the Galactic absorption and a more accurate ICM and background modeling. The details of the X-ray processing are provided in Bulbul et al. (2024), Liu et al. (2022). The X-ray observable count rate in the observer frame, which demonstrates a tight correlation with weak lensing shear data and is related to selection, is used as a mass proxy for this work, similar to the method adopted in G24.

3.2. Weak lensing data

In order to minimize the bias on galaxy cluster masses, we use weak gravitational lensing data to obtain reliable mass measurements for the calibration of the mass-observable scaling relations. We used weak lensing measurements in the form of tangential reduced shear profiles extracted from the shapes of background galaxies in the Dark Energy Survey year three (DES Y3) data (Sevilla-Noarbe et al. 2021; Gatti et al. 2021) and the Kilo Degree Survey (KiDS), and the Hyper Suprime-Cam (HSC) Strategic Survey program (Hildebrandt et al. 2021; Li et al. 2022). The details of shape measurements around the detected eROSITA clusters and the weak lensing mass calibration details are provided in Grandis et al. (2024), Kleinebreil et al. (2025) and Chiu et al. (2025). Besides the reduced tangential shear profile, we also use estimates of the redshift distribution for the employed background sources, estimates of the measurement uncertainty on the reduced shear profile, and calibrations on the contamination of the background sample by cluster member galaxies.

The Dark Energy Survey shares a common sky area of 4,060 deg² with eROSITA, and its observations have been conducted in the r, i, and z bands (Sevilla-Noarbe et al. 2021). The details of the application of the DES Y3 shear maps to the eRASS1 galaxy clusters are presented in Grandis et al. (2024). Tangential shear profiles were derived from the DES Y3 shape catalog (Gatti et al. 2021). The analysis yielded 2201 tangential shear profiles for eRASS1 galaxy clusters with a signal-to-noise ratio of 65.

The Hyper Suprime-Cam Subaru Strategic Program (Aihara et al. 2018) is a wide and deep optical survey in the g, r, i, z, and y bands. The three-year shape catalog (HSC Y3) (Li et al. 2022) was used to analyze galaxy cluster shear profiles in an overlapping sky area of ≈500 deg² between HSC Y3 and eRASS1. This analysis yielded tangential shear profiles, lensing covariance matrices, and redshift distributions for 96 eRASS1 galaxy clusters. The total signal-to-noise is 40 (for details see Chiu et al. 2025; Okabe et al. 2025; Chiu et al. 2022).

The Kilo-Degree Survey is an optical wide-field survey in the u, g, r, and i bands, dedicated to delivering weak gravitational lensing as well as photometric redshift measurements (de Jong et al. 2013). We use the shear maps and photometric redshifts of the fourth data release of the Kilo-Degree Survey (KiDS-1000) over a sky area of ≈1000 deg² (Kuijken et al. 2019; Giblin et al. 2021; Hildebrandt et al. 2021; Wright et al. 2020). The joint sky coverage with eROSITA yielded reduced tangential shear maps for 237 eRASS1 galaxy clusters (101 in the KiDS-North field, 136 in the KiDS-South field) with a total signal-to-noise of 19 (for details see Kleinebreil et al. 2025).

We highlight the importance of the overlap between eRASS1 clusters in KiDS and HSC (125 clusters) and in KiDS and DES (25 clusters). There is no overlap in the sky coverage of HSC and DES, so the reduced tangential shear maps deduced from KiDS-1000 provide an internal consistency check (Kleinebreil et al. 2025, G24).

4. Methodology

This section briefly presents the inference pipeline developed in G24 and explains the applied adaptations to consider the cosmological effects of introducing ULAs. To this end, we reviewed the observables, the HMF, the mass-observable scaling relations, and the weak gravitational lensing calibration.

The forward modeling cosmological pipeline is described in G24. It takes into account selection effects through the selection function presented in Clerc et al. (2024) to model the completeness and mass-observable scaling laws calibrated with weak gravitational lensing data by DES Y3, KiDS, and HSC Y3, presented in Grandis et al. (2024), Kleinebreil et al. (2025), Chiu et al. (2025), as well as a mixture model for eliminating care of possible contaminants (Kluge et al. 2024). The Bayesian analysis assumes a Poissonian likelihood for the galaxy cluster number counts. The pipeline is used for constraining a large set of cosmological parameters as well as scaling parameters simultaneously.

4.1. Statistical inference

This section summarizes the framework presented in G24 for the computation of the cluster count likelihood. The cluster abundance, n, per units of true mass, M; true redshift, z; and solid angle (true sky positions are noted with ℋ) follows

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}lnM\mathrm{d}z\mathrm{d}\mathcal{H}}=\frac{{\rho }_{\mathrm{m},0}}{M}\frac{\mathrm{d}ln{\sigma }^{-1}}{\mathrm{d}lnM}f(\sigma )\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathcal{H}},\end{array}$$(8)

where ρ_m, 0 is the matter density at present, σ(R, z) is the root mean square density fluctuation defined below, d V/d zd ℋ is the differential comoving volume per redshift per steradian. Further, f(σ) is the multiplicity function introduced by Tinker et al. (2008), with its parameters assumed to be fixed and universal. The root mean square density fluctuation, σ, is defined as the variance of a convolution of the square of a window function, W, with the matter power spectrum, P:

$$\begin{array}{c}\hfill {\sigma }^2(z,R)={\displaystyle {\int }_0^{\infty }}\mathrm{d}k{k}^{2}P(k,z){|W(kR)|}^2,\end{array}$$(9)

where k is the wave number of the considered mode and R is the defining radius of the corresponding overdensity. For W, the standard top-hat window function is used, which defined as the Fourier transform of a radially symmetrical uniform density distribution up to a certain radius:

$$\begin{array}{c}\hfill W(kR)={\displaystyle {\int }_{{\mathbb{R}}^3}}{\mathrm{d}}^3\mathbf{r}\frac{3\mathrm{\Theta }(R-r)}{4\pi R^3}{e}^{-i\mathbf{k}\mathbf{r}}=3\frac{sin(kR)-(kR)cos(kR)}{{(kR)}^3}.\end{array}$$(10)

The observed cluster number density is assumed to follow a Poisson statistic, as commonly performed in modern cluster surveys (e.g., Bocquet et al. 2024). We do not consider the sample variance (Hu & Kravtsov 2003), as its impact is negligible given the characteristics of the eRASS1 survey (i.e., a catalog of massive halos and a large volume, see Fumagalli et al. 2021). We note as λ the intensity of the Poisson process (the number density of objects per unit of observable) and x as the vector of the observables. This intensity depends on the background cosmological model. We note with Θ the set of cosmological parameters used. By definition, the expected number of objects whose observable properties belong to the subsample Ω of the observable parameter space follows

$$\begin{array}{c}\hfill N_{\{x\in \mathrm{\Omega }\}}(\mathrm{\Theta })={\displaystyle {\int }_{\mathrm{\Omega }}}\lambda (x|\mathrm{\Theta })\mathrm{d}x.\end{array}$$(11)

In the case of eRASS1, combined with the observation of the Legacy Survey DR10 (Kluge et al. 2024), the observable vector incorporates the observed X-ray count rate, $\widehat{C_{\mathrm{R}}}$; the observed optical richness, $\widehat{\lambda }$; the photometric redshift, $\widehat{z}$; the observed sky position, $\widehat{\mathcal{H}}$; and the reduced tangential shear profile, $\widehat{g_{\mathrm{t}}}$, for clusters that belong to the overlapping weak lensing surveys used for the mass calibration. The vector of observables consequently reads as

$$\begin{array}{c}\hfill x=\{{\widehat{C}}_{\mathrm{R}},\widehat{z},\widehat{\lambda },\widehat{\mathcal{H}},\widehat{g_t}\}.\end{array}$$

We can thus rewrite the intensity of the Poisson process as $\lambda (x|\theta )=(\mathrm{d}\widehat{n})/(\mathrm{d}\widehat{C_{\mathrm{R}}}\mathrm{d}\widehat{\lambda }\mathrm{d}\widehat{g_{\mathrm{t}}}\mathrm{d}\widehat{z}\mathrm{d}\mathcal{H})$, with

$$\begin{array}{c}\hfill \frac{\mathrm{d}\widehat{n}}{\mathrm{d}\widehat{C_{\mathrm{R}}}\mathrm{d}\widehat{\lambda }\mathrm{d}\widehat{g_{\mathrm{t}}}\mathrm{d}\widehat{z}\mathrm{d}\mathcal{H}}={\displaystyle \int \frac{\mathrm{d}n}{\mathrm{d}\ln M\mathrm{d}z\mathrm{d}\mathcal{H}}}\mathcal{P}(x|M,z)\mathcal{P}(I|x)\mathrm{d}\ln M\mathrm{d}z,\end{array}$$(12)

where 𝒫(x|M, z) is the probability distribution function related to the observables at a given mass and redshift (see Sect. 4.3), and 𝒫(I|x) is the selection function, describing the expected fraction of objects that we detect for a given set of observable x, described in Sect. 4.2.

4.2. Selection function

The selection function model 𝒫(I|x) was produced using the state-of-the-art eRASS1 digital twin simulation (Seppi et al. 2022; Comparat et al. 2020). The simulation reproduces the AGN and cluster distributions with high accuracy and produces synthetic eRASS1 observations. A detection probability is assigned to each set of observables using a Gaussian process classifier described in Clerc et al. (2024). In practice, the observables used to assign the probability are the true redshift z, the galactic column density n_H-corrected count rate C_R, and the sky position ℋ. The selection function depends on the sky position through the local background surface brightness, the local hydrogen column density, and the exposure time, which are not uniform across the eRASS1 sky. The selection function enters the statistical inference framework as described in Sect. 4.1.

4.3. Scaling relations

The eRASS1 cluster abundance pipeline primarily considers two mass proxies: the X-ray count rates, C_R, and the optical richness, λ. The scaling relations between these observables and the underlying dark matter halo mass are fitted together with the cosmological parameters. The count rate scaling relation is assumed to follow

$$\begin{array}{c}\hfill \langle ln\frac{C_{\mathrm{R}}}{C_{\mathrm{R},\mathrm{p}}}|M,z\rangle =ln{A}_{\mathrm{X}}+b_{\mathrm{X}}(z)ln\frac{M}{M_{\mathrm{p}}}+e_{\mathrm{x}}(z),\end{array}$$(13)

where the pivot count rate, pivot mass, and pivot redshift are fixed to C_R, p = 0.1 cts / s, M_p = 2 × 10¹⁴ M_⊙, and z_p = 0.35 by choice. The other terms follow

$$\begin{array}{c}\hfill b_{\mathrm{X}}(z)=B_{\mathrm{X}}+F_{\mathrm{X}}ln\frac{1+z}{1+z_{\mathrm{p}}}\end{array}$$(14)

and

$$\begin{array}{c}\hfill e_{\mathrm{x}}(z)=D_{\mathrm{X}}ln\frac{d_{\mathrm{L}}(z)}{d_{\mathrm{L}}(z_{\mathrm{p}})}+E_{\mathrm{X}}ln\frac{E(z)}{E(z_{\mathrm{p}})}+G_{\mathrm{X}}ln\frac{1+z}{1+z_{\mathrm{p}}}.\end{array}$$(15)

The richness scaling relation has a form similar to

$$\begin{array}{c}\hfill \langle ln\lambda |M,z\rangle =ln{A}_{\lambda }+b_{\lambda }(z)ln\left(\frac{M}{M_{\mathrm{p}}}\right)+C_{\lambda }ln\left(\frac{1+z}{1+z_{\mathrm{p}}}\right),\end{array}$$(16)

and the redshift dependence of the mass slope follows

$$\begin{array}{c}\hfill b_{\lambda }(z)=B_{\lambda }+D_{\lambda }ln\left(\frac{1+z}{1+z_{\mathrm{p}}}\right).\end{array}$$(17)

These scaling relations bridge the X-ray, optical, and shear data. They enter into the statistical inference framework via 𝒫(x|M, z), as described in Sect. 4.1.

4.4. Weak lensing calibration

This section demonstrates that the modeling introduced for the cluster mass calibration in G24 can be used in this work without further modifications. Future X-ray samples will explore lower-mass ranges and will require a deeper understanding of the lensing of dark matter halos in the presence of ULA. In the context of eRASS1 statistical precision, we demonstrate that adopting a generic NFW profile is sufficient. We begin by showing that the weak lensing mass bias model, which characterizes the discrepancy between the true total mass of a cluster and the weak lensing-inferred mass employed in this analysis, is largely insensitive to the ULA mass fraction. In practice, to accurately predict the mapping between the halo mass and the weak lensing profiles, we employ the shear inferred cluster mass, which results from fitting a reduced shear profile with our custom shear profile model, as proposed by Grandis et al. (2021). However, the shear inferred masses are biased compared to the true total mass of a cluster (Becker & Kravtsov 2011). The relation between weak lensing mass and true dark matter halo mass is parametrized with the weak lensing bias and weak lensing scatter. The values of the latter, as well as their uncertainties, are calibrated on the Monte-Carlo realizations of synthetic shear profiles (Grandis et al. 2024). These synthetic shear profiles are based on mass maps from hydrodynamical simulations augmented by the eRASS1 miscentering distribution, by the cluster member contamination results, by the shape and photometric redshift calibration uncertainties of respective weak lensing surveys, and by an uncorrelated large-scale structure noise contribution. We assume that these calibrations will remain valid over the full ULA mass range for small ULA fractions Ω_a/(Ω_cdm + Ω_a)≲0.1, since simulations do not show evidence of ULA effects on halo density profiles below this threshold (Schwabe et al. 2020; Vogt et al. 2023). As discussed below, even for a dark matter sector comprised solely by ULAs (Ω_a/(Ω_cdm + Ω_a) = 1), the calibrations remain valid for m_a ≳ 10⁻²⁶ eV. Consequently, the bias model remains the same as the one used in G24.

Furthermore, we show that for an eRASS1-like sample and a cluster abundance-based analysis, the impact of ULAs on the dark matter density profile used in weak lensing mass estimation can be safely neglected. However, for deeper eROSITA surveys, this assumption may break down due to the increased statistical contribution from group-sized halos. In general, the presence of ULA dark matter impacts the structure of halos, as explored most completely in gravity-only simulations in May & Springel (2023). Those simulations confirm prior phenomenological findings of Schive et al. (2014b), Marsh (2016) that ALPs form a so-called soliton core where quantum degeneracy pressure compensates for the inward pull of gravity in the inner regions of halos. The density profile of the core can be analytically approximated as

$$\begin{array}{c}\hfill {\rho }_{\text{soli}}(r)=\frac{1.9\times {10}^9(1+z){\left(\frac{m_{\mathrm{a}}}{{10}^{-23}\text{eV}}\right)}^{-2}{\left(\frac{r_{\text{soli}}}{\text{kpc}}\right)}^{-4}}{{(1+0.091{\left(\frac{r}{r_{\text{soli}}}\right)}^2)}^8}{\text{M}}_{\odot }{\text{kpc}}^{-3},\end{array}$$(18)

where z is the halo redshift. The characteristic radius of the soliton cores r_soli marks the scale at which the density profile transitions from a flat inner part to a steeply declining outer part. Such core profiles for different ULA particle masses m_a, and halo masses M at redshift z = 0.1 are shown in the left panel of Fig. 3 as dotted lines. Crucial to the phenomenology of halo profiles in ULA cosmology is the halo mass dependence of the core radius, given by Schive et al. (2014b) as

$$\begin{array}{c}\hfill r_{\text{soli}}=1.6{\left(\frac{m_{\mathrm{a}}}{{10}^{-22}\text{eV}}\right)}^{-1}{(1+z)}^{1/2}{\left(\frac{\zeta (z)}{\zeta (0)}\right)}^{1/6}{\left(\frac{M}{{10}^9{\text{M}}_{\odot }}\right)}^{-1/3}\text{kpc},\end{array}$$(19)

where ζ(z) = (18π² + 82(Ω_m(z)−1)−39(Ω_m(z)−1)²)/Ω_m(z), which we evaluate with the present day matter density Ω_m(z = 0) = 0.3. As shown in Fig. 3, smaller ULA particle masses lead to larger cores, which are also larger in lower-mass halos.

Fig. 3. Density profiles (left column) and projected density contrasts (right column) of halos of different mass (color coded) at redshift z = 1 for different ULA particle masses (rows) in the extreme case of a ULA-only dark matter sector. Full lines show the total profiles, dashed lines show the corresponding cold dark matter profiles (assumed to follow a Navarro-Frenk-White profile), and dotted lines show the soliton core profile. In the right column, we highlight the scales we use for the weak lensing measurement in white.

In this section and for this demonstration, we use the Navarro-Frenk-White profile for the collisionless part of the profile (Navarro et al. 1996) given by

$$\begin{array}{c}\hfill {\rho }_{\text{NFW}}(r)=\frac{{\rho }_0}{\frac{r}{r_{\text{s}}}{(1+\frac{r}{r_{\text{s}}})}^2},\end{array}$$(20)

with scale radius r_s and normalization ρ₀. We set these values as follows: For a given halo mass and redshift, we computed the radius r_500c using the definition for spherical overdensity masses and the mean concentration c_500c following the relation by Ragagnin et al. (2021). This sets r_s = r_500c/c_500c. Using the profile to compute the mass enclosed in r_500c, we set the normalization constant ρ₀ in accordance with Eq. 4 in Navarro et al. (1996) adjusted to the overdensity 500.

We constructed the composite model following the recommendations by May & Springel (2023). The transition radius r_t from a soliton core to the collisionless profiles is defined by the condition ρ_NFW(r_t) = ρ_soli(r_t). If for small ULA masses (m_a < 10⁻²⁴ eV) no r_t satisfies this condition, we pick r_t = max(r²ρ_soli). If there are two solutions, we picked the larger of the two. The total profile was constructed as

$$\begin{array}{c}\hfill \rho (r)=\{\begin{array}{c}{\rho }_{\text{soli}}(r)\text{if}r\le r_{\text{t}}\hfill \\ {\rho }_{\text{NFW}}(r)\text{if}r>r_{\text{t}}.\hfill \end{array}\end{array}$$(21)

May & Springel (2023) discusses that this model is only qualitatively correct for spherically averaged profiles. As such, it is not suited for our precise weak lensing measurements and can only be used to understand the phenomenological effect of ULA particles on massive halo profiles. In the left column of Fig. 3, we show the resulting total mass profiles for different ULA masses m_a and halo masses M.

However, the density profile is not directly observable by the weak gravitational lensing signal, which is sourced by the projected density contrast

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R)=⟨\mathrm{\Sigma }(<R)⟩-\mathrm{\Sigma }(R),\end{array}$$(22)

where R is the 2D cluster centric distance, and $\mathrm{\Sigma }(R)={\displaystyle \int \text{d}z}\rho (\sqrt{R^2+z^2})$, the projected matter density. The transformation from a 3D density profile to a density contrast is nonlocal, as weak lensing depends on the enclosed mass. Flat density profiles thus provide very little signal, while the effect of central overdensity can be observed well outside of their physical extent. These density contrasts are shown in the right panel of Fig. 3.

In this extreme case of a ULA-only dark matter sector, two interesting regimes can be observed. For intermediate ULA masses around m_a = 10⁻²⁴ eV pronounced soliton cores in the innermost part of the halos appear. Nonetheless, this effect is not observable by our weak lensing measurements, which are restricted to scales around ∼1 Mpc (highlighted in white in the right panel of Fig. 3). Small ULA particle masses m_a ≲ 10⁻²⁶ eV are a second interesting regime. For these ULA masses, cores become so extended that they dominate the density profile for group-scale halos. At m_a ∼ 10⁻²⁷ eV, also cluster scale halos with M < 10¹⁴ M_⊙ are dominated by cores. This observation can be interpreted in terms of the soliton radius. If it becomes larger than the cluster size, axions no longer cluster in the halo. This is the ULA mass range we refer to as the dark energy regime. In this context, the light axions, because of their flat inner profile shape, provide a vanishing weak lensing signal. In all cases, the excess surface mass density described by Equation (22) is mostly unaffected. Furthermore, the effect on ΔΣ is, as expected, more noticeable for low-mass halos (M ∼ 10¹³ M_⊙). The eRASS1 cosmological sample has a median mass of M_med = 10^14.39 M_⊙ while 90% of the halos have a mass larger than 10¹⁴ M_⊙. We therefore conclude that even a pure ULA dark matter sector for m_a ≳ 10⁻²⁶ eV does not significantly impact the prediction for the weak lensing signal of group and cluster scale halos, as compared to a Navarro-Frenk-White profile. Quantitative predictions for the weak lensing signal of groups and clusters would require large-volume hydro-dynamical simulations in ULA cosmology, which are not available to date, and therefore cannot be incorporated in this work. These considerations are left for future works.

Finally, in a mixed ULA-cold dark matter model, N-body simulations by Schwabe et al. (2020) show that soliton cores are formed for ULA fractions higher than Ω_a/(Ω_m + Ω_a)≳0.1. However, below this fraction, these simulations show no evidence of the formation of soliton cores. Below ULA masses of m_a ≲ 10⁻²⁸ eV, ULAs are in the dark energy regime and do not contribute to the WL signal (see Sect. 2.4).

We conclude that the weak lensing mass calibration remains valid throughout the full ULA mass range as long as the ULA fraction remains below the aforementioned threshold. Hence, we use the profile described in Grandis et al. (2024).

4.5. Mixture model

The eRASS1 cluster catalog contains three classes of objects: galaxy clusters (C), which are of interest, and active galactic nuclei (AGN), as well as background fluctuations misclassified as clusters (NC), the latter two considered contaminants. Our model simultaneously accounts for the cluster counts and the contaminant fractions through the Poisson mixture model, as described in detail in Kluge et al. (2024). The total density is the sum of the three-component model:

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{tot}}(x|\mathrm{\Theta })={\lambda }_{\mathrm{C}}(x|\mathrm{\Theta })+{\lambda }_{\mathrm{AGN}}(x)+{\lambda }_{\mathrm{NC}}(x).\end{array}$$(23)

In terms of the total number of objects in each class, we obtained

$$\begin{array}{c}\hfill N_{\mathrm{tot}}(\mathrm{\Theta })=N_{\mathrm{C}}(\mathrm{\Theta })+f_{\mathrm{AGN}}{N}_{\mathrm{tot}}(\mathrm{\Theta })+f_{\mathrm{NC}}{N}_{\mathrm{tot}}(\mathrm{\Theta }),\end{array}$$(24)

where f_AGN and f_NC are the respective fractions of contaminants. Consequently, we can express the number of AGN, false detections, and the total number of objects in the catalog as a function of the cosmology-dependent number of clusters N_C:

$$\begin{array}{c}\hfill \{\begin{array}{c}{N}_{\mathrm{tot}}(\mathrm{\Theta })=(1/(1-f_{\mathrm{AGN}}-f_{\mathrm{NC}}))N_{\mathrm{C}}(\mathrm{\Theta })\hfill \\ N_{\mathrm{AGN}}(\mathrm{\Theta })=(f_{\mathrm{AGN}}/(1-f_{\mathrm{AGN}}-f_{\mathrm{C}}))N_{\mathrm{C}}(\mathrm{\Theta })\hfill \\ N_{\mathrm{NC}}(\mathrm{\Theta })=(f_{\mathrm{NC}}/(1-f_{\mathrm{AGN}}-f_{\mathrm{NC}}))N_{\mathrm{C}}(\mathrm{\Theta })\hfill \end{array}.\end{array}$$(25)

Starting from the formalism above, we can write the total number of objects as

$$\begin{array}{c}\hfill N_{\mathrm{tot}}(\mathrm{\Theta })=\frac{1}{1-f_{\mathrm{AGN}}-f_{\mathrm{NC}}}{\displaystyle {\int }_x}{\lambda }_{\mathrm{C}}(x|\mathrm{\Theta })\mathrm{d}x,\end{array}$$(26)

and that the number density of contaminants follows λ_AGN(x|θ) = N_AGN(Θ)𝒫_AGN(x), and λ_NC(x|θ) = N_NC(Θ)𝒫_NC(x). In this equation, 𝒫 describes the probability distribution function of the respective object, depending on the observables. Finally, the likelihood becomes

$$\begin{array}{c}\hfill \begin{array}{cc}ln\mathcal{L}(\mathrm{\Theta })=\hfill & {\displaystyle {\sum }_{i}ln({\lambda }_{\mathrm{C}}(x_i|\mathrm{\Theta })}\hfill \\ \hfill & +N_{\mathrm{AGN}}(\mathrm{\Theta }){\mathcal{P}}_{\mathrm{AGN}}(x_i)+N_{\mathrm{NC}}(\mathrm{\Theta }){\mathcal{P}}_{\mathrm{NC}}(x_i))\hfill \\ \hfill & {\displaystyle -\frac{1}{1-f_{\mathrm{AGN}}-f_{\mathrm{NC}}}{\displaystyle {\int }_x}{\lambda }_{\mathrm{C}}(x|\mathrm{\Theta })\mathrm{d}x.}\hfill \end{array}\end{array}$$(27)

The contaminant fractions f_AGN and f_NC are fit for during the parameter inference.

4.6. Integration of ultralight axions

Constraining the ULA fraction requires us to predict the massive dark matter halo number density in the corresponding model. We introduce ULAs in the cluster abundance cosmological pipeline by implementing a modified version of the power spectrum used to produce the HMF (Diehl & Weller 2021). This framework provides us with the proper transfer function and growth of structure models. In practice, we use axionCAMB (Grin et al. 2022), a modified version of CAMB that adds the evolution of primordial density fluctuations caused by an axion field to the Boltzmann solver. The effect of ULAs on the HMF is illustrated in Fig. 4 for different ULA masses m_a and ULA relic densities Ω_a. We use axionCAMB to modify the Hubble rate as a function of redshift H(z) for a universe with dark matter composed of ULAs. Using an updated Hubble rate ensures that we properly account for the redshift dependence of the Hubble rate on the effective properties of ULAs (dark energy or dark matter behavior in the Hubble rate) when computing the distances between background galaxies and the lens (or the observer) needed for the weak gravitational lensing calibration of the scaling laws a described in Sect. 4.4. Fig. 5 shows the influence of ULA mass on the Hubble rate in a cosmology with a high ULA abundance as dark matter. We also use axionCAMB to compute the matter power spectrum. We then compute the r.m.s. density fluctuations (Eq. (9)) and use them as input to the multiplicity function f(σ) from Tinker et al. (2008) to compute the halo number density.

Fig. 4. Halo mass function with Planck Collaboration VI (2020) cosmological values (ΩΛ ≈ 0.7, Ωm + Ωa ≈ 0.3) at z = 0.1 with a fixed relative ULA density of Ωa/(Ωcdm + Ωa) = 0.5 and varying ULA mass (left) and a fixed ULA mass of ma = 10−25 eV and varying ULA abundance (right) compared to the ΛCDM HMF (black). Low-mass halo suppression becomes stronger with decreasing ULA mass and increasing ULA abundance. For higher ULA masses of ma ≳ 10−24 eV, the HMF becomes indistinguishable from the one in a ΛCDM cosmology, as is the case for Ωa → 0.

Fig. 5. Hubble rate $H(z)=\dot{a}(z)/a(z)$ as a function of redshift. The dotted black line assumes a flat ΛCDM cosmology with Ωm = 0.6 and ΩΛ = 0.4. The dashed black line assumes a flat ΛCDM cosmology with Ωm = 0.3 and ΩΛ = 0.7. The colored lines represent Hubble rates for a ULA cosmology with fixed relative densities Ωa = 0.3, Ωm = 0.3, and ΩΛ = 0.4. The colored Hubble rates differ only in the ULA mass, as shown by the color bar. This plot demonstrates that ULAs undergo a domain change in the low-mass regime between log10(ma[eV]) = −32 and log10(ma[eV]) = −33. On the low-mass end, the Hubble rate approaches that of a flat LCDM cosmology but explains parts of the dark energy component. In the modifications to the pipeline used in G24, we use this accurate Hubble rate extracted from axionCAMB in the distance measure computations.

Consequently, in this work, we choose to model the effect of ULAs solely by implementing the power spectrum produced by the corresponding ULA model in the HMF obtained by Tinker et al. (2008). The impact of ULAs is thus obtained through the transfer function and the growth of structure used to compute the r.m.s. density fluctuations in Equation (9). However, accounting for ULAs also requires recalculating the critical overdensity for halo collapse, given that they do not form structures below the Jeans mass. Different studies in the literature have provided fitting functions to model this effect, with smoothed k window functions (Bohr et al. 2021), or analytical formulas (Dome et al. 2025). However, these methods have not been tested in the framework of cosmological parameter inference. Additionally, most of the correcting effects are significant at masses much below the mass range of clusters below the detection limits of the galaxy cluster catalogs of the eROSITA All-Sky Survey (∼10¹⁰ M_⊙ for the typical ranges of the corrections, compared with ∼10¹⁴ M_⊙ for a typical cluster detected in eRASS1). Consequently, it is not necessary to apply these corrections for the eRASS1 sample. We note that better modeling will be required to extend the analysis to a lower-mass range below M ≲ 10¹³ M_⊙.

Below an ULA mass of 10⁻²⁸ eV, we exclude ULAs from the matter density in the HMF formalism to account for the fact that their de Broglie wavelength exceeds the Jeans length and ULAs are not expected to contribute to halo formation in this regime. For ULA masses of log₁₀(m_a[eV]) ≲ −32.5, the de Broglie wavelength of the ULA BEC exceeds the size of the universe, making ULAs indistinguishable from a cosmological constant Λ. As such, a model is indistinguishable from ΛCDM, we exclude this regime from our analysis.

For ULA masses −32.5 ≲ log₁₀(m_a[eV]) ≲ −31.5, despite not contributing to halo formation, ULAs affect the Hubble rate, which in turn influences the distance measurements relevant for the weak lensing calibration. Above log₁₀(m_a[eV]) ∼ −32.5, only objects at redshifts higher than z = 0.8 (maximum redshift for galaxy clusters in the eRASS1 cosmology sample) are affected by the regime change in the Hubble rate. This involves the galaxies used for the gravitational weak lensing calibration of the scaling relations (see Sect. 4.4 for details) that extend up to z ≲ 6 in redshift space. We consider this effect by computing the distances based on the ULA Hubble rate extracted from axionCAMB in this regime. All ULA mass bins above the lowest-mass bin are not affected by the regime change in the Hubble rate, i.e., ULA masses above log₁₀(m_a[eV]) ≲ −31.5 (see Fig. 2 for further details). Fig. 5 illustrates how the regime affects the Hubble rate in a cosmology with a high ULA abundance as a function of ULA mass (color coded).

5. Results and discussion

In this section, we present the results of the Bayesian inference fitting of cosmological and ULA parameters through galaxy cluster number counts. We used the X-ray count rate, the optical richness, and (where available) the reduced tangential shear profiles given in the eRASS1 galaxy cluster cosmology catalog.

5.1. Consistency checks

We implement several changes to the base cosmological pipeline of G24 as described in detail in Sect. 4. To ensure the adaptations do not introduce inconsistencies, we reproduce the main cosmological results in the new ultralight axion framework by choosing fixed ULA parameters, which are indistinguishable from a ΛCDM cosmology. This limit can be achieved either by setting the ULA relic density to Ω_a = 0 or for values log₁₀(m_a[eV]) ≳ −18 (see Fig. 1 and Fig. 4). The code axionCAMB only allows positive values for Ω_a, so we chose Ω_a = 10⁻⁹, several orders of magnitude below any other cosmological fluid component in the relevant redshift range. We additionally fix the ULA mass to log₁₀(m_a[eV]) = −12. The matter power spectrum and hence σ(M, z) with this choice of ULA parameters is in perfect agreement with the ΛCDM quantities. We checked that running the adapted pipeline leads to almost identical results as compared to the ones obtained in the ΛCDM case shown in G24. A subset of the obtained cosmological parameters is presented in Fig. 6. We conclude that the applied changes did not introduce inconsistencies with the main cosmological results from eRASS1.

Fig. 6. Consistency check: Posteriors for the cosmological parameters Ωm, σ8, and S8 in the ΛCDM (ma = 10−12 eV and Ωa = 10−9). These values were chosen arbitrarily such that they satisfy ma ≫ 10−18 eV and Ωa ≪ 10−3. The results are consistent with G24. The inner and outer contours represent the 68% and 95% contour levels, respectively.

5.2. Cosmological constraints

Following the approach developed by Hložek et al. (2018) and Rogers et al. (2023), we perform a binned analysis in the logarithmic ULA mass log₁₀(m_a[eV]). We consider a total of nine bins, from log₁₀(m_a[eV]) = −24 in the highest-mass bin to log₁₀(m_a[eV]) = −32. Independent cosmological and scaling relation parameters are considered within each ULA mass bin. The analysis is not sensitive to how the ULA mass is treated within a mass bin, as drawing the ULA mass randomly within a bin or fixing it to the central bin value does not influence the posterior distribution.

We assumed a flat (ΛCDM+ ULA) cosmology with massless neutrinos to get a conservative estimate of the effect of ULAs. We leave the matter energy density fraction at redshift z = 0, Ω_m, the baryon energy density fraction at redshift z = 0, Ω_b, the scalar spectral index, n_s, the normalization of the matter power spectrum, log(A_s), and the Hubble constant, H₀, free. We leave the same parameters free for the scaling relations as done in G24. The priors on the cosmological parameters are shown in Table 1.

For all mass bins except the one with m = 10⁻²⁵ eV, we find no deviations from the ΛCDM cosmological posteriors obtained by G24, and the posteriors are in good agreement between the different ULA mass bins. We show the cosmological posterior distributions for the cosmological parameters Ω_m, σ₈, and S₈ for a subset of ULA mass bins covering the full mass range in Fig. 7.

Fig. 7. Posterior distributions for Ωm, σ8, and S8 in a subset of ULA mass bins. We find consistent cosmological constraints across all mass bins within 1σ (except for the excluded ma = 10−25 eV mass bin). The results are furthermore compatible with G24. The inner and outer contours represent the 68% and 95% contour levels, respectively.

In the mass bin around m_a = 10⁻²⁵ eV (i.e., covering the range 10^−25.5 eV ≤ m_a ≤ 10^−24.5 eV), we find significant discrepancies compared to G24 in the X-ray and optical scaling relation parameters (e.g., the discrepancy between A_X and G_X is 10.9σ, the discrepancy between A_λ and C_λ is 8.0σ). Figure 8 shows the posterior distributions of both the X-ray and optical scaling relations, comparing a representative selection of ULA mass bins with G24. The scaling relations, together with the gravitational weak lensing calibration, translate the observed count-rates and optical richnesses to the total masses of the clusters in the eRASS1 sample. Together with the contamination model and the selection function, these relations include all astrophysical modeling between observation and theory. As demonstrated in Sect. 4.4, ULAs do not alter the weak lensing signal compared to the ΛCDM scenario in the m_a = 10⁻²⁵ eV mass bin. Consequently, the astrophysical modeling should be independent from the abundance of ULAs with a mass around m_a = 10⁻²⁵ eV, i.e., the scaling relation parameter posteriors should be comparable with the ones reported by G24. However, in contrast to all other ULA mass bins, this is not the case for the m_a = 10⁻²⁵ eV bin (see Fig. 8). Thus, the cosmological parameters and the ULA relic density in this bin are not reliable, which is why we opt to exclude the m_a = 10⁻²⁵ eV mass bin from the analysis.

Fig. 8. Posterior distributions in a representative subset of ULA mass bins for the X-ray scaling relations (top) and the optical scaling relations (bottom). The large discrepancies in the astrophysical scaling relation model parameters between the ma = 10−25 eV mass bin on the one hand and G24 and the remaining mass bins on the other hand indicates an unphysical posterior distribution in the ma = 10−25 eV ULA mass bin. Consequently, the concerned mass bin is excluded from the analysis.

Fig. 9. Upper bounds on the relative ULA density Ωa in each logarithmic ULA mass bin. The left plot shows the upper bounds on Ωa on a linear scale, and the right plot shows the upper bounds on a logarithmic scale. The bar height indicates the 95% confidence interval. The orange bars indicate the bounds obtained by eRASS1 galaxy clusters only, and the gray hatched bars show the bounds obtained from Planck 2015 CMB data by Hložek et al. (2018). The red bars show the bounds obtained from combining the bounds from eRASS1 galaxy clusters with those from Planck 2015 CMB data.

Furthermore, we find that the posterior distribution is highly sensitive to the choice of priors, indicating a lack of statistical robustness in the m_a = 10⁻²⁵ eV ULA mass bin. However, we report no constraining power on the ULA relic density in this bin for prior choices comparable to those used by G24, similar to the results for the m_a = 10⁻²⁴ eV mass bin. Consequently, the m_a = 10⁻²⁵ eV mass bin does not contribute to constraining the ULA fraction. Possible reasons for the reported discrepancies include a breakdown of the HMF fitting function or the modeling of the shear signal for this specific ULA mass range. A detailed understanding of this effect will require further investigation beyond the scope of this work. In the figures with an interpolated curve between the bin values (Figs. 11, 13), we assume the Planck upper confidence level (Hložek et al. 2018) as a reasonable reference point for a bin with lacking constraining power for interpolation purposes only.

5.3. Constraints on the ULA parameter space

We found an exclusion region for ULAs by performing a binned analysis, providing upper bounds on the relic density within logarithmic mass bins. Since ULAs form BECs with characteristic extensions given by their thermal de Broglie wavelengths λ_dB ∼ 1/m_a, ULAs with different masses have imprints on different corresponding scales. Ultralight axions with masses of order 𝒪(10⁻²² eV) have de Broglie wavelengths of order 𝒪(kpc), while ULAs with masses of order 𝒪(10⁻³³ eV) have de Broglie wavelengths comparable to the size of the universe. As galaxy clusters are probing the 𝒪(Mpc) scales and above, their abundance can constrain ULAs with masses of order 𝒪(10⁻²⁶ eV) and below. The inverse relation between scale and ULA mass highlights the importance of galaxy groups in the sample to explore higher ULA mass regimes. Very light ULAs with masses of order 𝒪(10⁻³³ eV) mimic a dark energy component, as BECs on the scale of the universe add a global energy component to the universe. Ultralight axions with masses slightly above this threshold undergo a transition from an effective dark energy component to an effective dark matter component at some redshift where the size of the universe exceeds the de Broglie wavelength of the ULA condensate. For this regime, no models exist for halo collapse in the presence of ULAs. Hence, we are not able to probe this region and exclude ULA masses log₁₀(m_a[eV]) < −32.5 from our analysis.

We find upper bounds on the ULA relic density Ω_a obtained from binning the logarithmic ULA mass with a bin width of Δlog₁₀(m_a[eV]) = 1 around powers of ten. We treat each bin independently. We find upper bounds on the ULA relic density in each ULA mass bin with log₁₀(m_a[eV]) ≤ −26, as presented in Fig. 9. We confirm the exclusion region below the log₁₀(m_a[eV]) = −25 mass bin, which has been reported by Hložek et al. (2015, 2018), and Rogers et al. (2023). eRASS1 galaxy cluster abundance alone yields the tightest constraints on the ULA relic density in the mass bins log₁₀(m_a[eV]) ∈ { − 27, −26} with

$$\begin{array}{cc}\hfill {\mathrm{\Omega }}_{\mathrm{a}}(-26.5<{log}_{10}(m_{\mathrm{a}}[\mathrm{eV}])\le -25.5)& <0.0079,\hfill \\ \hfill {\mathrm{\Omega }}_{\mathrm{a}}(-27.5<{log}_{10}(m_{\mathrm{a}}[\mathrm{eV}])\le -26.5)& <0.0035,\hfill \end{array}$$(28)

both at a 95% confidence level. Fig. 10 shows that these are the tightest upper bounds that have been presented in the literature. Combining our results with the available chains from Planck 2015 CMB data (Hložek et al. 2018), we find even tighter constraints in both mass bins:

$$\begin{array}{cc}\hfill {\mathrm{\Omega }}_{\mathrm{a}}(-26.5<{log}_{10}(m_{\mathrm{a}}[\mathrm{eV}])\le -25.5)& <0.0058,\hfill \\ \hfill {\mathrm{\Omega }}_{\mathrm{a}}(-27.5<{log}_{10}(m_{\mathrm{a}}[\mathrm{eV}])\le -26.5)& <0.0030.\hfill \end{array}$$(29)

Fig. 10. Comparison of ULA constraints obtained from eRASS1 galaxy cluster abundance (this work), Planck 2015 (Hložek et al. 2018), and 2018 CMB data; BOSS galaxy clustering data (both by Rogers et al. 2023); and combined analyses of eRASS1 and Planck 2015 (this work) and Planck 2018 and BOSS (Rogers et al. 2023). The bar height indicates the 95% exclusion limits of the ULA density Ωa times h2 = (H0/100)2 in each log10(ma[eV]) bin. eRASS1 yields the tightest upper bounds on dark matter ultralight axions in the log10(ma[eV]) ∈ { − 27, −26} bins.

From eRASS1 galaxy cluster number counts alone, we find that the ULA relic density Ω_a cannot exceed 2.7% of the total energy density of the universe for m_a ≤ 10^−27.5 eV and is bound below 0.84% for ULA masses 10^−27.5 eV ≤ m_a ≤ 10^−26.5 eV and below 0.36% for ULA masses 10^−26.5 eV ≤ m_a ≤ 10^−25.5 eV. A significant contribution of ULAs to the total dark matter density can thus be ruled out by the growth of structure as observed through eRASS1 galaxy cluster number counts.

It is particularly interesting to compare the eRASS1 cluster number count constraints with those found using BOSS galaxy clustering data as another large-scale structure probe. Rogers et al. (2023) studied the cosmological implications of ULAs using the twelfth data release of the Baryon Oscillation Spectroscopic Survey (BOSS) catalog (Dawson et al. 2013; Alam et al. 2017) and found upper bounds on the ULA fraction in a binned analysis in the range −32 ≤ log₁₀(m_a[eV]) ≤ −26. As expected, galaxy clustering and galaxy cluster number counts constrain different ULA mass regimes. This is due to the relation between the ULA mass and the scale it affects: The heavier the ULA, the larger the affected scale in the matter power spectrum, for instance, a ULA with a mass of ∼10⁻²⁴ eV is only observable on scales below ∼200 kpc, while a ULA with a mass of ∼10⁻²⁵ eV is already observable on scales ∼500 kpc and below (see Fig. 1). While, in comparison to other probes, Rogers et al. (2023) find the tightest constraints (Ω_ah² < 0.0418) in the ULA mass bin around m_a = 10⁻²⁵ eV, corresponding to an effect on smaller scales, eRASS1 cluster number counts yield the tightest constraints among complementary in the mass bins m_a ∈ {10⁻²⁷ eV, 10⁻²⁶ eV}. Ultralight axions in this mass range affect larger scales, which is constrainable using galaxy cluster number counts. Similar to galaxy clustering data, eRASS1 cluster number counts have reduced constraining power at scales extending the typical probe scales, as seen in the lowest-mass bins. On these largest scales (corresponding to ULA masses m_a ≲ 10⁻²⁸ eV), CMB measurements by the Planck collaboration (Planck Collaboration XIII 2016; Planck Collaboration VI 2020) have the highest constraining power and yield the tightest constraints (Hložek et al. 2015, 2018; Rogers et al. 2023). For better comparison with the results obtained using other probes (see Fig. 10), we also provide the bounds of each constraining bin for Ω_a and Ω_ah² in Table 2.

We emphasize that the weak lensing mass calibration as described in Sect. 4.4 remains valid over the full considered ULA mass range as the 95% confidence levels do not exceed ULA fractions of Ω_a/(Ω_m + Ω_a)≈0.1 for m_a ≲ 10⁻²⁶ eV, and below m_a ≈ 10^−28.5 eV, ULAs are not included in the HMF model and thus the modeling used in G24 remains valid. For higher ULA masses, Sect. 4.4 justifies the modeling choice of keeping NFW profiles for the weak lensing mass calibration up to any ULA fraction.

Our current pipeline is limited by the model of the halo abundance in ULA cosmology, where we use a modified power spectrum in the Tinker et al. (2008) fitting function. Although this is justified by numerical simulations (see Sect. 4.6), our results will be improved if suitable emulators are developed for this model. Additionally, for this pioneering work, we fix the masses of the neutrinos, while the effect of neutrinos is degenerate with ULAs. This modeling choice results in our constraints being conservative upper limits. Including massive neutrinos in the next analyses significantly tightened our constraints.

5.4. Forecasts for the deepest eROSITA survey

This section presents an estimation quantifying the potential of future galaxy cluster samples, illustrated by the deepest eROSITA All-Sky Surveys. The constraints presented in this work are based on eRASS1, which represents ∼22% of the stacked consequent eROSITA All-Sky Survey (eRASS:5, hereafter). In the near future, the release of deeper eROSITA data and associated catalogs will significantly increase the number of available galaxy clusters and groups. In particular, we expect the number of low-mass galaxy groups to increase significantly with higher exposure time. To assess the expected evolution of cluster abundance in future eROSITA data releases, we generate a mock sample representing a catalog for eRASS:5. This mock catalog is constructed using the methodology outlined in G24 to validate the pipeline.

We generated a dark matter halo catalog with the cosmology inferred in G24, using the halo mass function from Tinker et al. (2008). We then create mock X-ray and optical observables following the scaling relations described in Section 4.1. We then select the fraction of objects that would be detected using the selection function described in Clerc et al. (2024). The X-ray selection function depends on X-ray count rates, hydrogen column density, background, and exposure time. We scale the exposure time by a factor of 4.6 to simulate a survey with 4.6 times the depth of eRASS1. The result is a catalog containing ∼15 000 clusters with 0.1 < z < 0.8. Figure 12 presents the cluster abundance in eRASS1 and in the mock eRASS:5 catalog. In the deeper eROSITA survey, the statistical power of cosmological constraints will be driven by the detection of a large and significant sample of galaxy groups. At 10^13.8 M_⊙, we detect approximately five times more massive clusters than in eRASS1, where the constraints on the parameter space are mostly constrained by the number of lower-mass galaxy clusters and galaxy groups (see Fig. 4).

Fig. 11. Exclusion regions (95% confidence level) in the log10(ma[eV]) − Ωah2 space as constrained by different probes. Colored regions represent the part of the parameter space where one ULA species can be excluded at 95% confidence level, whereas white regions represent the unconstrained parameter space. eRASS1 galaxy cluster number counts (orange) exclude a new region in the parameter space between 10−27.5 eV ≲ ma ≲ 10−26 eV, corresponding to galaxy cluster scales. Combining eRASS1 galaxy cluster number counts with the Planck 2015 CMB data (red) further improves the eRASS1 bounds on the ULA relic density. For ULA masses below ma ≲ 10−28 eV, Planck 2018 CMB data yields the tightest bounds. Above ma ≳ 10−25 eV, galaxy clustering data from BOSS yields the tightest constraints. The curves are interpolated between the central bin values, which represent the 95% confidence levels in the respective bins. In the ma = 10−25 eV mass bin, the Planck upper confidence level has been assumed as a reasonable choice for the interpolation (for details on that bin, see Sect. 5.2).

Fig. 12. Forecasts for an eRASS:5 survey. Left: Number count of X-ray detected clusters for an eRASS:5 survey (red) compared with an eRASS1 catalog (black). The Poisson error bars are represented by filled lines. Right: Relative difference of eRASS:5 detected clusters and eRASS1 detected objects. eRASS:5 will contain five times more clusters at M = 1013.8 M⊙.

We assumed that the Poisson error of the cluster abundance is the main statistical limiting factor. This means that we assume that the additional low-mass objects possess X-ray, optical, and weak lensing quantities measured as precisely as the ones obtained for massive clusters. We provide an optimistic and a pessimistic scenario based on the potential inclusion of low-mass clusters in Fig. 13. In the optimistic case, our lower limit decreases by a factor of 2.2, while it is lowered by a factor of 1.7 in the pessimistic case. The reliability of our forecasts for the eRASS:5 data is influenced by the accuracy of the selection function; the robustness of the model used to predict the abundance, including the HMF model; the universality of the scaling relations over the full survey mass range; and the precision of forthcoming weak lensing measurements for the mass calibration of low-mass galaxy groups.

Fig. 13. Forecast for eRASS:5 data. The exclusion region of the ULA parameter space obtained from eRASS:5 mock catalog is presented with the orange dotted line for a pessimistic scenario and with the orange dashed line for an optimistic scenario, and can be compared with the plain orange region corresponding to the exclusion region obtained by eRASS1. For the missing bin, the Planck upper confidence level was assumed for the interpolation (see Sect. 5.2).

6. Conclusions

The unique properties of ULAs lead to observable phenomena such as solitonic cores in dark matter halos and interference patterns on cosmological scales. These features can be probed through astrophysical observations, providing potential evidence of their existence.

Galaxy cluster number counts inferred from the eRASS1 cluster cosmology catalog (Bulbul et al. 2024; Kluge et al. 2024) can be used to test and constrain the viable dark matter candidates, such as ULAs, through the evolution of the structure growth. Indeed, eRASS1, being the largest ICM-selected sample to date, has a significant constraining power. Additionally, it has the potential to probe the group regime, which is very sensitive to dark matter signatures.

In this work, we performed Bayesian inference on cosmological parameters in a universe with an additional ULA component parameterized by the logarithmic ULA mass and ULA relic density {log₁₀(m_a[eV]), Ω_a}. The matter power spectrum was computed by axionCAMB (Hložek et al. 2015), from which we computed the HMF using the multiplicity function by Tinker et al. (2008). The forward modeling process requires the selection function developed by Clerc et al. (2024) and additional parameters from the mixture model that take care of possible contaminants, the count-rate-mass scaling relation, and the richness-mass scaling relation, which are calibrated by gravitational weak lensing shear maps (Grandis et al. 2024). Apart from the adaptations described in Sect. 4.6, we used the framework established in G24 for the parameter inference.

We performed a binned analysis with logarithmic ULA mass bins covering the whole mass spectrum of ultralight axions from the dark matter regime starting at log₁₀(m_a[eV]) = −24 down to the dark energy regime with the lowest ULA mass at log₁₀(m_a[eV]) = −32 including an independent full cosmological inference in each bin while assuming a flat universe. We constrained the following set of cosmological parameters: the matter relic density, Ω_m, and the normalization of the matter power spectrum, log(A_s). We left the Hubble constant, H₀; the scalar spectral index, n_s; and the baryon relic density, Ω_b, free but could constrain them as done in G24. We chose a uniform prior on Ω_a between 0 and 1 in both scenarios. All common parameters have the same priors as in G24.

We performed a consistency check with the eRASS1 cosmological results from G24 by fixing ULA parameter values to log₁₀(m_a[eV]) ≫ −18 and Ω_a → 0, effectively switching off the effect of ULAs on structure growth. The results are consistent with G24, indicating that including ULAs does preserve the correctness and reliability of the inference framework.

The binned analysis of eRASS1 galaxy cluster number counts alone yields tight bounds on the ULA relic density, especially in the two highest ULA mass bins within the exclusion region log₁₀(m_a[eV]) ∈ { − 27, −26}. This ULA mass range corresponds to physical scales (de Broglie wavelengths) of ULA BECs comparable to the typical galaxy cluster scales around 1 Mpc to 10 Mpc. The parameter space excluded by eRASS1 galaxy clusters in comparison to complementary probes is shown in Fig. 11. In the mass range of 10^−27.5 eV ≤ m_a ≤ 10^−26.5 eV, the upper limit on the ULA relic density is Ω_a < 0.0036 of the total energy density of the universe. Similarly, in the mass range of 10^−26.5 eV ≤ m_a ≤ 10^−25.5 eV, the ULA relic density is constrained to Ω_a < 0.0084 of the total energy density of the universe. The m_a = 10⁻²⁵ eV mass bin was excluded from the analysis due to a lack of modeling and statistical robustness in that particular bin. In the mass regime 10^−32.5 eV ≤ m_a ≤ 10^−27.5 eV, the ULA relic density is bound below Ω_a ≲ 0.055. We conclude that the evolution of structure growth as observed by eRASS1 galaxy cluster number counts rules out ULAs as a significant contribution to the total dark matter density of the Universe. The constraints on the abundance at different confidence levels concerning a logarithmic ULA mass bin are presented in Fig. 9. Overall, these results highlight the potential of cluster abundance measurements and the evolution of the growth of structure as a probe of ULAs. For the first time, we have obtained constraints using cluster abundance in a redshift and scale regime complementary to the one probed by other cosmological probes.

In the near future, statistically more powerful cosmological simulations, including ULAs (or ultralight scalar fields, fuzzy dark matter) that predict an HMF model and shear profiles for a sufficiently large cosmological parameter space, will be required to provide unbiased constraints. We note that the constraints presented in this work were obtained using the eRASS1 cosmology sample of 5 259 X-ray-selected massive clusters. The primary sample contains a higher fraction of low-mass objects, so we are more sensitive to ULA signatures, provided that we carefully model the baryonic systematic effects that affect groups. Future analyses of deeper surveys with eROSITA at eRASS:5 depth will provide larger samples of low-mass halos with a greater constraining power in the ULA parameter space, improving the overall constraints by a factor of 1.7.

Acknowledgments

The authors thank the referee for the helpful and constructive comments on the draft. We further thank Christian Garrel for his contributions to the codebase and the valuable exchange with him. Thanks to Johannes Diehl, Julia Sisk-Reynes, and Luis A. Ureña-López for the helpful discussions. This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS software system developed by the German eROSITA consortium. S. Zelmer, E. Bulbul, V. Ghirardini, A. Liu, and X. Zhang acknowledge financial support from the European Research Council (ERC) Consolidator Grant under the European Union’s Horizon 2020 research and innovation program (grant agreement CoG DarkQuest No 101002585). N. Clerc was financially supported by CNES. T. Schrabback and F. Kleinebreil acknowledge support from the German Federal Ministry for Economic Affairs and Energy (BMWi) provided through DLR under projects 50OR2002, 50OR2106, and 50OR2302, as well as the support provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 415537506. M. Brüggen acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. S. Krippendorf’s work has been partially supported by STFC consolidated grants ST/T000694/1 and ST/X000664/1.

References

All Tables

All Figures

Fig. 1. Effect of ULAs on the matter power spectrum recovered assuming the Planck Collaboration VI (2020) cosmological values (ΩΛ ≈ 0.7, Ωm + Ωa ≈ 0.3) at z = 0.1. On the left, the change in the matter power spectrum with varying ULA mass is shown when the relative ULA density is fixed to Ωa/(Ωcdm + Ωa) = 0.5. The panel on the right displays the dependence of the power spectrum on varying ULA abundance at a fixed ULA mass of ma = 10−25 eV. The concordance ΛCDM model is shown in a black curve on both panels. Small-scale suppression becomes stronger with decreasing ULA mass and increasing ULA abundance. For higher ULA masses of ma ≫ 10−22 eV, the matter power spectrum becomes indistinguishable from the one in a ΛCDM cosmology, as is the case for Ωa → 0.

In the text

Fig. 2. Demonstration of the ULA mass and redshift dependent regime change in the Hubble rate. In the pink regime, ULAs contribute to the Hubble rate e in a manner similar to dark matter; in the blue region, ULAs contribute to the Hubble rate in a dark energy-like fashion. The orange shaded region shows the eRASS1 cluster redshift range, while the gray shaded region shows the redshift range of lensed galaxies from the DES, HSC, and KiDS data, described in Sect. 3.2. Above a ULA mass of ma ∼ 10−31.5 eV, ULAs should be modeled as dark matter in the Hubble rate, resulting in correct distance measurements during the parameter inference. For ULAs with masses below ma ∼ 10−31.5 eV, the distance measurements are affected by the regime change via the Hubble rate, where ULAs are treated as dark energy.

In the text

Fig. 3. Density profiles (left column) and projected density contrasts (right column) of halos of different mass (color coded) at redshift z = 1 for different ULA particle masses (rows) in the extreme case of a ULA-only dark matter sector. Full lines show the total profiles, dashed lines show the corresponding cold dark matter profiles (assumed to follow a Navarro-Frenk-White profile), and dotted lines show the soliton core profile. In the right column, we highlight the scales we use for the weak lensing measurement in white.

In the text

Fig. 4. Halo mass function with Planck Collaboration VI (2020) cosmological values (ΩΛ ≈ 0.7, Ωm + Ωa ≈ 0.3) at z = 0.1 with a fixed relative ULA density of Ωa/(Ωcdm + Ωa) = 0.5 and varying ULA mass (left) and a fixed ULA mass of ma = 10−25 eV and varying ULA abundance (right) compared to the ΛCDM HMF (black). Low-mass halo suppression becomes stronger with decreasing ULA mass and increasing ULA abundance. For higher ULA masses of ma ≳ 10−24 eV, the HMF becomes indistinguishable from the one in a ΛCDM cosmology, as is the case for Ωa → 0.

In the text

Fig. 5. Hubble rate $H(z)=\dot{a}(z)/a(z)$ as a function of redshift. The dotted black line assumes a flat ΛCDM cosmology with Ωm = 0.6 and ΩΛ = 0.4. The dashed black line assumes a flat ΛCDM cosmology with Ωm = 0.3 and ΩΛ = 0.7. The colored lines represent Hubble rates for a ULA cosmology with fixed relative densities Ωa = 0.3, Ωm = 0.3, and ΩΛ = 0.4. The colored Hubble rates differ only in the ULA mass, as shown by the color bar. This plot demonstrates that ULAs undergo a domain change in the low-mass regime between log10(ma[eV]) = −32 and log10(ma[eV]) = −33. On the low-mass end, the Hubble rate approaches that of a flat LCDM cosmology but explains parts of the dark energy component. In the modifications to the pipeline used in G24, we use this accurate Hubble rate extracted from axionCAMB in the distance measure computations.

In the text

	Fig. 6. Consistency check: Posteriors for the cosmological parameters Ωm, σ8, and S8 in the ΛCDM (ma = 10−12 eV and Ωa = 10−9). These values were chosen arbitrarily such that they satisfy ma ≫ 10−18 eV and Ωa ≪ 10−3. The results are consistent with G24. The inner and outer contours represent the 68% and 95% contour levels, respectively.
In the text

	Fig. 7. Posterior distributions for Ωm, σ8, and S8 in a subset of ULA mass bins. We find consistent cosmological constraints across all mass bins within 1σ (except for the excluded ma = 10−25 eV mass bin). The results are furthermore compatible with G24. The inner and outer contours represent the 68% and 95% contour levels, respectively.
In the text

Fig. 8. Posterior distributions in a representative subset of ULA mass bins for the X-ray scaling relations (top) and the optical scaling relations (bottom). The large discrepancies in the astrophysical scaling relation model parameters between the ma = 10−25 eV mass bin on the one hand and G24 and the remaining mass bins on the other hand indicates an unphysical posterior distribution in the ma = 10−25 eV ULA mass bin. Consequently, the concerned mass bin is excluded from the analysis.

In the text

Fig. 9. Upper bounds on the relative ULA density Ωa in each logarithmic ULA mass bin. The left plot shows the upper bounds on Ωa on a linear scale, and the right plot shows the upper bounds on a logarithmic scale. The bar height indicates the 95% confidence interval. The orange bars indicate the bounds obtained by eRASS1 galaxy clusters only, and the gray hatched bars show the bounds obtained from Planck 2015 CMB data by Hložek et al. (2018). The red bars show the bounds obtained from combining the bounds from eRASS1 galaxy clusters with those from Planck 2015 CMB data.

In the text

Fig. 10. Comparison of ULA constraints obtained from eRASS1 galaxy cluster abundance (this work), Planck 2015 (Hložek et al. 2018), and 2018 CMB data; BOSS galaxy clustering data (both by Rogers et al. 2023); and combined analyses of eRASS1 and Planck 2015 (this work) and Planck 2018 and BOSS (Rogers et al. 2023). The bar height indicates the 95% exclusion limits of the ULA density Ωa times h2 = (H0/100)2 in each log10(ma[eV]) bin. eRASS1 yields the tightest upper bounds on dark matter ultralight axions in the log10(ma[eV]) ∈ { − 27, −26} bins.

In the text

Fig. 11. Exclusion regions (95% confidence level) in the log10(ma[eV]) − Ωah2 space as constrained by different probes. Colored regions represent the part of the parameter space where one ULA species can be excluded at 95% confidence level, whereas white regions represent the unconstrained parameter space. eRASS1 galaxy cluster number counts (orange) exclude a new region in the parameter space between 10−27.5 eV ≲ ma ≲ 10−26 eV, corresponding to galaxy cluster scales. Combining eRASS1 galaxy cluster number counts with the Planck 2015 CMB data (red) further improves the eRASS1 bounds on the ULA relic density. For ULA masses below ma ≲ 10−28 eV, Planck 2018 CMB data yields the tightest bounds. Above ma ≳ 10−25 eV, galaxy clustering data from BOSS yields the tightest constraints. The curves are interpolated between the central bin values, which represent the 95% confidence levels in the respective bins. In the ma = 10−25 eV mass bin, the Planck upper confidence level has been assumed as a reasonable choice for the interpolation (for details on that bin, see Sect. 5.2).

In the text

	Fig. 12. Forecasts for an eRASS:5 survey. Left: Number count of X-ray detected clusters for an eRASS:5 survey (red) compared with an eRASS1 catalog (black). The Poisson error bars are represented by filled lines. Right: Relative difference of eRASS:5 detected clusters and eRASS1 detected objects. eRASS:5 will contain five times more clusters at M = 1013.8 M⊙.
In the text

Fig. 13. Forecast for eRASS:5 data. The exclusion region of the ULA parameter space obtained from eRASS:5 mock catalog is presented with the orange dotted line for a pessimistic scenario and with the orange dashed line for an optimistic scenario, and can be compared with the plain orange region corresponding to the exclusion region obtained by eRASS1. For the missing bin, the Planck upper confidence level was assumed for the interpolation (see Sect. 5.2).

In the text