© The Authors 2026

1. Introduction

The Universe, with its galaxies organized into clusters and filaments, is the result of approximately 13.8 billion years of evolution. Dark energy and dark matter (DM) shape spacetime and weave a cosmic web that connects astrophysical objects through structures governed by gravitational effects. The nature and existence of DM remain among the most significant unresolved questions in astrophysics. According to the prevailing Lambda cold dark matter (ΛCDM) model, DM consists of cold–moving at velocities much lower than the speed of light–and collisionless particles (e.g., Blumenthal et al. 1984; Davis et al. 1985; Frenk et al. 1988). The ΛCDM model, despite its simplicity, has achieved significant success. For example, it explains the power spectrum of the cosmic microwave background and accurately reproduces the large-scale distribution of matter throughout the Universe. When it comes to smaller scales, such as galaxies and subgalactic scales, the ΛCDM model encounters several challenges, including the core-cusp (Moore 1994; de Blok & McGaugh 1997), the too-big-to-fail (Boylan-Kolchin et al. 2011), the missing satellites (Moore et al. 1999; Klypin et al. 1999), and the diversity (Oman et al. 2015) problems. Recent advancements in simulations incorporating hydrodynamics and galaxy formation physics (Vogelsberger et al. 2020) have addressed some of these discrepancies within the ΛCDM paradigm, such as the core-cusp and the missing satellites problems (Pontzen & Governato 2012; Governato et al. 2012; Vogelsberger et al. 2012, 2014; Brooks et al. 2013; Chan et al. 2015; Wetzel et al. 2016; Engler et al. 2021). Moreover, baryonic feedback and tidal effects can alleviate the too-big-to-fail problem in satellite galaxies (Zolotov et al. 2012; Brooks & Zolotov 2014). However, in isolated dwarfs, where such processes are less effective, the too-big-to-fail and diversity problems remain significant challenges for the ΛCDM model (Papastergis et al. 2015; Callingham et al. 2020). Warm dark matter (WDM; e.g., Bode et al. 2001) stands out as an alternative to CDM, proposing particles with non-negligible peculiar velocities. Recent hydrodynamic simulations have shown that WDM models can delay star formation in low-mass haloes and lead to observable differences in stellar ages and central densities of dwarf galaxies (Herpich et al. 2014; Macciò et al. 2019; Shen et al. 2024). Other compelling alternative candidates include self-interacting dark matter (e.g., Spergel & Steinhardt 2000; Burkert 2000; Zavala et al. 2013), where particles interact significantly with one another, and scalar field or fuzzy dark matter (e.g., Hui et al. 2017; Zhang et al. 2019; Niemeyer 2020), which involve ultralight bosons.

In this work, we focus on WDM models by analyzing cosmological simulations from the DREAMS project (DaRk mattEr with AI and siMulationS, Rose et al. 2025). The DREAMS project integrates astrophysics, particle physics, and machine learning to investigate the nature of DM, with a particular emphasis on small-scale physics in the Universe. The simulations model our Universe across various volumes and resolutions, ranging from large-scale uniform-boxes and Milky Way (MW) zoom-ins to dwarfs zoom-in configurations, tailored for specific DM models. Despite the existence of various experimental constraints on the lower mass limit of the WDM particle – up to 9.7 keV, as indicated by analyses of gravitationally lensed systems and MW satellites (Nadler et al. 2021) – WDM remains a crucial benchmark for testing alternative DM models in cosmological simulations. In fact, WDM models represent a simple alternative to CDM and serve as a natural starting point for exploring scenarios beyond the ΛCDM paradigm, as demonstrated in the DREAMS project.

In this study, we focus on key galactic scaling relations – specifically those linking galaxy size, DM content, and stellar mass – and investigate how varying WDM particle masses affect these relations, with particular attention given to low-mass galaxies. Indeed, galactic scaling relations arise from the fundamental physics governing galaxy formation and evolution, making them inherently sensitive to the nature of DM. Our goal is to evaluate how well WDM models can be constrained through these relations. In addition to the nature of DM, astrophysical feedback processes – such as those from supernovae (SNe), active galactic nuclei (AGNe), or gravitational heating – can significantly modulate galactic scaling relations (Dekel & Birnboim 2006; Governato et al. 2007; Moster et al. 2010; Tortora et al. 2019, 2025; Busillo et al. 2023, 2025; Vogelsberger et al. 2013; Torrey et al. 2014; Weinberger et al. 2016; Pillepich et al. 2018; Roca-Fàbrega et al. 2021; Ni et al. 2023). To this end, we analyzed simulations at different resolutions and volumes, including both large-scale uniform-box simulations and high-resolution MW zoom-in simulations, and we compared the scaling relations obtained from these simulations with observational trends inferred from nearby dwarf galaxies in the Local Volume Database (LVDB; Pace 2025) and from the Spitzer Photometry and Accurate Rotation Curves (SPARC; Lelli et al. 2016).

This paper is the fourth in the CASCO project: Cosmological and Astrophysical parameters from Cosmological simulations and Observations. Busillo et al. (2023, hereafter Paper I) and Busillo et al. (2025, hereafter Paper II) used the CAMELS simulations within a ΛCDM cosmology to constrain cosmological and astrophysical parameters – such as SN and AGN feedback – by comparing simulated and observed scaling relations, considering both early- and late-type galaxies. Tortora et al. (2025, Paper III) focused on the emergence of the so-called golden mass, a characteristic stellar mass scale where star formation efficiency peaks, and studied how its value evolves across cosmic time. In this work (Paper IV), we extend the CASCO framework to a WDM cosmology by using DREAMS simulations to investigate how feedback and the nature of DM shape the internal structure of galaxies.

In Section 2, we present the DREAMS simulations; the observational datasets used for comparison, namely the LVDB and SPARC catalogs; and the fitting procedure adopted to constrain cosmological, astrophysical, and WDM model parameters based on the agreement between simulations and observations. In Section 3, we discuss our results regarding the impact of these parameters on galaxy scaling relations and galaxy counts, and we assess the effectiveness and limitations of the fitting approach. We then compare our findings with observational data. Finally, in Section 4, we summarize our conclusions and outline directions for future research.

2. Data and methodology

2.1. WDM DREAMS simulations

The DREAMS project¹ consists of an extensive suite of hydrodynamic simulations, complemented by N-body simulations that serve as their purely gravitational counterparts. In this section, we provide a description of the simulations used in this study, specifically those based on the WDM model, but we refer the reader to Rose et al. (2025) for a more detailed description of the DREAMS project and the simulation setup.

Depending on the simulated environment, the simulations are organized into two different suites, each conducted in regions of different spatial scales and resolutions: cosmological boxes and MW zoom-in. The uniform-box simulations span a volume of (25 h⁻¹ Mpc)³, adopt periodic boundary conditions, and include a 2 × 256³ grid of DM and gas particles. The DM mass resolution depends on the matter density parameter Ω_m and is given by 7.81 × (Ω_m/0.302) × 10⁷ h⁻¹ M_⊙, while the baryonic mass resolution is fixed at 1.27 × 10⁷ h⁻¹ M_⊙. The spatial resolution, set by the gravitational softening length, reaches 1.0 h⁻¹ kpc at z = 0. These simulations offer good statistical sampling thanks to their volume. However, the relatively low resolution of these runs, comparable to TNG300-1, limits the efficiency of baryonic feedback processes, especially those linked to star fformation and SNe. As a result, the stellar-DM relation tends to deviate from the results of higher-resolution simulations such as TNG100-1 and TNG50-1 (Busillo et al. 2025), where feedback is better resolved and more effective. To mitigate these resolution-dependent biases, we calibrate the stellar-DM mass relation of the uniform-box simulations by aligning them with the high-resolution TNG100-1 simulation, which serves as our reference. This empirical calibration improves the agreement with observational data and enhances the predictive reliability of the simulated galaxy populations. The calibration introduces a stellar-mass shift of approximately c ∼ 0.27 dex in log₁₀M_*/M_⊙. A detailed description of the calibration procedure, including a discussion of its limitations, is provided in Appendix A.

Although the relatively low mass resolution of the uniform-box runs limits their ability to fully resolve small-scale baryonic processes, recent analyses have demonstrated that simulations at comparable resolution can still retain sensitivity to the WDM particle mass. In particular, Rose et al. (2024) showed that, for fixed cosmology, convolutional neural networks trained on DM density fields can recover WDM masses up to ∼7 keV from N-body simulations with 256³ particles in a (25 h⁻¹ Mpc)³ box–identical to the DREAMS setup. When cosmological parameters are varied but the simulation resolution is kept fixed, their model remains accurate for WDM masses in the range 3.5–5 keV, though with slightly larger uncertainties. Furthermore, Costanza et al. (2025) used the DREAMS hydrodynamical uniform boxes to show that galaxy population statistics, such as stellar and gas mass distributions, preserve strong discriminatory power between WDM and CDM models up to ∼6 keV. These results motivate the inclusion of the uniform-box suite in our analysis as a statistically representative yet still physically sensitive component of the DREAMS framework.

Complementing this setup, the MW zoom-in simulations simulate MW-mass analogs, with no nearby massive neighbor, and target smaller regions with a characteristic scale of approximately 200 h⁻¹ kpc. These simulations achieve higher resolution, only a factor of two lower in mass resolution than TNG50-1, with DM particle masses of 1.2 × 10⁶ h⁻¹ M_⊙, average baryonic masses of 1.9 × 10⁵ h⁻¹ M_⊙, and a spatial resolution of 0.31 h⁻¹ kpc at z = 0. For this reason, no empirical calibration is applied to the MW zoom-in simulations, which are expected to capture feedback processes more accurately.

The construction of the MW zoom-in simulations involves a multistep procedure (see Appendix A in Rose et al. 2025), which includes intermediate DM-only runs at progressively increasing resolution. As a result of this process, the final high-resolution region can be surrounded by particles from lower-resolution region, and some sub-halos – particularly in the outskirts – can be partially or predominantly composed of low-resolution DM particles. However, by design, the high-resolution region extends to ∼5 times the virial radius of the MW at z = 0, ensuring that the main halo and most of its satellites lie well within the high-resolution volume, and are therefore minimally affected by low-resolution contamination. In cases where contamination could be significant, key baryonic processes such as gas cooling and star formation cannot be reliably modeled. To ensure physical consistency, we exclude from our analysis all galaxies in which more than 5% of the total DM mass originates from low-resolution particles.

It is also important to highlight that, although the central galaxy in each MW zoom-in simulation resides within a halo carefully selected to have a virial mass comparable to that of the MW at z = 0, the resulting galaxy properties are not guaranteed to resemble those of the MW itself. The diversity in feedback strength, and therefore in star formation efficiency, can lead to significant variations in stellar mass and morphology across the sample (see Appendix E). Thus, the MW-like host galaxies should not be interpreted as strict analogs of the MW, but rather as realizations within the broader class of galaxies hosted by MW-mass halos.

In the simulations of interest, cosmological parameters can either be variable (for the uniform-box runs) or fixed (for MW zoom-in). In the first case, Ω_m varies within the range [0.1, 0.5] and σ₈ within [0.6, 1.0], both sampled from a uniform distribution. For the MW zoom-in, the reference values from the Planck Collaboration (Planck Collaboration XI 2016) are adopted: Ω_m = 0.302 and σ₈ = 0.839. Other cosmological parameters, fixed in both suites according to the Planck Collaboration results, are Ω_Λ = 0.698 and H₀ = 100 h km s⁻¹ with h = 0.691. N-body simulations assume Ω_b = 0, while in hydrodynamic simulations, Ω_b = 0.046. This study focuses exclusively on hydrodynamic simulations. In this context, we note that because Ω_m is varied while Ω_b is fixed in the hydrodynamical runs, the ratio f_b = Ω_b/Ω_m, which sets the cosmic baryon fraction, also varies considerably. This ratio is relatively well constrained by cosmological observations (e.g., from the cosmic microwave background), typically around 0.15. In the simulations, however, it can depart significantly from that value, depending on the choice of Ω_m. These variations can influence the overall efficiency of galaxy formation and lead to increased scatter in quantities such as the stellar-to-halo mass ratio, beyond what would be expected if Ω_b and Ω_m were varied consistently within observational bounds.

2.1.1. Astrophysical parameters and WDM mass

The WDM simulations employed in this study are based on the TNG galaxy formation model (Springel et al. 2018; Pillepich et al. 2018; Nelson et al. 2018; Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2019; Pillepich et al. 2019). They start from redshift z = 127 and evolve down to z = 0 using the moving-mesh code AREPO (Springel 2010, Springel et al. 2019, Weinberger et al. 2020). These simulations adopt the fiducial parameters of the TNG model, except for those governing SN and AGN feedback, and the matter power spectrum.

The TNG model incorporates a feedback mechanism to describe the interactions between star formation, supermassive black holes, and the circumgalactic medium. To understand how these processes affect DM distribution and galaxy growth, three key parameters related to stellar winds (A_SN1 and A_SN2) and AGN activity (BH_FF) are varied in the simulations².

The first parameter is the wind energy, A_SN1, which enters the expression for the wind mass-loading factor, η_w, via e_w, the specific energy available to driving winds:

$$\begin{array}{c}\hfill {\eta }_w=\frac{2}{{{\nu }_w}^2}{e}_w(1-{\tau }_w).\end{array}$$(1)

Here, ν_w denotes the wind velocity (see Eq. 3) and τ_w is the fraction of energy released thermally, fixed at its fiducial value of 0.1. The parameter e_w is expressed as follows:

$$\begin{array}{cc}\hfill e_w& =A_{\mathrm{SN}1}\times {\overline{e}}_w[f_{w,Z}+\frac{1-f_{w,Z}}{1+{(Z/Z_{w,\mathrm{ref}})}^{{\gamma }_{w,Z}}}]\hfill \\ \hfill & \times N_{\mathrm{SNII}}{E}_{\mathrm{SNII},51}{10}^{51}{\mathrm{erg}\mathrm{M}}_{\odot }^{-1},\hfill \end{array}$$(2)

where Z is the metallicity of gas cells, ${\overline{e}}_w$ is the wind energy factor, and f_w,  Z accounts for the metallicity-dependent reduction. The reference metallicity is Z_w,ref, while γ_w,  Z sets the power of the metallicity-dependent reduction. Finally, N_SNII represents the number of Type II SNe per unit stellar mass formed, and E_{SNII,  51} is the energy available per core-collapse SN in units of 10⁵¹ erg.

The second parameter is the galactic wind velocity, A_SN2, which depends on the properties of the DM halo and redshift. It is governed by ν_w, determining the efficiency of gas ejection:

$$\begin{array}{c}\hfill {\nu }_w=max[A_{\mathrm{SN}2}{\kappa }_{w,\mathrm{ref}}\times {\sigma }_{\mathrm{DM}}{\left(\frac{H_0}{H(z)}\right)}^{1/3},{\nu }_{w,\min }].\end{array}$$(3)

In this expression, the parameter A_SN2 = κ_w/κ_w,  ref, where κ_w is a dimensionless normalization factor with fiducial TNG value κ_w,  ref = 7.4. The remaining quantities include σ_DM, the DM velocity dispersion, H(z), the Hubble function, and ν_w,min = 350 km s⁻¹, the minimum allowed wind speed. Increasing ν_w suppresses star formation at later cosmic times.

The third key parameter is the AGN feedback efficiency, BH_FF. Supermassive black holes influence gas heating and ejection, regulating star formation in massive galaxies. The model includes two feedback modes: a low-accretion state, which heats gas more uniformly, and a high-accretion state, which drives energetic outflows. The WDM DREAMS simulations vary only the high-accretion mode. Here, BH_FF = ϵ_f,  high/ϵ_{f,  high,  ref}, where ϵ_f,  high determines the fraction of accreted energy transferred to the surrounding gas, and ϵ_{f,  high,  ref} = 0.1 is its reference value:

$$\begin{array}{c}\hfill \mathrm{\Delta }\dot{E}={\mathrm{BH}}_{\mathrm{FF}}{ϵ}_{f,\mathrm{high},\mathrm{ref}}\times {ϵ}_r{\dot{M}}_{\mathrm{BH}}{c}^2,\end{array}$$(4)

with ϵ_r is the radiative efficiency and Ṁ_BH is the black hole accretion rate. Higher values of BH_FF enhance AGN thermal feedback, suppressing star formation in massive galaxies. This parameter, also known as BH feedback factor, is implemented – together with other parameters regulating AGN feedback – in the latest suite of CAMELS simulations (Ni et al. 2023; Tortora et al. 2025).

The TNG model assumes fiducial values for A_SN1, A_SN2 and BH_FF, equal to 1, but DREAMS simulations explore a broader range to assess their impact on galaxy evolution. Specifically, the parameter A_SN1 and BH_FF varies within the interval [0.25, 4], while A_SN2 spans the range [0.5, 2]. The astrophysical parameters are selected from a uniform distribution in logarithmic space.

In the WDM DREAMS simulations, the WDM is treated as a collisionless and pressureless fluid, similarly to CDM, and is modeled by a single parameter, P_WDM, defined as the inverse of the WDM particle mass (P_WDM = 1/M_WDM). In the case of uniform-box simulations, P_WDM is uniformly chosen in the interval [0.062, 0.556] keV⁻¹, corresponding to WDM particle masses between ∼1.8 and 16 keV. In the MW zoom-in simulations, P_WDM spans a broader range of [0.033, 0.556] keV⁻¹, meaning that the WDM mass varies between ∼1.8 and 30 keV. P_WDM, as well as the cosmological and astrophysical parameters, is sampled using a Sobol sequence (Sobol’ 1967), which ensures homogeneous coverage of the multidimensional parameter space.

The parameter P_WDM fully determines the distribution of the WDM thermal velocities. Operationally, in the simulations, the WDM models differ from the CDM models in the initial conditions. Specifically, following Bode et al. (2001), the particle mass parameter enters the linear power spectrum through a suppression factor β(k). As a consequence, the WDM power spectrum can be written as P_WDM(k) = β(k)P_CDM(k), which is damped at large comoving wave numbers, i.e., on small scales,

$$\begin{array}{cc}\hfill \beta (k)& =({(1+{(\alpha k)}^{2.4})}^{-5.0/1.2}{)}^2,\hfill \\ \hfill \alpha & =0.048h^{-1}\mathrm{Mpc}{\left(\frac{M_{\mathrm{WDM}}}{\mathrm{keV}}\right)}^{-1.15}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}}-{\mathrm{\Omega }}_{\mathrm{b}}}{0.4}\right)}^{0.15}{\left(\frac{h}{0.65}\right)}^{1.3}.\hfill \end{array}$$(5)

It is worth noting that in the limit where M_WDM → ∞, β(k)→1, and the power spectrum of the WDM coincides with that of the CDM. In this way, the DREAMS simulations cover a wide range of WDM models using a single parameter.

2.1.2. Simulations employed, astrophysical quantities, and selection criteria

The WDM DREAMS simulation suite provides 1024 distinct configurations of astrophysical, cosmological, and WDM mass parameters for the uniform-box runs. Additionally, an equivalent number of simulations is available for the MW zoom-in runs, where the cosmology is fixed.

Sub-halo identification is performed using the friends-of-friends and SUBFIND algorithms (Springel et al. 2001). In this framework, each sub-halo is defined as a gravitationally self-bound structure, and its total mass is computed as the sum of all bound particles, including DM, gas, stars, and black holes, without relying on a spherical overdensity threshold. In this work, we compare simulated galaxy scaling relations with observational counterparts by analyzing several key astrophysical quantities derived from SUBFIND at redshift z = 0:

• Stellar half-mass radius, R_*, 1/2: The radius enclosing half of the total stellar mass of the galaxy.

• Total stellar mass, M_*: Obtained as the sum of the stellar particle masses within a given sub-halo.

• Total mass, M_tot: The sum of all mass components (DM, stars, gas, and black holes) in a sub-halo.

• Mass components within the stellar half-mass radius: Stellar mass (M_*, 1/2), DM mass (M_DM, 1/2), and total mass (M_tot, 1/2).

• The DM fraction within the stellar half-mass radius, f_DM, 1/2 ≡ f_DM(< R_*, 1/2)≡M_DM, 1/2/M_tot, 1/2.

• Total-to-stellar mass ratio within the half-mass radius, M_tot, 1/2/M_*, 1/2.

• Number of star particles within the stellar half-mass radius, N_*, 1/2.

• Star formation rate (SFR).

• Number of galaxies, N_Galaxy: The total number of sub-halos.

To ensure a robust comparison between simulations and observational datasets, we imposed resolution-dependent selection criteria:

• The value of R_*, 1/2 must exceed the spatial resolution, which is 1 kpc for uniform-box simulations and 0.305 kpc for MW zoom-in runs.

• The DM fraction within the half-mass radius, f_DM, 1/2, must be positive for both simulation suites.

• The number of stellar particles within the half-mass radius, N_*, 1/2, must exceed 50 for both simulation suites, corresponding to stellar mass thresholds of approximately log₁₀M_*/M_⊙ ∼ 9.3 for the uniform-box simulations and ∼7.4 for the MW zoom-in simulations. This selection is applied when comparing the simulations to observational catalogs.

These criteria ensure that only well-resolved galaxies are included in the analysis and that a consistent comparison with observations can be made across all simulation sets. WDM simulations, however, may still be affected by numerical artifacts known as “beads-on-a-string,” which arise from discreteness noise along filaments and can produce spurious small-scale clumps (e.g. Wang & White 2007). Such effects become relevant near the cutoff scale of the WDM power spectrum, where the half-mode mass approaches the particle-mass resolution. As discussed by Rose et al. (2024), numerical fragmentation becomes significant only for very warm models (M_WDM ≲ 1.8 keV). Since this value marks the lower bound of our explored range, we expect that the vast majority of our simulations are unaffected by such artifacts and that any residual spurious clumping at the lowest WDM mass is minimal and does not impact our results. The scaling relations analyzed in this work are consistently measured for both the uniform-box and MW zoom-in simulations, focusing on well-resolved galaxies above our stellar particle-number threshold (i.e., N_*, 1/2 > 50). Spurious haloes generated by numerical fragmentation are typically dark and therefore excluded by these cuts. Moreover, the abundance of satellites – where such artifacts could, in principle, be more apparent – is examined only in the higher-resolution MW zoom-in suite, which further mitigates possible contamination. Together with the findings of Rose et al. (2024), these considerations indicate that numerical fragmentation does not materially affect the trends discussed in this paper.

2.2. Observational catalogs

In the following, we introduce the two observational catalogs, the SPARC and the dwarf galaxy catalog, which serve as benchmarks for comparing the trends obtained from simulations. We specify the astrophysical quantities extracted from these datasets that contribute to the analyzed scaling relations and detail the methods used to derive them.

2.2.1. SPARC catalog overview

The SPARC catalog (Lelli et al. 2016) provides a dataset of 175 disk galaxies, serving as a valuable reference for studying the interplay between baryonic matter and DM in galaxy dynamics. While it does not constitute a volume-limited sample, it includes a wide variety of late-type galaxy morphologies, rotational properties, and sizes, offering a representative view of disk galaxies in the local Universe. The SPARC galaxies span a stellar mass range of 6.9 ≲ log₁₀M_*/M_⊙ ≲ 11.4, covering a wide range of galaxy masses from low-mass disks to massive spirals.

The near-infrared photometry at 3.6 μm of SPARC was analyzed to estimate stellar masses, adopting a standard stellar mass-to-light ratio (Υ_* = 0.6 Υ_⊙). Busillo et al. (Paper I) tested the impact of this assumption by repeating the fitting-like analysis – which is introduced in Sect. 2.3 – with stellar mass-to-light ratios in the range Υ_* = 0.5 Υ_⊙ to Υ_* = 0.7 Υ_⊙. While this range leads to differences in stellar mass estimates of about 0.06–0.08 dex, the resulting variations in the astrophysical and cosmological parameters were found to be negligible, remaining well within the statistical uncertainties. The catalog’s high-resolution HI 21 cm rotation curves allow for the reconstruction of the radial distribution of mass within each galaxy. The total enclosed 3D mass at a given radius r is calculated using the relation M(r) = v²r/G, allowing a direct comparison between the observed baryonic components and the inferred DM contributions.

Selection criteria further refine our working sample: galaxies with inclinations below 30° are excluded to minimize projection effects, and only systems with well-measured rotation curves extending beyond the effective radius are retained, leading to a final selection of 152 galaxies.

To compare our results with the SPARC dataset, we select simulated galaxies based on their specific star formation rate (sSFR). We classify galaxies with sSFR > 10^−10.5 yr⁻¹ as star-forming galaxies (SFGs), while we consider those with sSFR < 10^−10.5 yr⁻¹ passive galaxies (PGs). This selection is applied consistently throughout Papers I, II, and III, providing a practical way to isolate actively star-forming systems in the simulations. Although the SPARC sample targets late-type galaxies, these are generally star-forming in the local Universe. Selecting simulated galaxies by sSFR thus offers a practical proxy to identify SPARC analogs. This allows a meaningful comparison between simulations and observations.

In this study, we focus on several key parameters, including the total stellar mass, M_*; the stellar half-mass radius, R_*, 1/2, defined multiplying the respective effective radii (R_e), which encloses half of the total 3.6 μm luminosity, by a constant factor of ∼1.35 (Wolf et al. 2010); the stellar, gas and total mass within this radius, M_*, 1/2, M_g, 1/2 and M_tot, 1/2, respectively. From these, we derive the enclosed DM mass, M_DM, 1/2, and compute the DM fraction, f_DM, 1/2. To complement these calculations, we incorporate virial mass estimates, M_vir, from Posti et al. (2019), obtained by fitting rotation curves using a Navarro-Frenk-White (NFW) DM halo profile (Navarro et al. 1996).

2.2.2. Dwarfs in the local volume

We also exploited the LVDB (Pace 2025) sample that is publicly available as a GitHub repository³. LVDB is an extensive compilation of observed properties for astrophysical objects within the Local Volume, encompassing dwarf galaxies, globular clusters, and other stellar systems. The database primarily focuses on resolved stellar systems within approximately 10 Mpc, though its completeness is highest for objects within 3 Mpc, and it provides a valuable resource for studying galaxy formation, structure, kinematics, and chemical evolution. The catalog covers a stellar mass range of 2.7 ≲ log₁₀M_*/M_⊙ ≲ 9.7, offering a broad view of low-mass galaxy populations in the nearby Universe. The LVDB incorporates data from multiple sources, including the McConnachie (2012) catalog, the Catalog and Atlas of Local Volume Galaxies (Karachentsev et al. 2013), and the Extragalactic Distance Database (Tully et al. 2009, Anand et al. 2021) and it is regularly updated to include new discoveries.

The LVDB contains key astrophysical parameters for each object, including structural, kinematic, and dynamical properties. For this study, we extract specific quantities necessary for establishing scaling relations, which are subsequently compared to numerical simulations. The primary parameters used in our analysis are:

• Stellar mass: The parameter mass_stellar stored in the catalog is defined as the logarithm of the stellar mass, estimated assuming a mass-to-light ratio of two, and derived from the absolute V-band magnitude (M_V) [log₁₀ M_⊙].

• Half-light radius: rhalf_sph_physical is the azimuthally averaged 2D half-light radius, computed as the geometric mean incorporating ellipticity corrections, expressed as $\mathtt{rhalf}\times \mathtt{distance}\times \sqrt{(1-\mathtt{ellipticity})}$ in parsecs. Here, rhalf represents the major axis of the half-light radius (or Plummer radius) in arcminutes, distance is the heliocentric distance, and ellipticity = 1 – minor/major axis.

• Dynamical mass within the half-light radius: mass_dynamical_wolf is the dynamical mass enclosed within the 3D half-light radius [M_⊙].

The 3D half-mass radius is determined following the method outlined in Wolf et al. (2010), as done for the SPARC catalog:

$$\begin{array}{c}\hfill R_{\ast ,1/2}[\mathrm{kpc}]=4/3\times \frac{\mathtt{rhalf}\_\mathtt{sph}\_\mathtt{physical}}{1000}.\end{array}$$(6)

The DM mass enclosed within the half-mass radius is computed as:

$$\begin{array}{c}\hfill M_{\mathrm{DM},1/2}[{\mathrm{M}}_{\odot }]=\mathtt{mass}\_\mathtt{dynamical}\_\mathtt{wolf}-0.5\times {10}^{\mathtt{mass}\_\mathtt{stellar}}.\end{array}$$(7)

Due to the lack of robust gas mass estimates for the majority of the dwarf galaxies in the LVDB, we do not include a gas mass contribution in our analysis. The fraction of DM within the half-light radius was subsequently derived as

$$\begin{array}{c}\hfill f_{\mathrm{DM},1/2}=\frac{M_{\mathrm{DM},1/2}}{\mathtt{mass}\_\mathtt{dynamical}\_\mathtt{wolf}}.\end{array}$$(8)

Total virial masses of the dwarf galaxies are not included in the catalog and therefore will not be included in the next analysis.

To ensure a meaningful comparison with simulations, we applied specific selection criteria to the LVDB. The subset of 14 dwarf galaxies used in this study satisfies the following conditions:

• The system must be classified as a confirmed physical galaxy (confirmed_real = 1) and a dwarf galaxy (confirmed_dwarf = 1).

• The stellar mass must satisfy mass_stellar > 7.4 [log₁₀ M_⊙], ensuring consistency with the resolution limit (N_*, 1/2 > 50) adopted in the WDM zoom-in simulations.

• The stellar half-mass radius must satisfy R_*, 1/2[kpc]> 0.305 kpc, ensuring that the selected objects have spatial resolutions above the simulation limits.

• The dynamical mass estimate must be positive and exceed the enclosed stellar mass (mass_dynamical_wolf > 0.5 × 10mass_stellar).

• The catalog must simultaneously contain all the quantities described above for the selected galaxies to ensure a complete dataset for comparison with simulations.

2.3. Analysis approach

To compare the simulated trends with observational data, we adopt a fitting-like procedure, implemented in Mathematica, and designed to account for statistical uncertainties and parameter degeneracies. This method was first introduced in Paper I and Paper II, and is adapted here to the current simulations and data. In this context, we define the distance estimator D² as in Papers I and II, but with a minor change in notation. This modification is intentional: we aim to emphasize that our fitting-like method and the D² estimator do not correspond to the standard definitions of fitting or the classical χ² statistic. Instead, they represent an alternative approach developed in the context of the CASCO paper series, specifically designed to assess how well the simulations reproduce the observed trends in selected galaxy catalogs. Despite the differences, our method shares several features with the classical χ². Below, we outline the main points of our analysis.

In this work, we focus on four key scaling relations used in the fit: the size–mass relation (R_*, 1/2 versus M_*); the DM mass and fraction within the stellar half-mass radius (M_DM, 1/2 versus M_* and f_DM, 1/2 versus M_*, respectively); and the relation between total sub-halo mass and stellar mass (M_tot versus M_*). In addition, we considered one relation not included in the fit but used elsewhere in our analysis. Specifically, we included the total-to-stellar mass ratio within the half-mass radius (M_tot, 1/2/M_*, 1/2) in our plots, as it offers improved visual clarity compared to f_DM, 1/2. Since both quantities carry essentially the same information, we use f_DM, 1/2 in the fit to remain consistent with our previous works, where this methodology was first implemented.

For the four scaling relations included in the fit, we first derive the observational trends by binning stellar mass and computing the 16th, 50th, and 84th percentiles in each bin. This provides the median trend and the scatter of the observational distribution. As discussed in Paper II, different binning choices can slightly influence this step. We then interpolate these percentiles linearly to define the observational function f_obs(x). For each simulation, we extract the corresponding scaling relations, assigning each data point coordinates (x_sim,  i,  y_sim,  i), where i indexes the simulation data and n is the total number of points in that simulation. We computed the distance estimator as

$$\begin{array}{c}\hfill D^2={\displaystyle \sum _{i=1}^n}\frac{{[y_{\mathrm{sim},\mathrm{i}}-{\mathcal{N}}_{\mathrm{i}}(f_{\mathrm{obs}}(x_{\mathrm{sim},\mathrm{i}},{\sigma }_{\mathrm{obs},\mathrm{i}}))]}^2}{{\sigma }_{\mathrm{obs},\mathrm{i}}^2},\end{array}$$(9)

where the numerator is the squared difference between the ordinate of the simulation point i, y_sim,  i, and a randomly extracted point from a Gaussian distribution, 𝒩, centered on the median of the observed scaling relation with a standard deviation equal to σ_obs, i. The quantity σ_obs,  i is defined as the average of σ₊ and σ₋, which are the absolute differences between the interpolated 16th and 84th percentile trends and the interpolated median trend at x_sim,  i. At this stage, a few clarifications are necessary. First, we emphasize that we used the Gaussian distribution to model observational uncertainties, as it offers a simple yet effective way to account for the observed scatter in the data, a choice that, as we demonstrate, proves useful later in the analysis. It should be noted, however, that this approach does not explicitly account for potential selection effects in the observational data. Specifically, it assumes that the data points are independent and identically distributed samples from the underlying galaxy population. A more accurate treatment of these biases would require advanced modeling of the selection effects, which we leave for future work. Second, the quantity σ_obs,  i accounts for the total observed scatter, encompassing both measurement errors and intrinsic astrophysical dispersion. Rather than comparing simulations to a single median trend, our analysis assesses how well they reproduce the full observed distribution.

Since each simulation contains a different number of data points, it is essential to normalize the D² estimator by this number, ensuring a fair comparison across simulations. We therefore define the reduced distance estimator as ${\stackrel{\sim }{D}}^2=D^2/(n-1)$. Furthermore, to avoid extrapolations, we exclude from the fit any simulation points where x_sim,  i falls outside the range of the observational median trend. The best-fitting simulation is the one with the minimum reduced distance estimator.

A simple fitting procedure, such as that described above, does not provide any uncertainty estimates for the best-fit parameters. Moreover, parameter degeneracy complicates the analysis: different parameter combinations can yield similarly good fits or produce very similar trends, making it difficult to identify a unique optimal configuration. This issue, which we discuss further in our work, motivates the use of a more robust approach: bootstrap resampling. By applying the bootstrap method to the observational dataset, the simulations, and the simulation output, we obtain a clearer statistical interpretation of the fit. Notably, since we have introduced a degree of randomness in the observational trends – by sampling from Gaussian distributions around the observed percentiles – this approach also allows us to estimate an uncertainty on the D² values themselves at the end of the bootstrap procedure.

We followed the steps described below. Using Mathematica ResourceFunction["BootstrapStatistics"]⁴, we generated bootstrap samples, subsets of the original data with replacement. We performed resampling with N ≈ 100 iterations, creating a larger statistical sample. Each resampled dataset undergoes the same D² analysis, producing a distribution of N best-fit simulations. This distribution allows us to estimate the best-fit parameters using the median and quantify their scatter as the absolute difference between the median and the 16th/84th percentiles. The statistical reliability of the bootstrap resampling procedure, and its potential impact on parameter recovery and uncertainty estimates, have been explicitly tested through a series of dedicated analyses presented in Appendix F.2. We note that owing to the discrete nature of the sample, the 16th and 84th percentiles may occasionally coincide with the median, especially for sharply peaked distributions. This analysis provides both estimates of parameter uncertainties and insights into possible degeneracies. By examining the distribution of best-fit simulations obtained from the bootstrap resamples, each associated with a specific set of model parameters, we can assess the stability of the inferred results. A broad distribution of best-fit simulations indicates that multiple parameter configurations yield comparably good fits to the data, revealing potential degeneracies. Conversely, a narrow distribution may suggest that the best-fit parameters are robustly constrained. This approach thus allows for a more reliable evaluation of the model’s sensitivity to the data and the significance of the resulting parameter estimates.

3. Results

In this section, we present the main results of our study. In Section 3.1, we show the impact of cosmological and astrophysical parameters, as well as the WDM mass, on the scaling relations. We perform this analysis using both uniform-box and MW zoom-in simulations. In addition, for the MW zoom-in simulations, we also evaluate the effect of these parameters on galaxy abundance. In Section 3.2, we assess the ability of our fitting method, supported by bootstrap resampling, to accurately estimate the parameters for both simulation suites. Finally, in Section 3.3, we compare the simulation trends with observational data from the SPARC and LVDB catalogs, determining the set of parameters that best reproduces the observations.

Fig. 1. Median scaling relations as a function of total stellar mass illustrating the dependence on cosmological parameters. The relations for stellar half-mass radius, total-to-stellar mass ratio within the stellar half-mass radius, DM mass within the stellar half-mass radius, and total galaxy mass are shown for both uncalibrated (solid line) and calibrated (dashed line) cases. Trends were calculated using galaxies from 1024 uniform-box simulations divided into the following stellar mass intervals: log10M*/M⊙ < 9.3, 9.3 ≤ log10M*/M⊙ < 9.6, 9.6 ≤ log10M*/M⊙ < 10.4, 10.4 ≤ log10M*/M⊙ < 11.2, and log10M*/M⊙ ≥ 11.2. Red and blue indicate trends for high and low parameter values, respectively, while the dashed gray line marks the threshold N*, 1/2 = 50. The thresholds and median values for the parameters are provided in Table 1. Error bars represent the scatter (16th–84th percentile range), while shaded areas indicate uncertainty estimates on the median values.

3.1. Galactic scaling relations in WDM cosmology

To investigate how cosmological and astrophysical parameters, including the WDM particle mass, influence galactic scaling relations, we use the full set of simulations. To clarify our methodology, we take the parameter Ω_m as an example and describe our approach step by step; the same procedure is then applied to each of the other parameters.

First, we aim to show how scaling relations vary for extreme values of the parameter. We select two threshold values for Ω_m, labeled “Low Ω_m” = 0.25 and “High Ω_m” = 0.35 (all threshold values are listed in Table 1). From the full simulation set, we construct two subsamples containing simulations with Ω_m values below the low threshold and above the high threshold, respectively. We then ensure that the distributions of the remaining parameters in each subsample remain similar to their original distributions in the full sample and close to their fiducial values. Next, we divide the stellar mass range of the simulations into intervals and, for each interval, compute the median of the quantities of interest, along with the 16th and 84th percentiles. These values are then used to reconstruct the median scaling relations by combining the points from each mass interval. Median statistics are robust against the influence of outliers, and evaluating all simulations together within specific mass ranges allows us to assess the impact of the parameters on galaxies of different masses while maximizing the statistical sample size.

On the other hand, the scatter does not represent an uncertainty on the median trend itself, but rather captures the statistical variation in the trends due to changes in the remaining parameters. As sources of uncertainty, we first consider the uncertainty on the medians, which is generally small and increases only moderately at higher stellar masses due to the reduced number of galaxies in those bins, reaching at most 0.08 dex. Second, we consider the error in estimating the stellar half-mass radius, which propagates into the measurements of the quantities computed within that radius. Details of this calculation are provided in Appendix C. The uncertainty tied to the R_*, 1/2 calculation depends on the number of stellar particles in the simulated galaxy, and can reach approximately 0.04–0.05 dex when N_*, 1/2 = 20. For each scaling relation, we add in quadrature the sources of uncertainties on the y-axis: namely, the error of the median, and that on the plotted quantity (discussed in Appendix C), if the latter is – or is computed within – R_*, 1/2. This procedure is applied consistently to both the low and high trends. We show the total uncertainties in the trends using shaded regions, while representing the scatter of the medians as error bars. In this section, we present the trends for galaxies with a minimum contribution of N_*, 1/2 = 20, along with the corresponding uncertainty estimates. However, when comparing with observational catalogs, we adopted a more conservative cut of N_*, 1/2 > 50. Overall, the use of the full simulation set allows for precise estimates of median values, thanks to the large number of galaxies in each stellar mass bin.

3.1.1. Uniform-box suite

We begin by analyzing the uniform-box simulations, to which we apply the selection criteria described in Sect. 2.1.2. The parameter thresholds are listed in Table 1, which shows, for each parameter, the threshold ranges along with the corresponding median values computed within each subsample. These median values remain consistent across the different stellar mass intervals.

Figure 1 illustrates the impact of cosmological parameters and the WDM particle mass on the galactic scaling relations discussed in Sect. 2.3, and also shows the corresponding trends obtained from the calibrated models (dashed lines) for comparison. Figure 1 shows that Ω_m has a significant impact on the total-to-stellar mass ratio and the DM mass within the stellar half-mass radius, shifting these relations toward higher DM fractions and greater DM content. Moreover, the total galaxy mass is higher for larger values of Ω_m, with this effect becoming more significant as the stellar mass increases and, in the lower stellar mass bins, higher values of Ω_m correspond to larger R_*, 1/2. The influence of σ₈ is generally weaker than that of Ω_m. The DM content and the total-to-stellar half-mass ratio within R_*, 1/2 show a mild increase at low stellar masses for higher values of σ₈, but this behavior reverses around log₁₀M_*/M_⊙ ∼ 10.6. Conversely, the size–mass relation appears to flatten for elevated σ₈ values, starting from galaxies with log₁₀M_*/M_⊙ ∼ 10. In addition, the total mass decreases with increasing parameter values, with little variation across stellar mass bins. Moreover, we find that the σ₈ also induces a slight horizontal shift in the highest mass bin, where higher values of σ₈ are associated with higher total stellar masses. These results are consistent with the findings of Papers I and II, where the effects of cosmological parameters in a ΛCDM framework were analyzed separately for SFGs and PGs, respectively. In our analysis, we do not differentiate between these morphological types. However, as discussed in Paper II, the influence of cosmological parameters is independent of galaxy type. Moreover, we find that the WDM particle mass does not significantly alter the trends in the scaling relations. However, we do observe a horizontal shift, with higher WDM particle masses associated with slightly lower stellar masses in the most massive bin. Finally, the comparison with the calibrated trends indicates that calibration effects become more pronounced at higher stellar masses, resulting in consistently lower values for all measured quantities compared to the uncalibrated case. This behavior is expected, and the same effect is naturally observed in Fig. 2, where we show how the scaling relations vary as a function of the astrophysical parameters.

Fig. 2. Same as Fig. 1, but for the astrophysical parameters ASN1, ASN2, and BHFF.

Fig. 3. Median scaling relations as a function of total stellar mass showing the dependence on astrophysical parameters and WDM mass for both MW zoom-in simulations (solid line) and calibrated uniform-box simulations (dashed line) for comparison. Trends are derived from 1024 MW zoom-in simulations divided into the following stellar mass intervals: log10M*/M⊙ < 7.5, 7.5 ≤ log10M*/M⊙ < 8.5, 8.5 ≤ log10M*/M⊙ < 10, and log10M*/M⊙ ≥ 10. Red and blue indicate trends for high and low parameter values, respectively. The dashed gray line indicates N*, 1/2 = 50.

Fig. 4. Differences between the median trends corresponding to the high and low values of two parameters: ASN1 (green) and MWDM (brown). Solid lines represent MW zoom-in simulations; dashed lines represent uniform-box runs. Shaded regions show propagated uncertainties. The dashed gray line indicates N*, 1/2 = 50.

Fig. 5. Median number of galaxies as a function of median total stellar mass computed in the stellar mass bins log10M*/M⊙ < 6.5, 6.5 ≤ log10M*/M⊙ < 7.5, 7.5 ≤ log10M*/M⊙ < 8.5, and log10M*/M⊙ ≥ 8.5. The four panels, from top to bottom, correspond to variations in ASN1, ASN2, BHFF, and MWDM, respectively. In each case, the red and blue curves represent parameter values above and below the thresholds listed in Table 1, respectively. The plots are based on the 1024 MW zoom-in simulations from DREAMS.

Regarding the effect of astrophysical parameters on scaling relations, the wind energy–parameterized by A_SN1–has the strongest impact, systematically shifting all four scaling relations upward (Fig. 2). In addition, in the highest stellar mass bin, larger values of A_SN1 suppress star formation, leading to a slight reduction in the median stellar mass of galaxies. Notably, the A_SN1 parameter operates almost independently of the total stellar mass of the galaxy within the analyzed mass range, a result consistent with Papers I, II and III, where similar trends were observed for SFGs and PGs, respectively, using CAMELS simulations.

The effects of wind velocity (i.e., A_SN2) on stellar mass accretion are primarily evident at high stellar masses. In particular, higher values of the parameter increase, at fixed stellar mass, the DM content and the total-to-stellar mass ratio within the half-mass radius, as well as the total sub-halo mass. The size–mass relation, on the other hand, exhibits a turning point around log₁₀M_*/M_⊙ ∼ 10.6: above this threshold, larger A_SN2 values are associated with larger sizes, whereas below it, the trend reverses and sizes decrease. Since these effects are most pronounced at high masses, we first compare our results with Paper II, which focuses on PGs, typically more massive than SFGs. There, high A_SN2 values consistently increase all four scaling relations, in agreement with our findings. The turning point we observe in the size–mass relation is also seen in Paper III, which reports a similar mass scale of log₁₀M_*/M_⊙ ∼ 10.5. Below this mass scale, the trend inverts: stronger feedback reduces galaxy sizes. This behavior aligns with Paper I, which examines SFGs, generally less massive systems, and finds that large A_SN2 values lead to systematically smaller galaxy sizes.

Finally, the effects of the black hole feedback factor, BH_FF, are broadly consistent with the trends observed in the CAMELS simulations, as presented in Paper III, which showed that, at fixed stellar mass, higher values of the BH_FF parameter reduce the total-to-stellar mass ratio already at stellar masses of log₁₀M_*/M_⊙ ∼ 10. Here, we also examine other scaling relations and find that high BH_FF values are associated with smaller galaxy sizes, lower DM fractions within R_*, 1/2, and lower total galaxy masses. As observed for the A_SN1 parameter, the BH_FF parameter also induces a slight reduction in the total stellar mass of galaxies within the highest stellar mass bin.

Among the scaling relations analyzed in this section, the ratio M_tot, 1/2/M_*, 1/2 is particularly insightful, as it exhibits a minimum that defines the so-called golden mass (Dekel & Birnboim 2006; Moster et al. 2010; Tortora et al. 2010, 2019, 2025). This characteristic mass scale appears as a bend or extremum in several scaling relations and is especially associated with the turnover in the stellar-to-halo mass relation (i.e., the M_tot–M_* relation), or equivalently with the peak in the correlation between star formation efficiency, defined as ϵ_SF = M_*/(M_hf_b), where M_h is the total halo mass, and mass. Physically, the golden mass is thought to arise from the interplay of various processes – including SN and AGN feedback, virial shocks, and cold gas accretion via streams (see Paper III and references therein). A focused analysis of how this feature varies across different cosmological and feedback parameters, and between CDM and WDM models, is presented in Appendix D. This analysis, based on the DREAMS and CAMELS uniform-box simulations at comparable resolution, confirms the presence of a well-defined golden mass in both cosmologies. Notably, WDM simulations tend to exhibit more regular and pronounced minima in the M_tot, 1/2/M_*, 1/2 ratio compared to CDM at the same resolution.

3.1.2. MW zoom-in suite

Figure 3 shows the dependence of the median scaling relations for astrophysical parameters and M_WDM, derived from MW zoom-in simulations, alongside those from the calibrated uniform-boxes, for comparison. The consistency between the two sets of simulations is confirmed as they share the same parameter thresholds (Table 1), and the median values are very similar, except for the WDM mass when M_WDM exceeds 10 keV. In this case, the WDM mass is approximately 12.3 keV in the uniform-boxes compared to 15.4 keV in the MW zoom-in simulations; however, this discrepancy does not impact our analysis, as the scaling relations are largely insensitive to variations in M_WDM at high values of this parameter. From the comparison between the two simulation suites, we observe discrepancies in some of the proposed scaling relations, most notably in the size–mass relation. In particular, higher-resolution simulations tend to yield significantly smaller stellar half-mass radii compared to the uniform-box runs. This likely contributes to the lower trends observed in other relations measured within R_*, 1/2. By contrast, the total stellar mass is broadly consistent with the trends seen in the lower-resolution suite.

These discrepancies can be attributed to a combination of factors. A first and obvious contribution arises from resolution differences between the uniform-box and the MW zoom-in suites, which cannot be fully compensated by a simple stellar-mass calibration. Indeed, such a calibration only corrects for global offsets, without accounting for more complex effects related to gravitational softening, energy transfer between baryons and DM, or the detailed structure of the interstellar medium. As shown in Appendices A and B, calibrating the 1P uniform-box run against TNG100-1 or TNG50-1 significantly improves the agreement of the scaling relations with the reference runs, leaving only moderate residual offsets. It is worth noting that TNG50-1 has a resolution comparable to that of the MW zoom-in simulations, while TNG100-1 is coarser. By contrast, when the same 1P simulation is calibrated against the satellite population extracted from the MW zoom-in suite, the resulting scaling relations remain systematically above the MW satellite trends, despite exhibiting a stellar-mass shift comparable to that obtained when calibrating against TNG100-1 or TNG50-1. This behavior indicates that resolution alone cannot explain the observed offsets: differences of this magnitude are not seen when comparing to TNG50-1 or TNG100-1, which suggests that environmental effects specific to the MW zoom-in simulations play a major role. We verified in Appendix B.2 that an additional calibration based on the stellar half-mass radius, R_*, 1/2, further reconciles the two datasets, substantially improving the match with the MW zoom-in trends, however we do not pursue this approach here. The MW zoom-in galaxies, by construction, probe a highly specific environment – that of satellites within MW-mass haloes – and therefore cannot be considered representative of the general galaxy population at similar stellar masses.

Additional factors and motivations should also be considered. In particular, part of the discrepancy between the trends of the two simulation suites may be attributed to cosmic variance – i.e., the natural statistical variations caused by changing the random seed used to generate the initial conditions in cosmological simulations. Although our DREAMS simulations do not include multiple fiducial realizations (where only the seed is varied) to directly measure this effect, we expect its impact to be small. As a reference, Paper I provides an estimate of the effect of cosmic variance on the CAMELS simulations, finding the scatter in relevant galaxy properties to be at the percent level⁵.

Moreover, we observe a flattening of the size–mass relation that amplifies the discrepancies between the two simulation suites. This flattening appears in the uniform-box simulations at total sub-halo stellar masses below log₁₀M_*/M_⊙ ∼ 10.5, while the MW zoom-in simulations also show some flattening at stellar masses below log₁₀M_*/M_⊙ ∼ 9, maintaining a relatively shallow trend over the full mass range. Such flattening can be attributed to a combination of physical and numerical factors. In particular, Ludlow et al. (2019b) demonstrate that numerical artifacts may arise due to resolution limits, notably two-body scattering effects. In cosmological simulations, the spurious transfer of kinetic energy from massive DM particles to lighter stellar particles leads to artificial size inflation over time. This phenomenon becomes significant when the particle mass ratio μ = m_DM/m_gas exceeds unity, where m_DM and m_gas are the masses of individual DM and gas particles, respectively. Specifically, Ludlow et al. (2019b) demonstrate that simulations from EAGLE project (Schaye et al. 2015) with μ ≳ 5.4 show an artificial increase in galaxy sizes, which occurs at stellar masses below log₁₀M_*/M_⊙ ∼ 9.5, which leads to a systematic flattening of the size–mass relation⁶. In DREAMS uniform box simulations, Ω_m varies from 0.1 to 0.5, causing the particle mass ratio μ to range approximately from 2.0 to 10.2. In contrast, the MW zoom-in simulations have a roughly constant μ ∼ 6.3. Although these numerical effects likely contribute to the observed flattening, they do not fully account for the discrepancy between the two simulation suites. Even if the size–mass relation were extrapolated without flattening to lower stellar masses, a significant offset would remain.

Trends in Fig. 3 confirm that A_SN1 plays a dominant role down to log₁₀ M_∗/M_⊙ ∼ 7, driving an increase in R_*, 1/2 and shifting the total-to-stellar mass ratio and DM mass upward. In other words, at fixed stellar mass, higher values of A_SN1 lead to larger galaxies with lower stellar fractions and a higher DM content, highlighting the efficiency of SN feedback in regulating star formation. On the other hand, the effects of A_SN2 and BH_FF appear to be weak. For BH_FF this is expected in MW zoom-in simulations, which primarily probe lower-mass galaxies, as they tend to regulate star formation and black hole growth more significantly in higher-mass systems. Nevertheless, we find that the trends of R_*, 1/2 with respect to A_SN2 and BH_FF are qualitatively consistent with those observed in uniform-box simulations. Notably, the turning point where the size–mass relation inverts for A_SN2 now occurs at log₁₀M_*/M_⊙ ∼ 9.6. Moreover, we observe that an increase in A_SN2 leads to a reduction in the median stellar mass within the most massive bin, consistent with the expected effects of intensified stellar feedback. Conversely, higher values of BH_FF correlate with increased stellar masses even in lower-mass bins. This trend likely reflects the fact that AGN feedback predominantly impacts galaxies exceeding a certain mass threshold, rather than indicating a direct causal enhancement of stellar mass. Notably, the turning point where the size–mass relation inverts for A_SN2 now occurs at log₁₀M_*/M_⊙ ∼ 9.6. Finally, M_WDM shows a weak influence on the total-to-stellar half-mass ratio within the half-mass radius at low masses. This effect suggests that higher-resolution simulations could help clarify the impact of this parameter.

In Fig. 4, we show the difference between the median trends obtained for the high and low settings of two representative parameters: A_SN1 and M_WDM. The associated uncertainties are computed through standard error propagation, based on the total uncertainties estimated for the individual high and low trends, and are displayed as shaded regions. The figure shows results from both the uniform-box simulations (dashed lines) and the MW zoom-in runs (thick lines). Overall, the impact of these parameters is noticeably stronger in the MW zoom-ins. For A_SN1, the difference between the extreme configurations reaches 0.5–0.7 dex, while in the uniform-box simulations the effect remains more modest, around 0.2–0.4 dex. This suggests that feedback processes regulated by A_SN1 are more effective or more easily captured in the zoom-in setup. Regarding the WDM mass, the resulting variations are smaller in magnitude – typically on the order of 0.10-0.15 dex – but still detectable, particularly in quantities such as the ratio M_tot, 1/2/M_*, 1/2 and M_DM, 1/2. While subtle, these shifts imply that WDM mass can leave a non-negligible imprint on the inner structure of galaxies. Although the WDM particle mass exhibits a weak signature on the galaxy scaling relations analyzed here, the high-resolution MW zoom-in simulations reveal a clear additional effect: a systematic reduction in the number of galaxies as a function of total stellar mass with decreasing WDM particle mass. This trend reflects the expectation that WDM suppresses the formation of low-mass sub-haloes (e.g., Lovell et al. 2014; Read et al. 2017), thereby limiting galaxy formation at small scales. DREAMS simulations confirm this behavior: the total number of galaxies decreases as the WDM particle mass decreases, with the suppression becoming increasingly pronounced at lower stellar masses.

Figure 5 illustrates the mass-dependent suppression by showing the median number of galaxies as a function of stellar mass, comparing models with M_WDM < 2 keV and M_WDM > 10 keV. The figure also reports the median trends for different values of the other feedback parameters. The threshold values used to define the parameter subsets are the same as those adopted in this section for the other scaling relations (see Table 1). Among all feedback parameters, A_SN1 has the strongest influence on galaxy abundance, as expected. Its effect manifests primarily as a mass-independent vertical shift in the number of galaxies. In contrast, the WDM particle mass introduces a mass-dependent suppression: the number of galaxies is strongly reduced at log₁₀M_*/M_⊙ ≲ 8, with the effect gradually weakening and becoming negligible at log₁₀M_*/M_⊙ ∼ 9. Regarding the other parameters, BH_FF has no measurable impact on galaxy abundance in this stellar mass range, consistent with AGN feedback becoming relevant only in more massive systems. A_SN2 produces a minor increase in the number of galaxies. The overall sensitivity of galaxy abundance to the feedback parameters and on the WDM mass is qualitatively consistent with the findings of Rose et al. (2025). While the suppression of the stellar mass function by WDM is theoretically expected, our results provide a direct, quantitative realization of this effect within a fully hydrodynamical framework that systematically explores baryonic feedback. The stellar mass dependence of galaxy abundance thus emerges as a promising diagnostic to help disentangle DM properties from feedback processes.

In this section, we have focused on the cumulative median trends of all galaxies in the MW zoom-in simulation, including both the central (host) galaxy and its satellites, with particular attention given to the lowest-mass systems. However, it is also informative to isolate the host galaxies and analyze their properties separately. A detailed discussion of the host galaxy properties – as a function of the astrophysical feedback parameters and the WDM particle mass – is presented in Appendix E, where we find that their overall trends are qualitatively consistent with those of the full galaxy population.

3.2. Evaluation of the fitting procedure

In this section we assess the effectiveness of the fitting procedure described in Sect. 2.3 in recovering astrophysical and cosmological parameters, including the mass of the WDM, through an internal validation process. This methodology closely follows the framework established by Busillo et al. (2023). As a first step, we select a subset of simulations to serve as reference datasets for evaluation – that is, datasets with known input values used as benchmarks to test the quality of the fitting procedure. We then apply the fitting algorithm to infer astrophysical and cosmological parameters as well as the WDM mass. Since the true input values are known in these reference datasets, we can directly evaluate the accuracy of the inference by comparing the recovered parameters to the ground truth. To quantify the level of agreement, we compute the Pearson correlation coefficient and determine its statistical significance⁷. Once a significant correlation between true and inferred parameters is established, we assess whether the relationship is consistent with an ideal one-to-one correspondence, performing a weighted linear fit and computing the slope m and intercept c. The recovery is considered successful if m is statistically compatible with 1 and c with 0.

3.2.1. Parameter recovery in uniform-box simulations

We assess here the performance of the fitting procedure when applied to uniform-box simulations. To ensure consistency, we apply the same selection criteria, used for the comparison with the SPARC catalog and described in Sect. 2.1.2, to all simulations, including the subset serving as the reference. Correlation coefficients and linear fit parameters, are reported in Appendix F.1 (Table F.1).

Figure 6 exhibits the recovery plots for cosmological and astrophysical parameters, as well as for the WDM mass. Strong positive correlations are observed for Ω_m, σ₈, A_SN1, and A_SN2, with statistically significant results at the α = 0.05 level. The linear fits support these findings: for Ω_m, A_SN1, and σ₈, the slope and intercept are statistically consistent with a one-to-one relation, indicating a successful recovery. For A_SN2, the slope is slightly underestimated and the intercept is overestimated.

Fig. 6. Recovery plot for the uniform-box simulations showing the bootstrap results with uncertainties for the estimated parameters: Ωm, σ8, MWDM, ASN1, ASN2, and BHFF against true values. These parameters were obtained by comparing the simulations with the simulation outputs mimicking the SPARC catalog. The dashed gray line indicates the ideal one-to-one relation between estimated and true parameter values.

By contrast, M_WDM and BH_FF exhibit negative correlations with the ground truth and large uncertainties, suggesting that the fitting procedure fails to constrain these parameters. The failure to recover M_WDM is expected, as the median trends show a lack of sensitivity to this parameter (see Fig. 1), whereas for BH_FF, it can be attributed to the sample selection. Indeed, the median trends discussed in Sect. 3.1.1 are computed over the full galaxy population, while here the fitting procedure is restricted to SFGs, which underrepresent massive galaxies where the effect of BH_FF is more prominent. Interestingly, while our analysis does not constrain M_WDM, other approaches have reported sensitivity to WDM masses of several keV (e.g. Rose et al. 2024; Costanza et al. 2025). These constraints rely on methodological choices that differ from ours – such as deep-learning models applied to DM-only density fields or machine-learning emulators trained on galaxy-population statistics – which extract information in ways complementary to the global scaling-relation framework adopted here.

These results are in agreement with Paper I, where a similar analysis on CDM CAMELS simulations identified Ω_m and A_SN1 as the most robustly recovered parameters. Moreover, the parameter σ₈ was reasonably well recovered, although with larger associated uncertainties. In contrast, A_SN2 exhibited less consistent behavior, passing the consistency test in some cases but failing in others, which indicates greater difficulty in obtaining reliable constraints. The DREAMS and CAMELS simulations change different AGN feedback parameters, so a direct comparison is not possible.

3.2.2. Parameter recovery in MW zoom-in simulations

To assess how well the fitting procedure performs when applied to the MW zoom-in simulations, we used simulated galaxy samples filtered according to the selection criteria described in Sect. 2.1.2, ensuring consistency with the SPARC and dwarf galaxy catalogs.

Figure 7 shows parameter estimates with uncertainties from scaling relation fits, compared to the ground truth provided by the simulation. In the comparison with the dwarf-like simulated samples, we consider two cumulative fit cases: one including all the scaling relations, and one excluding the M_tot versus M_* relation. In the figure, we show the results obtained by considering all scaling relations. Correlation coefficients and fit parameters, are reported in Table F.1 in the Appendix F.1.

Fig. 7. Recovery plot for the MW zoom-in simulations. The figure shows the bootstrap results with uncertainties for the estimated parameters: MWDM, ASN1, ASN2, and BHFF (from top to bottom). These parameters were obtained by comparing the simulations with the simulation outputs mimicking SPARC and the dwarf catalog. The dashed gray line indicates the ideal one-to-one relation between estimated and true parameter values.

The parameter A_SN1 passes the correlation test at a significance level of α = 0.05 for both the SPARC-matched and dwarf-matched simulation samples. In both samples, however, the linear fit systematically underestimates the slope and overestimates the intercept. For M_WDM, a strong correlation is observed with the SPARC analog; however, the linear fit deviates significantly from the ideal one-to-one relation, with the slope close to zero, in clear disagreement with the expected value of one. Overall, in both cases considered here, the SPARC and dwarf analogs, these results indicate that we are not able to constrain the WDM particle mass using scaling relations at the current resolution. Finally, both A_SN2 and BH_FF exhibit no significant correlation, consistent with the findings reported in Sect. 3.1.2.

3.3. Comparing simulated and observed scaling relations

After discussing the impact of each parameter on the scaling relations and validating our fitting procedure using simulation output, we now compare the scaling relations derived from the SPARC and dwarf galaxy catalogs (see Sect. 2.2) with those obtained from the WDM simulations. We apply both a standard fitting procedure and the bootstrap method described in Sect. 2.3, to constrain the variable parameters adopted in the simulations.

3.3.1. Fitting uniform-box simulations to observations

In Fig. 8 we present the simulations that best reproduce the observational catalogs, displaying both the trends from the calibrated uniform-box suite and from the MW zoom-in simulations within the same panel. Here, we focus exclusively on the uniform-box runs, leaving the discussion of the MW zoom-in simulations for the next section.

Fig. 8. Best-fit simulations (points) compared to median observational trends (shaded areas) from SPARC (green) and the dwarf catalog (orange). From top to bottom: Stellar half-mass radius, total-to-stellar mass ratio within R*, 1/2, DM mass within R*, 1/2, and total versus stellar mass. Green points indicate the best-fit uniform-box versus SPARC, white points represent the MW zoom-in best matching SPARC, and orange points show the MW zoom-in best matching the dwarf catalog.

The simulated trends (green points) show good agreement with the SPARC observational data (green shaded area). The best-fit parameters are Ω_m = 0.26, σ₈ = 0.99, S₈ = 0.93, A_SN1 = 1.75, A_SN2 = 1.42, BH_FF = 1.45 and ${\mathrm{P}}_{\mathrm{WDM}}=0.38{\mathrm{keV}}^{-1}$ (i.e., M_WDM = 2.6 keV). We include S₈ as a derived parameter, defined as $S_8:={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$. The cumulative reduced D-square is ${\stackrel{\sim }{D}}^2=1.09$, while all individual and cumulative D² and ${\stackrel{\sim }{D}}^2$ values are reported in Table G.1. The fitting results indicate significantly higher values for A_SN1, A_SN2, σ₈, and S₈ compared to those reported in Paper I, where A_SN1 = 0.48, A_SN2 = 1.24, σ₈ = 0.83, and S₈ = 0.78, while the matter density is similar: Ω_m = 0.27. We remind the reader that Paper I, which analyzes CDM CAMELS simulations, included the same cosmological and astrophysical parameters, except for AGN feedback and the WDM mass.

Bootstrap resampling analysis further supports the results of the previous fit; the parameter estimates and their associated uncertainties are reported in Table G.3. In the text below, we also report the corresponding percentage uncertainties, computed relative to the full prior range of each parameter – that is, each uncertainty is expressed as a fraction of the allowed variation range. The results are as follows: ${\mathrm{\Omega }}_{\mathrm{m}}=0.{28}_{-0.04}^{+0.08}$ (±10/20%), ${\sigma }_8=0.{89}_{-0.18}^{+0.10}$ (±25/45%), $S_8=0.{85}_{-0.12}^{+0.09}$ (±10/13%), $A_{\mathrm{SN}1}=1.{67}_{-0.84}^{+0.49}$ (±13/22%), $A_{\mathrm{SN}2}=1.{46}_{-0.21}^{+0.32}$ (±14/21%), ${\mathrm{BH}}_{\mathrm{FF}}=1.{40}_{-0.66}^{+0.59}$ (±16/18%), and ${\mathrm{P}}_{\mathrm{WDM}}=0.{37}_{-0.13}^{+0.10}{\mathrm{keV}}^{-1}$ (±20/27%). Applying the inverse transformation, M_WDM = 1/P_WDM, we obtain a corresponding WDM particle mass of $M_{\mathrm{WDM}}=2.7_{-0.6}^{+1.3}$ keV. Since this transformation is nonlinear, the propagated uncertainties become asymmetric. In this case, it is more appropriate to express the uncertainty relative to the central value (i.e., 21/54%) rather than to the range of variation of the original parameter. The reduced D-square from the bootstrap is ${\stackrel{\sim }{D}}^2=5.{08}_{-0.36}^{+0.32}$. Next, we compare these results with those from Paper I, noting that, at these resolutions, the WDM mass has no significant influence on the scaling relations analyzed in our study. While the cosmological parameters are consistent within 1σ with those recovered in Paper I–${\mathrm{\Omega }}_{\mathrm{m}}=0.{27}_{-0.05}^{+0.01}$, ${\sigma }_8=0.{83}_{-0.11}^{+0.08}$, and $S_8=0.{78}_{-0.09}^{+0.03}$–the supernova feedback parameters are considerably higher with respect to Paper I, as they found $A_{\mathrm{SN}1}=0.{48}_{-0.16}^{+0.25}$ and $A_{\mathrm{SN}2}=1.{21}_{-0.34}^{+0.03}$. Furthermore, the parameter estimates here are generally less precise with respect to Paper I. These discrepancies are attributed to the introduction of the BH_FF parameter, which has a more significant impact in this study compared to the AGN feedback parameters used in Paper I (see Sect. 3.1). The higher inferred values of A_SN1 and A_SN2 compared to Paper I, and of BH_FF relative to its fiducial value of 1, likely reflect their distinct effects on the model’s scaling relations: A_SN1 tends to shift them upward, while BH_FF generally shifts them downward, with a dependence on stellar mass. The effect of A_SN2 is more intricate, but when focusing on SFGs, its impact resembles that of BH_FF: higher values of A_SN2 are associated with smaller galaxy sizes and lower enclosed quantities, although this does not hold for the total stellar mass. Therefore, instead of converging to very low values of A_SN1 to fit SPARC data, as seen in Paper I, the fit identifies a new combination of parameters that provides a good match to the observations. Moreover, the inclusion of the calibration procedure – introduced to mitigate the tendency of simulations to overestimate the scaling relations – prevents the fit from favoring artificially low values of A_SN1, which would otherwise shift the predicted trends downward. Indeed, in the absence of calibration, the simulations yield cosmological parameter values and A_SN1 that are systematically lower, as shown in Table G.4: ${\mathrm{\Omega }}_{\mathrm{m}}=0.{26}_{-0.04}^{+0.04}$, ${\sigma }_8=0.{87}_{-0.16}^{+0.11}$, $S_8=0.{80}_{-0.12}^{+0.11}$, and ${\mathrm{A}}_{\mathrm{SN}1}=1.{04}_{-0.56}^{+0.78}$. Compared to the calibrated runs, we also find a lower value of ${\mathrm{A}}_{\mathrm{SN}2}=1.{27}_{-0.48}^{+0.47}$, and a higher value of ${\mathrm{BH}}_{\mathrm{FF}}=1.{87}_{-0.95}^{+1.38}$. This helps explain the broader uncertainties seen in the results, which are likely due to stronger parameter degeneracies caused by the influence of BH_FF. In the future, a dedicated study incorporating PGs, where the impact of BH_FF and A_SN2 is more significant, will be performed, possibly helping to break these degeneracies and lead to tighter constraints on the model parameters.

3.3.2. Fitting MW zoom-in simulations to observations

In this section, we compare the MW zoom-in simulations with the observational catalogs from SPARC and dwarf galaxies. Before going into the details of the analysis, a few considerations are in order.

The limited number of galaxies involved reduces the statistical robustness of the fit. We computed the median number of galaxies used in the comparison between simulations and observations, along with the 16th and 84th percentiles. Taking into account the stellar mass range covered by the catalogs and applying a possible cut in sSFR, we compare ${12}_{-7}^{+14}$ galaxies with SPARC and ${10}_{-6}^{+11}$ galaxies with the dwarf catalog. When the sample size is very small, some simulations may appear to agree well with the data due to statistical fluctuations rather than an accurate description of the underlying physics. Nevertheless, we retain such cases to avoid introducing selection biases. Importantly, a low number of galaxies can sometimes be associated with strong supernova feedback – particularly high values of the A_SN1 parameter, as shown in Fig. 5. However, this is not a strict correlation: some simulations with weak feedback also yield few galaxies, and vice versa. Moreover, within the stellar mass range probed by the dwarf catalog (approximately log₁₀M_*/M_⊙ ∼ 7.4 to 8), the actual number of satellite galaxies observed in the MW is only about two. As such, models producing few galaxies in this regime remain physically plausible and cannot be excluded from the analysis on the basis of low statistics alone.

As shown in Sects. 3.1 and 3.2, our analysis of the MW zoom-in simulations allows us to constrain only the feedback parameter A_SN1. Although these simulations contain a relatively small number of galaxies, the fixed cosmological background reduces parameter degeneracies, thereby improving the robustness of the constraints on A_SN1. Next, we present our results, starting with the comparison with the SPARC catalog, and then compare them with those obtained from the uniform-box simulations.

In Fig. 8 we show the trends from the best-fit simulation (white points) compared to the SPARC catalog. The best-fit simulation parameter configuration is the following: A_SN1 = 0.85, A_SN2 = 0.62, BH_FF = 0.67, and ${\mathrm{P}}_{\mathrm{WDM}}=0.22{\mathrm{keV}}^{-1}$ (i.e., M_WDM = 4.50 keV), yielding a reduced D-squared value of ${\stackrel{\sim }{D}}^2$ = 1.27 (see Table G.1). The fit obtained via bootstrap resampling is consistent with our initial parameter estimates, with the exception of the BH_FF. In fact, we find that (Table G.2) $A_{\mathrm{SN}1}=1.{31}_{-0.52}^{+0.83}(\pm 14/22\%)$, $A_{\mathrm{SN}2}=0.{89}_{-0.32}^{+0.67}(\pm 21/45\%)$, ${\mathrm{BH}}_{\mathrm{FF}}=1.{08}_{-0.72}^{+1.23}(\pm 19/33\%)$, and ${\mathrm{P}}_{\mathrm{WDM}}=0.{31}_{-0.20}^{+0.18}{\mathrm{keV}}^{-1}(\pm 37/40\%)$, which corresponds to M_WDM = 3.56 keV ( ± 37/181%), with ${\stackrel{\sim }{D}}^2=2.{71}_{-0.70}^{+0.63}$. This estimate of A_SN1 is consistent with that derived from the uniform-box analysis, with comparable uncertainties.

The best-fit simulation with respect to the dwarf catalog is shown in Fig. 8 (orange points). The parameters of the best-fit simulation are: A_SN1 = 0.74, A_SN2 = 1.84, BH_FF = 0.56 and P_WDM = 0.39 keV⁻¹ (i.e., M_WDM = 2.55 keV), with a reduced D-squared value of ${\stackrel{\sim }{D}}^2$ = 0.96. As usual, the individual D-squared values for the different scaling relations are listed in Table G.1. The bootstrap analysis shows (Table G.2): $A_{\mathrm{SN}1}=0.{68}_{-0.34}^{+0.50}(\pm 9/13\%)$, $A_{\mathrm{SN}2}=0.{91}_{-0.34}^{+0.77}(\pm 23/51\%)$, ${\mathrm{BH}}_{\mathrm{FF}}=0.{78}_{-0.46}^{+1.14}(\pm 12/30\%)$, and ${\mathrm{P}}_{\mathrm{WDM}}=0.{37}_{-0.17}^{+0.10}{\mathrm{keV}}^{-1}(\pm 20/35\%)$, i.e., M_WDM = 7.04 keV (±22/85%), with ${\stackrel{\sim }{D}}^2=2.8_{-0.9}^{+1.7}$. The estimate of A_SN1 obtained from the comparison with the dwarf galaxy catalog is more precise than those derived from the SPARC catalog, uniform-box, and MW zoom-in simulations. However, the dwarf catalog tends to favor models with lower A_SN1 values compared to the overall results from the SPARC analysis. Despite these differences, the estimates remain consistent within their respective uncertainties.

4. Conclusions

In this work, we extended the CASCO project by analyzing DREAMS hydrodynamic simulations in a WDM cosmology (Rose et al. 2025) and comparing them to observational data from the SPARC and LVDB catalogs. Our goal was to evaluate how cosmological and astrophysical parameters, especially the WDM particle mass, affect galactic scaling relations. We considered two complementary simulation suites – uniform-box and MW zoom-in runs – characterized by different numerical resolutions and volumes, focusing on the small-scale signatures of WDM. We began by analyzing how variations in simulation parameters affect the median trends of the scaling relations between stellar mass and four key quantities: the stellar half-mass radius, the enclosed DM mass, the enclosed total-to-stellar mass ratio, and the total mass. We then constrained simulation parameters by comparing DREAMS WDM outputs to observations from SPARC and LVDB by using a bootstrap-based fitting method developed in Papers I and II. This approach was validated across both simulation types for its ability to constrain cosmological, astrophysical, and WDM parameters. Below, we summarize our main results:

• Simulated scaling relations in the uniform-box simulations, where cosmological parameters were varied, demonstrate a clear dependence on the matter density parameter. Higher values of Ω_m systematically increase the median values across all four scaling relations. In contrast, σ₈ has milder effects: It reduces stellar sizes and total masses while increasing the DM content and total-to-stellar mass ratio below log₁₀M_*/M_⊙ ∼ 10.8, where the trend reverses.

• Astrophysical parameters were varied in both simulation suites, allowing for a comparative analysis. The SN wind energy (i.e., A_SN1) has the most significant impact: Increasing A_SN1 leads to a systematic rise in galaxy sizes, higher total-to-stellar mass ratios, and DM masses within the stellar half-mass radius, as well as an overall increase in the total galaxy mass. This effect is observed in both the uniform-box and MW zoom-in simulations, but it is more pronounced in the latter, where differences between the high and low median trends reach 0.5–0.7 dex, compared to 0.2–0.4 dex in the uniform-box suite. The SN wind velocity (A_SN2) and black hole feedback factor (BH_FF) generally have weaker effects than A_SN1 and primarily affect galaxies with higher stellar masses. As a result, the MW zoom-in simulations, which include mostly low-mass galaxies, are less sensitive to these parameters. The BH_FF parameter is found to suppress all four scaling relations at the high-mass end, as consistently observed across both simulation suites.

• Regarding the WDM particle mass, previous studies have shown that simulations at resolutions comparable to the DREAMS uniform boxes can retain constraining power on WDM models up to ∼6 keV, particularly when analyzed using population-level statistics or machine-learning techniques (e.g., Rose et al. 2024; Costanza et al. 2025). In our analysis based on scaling relations that are expected to be most sensitive to the underlying DM physics, we find no significant differences in the uniform-box simulations across the explored WDM mass range. By contrast, slight trends in M_tot, 1/2/M_*, 1/2 and M_DM, 1/2 emerge in the MW zoom-in simulations, suggesting that higher-resolution simulations could offer valuable insights and represent a promising direction for future study.

• Moreover, although the MW zoom-in and uniform-box simulations exhibit similar overall behaviors, a direct comparison between the MW zoom-in trends and those of the calibrated uniform-box suite also indicates that high-resolution simulations produce galaxies with smaller stellar half-mass radii as well as lower values of M_DM, 1/2 and M_tot, 1/2/M_*, 1/2. The total mass versus stellar mass trend, however, remains smoothly aligned across both simulation types. These systematic shifts are likely driven by resolution effects and environmental factors, among other influences, that can impact the underlying astrophysical processes.

• We also analyzed galaxy abundance as a function of total stellar mass, showing that this relation may help disentangle the effects of WDM particle mass from those of other parameters due to its strong dependence on total stellar mass.

• We confirm the presence of a golden mass, a characteristic stellar mass scale where star formation efficiency peaks, in the uniform-box WDM simulations, which show more clearly defined minima compared to CAMELS CDM simulations at the same resolution.

• Internal validation of the bootstrap-based fitting strategy was performed by accounting for the characteristics of the observational catalogs used for comparison. While this approach does not reflect the method’s full potential when applied to the two simulation suites, it offers a realistic estimate of its ability to constrain the model parameters. In the uniform-box suite, the most robustly constrained parameters are Ω_m and A_SN1. The parameter σ₈ is reasonably well recovered, though with larger uncertainties, while A_SN2 passes the correlation test but fails to reproduce an ideal one-to-one relation with the ground truth. In the MW zoom-in simulations, A_SN1 shows significant correlation with the target values in both the SPARC-analog and dwarf-analog cases, though the one-to-one relation is only closely approached in the former. This reduced performance is likely due to the limited statistics of the MW zoom-in sample, which hinders tighter constraints. Finally, the WDM particle mass cannot be constrained.

• Comparison between the simulations and observational catalogs yielded the following results. In the uniform-box simulations, ${\mathrm{\Omega }}_{\mathrm{m}}=0.{28}_{-0.04}^{+0.08}$ is in agreement with the value reported for the CAMELS simulations analyzed in Paper I ($0.{27}_{-0.05}^{+0.01}$), with uncertainties at the 10–20% level. The estimate for ${\sigma }_8=0.{89}_{-0.18}^{+0.10}$ is slightly higher than in the earlier analysis ($0.{83}_{-0.11}^{+0.08}$), though still compatible within 1σ, resulting in a comparable S₈. The feedback parameter A_SN1 is significantly higher than in Paper I ($1.{67}_{-0.84}^{+0.49}$ versus $0.{48}_{-0.16}^{+0.25}$), although it is more loosely constrained. In contrast, $A_{\mathrm{SN}2}=1.{46}_{-0.21}^{+0.32}$ shows a milder increase compared to the previous estimate ($1.{21}_{-0.34}^{+0.03}$). This discrepancy may stem from the fact that the CAMELS and DREAMS simulations change different AGN feedback parameters–resulting in negligible effects in the former and more significant impacts in the latter–making a direct comparison unfeasible. For the MW zoom-in simulations, only A_SN1 is meaningfully constrained. Its estimates – $1.{31}_{-0.52}^{+0.83}$ from the SPARC catalog and $0.{68}_{-0.34}^{+0.50}$ from the dwarf catalog – are consistent with the uniform-box result of $1.{67}_{-0.84}^{+0.49}$. The dwarf catalog provides the tightest constraint, although it tends to prefer lower values, closer to the estimate from Paper I.

Our analysis highlights the value of combining simulations with different resolutions and observational catalogs that probe complementary stellar mass ranges. While we find that the current constraints are insufficient to determine the WDM particle mass, our results in a WDM framework reveal subtle mass-dependent effects, particularly in the galaxy abundance as a function of stellar mass, that may help disentangle the influence of WDM from other parameters. This suggests that a multiscale approach, leveraging both zoom-in and uniform-box simulations, is a promising avenue for future studies aiming to constrain the physics of structure formation in nonstandard cosmologies. In future work, we plan to extend our analysis by examining the DM density profiles of different components and by quantifying enclosed masses within radii smaller than the effective radius (r < R_*, 1/2). These inner regions are expected to be more sensitive to the suppression of small-scale power induced by WDM, especially in low-mass (dwarf) galaxies. Furthermore, we will investigate whether the core-like structures predicted in DM-only WDM simulations (Macciò et al. 2012, 2013) are preserved or altered in the presence of baryonic feedback. We also aim to explore the redshift and environmental dependence of these trends, which may further enhance the discriminatory power between WDM and CDM models. In this context, upcoming large-scale surveys, such as Euclid (Marleau et al. 2025; Euclid Collaboration: Marleau et al. 2026), with their ability to map the faint end of the luminosity function with unprecedented precision, will be instrumental. When combined with future improvements in hydrodynamic simulations – including enhanced resolution and refined baryonic physics – these observational datasets will allow us to better probe the role of WDM and its viability as an alternative to CDM.

Acknowledgments

All the calculations underlying this work have been performed via the use of Wolfram Mathematica ver. 13.1. A.M.B. acknowledges support from grant FI-CCA-Research-00011826 from the Simons Foundation. A.F. acknowledges support from the National Science Foundation under Cooperative Agreement 2421782 and the Simons Foundation award MPS-AI0001051. N.R.N. acknowledges support from the Guangdong Science Foundation grant (ID: 2022A1515012251). P.T. acknowledges support from NSF-AST 2346977 and the NSF-Simons AI Institute for Cosmic Origins which is supported by the National Science Foundation under Cooperative Agreement 2421782 and the Simons Foundation award MPS-AI-00010515.

References

Appendix A: Calibration of the uniform-box simulation suite with respect to TNG100-1

In this section, we describe the calibration procedure adopted for the lower-resolution uniform-box DREAMS simulations, based on the fiducial TNG model. This calibration is essential to reproduce observational trends more accurately.

Numerical resolution significantly impacts the effectiveness of astrophysical feedback processes. Lower resolution implies more massive gas particles, leading to a less accurate representation of the gas distribution. This limits the ability to resolve star formation, which occurs in dense, small-scale environments. Furthermore, feedback models, such as those associated with SN energy release, may become ineffective or poorly modeled, as the energy is distributed over excessively large masses. Higher resolution enables better physical modeling and more efficient feedback implementation. To account for these limitations, we calibrated the WDM uniform-box simulations of DREAMS using the TNG100-1 fiducial run as a reference. TNG100-1 features a DM particle mass of 5.1 × 10⁶ h⁻¹ M_⊙, a baryonic mass resolution of 9.4 × 10⁵ h⁻¹ M_⊙, and a comoving box size of 75h⁻¹Mpc per side. As no WDM uniform-box simulation is available for direct comparison, we use the 1P ⁸ WDM 16 keV simulation as a proxy. These simulations are specifically designed to isolate the impact of individual parameters by varying one at a time while keeping all others fixed. In this case, the only difference from the reference setup is the WDM particle mass. At this resolution, we verified that the WDM mass parameter has a negligible impact on the relevant physical processes (see Sect. 3.1). The calibration is based on the stellar-to-DM mass relation (M_* versus M_DM) and proceeds as follows:

• TNG100-1 interpolation: we compute the median, 16th, and 84th percentiles of the M_* versus M_DM relation from TNG100-1 and generate an interpolating function of M_* as a function of M_DM.

• Cost function definition: for each object in the WDM simulation (1P, 16 keV), we evaluate the squared difference between the stellar mass (shifted by a constant c_M* in log) and the interpolated TNG100-1 prediction, normalized by the TNG100-1 uncertainty (average between 16th and 84th percentiles) at the same M_DM.

• Filtering: objects outside the interpolation validity range are excluded.

• Minimization: we determine the optimal logarithmic shift c_M* by minimizing the sum of the normalized residuals. The resulting offset is c_M* ∼0.27 in log₁₀M_*/M_⊙.

Finally, we rescale stellar-mass-dependent astrophysical quantities as follows:

• [Inline555${\text{log}}_{10}{M}_{*}/M_{\odot }\to {\log }_{10}{M}_{*}/M_{\odot }+c_{M_{*}}$;

• ${\text{log}}_{10}{M}_{\text{tot}}/M_{\odot }\to {\log }_{10}(({10}^{c_{M_{*}}}{M}_{*}+M_{\text{DM}}+M_{\text{gas}}+M_{\text{BH}})/M_{\odot }$);

• $f_{\text{DM},1/2}\to \frac{M_{\text{DM},\text{1/2}}}{{10}^{{\text{c}}_{{\text{M}}_{*}}}{M}_{*,1/2}+M_{\text{DM},1/2}+M_{\text{BH},1/2}+M_{\text{gas},1/2}}$.

Figure A.1 shows the stellar-DM mass relation before and after applying the calibration. While this procedure effectively mitigates resolution effects on the M_* versus M_DM relation, it is important to note that the calibration is performed a posteriori and tailored to the TNG100-1 reference run. Nevertheless, it consists solely of a rigid shift in stellar mass and does not involve tuning of the simulation parameters themselves. Therefore, it does not trivially favor the reference configuration when analyzing other scaling relations. Indeed, the effects of physical parameters in the simulations are generally non linear and mass-dependent, whereas the applied correction is uniform across all stellar masses. In Fig. A.2, we present the median values together with the 16th and 84th percentiles for the same set of scaling relations adopted in our analysis, for both the reference simulation TNG100-1 and the 1P run at 16 keV, before and after the stellar-mass calibration. We note that the calibration described in this section brings the DREAMS uniform-box simulation into much closer agreement with the reference model, with a substantial overlap between the two. This indicates that a systematic shift in M_* has a generally beneficial effect on the other reference scaling relations as well.

Fig. A.1. Stellar-DM mass relation for the WDM uniform-box 1P DREAMS simulation with a 16 keV WDM particle mass (blue), its calibrated counterpart (yellow), and the TNG100-1 fiducial run (light gray), used as reference.

Fig. A.2. Median scaling relations with 16th and 84th percentiles for the reference simulation TNG100-1 (light gray), the uncalibrated 1P run at 16 keV (blue), and the stellar-mass–calibrated 1P run (yellow).

Appendix B: Calibration of the uniform-box simulation suite with respect to the MW zoom-ins

Fig. B.1. Stellar–DM mass relation for the WDM uniform-box 1P DREAMS simulation with a 16 keV WDM particle mass (blue), its stellar-mass–calibrated counterpart (yellow), and the corresponding reference samples (gray). The upper panel shows the calibration relative to the satellite galaxies from the MW zoom-in simulations, while the lower panel shows the calibration relative to TNG50-1.

In this section we extend the analysis presented in Appendix A by performing a two-step calibration of the standard uniform-box 1P simulation with a 16 keV WDM particle mass, using as reference the population of satellite galaxies extracted from the 1024 MW zoom-in simulations.⁹ We refer to this procedure as a double calibration, as it is carried out sequentially on the stellar mass, M_*, and on the stellar half-mass radius, R_*, 1/2.

B.1. Stellar-mass calibration

We calibrate the stellar mass with respect to the satellites of the MW zoom-in simulations, following the same procedure outlined in Appendix A. We find a best-fitting offset of c_M* = 0.32 in log₁₀M_*/M_⊙ (see upper panel of Fig. B.1), close to the value obtained in Appendix A (i.e., c_M* = 0.27). After applying this correction, the median scaling relations of the calibrated uniform-box simulation (yellow) remain systematically above those of the MW zoom-in satellites (light gray), especially for relations involving the stellar half-mass radius R_*, 1/2 (see Fig. B.2). This residual discrepancy may arise from a combination of factors, including resolution effects that can limit the modeling of small-scale baryonic structure, as well as physical or environmental differences between the two simulation suites.

To further test the impact of resolution, we repeated the calibration using the TNG50-1 run as reference. TNG50-1 has a spatial and mass resolution comparable to that of the MW zoom-in simulations, and therefore provides an appropriate reference for this comparison. Using TNG50-1 as reference, we calibrated the stellar mass of the uniform-box run (see lower panel of Fig. B.1); the resulting shift, c_M* = 0.45, is slightly higher but of the same order as those obtained for TNG100-1 and the MW zoom-ins. The resulting scaling relations of the calibrated uniform-box simulation (red) closely follow those of TNG50-1 (dark gray; see Fig. B.3). It is plausible that environmental factors, inherent to the MW zoom-in sample of MW-mass haloes, also contribute to the residual offsets at fixed stellar mass. To further explore this aspect, we examined the uniform-box runs and found that regions of higher local density–for example, spheres of 350 kpc radius containing at least ten neighboring galaxies–tend to host systems with lower total mass at fixed stellar mass. Although only indicative, this behavior points to a possible environmental contribution to the differences observed between the two simulation suites.

B.2. Stellar half-mass radius calibration

Building on the stellar mass correction (c_M* = 0.32; see Appendix B.1), we identified residual offsets in the size–mass relation and in quantities measured within R_*, 1/2. To mitigate these, we used the R_*, 1/2 versus M_DM relation from the MW zoom-in satellites (see Fig. B.2), which yields a logarithmic shift of c_R = −0.30, corresponding to:

$$\begin{array}{c}\hfill {log}_{10}{R}_{\ast ,1/2}/\mathrm{kpc}\to {log}_{10}{R}_{\ast ,1/2}/\mathrm{kpc}+c_{\mathrm{R}}.\end{array}$$(B.1)

A calibration on R_*, 1/2 implies corresponding adjustments to all quantities defined within this radius. The stellar mass within R_*, 1/2 is easily corrected, since $M_{\ast ,1/2}\equiv \frac{1}{2}{M}_{\ast }$ and therefore follows the same scaling as M_*. For the DM mass within the stellar half-mass radius, we assume a simple power-law dependence around R_e, M_DM ∝ r^η, where η = 2 corresponds to an NFW profile and η = 1.2 to a contracted NFW profile (Boylan-Kolchin et al. 2005; Tortora et al. 2018). Accordingly,

$$\begin{array}{c}\hfill {log}_{10}{M}_{\mathrm{DM},1/2}/M_{\odot }\to {log}_{10}{M}_{\mathrm{DM},1/2}/M_{\odot }+\eta c_R.\end{array}$$(B.2)

Throughout the analysis we adopt η = 2, noting that the differences with a contracted NFW assumption are marginal. Finally, when computing the total mass enclosed within R_*, 1/2, we do not rescale the black hole and gas masses, as the resulting error is negligible. For completeness, we note that recalibrating these components would have only a minor effect, since their median contributions to the total mass within R_*, 1/2 in the simulations are about 0.2% for the black hole and 6% for the gas. The 1P simulation calibrated on both M_* and R_*, 1/2 is represented by the orange trends in Fig. B.3, which closely follow those of the MW zoom-in satellites (light gray trends). For comparison with the other results in this paper, we performed a bootstrap fit of the uniform-box simulations to the SPARC catalog after the double calibration. In this case, the resulting A_SN1 value is significantly higher, around 3.5. This occurs because the double calibration causes the uniform box trends to be systematically underestimated, requiring a higher A_SN1 to raise the trends and match the SPARC data.

Fig. B.2. Stellar half-mass radius–DM mass relation for the WDM uniform-box 1P DREAMS simulation with a 16 keV WDM particle mass (blue), its calibrated counterpart (yellow), and the satellite galaxies from the MW zoom-in simulations used as reference (gray). The equivalent plot for the stellar-mass calibration is shown in the left panel of Fig. B.1.

Fig. B.3. Median scaling relations with 16th and 84th percentiles for the MW zoom-in satellites (light gray), TNG50-1 (dark gray), the uncalibrated 1P run (blue), the 1P run calibrated on stellar mass with respect to the MW zoom-ins (yellow), the same run further calibrated on stellar half-mass radius (orange), and the 1P run calibrated on stellar mass with respect to TNG50-1 (red).

Appendix C: Determination of uncertainties on simulated galaxy quantities

In this section, we evaluate the typical uncertainties of relevant simulated galaxy quantities (e.g., effective radius) as a function of galaxy stellar mass. We apply this procedure for the galaxies in both MW zoom-in and uniform-box simulations, in order to evaluate quantitatively how the low particle number affects the knowledge of these quantities.

For this analysis, we both assume that light follows mass and spherical symmetry. We consider a galaxy as a collection of randomly extracted point particles from the following deprojected Sérsic distribution (e.g., Prugniel & Simien 1997):

$$\begin{array}{c}\hfill L(<x)=L_{\mathrm{tot}}\frac{\mathrm{\Gamma }(n(3-p(n)),0,x^{1/n})}{\mathrm{\Gamma }(n(3-p(n)))},\end{array}$$(C.1)

where p(n)≊1 − 0.6 n⁻¹ + 0.05 n⁻², n is the Sérsic index, L_tot is the total luminosity of the galaxy,

$$\begin{array}{c}\hfill x=\frac{r}{R_{\mathrm{e}}}{b}^n(n),\end{array}$$(C.2)

is an adimensional radius (with b(n)≊2 n − 1/3 + 9.9  × 10⁻³ n⁻¹ and R_e projected effective radius), Γ(z₀, z₁, z₂):=Γ(z₀, z₁)−Γ(z₀, z₂) is the generalized incomplete Gamma function, Γ(z₀, z₁) the incomplete Gamma function, and Γ(z) is the Euler Gamma Function. The deprojected effective radius can thus be obtained by solving L(< x) = 0.5 for x and then substituting the obtained value in Eq. (C.2).

To find the uncertainties on the deprojected effective radius, we start by considering fixed values for the projected effective radius and Sérsic index equal to log₁₀(R_e/kpc) = 0.4 (which is roughly equal to the median of the dwarfs’ effective radii in Fig. 3) and n = 1 (noting that the Sérsic index of dwarf galaxies is typically between 0.5 and 1.5, see for example Lazar et al. (2024, Fig. 13) while for the SFGs is ∼1, Ward et al. 2024).

In the case of the uniform-box simulations, for example (the procedure is analogous for the MW zoom-in simulation), we select four different fixed total numbers of stellar mass particles (N_par = [40, 150, 600, 2000]), which correspond to a mass range between log₁₀M_*/M_⊙ ∼ 8.8 and 10.6. For each value of N_par, we produce N_gal = 400 random galaxies. For all these galaxies, we then compute the respective normalized growth curve, obtained by counting the number of particles within a generic radius r and then dividing the count by the total number of particles. We can then estimate the median deprojected effective radius for a galaxy with number of particles equal to a given N_par by finding the empirical deprojected effective radius (i.e., the nearest r value such that the normalized growth curve is equal to 0.5) for each of the N_gal randomly produced galaxies, and evaluating the median and quantiles of these radii. The trend of the effective radii estimated from this procedure as a function of stellar mass is shown in Fig. C.1. As expected, the uncertainties are larger when the number of stellar mass particles is lower. The maximum error on the effective radius, when the total number of stellar mass particles is 40, does not exceed 0.05 dex.

Fig. C.1. Estimate and uncertainty of the half-light radius for galaxies in the DREAMS uniform-box simulations, computed for four different fixed total numbers of stellar mass particles (Npar = [40, 150, 600, 2000]). The dashed red line marks the true half-light radius.

Given the uncertainty on the deprojected effective radius, one can determine the uncertainties on all other derivated quantities: for example, the uncertainty on the stellar mass within the deprojected effective radius can be obtained by considering for each simulated galaxy the number of particles within the empirical deprojected effective radius, and then converting to stellar mass by multiplying the results for the resolution of the star particles. Uncertainty on the DM mass within the deprojected effective radius can then be obtained by considering the M_DM, 1/2/M_*, 1/2 values corresponding to a certain number of particles, and proceed by error propagation by multiplying this ratio for the estimated M_*, 1/2. Similarly, the uncertainties on the total mass within the deprojected half-light radius can be obtained from the expected ratio, M_tot, 1/2/M_*, 1/2, for a given number of particles. Since we already have an estimate of the uncertainty on M_*, 1/2, we can then propagate this to determine the uncertainty on the total-to-stellar half-mass ratio. We observe that the propagated uncertainties on M_DM, 1/2 and on M_tot, 1/2/M_*, 1/2 are also very small and do not exceed 0.04 dex. The uncertainty estimates derived in this section are incorporated in the analysis presented in Sect. 3.1, where they are used to propagate errors on quantities computed within the stellar half-mass radius. In that section, the resulting uncertainties are combined in quadrature with the statistical uncertainty on the medians and shown as shaded regions in Figs. 1, 2, and 3.

Appendix D: Analysis of the golden mass in DREAMS and CAMELS

Fig. D.1. Median trends of the total-to-stellar half-mass ratio as a function of total stellar mass, for different cosmological and astrophysical parameters in the DREAMS (left, thick lines) and CAMELS (right, dashed lines) simulations. Colored curves correspond to subsets of simulations grouped by parameter value; the legend indicates the median value for each group.

In this section, we investigate the so-called golden mass, the characteristic stellar mass scale at which the ratio M_tot, 1/2/M_*, 1/2 reaches a minimum. This minimum marks the peak of star formation efficiency and a transition in galaxy formation physics (Tortora et al. 2025). The analysis presented in this work is based on the DREAMS and CAMELS uniform-box simulations, which have a resolution comparable to that of TNG300 within a volume of (25 h⁻¹Mpc)³. These results are then compared to other CAMELS simulations with the same resolution but covering a larger volume of (50 h⁻¹Mpc)³, as presented in Paper III.

We note that, unlike Paper III, which defines the golden mass using the ratio M_DM, 1/2/M_*, 1/2, we adopt M_tot, 1/2/M_*, 1/2 throughout our analysis. This choice includes the contributions from baryons in the numerator, yielding slightly higher values for the ratio. However, as explicitly stated in footnote 4 of Paper III, the effect of using the total mass instead of the DM mass within the same aperture leads to an increase in the golden mass of at most 0.05–0.1 dex and does not affect the physical conclusions.

Figure D.1 illustrates how the total-to-stellar half-mass ratio responds to variations in cosmological and astrophysical parameters, with results shown for DREAMS (left) and CAMELS (right). Since the parameter distributions differ between DREAMS and CAMELS, the comparison between CDM and WDM is not one-to-one, and differences in the resulting curves are partly driven by sampling choices. Moreover, cosmic variance may also contribute to the observed discrepancies. The golden mass is identified by fitting a parabola to the median trend of the ratio as a function of stellar mass. However, not all parameter choices yield a clear minimum. In several cases, especially for extreme parameter values, the trend is either monotonic or only mildly curved, indicating that the golden mass is not always well-defined. When the minimum is identifiable, its location varies systematically with the physical or cosmological parameter under consideration. For the DREAMS simulations, minima are found around log₁₀M_*/M_⊙∼11.1 for variations in both Ω_m and σ₈, and approximately 11.0 for moderate values of A_SN1. Higher values of A_SN1 often do not show a detectable minimum. A_SN2 shows a gradual shift in the minimum, ranging from 11.2 (stronger winds) to 10.8 (weaker winds). AGN feedback, parametrized by BH_FF, produces minima between 10.8 and 11.2 when present. In the CAMELS simulations, the minima are located at even higher stellar masses than in Paper III. The golden mass occurs around log₁₀M_*/M_⊙∼11.2 for varying Ω_m, and between 11.3 and 11.4 for changes in σ₈. For A_SN2, the minimum shifts from 11.4 at weak feedback to 10.6 for the strongest winds. In the case of A_SN1, when a minimum is present, it lies between 11.4 and 11.5.

Compared to Paper III, which consistently reports a golden mass near log₁₀M_*/M_⊙∼10.6, the values found in both DREAMS and CAMELS uniform boxes are systematically higher. Our analysis confirms the presence of a golden mass in both CDM and WDM scenarios, with trends that qualitatively agree with those reported in Paper III. Interestingly, WDM models tend to show more regular and well-defined minima in M_tot, 1/2/M_*, 1/2, possibly reflecting a smoother buildup of central mass in the absence of small-scale structure. Despite the differing parameter distributions in DREAMS and CAMELS the median trends consistently highlight the golden mass as a robust feature of galaxy formation, shaped by both feedback and the DM model. The higher golden masses observed in smaller-volume simulations relative to Paper III is likely due to volume-dependent sampling effects. Larger volumes sample intermediate- and high-mass galaxies more thoroughly, allowing a more accurate determination of the minimum in M_tot, 1/2/M_*, 1/2 at slightly lower stellar masses. This may explain the systematic shift without invoking differences in underlying physics.

Appendix E: Properties of host galaxies in the MW zoom-in simulations

In this section, we summarize the main properties of the MW-like host galaxies from the MW zoom-in simulations of the DREAMS project. These simulations are constructed by selecting DM halos with masses comparable to that of the MW from a larger cosmological volume, which are then re-simulated at higher resolution using a zoom-in technique (see Rose et al. 2025). As a result of this setup, each simulation contains a single, central galaxy–the host galaxy–that forms at the center of the high-resolution region. The host galaxies generally exhibit a contamination by low-resolution DM particles that remains below 2.5%, with a median value of 0%. An analysis of the galaxy population in the MW zoom-in simulations shows that the level of contamination systematically increases with distance from the central galaxy, whereas the innermost regions, up to ∼ 1000 kpc (far beyond the virial radius), display a median contamination of 0% and only minimal scatter, increasing rapidly at larger distances. These tests confirm that the impact of low-resolution DM particle contamination is minimal for the host galaxies and the satellite populations considered in our analysis.

Table E.1 presents the minimum, maximum, and median values of the key astrophysical quantities analyzed throughout this work. Since each host galaxy evolves within a different realization of baryonic feedback parameters, its physical properties can vary significantly and often differ from those of the MW itself.

In Fig. E.1, we show the trends of R_*, 1/2, M_DM, 1/2, M_tot, 1/2/M_*, 1/2, M_tot and M_* as functions of the astrophysical parameters and the inverse WDM particle mass. From Fig. E.1 we observe that increasing A_SN1 leads to larger stellar half-mass radii, higher enclosed DM masses, and lower stellar masses, resulting in a steep increase in the ratio M_tot, 1/2/M_*, 1/2. The trends with A_SN2 are non-monotonic: both R_*, 1/2 and M_DM, 1/2 rise mildly up to A_SN2 ∼ 1.1 and decline thereafter, while the stellar mass decreases more steadily. This produces a moderate, but continuous, increase in M_tot, 1/2/M_*, 1/2. Higher BH_FF values are associated with a slight increase in stellar mass, a decrease in M_DM, 1/2, and a mild reduction in R_*, 1/2, leading to a lower mass ratio. No significant variations are observed as a function of the inverse WDM particle mass across any of the plotted quantities. Finally, the total mass of the host galaxies does not appear to be significantly affected by any of the simulation parameters. Overall, the qualitative influence of astrophysical and WDM parameters on the host galaxies in the MW zoom-in simulations is consistent with the trends observed in Sect. 3.1.2 for the broader galaxy population. The only notable difference concerns the total galaxy mass, which shows no clear dependence on simulation parameters for the hosts alone, whereas a significant trend is observed when considering the full galaxy sample.

Fig. E.1. Trends of stellar half-mass radius, total-to-stellar mass ratio within the stellar half-mass radius, DM mass within the stellar half-mass radius, total galaxy mass and total stellar mass for the host galaxies in the MW zoom-in simulations, as a function of astrophysical parameters and the inverse of the WDM mass parameter. The light blue points represent individual galaxies; red lines indicate the median trends and their scatter.

Appendix F: Consistency and robustness tests

Fig. F.1. Recovery plot for the uniform-box simulations, showing the bootstrap results with uncertainties for the estimated parameters: Ωm, σ8, MWDM, ASN1, ASN2 and BHFF against their true values. These parameters are obtained by comparing the simulations with the corresponding simulation outputs mimicking the SPARC catalog. The dashed gray line indicates the ideal one-to-one relation between estimated and true parameter values. Blue points correspond to the reference analysis with N = 100, while grey points show the test performed with N = 10.

In this section, we present a set of complementary validation tests aimed at assessing the overall reliability and robustness of our fitting procedure based on the bootstrap resampling method. We first examine the internal consistency of the recovered parameter correlations (Appendix F.1), and then we perform dedicated tests to verify that the bootstrap approach does not introduce biases or underestimate uncertainties (Appendix F.2).

F.1. Consistency tests on parameter correlations

In this section, we present the results of the correlation tests described in Sect. 3.2. Specifically, Table F.1 reports, for both the uniform-box and MW zoom-in samples, the outcomes of the Pearson correlation test along with the estimated slope m and intercept c for the parameters exhibiting statistically significant (i.e., α < 0.05) correlations.

As discussed in the main text, we consider the fitting procedure described in Sect. 2.3 to be effective for estimating a given parameter when the correlation is significant, the slope m is consistent with 1, and the intercept c is consistent with 0.

F.2. Bootstrap resampling robustness

To further validate the reliability of the bootstrap–based fitting method described in Sect. 2.3, we performed additional tests aimed at verifying that the resampling procedure does not bias parameter recovery or underestimate the associated uncertainties, particularly in the presence of parameter degeneracies. These tests are carried out using the uniform-box simulation suite, and for computational reasons we combine bootstrap estimates obtained with N = 100 (used in Sect. 3.2.1) and with N = 10, which allows us to extend the statistical sampling while keeping the analysis tractable. Although a smaller N increases the scatter in individual realizations, it provides a useful stress test for the stability and statistical validity of the method.

Figure F.1 presents the parameter–recovery plots, similarly to Fig. 6 but including the N = 10 points in gray. The parameters Ω_m and A_SN1 closely follow the one–to–one relation, indicating accurate recovery and no systematic bias. A_SN2 also follows the general trend but with larger dispersion, while σ₈ shows a noisier behavior. The WDM mass parameter and the AGN feedback efficiency are not well recovered, as expected given the limited sensitivity of the adopted observables to these quantities. We notice a mild systematic trend across several parameters: low true values tend to be slightly overestimated, whereas high true values are mildly underestimated. This pattern is consistent with a boundary bias, which occurs when the true parameter lies near the edges of a finite grid, causing the recovered distribution to be compressed toward the center of the sampled range (see, e.g., Davison & Hinkley 1997; Silverman 2018). This behavior reflects the finite coverage of the uniform-box simulations rather than a bias intrinsic to the bootstrap resampling.

The normalized residual distributions, Δ_p = (p_true − p_rec)/σ_p, where p_true denotes the value of the parameter set in the simulation, p_rec is the recovered estimate, and σ_p = (σ₊ + σ₋)/2 is the mean bootstrap uncertainty from the 16th and 84th percentiles, are shown in Fig. F.2 together with their Gaussian fits; the fitted mean (μ) and standard deviation (σ) are indicated in each panel. The parameters Ω_m and A_SN1 exhibit nearly Gaussian residuals (μ ≃ 0, σ ≃ 1), confirming unbiased recovery and realistic uncertainty estimates. A_SN2 and σ₈ display broader distributions (σ > 1), consistent with the larger intrinsic scatter already visible in their recovery plots. In contrast, M_WDM and BH_FF show clearly non-Gaussian and asymmetric residuals (σ ≫ 1, μ ≠ 0), reflecting the limited constraining power of the available observables rather than any bias introduced by the bootstrap procedure. Overall, these results demonstrate that the bootstrap resampling method yields statistically robust and unbiased estimates for parameters effectively constrained by the data, while correctly propagating the increased uncertainty for less constrained ones.

Fig. F.2. Normalized residual distributions (Δp) for all parameters in the uniform-box suite. Each panel shows the empirical distribution (histogram) together with a Gaussian fit (red solid curve); the fitted mean (μ) and standard deviation (σ) are indicated in each panel. Histograms are normalized to unit area using Mathematica’s "PDF" option, allowing direct comparison with a standard normal.

Appendix G: Fit results and parameter constraints

In this section, we present the results of the fits obtained through bootstrap resampling, as discussed in Sect. 3.3.

All Tables

All Figures

Fig. 1. Median scaling relations as a function of total stellar mass illustrating the dependence on cosmological parameters. The relations for stellar half-mass radius, total-to-stellar mass ratio within the stellar half-mass radius, DM mass within the stellar half-mass radius, and total galaxy mass are shown for both uncalibrated (solid line) and calibrated (dashed line) cases. Trends were calculated using galaxies from 1024 uniform-box simulations divided into the following stellar mass intervals: log10M*/M⊙ < 9.3, 9.3 ≤ log10M*/M⊙ < 9.6, 9.6 ≤ log10M*/M⊙ < 10.4, 10.4 ≤ log10M*/M⊙ < 11.2, and log10M*/M⊙ ≥ 11.2. Red and blue indicate trends for high and low parameter values, respectively, while the dashed gray line marks the threshold N*, 1/2 = 50. The thresholds and median values for the parameters are provided in Table 1. Error bars represent the scatter (16th–84th percentile range), while shaded areas indicate uncertainty estimates on the median values.

In the text

	Fig. 2. Same as Fig. 1, but for the astrophysical parameters ASN1, ASN2, and BHFF.
In the text

Fig. 3. Median scaling relations as a function of total stellar mass showing the dependence on astrophysical parameters and WDM mass for both MW zoom-in simulations (solid line) and calibrated uniform-box simulations (dashed line) for comparison. Trends are derived from 1024 MW zoom-in simulations divided into the following stellar mass intervals: log10M*/M⊙ < 7.5, 7.5 ≤ log10M*/M⊙ < 8.5, 8.5 ≤ log10M*/M⊙ < 10, and log10M*/M⊙ ≥ 10. Red and blue indicate trends for high and low parameter values, respectively. The dashed gray line indicates N*, 1/2 = 50.

In the text

	Fig. 4. Differences between the median trends corresponding to the high and low values of two parameters: ASN1 (green) and MWDM (brown). Solid lines represent MW zoom-in simulations; dashed lines represent uniform-box runs. Shaded regions show propagated uncertainties. The dashed gray line indicates N*, 1/2 = 50.
In the text

Fig. 5. Median number of galaxies as a function of median total stellar mass computed in the stellar mass bins log10M*/M⊙ < 6.5, 6.5 ≤ log10M*/M⊙ < 7.5, 7.5 ≤ log10M*/M⊙ < 8.5, and log10M*/M⊙ ≥ 8.5. The four panels, from top to bottom, correspond to variations in ASN1, ASN2, BHFF, and MWDM, respectively. In each case, the red and blue curves represent parameter values above and below the thresholds listed in Table 1, respectively. The plots are based on the 1024 MW zoom-in simulations from DREAMS.

In the text

Fig. 6. Recovery plot for the uniform-box simulations showing the bootstrap results with uncertainties for the estimated parameters: Ωm, σ8, MWDM, ASN1, ASN2, and BHFF against true values. These parameters were obtained by comparing the simulations with the simulation outputs mimicking the SPARC catalog. The dashed gray line indicates the ideal one-to-one relation between estimated and true parameter values.

In the text

Fig. 7. Recovery plot for the MW zoom-in simulations. The figure shows the bootstrap results with uncertainties for the estimated parameters: MWDM, ASN1, ASN2, and BHFF (from top to bottom). These parameters were obtained by comparing the simulations with the simulation outputs mimicking SPARC and the dwarf catalog. The dashed gray line indicates the ideal one-to-one relation between estimated and true parameter values.

In the text

Fig. 8. Best-fit simulations (points) compared to median observational trends (shaded areas) from SPARC (green) and the dwarf catalog (orange). From top to bottom: Stellar half-mass radius, total-to-stellar mass ratio within R*, 1/2, DM mass within R*, 1/2, and total versus stellar mass. Green points indicate the best-fit uniform-box versus SPARC, white points represent the MW zoom-in best matching SPARC, and orange points show the MW zoom-in best matching the dwarf catalog.

In the text

	Fig. A.1. Stellar-DM mass relation for the WDM uniform-box 1P DREAMS simulation with a 16 keV WDM particle mass (blue), its calibrated counterpart (yellow), and the TNG100-1 fiducial run (light gray), used as reference.
In the text

	Fig. A.2. Median scaling relations with 16th and 84th percentiles for the reference simulation TNG100-1 (light gray), the uncalibrated 1P run at 16 keV (blue), and the stellar-mass–calibrated 1P run (yellow).
In the text

	Fig. B.1. Stellar–DM mass relation for the WDM uniform-box 1P DREAMS simulation with a 16 keV WDM particle mass (blue), its stellar-mass–calibrated counterpart (yellow), and the corresponding reference samples (gray). The upper panel shows the calibration relative to the satellite galaxies from the MW zoom-in simulations, while the lower panel shows the calibration relative to TNG50-1.
In the text

	Fig. B.2. Stellar half-mass radius–DM mass relation for the WDM uniform-box 1P DREAMS simulation with a 16 keV WDM particle mass (blue), its calibrated counterpart (yellow), and the satellite galaxies from the MW zoom-in simulations used as reference (gray). The equivalent plot for the stellar-mass calibration is shown in the left panel of Fig. B.1.
In the text

	Fig. B.3. Median scaling relations with 16th and 84th percentiles for the MW zoom-in satellites (light gray), TNG50-1 (dark gray), the uncalibrated 1P run (blue), the 1P run calibrated on stellar mass with respect to the MW zoom-ins (yellow), the same run further calibrated on stellar half-mass radius (orange), and the 1P run calibrated on stellar mass with respect to TNG50-1 (red).
In the text

	Fig. C.1. Estimate and uncertainty of the half-light radius for galaxies in the DREAMS uniform-box simulations, computed for four different fixed total numbers of stellar mass particles (Npar = [40, 150, 600, 2000]). The dashed red line marks the true half-light radius.
In the text

	Fig. D.1. Median trends of the total-to-stellar half-mass ratio as a function of total stellar mass, for different cosmological and astrophysical parameters in the DREAMS (left, thick lines) and CAMELS (right, dashed lines) simulations. Colored curves correspond to subsets of simulations grouped by parameter value; the legend indicates the median value for each group.
In the text

Fig. E.1. Trends of stellar half-mass radius, total-to-stellar mass ratio within the stellar half-mass radius, DM mass within the stellar half-mass radius, total galaxy mass and total stellar mass for the host galaxies in the MW zoom-in simulations, as a function of astrophysical parameters and the inverse of the WDM mass parameter. The light blue points represent individual galaxies; red lines indicate the median trends and their scatter.

In the text

Fig. F.1. Recovery plot for the uniform-box simulations, showing the bootstrap results with uncertainties for the estimated parameters: Ωm, σ8, MWDM, ASN1, ASN2 and BHFF against their true values. These parameters are obtained by comparing the simulations with the corresponding simulation outputs mimicking the SPARC catalog. The dashed gray line indicates the ideal one-to-one relation between estimated and true parameter values. Blue points correspond to the reference analysis with N = 100, while grey points show the test performed with N = 10.

In the text

	Fig. F.2. Normalized residual distributions (Δp) for all parameters in the uniform-box suite. Each panel shows the empirical distribution (histogram) together with a Gaussian fit (red solid curve); the fitted mean (μ) and standard deviation (σ) are indicated in each panel. Histograms are normalized to unit area using Mathematica’s "PDF" option, allowing direct comparison with a standard normal.
In the text