2025ApJ...986..127N

We present 72 cosmological dark-matter–only N-body zoom-in simulations with initial conditions beyond cold, collisionless dark matter (CDM), as the first installment of the COZMIC suite. We simulate Milky Way (MW) analogs with linear matter power spectra P(k) for (i) thermal-relic warm dark matter (WDM) with masses m_WDM ∈ [3, 4, 5, 6, 6.5, 10]keV, (ii) fuzzy dark matter (FDM) with masses m_FDM ∈ [25.9,69.4, 113, 151, 185, 490] × 10^–22eV, and (iii) interacting dark matter (IDM) with a velocity-dependent elastic proton scattering cross section σ = σ₀v^n^, relative particle velocity scaling n ∈ [2, 4], and dark matter mass m_IDM ∈ [10^–4, 10^–2, 1] GeV. Subhalo mass function (SHMF) suppression is significantly steeper in FDM versus WDM, while dark acoustic oscillations in P(k) can reduce SHMF suppression for IDM. We fit SHMF models to our simulation results and derive new bounds on WDM and FDM from the MW satellite population, obtaining m_WDM > 5.9 keV and m_FDM > 1.4 × 10^–20 eV at 95% confidence; these limits are ≈10% weaker and 5× stronger than previous constraints owing to the updated transfer functions and SHMF models, respectively. We estimate IDM bounds for n = 2 (n = 4) and obtain σ₀ < 1.0 × 10^–27 cm², 1.3 × 10^–24 cm², and 3.1 × 10^–23 cm² (σ₀ < 9.9 × 10^–27 cm², 9.8 × 10^–21 cm², and 2.1 × 10^–17 cm²) for m_IDM = 10^–4, 10^–2, and 1 GeV, respectively. Thus, future development of IDM SHMF models can improve IDM cross section bounds by up to a factor of ∼20 with current data. COZMIC presents an important step toward accurate small-scale structure modeling in beyond-CDM cosmologies, critical to upcoming observational searches for dark matter physics.