TangoSIDM Project: Is the Stellar Mass Tully-Fisher relation consistent with SIDM?

Abstract

Self-interacting dark matter (SIDM) has the potential to significantly influence galaxy formation in comparison to the cold, collisionless dark matter paradigm (CDM), resulting in observable effects. This study aims to elucidate this influence and to demonstrate that the stellar mass Tully-Fisher relation imposes robust constraints on the parameter space of velocity-dependent SIDM models. We present a new set of cosmological hydrodynamical simulations that include the SIDM scheme from the TangoSIDM project and the SWIFT-EAGLE galaxy formation model. Two cosmological simulations suites were generated: one (Reference model) which yields good agreement with the observed $z=0$ galaxy stellar mass function, galaxy mass-size relation, and stellar-to-halo mass relation; and another (WeakStellarFB model) in which the stellar feedback is less efficient, particularly for Milky Way-like systems. Both galaxy formation models were simulated under four dark matter cosmologies: CDM, SIDM with two different velocity-dependent cross sections, and SIDM with a constant cross section. While SIDM does not modify global galaxy properties such as stellar masses and star formation rates, it does make the galaxies more extended. In Milky Way-like galaxies, where baryons dominate the central gravitational potential, SIDM thermalises, causing dark matter to accumulate in the central regions. This accumulation results in density profiles that are steeper than those produced in CDM from adiabatic contraction. The enhanced dark matter density in the central regions of galaxies causes a deviation in the slope of the Tully-Fisher relation, which significantly diverges from the observational data. In contrast, the Tully-Fisher relation derived from CDM models aligns well with observations.

1 Introduction

The self-interacting dark matter paradigm (SIDM) postulates that dark matter particles engage in gravitational interactions with ordinary particles while exhibiting non-gravitational interactions among themselves. Arising as a natural prediction of dark sector models beyond the Standard Model (e.g. Spergel & Steinhardt 2000; Tulin & Yu 2018), SIDM is expected to manifest detectable astrophysical signatures (e.g. Adhikari et al. 2022). Moreover, it offers a potential explanation for the most challenging discrepancy between $\Lambda$ cold dark matter ($\Lambda$CDM) numerical simulations and observations: the diverse distribution of dark matter within dwarf galaxies (see e.g. Oman et al. 2015; Santos-Santos et al. 2020; Hayashi et al. 2021; Sales et al. 2022; Borukhovetskaya et al. 2022).

Within the SIDM framework, interactions among dark matter particles dynamically alter the internal structure of dark matter halos. This modification involves the transfer of heat from the outer parts to the inner halo, resulting in an increase in the velocity dispersion, and a reduction of dark matter densities in the central regions (e.g. Davé et al. 2001; Colín et al. 2002; Vogelsberger et al. 2012; Rocha et al. 2013; Dooley et al. 2016; Vogelsberger et al. 2016). The crucial parameter governing the rate of dark matter particle interactions is the cross section per unit mass, denoted as $\sigma/m_{\chi}$ (e.g. Robertson et al. 2017; Kahlhoefer et al. 2019; Kummer et al. 2019; Vogelsberger et al. 2019; Banerjee et al. 2020; Shen et al. 2021). Measurements derived from the shape and collision of nearby galaxy clusters constrain this parameter to be ${<}1~{}\rm{cm}^{2}\rm{g}^{-1}$ (e.g. Randall et al. 2008; Dawson et al. 2013; Massey et al. 2015; Harvey et al. 2015; Wittman et al. 2018; Harvey et al. 2019; Sagunski et al. 2021; Andrade et al. 2022).

While various studies have explored the impact of SIDM under a small and constant cross section, prevailing particle physics models advocate for a velocity-dependent framework, where $\sigma/m_{\chi}$ allows dark matter to behave as a collisional fluid on small scales while remaining essentially collisionless over large scales (e.g. Pospelov et al. 2008; Arkani-Hamed et al. 2009; Buckley & Fox 2010; Feng et al. 2010; Boddy et al. 2014; Tulin & Yu 2018). Under this velocity-dependent scheme, $\sigma/m_{\chi}$ can be ${<}1~{}\rm{cm}^{2}\rm{g}^{-1}$ for high dark matter velocities at large scales, aligning with the constrains of cluster-size haloes, and exceed ${>}100~{}\rm{cm}^{2}\rm{g}^{-1}$ for low dark matter velocities in order to explain the diverse dark matter distribution within dwarf galaxies (e.g. Correa 2021; Gilman et al. 2021; Correa et al. 2022; Yang et al. 2023; Silverman et al. 2023; Nadler et al. 2023; Shah & Adhikari 2023; Gilman et al. 2023). Although SIDM has been robustly constrained on galaxy cluster scales, uncertainties persist in the lower-mass galaxy regime due to the difficulty in isolating the impact of baryonic physics from dark matter interactions.

Recent studies exploring the co-evolution of baryons and SIDM in isolated systems indicate that non-bursty stellar feedback may not significantly alter SIDM density profiles in dwarf galaxies (e.g. Vogelsberger et al. 2014; Robles et al. 2017; Sameie et al. 2021). Conversely, hydrodynamical simulations incorporating SIDM and a bursty stellar feedback model reveal distinctions in velocity dispersion profiles between SIDM and CDM haloes (Burger et al. 2022), suggesting the need for more detailed investigations into the interplay between SIDM and various feedback models. In more massive systems, the intricate interplay between SIDM and baryons is even more challenging. Studies that modelled the evolution of Milky Way-like systems and galaxy clusters (e.g. Robertson et al. 2019; Despali et al. 2019; Sameie et al. 2021; Rose et al. 2022) found that baryon contraction results in the formation of denser and cuspier central density profiles under SIDM compared to CDM. Analytical studies focusing on the gravitational contribution of a baryonic disc and bulge reached similar conclusions (Robles et al. 2019; Silverman et al. 2023; Jiang et al. 2023). However, uncertainties persist regarding how the increased cuspiness of SIDM haloes depends on the specific SIDM model parameters or the strength of galaxy feedback models.

This paper seeks to address this knowledge gap by introducing a new set of cosmological hydrodynamical simulations. These simulations integrate the SIDM model derived from the TangoSIDM project, with the baryonic physics from the SWIFT-EAGLE galaxy formation model. The goals of the TangoSIDM project are to derive robust constraints on the dark matter cross section from observations of dwarf and Milky Way-type galaxies. In this study, we take a pivotal first step by demonstrating how the stellar mass Tully-Fisher relation, a well-established galaxy scaling relation, can be leveraged to derive robust constraints on the parameter space of velocity-dependent SIDM models. The structure of this paper is organized as follows. Section 2 describes the SIDM and baryonic subgrid models employed in our simulations. In Section 3, we show how SIDM influences key galaxy properties, including stellar masses, sizes, and star formation rates. Section 4 compares the dark matter density profiles of haloes between CDM and various SIDM models. Section 5 undertakes an in-depth analysis of the stellar mass Tully-Fisher relation, and demonstrates it rules out the velocity-dependent SIDM models studied in this work. Section 6 discusses the SIDM parameter space, and Section 7 summarizes the paper’s findings.

2 Simulation setup

TangoSIDM is a simulation project dedicated to modelling cosmological simulations that capture the intricacies of structure formation within a $\Lambda$SIDM universe. This work introduces the first realization of hydrodynamical cosmological volumes, each spanning 25 Mpc on a side, as integral compontents of the TangoSIDM project. To produce these simulations, the SWIFT code (Schaller et al. 2023) was employed. SWIFT includes advanced hydrodynamics and gravity schemes. The gravity solver employs the Fast Multiple Method (Greengard & Rokhlin 1987) with an adaptive opening angle, while for hydrodynamics the SPHENIX SPH scheme (Borrow et al. 2022), specifically designed for galaxy formation sub-grid models, was utilized.

The simulations follow the evolution of 376${}^{3}$ dark matter particles and 376${}^{3}$ gas particles to redshift $z=0$. The softening is set to 2.66 comoving kpc at early times, but is frozen a physical value of 700 pc at $z=2.8$. The dark matter particle mass is $9.70\times 10^{6}~{}\rm{M}_{\odot}$ and the gas initial particle mass is $1.81\times 10^{6}~{}\rm{M}_{\odot}$. The starting redshift of the simulations is $z=127$. The initial conditions were calculated using second-order Lagrangian perturbation theory with the method of Jenkins (2010, 2013). The adopted cosmological parameters are $\Omega_{\rm{m}}=0.307$, $\Omega_{\Lambda}=0.693$, $h=0.6777$, $\sigma_{8}=0.8288$ and $n_{\rm{s}}=0.9611$ (Planck Collaboration et al. 2014).

2.1 TangoSIDM model

Figure 1: Momentum transfer cross section as a function of the relative scattering velocity among dark matter particles for the SIDM models featured in this work (Table 1). The figure shows two velocity-dependent models, namely SigmaVel60 (light blue line) and SigmaVel30 (dark blue line), alongside SigmaConstant10 (orange line), which uses a constant cross section, $\sigma_{T}/m_{\chi}=10$ cm${}^{2}$g${}^{-1}$. The top x-axis indicates the typical halo mass that hosts orbits of the velocities indicated on the bottom x-axis.

The TangoSIDM project, encompassing its models and SIDM implementation, was presented in Correa et al. (2022). In this section we briefly summarise the key elements of the SIDM model, with further details available in the aforementioned reference.

Four dark matter models were generated for this study: the cold collisionless dark matter model (hereafter CDM); a SIDM model with a constant scattering cross section of 10 cm${}^{2}\rm{g}^{-1}$ (hereafter SigmaConstant10); and two SIDM models featuring velocity-dependent cross sections (see Fig. 1). Although the SigmaConstant10 model has been ruled out by observations of galaxy clusters (e.g. Harvey et al. 2015, 2019), it serves as a control model for comparative analysis. Among the velocity-dependent models, one has a cross section that is below 1 cm${}^{2}\rm{g}^{-1}$ at high velocities ($v{>}150$ km s${}^{-1}$) and increases with decreasing velocity, reaching 60 cm${}^{2}\rm{g}^{-1}$ at 10 km s${}^{-1}$ (hereafter SigmaVel60 model). The other velocity-dependent model has a cross section smaller than 8 cm${}^{2}\rm{g}^{-1}$ at velocities surpassing $200$ km s${}^{-1}$ (dropping below 1 cm${}^{2}\rm{g}^{-1}$ at ${\approx}1000$ km s${}^{-1}$) and increases with decreasing velocity, reaching 30 cm${}^{2}\rm{g}^{-1}$ at 10 km s${}^{-1}$ (hereafter SigmaVel30 model).

The SigmaVel60 and SigmaVel30 models represent two extreme scenarios for the rate of dark matter interactions in Milky Way-mass systems. Despite both models adhering to the SIDM constraints derived from cluster-size haloes, there are important differences. In SigmaVel60, interactions reach 1-2 cm${}^{2}\rm{g}^{-1}$ around 100 km s${}^{-1}$, therefore this model produces a low rate of interactions in the center of Milky Way-like haloes. In contrast, SigmaVel30 exhibits a cross section of 10-20 cm${}^{2}\rm{g}^{-1}$ at 100 km s${}^{-1}$, imposing a stronger rate of interaction.

The velocity-dependent cross sections are modelled under the assumption that dark matter particle interactions are mediated by a Yukawa potential dependent on three parameters: the dark matter mass $m_{\chi}$; the mediator mass $m_{\phi}$; and the coupling strength $\alpha_{\chi}$. While there is no analytical form for the differential scattering cross-section due to a Yukawa potential, the Born-approximation (Ibe & Yu 2010)$-$applicable when treating the scattering potential as a small perturbation$-$yields the differential cross-section of the dark matter-dark matter interactions

$$\frac{{\rm{d}}\sigma}{{\rm{d}}\Omega}=\frac{\alpha_{\chi}^{2}}{m_{\chi}^{2}(m_{\phi}^{2}/m_{\chi}^{2}+v^{2}\sin^{2}(\theta/2))^{2}}.$$ (1)

While in the model with a constant cross section the dark matter scattering is isotropic, in the velocity-dependent cross section models the scattering is anisotropic. For anisotropic scattering the momentum transfer cross section, defined as

$$\sigma_{T}/m_{\chi}=2\int(1-|\cos\theta|)\frac{{\rm{d}}\sigma}{{\rm{d}}\Omega}{\rm{d}}\Omega,$$ (2)

is useful to consider, because it is weighted by the scattering angle and therefore it does not overestimate the scattering with $\theta>\pi/2$ (Kahlhoefer et al. 2015). Table 1 shows the SIDM model parameters adopted in this work and Fig. 1 shows the momentum transfer cross sections. The figure shows the velocity-dependent models (light blue and dark blue lines) and the constant cross section model (orange line). While the bottom x-axis shows the relative velocity between dark matter particles, the top x-axis indicates the typical halo mass that hosts circular orbits of such velocities.

Table 1: SIDM models analysed in this work. Form left to right: Model name, SIDM parameters for each model (dark matter mass, $m_{\chi}$, mediator mass, $m_{\phi}$, and coupling strength, $\alpha$) and type of dark matter interaction.

	SIDM parameters			DM interaction
Model	$m_{\chi}$	$m_{\phi}$	$\alpha$
Name	[GeV]	[MeV]
CDM	-	-	-	No interaction
SigmaConstant10	-	-	-	Isotropic
SigmaVel30	2.227	0.778	$4.317\times 10^{-5}$	Anisotropic
SigmaVel60	3.855	0.356	$1.027\times 10^{-5}$	Anisotropic

2.2 SWIFT-EAGLE model

The SWIFT-EAGLE model, an open-source galaxy formation model implemented in SWIFT, is derived from the original EAGLE model (Schaye et al. 2015; Crain et al. 2015). While it has common modules to those of EAGLE, SWIFT-EAGLE includes new developments and improvements. A detailed model description can be found in Bahé et al. (2022) and Borrow et al. (2023). Below we provide a summary.

SWIFT-EAGLE incorporates the element-by-element sub-grid radiative gas cooling and photoheating prescription from Ploeckinger & Schaye (2020), which accounts for the inter-stellar radiation field and self-shielding of dense gas, as well as the UV/X-ray background from galaxies and quasars according to Faucher-Giguère (2020). Star formation is implemented stochastically, following the Schaye & Dalla Vecchia (2008) pressure law, as in the original EAGLE model. A polytropic equation of state, $P\propto\rho^{4/3}$, sets a minimum limit on the gas pressure. The star formation rate per unit mass is calculated from the gas pressure, employing an analytical formula designed to reproduce the observed Kennicutt–Schmidt law (Kennicutt 1998) in disc galaxies. A gas particle is star-forming if its subgrid temperature $T<1000$ K, or if its density (expressed in units of hydrogen particles per cubic cm, $n_{\rm{H}}$) is $n_{\rm{H}}>10~{}\rm{cm}^{-3}$ and temperature $T<10^{4.5}$ K.

The stellar initial mass function assumes the form of Chabrier (2003) within the range 0.1-100 M${}_{\odot}$, with each particle representing a simple age stellar population. Stellar feedback is implemented stochastically, following the prescription of Dalla Vecchia & Schaye (2012), where stars with masses between 8 $\rm{M}_{\odot}$ and 100 $\rm{M}_{\odot}$ explode as core-collapse supernovae. The resulting energy is transferred as heat to the surrounding gas, following Chaikin et al. (2022).

The energy injected into the gas corresponds to $10^{51}$ erg per supernova times a dimensionless coupling efficiency factor, $f_{\rm{E}}$, that follows the same scaling function as in EAGLE,

$$f_{\rm{E}}=f_{\rm{E,max}}-\frac{f_{\rm{E,max}}-f_{\rm{E,min}}}{1+\rm{exp}\left(\frac{-\log_{10}Z/Z_{0}}{\sigma_{Z}}\right)\rm{exp}\left(\frac{\log_{10}n_{H}/n_{H,0}}{\sigma_{n}}\right)}.$$ (3)

As can be seen, $f_{\rm{E}}$ depends on a number of free parameters: $f_{\rm{E,min}}$ and $f_{\rm{E,max}}$, which set the minimal and maximal feedback energies, $n_{\rm{H},0}$ and $Z_{0}$ defined as the density and metallicity pivot point around which the feedback energy fraction plane rotates, and $\sigma_{Z}$ and $\sigma_{n}$, the width of the feedback energy fraction sigmoids in the metallicity and density dimensions.

In addition to the energy released through star formation, star particles also release metals into the inter-stellar medium (ISM) through four evolutionary channels: AGB stars, winds from massive stars, core-collapse supernvae and Type Ia supernovae. This process follows the methodology discussed in Wiersma et al. (2009) and Schaye et al. (2015). The abundances of 9 elements (H, He, C, N, O, Ne, Mg, Si, Fe) are tracked.

The formation and growth of supermassive black holes are modelled following Bahé et al. (2022). Initially seeded within friends-of-friends dark matter groups of mass $10^{10}$ M${}_{\odot}$, black holes accretion rates follow the Eddington-limited Bondi accretion rate. The feedback mechanism from active galactic nucleus (AGN) activity is implemented following Booth & Schaye (2009). The energy depends on the accreted mass, $\Delta m$, onto the black hole as, $\Delta E=\epsilon_{\rm{r}}\epsilon_{\rm{f}}\Delta mc^{2}$, where $\epsilon_{\rm{r}}=0.1$ is the default value. This energy is stored in a reservoir carried by each black hole particle until it can be utilized to heat the nearest gas particle, inducing a temperature increase of $\Delta T_{\rm{AGN}}$. The coupling efficiency, $\epsilon_{\rm{f}}$, and the heating temperature of AGN feedback are free parameters.

2.3 Reference & WeakStellarFB SWIFT-EAGLE models

This work investigates the evolution of galaxies for two distinct SWIFT-EAGLE models. In the first, referred to as the Reference model, the free parameters described in the previous subsection were calibrated in a (25 Mpc)${}^{3}$ volume to reproduce the galaxy stellar mass function and galaxy mass-size relation. The second, named the WeakStellarFB model, adopts parameters that produce Milky Way-mass galaxies with very weak stellar feedback. Table 2 provides a comprehensive listing of the subgrid parameter values for both models.

The parameters for the Reference model were derived within the CDM framework using emulators that employed the Gaussian Process Regression-based python module SWIFTEmulator (Kugel & Borrow 2022). Further details on the calibration and emulation technique can be found in Borrow et al. (2023). Note that the SIDM simulations with the SWIFT-EAGLE Reference model adopt the parameters listed in Table 2, no re-calibration was performed to account for the SIDM effects.

The original parameters from the EAGLE simulations were calibrated to reproduce the $z=0.1$ galaxy stellar mass function, the relation between galaxies stellar mass and galaxies’ central black hole masses, as well as disc galaxy sizes (Crain et al. 2015). While the SWIFT-EAGLE model was inspired by EAGLE, significant differences exist, such as the gravity and hydrodynamics solver, cooling rates, supernovae and AGN feedback energy deposition into the ISM. Because of these differences, applying the original EAGLE parameter values in the SWIFT-EAGLE model yields different results. Relative to the Reference model, the WeakStellarFB model exhibits a weaker stellar feedback at the specific mass scale of $10^{12}~{}\rm{M}_{\odot}$ haloes, attributed to the lower value of $f_{\rm{E,max}}$ and higher $n_{\rm{H,0}}$. This combination results in a lower coupling efficiency factor $f_{\rm{E}}$ at fixed hydrogen number density, justifying its nomenclature “WeakStellarFB”.

In Section 3 and Appendix A, we show that both the Reference and WeakStellarFB models yield stellar mass functions, specific star formation rates, and stellar-to-halo mass relations that closely align with observational data. However, the stellar feedback in the WeakStellarFB model is less efficient in Milky Way-mass systems, making them more compact by redshift zero. The primary objective of exploring SIDM under these two galaxy models is to understand the impact of dark matter collisions in the central regions of galaxies. We aim to discern how SIDM coevolves with the dynamical heating from supernova explosions and evaluate whether our conclusions regarding the impact of SIDM on galaxies remain robust in the face of variations in feedback models.

Table 2: Subgrid parameter values of the SWIFT-EAGLE galaxy formation model that regulate stellar and AGN feedback. The left column identifies each parameter, with detailed descriptions provided in the text. The middle and right columns list the parameter values adopted in the Reference and WeakStellarFB models, respectively.

Parameters	Reference	WeakStellarFB
$f_{\rm{E,min}}$	0.388	0.5
$f_{\rm{E,max}}$	7.37	5.0
$n_{\rm{H,0}}$ [cm${}^{-3}$]	0.412	1.46
$\sigma_{\rm{Z}}$	0.311	0.275
$Z_{0}$	0.00134	0.00134
$\sigma_{\rm{n}}$	0.428	1.77
$\epsilon_{\rm{f}}$	0.035	0.1
$\Delta T_{\rm{AGN}}$ [K]	10${}^{8.62}$	10${}^{8.5}$

Figure 2: Galaxy scaling relations at redshift $z=0$ for the Reference (top panels) and WeakStellarFB model (bottom panels). The columns show the stellar-to-halo mass ratio ($M_{*}/M_{200\rm{c}}$) as a function of halo mass (left), the projected stellar half-mass radius as a function of stellar mass (middle) for all galaxies, and the specific star formation rate (sSFR, SFR/$M_{*}$) as a function of stellar mass for actively star forming-galaxies (right). In all panels the curves correspond to the median for the Reference and WeakStellarFB models produced in CDM (blue lines), SigmaVel60 (purple lines) and SigmaConstant10 (orange lines) frameworks. SigmaVel30, though not shown, follows a similar trend as SigmaVel60. The shaded regions mark the 16th-84th percentiles of the relations. These models are contrasted with various observational datasets and the EAGLE simulations (black lines). The left panels show the stellar-to-halo mass relation from Behroozi et al. (2019). In the middle panels, galaxy sizes are compared with citetLange15 (green circles) dataset. The right panels compare the sSFR with those reported by Bauer et al. (2013) and Chang et al. (2015). SIDM appears to have minimal impact on galaxy masses and star formation rates. However, it significantly alters galaxy sizes, leading to increases by up to a factor of 2 for SigmaConstant10.

2.4 Halo catalogue and definitions

Halo catalogues were generated using the VELOCIraptor halo finder (Elahi et al. 2011; Elahi et al. 2019; Cañas et al. 2019). VELOCIraptor uses a 3D-friends of friends (FOF) algorithm to identify field haloes, and subsequently applies a 6D-FOF algorithm to separate virialised structures and identify sub-haloes of the parent haloes (Elahi et al. 2019). Throughout this work, virial halo masses ($M_{200\rm{c}}$) are defined as all matter within the virial radius $R_{200\rm{c}}$, for which the mean internal density is 200 times the critical density, $\rho_{\rm{crit}}$, which is $127.5\rm{M}_{\odot}\rm{kpc}^{-3}$ at $z=0$. In each FOF halo, the ‘central’ subhalo is the one that is most likely the core in phase space, which is nearly always the most massive. The remaining subhaloes within the FOF halo are its satellites. The resolution of the simulations is sufficient to resolve (sub-)haloes down to ${\sim}10^{10}~{}\rm{M}_{\odot}$ with $10^{3}$ particles within $R_{200}$. Galaxy stellar masses, sizes and star formation rates are always defined within an aperture of 50 kpc.

3 Galaxy properties

In this section we analyse key galaxy properties from the Reference and WeakStellarFB models: the $z=0$ stellar-to-halo mass relation, projected galaxy sizes, and star formation rates, and we compare them against observational data. It is important to point out that during the calibration of the subgrid parameters for feedback under CDM, the $z=0$ galaxy stellar mass function and the stellar mass-size relation were considered, and as a result, the simulations do not provide predictions for these. We remind the reader that the subgrid parameters from the Reference model were only calibrated under the CDM framework and not under SIDM. The SIDM simulations use the same subgrid parameter values as CDM for both the Reference and WeakStellarFB models. The $z{=}0$ galaxy stellar mass function is presented in Appendix A.

Fig. 2 illustrates three galaxy scaling relations from the Reference (the top panels) and WeakStellarFB models (bottom panels). In the left panels, the ratio between the galaxy stellar mass and halo mass ($M_{*}/M_{200\rm{c}}$) is plotted as a function of the host halo mass. Coloured curves represent the median relations for central galaxies, with shaded regions indicating the 16-84th percentiles. A comparison is made with the stellar-to-halo mass relation from the EAGLE simulation and from UNIVERSEMACHINE (Behroozi et al. 2019). Notably, the WeakStellarFB model aligns best with the original EAGLE data (McAlpine et al. 2016). At fixed halo mass, galaxies from the WeakStellarFB model are more massive than those from the Reference model (consistent with a comparison of the stellar mass functions). The dark matter framework does not significantly alter the stellar-to-halo mass relation. For clarity, the SigmaVel30 model is not shown, as it follows a trend similar to SigmaVel60.

Moving to the middle panels of Fig. 2, the stellar half-mass radius is shown as a function of stellar mass. The half-mass radius is defined as the radius that encloses 50 per cent of the stellar mass, and is computed from all bound star particles within a projected 2D circular aperture of 50 kpc radius. The simulations are compared against the GAMA survey (Lange et al. 2015), and the EAGLE simulation (McAlpine et al. 2016). An interesting feature emerges in the bottom middle panel, revealing a U-shape trend in the galaxy size-mass relation. Galaxies within the mass range of $10^{9}$ to $10^{11}~{}\rm{M}_{\odot}$ become too compact due to excessive radiative losses at high gas densities. To counteract this issue, the EAGLE model introduced a dependence of the stellar feedback energy on the gas density (eq. 3), so that higher density gas receives a larger amount of energy from stellar explosions (Crain et al. 2015). The WeakStellarFB model incorporates the density-dependent stellar feedback energy, but its parameter values are such that the feedback strength remains inadequate. The coupling efficiency factor applied to the supernova energy that is injected into that gas is smaller than in the Reference model. Therefore, while stellar and AGN feedback in the WeakStellarFB model can prevent the formation of excessively massive galaxies, it does not guarantee the formation of extended galaxies with realistic sizes. A more careful approach, or tuning of the energy parameters, is required for feedback to effectively eject low-angular momentum gas, increase the median angular momentum of the ISM gas that remains to form stars, and form more extended galaxies (e.g. Brook et al. 2012).

For stellar masses ${\approx}10^{9}~{}\rm{M}_{\odot}$ the WeakStellarFB model predicts galaxies with sizes that agree with EAGLE, and do not seem to suffer from overcooling and compactness. However, these sizes appear large when compared to the dataset of Lange et al. (2015). The Reference model, calibrated to match the size-mass relation from Lange et al. (2015), yields galaxies that are still overly extended, partly due to the sampling noise in gravitational interactions between stars and dark matter, that leads to spurious size growth (Ludlow et al. 2019a, 2023).

The middle panels also show the evident impact of SIDM on galaxy sizes. Dark matter particle interactions heat the inner halo, leading to core formation in the central regions and dynamically heating the surrounding gas and stars, promoting the formation of more extended galaxies. However, this is insufficient to counteract the overcooling and compactness observed in the WeakStellarFB model for galaxies more massive than $10^{10}~{}\rm{M}_{\odot}$. The top middle panel shows that the SigmaVel60 model, characterized by a large cross section for galaxies less massive than $10^{9}M_{\odot}$, produces sizes that are close to those for the SigmaConstant10 model. For these masses, as the cross section decreases, the galaxy sizes from SigmaVel60 decrease relative to those for SigmaConstant10, and become similar to those of CDM.

Ludlow et al. (2023) found than in CDM hydrodynamical simulations like EAGLE, which share the same numerical resolution as the TangoSIDM simulations, the galaxies’ half-mass radius remains robust against spurious collisional heating only for halo masses $M_{200\rm{c}}{\gtrsim}10^{11.7}~{}\rm{M}_{\odot}$. This suggests that our galaxies’ sizes are free from spurious heating if the galaxies are more massive than $M_{*}{\gtrsim}10^{10}~{}\rm{M}_{\odot}$. We note, however, that resolution effects may have a stronger impact on CDM simulations than on SIDM simulations, in which case the effect of SIDM on sizes relative to CDM may be underestimated.

The right panels of Fig. 2 display the median specific star formation rates (sSFR) for actively star forming-galaxies, with galaxies classified as star-forming if their sSFR ${>}10^{-11}$ yr${}^{-1}$. The panels reveal that the $z=0$ sSFR from the Reference model are in agreement with the sSFR from the EAGLE simulations and the dataset from Chang et al. (2015), and are within a factor of 5 from the Bauer et al. (2013) data. The WeakStellarFB model has sSFR lower than Reference. Interestingly, there are no differences in the median sSFR trends between simulations with CDM vs. SIDM.

Figure 3: Dark matter density profiles, $\rho_{\rm{DM}}$, of $10^{11}$ M${}_{\odot}$, $10^{11.5}$ M${}_{\odot}$ and $10^{12}$ M${}_{\odot}$ haloes from the Reference (purple solid lines) and WeakStellarFB (blue dot-dashed lines) models under CDM (left panels), SigmaConstant10 (second panels from the left), SigmaVel30 (third panels from the left) and SigmaVel60 (right panels). The panels compare $\rho_{\rm{DM}}$ between hydrodynamical (CDM and SIDM) simulations and dark matter-only simulations (orange dashed lines). The coloured lines highlight the median values and the shaded regions the 16-84th percentiles. Additionally, the black solid line corresponds to the NFW profile (estimated using the concentration-mass relation from Correa et al. 2015), and the black dotted and dashed-dotted lines indicate the convergence radii (see text for definition). The differences in the profiles between haloes of the same mass highlights the impact of baryonic effects and dark matter particle interactions on the central haloes densities.

Figure 4: Stacked dark matter density profiles, $\rho_{\rm{DM}}$, of the 32 most massive haloes in the box (with masses larger than $10^{12}$ M${}_{\odot}$) at $z=0$ from the Reference model under CDM (left panel) and SigmaVel30 (right panel). The panels show the median density evolution between redshifts 0 and 2. The coloured lines highlight the median values and the black solid line shows the NFW profile of the haloes at redshift zero (estimated using the concentration-mass relation from Correa et al. 2015). The black dotted and dashed-dotted lines indicate the convergence radii (see text for definition). While there is no large difference in the median density profiles of haloes over time in the CDM, the SigmaVel30 model shows that the central density increases.

4 Dark Matter Density profile

In the following analysis, we compare our findings with prior studies on SIDM. Specifically, we examine the dark matter density profiles of central haloes with masses in the range of $10^{10.9}-10^{11.1}$ M${}_{\odot}$, $10^{11.4}-10^{11.6}$ M${}_{\odot}$ and $10^{11.9}-10^{12.1}$ M${}_{\odot}$ from both the Reference and WeakStellarFB models under CDM, SigmaConstant10, SigmaVel30 and SigmaVel60.

The panels in Fig. 3 compare the dark matter density profiles, $\rho_{\rm{DM}}$, between hydrodynamical simulations (CDM and SIDM) of the Reference model (purple solid lines) and the WeakStellarFB model (blue dot-dashed lines). Additionally, dark matter-only simulations are presented as orange dashed lines. To facilitate the comparison, the NFW density profile (black solid lines) is included, estimated using the concentration-mass relation from Correa et al. (2015). We also include convergence radii defined as the minimum radius where the mean density converges at the 20 and 10 per cent level, $r_{\rm{c},20}$ (dash-dotted lines) and $r_{\rm{c},10}$ (dotted lines) respectively, relative to a simulation of higher resolution. The convergence criterion $r_{\rm{c},10}$, presented by Ludlow et al. (2019b), is defined as $r_{\rm{c},10}=0.055l(z)$, where $l(z)$ is the (comoving) mean inter-particle separation. At $z=0$ this separation is $l=L_{\rm{b}}/N_{\rm{p}}^{1/3}=52.7$ kpc, given $L_{\rm{b}}=25$ cMpc and $N_{\rm{p}}=2\times 376^{3}$ particles. In addition to Ludlow et al. (2019b) criterion, we include a relaxed convergence criterion given by $r_{\rm{c},20}=0.034l(z)$. This is motivated by the findings of Schaller et al. (2015), who showed that the differences in the mean density profiles from the EAGLE hydrodynamical and DM-only simulations are significantly larger than 10%. The value of 0.034 is obtained from eq. (18) of Ludlow et al. (2019b) after decreasing $\kappa_{\rm{P}03}\equiv\frac{t_{\rm{relax}}}{t_{200}}$ by a factor of 2 (see Power et al. 2003).

The bottom panels of Fig. 3 show that, in $10^{11}$ M${}_{\odot}$ haloes, baryons do not affect $\rho_{\rm{DM}}(r)$ beyond $r_{\rm{c},20}$. Under both CDM and SIDM, the hydrodynamical and DM-only simulations yield consistent $\rho_{\rm{DM}}$. In the CDM models, $\rho_{\rm{DM}}(r)$ agrees with the NFW prediction for $r>r_{\rm{c},20}$, while in SIDM models, dark matter particle interactions create the expected constant-density isothermal cores (see also e.g., Colín et al. 2002; Vogelsberger et al. 2012; Peter et al. 2013; Rocha et al. 2013, and Correa et al. 2022 for DM-only TangoSIDM density profiles). This cored $\rho_{\rm{DM}}$ corresponds to the median profile of the central $10^{11}$ M${}_{\odot}$ halo population. However, note that since the velocity-dependent SIDM models under consideration exhibit large cross sections at the $10^{11}~{}\rm{M}_{\odot}$ mass-scale, some SIDM haloes may potentially undergo core-collapse and form a cuspy central density profile.

The bottom panels of Fig. 3 reveal that galaxies with stellar masses as high as $10^{9}$ M${}_{\odot}$ do not produce sufficiently strong feedback to affect the underlying dark matter distribution. In agreement with our results, Robles et al. (2017) modelled dwarf galaxies within $10^{10}~{}\rm{M}_{\odot}$ haloes under both CDM and SIDM using the zoom-in FIRE cosmological model. They concluded that, for these low-mass systems, the final density profile of SIDM haloes was not strongly influenced by the stellar mass of the galaxy, exhibiting cored density profiles regardless of hosting galaxies with stellar masses ranging from $10^{5}$ to $10^{7}~{}\rm{M}_{\odot}$. Furthermore, Burger et al. (2022) showed that both CDM and SIDM can yield haloes with cored density profiles. The difference lies in the fact that, under SIDM, galaxies can be embedded in haloes with cored central dark matter profiles, irrespective of whether they have a smooth star formation history and non-bursty supernova feedback. In contrast, under CDM, galaxies would require a bursty star formation rate to generate strong supernova feedback that leads the impulsive cusp-core transformation.

Back to our results, the middle panels of Fig. 3 show that, in $10^{11.5}$ M${}_{\odot}$ haloes, baryons impact on the dark matter distribution from the SIDM models. Under CDM, $\rho_{\rm{DM}}(r)$ from the Reference hydrodynamical and DM-only simulations agree, but under SIDM, they diverge. The SIDM DM-only simulations produce lower-density and larger-core profiles compared to the SIDM hydrodynamical simulations, which more closely follow the NFW prediction. The WeakStellarFB model generates cuspier density profiles than the Reference model, both under CDM and SIDM. This result suggests that the increased baryonic concentration in the WeakStellarFB model, relative to the Reference model, enhances the central concentration of the dark matter distribution.

The influence of baryons becomes more pronounced in $10^{12}$ M${}_{\odot}$ haloes, as shown in the top panels of Fig. 3. In all dark matter models (CDM and SIDM), the density profiles between hydrodynamical and DM-only simulations no longer agree. Hydrodynamical-CDM models produce a cuspier $\rho_{\rm{DM}}(r)$ than the NFW profile (in line with predictions from adiabatic contraction models (e.g. Gnedin 2006). Similarly, hydrodynamical-SIDM models produce a very cuspy $\rho_{\rm{DM}}(r)$, in contrast to the cored $\rho_{\rm{DM}}$ profiles produced in the DM-only SIDM models. Consistent with our results, previous works by Elbert et al. (2018), Sameie et al. (2021) and Rose et al. (2022) showed that, under SIDM, dark matter density profiles can be either cuspy or even cuspier than their CDM counterparts, depending on the baryonic concentration. Sameie et al. (2021) analysed the density profiles of $10^{12}~{}\rm{M}_{\odot}$ haloes, modelled in high-resolution zoom-in simulations of SIDM within the FIRE galaxy formation scheme (Hopkins et al. 2018). Their study showed that SIDM haloes can reach higher and steeper central densities than their CDM counterparts. In a similar approach, Rose et al. (2022) presented zoom-in SIDM simulations of Milky Way-like galaxies with the IllustrisTNG (Pillepich et al. 2018) galaxy formation model. They concluded that baryon contraction begins to have an impact on the density profiles of haloes when their embedded galaxies reach stellar masses of $10^{8}~{}\rm{M}_{\odot}$. For higher-mass systems such as groups and clusters, the work of Robertson et al. (2021) concluded that the haloes profile strongly depends on the final baryonic distributions. They showed this from the analysis of dark matter halo densities modelled with SIDM and the baryonic physics model of EAGLE (Schaye et al. 2015) in a zoom-in simulation. Similarly, Despali et al. (2019), employing zoom-in SIDM simulations of galaxies with the IllustrisTNG model, showed that smaller-size galaxies were embedded in cuspy SIDM haloes, while more extended galaxies resided in cored-profile haloes.

The gravitational influence of baryons not only increases central dark matter densities in SIDM models, but also diversify the haloes’ dark matter distribution, which can be seen from the increased scatter around the median density profiles (shaded region in Fig. 3). This diversity could be attributed to variations in the assembly history of galaxies, influencing whether baryons dominate the central gravitational potential sooner or later.

In Fig. 4, we investigate how the haloes assembly history shapes the evolution of the DM density profile. We select the 32 most massive haloes (with masses larger than $10^{12}$ M${}_{\odot}$) at $z=0$ in the cosmological box from the Reference model under CDM (left panel) and SigmaVel30 (right panel). The panels show the evolution of the median $\rho_{\rm{DM}}$ between redshifts 0 and 2 (coloured lines). The black solid line shows the NFW profile of the haloes at redshift zero (estimated using the concentration-mass relation from Correa et al. 2015). The left panel shows that except for the inner few kpc, there is minimal evolution of $\rho_{\rm{DM}}(r)$ under CDM. In this case, haloes formed a cuspy profile by redshift two, the subsequent impact of the central galaxies, through ejection of energy via supernova- and AGN-driven winds, leads to the formation of small cores in the center.

Under SIDM, the haloes’ density evolves. At redshift two, the central dark matter density of SIDM haloes is lower than for their CDM counterparts. However, as galaxies in SIDM haloes grow in mass, baryons start to dominate the central potential. In response dark matter particles thermalise through frequent interactions and accumulate in the center of the baryonic-dominated potential. Over time, this results in an overconcentration of dark matter, manifesting as a highly cuspy density profile. This can be seen in the increasing central density of SIDM haloes in the right panel of Fig. 4. In the WeakStellarFB models, however, the situation is slightly different. Due to the early domination of baryons of the central potential, SIDM haloes quickly formed highly cuspy density profiles, with minimal evolution in the redshift range zero to two. Further details are presented in Appendix B.

In this section we have shown how baryons impact on the dark matter distribution under SIDM and CDM. While not an entirely novel result, this study presents the first cosmological simulations of a galaxy population under different velocity-dependent SIDM models and baryonic feedback schemes. The resulting features of the galaxy population have important implications for studies aiming to constrain SIDM by directly comparing to observational datasets. This is shown and discussed in the next section.

5 Tully-Fisher Relation

The galaxy sample from the TangoSIDM simulations is characterized by distinct sizes, varying from highly extended to compact, depending on the stellar feedback model (Reference versus WeakStellarFB, as illustrated in Fig. 2). Simultaneously, the sample includes haloes with distinct dark matter distributions, with SIDM haloes having densities that deviate significantly from the NFW profile, as shown in Fig. 3. The sample’s stellar-to-halo mass relation is consistent with observations (as depicted in the left panels of Fig. 2), and therefore haloes of a given mass host galaxies of the correct mass range. In this section, we test our galaxy sample with the stellar-mass Tully-Fisher relation, which establishes a correlation between the stellar mass and circular speed at a characteristic radius of spiral galaxies. First investigated by Tully & Fisher (1977), the relation has since become one of the best studied galaxy scaling relations (see e.g. Bell & de Jong 2001; Ziegler et al. 2002; Pizagno et al. 2007; Avila-Reese et al. 2008; Reyes et al. 2011; Catinella et al. 2023; Ristea et al. 2024), so that numerous studies have delved into its cosmological origin using both semi-analytical approaches and simulations (see e.g. Steinmetz & Navarro 1999; Dutton & van den Bosch 2012; Cattaneo et al. 2014; Desmond & Wechsler 2015; Ferrero et al. 2017).

Our analysis in this section demonstrates that when TangoSIDM galaxies are too compact or when dark matter is overly concentrated in the center, their rotation curves peak at much higher velocities than observed. This poses a powerful challenge for the validity of SIDM models. To quantify the significance of this constraint, Section 5.3 assesses which simulated galaxy samples, drawn from the Reference versus WeakStellarFB models under the various dark matter scenarios (presented in Section 5.1), are consistent with the observational sample (introduced in Section 5.2). This consistency test implies assessing the likelihood that the two sets of samples (simulated and observational) were drawn from the same, albeit unknown, probability distribution. Following this, in Section 5.4, we analyse the deviation of TangoSIDM galaxies from the observed Tully-Fisher relation. Subsequently, we evaluate the statistical significance of this deviation in Section 5.5.

Figure 5: Effective radius as a function of stellar mass for $z=0$ disc-type galaxies from the Reference (orange solid line) and WeakStellarFB (blue dashed line) models under CDM (left panel), SigmaConstant10 (second panel from the left), SigmaVel30 (third panel from the left) and SigmaVel60 (right panel). The observational sample is shown in grey symbols, with crosses corresponding to the SPARC dataset, triangles to the Pizagno et al. (2007) catalog and stars to the Reyes et al. (2011) data. The solid lines indicate the median relations for both the compiled observational sample (black) and the simulated sample (in color), while the shaded regions highlight the 16-84th percentiles. A visual inspection suggests that the simulated samples from the Reference model closely agree with the observational data, whereas the samples from the WeakStellarFB model do not. A statistical analysis using the Kolmogorov-Smirnov test reveals that only the massive ($M_{*}\geq 10^{10}~{}\rm{M}_{\odot}$) simulated galaxies from the Reference model under CDM, SigmaVel30 and SigmaVel60 are not significantly different from the observational sample.

5.1 Simulated sample

We create a subsample of disc-type galaxies using the fraction of stellar kinetic energy invested in ordered co-rotation, $\kappa_{\rm{co}}$, defined as

$$\kappa_{\rm{co}}=\frac{K_{\rm{co-rot}}}{K}=\frac{1}{K}\sum_{i}^{r<50\rm{kpc}}\frac{1}{2}\,m_{i}\left[L_{z,i}/(m_{i}\,R_{i})\right]^{2},$$ (4)

to quantify morphology (see e.g. Correa et al. 2017). In eq. (4), the sum is over all stellar particles within a spherical radius of 50 kpc centered on the minimum of the potential, $m_{i}$ is the mass of each stellar particle, $K(=\sum_{i}^{r<50\rm{kpc}}\frac{1}{2}m_{i}v_{i}^{2})$ the total kinetic energy, $L_{z,i}$ the particle angular momentum along the direction of the total angular momentum of the stellar component of the galaxy and $R_{i}$ is the projected distance to the axis of rotation. See also Sales et al. (2010) and Correa & Schaye (2020) for more details on $\kappa_{\rm{co}}$.

To create a disc-type galaxy subsample within each simulation, we use the criterion $\kappa_{\rm{co}}{>}0.3$ following Correa & Schaye (2020), who showed that values in the range $\kappa_{\rm{co}}{=}0.3-0.35$ select disc-type galaxies from the EAGLE simulations that agree with the distribution of disc galaxies from SDSS in the morphology-stellar mass-halo mass plane. This results in a selection of 61 disc-type galaxies per simulation with stellar masses ranging from $10^{9}~{}\rm{M}_{\odot}$ to $1.2\times 10^{11}~{}\rm{M}_{\odot}$ and effective sizes, denoted as $R_{\rm{eff}}$, ranging from 1.4 kpc to 17.3 kpc. Note that $R_{\rm{eff}}$ is defined as the 2D projected size enclosing 50 per cent of the total K-band luminosity. The total luminosity is computed from all bound star particles within a projected 2D circular aperture of 50 kpc radius. The luminosities are intrinsic (i.e. dust-free) and are calculated at each output time and for each star particle, accounting for its age, mass, and metallicity. This calculation is performed by the SWIFT code using the photometric tables from Trayford et al. (2015). Finally, we estimate the circular velocity at the effective radius, $V_{\rm{circ}}(R_{\rm{eff}})$, as follows $V_{\rm{circ}}(R_{\rm{eff}})=\sqrt{GM(<R_{\rm{eff}})/R_{\rm{eff}}}$, where the sum $M(<R_{\rm{eff}})$ considers the total mass of baryons (stars and gas) and dark matter enclosed within $R_{\rm{eff}}$.

5.2 Observational sample

We compile an observational sample by joining the catalogs of disk galaxies from Lelli et al. (2016), Pizagno et al. (2007) and Reyes et al. (2011), resulting in a dataset of 429 disc galaxies with stellar masses within the range of $10^{9}~{}\rm{M}_{\odot}$ to $2\times 10^{11}~{}\rm{M}_{\odot}$, and effective radii, $R_{\rm{eff}}$, spanning from 1.2 kpc to 18.5 kpc. Note that $R_{\rm{eff}}$ is defined as the radius encompassing half of the total galaxy luminosity. Rotational curves at $R_{\rm{eff}}$ were either directly extracted or estimated from each catalog. In the following, we provide a more detailed overview of these datasets.

Lelli et al. (2016) presented the Spitzer Photometry and Accurate Rotation Curves (SPARC) dataset, a galaxy catalog of 175 disc galaxies with near-infrared photometry at 3.6 $\mu$m and well-defined, high-quality HI rotation curves. For our analysis, we extracted inclination-corrected circular velocities, total luminosity at 3.6 $\mu$m, and effective radii directly from SPARC. We followed Lelli et al. (2017) and determined stellar masses using a constant mass-to-light ratio of $\Gamma=0.5~{}\rm{M}_{\odot}/\rm{L}_{\odot}$, which was motivated by stellar population synthesis models (Schombert & McGaugh 2014) using a Chabrier IMF. The total circular velocity at the effective radius was computed by interpolating the rotational curves.

The catalog derived by Pizagno et al. (2007) consists of 163 spiral galaxies featuring resolved $H_{\alpha}$ rotation curves. We utilized the effective radius and circular velocity at the effective radius directly from this catalog and estimated stellar masses using the $i$-band magnitudes, assuming a constant I-band mass-to-light ratio of 1.2, $M_{*}=1.2\times 10^{0.4(i_{\odot}-i)}~{}\rm{M}_{\odot}$ with $i_{\odot}=4.11$. The mass-to-light-ratio is adopted for a Chabrier IMF and it assumes the contribution of disc+bulge (Portinari et al. 2004). The effective radius for this sample is defined as the radius at $2.2\times R_{\rm{disk}}$, where $R_{\rm{disk}}$ is the disc exponential scale length.

Finally, Reyes et al. (2011) provided an improved estimate of disk rotation velocities for a subset of SDSS galaxies. This dataset includes the $i$-band Petrosian half-light radius, $r$-band Petrosian absolute magnitude ($M_{r}$), and $g-r$ colour, all $k$-corrected to $z=0$ and corrected for Galactic and internal extinction. Stellar masses were estimated following Bell et al. (2003),

$$M_{*}=10^{[\log_{10}(L_{r}/L_{r,\odot})+\log_{10}(M_{*}/L_{r})+\log_{10}h^{2}]}~{}\rm{M}_{\odot},$$ (5)

where $\log_{10}(L_{r}/L_{r,\odot})=-0.4(M_{r}-M_{r,\odot}+1.1z)$ with $M_{r,\odot}=4.76$, and $\log_{10}(M_{*}/L_{r})=-0.306+1.097\cdot(g-r)-0.093$, where the last term, $-0.093$, corresponds to the conversion from a modified Salpeter IMF to a Chabrier IMF (as indicated in Gallazzi et al. 2008).

To estimate the rotation velocity at $R_{\rm{eff}}$, we used the arctangent model

$$V_{\rm{circ,obs}}(R^{\prime})=V_{0}+\frac{2}{\pi}V_{\rm{c,obs}}\arctan\left(\frac{R^{\prime}-R_{0}}{R_{\rm{TO}}}\right).$$ (6)

Reyes et al. (2011) fitted this model to each rotational curve from the sample and provided the four free parameters: the systemic velocity $V_{0}$, the asymptotic circular velocity $V_{\rm{c,obs}}$, the spatial center $R_{0}$, and the turn-over radius $R_{\rm{TO}}$, at which the rotation curve starts to flatten out. We use the above expression for $V_{\rm{circ,obs}}(R^{\prime})$ and estimate it at $R_{\rm{eff}}$, by converting $R_{\rm{eff}}$ into arcsecond units and correcting for inclination as follows $V_{\rm{circ}}=V_{\rm{circ,obs}}(R^{\prime})/\sin(i)$.

Figure 6: The stellar mass Tully-Fisher relation, i.e. total circular velocity at the effective galactic radius as a function of stellar mass for disc-type galaxies. The median relations for disc galaxies in the Reference and WeakStellarFB models are represented by orange and blue lines, respectively, with shaded regions highlighting the $1{-}99$th percentiles. The panels show the Tully-Fisher relation for simulated galaxies under CDM (left panel), SigmaConstant10 (second panel from the left), SigmaVel30 (second panel from the right) and SigmaVel60 (right panel). Similar to Fig. 5, the panels also display the observational sample in grey symbols, with crosses corresponding to the SPARC dataset, triangles to the Pizagno et al. (2007) catalog and stars to the Reyes et al. (2011) data. The solid and dashed lines depict the best-fitting linear relations to the observational sample and simulated samples, respectively. The figure reveals a close agreement between observations and the Reference and WeakStellarFB models under CDM. However, this agreement is not maintained for SIDM. Under SIDM, the slope of the Tully-Fisher relation from the simulated sample deviates from the observed relation, with the most significant deviation occurring in the SigmaVel30 model, followed by SigmaVel60. The deviation between the relations (observational vs. simulated) is statistically significant at the $98\%$ level in the SigmaVel30 model for galaxies with masses ${\geq}10^{10}~{}\rm{M}_{\odot}$, and at the $95\%$ confidence level in the SigmaVel60 model for galaxies with masses ${\geq}1.3{\times}10^{10}~{}\rm{M}_{\odot}$.

5.3 The mass-size plane

We compare the observational sample in the mass-size plane with a sample of disc-type galaxies taken from the simulations. Fig. 5 shows the effective radius as a function of stellar mass for disc-type galaxies from the Reference (orange solid line) and WeakStellarFB (blue dashed line) models under CDM (left panel), SigmaConstant10 (second panel from the left), SigmaVel30 (third panel from the left) and SigmaVel60 (right panel). The coloured lines represent the median relations, and the shaded regions depict the 16-84th percentiles. The observational sample is shown in grey symbols, and its median relation is depicted in solid black line.

The panels in Fig. 5 indicate that the median trend of the simulated samples from the Reference model agrees with the observational data, while the simulated galaxies from the WeakStellarFB do not, as they become quite compact around a stellar mass of $10^{10}~{}\rm{M}_{\odot}$. To determine the statistical significance of the differences in the mass-size plane between the simulated and observational samples, we perform a Kolmogorov-Smirnov test (KS) for two samples. Given that the observed sample is not a volume-limited sample, we opt not to account for the mass distribution. Instead, we make a quantitative analysis by dividing the samples into bins of stellar mass and comparing the size distributions. For each stellar mass bin, we test the null hypothesis that the two samples—observational and simulated—were drawn from the same distribution. A confidence level of 95% is chosen, implying that we reject the null hypothesis if the $p$-value is less than 0.05. The aim of the KS test is to identify the simulated galaxy sample that is most likely drawn from the distribution function of the observational sample, making it statistically equivalent.

We separate the samples into three stellar mass bins ([$10^{9}$, $3{\times}10^{9}$], [$3{\times}10^{9}$, $10^{10}$], and [$10^{10}$, $3{\times}10^{10}~{}\rm{M}_{\odot}$]), and compare the observational sample and simulated galaxies from the Reference model under CDM. The KS statistical analysis returns $p$-values of 0.36, 0.02 and 0.13, respectively under each mass bin. The low $p$-value of 0.02 for simulated galaxies with stellar masses between $3{\times}10^{9}$ and $10^{10}~{}\rm{M}_{\odot}$ indicates that those galaxies do not conform to the observed size distribution, whereas galaxies in the other mass bins do. We further analyse the samples of galaxies from Reference + SigmaConstant10, SigmaVel30 and SigmaVel60, contrasting them with the observational sample. For Reference + SigmaConstant10, the analysis returns the following $p$-values of $8{\times}10^{-4}$, 0.78 and 0.15. Similarly, Reference + SigmaVel30 returns $p$-values of $8{\times}10^{-4}$, 0.31 and 0.09, whereas Reference + SigmaVel60 yields $p$-values of 0.59, $5{\times}10^{-3}$ and 0.06. For all Reference models (CDM + SIDM), the large $p$-values in the stellar mass bins $10^{10}-3{\times}10^{10}~{}\rm{M}_{\odot}$ indicate that we cannot reject the null hypothesis. Therefore, at the high mass end the samples, both observational and simulated, are not significantly different at the 95$\%$ confidence level and could be drawn from the same size distribution.

Differently, the WeakStellarFB model (under CDM or SIDM) fails to produce galaxies with sizes that agree with the observations. The KS test returns small $p$-values (${<}0.01$) for galaxies more massive than $3{\times}10^{9}~{}\rm{M}_{\odot}$. For lower mass galaxies, in the regime where the overcooling of the model has a lesser impact (as discussed in Section 3), the KS test yields $p$-values of 0.06, 0.1 and 0.98 for WeakStellarFB + SigmaConstant10, + SigmaVel30, and + SigmaVel60, respectively. From what we conclude that in the low mass end, the WeakStellarFB model under SIDM, produces galaxies whose sizes are not statistically different from the observations.

5.4 Tully-Fisher relation

The Tully-Fisher relation is shown in Fig. 6, where the y-axis corresponds to the total circular velocity at the effective galactic radius and the x-axis corresponds to the stellar mass. The left panel displays the Tully-Fisher relation for disc galaxies from the Reference (orange solid line) and WeakStellarFB (blue solid line) models under CDM. Moving from left to right, the subsequent panels show the relation for disc galaxies under the SigmaConstant10 model, SigmaVel60 model, and SigmaVel30 model. Similar to Fig. 5, the panels also show the observational sample in grey symbols. Coloured lines highlight the median relations from the simulations, while shaded regions represent the 1-99th percentiles.

The figure shows a tight correlation between circular velocity and stellar mass, as expected. This correlation is further highlighted by the best-fitting linear relation to the observational sample (black solid lines). The best-fitting parameters of the relation, $\log_{10}(V_{\rm{circ}}/{\rm{km~{}s^{-1}}})=a\log_{10}(M_{*}/10^{10}{\rm{M}}_{\odot})+b$, are $a=0.34{\pm}0.01$ and $b=2.07{\pm}0.01$. The parameters and $5-95\%$ confidence intervals were estimated by bootstrapping the observational sample. Similarly, we created a joint sample of galaxies from both the Reference and WeakStellarFB models, and via the bootstrap method we estimated the best-fitting linear relations from the simulations, which are depicted by black dashed lines in the panels.

Fig. 6 shows the close agreement between the Tully-Fisher relation derived from the observational sample and that of the simulated sample of disc galaxies from both the Reference and WeakStellarFB models under CDM. However, this agreement is not maintained when considering the SIDM models. Under the SIDM framework, the slope of the Tully-Fisher relation from the simulated sample begins to deviate relative to the observed Tully-Fisher relation. The largest deviation occurs in the SigmaVel30 model, followed by the SigmaVel60 model. The shift in $V_{\rm{circ}}(R_{\rm{eff}})$ found in galaxies within the SIDM models is attributed to the large central dark matter densities that result from the dark matter particle interactions. Consequently, at constant $M_{*}$, the increased enclosed dark matter mass drives a higher $V_{\rm{circ}}(R_{\rm{eff}})$ compared to CDM, thereby altering the slope of the relation.

Our analysis in the previous subsection established that only the disc galaxies from the Reference model under CDM, SigmaVel30 and SigmaVel60 were statistically comparable to the observational sample. This was not found for galaxies from the WeakStellarFB model under any dark matter model. Nonetheless, in Fig. 6, we intentionally include the trend from the WeakStellarFB model to highlight how the deviation from the observed Tully-Fisher relation increases under SIDM, particularly when galaxies become more compact. Notably, under CDM, the Tully-Fisher relations from both the Reference and WeakStellarFB models closely agree. This finding appears to contradict the conclusions drawn by Ferrero et al. (2017), who posited that $\Lambda$CDM models should be capable of matching the observed Tully-Fisher relation, provided that galaxy sizes are well reproduced, and that halos respond approximately adiabatically to galaxy assembly. This will be further addressed in future work, with more variations of the stellar feedback model and larger number statistics from the simulated sample.

5.5 Statistical analysis

The panels in Fig. 6 reveal a discernible departure of the Tully-Fisher relation from disc galaxies under SIDM relative to the observed Tully-Fisher relation. To quantify the significance of this deviation and to assess the likelihood of a similar deviation in the observational sample, we perform a statistical analysis focusing on the observational sample and the simulated samples from the Reference model under SigmaVel30 and SigmaVel60. Reference+SigmaConstant10 is not considered in the analysis because this particular SIDM model has already been ruled out by observations of galaxy clusters. Additionally, the WeakStellarFB models are excluded from the analysis due to their significant difference from the observations (as established in Section 5.3).

Given that the deviation in the SIDM models is prominent in massive galaxies, as demonstrated in Section 4 for haloes more massive than $10^{12}~{}\rm{M}_{\odot}$, and considering that these galaxies are statistically equivalent to the observational sample in terms of their size distribution, as demonstrated in Section 5.3, we apply a selection cut in stellar mass of $10^{10}~{}\rm{M}_{\odot}$. This allows us to analyse the Tully-Fisher relation for only massive galaxies, while disregarding the influence of low-mass systems that do not present significant changes in their central density profiles relative to CDM. This stellar mass cut results in a subsample of 287 real galaxies and 33 simulated galaxies from each SIDM model.

We perform a bootstrap analysis around these subsamples using 10,000 iterations. In each iteration, we create random samples (with replacement) for both the observational and simulated datasets, and calculate the slopes of their respective Tully-Fisher relations. After all iterations, we calculate the mean value and confidence intervals of the slopes. For the observational sample, a slope of $a=0.34{\pm}0.02$ is obtained. When comparing this slope with its value for the entire sample (calculated in Section 5.4), we find that the observed Tully-Fisher relation does not change when we consider only the subsample of massive galaxies.

For the simulated sample from the SigmaVel30 model, we find a slope of $a=0.48{\pm}0.08$, and for the SigmaVel60 model, a slope of $a=0.41{\pm}0.07$. The SigmaVel30 model exhibits a strong deviation from the observed Tully-Fisher relation, as indicated by the different slope. When we assess the differences between these slopes, we obtain a $p$-value of 0.012, which indicates the frequency that each random bootstrap sample from the simulations had a slope lower than the slope from the observational bootstrap sample. Consequently, we conclude that the Reference+SigmaVel30 model, despite producing galaxies with stellar masses and sizes in good agreement with the observations (Fig. 5), produces a Tully-Fisher relation that deviates from the observed one at the $98\%$ confidence level.

The Tully-Fisher relation from the SigmaVel60 model also deviates from the observed relation, although it is not as pronounced as in the SigmaVel30 case. The difference between these slopes yields a $p$-value of 0.13, signifying that in ${\sim}13\%$ of the bootstrap samples from the SIDM model, a slope equivalent or lower than the one derived from the observations arises. Therefore, we cannot reject the null hypothesis that both Tully-Fisher relations, from the observations and simulations, are drawn from the same distribution.

We further investigate this and calculate the minimum stellar mass above which the simulated galaxy sample from SigmaVel60 produces a Tully-Fisher relation that deviates significantly from the observed one. This cut is identified for galaxies with $M_{*}\geq 1.3{\times}10^{10}~{}\rm{M}_{\odot}$. For these refined subsamples, the observational sample yields a slope of $a=0.32{\pm}0.03$, while the SigmaVel60 model produces a slope of $a=0.45{\pm}0.09$. The bootstrap analysis yields a low $p$-value of 0.042, indicating that over this mass range, the SigmaVel60 model produces a Tully-Fisher relation that deviates from the observed one at the $95\%$ confidence level. As a control test, we assess the difference in the Tully-Fisher relations from the Reference+CDM model and observations over this mass range of ${>}1.3\times 10^{10}~{}\rm{M}_{\odot}$. We obtain a $p$-value of 0.26. Thus, we affirm that the Reference+CDM model maintains a good agreement with the observations.

The deviation of the Tully-Fisher relation from disc galaxies (under the Reference+SigmaVel30 models) relative to the observed Tully-Fisher relation is non-negligible. In this section, we have shown that it is statistically significant, which indicates that we can rule out the SigmaVel30 model over the mass range ${\gtrsim}10^{10}~{}\rm{M}_{\odot}$ with $98\%$ confidence. The rejection of this SIDM model relies on the assumption that the Reference galaxy formation model is valid, as we have demonstrated through the good agreement of the stellar mass-halo mass relation (Fig. 2), stellar masses (Fig. 8) and galaxy sizes (Fig. 5, supported by statistical analysis of Section 5.3) with observations.

The deviation of the Tully-Fisher relation, relative to observations, is driven by the impact of SIDM, which produces haloes with high central dark matter densities (Fig. 3). SIDM therefore raises $V_{\rm{circ}}(R_{\rm{eff}})$ at fixed stellar mass, and increases the slope of the relation as shown in this section. This physical effect, constrained by the Tully-Fisher relation, is ruled out. Note, however, that the rejection of the SigmaVel30 model is specific to a certain “mass range" (i.e. ${\gtrsim}10^{10}~{}\rm{M}_{\odot}$), because halo density evolution depends on the value of the cross section, which in turn depends on halo mass (Fig. 1). Therefore, only cross sections influencing the steepness of the haloes’ density throughout their evolution are ruled out. In a similar manner, we argue that the Tully-Fisher relation can be utilized to rule out the SigmaVel60 model over the mass range ${\gtrsim}1.3\times 10^{10}~{}\rm{M}_{\odot}$ with $95\%$ confidence. These findings have significant implications for the SIDM parameter space, which are discussed in the next section.

Figure 7: Momentum transfer cross section, $\sigma_{\rm{T}}/m_{\chi}$, plotted as a function of relative scattering velocity of dark matter particles. The blue solid lines shows the velocity-depenent SIDM models presented in this work, SigmaVel30 and SigmaVel60 (see Table 1 for details). The bottom x-axis indicates the relative velocity between dark matter particles, while the top x-axis indicates the typical halo mass that hosts orbits of such velocities. The shaded regions demarcate areas of the SIDM plane excluded by this work. The dark green and orange shaded regions highlight the excluded parameter space that is directly extracted from the simulations. The lighter green and orange shaded regions mark larger regions that are excluded based on the assumption that higher mass haloes under SIDM models with larger cross sections would exhibit a large deviation in the Tully-Fisher plane from the observations.

6 Discussion

6.1 SIDM parameter space

The SIDM parameter space, characterized by the self-interaction cross section as a function of the relative velocities between dark matter particles, is a topic of extensive debate. Robust constraints on the cross section on large scales (high dark matter particle velocities) have been established by studies of galaxy clusters (e.g. Randall et al. 2008; Dawson et al. 2013; Massey et al. 2015; Harvey et al. 2015; Wittman et al. 2018; Harvey et al. 2019; Sagunski et al. 2021; Andrade et al. 2022). However, the cross section for Milky Way-size galaxies and lower-mass systems remains highly uncertain. Recent proposals suggest that the cross section in dwarf-size galaxies should be as large as $100~{}\rm{cm}^{2}\rm{g}^{-1}$ (e.g. Correa 2021; Turner et al. 2021; Silverman et al. 2023; Slone et al. 2023; Yang et al. 2023), in order to address the diversity problem through halo core expansion and core collapse. At the scale of Milky Way-mass galaxies, Correa (2023) argues that the cross section should be lower than $10~{}\rm{cm}^{2}\rm{g}^{-1}$. Otherwise the frequent interactions between the Milky Way-mass systems and their satellites would lead to excessive mass loss and destruction of satellites, giving rise to unrealistic satellite populations.

This section discusses the new constraints on the SIDM parameter space presented in Section 5.5. Our work has shown that the co-evolution of baryons and dark matter self-interactions strongly impacts the evolution of galaxies. Compared with CDM, galaxies in SIDM hydrodynamical simulations not only tend to grow more extended (Section 3), but also contain enhanced dark matter central densities (Section 4). This behavior results in a deviation in the Tully-Fisher relation relative to an observational dataset (Section 5.4). In Section 5.5, we found that this deviation is statistically significant in the SigmaVel30 model for galaxies more massive than ${\gtrsim}10^{10}~{}\rm{M}_{\odot}$, and in the SigmaVel60 model for ${\gtrsim}1.3\times 10^{10}~{}\rm{M}_{\odot}$. Next, we place these constraints on the velocity-$\sigma_{\rm{T}}/m_{\chi}$ plane.

In what follows we argue that we can rule out velocity-cross section pairs that govern the evolution of haloes hosting the massive disc galaxies that significantly deviate from the observations in the Tully-Fisher plane. To identify the velocity-cross section pairs, we therefore select all disc galaxies from the Reference + SigmaVel30 and Reference + SigmaVel60 models with stellar masses larger than $10^{10}~{}\rm{M}_{\odot}$ and $1.3{\times}10^{10}~{}\rm{M}_{\odot}$, respectively. We follow the assembly histories of the haloes hosting these galaxies across the simulation snapshots until redshift 2, the redshift below which the haloes’ density profiles are well resolved and commence substantial evolution (refer to Fig. 4, Section 4 and Appendix B). We determine the median and 16-84th percentiles of their mass accretion histories, $M_{200}(z)$, and convert these to the circular velocity, $V_{\rm{circ}}(z)$. We assume that $V_{\rm{circ}}(z)$ is the average velocity of the dark matter particles within these haloes over the redshift range 0-2, and using eqs. (1) and (2) we estimate the haloes’ average dark matter cross section. The individual haloes’ $M_{200}(z)$, $V_{\rm{circ}}(z)$ and cross sections are shown in Appendix C.

The derived velocity-cross section pairs establish the limits above which the SigmaVel30 and SigmaVel60 models produce overly enhanced central dark matter densities in massive disc galaxies. Therefore, we mark these limits as regions where the SigmaVel30 and SigmaVel60 models are ruled out with $98\%$ and $95\%$ confidence, as shown in the green and orange shaded areas in Fig. 7. The figure depicts the momentum transfer cross section, $\sigma_{\rm{T}}/m_{\chi}$, as a function of relative dark matter particle scattering velocity. The curves show the velocity-dependent SIDM models presented in this work, SigmaVel30 and SigmaVel60 (see Table 1). While the bottom x-axis highlights the relative velocity between dark matter particles, the top x-axis indicates the typical halo mass that hosts orbits of such velocities. The dark green and orange shaded regions highlight the newly excluded parameter space that is directly extracted from the simulations. The lighter green and orange shaded regions mark further excluded regions under the assumption that higher mass haloes under SIDM models with higher cross sections would exhibit a large deviation in the Tully-Fisher plane from the observations.

While this finding imposes strong constraints on velocity-dependent models, it does not entirely rule them out. There is still room for models where the cross section reaches $100~{}\rm{cm}^{2}\rm{g}^{-1}$ at 10 km s${}^{-1}$, provided that it decreases to less than $1~{}\rm{cm}^{2}\rm{g}^{-1}$ at 150 km s${}^{-1}$. In Fig. 7, we refrain from extending the SIDM parameter space to velocities larger than 500 km s${}^{-1}$, since those are not covered by the simulations. Our future plans include expanding this analysis to larger scales, employing larger cosmological boxes and more statistical power through increased numerical resolution to model a more extensive sample with lower mass disc galaxies. Additionally, we aim to explore the circular velocities of dwarf galaxies in more detail. We anticipate that with sufficient resolution and statistics, the modelling of dwarf galaxies, even with the inclusion of baryons as shown in the bottom panels of Fig. 3, may yield lower values of $V_{\rm{circ}}(R_{\rm{eff}})$ relative to an observational sample at fixed stellar mass, consequently resulting in a deviation of the Tully-Fisher relation. This analysis, coupled with methodology improvements such as mock observations of HI discs for extracting rotational curves, will be the focus of future work.

7 Conclusions

The SIDM parameter space, while extensively explored in recent years, remains notably uncertain, particular for Milky Way-size galaxies and smaller systems. This uncertainty arises due to the inherent challenge of isolating the impact of baryonic physics from dark matter interactions. Recent studies (e.g. Robertson et al. 2019; Despali et al. 2019; Sameie et al. 2021; Rose et al. 2022; Burger et al. 2022; Jiang et al. 2023) have reported that the prevalence of baryons in the central gravitational potential leads to the formation of denser and more cusp-like central density profiles under SIDM compared to CDM. Nevertheless, uncertainties persist regarding how the increased cuspiness of SIDM haloes correlates with the specific SIDM model parameters and the strength of galaxy feedback. To address these uncertainties, this study introduces a new set of cosmological hydrodynamical simulations. These simulations include the SIDM model derived from the TangoSIDM project (Correa et al. 2022) and leverage the baryonic physics from the SWIFT-EAGLE galaxy formation model (Borrow et al. 2023; Schaller et al. 2023).

Two cosmological simulation suites were generated: The Reference model, calibrated in a (25 Mpc)${}^{3}$ volume to reproduce the galaxy stellar mass function and galaxy mass-size relation; and the WeakStellarFB model, featuring less efficient stellar feedback around Milky Way-like systems. Each galaxy formation model (Reference and WeakStellarFB) was simulated under four dark matter cosmologies: CDM, SigmaConstant10 (a SIDM model with a constant cross section of 10 cm${}^{2}\rm{g}^{-1}$), and SigmaVel30 and SigmaVel60, two SIDM models with velocity-dependent cross sections (see Fig. 1). SigmaVel60 has a cross section smaller than 1 cm${}^{2}\rm{g}^{-1}$ at high velocities ($v{>}150$ km s${}^{-1}$) and increases with decreasing velocity, reaching 60 cm${}^{2}\rm{g}^{-1}$ at 10 km s${}^{-1}$. SigmaVel30 has a cross section smaller than 8 cm${}^{2}\rm{g}^{-1}$ at velocities surpassing $200$ km s${}^{-1}$ (dropping below 1 cm${}^{2}\rm{g}^{-1}$ at ${\approx}1000$ km s${}^{-1}$) and it also increases with decreasing velocity. These SIDM models we selected to represent two extreme scenarios for the rate of dark matter interactions in Milky Way-mass systems. The SWIFT-EAGLE models were selected to determine whether the impact of SIDM on galaxies remains robust when subjected to variations in feedback models.

Our findings indicate that SIDM does not significantly alter global galaxy properties such as stellar masses and star formation rates, but it does impact galaxy sizes, making galaxies more extended (Fig. 2). Dark matter particle interactions heat the inner halo, leading to core formation in the central regions of haloes less massive than $10^{11}~{}\rm{M}_{\odot}$ and dynamically heating the surrounding gas and stars, promoting the formation of more extended galaxies. However, we have found that the impact of SIDM is insufficient to counteract the gas overcooling and size compactness in galaxies from the WeakStellarFB model.

In massive haloes (${\sim}10^{12}~{}\rm{M}_{\odot}$), baryonic influence on SIDM distributions result in steeper dark matter density profiles than those produced in CDM from adiabatic contraction (Fig. 3). This feature is enhanced in the WeakStellarFB model, suggesting that the increased baryonic concentration in the model, relative to the Reference model, enhances the central concentration of the dark matter distribution. Under SIDM, the haloes density profile evolved differently (Fig. 4). As galaxies grow in mass, baryons begin to dominate the central gravitational potential, causing dark matter particles to thermalise through frequent interactions and accumulate in the center, resulting in cuspy dark matter density profiles.

The enhanced dark matter density at the centers of galaxies results in a notable deviation in the slope of the Tully-Fisher relation, significantly diverging from observations. We assembled an observational sample of $z\approx 0$ disc galaxies by combining the catalogs from Pizagno et al. (2007), Reyes et al. (2011) and Lelli et al. (2016). Our analysis reveals that while the simulated massive galaxies from the Reference model under SigmaVel30 and SigmaVel60 are not significantly different from the observational sample in the galaxy mass-size plane (Fig. 5), they strongly deviate in the Tully-Fisher plane (Fig. 6). This is due to a shift in $V_{\rm{circ}}(R_{\rm{eff}})$ found in galaxies within the SIDM models, driven by the large central dark matter densities that result from the dark matter particle interactions. Consequently, at constant $M_{*}$, the increased enclosed dark matter mass leads to a higher $V_{\rm{circ}}(R_{\rm{eff}})$ compared to CDM, altering the slope of the relation. In contrast, the Tully-Fisher relation derived from CDM models aligns well with observations.

We have conducted a statistical analysis to assess the significance of the discrepancy between the SIDM models and observations in the Tully-Fisher plane. Our findings indicate that galaxies from the Reference+SigmaVel30 model more massive than $10^{10}~{}\rm{M}_{\odot}$ deviate from the observational sample at the 98% confidence level, while galaxies with masses exceeding $1.3{\times}10^{10}~{}\rm{M}_{\odot}$ from the Reference+SigmaVel60 model deviate at the 95% confidence level. These constraints, when translated into the velocity-$\sigma_{\rm{T}}/m_{\chi}$ plane (Fig. 7), reveal that the cross section should be smaller than 0.5 cm${}^{2}\rm{g}^{-1}$ for velocities of $\sim$150-200 km s${}^{-1}$ and smaller than 10 cm${}^{2}\rm{g}^{-1}$ for velocities of 110-180 km s${}^{-1}$.

Our study reveals that the Tully-Fisher plane, encompassing galaxy sizes, stellar masses, and circular velocities, serves as a powerful observable for discerning and excluding velocity-dependent SIDM models. In future work we will focus on improving the datasets (higher numerical resolution and larger cosmological box size for the simulations, as well as larger data compilation from observational surveys) and refining the methodology, including the creation of mock HI rotation curves, with the goal of carrying out more accurate and precise comparisons.

Data availability

The data supporting the plots within this article will be made publicly available in http://www.tangosidm.com. See example in Table 3.

Appendix A Galaxy Stellar Mass Function

Figure 8: The galaxy stellar mass function at $z=0$ for the ‘WeakStellarFB’ (left panel) and ‘Reference’ (right panel) galaxy formation models under the CDM (blue lines), SigmaVel60 (purple line), SigmaVel30 (red line) and SigmaConstant10 (orange line) schemes. The black line corresponds to the galaxy stellar mass function of the original EAGLE 25 REF model and the green line corresponds the measurements from the DR4 GAMA survey at $z<0.1$ (Driver et al. 2022). Both models, Reference and WeakStellarFB, produce a galaxy number density in the stellar mass range $10^{8}-10^{11}~{}\rm{M}_{\odot}$ that is in agreement with EAGLE and the observational data within 0.2 dex. Note that the models considered here were calibrated to reproduce the $z=0$ galaxy stellar mass function.

Fig. 8 shows the $z=0$ galaxy stellar mass function for the WeakStellarFB (left panel) and Reference (right panel) galaxy formation models under the CDM (blue lines), SigmaVel60 (purple line), SigmaVel30 (red line) and SigmaConstant10 (orange line) schemes. The simulation results are compared to the original EAGLE REF model (Schaye et al. 2015), and to the DR4 Galaxy And Mass Assembly (GAMA) survey (Driver et al. 2022). The EAGLE data shown throughout this section is taken from the EAGLE reference model run in a (25 Mpc)${}^{3}$ box with the same resolution as the TangoSIDM simulations, which were also run in a (25 Mpc)${}^{3}$ volume.

Both Reference and WeakStellarFB produce a galaxy number density in the stellar mass range $10^{8}-10^{11}~{}\rm{M}_{\odot}$ that is in close agreement with EAGLE and within 0.2 dex of the observational data. While Fig. 8 seems to indicate that SIDM does not strongly affect the galaxy stellar mass function, it does decrease the number of satellites (as shown in Vogelsberger et al. 2012; Nadler et al. 2020; Correa et al. 2022). SIDM interactions enhance the disruption of subhaloes by tidal stripping from the host. We find that from the 685 satellite galaxies in the stellar mass range $10^{7-10}\rm{M}_{\odot}$ from the Reference/CDM model, 639 (93%) survive in the Reference/SigmaVel60 model and 544 (79%) survive in the Reference/SigmaConstant10 model.

Appendix B Density evolution

This appendix expands the discussion presented in Section 3, where we showed that under SIDM halo dark matter density profiles evolve differently than under CDM (Fig. 4). We have found that as galaxies within SIDM haloes grow in mass, baryons assume a dominant role in the galaxies’ central gravitational potential. Consequently, dark matter particles thermalise through frequent interactions, accumulating in the center of the baryon-dominated potential. Fig. 9 shows the density evolution of the 32 most massive haloes from the WeakStellarFB model under CDM (left panel) and SigmaVel60 (middle panel), and from the Reference model under SigmaVel60 (right panel). The coloured lines represent the median density evolution between redshifts 0 and 2. In the WeakStellarFB models, the early dominance of baryons in the central potential results in the rapid formation highly cuspy density profiles, which for SIDM remains with minimal evolution in the redshift range zero to two. In contrast, under CDM haloes there is a slight decrease in cuspiness by redshift zero. The right panel of Fig. 9 demonstrates that, in the SigmaVel60/Reference model, the median central density of haloes slightly increases over time, as was the case for the SigmaVel30/Reference model.

Figure 9: Stacked median dark matter density profile, $\rho_{\rm{DM}}$, for the 32 most massive haloes from the WeakStellarFB model under CDM (left panel) and SigmaVel60 (middle panel), and from the Reference model under SigmaVel60 (right panel). The coloured lines show the median values at different redshifts, and the black solid line shows the NFW profile of the haloes at redshift zero. The black dashed-dotted lines indicate the convergence radius (see Section 4 for the definition). The median density profiles of haloes over time in WeakStellarFB/CDM are cuspy by redshift 2 and slightly decrease by redshift zero. For the WeakStellarFB/SigmaVel60 model, the median density profiles of the haloes is cuspier by redshift 2, and they do not largely evolve towards redshift zero.

Appendix C Assembly history

Section 6.1 reported an important discrepancy found in massive disc galaxies within the SIDM framework when compared to observations in the Tully-Fisher plane. This discrepancy was translated into an exclusion zone within the SIDM parameter space. Our approach involved identifying velocity-cross section pairs that lead to the formation of galaxies with exceedingly large $V_{\rm{circ}}(R_{\rm{eff}})$. In this section, we provide further details on the methodology employed to determine the lower limits for velocity and cross section, above which the SigmaVel30 and SigmaVel60 models are ruled out.

To identify these velocity-cross section pairs, we select all disc galaxies from the Reference + SigmaVel30 and Reference + SigmaVel60 models with stellar masses larger than $10^{10}~{}\rm{M}_{\odot}$ and $1.3{\times}10^{10}~{}\rm{M}_{\odot}$, respectively. We follow the assembly histories of the haloes hosting these galaxies across the simulation snapshots until redshift 2 (the redshift below which the haloes’ density profiles are well resolved and commence substantial evolution). The left panel of Fig. 10 shows the mass accretion history $M_{200}(z)$, of the haloes from the SigmaVel30 (dark blue lines) and SigmaVel60 (light blue lines) models under the Reference galaxy model. Converting $M_{200}(z)$ into circular velocity, $V_{\rm{circ}}(z)$, we show these values in the second form the left panel of Fig. 10.

We assume that $V_{\rm{circ}}(z)$ corresponds to the average velocity of the dark matter particles within these haloes. Therefore, to estimate the corresponding average dark matter particle cross sections of these haloes, we use eq. (1), assume $v=V_{\rm{circ}}(z)$ and integrate over the scattering angle (as done in eq. 2). The evolution of the dark matter haloes’ average cross sections, $\sigma_{T}/m_{\chi}$, as a function of redshift is shown in the second panel from the right. The right most panel displays the cross section as a function of the haloes circular velocities for the SigmaVel30 and SigmaVel60 models (grey lines). As expected, all the velocity-cross section pairs that were obtained from the haloes’ evolution align with the velocity-$\sigma_{T}/m_{\chi}$ relation from the models (eq. 2).

In the last step, at each redshift we determine the 16-84th percentiles in the distribution of the haloes circular velocities. We highlight these percentage ranges in red in the right panel and mark them as the limits above which the SigmaVel30 and SigmaVel60 models produce overly enhanced central dark matter densities in massive disc galaxies. Therefore, these limits represent the lower bounds above which the SigmaVel30 and SigmaVel60 models are ruled out with $98\%$ and $95\%$ confidence, respectively.

Figure 10: Left panel: Mass assembly, $M_{200}(z)$, illustrating the evolution of haloes hosting galaxies that exhibit significant deviations from observations in the $z=0$ Tully-Fisher plane. Dark blue lines represent the mass evolution of individual haloes from the SigmaVel30 model, while light blue lines correspond to the SigmaVel60 model. Second panel from the left: Circular velocity, $V_{\rm{circ}}(z)$, as a function of redshift for the same haloes as in the left panel. Second panel from the right: Cross section, computed from the circular velocity, as a function of redshift. Right panel: Velocity-cross section plane. The grey dashed lines mark the cross section dependence for the SigmaVel30 and SigmaVel60 models. The red limits mark the 16th-84th percentiles of the haloes’ evolution in the cross section-velocity plane.

Table 3: Observational data used in this work. Column 2 provides the sample the galaxy belongs to: ‘S’ (SPARC, Lelli et al. 2016), ‘R’ (Reyes et al. 2011), ‘P’ (Pizagno et al. 2007). Note that for the Reyes et al. and Pizagno et al. datasets, the galaxy names correspond to their SDSS names. The complete table can be found online in http://www.tangosidm.com.

Name	Sample	$M_{*}$ [M${}_{\odot}$]	$R_{\rm{eff}}$ [kpc]	$V_{\rm{circ}}(R_{\rm{eff}})$ [km s${}^{-1}$]
ESO079-G014	S	2.59e+10	7.23	140.99
ESO116-G012	S	2.15e+09	2.75	80.63
ESO563-G021	S	1.56e+11	10.59	294.74
F568-3	S	4.17e+09	7.47	91.87
F568-V1	S	1.91e+09	4.40	101.01
J001006.61-002609.7	R	9.64e+09	2.42	94.86
J001708.75-005728.9	R	4.57e+09	3.13	107.83
J002844.82+160058.8	R	2.91e+10	6.08	106.65
J003112.09-002426.4	R	1.53e+10	2.00	138.94
J004916.23+154821.0	R	7.49e+09	5.65	107.57
J004935.71+010655.2	R	3.72e+10	5.40	117.47
J011750.26+133026.3	R	3.80e+09	3.83	65.96
J012317.00-005421.6	R	1.71e+10	2.33	137.88
J012340.12+004056.4	R	2.14e+10	3.01	156.31
J012438.08-000346.4	P	2.07e+10	6.69	161.71
J013142.14-005559.9	P	6.72e+10	13.90	225.19
J013600.15+003948.6	P	3.10e+10	6.15	179.73
J013752.69+010234.8	P	5.29e+10	8.18	277.10
J014121.94+002215.7	P	1.88e+10	3.27	195.82
J015746.24-011229.9	P	8.53e+10	7.19	310.80
J015840.93+003145.2	P	4.78e+10	9.46	189.52
…	.	…	…	…