Romulus25: Cosmological Simulation Overview

• Romulus25 is a flagship 25 Mpc cosmological hydrodynamical simulation featuring high mass resolution, detailed SMBH seeding, and sub-kpc force resolution.
• It leverages an optimized ChaNGa tree+SPH code with explicit subgrid dynamical friction and angular-momentum-aware accretion to capture galaxy and SMBH co-evolution.
• Romulus25 supports diverse analyses including off-center SMBH dynamics, dual AGN formation, gravitational-wave host identification, and ultra-diffuse galaxy evolution.

Romulus25, often stylized as ROMULUS25, is the flagship 25 Mpc-on-a-side cosmological hydrodynamical simulation of the Romulus suite, evolved with the ChaNGa tree+SPH code to study galaxy formation together with supermassive black hole (SMBH) formation, dynamics, accretion, and feedback in a single uniform, periodic volume. In the literature, its distinctive role is the combination of $\sim 10^5\,M_\odot$ mass resolution, sub-kpc force resolution, local $10^6\,M_\odot$ black-hole seeding, angular-momentum-aware Bondi-like accretion, and explicit subgrid dynamical friction, which together enable analyses ranging from SMBH–galaxy co-evolution and wandering SMBHs to dual AGN, gravitational-wave progenitors, ultra-diffuse galaxies, low-surface-brightness systems, and Milky Way analog satellites (Tremmel et al., 2016).

1. Numerical architecture

Romulus25 is described as a 25 Mpc cosmological volume run to low redshift with ChaNGa, using modern SPH, metal and thermal diffusion, a cosmic UV background, low-temperature metal cooling, star formation, stellar feedback, and SMBH physics. Across the Romulus25 literature, the numerical parameters most consistently quoted are a dark-matter particle mass of $3.39\times10^5\,M_\odot$, a gas particle mass of $2.12\times10^5\,M_\odot$, and a gravitational softening reported as 350 pc spline or 250 pc Plummer-equivalent, with gravity becoming Newtonian at $\sim 700$ pc (Tremmel et al., 2016).

Quantity	Reported value	Context
Box size	25 Mpc on a side	Uniform cosmological volume
Dark matter particle mass	$3.39\times10^5\,M_\odot$	Standard Romulus25 resolution
Gas particle mass	$2.12\times10^5\,M_\odot$	Standard Romulus25 resolution
Star particle mass	$6\times10^4\,M_\odot$	$0.3\,m_{\rm gas}$
Force softening	350 pc spline; 250 pc Plummer-equivalent	Newtonian beyond $\sim700$ pc

The baryonic model couples a stochastic star-formation law to blastwave supernova feedback. Gas forms stars when $10^6\,M_\odot$0 K and $10^6\,M_\odot$1, with probability

$10^6\,M_\odot$2

where $10^6\,M_\odot$3, $10^6\,M_\odot$4, and $10^6\,M_\odot$5 yr. Type II supernovae inject $10^6\,M_\odot$6 erg with 75% coupling in the blastwave prescription; Type Ia supernovae, stellar winds, metal diffusion, and UV-background heating are also included (Wright et al., 28 Jul 2025).

A notable methodological element is the parameter optimization described for the Romulus suite: the free subgrid parameters regulating star formation and feedback were tuned through a multidimensional Kriging-based search on zoom-in runs, scored against several $10^6\,M_\odot$7 observational constraints including stellar mass–halo mass, HI gas fraction, specific angular momentum, and SMBH mass versus stellar mass (Tremmel et al., 2016). This calibration framework is central to how Romulus25 was positioned as a simulation intended simultaneously for galaxy structure and SMBH demography.

2. SMBH formation, accretion, feedback, and dynamics

Romulus25’s SMBH model is defined by local gas conditions rather than a global halo-mass seeding threshold in its core implementation. A gas particle eligible to form a star becomes a $10^6\,M_\odot$8 SMBH seed if it is metal-poor, dense, and warm, typically quoted as $10^6\,M_\odot$9, $3.39\times10^5\,M_\odot$0, and $3.39\times10^5\,M_\odot$1 K. The literature further notes that the vast majority of seeds form at $3.39\times10^5\,M_\odot$2, within the first Gyr of cosmic time (Tremmel et al., 2016).

Accretion follows a modified Bondi–Hoyle–Lyttleton prescription with angular-momentum support and a density boost. In one commonly quoted form,

$3.39\times10^5\,M_\odot$3

with

$3.39\times10^5\,M_\odot$4

and accretion capped at the Eddington rate (Ricarte et al., 2019).

Feedback is thermal and isotropic. The injected energy over a SMBH timestep is

$3.39\times10^5\,M_\odot$5

with $3.39\times10^5\,M_\odot$6 and $3.39\times10^5\,M_\odot$7, typically distributed to the nearest 32 gas particles together with a brief cooling shut-off (Sharma et al., 2019). Romulus25 therefore adopts a single-mode thermal AGN-feedback channel in its foundational formulation.

The other defining feature is orbital dynamics. Romulus25 does not artificially pin SMBHs to galaxy centers; instead it applies an explicit Chandrasekhar-like subgrid dynamical-friction acceleration,

$3.39\times10^5\,M_\odot$8

with the Coulomb logarithm set by impact parameters tied to the force softening and the local deflection scale. This choice is what enables self-consistent tracking of long-lived off-center SMBHs, delayed pair formation, and sub-kpc orbital decay (Tremmel et al., 2016).

Within the cosmological run, two SMBHs merge numerically when their separation is below $3.39\times10^5\,M_\odot$9 kpc and they satisfy a binding criterion written as

$2.12\times10^5\,M_\odot$0

This marks the endpoint of the galaxy-scale orbital-decay phase in Romulus25; binary hardening below $2.12\times10^5\,M_\odot$1 pc is not directly resolved in the cosmological simulation itself (Tremmel et al., 2017).

3. Black hole–galaxy co-evolution

A major scientific use of Romulus25 is the study of average SMBH growth relative to host-galaxy assembly. In star-forming galaxies with stellar masses between $2.12\times10^5\,M_\odot$2 and $2.12\times10^5\,M_\odot$3, the black-hole accretion rate (BHAR) tracks the star-formation rate (SFR). When averaged over $2.12\times10^5\,M_\odot$4 Myr, the relation is reported as

$2.12\times10^5\,M_\odot$5

with $2.12\times10^5\,M_\odot$6 scatter $2.12\times10^5\,M_\odot$7 dex, equivalent to $2.12\times10^5\,M_\odot$8. The same study states that this relation shows no measurable evolution in slope or normalization from $2.12\times10^5\,M_\odot$9 out to $\sim 700$0–6 and no systematic offset across halo masses $\sim 700$1–$\sim 700$2 or between field and cluster/proto-cluster environments (Ricarte et al., 2019).

Romulus25 has also been used to test whether mergers are the dominant trigger of SMBH growth. In that analysis, cross-correlations of BHAR with time since or to major mergers show no enhancement or delay on $\sim 700$3 up to 0.5 Gyr, leading to the stated result that there is no evidence for a connection between black-hole growth and galaxy mergers, on any timescale and at any redshift, within the limitations of the simulation volume and resolution (Ricarte et al., 2019). The same work gives a three-parameter regression involving SFR, $\sim 700$4, and cold-gas fraction that explains $\sim 700$5 of the BHAR variance.

At dwarf-galaxy masses, Romulus25 predicts large scatter in $\sim 700$6. In isolated central dwarfs with $\sim 700$7, “overmassive” and “undermassive” SMBHs follow different evolutionary tracks. Overmassive BHs grow through a mixture of BH mergers and relatively high average accretion rates, begin suppressing star formation around $\sim 700$8, and by $\sim 700$9 are associated with lower central stellar mass density, lower HI content, and lower SFR than undermassive counterparts, especially above $3.39\times10^5\,M_\odot$0. The quoted quenched fraction is $3.39\times10^5\,M_\odot$1 for overmassive-BH hosts, versus 0% for the undermassive bin in the summarized table (Sharma et al., 2019).

The dwarf-BH census is extended by x-ray selection studies. Romulus25 predicts that the local MBH occupation fraction drops below unity at $3.39\times10^5\,M_\odot$2 and $3.39\times10^5\,M_\odot$3, and that 76% of all MBHs in local dwarf galaxies are “hidden” because their x-ray luminosities are too low or too contaminated by x-ray binaries and hot ISM emission for current surveys (Sharma et al., 2022). This result is paired with the claim that local dwarf AGN in Romulus25 follow observed scaling relations between AGN x-ray luminosity, stellar mass, and SFR, while showing slightly higher active fractions and number densities than comparable x-ray observations.

Romulus25 has also been used to connect SMBH mass scatter to circumgalactic metal transport. For $3.39\times10^5\,M_\odot$4-like galaxies, SMBHs with accreted mass above the empirical $3.39\times10^5\,M_\odot$5 relation are reported to be about 15% more effective at removing metals from the ISM than under-massive SMBHs in star-forming galaxies; they suppress overall star formation and more effectively move metals from the ISM into the CGM, while showing little evidence for evacuating gas from their halos (Sanchez et al., 2023).

4. Off-center SMBHs, close pairs, dual AGN, and triple systems

Because Romulus25 explicitly models SMBH orbital decay, it predicts a substantial off-center SMBH population. In a sample of 26 isolated Milky Way–mass halos at $3.39\times10^5\,M_\odot$6, the simulation yields an average of $3.39\times10^5\,M_\odot$7 SMBHs within the virial radius, excluding satellite-halo SMBHs, and $3.39\times10^5\,M_\odot$8 SMBHs within the central 10 kpc, including the central SMBH. Wandering SMBHs are long-lived: 90% of SMBHs now at $3.39\times10^5\,M_\odot$9 kpc first crossed that boundary more than 2 Gyr ago, and in disk-dominated hosts only $2.12\times10^5\,M_\odot$0 of non-central SMBHs within 10 kpc lie within $2.12\times10^5\,M_\odot$1 of the stellar-disk plane, compared with 50% for an isotropic distribution, a deficit significant at $2.12\times10^5\,M_\odot$2 (Tremmel et al., 2018).

Romulus25 also provides the first self-consistent cosmological prediction for close SMBH pair formation timescales after galaxy mergers. The average comoving formation-rate density of close SMBH pairs is quoted as $2.12\times10^5\,M_\odot$3. The pair-formation delay is defined as $2.12\times10^5\,M_\odot$4, where the endpoint corresponds to a bound pair with separation $2.12\times10^5\,M_\odot$5 pc. The summarized distribution gives a median $2.12\times10^5\,M_\odot$6 Gyr and a 75th percentile of $2.12\times10^5\,M_\odot$7 Gyr, while only $2.12\times10^5\,M_\odot$8one in eight galaxy mergers produces a sub-kpc BH pair within a Hubble time. Low stellar-mass-ratio, low-core-density mergers are especially inefficient and are a major source of long-lived wanderers (Tremmel et al., 2017).

The dual-AGN population has been analyzed directly in Romulus25 for $2.12\times10^5\,M_\odot$9. Defining AGN as SMBHs with $6\times10^4\,M_\odot$0 and duals as pairs of such SMBHs within 30 pkpc, the study reports that the number of both single and dual AGN increases from lower to higher redshift. Dual AGN with separations of 0.5–4 kpc are twice as numerous as duals with separations of 4–30 kpc, and all dual AGN in the sample result from major mergers. When followed forward from halo infall, the outcomes are summarized as $6\times10^4\,M_\odot$1 successful SMBH mergers, $6\times10^4\,M_\odot$2 stalled binaries, and $6\times10^4\,M_\odot$3 ejections (Saeedzadeh et al., 2024).

Triple-SMBH dynamics have been explored with high-resolution follow-up $6\times10^4\,M_\odot$4-body integrations initialized from Romulus25 systems. These studies agree that the two most massive SMBHs form an efficiently hardening binary that coalesces within fractions of the Hubble time, while the lightest SMBH is ejected, remains on a wide orbit, forms a stable hierarchical triple, or later binds to the post-merger remnant. In a triaxial follow-up using ROMULUS25 initial conditions, the hardening rates $6\times10^4\,M_\odot$5–$6\times10^4\,M_\odot$6 show no strong systematic variation with the host triaxiality parameter $6\times10^4\,M_\odot$7, and the effect of triaxiality on coalescence time is described as mild, within a factor $6\times10^4\,M_\odot$8 (Koehn et al., 2023, Saha et al., 23 Oct 2025).

5. Gravitational-wave source hosts and host-galaxy identification

Romulus25 has become a key theoretical dataset for host-galaxy studies of low-frequency gravitational-wave sources. For the nano-Hz background relevant to pulsar timing arrays, the simulation is used to argue that the dominant contribution arises from sources with high chirp masses that are likely to reside in low-redshift early-type galaxies with high stellar mass, largely old stellar populations, and low star-formation rate, located at the centers of galaxy groups and showing evidence of recent mergers. The source masses are also stated to correlate with the halo mass and stellar mass of the host galaxies (Saeedzadeh et al., 2023).

Electromagnetic host identification has been developed in two complementary directions. In imaging, mock telescope images were generated for 201 simulated galaxies in Romulus25 hosting recent MBH mergers, and a linear discriminant analysis combining multiple morphological statistics was reported to identify MBH-merger host galaxies with accuracies that increase with chirp mass and mass ratio. For mergers with $6\times10^4\,M_\odot$9 and mass ratio $0.3\,m_{\rm gas}$0, the accuracy reaches $0.3\,m_{\rm gas}$1 and does not decline for at least $0.3\,m_{\rm gas}$2 Gyr after numerical merger. The proposed explanation is that the most distinctive host signature is a prominent classical bulge rather than short-lived tidal disturbance (Bardati et al., 2023).

In integral-field spectroscopy, the same Romulus25 host-galaxy sample was post-processed into synthetic optical IFU datacubes, and stellar kinematic maps were analyzed with another linear discriminant model. For the high-chirp-mass, high-mass-ratio subsample, the reported performance reaches $0.3\,m_{\rm gas}$3 accuracy and $0.3\,m_{\rm gas}$4 precision, with the hosts characterized by lower $0.3\,m_{\rm gas}$5, larger kinematic misalignment $0.3\,m_{\rm gas}$6, and higher mean stellar dispersion. These signatures are described as long-lived for at least $0.3\,m_{\rm gas}$7 Gyr after the numerical merger and are associated with massive early-type galaxies that have undergone major mergers (Bardati et al., 2024).

Taken together, these Romulus25 results define a coherent host-selection picture in which bulge prominence, slow rotation, and stellar kinematic misalignment act as signposts for the most easily classifiable massive-black-hole binaries and mergers. This suggests that, within the assumptions of the mock-observation pipelines, galaxy structure and stellar kinematics can materially reduce the candidate-host set inside a gravitational-wave localization region.

6. Diffuse, low-surface-brightness, and ultra-diffuse galaxies

Romulus25 has also been used extensively to study low-surface-brightness and ultra-diffuse galaxies as products of merger-driven angular-momentum evolution. In isolated field dwarfs, the simulation identifies 134 UDGs among 613 isolated dwarfs with $0.3\,m_{\rm gas}$8, corresponding to 22% of the sample, or “$0.3\,m_{\rm gas}$9” in the abstract formulation. These isolated UDGs are described as products of major mergers that typically occur early, almost never after $\sim700$0, produce a temporary boost in spin, and redistribute star formation to larger radii, suppressing central star formation and steepening negative $\sim700$1 color gradients while maintaining global SFRs (Wright et al., 2020).

The methodological definition of a UDG is itself shown to matter strongly in Romulus25. Using multiple surface-brightness and size cuts across isolated, satellite, and cluster environments, one study concludes that the number of galaxies classified as UDGs can change by up to approximately 45% of the dwarf population depending on the adopted definition. The same work shows that UDG classification depends on viewing orientation and that this dependence decreases as environmental density increases. Under the fiducial definition used there, isolated UDGs are more oblate, or diskier, than non-UDG dwarfs, consistent with a high-spin merger channel; under more permissive definitions, that link is erased (Nest et al., 2021).

A related Romulus25 result concerns classical low-surface-brightness (LSB) galaxies. For late-type central galaxies with $\sim700$2, classical LSBs account for $\sim700$3 of the population. They are reported to form predominantly through major mergers in which the secondary galaxy is co-rotating and aligned with the primary gas disk and/or has above-average orbital angular momentum at infall. The merger product is a high-spin galaxy with star formation spread out and inefficient, producing large HI reservoirs, slow enrichment, and low central surface brightness; the quoted HI depletion time is $\sim700$4 Gyr, about $\sim700$5 longer than for the HSB comparison sample (Wright et al., 28 Jul 2025).

Romulus25 has also been compared with NIHAO and with HI follow-up of SMUDGes UDG candidates. In that context, the simulation is described as producing field UDGs through major-merger-driven high-spin halos, yet yielding present-day HI masses, stellar masses, and SFRs that are qualitatively and quantitatively consistent with both NIHAO and the observed sample. When controlling for $\sim700$6, both Romulus25 and the observations show a positive correlation between gas richness $\sim700$7 and effective radius $\sim700$8, without evidence that gas-rich UDGs and classical dwarfs form distinct populations (Motiwala et al., 20 Feb 2025).

7. Satellite populations, environmental applications, and recognized limitations

Beyond SMBH and diffuse-galaxy science, Romulus25 has been used to place Milky Way analogs in a statistical cosmological context. Depending on the analog definition, one study identifies between 66 and 97 Milky Way analogs in the 25 Mpc box and finds that the number of satellites hosted by an analog increases predominantly with host stellar mass. The satellite quenched fraction also increases with host stellar mass and potentially in higher-density environments; paired systems within 1 Mpc of another Milky Way-mass-or-larger halo may have higher quenched fractions than isolated systems, although the evidence is described as not conclusive (Nest et al., 2023).

These applications are accompanied by explicit caveats. One co-evolution study states that the $\sim700$9 volume lacks rare luminous quasars with $10^6\,M_\odot$00 and the most massive $10^6\,M_\odot$01 SMBHs, while the $10^6\,M_\odot$02–250 pc spatial resolution cannot follow sub-100 pc angular-momentum losses, molecular hydrogen, or cold clumps in detail. The same source notes that the adopted AGN feedback is single-mode thermal and that the simulation does not include an explicit multi-phase ISM or high-temperature metal-line cooling (Ricarte et al., 2019).

At the SMBH-merger level, Romulus25 itself stops at the numerical-pair stage near 700 pc, so binary hardening through three-body scattering, gas torques, gravitational recoil, and triple-BH interactions is not directly modeled in the cosmological run (Tremmel et al., 2017). This is why later triple-system and triaxial-host studies reinitialize Romulus25 systems in higher-resolution direct or hybrid $10^6\,M_\odot$03-body solvers (Koehn et al., 2023, Saha et al., 23 Oct 2025).

A more recent structural critique concerns thin dwarfs. Comparing GAMA, DESI, ALFALFA, and Nearby Galaxy catalogs with TNG50, FIREbox, and Romulus25, one study finds no $10^6\,M_\odot$04 simulated galaxies flatter than $10^6\,M_\odot$05 in Romulus25, in clear contrast with the observational samples. The paper attributes this discrepancy to limitations in resolution and in the treatment of baryonic physics, and presents it as a challenge for current cosmological hydrodynamical simulations of low-mass disk formation (Benavides et al., 11 Dec 2025).

Romulus25 therefore occupies a specific methodological position in the literature: large enough to support demographic and environmental statistics, yet detailed enough to model SMBH seeding and galaxy-scale orbital decay; sufficiently predictive to support multi-messenger host studies, yet still limited at the sub-700-pc binary scale and in some aspects of dwarf-galaxy structure. Its continuing significance lies in this combination of breadth and explicit SMBH dynamics, which has made it a recurrent reference point across galaxy evolution, SMBH demographics, and low-frequency gravitational-wave astrophysics.