We study the implications of relaxing the requirement for ultralight axions to account for all dark matter in the Universe by examining mixed dark matter MDM cosmologies with axion fractions f leq 0.3 within the fuzzy dark matter FDM window 1025 eV lesssim m lesssim 1023 eV. Our simulations using a new MDM gravity solver implemented in AxiREPO capture wave dynamics across various scales with high accuracy down to redshifts zapprox 1. We identify halos with Rockstar using the CDM component and find good agreement of inferred halo mass functions HMFs and concentrationmass relations with theoretical models across redshifts z110. This justifies our halo finder approach a posteriori as well as the assumptions underlying the MDM halo model AxionHMcode. Using the inferred axion halo mass  cold halo mass relation MtextaMtextc and calibrating a generalised smoothing parameter alpha to our MDM simulations we present a new version of AxionHMcode. The code exhibits excellent agreement with simulations on scales klt 20  h cMpc1 at redshifts z13.5 for fleq 0.1 around the fiducial axion mass m  1024.5 eV   3.16times 1025 eV with maximum deviations remaining below 10. For axion fractions fleq 0.3 the model maintains accuracy with deviations under 20 at redshifts zapprox 1 and scales klt 10  h cMpc1 though deviations can reach up to 30 for higher redshifts when f0.3. Reducing the runtime for a single evaluation of AxionHMcode to below 1 minute these results highlight the potential of AxionHMcode to provide a robust framework for parameter sampling across MDM cosmologies in Bayesian constraint and forecast analyses.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Improved Halo Model Calibrations for Mixed Dark Matter Models of Ultralight Axions
{'Improved Halo Model Calibrations for Mixed Dark Matter Models of Ultralight Axions': "Tibor Dome 1 , 2 ⋆ , Simon May 3 , 4 , Alex Laguë 5 , David J. E. Marsh 7 , Sarah Johnston 6 , Sownak Bose 6 , Alex Tocher 1 , 2 , Anastasia Fialkov 1 , 2  \n1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK  \n- 2 Kavli Institute for Cosmology, Madingley Road, Cambridge, CB3 0HA, UK\n- 3 Perimeter Institute for Theoretical Physics, 31 Caroline Street North, Waterloo, ON, N2L 2Y5, Canada\n- 4 Department of Physics, North Carolina State University, Raleigh, NC, 27695-8202, USA\n- 5 Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street, Philadelphia, PA, USA 19104\n- 6 Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK\n- 7 King's College London, Strand, London, WC2R 2LS, United Kingdom  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': 'We study the implications of relaxing the requirement for ultralight axions to account for all dark matter in the Universe by examining mixed dark matter (MDM) cosmologies with axion fractions f ≤ 0 . 3 within the fuzzy dark matter (FDM) window 10 -25 eV ≲ m ≲ 10 -23 eV. Our simulations, using a new MDM gravity solver implemented in A xi REPO, capture wave dynamics across various scales with high accuracy down to redshifts z ≈ 1. We identify halos with R ockstar using the CDM component and find good agreement of inferred halo mass functions (HMFs) and concentration-mass relations with theoretical models across redshifts z = 1 -10. This justifies our halo finder approach a posteriori as well as the assumptions underlying the MDM halo model A xion HM code . Using the inferred axion halo mass-cold halo mass relation M a( M c) and calibrating a generalised smoothing parameter α to our MDM simulations, we present a new version of A xion HM code . The code exhibits excellent agreement with simulations on scales k < 20 h cMpc -1 at redshifts z = 1 -3 . 5 for f ≤ 0 . 1 around the fiducial axion mass m = 10 -24 . 5 eV = 3 . 16 × 10 -25 eV, with maximum deviations remaining below 10 %. For axion fractions f ≤ 0 . 3, the model maintains accuracy with deviations under 20 % at redshifts z ≈ 1 and scales k < 10 h cMpc -1 , though deviations can reach up to 30 % for higher redshifts when f = 0 . 3. Reducing the run-time for a single evaluation of A xion HM code to below 1 minute, these results highlight the potential of A xion HM code to provide a robust framework for parameter sampling across MDMcosmologies in Bayesian constraint and forecast analyses.  \nKey words: dark matter - cosmology: large-scale structure of Universe - methods: numerical', '1 INTRODUCTION': "Ideal candidates for ultralight bosonic dark matter (DM) - also called fuzzy dark matter (FDM) - from particle physics are ultralight axions which can emerge from both field theory and string theory. In field theory, these candidates are typically axion-like particles (ALPs), which share their conceptual origins with quantum chromodynamics (QCD) axions (Preskill et al. 1983; Abbott & Sikivie 1983; Dine & Fischler 1983) but can have much lighter masses with couplings to the Standard Model typically taken as free parameters (although see Dias et al. 2014; Kim & Marsh 2016). In string theory, pseudoscalar fields with axion-like properties arise in compactifications as Kaluza-Klein (KK) zero modes of antisymmetric tensor fields, with their masses and couplings determined by the topology  \nand geometry of the compact manifold (Svrcek & Witten 2006; Arvanitaki et al. 2010; Visinelli & Vagnozzi 2019; Gendler et al. 2023).  \nIn both field theory and string theory, the lightest of these axions could have a mass lower than 10 -19 eV, which implies the existence of particles with de Broglie wavelengths on galactic or even cosmological scales. In the last decade, several extensions of the simplest one-field FDM model have been explored in the literature. Some notable studies include vector FDM where FDM is a higher-spin field (Amin et al. 2022), multi-field FDM (Gosenca et al. 2023; Luu et al. 2024), and self-interacting FDM (Mocz et al. 2023) where the quartic self-coupling term is included. In the self-interacting model, soliton cores can undergo a phase transition from dilute to denser states, and small-scale structure formation may be enhanced.  \nString theory compactifications typically yield a plenitude of N = O (100) axion fields with logarithmically distributed masses (Mehta et al. 2021). Importantly, the relic density of ultralight parti-  \ncles, Ω FDM, does not necessarily match the total DM density of the Universe, Ω m -Ω b. For example, Bachlechner et al. (2018) show that in theories with N = O (100) axions and a lightest axion of mass m ∼ 10 -22 eV, axion misalignment can generate a DM abundance of Ω FDM ∼ 0 . 1. Models with a small number N = O (1) of axions can also give rise to Ω FDM ∼ 0 . 1 via axion misalignment (Cicoli et al. 2022) if the product S f a of the instanton action S (giving rise to the axion potential) and the axion decay constant f a satisfies S f a ≳ 1, slightly violating the Weak Gravity Conjecture (Hebecker et al. 2018; Alonso & Urbano 2017). This makes phenomenological multi-field FDM models a natural extension within the framework of string theory.  \nThe phenomenology of the two-field FDM model can be complex, novel, and compelling, as shown by several studies. Huang et al. (2023) conducted the first two-field cosmological simulations with cold dark matter (CDM) initial conditions and showed that minor admixtures of heavier axions to a lighter-mass major component can neither significantly disturb stable major-component solitons nor form minor-component solitons. Gosenca et al. (2023) studied stellar heating in FDM halos with two or more fields, while Glennon et al. (2023) simulated soliton collision of two-field FDM with inter-field and self-interactions. Jain et al. (2024) investigated Bose-Einstein condensation (BEC) formation in the kinetic regime, and Luu et al. (2024) showed that in the regime where the mass ratio m 2 / m 1 is in the range 2 -7, a minimal model with uncoupled fields predicts a central soliton with a nested structure, distinguishable from the typical flat-core soliton found in one-field halos. Importantly, multi-field FDM circumvents the Catch 22 problem that challenges the use of the simplest non-interacting one-field FDM as a solution to the small scale problems of CDM (Marsh & Silk 2014).  \nIn the context of mixed dark matter (MDM) models, we make the simplifying assumption that only one of the DM particle species is ultralight while the others have negligible de Broglie wavelengths. In this case, we group the combined relic density of the (near-)collisionless components into Ω CDM since the species can be modelled as CDM by virtue of the correspondence between the Schrödinger-Poisson and Vlasov-Poisson equations (Mocz et al. 2018). How do observations constrain MDM models? While the FDM constraints compiled in Dome et al. (2022) assume a pure FDM cosmology, many constraints were recently reported assuming axions do not constitute all of the DM in the Universe. In the mass range 10 -32 eV ≲ m ≲ 10 -26 eV, FDM can only comprise a few percent of the total DM (Hlozek et al. 2015; Hložek et al. 2018; Laguë et al. 2022; Rogers et al. 2023). For higher values of m , the bound becomes weaker since the FDM power suppression moves to smaller scales, while for smaller values of m , FDM behaves as dark energy and is strongly degenerate with ΩΛ .  \nAt the higher-mass end, Lyα forest and ultra-faint dwarf (UFD) data indicates that DM cannot be fully described by pure FDM models in the mass range 10 -21 eV ≲ m ≲ 3 × 10 -19 eV (Armengaud et al. 2017; Marsh & Niemeyer 2019; Zimmermann et al. 2021; Dalal & Kravtsov 2022; Iršiˇc et al. 2024), with a stronger bound for 10 -23 eV ≲ m ≲ 10 -21 eV where FDM must not comprise more than O (10 %) of the total DM (Kobayashi et al. 2017). The Dark Energy Survey year 1 was used by Dentler et al. (2022) to search for shear-correlation suppressions caused by FDM (assuming pure FDM), ruling out the existence of pure FDM in the mass range 10 -25 eV ≲ m ≲ 10 -23 eV (also see Preston et al. 2024). However, current constraints suggest that it is possible for FDM to exist in fairly large portions (albeit not the entirety of the DM) in the region 10 -25 eV ≲ m ≲ 10 -23 eV, referred to as the FDM window .  \nRecent studies have demonstrated the potential of 21 cm inter-  \nferometers, particularly the Hydrogen Epoch of Reionisation Array (HERA), to detect FDM fractions of f ≳ O (1 %) within the FDM window and beyond (Jones et al. 2021; Flitter & Kovetz 2022). This potential was further supported by Lazare et al. (2024), who in addition established an upper bound on FDM with a particle mass m = 10 -23 eV, restricting it to 16 % of the total DM. Their constraints are based on a range of observations, including UV luminosity function (LF) data from the Hubble Space Telescope, constraints on the neutral hydrogen fraction from high-redshift quasar spectroscopy, cosmic microwave background (CMB) optical depth measurements from Planck, and upper bounds on the 21 cm power spectrum from HERA. These constraints tighten for smaller masses, reaching down to 1 % for m = 10 -26 eV. However, their machine learning-based emulator does neither account for the 'quantum pressure' term in the evolution of the axion field, nor the impact of FDM on star formation, and 21 cm FAST relies on the excursion set formalism (Mesinger et al. 2011), which may not adequately capture small-scale physics. Using UVLF and Planck CMB data only, Winch et al. (2024) report slightly weaker constraints, restricting FDM to less than ≈ 22 % of the DM at m = 10 -23 eV and to less than ≈ 5 % at m = 10 -26 eV.  \nThe future Square Kilometre Array (SKA) is expected to provide strong constraints on FDM and MDM models, at few per cent level for masses m ∼ 10 -22 eV, and improving to below 1 % for lighter masses (Bauer et al. 2021; Dome et al. 2024). For pure FDM, Hotinli et al. (2022) showed that by means of velocity acoustic oscillations in the large-scale 21 cm power spectrum, HERA may be sensitive to axions with masses up to m ≈ 10 -18 eV (see also Marsh 2015, for FDM relative velocity e ff ects).  \nShevchuk et al. (2023) argued that FDM with a particle mass of m ≲ 10 -24 eV is inconsistent with the observed Einstein radii of several strong lensing systems. However, their analysis is based on simplified models, such as expressing total DM density profiles as ρ DM( r ) = f ρ FDM( r ) + (1 -f ) ρ CDM( r ), and assuming idealised solitons derived solely from FDM simulations, without accounting for contributions from cold dark matter or baryons.  \nWe thus believe that the FDM window remains an appealing regime to explore in MDM simulations, for a number of reasons. First, Blum & Teodori (2021) proposed that O (10 %) of DM in the form of FDM with particle mass in the window, m ∼ 10 -25 eV, could explain the tension between inferences of H 0, the current expansion rate of the Universe, which are based on the time delay in lensed quasar measurements, and those based on CMB observations. Note that both Shevchuk et al. (2023) and Blum & Teodori (2021) assume a soliton core in the (lens) density profile for FDM fractions of O (10 %), the validity of which we will address in more detail later (see Secs. 4.2 and 4.3). Second, the 5 σ tension between Planck CMB and Lyα forest data might prefer a fraction ≈ (1 -5) % of FDM with particle mass m ∼ 10 -25 eV, which is very close to the axion mass that we probe (Rogers et al. 2023). Third, a fraction f ∼ 10 % of FDMwas suggested to explain the suppressed amplitude of the matter power spectrum at late times (also known as the σ 8 / S8 tension, Allali et al. 2021; Ye et al. 2023).  \nIn this work, our aim is to push both the numerical and (semi-) analytical frontiers of MDM modelling. Alongside implementing a new MDM gravity solver and conducting a suite of simulations within the FDM window for a particle mass of m = 10 -24 . 5 eV = 3 . 16 × 10 -25 eV, we refine the halo model framework A xion HMcode (Vogt et al. 2023), calibrating some of its key parameters to the MDM simulations. We will achieve this mainly by modelling the axion halo mass-cold halo mass relation M a( M c) as a broken power law below the linear Jeans mass and calibrating a generalised transition smoothing parameter α .  \nThis paper is organised as follows. We start with the physics of MDM in Sec. 2.1 and describe the numerical methodology in Sec. 2.5, which we implement in A xi REPO (May & Springel 2021, 2023). Key findings regarding large-scale structure and halo statistics are presented in Sec. 3. We detail our improvements to MDM halo model calibrations in Sec. 4 before we conclude in Sec. 5.", '2.1 Mixed Dark Matter Physics': 'We define the FDM fraction (of the total matter density) as  \nf ≡ Ω FDM Ω FDM + Ω CDM + Ω b , (1)  \nwhere we adopt Planck Collaboration et al. (2016) cosmological parameters. In this paper, the simulation-based analyses (up to Sec. 4) rely on DM-only modelling, where baryons are e ff ectively absorbed into the CDM component and assumed to behave similarly to CDM, such that Ω CDM ← Ω CDM + Ω b and Ω b ← 0. In other words, for the purposes of our simulations, we set Ω b = 0 and increase Ω CDM correspondingly.  \nIn DM-only models (as in our simulations) this definition of f is standard (e. g. Laguë et al. 2023), and f max = 1 corresponds to pure FDM. In models that distinguish between the CDM component and baryons, our definition (1) di ff ers from the more commonly used ˜ f = Ω FDM / ( Ω FDM +Ω CDM), which is also employed in the Vogt et al. (2023) MDM halo model central to this paper. In such cases, the maximum FDM fraction is f max = 0 . 843. Converting between the more common definition and Eq. (1) simply involves downscaling ˜ f by the factor ( Ω FDM + Ω CDM) / ( Ω FDM + Ω CDM + Ω b).  \nWe consider axions and axion-like particles, characterised by a periodicity defined by the energy scale f a and represented as an angular variable ϕ . The initial, random value of the misalignment angle θ a = ϕ/ f a (or the distribution thereof in our Hubble patch) determines the FDM abundance Ω FDM. In case of the QCD axion, the associated global U(1) symmetry is the Peccei-Quinn (PQ) symmetry whose angular pseudoscalar field ϕ couples to the strong sector and via non-perturbative QCD e ff ects (instantons) dynamically relaxes the θ term in the QCD Lagrangian to zero, thus \'solving\' the strong-CP problem (Peccei & Quinn 1977). Similar nonperturbative e ff ects also exist in string-theoretical realisations of axions (brane instantons, Svrcek & Witten 2006; Svrcek 2006; Blumenhagen et al. 2009) and in both cases lead to the spontaneous breaking of the exact shift symmetry ϕ → ϕ + const. at scale Λ ≪ f a, rendering axions and axion-like particles massive (hence the term pseudo-Nambu-Goldstone boson). The potential produced by the non-perturbative e ff ects respects the periodicity of ϕ and is usually written as  \nV ( ϕ ) = Λ 4 " 1 -cos ϕ f a !# , (2)  \nwhere the axion mass can be expressed as m = Λ 2 / f a in case of string theory models and using temperature-dependent terms in case of QCD axions (Niemeyer 2020).  \nIn the limit of small field displacements away from the potential minimum, ϕ ≪ f a, and ignoring the Hubble drag in a late-time Friedmann-Lemaître-Robertson-Walker (FLRW) background cosmology, the axion satisfies the Klein-Gordon equation (in natural units),  \n( -∂µ∂ µ + m 2 ) ϕ = 0 . (3)  \nHere we are interested in the case where ϕ has both spatial and temporal fluctuations. On scales much larger than the Compton wavelength (in natural units)  \nλ c ≡ 2 π m , (4)  \nbut much smaller than the particle horizon, one can employ a nonrelativistic approximation to the dispersion relation and a Newtonian approximation to the gravitational interaction embedded in the covariant derivatives (since ∂µ∂ µ = g µν ∂µ∂ν ) of the field equation.  \nIt is then convenient to define a complex scalar field  \nψ = √ ρ FDM e i γ , (5)  \nthe \'classical\' wavefunction, constructed from the amplitude and phase of the field ϕ ,  \nϕ = √ ρ FDM cos( mt -γ ) . (6)  \nThe wavefunction then obeys  \ni ∂ t + 3 2 ˙ a a ! ψ = -1 2 m ∇ 2 + m Ψ ! ψ, (7a)  \n∇ 2 Ψ = 4 π G ( ρ CDM + ρ FDM -¯ ρ tot) , (7b)  \nwhere Ψ is the Newtonian gravitational potential, the FDM density can be expressed as ρ FDM = | ψ | 2 and ¯ ρ tot is the mean of the total DM density. This is simply the non-linear Schrödinger-Poisson (SP) system of equations for a self-gravitating many-body field in a potential well, embedded in an expanding Universe. Note that the right hand side of the first equation of system (7) vanishes for the unperturbed background and the energy density in the axion field ρ FDM = | ψ | 2 ∝ a -3 redshifts like matter. In terms of comoving coordinates x ≡ r / a , the SP equations become  \ni ∂ψ c ∂ t = -a -2 1 2 m ∇ 2 c ψ c + a -1 m Ψ c ψ c , (8a)  \n∇ 2 c Ψ c = 4 π G ( ρ CDM,c + ρ FDM,c -¯ ρ m,c) , (8b)  \nwhere we have defined comoving quantities, relating to physical quantities as  \nρ i , c ≡ a 3 ρ i , ψ c ≡ a 3 / 2 ψ, ∇ c ≡ a ∇ , Ψ c ≡ a Ψ , (9)  \nfor i ∈ { CDM , FDM , m } , where \'m\' denotes total matter.', '2.2 Characteristic Scales in MDM': "Afundamental length scale in the context of MDM is the linear Jeans scale (Hu et al. 2000; Marsh & Silk 2014),  \nk J = GLYPH<16> 16 π Ga 4 ¯ ρ m GLYPH<17> 1 / 4 GLYPH<18> m ℏ GLYPH<19> 1 / 2 ≃ 66 . 5 (1 + z ) -1 / 4 Ω m h 2 0 . 12 ! 1 / 4 GLYPH<18> m 10 -22 eV GLYPH<19> 1 / 2 cMpc -1 . (10)  \nThis scale emerges from the non-vanishing e ff ective sound speed of scalar fields, which introduces a Jeans-like dispersion relation for perturbation modes. Notably, the linear Jeans scale exhibits minimal dependence on redshift and is una ff ected by the axion fraction f . The corresponding linear Jeans mass can be written as  \nM J = 4 3 π π k J ! 3 ¯ ρ m ≈ (1 -2) × 10 11 M ⊙ / h for z = 1 -4 , (11)  \nwhere we have assumed spherical basis functions, and as usual we have taken m = 3 . 16 × 10 -25 eV. In pure FDM cosmologies ( f = f max), no halos can form with mass M < M J due to the pressure  \nTable 1. Scaling of several characteristic MDM scales with axion mass m , redshift z , and axion fraction f . The best-fit results for the scaling of the cuto ff scale k cut with axion mass, redshift, and axion fraction were obtained by evaluating Eq. (15) and comparing it to the virial radius in the redshift range z = 1 -4, axion mass range m = 10 -26 -10 -23 eV, and axion fraction range f = 0 . 01 -0 . 3, respectively. These length scales show no dependence on the axion fraction f and are useful scales to consider for both FDM and MDM cosmologies.  \nsupport from the scalar field. Some residual structures (not halos) may still form via fragmentation. However, this estimate is based on a purely linear e ff ect and indeed the growth rate is enhanced at second order (Li et al. 2019), hence heuristic reasoning based on experience from the baryonic Jeans length can only be applied with caution. Schive et al. (2014b) derived a slightly more accurate minimum halo mass from the properties of FDM solitons that can form, but their estimate is in close agreement with the linear Jeans mass M J.  \nThe Jeans mass can be generalised to non-linear structures, by adding the dependence on the halo profile ρ NFW of the CDM halo. Hu et al. (2000) and Marsh & Silk (2014) showed that in a pure axion cosmology, no virialised halo is expected to form if this socalled halo Jeans scale r hJ is larger than the virial radius R vir. The corresponding mass scale is called the cut-o ff mass M cut. In comoving units, the halo Jeans wavenumber is given by (cf. Eq. 10)  \nk hJ = 66 . 5 (1 + z ) -1 / 4 Ω m h 2 0 . 12 ! 1 / 4 GLYPH<18> m 10 -22 eV GLYPH<19> 1 / 2 ρ NFW( r hJ) ¯ ρ m ! 1 / 4 cMpc -1 . (12)  \nHere, k hJ is calculated by converting r hJ = π/ k hJ as in Marsh & Silk (2014) instead of r hJ = 2 π/ k hJ as in Hu et al. (2000). Since we assume that r hJ ≤ R vir, the Navarro-Frenk-White (NFW) density at r hJ can be written as (Navarro et al. 1997)  \nρ NFW( r hJ) = ¯ ρ m ∆ vir c 2 3 f ( c ) R vir r hJ(1 + r hJ c / R vir) 2 , (13)  \nwith f ( x ) = -x / (1 + x ) + ln(1 + x ) and c = R vir / r -2 being the halo concentration. The virial mass definition we adopt is  \nM vir = M tot = 4 π 3 ¯ ρ R 3 vir ∆ vir( z ) , (14)  \nwhere ∆ vir is the virial overdensity (Bryan & Norman 1998). Substituting Eq. (13) back into Eq. (12), we obtain  \nr 3 hJ = GLYPH<18> π 66 . 5 GLYPH<19> 4 (1 + z ) GLYPH<18> m 10 -22 eV GLYPH<19> -2 Ω m h 2 0 . 12 ! -1 3 f ( c ) ∆ vir c 2 (1 + r hJ c / R vir) 2 R vir , (15)  \nwhere r hJ and R vir are measured in ckpc and we have assumed spherical basis functions as we did for Eq. (11). Eq. (15) can then be solved for r hJ. The dependence on axion fraction f and redshift z only enters via the concentration c which we take to be the mean cosmic concentration at mass M vir as per Eqs. (43). In fact, after comparing to the virial radius, we find that the cut-o ff mass is e ff ectively independent of f , and attains values M cut ≈ 1 . 7 × 10 12 -3 . 1 × 10 13 M ⊙ / h for our MDM cosmologies across redshifts z = 1 -4.  \nHow do we interpret the cut-o ff mass M cut? In pure FDM, while the Jeans mass M J provides a fundamental lower bound on the mass  \nof any halo that can form based on the balance between gravity and 'quantum pressure' in the linear regime, the cut-o ff mass is a halospecific lower bound that describes the minimum mass of virialised halos that can form. In MDM cosmologies, virialised halos can form below M cut, but they will be primarily CDM-dominated.  \nAnother important scale is the characteristic mass at which the non-linear matter power spectra and HMFs begin to deviate from the predictions of the CDM model. This characteristic mass, denoted as M 0, is given by (Schive et al. 2016)  \nM 0 = 1 . 6 × 10 10 ( m / 10 -22 eV) -4 / 3 M ⊙ (16)  \nThe scaling of M 0 is expected to be almost independent of redshift, as it is primarily set during the radiation-dominated epoch (Hu et al. 2000).  \nFinally, we highlight a fourth mass scale, M conc ≈ 20 -50 M 0, at which the mean halo concentration-mass relation c ( M ) is significantly a ff ected by axion physics. This was demonstrated by Bose et al. (2016); Ludlow et al. (2016); Dentler et al. (2022) and can be attributed to the fact that halo concentration is determined when only a small fraction ( ≈ 0 . 01) of the mass has accumulated (Navarro et al. 1997; Bullock et al. 2001). These four mass scales are illustrated in Fig. 1, using a fiducial MDM cosmology with axion mass m = 3 . 16 × 10 -25 eV and axion fraction f = 0 . 1 at redshift z = 1.  \nFor reference, we show the scalings of these four characteristic length scales with axion mass m , redshift z and axion fraction f in Table 1. We translate between mass and inverse length scales via M = 4 π ¯ ρ R 3 / 3. Among the four scales, the cut-o ff scale k cut shows the strongest dependence on axion mass, scaling as k cut ∝ m 0 . 66 . Even though these scales were first introduced for pure FDM cosmologies and thus have no dependence on the axion fraction (while M cut is e ff ectively independent of f ), they remain valuable to consider in the context of MDM cosmologies as well. Specifically, we will show in Sec. 3.1 that despite its weak dependence on the axion fraction, M 0 serves as an e ff ective scale for characterising the suppression of the halo mass function (HMF) at the low-mass end. In Sec. 3.3, we will show that the weak dependence of M conc on the axion fraction is corroborated by inferred concentration-mass relations. In addition, we will establish in Sec. 4.3 that M cut provides a useful approximate scale below which halos are predominantly CDM-dominated.", '2.3 Initial Conditions': "The goal of an initial conditions (ICs) generator for a cosmological simulation is to faithfully reproduce the statistical properties of the density field in the early universe with a finite number of point particles (for the CDM component) or grid cells (for the FDM component). Given a pre-IC particle distribution representing a completely homogeneous universe (we use grid pre-ICs, see Dome et al. 2022), the next task is to perturb the particles to produce a density field that reproduces the cosmologically relevant expected statistical properties. Let us consider a (Gaussian) over-density field δ ( r ) that is completely described by its power spectrum P ( k ) = ⟨ δ k δ ∗ k ⟩ . It is customary to express the amplitude of density fluctuations in terms of the transfer function T ( k , a ),  \nP ( k , a ) ≈ Bk ns T ( k , a ) 2 , (17)  \nwhere n s is the spectral index, and B can be expressed in terms of the normalisation A s and the pivot scale k piv (Planck Collaboration et al. 2016).  \nSetting up ICs for cosmological simulations at a certain scale factor a thus involves generating a white noise sample of random val-  \nFigure 1. Characteristic scales in the MDM framework. We illustrate the linear Jeans mass M J (see Eq. 10), the cut-o ff mass M cut (see Eq. 15 and associated text), and the characteristic mass M 0, where deviations from CDM become apparent in the matter transfer function (TF) and HMF (see Eq. 16). Additionally, we identify M conc, the mass scale at which the concentration-mass relation c ( M ) shows notable changes. Results are presented for a fiducial MDM cosmology with axion mass m = 3 . 16 × 10 -25 eV and axion fraction f = 0 . 1 at redshift z = 1.  \n<!-- image -->  \n/circledot  \nues µ ( r ) (typically sampled from a Gaussian distribution with zero mean and unit variance) and requiring that their amplitudes follow a specific power spectrum P ( k , a ). This is achieved by multiplying the Fourier transformed white noise field µ k with the square root of the power spectrum, i. e. for all k representable on a grid of given resolution set  \nδ k = p P ( k , a ) µ k = √ Bk ns / 2 T ( k , a ) µ k . (18)  \nThe real-space over-density field δ ( r ) is then obtained by inverse Fourier transformation, and this procedure is typically called ' k -space sampling'.  \nNote that a product in Fourier space simply corresponds to a convolution in real space, i. e. Eq. (18) is equivalent to  \nδ ( r ) = T ( r , a ) ⋆ µ ( r ) , (19)  \nwhere T ( r , a ) is the real-space counterpart of T ( k , a ) = √ Bk ns / 2 T ( k , a ), and ' ⋆ ' denotes a convolution. It is hence mathematically equivalent whether Eq. (18) is evaluated in Fourier space, followed by an inverse transform ( k -space sampling), or whether Eq. (19) is evaluated using an inverse transform of T ( k , a ) followed by the convolution (real-space sampling). Most cosmological IC codes follow the first approach (see e. g. Bertschinger 2001), while e. g. Pen (1997) and Sirko (2005) use the second or variations thereof.  \nHowever, the discrete realisations of the density fields derived with the two approaches will have significant di ff erences. This has been demonstrated conclusively by Sirko (2005), who showed amongst others that employing Eq. (18) imposes periodicity of the real-space transfer function on box scales and leads to an underestimation of the two-point correlation function on large (sub-box) scales. This discrepancy arises because, with a finite number of particles in a finite box, achieving high accuracy for estimates of both the correlation function and the power spectrum is challenging due to finite box e ff ects and the limitations in resolving all relevant scales. Even with an infinite number of particles, the finite volume of the box introduces periodic boundary conditions and aliasing e ff ects that can distort the two-point correlation function, particularly at large scales. Importantly, in this work we use the public code M usic to generate ICs (Hahn & Abel 2011), which is likewise based on a convolution of Gaussian white noise with a real-space transfer function  \nkernel. The linear MDM power spectrum which serves as an input to M usic is calculated using A xion C amb (Hlozek et al. 2015).  \nHaving generated the Eulerian-grid seed density field δ ( r ) using real-space sampling while ensuring its consistency with both the input correlation function and the input power spectrum, we need to calculate the initial positions and velocities of the IC macroparticles. This is typically achieved using first (1LPT) or second order Lagrangian perturbation theory (2LPT), where δ ( r ) is used as the source field for Lagrangian perturbation theory. The displacement field in 1LPT only contains contributions from the gravitational potential, which is often called the Zel'dovich approximation (Zel'dovich 1970). For all the simulations analysed in this work, we employ a more accurate representation by incorporating secondorder e ff ects using M usic .  \nHow do we take account of the multi-fluid nature of MDM beyond calculating separate density transfer functions for each component and coupling them gravitationally in the 2LPT? The growth of density perturbations in a two-component fluid can only be correctly reproduced if besides the di ff erent initial amplitudes of density perturbations also the di ff erence in initial velocities between the two components are respected (Yoshida et al. 2003). ICs for the twocomponent fluid thus ought to reflect these important di ff erences between the two components. 1 The correct growth of fluctuations in both components consistent with the predictions from linear perturbation theories can only be achieved with a proper modelling of (relative) velocities (Somogyi & Smith 2010; Laguë et al. 2021).  \nTo obtain a velocity transfer function, recall that the curl of a peculiar velocity field, ∇× v ∝ a -1 , drops o ff with the expansion of the Universe and can be neglected at late times since there is no source for vorticity (although see Mocz et al. 2017; Hui et al. 2021). Hence we can write v as the gradient of a velocity potential, v = ∇V , and so  \nvk = i k V k . (20)  \nLinear Newtonian perturbation theory applied to a pure DM cosmol-  \n1  \nFigure 2. Velocity transfer functions (squared) in an MDM cosmology consisting of a majority CDM component with an f = 0 . 1 axion admixture of boson mass m = 10 -24 . 5 eV at z = 127. The unapproximated FDM (blue solid) and CDM (green solid) velocity transfer functions in Newtonian gauge (NG) are calculated using a modified version of A xion C amb . The vertical dashed line marks the FDM linear Jeans scale k J (Eq. 10). The approximated Newtonian gauge FDM velocity, derived from MDM matter transfer functions (see Eq. 24), is represented by the orange solid line. This L ( k ) model provides a poor approximation to the blue solid line on small scales k > 5 cMpc -1 . For completeness, we also show the (unapproximated) FDM velocity transfer function in synchronous gauge (SG), calculated similarly using a modified version of A xion C amb . Note that CDM SG is zero by definition, and so does not appear.  \n<!-- image -->  \nogy (either CDM or FDM) thus implies that  \nvk = ia k k 2 d δ k d t . (21)  \nFor the growing mode of δ , we obtain  \nvk = ia k k 2 Ha δ k F ( Ω m) , (22)  \nwhere  \nF ( Ω m) ≡ d ln( D + ( k , a )) d ln( a ) . (23)  \nNote that Eq. (22) remains valid even in MDM cosmologies, where we have one such equation for each DM component. In contrast to pure CDM, the FDM growth rate depends on scale k , which arises due to the non-vanishing sound speed of the axion fluid (Hwang & Noh 2009; Marsh 2016).  \nIt is nonetheless possible to approximately separate out the dependence on scale k following Laguë et al. (2021),  \nD FDM + ( k , a ) ≈ L ( k ) D CDM + ( a ) . (24)  \nNote that we can express L ( k ) as  \nL ( k ) ≈ D FDM + ( k , a ) D CDM + ( a ) = s ∆ FDM ( k , a ) ∆ CDM ( k , a ) , (25)  \nwhere ∆ ( k , a ) ≡ T ( k , a ) 2 is the squared transfer function of the respective DM component. In CDM, recall that T ( k , a ) = T ( k ) D + ( a ).  \nWe can now approximate the velocity field of the axion component in MDM cosmologies using the CDM velocity field as  \nv FDM k ≈ L ( k ) v CDM k . (26)  \nTo approximate the squared velocity transfer function ∆ FDM v ( k , a )  \nfor the axion component in an MDM cosmology, it should thus be su ffi cient to calculate ∆ CDM v ( k , a ) in the corresponding pure CDM cosmology and construct the squared ratio of FDM-to-CDM matter transfer functions ∆ FDM ( k , a ) / ∆ CDM ( k , a ). Assuming separability of D FDM + ( k , a ) based on Eq. (25), we have  \n∆ FDM v ( k , a ) ≈ ∆ FDM ( k , a ) ∆ CDM ( k , a ) ∆ CDM v ( k , a ) . (27)  \nNow, both ∆ CDM v ( k , a ) and the squared transfer function ratio can be generated easily using the public version of A xion C amb (Hlozek et al. 2015). We call Eq. (27) the L ( k ) model since it relies on the assumption of separability. To assess the fidelity of the L ( k ) model, we have modified A xion C amb such that we can retrieve the unapproximated FDM velocity transfer function ∆ FDM v ( k , a ) in Newtonian gauge directly from the Boltzmann code. In Fig. 2, we show velocity transfer functions in an MDM cosmology with a 10 % admixture of a m = 10 -24 . 5 eV axion field. We see that the L ( k ) model provides a poor approximation for ∆ FDM v ( k , a ) on small scales k > 5 cMpc -1 in Newtonian gauge, 2 justifying our more reliable, unapproximated approach to velocity transfer functions.  \nEquipped with the velocity and matter transfer functions of each component, we use the public code M usic to generate ICs (Hahn & Abel 2011). M usic takes the transfer functions at a particular redshift (in our case, z = 127) as input, and generates initial positions and velocities for macroparticles of each component.", '2.4 Madelung Formulation': "How do we initialise the FDM wavefunction on a grid from macroparticle positions and velocities? We use the Madelung change of variables (Madelung 1927), which can aid with physical intuition. We start with the decomposition (5) of the wavefunction into its amplitude and phase, ψ c = √ ρ FDM,c e i γ , and define the Madelung velocity as the gradient of the phase,  \nv M ≡ ∇ γ m . (28)  \nOn linear scales before shell crossing, the fluid velocity is a gradient flow, and it resembles that of a superfluid. The Schrödinger equation can then be written as  \n∂ρ FDM,c + ∇ c · ( ρ FDM,c v M) = 0 ,  \n∂ t (29a) ∂ v M ∂ t + a -2 v M · ∇ c v M = -a -1 ∇ c Ψ c + a -2 1 2 m 2 ∇ c ∇ 2 √ ρ FDM,c √ ρ FDM,c ! . (29b)  \nThe Schrödinger equation possesses a U(1) symmetry, which amounts to the rotation of ψ by a phase. With the identification of the fluid velocity, what is normally understood as probability conservation (i. e. conservation of the associated Noether current) in quantum mechanics is recast as mass conservation in Eq. (29a). The last term in the Euler equation (29b) is often referred to as the 'quantum pressure' term. It is a misnomer since we have a classical system. In addition, the term arises from a stress tensor rather than mere pressure:  \nΣ i j = 1 4 m 2 GLYPH<16> ρ -1 ∂ i ρ∂ j ρ -∂ i ∂ j ρ GLYPH<17> = -ρ 4 m 2 ∂ i ∂ j ln( ρ ) , (30)  \ni. e. ∂ i ( ∇ 2 √ ρ / √ ρ ) / (2 m 2 ) = -ρ -1 ∂ j Σ i j . We have dropped the subscripts on the FDM density field in Eq. (30), ρ = ρ FDM,c. The stress tensor Σ i j represents how the fluid description accounts for the underlying wave dynamics. It shows how the particle limit is obtained: for large m , the Euler equation reduces to that for a pressureless fluid, as is appropriate for particle DM. The insight that the wave formulation in the large m limit can be used to model particle CDM was exploited by Widrow & Kaiser (1993). The wave description e ff ectively reshu ffl es information in a phase-space Boltzmann distribution into a position-space wavefunction and o ff ers a number of insights that might otherwise be obscured (Uhlemann et al. 2019; Garny et al. 2020). This correspondence can be formalised (Mocz et al. 2018). It is worth noting that A xion C amb also employs the Madelung description.  \nTo initialise the FDM wavefunction, we calculate the FDM phase γ by constructing the Madelung velocity field v M and solving Eq. (28) in Fourier space. Equipped with CDM macroparticle positions and velocities as well as the FDM wavefunction, we can now evolve the joint MDM field.", '2.5 Pseudo-Spectral Method': 'We use a spectral method to simulate MDM structure formation implemented in the A xi REPO code (May & Springel 2021, 2023). The system of equations (8) is solved using a second-order symmetrised split-step pseudo-spectral Fourier method, colloquially called a \'kick-drift-kick\' leapfrog-like scheme. For a small time step ∆ t , the time evolution can be simplified using the following approximation (Edwards et al. 2018; May & Springel 2021):  \nψ c( t + ∆ t , x )  \n= T exp " -i Z t +∆ t t -ℏ 2 m 1 a ( t \' ) 2 ∇ 2 c + m ℏ 1 a ( t \' ) Ψ c( t \' , x ) ! d t \' # ψ c( t , x ) ≈ exp " i ∆ t 2 -ℏ m 1 a ( t ) 2 ∇ 2 c -m ℏ 1 a ( t ) Ψ c( t + ∆ t , x ) -m ℏ 1 a ( t ) Ψ c( t , x ) !# × ψ c( t , x ) ≈ exp " -i m ℏ 1 a ( t ) ∆ t 2 Ψ c( t + ∆ t , x ) # exp " i ℏ m 1 a ( t ) 2 ∆ t 2 ∇ 2 c # × exp " -i m ℏ 1 a ( t ) ∆ t 2 Ψ c( t , x ) # ψ c( t , x ) , (31)  \nwhere T is the time ordering operator and, using the BakerCampbell-Hausdor ff formula, the time evolution operator has been split into three unitary parts which do not mix functions of the position and derivative operators. This makes it natural to automatically couple the method to particle-based N -body techniques that evolve collisionless components such as CDM and stellar particles on the same sub-time step spacing. Coupling to gas cells is also straightforward and is achieved via the full gravitational potential Ψ c (including the baryonic contribution) in both the SP equations and the forces evolving the gas cells. Simulations involving mixed ultralight and baryonic physics are left for future work.  \nThe fields ψ c and Ψ c are discretised on a uniform Cartesian grid with N 3 mesh points in a periodic box of length L to allow for e ffi cient numerical computations using the Fast Fourier Transform (FFT). The pseudo-spectral method is summarised in Algorithm 1. The choice of the time step ∆ t is determined by the requirement that the phase di ff erence in the exponentials must not exceed 2 π , at which point the time step would be incorrectly \'aliased\' to a smaller time step corresponding to the phase di ff erence subtracted by a multiple of 2 π due to the periodicity of the exponential function. The kicks  \nAlgorithm 1 Pseudo-Spectral Method for MDM  \nand the drift yield separate constraints for ∆ t , both of which must be simultaneously fulfilled. The resulting time step criterion is  \n∆ t < min 4 3 π m ℏ a 2 ∆ x 2 , 2 π ℏ m a 1 | Ψ c , max | ! , (32)  \nwhere ∆ x = L / N is the spatial resolution and Ψ c , max is the maximum value of the potential. Note that Eq. (32) is essentially a CourantFriedrichs-Lewy (CFL) condition. The dependence ∆ t ∝ ∆ x 2 can be viewed as a reflection of the relation of the Schrödinger equation to di ff usion problems. Since N -body codes for gravity and Eulerian fluid solvers scale as ∝ ∆ x , this adds computational cost to the simulations. However, Algorithm (1) allows for machine precision control of the total kinetic energy and achieves spectral (exponential) convergence in space.  \nWhen considering the velocity field v c = v M = ℏ ∇ c γ/ m , another constraint on the validity of the discretisation becomes apparent. Since the di ff erence in the gradient of the phase between two points can be at most 2 π , it follows that the discretised velocity field cannot exceed a maximum value (depending on the concrete form of the discretised gradient operator) of about  \nv max = ℏ m 2 π ∆ x . (33)  \nVelocities v ≥ v max cannot be represented in a simulation with resolution ∆ x , which translates into a constraint on resolution, which should be good enough to resolve the de Broglie wavelength λ dB of the largest velocities:  \n∆ x < π ℏ mv max ≡ 1 2 λ dB( v max) . (34)  \nRequirements (32) and (34) exemplify why FDM simulations are computationally much more costly than traditional particle-based CDM simulations, where resolution can be set independent of velocities and time step constraints are less restrictive (e. g. Springel 2005).', '2.6 Simulation Setup': 'As mentioned, in the MDM model only one of the DM particle species - occupying a fraction f of the total matter content - is ultralight while all others are assumed to have negligible de Broglie wavelengths and are thus modelled as CDM. Schwabe et al. (2020) pointed out that a soliton may not form in cosmological simulations with f < 0 . 1. This result was confirmed by Laguë et al. (2023) who investigated the impact of a mixture of CDM and FDM in various proportions f = [0 , 1 , 10 , 50 , 100] % and for ultralight particle masses ranging over five orders of magnitude (2 . 5 × 10 -25 eV - 2 . 5 × 10 -21 eV) using A xio N yx , albeit mostly for relatively small box sizes, L box = 1 cMpc / h . The authors also implemented a modified friends-of-friends (FOF) halo finder and found good agreement between the inferred halo abundance and the predictions from the adapted halo model A xion HM code in a narrow mass range. Expanding upon their findings, we aim to identify distinctive characteristics  \nTable 2. Overview of the MDM simulation suite: (1) FDM fraction f ; (2) side length of the simulation box L box; (3) number of CDM N -body particles; (4) number of FDM grid cells; (5) mass per CDM N -body particle. The softening scale for CDM is fixed to ϵ = 1 . 78 ckpc / h in comoving units and capped at ϵ = 0 . 89 pkpc / h in physical units following Power et al. (2003). All simulations are DM-only and assume an axion mass of m = 10 -24 . 5 eV = 3 . 16 × 10 -25 eV.  \nthat set MDM apart from single-particle models, thereby enhancing our understanding of its cosmological implications.  \nFor our MDM simulations, we adopt an axion mass of m = 10 -24 . 5 eV = 3 . 16 × 10 -25 eV and vary the FDM fraction in the set f ∈ [0 . 01 , 0 . 1 , 0 . 2 , 0 . 3]. Note that due to the presence of a dominant CDM component ( f < 0 . 5), the potential wells in which the wavefunction evolves become steeper, which increases the FDM velocity dispersion and decreases its de Broglie wavelength λ dB, Eq. (34). MDM simulations thus require a higher resolution than their pure FDM ( f = f max) counterparts. Consequently, conducting large-scale cosmological simulations involving the solution of the full Schrödinger-Poisson system within the MDM framework is challenging. While we ascertain the necessary CDM and FDM resolution requirements according to the criteria outlined in Sec. 2.5 and Dome et al. (in prep), we thus ensure at the same time that our selection of simulation parameters is conservative. The final specifications of our MDM simulation suite are given in Table 2.  \nTo ensure the fidelity of the high redshift evolution of the simulations, we performed several tests: we verified that across several orders of magnitude the partial (FDM and CDM components) and total power spectra of the first snapshot ( z = 127) replicate the target power spectra obtained using the modified version of A xion -C amb (see Sec. 2.3); we also tested the linear growth prediction D + ( k , a ) ∝ a (from linear theory) which holds on large scales for both partial and total power spectra. Finally, we exclude snapshots with redshifts below z ≲ 1, as even with our conservative selection of simulation parameters, there is a risk that their highest velocity dispersions may remain unresolved by the FDM solver (see Eq. 34).  \nWe visualise the MDM density distribution for an axion fraction f = 0 . 1 in Fig. 3, juxtaposing the FDM and CDM fields. The clustering of the former is visibly suppressed below the axion Jeans scale k J, Eq. (10). The smallest-mass halos can remain completely CDMdominated. As known from the pure FDM case ( f = f max), the wave dynamics are reflected in interference fringes along filaments as well as granule structures in the DM halos.', '3 HALO MASS DISTRIBUTION AND DENSITY PROFILES': 'In the following, we present a study of HMFs and halo density profiles in MDM cosmologies. To identify halos, we use the R ockstar halo finder (Behroozi et al. 2013), which is based on adaptive hierarchical refinement of FOF groups in six dimensions (position and momentum space). This method provides robust tracking of substructure, being grid-independent, orientation-independent, and resilient to noise. However, R ockstar is a particle-based halo finder and thus  \ncan only be applied to particle distributions. Although it is theoretically possible to convert FDM grid cells into particles by concentrating the mass into points at the centre of each cell, we found 3 that this approach leads to unreliable FOF groupings due to artifacts introduced by the grid features in the particle distribution.  \nGiven that the dominant component in our MDM simulations is CDM ( f < 0 . 5), we base our halo identification exclusively on the CDM component. This is achieved by running R ockstar on the (equal-mass) CDM particle distribution after correcting the CDM particle mass m CDM such that the total mass of CDM particles present within the simulation box aligns with the anticipated total DM mass. 4 We will justify the usage of R ockstar on the CDM component a posteriori.', '3.1 Halo Mass Functions': "To measure the halo abundance in the simulations, we choose logarithmic mass bins of width ∆ log( M tot) = 0 . 03 and count the number of R ockstar halos that fall into each bin. We show the resulting HMFs in Fig. 4 at redshifts z = 1 -4. While CDM and the f = 0 . 01 cosmology follow the bottom-up structure formation paradigm in which small-mass halos form first and thus the small-mass end of their HMF barely changes from z = 4 to z = 1, MDM cosmologies with f ≳ 0 . 1 violate this picture and many small-mass halos are assembled at low redshift, with corresponding changes in the smallmass amplitude of their HMF.  \nWe now compare against theoretical predictions for the HMF by employing the Sheth-Tormen approach for the halo mass function (Press & Schechter 1974; Sheth & Tormen 2002),  \n1 M d n d ln( M ) = 1 2 ¯ ρ ( z ) M 2 f ( ν ) GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> d ln( σ 2 ) d ln( M ) GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> , (35)  \nwhere n is the halo number density, ν = δ crit /σ ( M , z ) is the peak height with δ crit = 1 . 686 the critical linear density threshold for halo collapse, and the multiplicity function (for ellipsoidal collapse) is given by  \nf ST( ν ) = A r 2 π √ q ν (1 + ( √ q ν ) -2 p ) e -q ν 2 2 , (36)  \nwith A = 0 . 3222, p = 0 . 3 and q = 0 . 707. We adopt a spherical top hat window function W in real space when calculating the variance of the linear MDM power spectrum, 5  \nσ ( R , z ) 2 = 1 2 Z ∞ P L ( k , z ) ˜ W ( Rk ) 2 k 2 d k ,  \n˜ W ( x ) = 3 3 (sin x -x cos x ) .  \n2 π 0 (37) x (38)  \nThe variance can be transformed into a function of the halo mass via M = 4 π ¯ ρ R 3 / 3.  \nWe show these model predictions as dashed curves in Fig. 4 and find a reasonable agreement between the model and the simulations,  \nFigure 3. MDMstructure formation for m = 3 . 16 × 10 -25 eV and f = 0 . 1. We plot projected (comoving) densities of the FDM (left column) and CDM (right column) components along the line of sight on a logarithmic scale. The redshift spacings (di ff erent rows) correspond to equal logarithmic spacings in the scale factor. The axion Jeans scale prevents FDM from clustering around small-mass halos, which thus remain CDM-dominated. The wave dynamics of FDM is reflected in interference fringes along filaments as well as the presence of granule structures in the DM halos, most visible towards lower redshift.  \n<!-- image -->  \nFigure 4. MDM halo mass function (HMF) across redshifts z = 1 -4. We show results for a m = 3 . 16 × 10 -25 eV MDM cosmology at various FDM fractions f (see legend) and the reference CDM cosmology. Solid lines with markers show the HMF inferred from the simulation. The dashed curves trace the Sheth-Tormen HMF (based on linear MDM power spectra), which show good agreement with the inferred HMF. Dotted curves are best-fit results against Eq. (39), where we keep α 1 = 1 . 1 and α 2 = 2 . 2 fixed during the fit. The vertical dashed line corresponds to the mass resolution 50 × m ' CDM , where m ' CDM is the rescaled CDM particle mass (see text).  \n<!-- image -->  \n/circledot  \n/circledot  \nTable 3. Best-fit values of β in Eq. (39) for our MDM simulation suite across a range of redshifts z = 1 -8. Note that the dependence on f is substantial at higher redshift z ≳ 3. For f = 0 . 01 the exact value of β is less relevant since only the combination β f enters the parametrisation (39).  \nat all redshifts shown z = 1 -4 and in fact out to z = 10 (not shown). The model curves for f ≳ 0 . 1 branch o ff from the CDM curves at the characteristic mass M 0 (see Eq. 16) below which the FDM HMF  \nis suppressed with respect to the CDM case. While M 0 has been introduced for pure FDM ( f = f max), our MDM HMFs also branch o ff at approximately M 0 = 2 . 34 × 10 13 M ⊙ / h , suggesting that M 0 remains a useful and accurate mass scale even when f < f max.  \nWe also aim to provide an analytic form for the HMFs as a function of the FDM fraction f . To that end, we generalise the twoparameter model of Schive et al. (2016) and write  \nn ( M tot) = n ST ( P c L )       1 -β f + β f 1 + M tot M 0 ! -α 1 ! -α 2       , (39)  \nwhere n ST ( P c L ) is the Sheth-Tormen HMF in CDM. The steepness of the suppression is controlled by the parameter α 2 (but also β for M significantly smaller than M 0) while the sharpness of the transition at M ≈ M 0 is controlled by α 1. Fitting this model to N -body simulations of pure FDM ( f = f max) results in α 1 = 1 . 1, α 2 = 2 . 2 (Schive et al. 2016). Schive et al. (2016) have taken special care when removing spurious halos and have found that the redshift-independent suppression (1 + ( M h / M 0) -1 . 1 ) -2 . 2 provides a good fit across a range of redshifts z > 4. Subsequent works simulating bona fide FDM including the axion wave dynamics (May & Springel 2023) have found that this parametrisation remains robust down to at least z = 3. Our  \nFigure 5. Spherically averaged density profiles in MDM at redshift z = 1. The mass bins range from 10 9 -10 10 M ⊙ / h in the top left to 10 12 -10 13 M ⊙ / h in the bottom right. The shaded areas delineate the standard error of the median while the solid curves trace the median profile in each mass bin. Arrows at the bottom of the panels indicate inner convergence radii r conv = 5 ϵ . The r -1 and r -2 power laws delineate two characteristic regimes of an NFW-like profile. Einasto best-fits are shown by circular markers, and we find that three-parameter Einasto model provides a good parametrisation of MDM density profiles across the entire range of FDM fractions f = 0 . 0 -0 . 3.  \n<!-- image -->  \nMDMHMFsarealso well fit using these reference values for α 1 and α 2, hence we adopt them throughout this work rather than refitting them.  \nThe only parameter that requires tuning is thus β . This parameter is crucial for achieving higher levels of suppression, 1 -β f , to the left of the characteristic mass M 0 than would be allowed by the value of f alone. Best-fit values of β are presented in Table 3. Our results indicate that β typically decreases as f increases and as z decreases. The dependence on f is notably weaker at low redshift ( z ≲ 3). Specifically, the values β = 3 . 1, β = 2 . 9 and β = 2 . 4 provide a good fit at redshifts z = 3, z = 2 and z = 1, respectively. The best-fit results at z = 1 -4 are depicted as dotted curves in Fig. 4. The overall agreement with inferred HMFs is good across axion fractions f = 0 . 01 -0 . 3, indicating that the mass scale M 0 accurately characterises the suppression of the HMF despite being independent of f . The deviation from the inferred HMF at the high-mass end in case of f = 0 . 3 reveals limitations of parametrisation (39) as well as a potential dependence of M 0 on f .", '3.2 Density Profiles': 'A well-known prediction of pure FDM models ( f = f max) is the formation of solitonic cores at the centre of halos (Schive et al. 2014a). The situation is more complex in MDM, where in case of CDMdominated models ( f < 0 . 5) we expect the soliton to be buried beneath the CDM central density cusp. Before we verify this intuitive picture, we first quantify the total DM density profile averages in various mass bins. To better sample the FDM field in the centre of halos, we interpolate FDM density values (which are defined on a Eulerian grid in the simulation code) trilinearly to inter-grid points, while CDMdensity profiles are obtained using simple C osmic P rofiles routines (Dome 2023). The density profiles are calculated around halo centres as identified by R ockstar , which neglects the FDM component. This approach is justified in our CDM-dominated models ( f < 0 . 5), where the CDM cusp dominates the gravitational potential in the centre and shapes its minimum. Note that the FDM core exhibits time-dependent behaviour with no stable centre, undergoing random walks by an amount of the order of the soliton radius on  \ntimescales τ ∝ ( m σ 2 FDM ) -1 , driven by perturbations away from the perfect stationary state (Chiang et al. 2021; Li et al. 2021).  \nThe resulting density profiles in four mass bins from M tot = 10 9 -10 13 M ⊙ / h are shown in Fig. 5, tracing the two characteristic powerlaw regimes of an NFW-like profile (Navarro et al. 1997). The solid curves show the median density profile in logarithmic radial bins of width ∆ log( r ) = 0 . 05 while the shaded area delineates the standard error of the median. Note that in CDM the scale radius rs ≡ r -2 at which the logarithmic slope has the isothermal value of -2, i. e. d ln( ρ ) / d ln( r ) | r -2 = -2, migrates towards larger normalised radii as the halo mass grows. In other words, the concentration  \nc ≡ R vir r -2 (40)  \nof CDM halos decreases as the halo mass increases. This wellknown result reflects the higher background density at earlier epochs when smaller mass halos form (Navarro et al. 1997; Bullock et al. 2001; Ludlow et al. 2014).  \nTo approximate the halo density profiles using an analytic form and to determine their concentration, we invoke the Einasto profile (Einasto 1965)  \nln ρ E( r ) ρ -2 ! = -2 ζ       r r -2 ! ζ -1       (41)  \nand assess the fidelity thereof a posteriori. Best-fit Einasto profiles are determined by adjusting the three parameters ( ζ , r -2 and ρ -2) of Eq. (41) in order to minimise a figure-of-merit defined as  \nψ 2 = 1 N bin N bin X i = 1 GLYPH<2> ln( ρ i ) -ln( ρ E( ri ; ρ -2; r -2; ζ )) GLYPH<3> 2 . (42)  \nWe choose an inner convergence radius of r conv = 5 ϵ , where ϵ is the gravitational softening length. This choice is based on our finding that the 3 ϵ recommendation by Power et al. (2003) is not su ffi ciently conservative, particularly for small-mass bins, potentially leading to an underestimation of halo concentration. In addition, radial bins that exceed the outer limit of 0 . 8 R vir are also discarded in the fit since they might correspond to radii where halos are not fully relaxed (Ludlow et al. 2016). Best-fit results are shown in Fig. 5 and reproduce MDM density profiles very well across the entire range of FDM fractions f = 0 . 0 -0 . 3 at the redshift shown ( z = 1) but in fact out to z = 10 (not shown), albeit at much lower statistical significance.', '3.3 Concentration-Mass Relation': 'We now turn to the concentration-mass relation c ( M , z ) which follows trivially from the Einasto fits. The MDM concentration-mass relation inferred from the simulations at z = 1 is shown in Fig. 6 for logarithmic mass bins of width ∆ log( M tot) = 0 . 2. While the CDM concentration follows the characteristic decrease toward higher mass as mentioned before, this decrease is weaker at higher FDM fractions f and even seems to reverse for f = 0 . 3. To model the concentration mass relation, we invoke the Ludlow et al. (2016) approach based on extended Press-Schechter (EPS) theory. It stipulates that in order to estimate the mean concentration of halos at a given redshift z 0, one needs to solve the following system of coupled non-linear equations,  \n⟨ ρ -2 ⟩ ρ 0 = C H ( z -2) H ( z 0) ! 2 , (43a)  \nM -2 M 0 = erfc        δ sc( z -2) -δ sc( z 0) p 2( σ 2 ( f coll M 0) -σ 2 ( M 0))        . (43b)  \nc  \nFigure 6. Concentration-mass relation c ( M , z ) in MDM cosmologies at redshift z = 1. Square markers indicate median concentrations obtained from Einasto best-fits. Analytic model predictions from Ludlow et al. (2016) for parameters f coll = 0 . 02 and C = 650 (see text) are traced by dashed curves, showing good agreement across f = 0 . 0 -0 . 3.  \n<!-- image -->  \n/circledot  \nThe first equation relates the mean inner density within the scale radius, ⟨ ρ -2 ⟩ , with the critical density of the Universe at the collapse redshift, z -2. The second equation expresses the collapse mass fraction using EPS theory (Lacey & Cole 1993), encapsulating the physical meaning of the collapse redshift as the redshift at which the characteristic mass, M -2, was first contained in progenitors more massive than a fraction f coll of the final halo mass M 0.  \nFollowing the recommendations of Ludlow et al. (2016), we set f coll = 0 . 02, C = 650 and solve Eqs. 43 in the respective MDM cosmologies. Note that the linear variance σ 2 entering the equations depends on the linear MDM matter power spectrum. We show the analytic model predictions in Fig. 6 and find good agreement with the inferred c ( M , z = 1) relation across the entire range of FDM fractions f = 0 . 0 -0 . 3.  \nAccording to the model, the concentration c exhibits a very flat U -shaped profile as a function of mass M tot for cosmologies with high values of f . It indicates that the strict bottom-up structure formation picture in which smaller-mass halos have necessarily higher concentration is invalidated in MDM cosmologies. Note that this result is well-known in pure FDM and warm DM (WDM) cosmologies (Schneider et al. 2012; Dome et al. 2022). Even though we cannot directly resolve the concentration mass scale M conc discussed around Fig. 1, a visual extrapolation of the inferred concentration-mass relation c ( M , z ) to higher masses suggests that the curves converge around M conc ≈ 4 . 7 × 10 14 M ⊙ / h . This implies that M conc ≈ 20 -50 M 0 represents a characteristic scale, independent of the axion fraction, at which the c ( M , z ) relation in MDM diverges from that in CDM.', '4.1 Halo Model Details': "Having established that identifying halos based solely on the CDM component produces reliable HMFs and c ( M , z ) relations, we now focus on calibrating the halo model framework A xion HM code (Vogt  \net al. 2023) to the MDM simulations and compare the resulting nonlinear power spectra. A xion HM code , an adaptation of HMC ode -2020 (Mead et al. 2021), incorporates ultralight particles as part of the DM. It builds on the halo model originally implemented in Dentler et al. (2022) to constrain axions using weak lensing shear statistics. The code takes as input a range of cosmological and DM parameters: Ω m, Ω b, f , m , H 0, ns , As , k piv. Here, the axion density parameter is Ω a ≡ Ω FDM = f Ω m, while the combined density of CDM and baryons ('cold' matter) is Ω c = Ω m -Ω a = Ω CDM + Ω b. Using this input, A xion HM code calculates the non-linear power spectrum at a desired redshift z . The linear FDM power spectrum, which is integral to the model, is calculated using the public version of A x -ion C amb . The total matter overdensity in the MDM cosmology can be decomposed into a sum of the cold matter, δ c, and axions, δ a,  \nδ m = Ω c Ω m δ c + Ω a Ω m δ a . (44)  \nThe non-linear power spectrum in A xion HM code is constructed as  \nP ( k ) = Ω c Ω m ! 2 P c( k ) + 2 Ω c Ω a Ω 2 m P c,a( k ) + Ω a Ω m ! 2 P a( k ) , (45)  \nwhere P c, P c,a( k ) ∝ δ c δ a and P a are the cold, cross and axion power spectrum, respectively. For P c( k ), the standard halo model (see Sec. 4.4 for improvements on top of this) is adopted, splitting the power spectrum into the one-halo term P 1h and two-halo term P 2h ,  \nP c( k ) = P 1h c ( k ) + P 2h c ( k ) . (46)  \nFor axions we adopt the biased tracer formalism (Massara et al. 2014) which assumes that a sub-component, δ L, cannot cluster and evolves approximately linearly while the remaining fraction  \nF h = 1 ¯ ρ a Z ∞ M cut d M c n ( M c) b ( M c) M a( M c) ∈ [0 , 1] (47)  \nis in halos, i. e.  \nδ a = F h δ h + (1 -F h) δ L . (48)  \nHere, n ( M c) denotes the cold HMF, and we have introduced the cold halo bias, b ( M c), as well as the axion halo mass to cold halo mass relation, M a( M c). The biased tracer formalism assumes that axion halos only form in and around cold matter halos and thus the halo mass function for axions is the same as for the cold field, n ( M a)d M a = n ( M c)d M c, and the linear axion halo bias corresponds to the cold halo bias, b ( M a) = b ( M c), i. e. M a is itself a function of M c, a relation we need to specify (see below). In total there are three new quantities we have to provide to complete the MDM halo model: the cut-o ff mass, M cut, the axion halo mass relation, M a( M c), and the axion halo density profile ρ a( r , M a , z ). For details see Vogt et al. (2023).  \nTo construct the cold halo density profile, we retain the Bullock et al. (2001)-inspired model for the concentration parameter of cold halos of mass M c,  \nc c( M c , z ) = B 1 + z f( M c , z ) 1 + z ! , (49)  \nwhere z f denotes the formation redshift, which is defined by  \nD + ( z f) D + ( z ) σ c(0 . 01 M c , z ) = δ crit , (50)  \nwith D + ( z ) denoting the total MDM linear growth rate (depending only on cosmology and redshift). The minimum halo concentration B = 5 . 196 is attained when the solution to Eq. (50) yields z f < z .  \nFigure 7. Spherically averaged density profile of the FDM component in the f = 0 . 1 MDM cosmology at redshift z = 1. The top panel covers a total halo mass range of M tot = 10 9 -10 13 M ⊙ / h while the lower panel focuses on the range M tot = 10 10 . 4 -10 12 . 4 M ⊙ / h with a finer mass resolution of ∆ log( M tot) = 0 . 2. Shaded areas delineate the standard error of the median, and solid curves trace the median profile in each mass bin. Arrows at the bottom of the panels indicate the FDM grid resolution scale of ∆ x / 2 = 14 . 6 ckpc / h , rescaled by the mean virial radius in each mass bin. For masses M tot > M cut = 1 . 7 × 10 12 M ⊙ / h , the axion profile displays a distinct core-like feature, indicating the presence of a soliton. The steep decline in central densities between neighboring mass bins becomes pronounced for M tot < M J = 1 . 0 × 10 11 M ⊙ / h , reflecting a diminished influence of the FDM component on structure formation.  \n<!-- image -->  \nEven though the Bullock et al. (2001) model has been shown to predict a sharp decline in the concentration at high mass that is inconsistent with N -body simulations and to insu ffi ciently capture e ff ects in non-CDM cosmologies such as WDM, the more accurate Ludlow et al. (2016) approach introduced in Sec. 3.2 is less suitable for the halo model since it typically underpredicts CDM small-scale power around k ≳ 5 cMpc -1 in both pure CDM and MDM simulations. Shifting to Ludlow et al. (2016) would also necessitate recalibrating HMC ode -2020 parameters (most notably the halo bloating parameter η ) without gaining in precision.", '4.2 Cut-O ff Mass': "We first perform a sanity check of the cut-o ff mass M cut adopted by A xion HM code . It is obtained by invoking the concept of the halo  \nFigure 8. Axion halo mass-cold halo mass relation M a( M c) in MDM cosmologies at various axion fractions f = 0 . 0 -0 . 3 across redshifts z = 1 -4. The mass resolution is ∆ log( M tot) = 0 . 2 as in Fig. 7 and we estimate the axion halo mass M a = 4 π R R vir 0 d r r 2 ρ FDM( r ) by integrating the axion density profile out to the virial radius R vir. For consistency we estimate the cold halo mass M c likewise via integration since the direct estimate M c = Ω c / Ω m M tot from the R ockstar virial mass can deviate therefrom at the high-mass end due to resolution e ff ects. Dashed curves indicate the cosmic mean relation M a = Ω a Ω c M c, currently implemented in A xion HM code . Colored arrows at the bottom of the panels indicate the cut-o ff mass M cut = 5 × 10 11 -10 12 M ⊙ / h while the black arrow denotes the Jeans mass M J. Note that we infer a steeper M a( M c) relation than the cosmic average for M c ≲ M J, with steepness increasing toward lower f . The transition range widens toward high axion fractions, covering M J -M cut in some cases. We fit a broken power law, Eq. (52), for M a > 10 7 M ⊙ / h (grey arrow) and show results as dotted curves.  \n<!-- image -->  \n/circledot  \n/circledot  \nJeans scale r hJ from Sec. 2.2. To verify the accuracy of Eq. (15) as the length scale below which axions fail to cluster inside DM halos, we calculate density profiles of the axion component in the MDM simulations. Fig. 7 illustrates these profiles for the f = 0 . 1 MDM cosmology at redshift z = 1 across various mass bins. Our analysis, which includes simulations across axion fractions f = 0 . 0 -0 . 3 and redshifts up to z = 10, consistently shows that halos with mass M tot ≳ M cut exhibit a distinct soliton core. Below the cut-o ff mass, and particularly below the Jeans mass M J, the central density of the core declines sharply, resulting in featureless axion profiles for halos with M tot ≪ M J.  \nThis steep decline in profile density supports M cut as an e ff ective estimate of the mass scale below which halos (slowly) transition to being predominantly CDM-dominated in MDM. It is important to note that in the transitional mass range M J -M cut, 'quantum pressure' continues to play a role but is insu ffi cient to fully counteract gravitational e ff ects (recall definition of M J in Sec. 2.2). The transition to CDM-dominated halos can occur over a large mass range exceed-  \ning 0.5 dex, and is not always complete for mass M tot ≈ M J. In the following, we study integrated axion density profiles (i. e. the axion halo mass relation) and we will see that the Jeans mass and cut-o ff mass can help describe the transition toward CDM-dominated halos.", '4.3 Axion Halo Mass Relation': 'We now turn to the relationship between the axion halo mass, M a, and the cold halo mass, M c. In the public version of A xion HM code , it is assumed that the axion mass M a follows the cosmic abundance fraction relative to the cold halo mass M c down to the cut-o ff mass M cut, i. e. that M a = ( Ω a / Ω c) M c for M c > M cut. However, as seen in Sec. 4.2, there are strong indications that this simplistic approach may not accurately capture the axion halo mass-cold halo mass relation found in simulations. Fig. 8 shows the M a( M c) relation inferred from our MDM simulations at various axion fractions f = 0 . 0 -0 . 3 across redshifts z = 1 -4. We estimate both the axion and cold halo  \nTable 4. Best-fit values of β 2 in Eq. (52) for our MDM simulation suite across a range of redshifts z = 1 -8. In the fitting process, we fix β 1 at 1, as allowing β 1 to vary only results in minor improvements. Note that β 2 increases with higher redshift and lower values of f .  \nmass by integrating the respective density profile out to the virial radius R vir,  \nM a,c = 4 π Z R vir 0 d rr 2 ρ a,c( r ) . (51)  \nThe cold halo mass could also be obtained directly from the R ock -star virial mass, M c = ( Ω c / Ω m) M tot, but the two estimates can deviate by more than 0 . 5 dex (at the high-mass end) due to resolution e ff ects, hence for consistency we determine M c via integration.  \nAs shown in Fig. 8, the cosmic average relation, M a = ( Ω a / Ω c) M c, is a valid approximation at the high-mass end. 6 Toward lower mass, the simulated M a( M c) relation becomes significantly steeper. Deviations from the cosmic mean exceed 1 dex in several mass bins. These deviations become particularly pronounced below the Jeans mass, M J, consistent with the axion density profiles discussed in Sec. 4.2. The break can span a wide mass range, especially at higher axion fractions, typically encompassing M J -M cut and exceeding 0.5 dex for high values of f . To better capture this behaviour, we propose to parametrise the axion halo mass-cold halo mass relation as a broken power law:  \nM a =       1 + M c M J ! -β 1       -β 2 Ω a Ω c M c , (52)  \nwhere the steepness of the suppression is controlled by the parameter β 2 while the sharpness of the transition at M ≈ M J is controlled by β 1, similar to Eq. (39). The linear Jeans mass is given by Eq. (11). Letting β 1 vary during the fit along with β 2 leads to only minor improvements, hence we fix β 1 = 1. Best-fit values of β 2 for M a > 10 7 M ⊙ / h are presented in Table 4, and the corresponding best-fit curves are illustrated in Fig. 8. We observe a weak dependence of β 2 on redshift, but a notable decrease of β 2 as the axion fraction f increases, indicating a weaker suppression. This analysis shows that MDM cosmologies dominated by CDM f < 0 . 5 are different from pure axion cosmologies ( f = f max), where axion halos do not exist below M J and virialised axion halos only form above M cut (see Sec. 2.2).  \nAnother implication of the M a( M c) relation deviating from the cosmic mean, where M a = ( Ω a / Ω c) M c, is that expressing total DM density profiles as ρ DM( r ) = f ρ FDM( r ) + (1 -f ) ρ CDM( r ) (see e. g. Shevchuk et al. 2023) can at best be a good approximation at  \nhigh masses above the cut-o ff mass M cut. A more detailed analysis of axion density profiles, including their parametrisation via soliton + NFWprofiles, their (non-)formation at low axion fractions (see e. g. Schwabe et al. 2020), and their use in axion forecasts and constraints, is beyond the scope of this paper and will be addressed in future work.', '4.4 A xion HM code Parameters': 'The updated version of A xion HM code has several improvements over the original implementation by Vogt et al. (2023), which we now summarise. As in Vogt et al. (2023), we adopt the HMC ode -2020 parameters (Mead et al. 2021), which have been introduced to improve the model in its fit to Λ CDM simulations over the standard halo model. The HMC ode -2020 parameters were calibrated using the M ira T itan matter power spectrum emulator of Heitmann et al. (2016); Lawrence et al. (2017). This cosmic emulator encompasses eight cosmological parameters and provides an accuracy of 4 % for k < 7 h cMpc -1 . In turn, the accuracy of HMC ode -2020 when compared to simulated Λ CDMdata is excellent with a RMS error of less than 2.5 % for k < 10 h cMpc -1 and z < 2 (Mead et al. 2021).  \nThe HMC ode -2020 parameters are only calibrated up to z = 2 and it is a priori not guaranteed that a non-linear power spectrum with these parameters at z > 2 is more accurate than without the parameters (i. e. standard halo model). We adopt the HMC ode -2020 parameters up to z = 3 . 5 and will show in Sec. 4.5 that not only is the agreement with simulations within the 10 % margin for pure Λ CDM, but that the improvements to A xion HM code (including generalising and recalibrating one of the HMC ode -2020 parameters) yield good agreement with simulated non-linear power spectra up to at least z = 3 . 5 for axion fractions that we can assess, f < 0 . 3.  \nOne of the HMC ode -2020 parameters is the one-halo term damping. In the standard halo model approach (see Eq. 46), the one-halo term is typically constant on large scales. However, this does not accurately reflect mass and momentum conservation. It was demonstrated by Smith et al. (2003) that the one-halo term should increase as P 1h( k ) ∝ k 4 at small k (i. e., it should dampen compared to a constant at small k ). To address this, HMC ode -2020 implements a modification:  \nP 1h c ( k ) → P 1h c ( k ) ( k / k ∗ ) 4 1 + ( k / k ∗ ) 4 . (53)  \nThis adjustment ensures that the one-halo term grows as expected and is suppressed on large scales. Consequently, on large scales, the non-linear power spectrum is primarily determined by the two-halo term, which aligns with the (perturbed) linear power spectrum. The suppression e ff ect is controlled by the free parameter k ∗ , which was fitted as:  \nk ∗ = 0 . 05618 × σ 8 , c( z ) -1 . 013 h cMpc -1 . (54)  \nWe can now list the improvements made in the new version of A xion HM code :  \n- · To construct the two-halo term for the cold matter component, we apply a perturbative damping to the linear power spectrum, so that  \nP 2h c ( k ) = P L c ( k ) 1 -q ( k / k d) n d 1 + ( k / k d) n d ! . (55)  \nThis formulation aligns with the recommendations of Mead et al. (2021), who suggest that, given the precision of modern halo models, it is crucial to account for the largest-scale non-linear e ff ects.  \nTable 5. A xion HM code parameters, built on top of the standard halo model. For each parameter, we provide its description, the equation defining it, the default value in the standard halo model, the fitted functional form or value, and an example of this function evaluated at z = 1 for a standard MDM cosmology with f = 0 . 1, m = 3 . 16 × 10 -25 eV and Planck Collaboration et al. (2016) cosmological parameters.  \nSpecifically, perturbation theory indicates that the most significant non-linear e ff ect on large scales is a small damping of power. The term introduced in Eq. (55) is designed to capture this e ff ect, with the three parameters fitted and given by Mead et al. (2021),  \nk d = 0 . 05699 × σ 8 , c( z ) -1 . 089 h cMpc -1 , q = 0 . 2696 × σ 8 , c( z ) 0 . 9403 , (56) n d = 2 . 853 .  \nNote that the di ff erence between perturbatively corrected linear theory and the standard two-halo term is tiny but significant gains in computational time are made when replacing the integral expression of the standard two-halo term by Eq. (55).  \n- · We achieve additional significant speed-ups through modularisation and in-memory storage of arrays. For example, the formation redshift z f for a given cold halo mass M c is computed once and subsequently stored. As a result, we reduce the runtime on a single-core machine to under 1 minute for a single evaluation of A xion HM code , with only a minor increase in memory overhead.\n- · We fix a bug in the implementation of the halo bloating e ff ect mediated by the parameter η which scales the wavenumber k in the Fourier transformation of the NFW profile of cold halos as  \n˜ u ( k , M , z ) → ˜ u ( ν η k , M , z ) . (57)  \nThe halo bloating parameter was fitted to (Mead et al. 2021)  \nη = 0 . 1281 × σ 8 , c( z ) -0 . 3644 , (58)  \nwhich we also adopt in A xion HM code .  \n- · We adjust the cosmic axion halo mass-cold halo mass relation M a( M c) = ( Ω a / Ω c) M c by incorporating a broken power law below the linear Jeans mass M J, in accordance with the findings from our MDM simulations which reveal a strong decrement in the inferred axion halo mass M a compared to the cosmic mean relation (see Eq. 52). Consequently, lower integration limits for various quantities, such as the clustered fraction of Eq. (47), are adjusted from M cut down to 0.\n- · Wecontinue to use the smoothing parameter α , which allows us to overcome the simplistic assumption of a purely additive behaviour  \nof one- and two-halo terms by modelling  \nP c( k ) = ( P 1h c ( k ) α + P 2h c ( k ) α ) 1 /α . (59)  \nMead et al. (2021) fitted parameter α to a general form,  \nα = 1 . 875 × (1 . 603) n e ff ,c ( z ) , (60)  \nwhere the e ff ective spectral index at the non-linear length scale is  \nn e ff ,c( z ) = -d ln( σ 2 c ( R , z )) d ln( R ) GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> σ c = δ crit -3 . (61)  \nWe first generalise the smoothing parameter to the cold matter component of MDM. To prevent strong smoothing e ff ects to spill over from the quasi-linear regime into the large-scale regime, the coldcold smoothing is modeled using a logistic function between α 1 , c and α 2 , c:  \nα c( k , z ) = α 2 , c( z ) + α 1 , c( z ) -α 2 , c( z ) 1 + exp(( k piv -k ) / ∆ k ) . (62)  \nThe width of the transition is controlled by ∆ k = 0 . 1 h cMpc -1 . In the small-scale limit, we use the following expression:  \nα 1 , c( z ) = 1 . 875 × (1 . 603) n e ff ,c ( z ) 1 + a GLYPH<16> 10 -24 m GLYPH<17> b GLYPH<16> Ω a Ω m GLYPH<17> d (1 + z ) e . (63)  \nThis parametrisation reduces to the default smoothing parameter for Λ CDM in the limit of a high axion mass m → ∞ or a low axion fraction Ω a → 0, given that b , d > 0. Due to our lack of MDM simulation for axion mass values other than m = 3 . 16 × 10 -25 eV, we fix the exponent b = 0 . 0450 based on the expectation of a weak dependence on m . In the large-scale limit, where k < k piv = 0 . 5 h cMpc -1 , the smoothing is well fit by  \nα 2 , c( z ) =        α 1 , c( z ) for f < 0 . 01 max GLYPH<16> α 1 , c( z ) , 1 . 10 1 + z GLYPH<17> for f ≥ 0 . 01 . (64)  \nWe fit parameters a , d , e against the simulated cold-cold power spectrum P c (see Eq. 59) using the Nelder & Mead (1965) simplex algorithm with equal logarithmic weights in k ∈ [0 . 2 -30] h cMpc -1 and  \nFig. 10 displays the non-linear matter power spectrum predictions from the updated version of A xion HM code at redshifts z = 1 -3 . 5. We use the default configuration, which includes large-scale damping of the one-halo term, perturbative damping of the two-halo term, halo bloating, and transition smoothing between one- and twohalo terms in the quasi-linear regime. These predictions are compared with MDM power spectra derived from our simulations using P ylians (Villaescusa-Navarro 2018).  \n<!-- image -->  \n2  \nFigure 9. A xion HM code vs MDM simulations. We compare the non-linear power spectrum for our MDM simulation (solid) with the halo model implementation using the new version of A xion HM code (dashed). Results are shown for f = 0 . 1 and m = 3 . 16 × 10 -25 eV at redshifts z = 1. We show the full power spectrum (blue) and the three contributing terms from Eq. (45) (green, orange, and red). For comparison, we show the total matter power spectrum in a pure Λ CDM cosmology (black solid). The dashed vertical line indicates the fundamental frequency of the simulation box k f = 2 π/ L box while the Nyquist frequency k Ny = π N 1 / 3 CDM / L box is shown as a vertical dotted line.  \nusing data from z = 1 , 2 , 3 and f = 0 . 0 , 0 . 01 , 0 . 1 , 0 . 2 , 0 . 3 with equal weight, minimising the figure-of-merit  \nψ 2 = 1 N X z , f , k | log(AxHM( k , z , f )) -log(Sim( k , z , f )) | 2 . (65)  \nIn addition, we find better agreement with simulations if we add an independent smoothing parameter to the cold-axion cross power spectrum P c,a (which depends on P c):  \nP c,a( k ) = F h GLYPH<16> P 1h c,a ( k ) α c,a + P 2h c,a ( k ) α c,a GLYPH<17> 1 /α c,a + (1 -F h) q P c( k ) P L a ( k ) (66)  \nWe parametrise α c,a in the exact same way as the small-scale coldcold smoothing term α 1 , c( z ) in Eq. (63), and fit parameters a , d , e against the simulated P c,a spectrum. Best-fit values for α 1 , c , α 2 , c and α c,a, along with all other A xion HM code parameters, are provided in Table 5.  \n· In the new version of A xion HM code , we activate the halo bloating parameter η (for the cold halo density profile), along with the one-halo damping (see Eq. 53), the two-halo damping (see Eq. 55), and the smoothing parameter (see Eqs. 63 and 62) by default. This configuration significantly enhances the agreement with MDM simulations, both qualitatively and quantitatively. These improvements validate the necessity of HMC ode -2020-like parameters for accurate MDMsimulations, as initially suggested by Vogt et al. (2023).', '4.5 Non-Linear Power Spectra': "We now validate the updated version of A xion HM code . Recall that the M a( M c) relation was fitted over the redshift range z = 1 -8, the smoothing parameter was fitted for z = 1 -3, and the remaining HMC ode -2020 parameters are calibrated up to z = 2. As a re-  \nsult, A xion HM code should be used with caution outside the range z ≈ 1 -3 and f ≈ 0 -0 . 3. In Fig. 9, we illustrate its predictions for the non-linear power spectrum, compared to results from our MDM simulations for f = 0 . 1 at redshift z = 1. We find very good agreement between the total predicted and simulated matter power spectrum (blue solid and dashed). The moderate suppression of power relative to Λ CDM on scales k > 1 h cMpc -1 is well captured. The agreement between A xion HM code and simulations extends to the cold matter power spectrum (green), though there is a slight discrepancy for the cross and FDM-FDM power spectra (orange and red) at k > 10 h cMpc -1 . Note that beyond k > 10 h cMpc -1 we expect strong e ff ects from baryonic physics (see e.g. Mead et al. 2021), which are not accounted for in either our MDM simulations or A x -ion HM code .  \nA xion HM code also captures the enhanced non-linear power on scales k ≈ (1 -10) h cMpc -1 relative to Λ CDM. The e ff ect has been reported by Vogt et al. (2023) who traced it to an enhancement in the cross and FDM-FDM power spectra, which in turn is caused by the coherence of the soliton in the axion halo density profile, increasing the correlation function of the axion field on such scales. The enhancement is stronger for axion masses around m ≈ 10 -22 eV (shifting toward higher k ) and was also observed in pure FDM cosmologies when accounting for 'quantum pressure' (Nori et al. 2018; May & Springel 2023).  \nAs anticipated, A xion HM code demonstrates excellent agreement (deviations less than 10 %) with simulations for pure CDM ( f = 0 . 0) across all scales and redshifts examined. For an axion fraction of f = 0 . 1, the code maintains high accuracy on both large (small k ) and small scales (large k ) at all considered redshifts, with maximum deviations also remaining below 10 %. When the axion fraction is increased to f = 0 . 2, maximum deviations rise to approximately 20 % for scales with k < 10 h cMpc -1 ; however, at redshifts around z ≈ 1, deviations on these scales remain under 10 %. For f = 0 . 3, the maximum deviations increase to about 30 % for scales with k < 10 h cMpc -1 , though at redshifts around z ≈ 1, deviations do not exceed approximately 20 % on these scales.  \nWe are unable to assess the agreement with MDM simulations for redshifts z < 1 due to the potential for unresolved high velocity dispersions in the FDM solver (see Sec. 2.6). Similarly, neither HMC ode -2020 nor A xion HM code are anticipated to perform well at redshifts above z ≈ 3 . 5, which defines the upper limit of our redshift window. Note that our MDM simulation data is limited to FDM fractions f ≤ 0 . 3. However, this limitation is not critical, as ultralight axions are permitted to exist in substantial portions within the FDM window but not at extremely high fractions close to f max, which are increasingly ruled out with greater significance beyond the FDM window. Additionally, our fitting of A xion HMcode parameters is constrained to Λ CDM and MDM cosmologies with axion mass m = 3 . 16 × 10 -25 eV, which lies within the FDM window 10 -25 eV ≲ m ≲ 10 -23 eV. Due to the lack of MDM simulation data for other axion mass values, we cannot validate the accuracy of A xion HM code for axion masses significantly di ff erent from m = 3 . 16 × 10 -25 eV, except in the limit as m → ∞ . The exponent value b = 0 . 0450 in Eq. (63) is thus an educated guess based on the  \nFigure 10. Comparison of non-linear matter power spectrum predictions from the updated version of A xion HM code and MDM simulations at across redshifts z = 1 -3 . 5. Solid lines represent the new model, while dotted lines show the previous version from Vogt et al. (2023). Horizontal dashed and dotted lines indicate 10 % and 20 % deviations, respectively. Arrows at the bottom mark comoving Jeans scale k J at each redshift (see Eq. 10). The axion fraction f increases from the top left panel ( f = 0 . 0) to the bottom right ( f = 0 . 3) along with a decrease in overall accuracy of A xion HM code .  \n<!-- image -->  \nexpectation of a weak dependence on m . Future MDM simulations will provide the opportunity to refine this estimate and validate the parameter b .", '5 CONCLUSIONS': 'In view of tighter constraints being put on pure, single-field ultralight axion DM (see Dome et al. 2022, for a compilation), it is promising to relax the requirement that ultralight axions must comprise all of the DM in the Universe. In this work, we focused on the FDM window 10 -25 eV ≲ m ≲ 10 -23 eV in which ultralight axions are allowed to exist in large portions (albeit not f = f max) (although see Shevchuk et al. 2023; Lazare et al. 2024; Winch et al. 2024).  \nMethods: We implemented an MDM gravity solver (see Sec. 2.5) and ran state-of-the-art simulations of mixed ultralight axion cosmologies dominated by CDM ( f < 0 . 5). Our MDM simulations were designed to capture the wave dynamics across small and intermediate length scales, with particular emphasis on achieving numerical convergence and resolution thresholds. ICs were set up carefully via second order Lagrangian perturbation theory (2LPT) using M u -sic (Hahn & Abel 2011) based on matter and velocity power spectra  \ncalculated using A xion C amb (Hlozek et al. 2015). By rigorously enforcing criteria such as velocity resolution and resolving the axion half-mode scale k 1 / 2, while also ensuring accurate representation of halo populations, our simulations faithfully reproduce internal halo structures above redshifts z ≈ 1, providing a robust platform for evaluating common (semi-)analytical techniques such as halo models.  \nHalo Mass Distribution: We identified halos using the R ockstar particle-based halo finder applied on the (equal-mass) CDM distribution and report total halo mass distributions ( M tot = M c + M a) in MDM. We found good agreement between the Sheth-Tormen model based on the linear MDM matter power spectrum and the inferred HMF across a very wide range of redshifts z = 1 -10 and axion fractions f = 0 . 0 -0 . 3, justifying the usage of R ockstar on the CDM component a posteriori. The HMF in MDM cosmologies branches o ff from the CDM one at the characteristic mass M 0 = 1 . 6 × 10 10 ( m / 10 -22 eV) -4 / 3 M ⊙ , but instead of plateauing toward low mass as in WDM or turning over as in pure FDM, the HMFcontinues to increase in a power-law fashion as M tot decreases. By providing best-fit results to a one-parameter model of the MDM HMF, Eq. (39), we hope to facilitate parameter sampling across MDMcosmologies in Bayesian constraint and forecast analyses.  \nDensity profiles: We fit the Einasto model against total DM density profiles and found the Einasto parametrisation to be reliable across the entire range of FDM fractions f = 0 . 0 -0 . 3 and redshifts z = 1 -10 studied. The resulting median concentration-mass relation c ( M tot) is in good agreement with the Ludlow et al. (2016) model based on extended Press-Schechter theory, suggesting that instead of turning over at around two decades above the half-mode mass, 100 × M 1 / 2, as in pure WDM and FDM, the concentration exhibits a decrease before recovering and increasing toward smaller halo mass M tot. This results in an e ff ective flat U-shaped c ( M tot) relation for high values of the axion fraction f , and is in agreement with insights from the halo mass distribution.  \nCalibrating AxionHMcode: We aimed at improving the calibration of the halo model code A xion HM code based on a biased tracer approach using insights from our MDM simulations. The aforementioned success of reproducing analytical total halo mass distributions and concentration-mass relations based on halos identified solely using the CDM component lends additional credence to the viability of the biased tracer approach. The modifications we introduce (apart from minor ones) are threefold: First, we model the axion halo mass-cold halo mass relation M a( M c) as a broken power law below the Jeans mass M J and retain the cosmic mean relation M a = ( Ω a / Ω c) M c above M J. Second, we generalise the transition smoothing parameter α to MDM with a dependence on m and Ω a while heeding the spill-over e ff ects of strong smoothing on the largescale regime using a logistic function for the wavenumber-dependent smoothing α c( k ). Third, we introduce various speed-ups by making sure numerical functions are not evaluated too often, leading to a slight increase in memory requirements while reducing run-time on a single-core machine to below 1 minute for a single evaluation of A xion HM code . The code exhibits excellent agreement with simulations for pure Λ CDM, with deviations under 10 % on scales below k < 20 h cMpc -1 and redshifts z = 1 -3 . 5. For axion fractions f ≤ 0 . 3, the model maintains accuracy with deviations under 20 % at redshifts z ≈ 1 and scales k < 10 h cMpc -1 , though deviations can reach up to 30 % for higher redshifts when f = 0 . 3.  \nOutlook: Mixed ultralight axion cosmologies dominated by CDM ( f < 0 . 5) have their clustering properties determined to first order by the CDM component. More precisely, the peak height distribution is largely shaped by the cold component as first suggested by Massara et al. (2014), motivating not only the very assumptions underlying A xion HM code but also the identification of halos based solely on the CDM component. MDM models might have the potential to reconcile observational constraints while providing a semiphenomenological route to understanding the nature of DM. In upcoming forecast and constraint analyses based on observational data, having good control of non-linear predictions in MDM models is key, and this work aims to contribute to that undertaking.', '6 ACKNOWLEDGEMENTS': 'We acknowledge useful discussions with Sophie Vogt. TD acknowledges support from the Isaac Newton Studentship and the UK Research and Innovation (UKRI) Science and Technology Facilities Council (STFC) under Grant No. ST / V50659X / 1. SM acknowledges support by the National Science Foundation under Grant No. 2108931. SB is supported by the UKRI Future Leaders Fellowship (Grant No. MR / V023381 / 1). AL acknowledges support from NASA grant 21-ATP21-0145. DJEM is supported by an Ernest Rutherford Fellowship from the STFC, Grant No. ST / T004037 / 1. This work used the DiRAC@Durham facility managed by the Insti-  \ntute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility, with equipment funded by BEIS capital funding via STFC capital grants ST / K00042X / 1, ST / P002293 / 1, ST / R002371 / 1 and ST / S002502 / 1, Durham University and STFC operations grant ST / R000832 / 1.', '7 DATA AVAILABILITY': 'MDMsnapshot data and post-processing scripts are made available upon reasonable request.', 'REFERENCES': 'Heitmann K., et al., 2016, ApJ, 820, 108'}