ROMULUS25 Cosmological Simulation

• ROMULUS25 is a large-scale, high-resolution cosmological simulation that models galaxy formation, SMBH dynamics, and cosmic structure evolution in a 25 Mpc volume.
• It utilizes advanced physical modules including star formation, supernova feedback, and detailed SMBH seeding and accretion models to reproduce observed galaxy scaling relations.
• The simulation robustly captures merger-driven morphological changes, environmental quenching of satellites, and SMBH triple dynamics, offering valuable insights for galaxy evolution research.

The ROMULUS25 cosmological simulation is a large-scale, high-resolution suite designed to model galaxy formation, supermassive black hole (SMBH) growth, and cosmic structure evolution in a 25 Mpc (comoving) cubic volume. It integrates advanced astrophysical prescriptions for star formation, SMBH dynamics, feedback, and baryonic processes using the ChaNGa Tree+SPH code, aiming to reproduce a broad range of observed galaxy scaling relations and the co-evolution of galaxies and their central SMBHs from high redshift to the present.

1. Simulation Framework and Core Physical Modules

ROMULUS25 utilizes the massively parallel ChaNGa code, which implements a combination of Tree gravity and Smoothed Particle Hydrodynamics (SPH). The code includes improvements such as thermal diffusion, geometric mean SPH force calculations to reduce numerical surface tension, and an optimized blastwave model for supernova (SN) feedback. The simulation tracks dark matter and gas at high mass resolution: DM particle mass is $\sim3.39\times10^5\, M_\odot$, gas particle mass is $\sim2.12\times10^5\, M_\odot$, with a Plummer-equivalent softening of 250 pc. The dark matter is oversampled relative to baryons (typically by a factor $\sim3.4$) for enhanced dynamical tracking.

Key baryonic physics modules:

• Star formation employs a stochastic, local dynamical time prescription,

$$p = \left(\frac{m_{\mathrm{gas}}}{m_{\mathrm{star}}}\right)\left[1 - \exp\left(-c_*\,\frac{\Delta t}{t_{\mathrm{dyn}}}\right)\right]$$

with $c_*=0.15$, and threshold density $n_* = 0.2$ cm$^{-3}$.

• SN feedback applies the blastwave formalism, coupling $75\%$ of SN energy ($\epsilon_{\mathrm{SN}} = 0.75$) to nearby gas.
• Metal cooling, UV background, and self-shielding are treated with up-to-date rate tables and cross-sections.

2. Supermassive Black Hole Seeding, Dynamics, and Accretion

ROMULUS25 advances previous approaches by coupling SMBH seeding and dynamics directly to local gas properties rather than imposing a fixed halo mass threshold. Seeds form when gas conditions satisfy:

• Metallicty $Z < 3\times10^{-4}$,
• Density $n > 3$ cm$^{-3}$ (15$\times$ SF threshold),
• Temperature $9500\,\mathrm{K} < T < 10000\,\mathrm{K}$.

The simulation incorporates a physically motivated dynamical friction module based on Chandrasekhar's formula,

$$a_{\mathrm{DF}} = -4\pi G^2 M\, \rho(<v_{\mathrm{BH}})\, \ln\Lambda\, \frac{v_{\mathrm{BH}}}{v_{\mathrm{BH}}^3}$$

with $\ln\Lambda$ set to the ratio of the gravitational softening to the BH influence radius. This allows for realistic SMBH orbital decay and delays in nuclear coalescence, crucial for dual and triple SMBH events.

Accretion rates are governed by a modified Bondi–Hoyle–Lyttleton formalism accounting for angular momentum:

$\dot{M} = \alpha \begin{cases} \frac{\pi G^2 M_{\mathrm{BH}}^2 \rho}{(v_{\mathrm{bulk}}^2 + c_s^2)^{3/2}} & v_{\mathrm{bulk}} > v_\theta\[2ex] \frac{\pi G^2 M_{\mathrm{BH}}^2 \rho c_s}{(v_\theta^2 + c_s^2)^2} & v_{\mathrm{bulk}} < v_\theta \end{cases}$

Here, $\alpha = (n / n_{*,SF})^\beta$ for $n \geq n_{*,SF}$, $\beta = 2$ (sub-grid boost), $v_\theta$ is the mean tangential velocity, and $n_{*,SF}$ is the star formation threshold density.

3. Feedback, Quenching, and Star Formation Histories

ROMULUS25 demonstrates that AGN (SMBH) feedback—implemented as isotropic thermal deposition of energy into the nearest $N=32$ gas particles [$E = \epsilon_r \epsilon_f \dot{M} c^2 \Delta t$ with $\epsilon_r = 0.1$, $\epsilon_f = 0.02$]—effectively regulates star formation, especially in massive galaxies ($M_{\rm halo} \gtrsim 10^{12} M_\odot$), and is significantly more effective compared to SNe feedback alone. Galaxies in ROMULUS25 reproduce both the observed cosmic star formation history (peak at $z \sim 2$), the SMHM relation, and the $M_{\mathrm{BH}}$–$M_{*}$ relation with low scatter at the high-mass end, while allowing for substantial variance in dwarfs ($M_{*} < 10^{10} M_\odot$) (Tremmel et al., 2016, Sharma et al., 2019).

Feedback from overmassive SMBHs in dwarf galaxies is shown to dominate the integrated energy injection compared to SNe, suppressing HI content, lowering central stellar densities, and quenching star formation in systems above $M_{*} \gtrsim 10^9 M_\odot$ (Sharma et al., 2019).

4. Merger-Driven Morphologies, Disk Spin, and Low Surface Brightness Galaxies

ROMULUS25 reveals that classical Low Surface Brightness (LSB) galaxies—comprising approximately 60% of central galaxies with $8 \leq \log_{10}(M_*/M_\odot) \leq 10$—are preferentially formed via major mergers in which the secondary is co-rotating and highly aligned with the primary's disk. These mergers impart excess orbital angular momentum such that the final disk exhibits large scale length $r_d$ and low central surface brightness, well-described by an exponential profile:

$\mu(r) = \mu_0 + 1.086 \frac{r}{r_d}$

LSB galaxies show nearly constant, but distributed, star formation rates, large HI gas reservoirs, and lower metallicities (Wright et al., 28 Jul 2025). Their bulges are more diffuse and optically redder than those of HSB galaxies, though lower extinction renders the integrated colors bluer.

Ultra-Diffuse Galaxies (UDGs) are also frequently realized via early major mergers (typically before $z \sim 1$), whose spin-up temporarily increases the host angular momentum by up to $\sim200\%$, displacing star formation to outer regions and aging the central stellar populations (Wright et al., 2020). Both LSB and UDG formation mechanisms are internal and merger-driven, not requiring atypical halo spins or environmental quenching.

5. Satellite Populations and Environmental Effects

Milky Way analogs selected from ROMULUS25 display satellite luminosity functions, HI content, and quenched fractions consistent with observed samples (SAGA, ELVES). Satellite number scales predominantly with host stellar mass, while environmental factors—such as proximity to another MW-mass halo—can increase quenched fractions, especially for paired systems ($<1$ Mpc separation). Quenching is determined by sSFR $< 10^{-11}$ yr$^{-1}$ and is sensitive to surface brightness limits (Nest et al., 2023).

Specific frequency definitions for satellites ($S_{n,mass}$, $S_{n,env}$) quantitatively link satellite counts to host properties and environmental isolation. ROMULUS25 further enables assessment of satellite quenching mechanisms and the role of halo–environment interactions in shaping the galaxy population.

6. SMBH Triple Dynamics and Triaxiality

High-resolution follow-up N-body studies of triple SMBH systems with ROMULUS25 initial conditions demonstrate that, in all configurations, the two most massive SMBHs consistently form rapidly hardening binaries that coalesce within $<H_0^{-1}$ via stellar hardening and GW emission. The third SMBH typically remains on a wide orbit or forms a stable hierarchical triple system (Koehn et al., 2023, Saha et al., 23 Oct 2025). The binary hardening and GW-driven merger rates are described by coupled differential equations:

$$\frac{d}{dt}\left(\frac{1}{a}\right) = s_{\mathrm{SH}} + \frac{64}{5} \frac{G^3\mu M^2}{c^5 a^5}\frac{1 + \frac{73}{24}e^2+\frac{37}{96}e^4}{(1-e^2)^{7/2}}$$

$$\frac{de}{dt} = -\frac{304}{15} \frac{G^3 \mu M^2}{c^5 a^4} \frac{e + \frac{121}{304}e^3}{(1-e^2)^{5/2}}$$

where $s_{\mathrm{SH}}$ is the stellar hardening rate extracted from N-body integration.

The simulations confirm that triaxiality (non-spherical potentials) does not qualitatively change the merger outcomes—binary formation and hardening rates remain robust, and hierarchical triple configurations are favored when the central stellar density is sufficiently cuspy. These dynamics underlie the assembly of massive SMBHs in the ΛCDM paradigm (Saha et al., 23 Oct 2025).

7. Numerical Robustness, Data Accessibility, and Community Relevance

ROMULUS25 exploits both SPH and AMR approaches (see GADGET, Enzo tests) to validate physical results. Initial conditions generators are required to self-consistently treat baryonic pressure, temperature fluctuations, and dark matter–baryon streaming velocities; imprecise initial setups can yield artificially enhanced central densities and premature star formation (O'Leary et al., 2012). The adopted methodology in ROMULUS25 is designed to mitigate such systemic biases.

By analogy with Skies and Universes data-sharing models (Klypin et al., 2017), ROMULUS25 simulation products—if published—would facilitate wide accessibility, cross-survey catalog construction, and comprehensive statistical analyses relevant for galaxy formation, SMBH growth, and cosmological inference.

---

ROMULUS25 stands as a technically advanced cosmological simulation supporting investigations across SMBH physics, galaxy evolution, and environment-dependent processes, whose detailed feedback models and merger-driven morphological mechanisms have produced results highly consistent with a wide range of observational constraints.