Shining Light on the Hosts of the Nano-Hertz Gravitational Wave Sources: A Theoretical Perspective

The formation of supermassive black holes (SMBHs) in the Universe and its role in the properties of the galaxies is one of the open questions in astrophysics and cosmology. Though, traditionally, electromagnetic waves have been instrumental in direct measurements of SMBHs, significantly influencing our comprehension of galaxy formation, gravitational waves (GW) bring an independent avenue to detect numerous binary SMBHs in the observable Universe in the nano-Hertz range using the pulsar timing array observation. This brings a new way to understand the connection between the formation of binary SMBHs and galaxy formation if we can connect theoretical models with multi-messenger observations namely GW data and galaxy surveys. Along these lines, we present here the first paper on this series based on {\sc Romulus25} cosmological simulation on the properties of the host galaxies of SMBHs and propose on how this can be used to connect with observations of nano-Hertz GW signal and galaxy surveys. We show that the most dominant contribution to the background will arise from sources with high chirp masses which are likely to reside in low redshift early-type galaxies with high stellar mass, largely old stellar population, and low star formation rate, and that reside at centers of galaxy groups and manifest evidence of recent mergers. The masses of the sources show a correlation with the halo mass and stellar mass of the host galaxies. This theoretical study will help in understanding the host properties of the GW sources and can help in establishing a connection with observations.


INTRODUCTION
The discovery of Gravitational Waves (GWs) by the LIGO-Virgo-KAGRA (LVK) Collaboration from coalescing compact object binaries of a few tens of solar masses inaugurated the era of gravitationalwave astronomy, enabling the observations of previously inaccessible astrophysical phenomena (Aasi et al. 2015;Abbott et al. 2016a;Acernese et al. 2014Acernese et al. , 2019;;Abbott et al. 2018;Akutsu et al. 2019Akutsu et al. , 2020)).Following this initial discovery, several more binary objects have been detected, one of which (GW170817) also had an electromagnetic (EM) counterpart and stands as the first multi-messenger measurement involving GW signal (Abbott et al. 2016b(Abbott et al. , 2017(Abbott et al. , 2021b,a;,a;Abbott et al. 2023a;The LIGO Scientific Collaboration et al. 2021;Abbott et al. 2023b).With the aid of ongoing and upcoming networks of GW detectors, several more detections of coalescing black hole binaries are likely over the frequency range of 10 Hz and above.
Along with the high-frequency GW signal, coalescing supermassive black holes (SMBHs) can also produce GW signals detectable at low-frequency bands, ranging from a few nano-Hertz to milli-Hertz range.In the milli-Hertz frequency range, upcoming GW detectorssuch as Laser Interferometer Space Antenna (LISA; Amaro-Seoane † E-mail:suvodip@tifr.res.inet al. 2017;Baker et al. 2019) -can probe signal from the SMBHs of masses in the range approximately 10 4 -10 7 M ⊙ .The nano-Hertz GW signal from sources with masses above 10 8 M ⊙ can be detected and characterized using the timing data from several extremely wellstudied millisecond pulsars (Sesana et al. 2008b;Foster & Backer 1990).These signals are the target of the International Pulsar Timing Array (IPTA) collaboration (Antoniadis et al. 2022), comprising the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; McLaughlin 2013), the Indian Pulsar Timing Array Project (InPTA; Joshi et al. 2018) and the Parkes Pulsar Timing Array (PPTA; Manchester et al. 2013).Along with the IPTA Collaboration, the Chinese Pulsar Timing Array (CPTA; Xu et al. 2023a) are also making measurements in this band.In the future with the operation of the Square Kilometer Array (SKA) (Sesana et al. 2008b;Terzian & Lazio 2006;Janssen et al. 2014) more accurate measurement of this signal will be possible (Burke-Spolaor et al. 2019).The recent evidence of the stochastic GW background (SGWB) in the nano-Hertz range by CPTA (Xu et al. 2023b), EPTA+InPTA (Antoniadis et al. 2023a), NANOGrav (Agazie et al. 2023) and PPTA (Reardon et al. 2023) promises to open an exciting new window onto the evolving population of binary supermassive black holes (SMBBHs) in the Universe.
The presence of nano-Hertz GW signal leads to several interesting questions such as: What are the astrophysical properties of the host galaxies of the source SMBBHs, and can one use conventional galaxy surveys to identify (if not uniquely detect) these host galaxies?Since SMBHs reside in the centers of galaxies, SMBBHs are expected to be byproducts of galaxy mergers.Consequently, SMBBHs-host galaxy identifications can potentially shed light on the pathways leading to the formation of the SMBBHs, including their dynamical evolution from the time of first encounter, and more generally on the astrophysics of galaxy mergers.They can also potentially provide insights into the growth of SMBHs and the implications of SMBH-SMBH mergers on galaxy formation (Cattaneo et al. 2009), including conditions leading to the transition between radiative versus kinetic feedback modes (eg Narayan & Quataert 2005;Merloni & Heinz 2008;Benson & Babul 2009;Babul et al. 2013).(See also O'Sullivan et al. 2012, Reynolds et al. 2014, Prasad et al. 2020and O'Sullivan et al. 2021 for rare examples of observed SMBHs in unexpected state.)This would be an important step towards a new paradigm of multi-messenger science capable of addressing a broad spectrum of questions related to astrophysics and cosmology.However, typical PTA localizations in the near term are expected to encompass several thousand (if not more) galaxies.Theoretical and computational modeling offers an opportunity to not only explore environments where SMBH -SMBH mergers take place but also a way to narrow the field of candidate host galaxies for more detailed observational scrutiny (Rosado & Sesana 2014;Goldstein et al. 2019;Volonteri et al. 2003Volonteri et al. , 2020;;Muhamed Kozhikkal et al. 2023).In this study, we use results from a high-resolution Romulus cosmological simulation (Tremmel et al. 2017(Tremmel et al. , 2019) ) to explore this possibility.As we discuss in §3, Romulus suite of simulations is especially suited for investigating SMBH/SMBBH-galaxy connections because of its unique approach to seeding, accretion, and especially the dynamics of supermassive black holes.Consequently, the simulations have previously been used to explore a variety of related topics, including the timescale for the formation of close SMBH pairs following galaxy mergers (Tremmel et al. 2018), the galaxy-SMBH coevolution (Ricarte et al. 2019a), the origin and demographics of wandering black holes (Ricarte et al. 2021), and the demographics of dual active galactic nuclei (Saeedzadeh et al. 2024).
The paper is organized as follows: In §2, we briefly discuss the motivation behind the present study and in §3, we discuss the Romulus simulation.The expected SGWB based on the simulation and the astrophysical properties of the galaxies hosting SGWB sources are discussed in §4 and §5.Among the properties we consider are: gas density ( gas ), star formation rate (  * or SFR), stellar mass (M * ), galaxy morphology, galaxy color, specific star formation rate (sSFR ≡  * / * ) and halo mass (M h , which we take to be M 500 ; see §3.4 for definition of M 500 ).Then, we discuss possible techniques to validate the connection between the SMBBHs and their host galaxies that we find in §6.Finally, we summarise our findings and discuss the future outlook in §7.

MOTIVATION
On one hand, we have the recently detected SGWB in the nano-Hertz range from coalescing SMBHs of mass M > 10 7 M ⊙ .On the other hand, we have spectroscopic/photometric galaxy surveys that are capable of detecting faint galaxies up to high redshifts, some of which will be hosts of the SMBBHs that are contributing to the nano-Hertz SGWB.The combination of these two opens up the prospect of a new multi-messenger science that can shed light on several key questions in astrophysics and cosmology.A limited list of these key questions are: (i) How do the SMBHs grows with time?(ii) How do SMBBHs form and is there a relationship between their formation and one or more properties of the host galaxies?(iii) Do the astrophysical properties of the host galaxies play a role in coalescing of the SMBBHs?(iv) What is the occupation number of the SMBHs in galaxies (or halos) of different masses?
We are interested in understanding the theoretical dimensions of these questions and in identifying whether the key astrophysical properties of the host galaxies can be predicted based on our current understanding of galaxy formation.In this paper, we explore the astrophysical "tells" of galaxies that host SMBBHs in the Romulus simulation volume.We also investigate the properties of the halos of these galaxies.Although the Romulus simulations can track black holes across nearly three orders of magnitude in mass (10 6 -10 9 M ⊙ ), in the present paper we focus primarily on coalescing binaries black holes that can contribute to the stochastic gravitational wave background in the frequency band accessible to PTA.We perform a simulation-based study of the correlations between the SMBBHs and their host galaxies.The specific galaxy properties we focus on include their morphology, star formation rate, galaxy color, stellar mass, gas density, and halo mass.Uncovering a theoretical connection between the properties of the host galaxy and its SMBBH will help motivate observational and data analysis strategies aimed at identifying the host galaxies of the GW sources from the photometric/spectroscopic galaxy catalogs.This, in turn, can contribute to building a data-driven understanding of the evolution of SMBBHs in galaxies.
In future papers in this series, we will consider black holes accessible to LISA, examine possible connections between these and SMBBHs detectable with PTA, and its implementation on the latest nano-Hertz observations (Xu et al. 2023b;Antoniadis et al. 2023b;Agazie et al. 2023;Zic et al. 2023) to identity the possible host candidates.For completeness, we note that there are several analytical and numerical simulation-based studies estimating SGWB signal in the PTA frequency range (Rajagopal & Romani 1995;Jaffe & Backer 2003;Kocsis & Sesana 2011;Chen et al. 2017;DeGraf et al. 2021;Izquierdo-Villalba et al. 2022, 2023).Additional studies are referenced throughout the text.

THE ROMULUS SIMULATIONS
In this work, we present results from the analysis of the Romulus25 simulation, which is a (25 cMpc) 3 cosmological volume simulation from the Romulus suite (Tremmel et al. 2017(Tremmel et al. , 2019;;Butsky et al. 2019;Jung et al. 2022;Saeedzadeh et al. 2023).
The simulation was run using the Tree+Smoothed Particle Hydrodynamics (Tree+SPH) code CHaNGa (Menon et al. 2015;Wadsley et al. 2017).CHaNGa incorporates the standard physics models previously employed in the simulation codes GASOLINE/GASOLINE2 and has been extensively tested (Stinson et al. 2006).Including modules for star formation, stellar feedback, turbulent diffusion (Shen et al. 2010), the UV background, low-temperature metal cooling, and an improved treatment of both weak and strong shocks.However, the models for SMBH formation, dynamics, growth, and feedback are novel (Tremmel et al. 2015(Tremmel et al. , 2017)).We will discuss them in more detail in the following subsections.
The Romulus25 simulation was run assuming a flat ΛCDM universe with cosmological parameters consistent with the Planck 2016 results (Ade et al. 2016): Ω m = 0.309, Ω Λ = 0.691, Ω b = 0.0486, H 0 = 67.8km s −1 Mpc −1 , and  8 = 0.82.The simulation has a Plummer equivalent gravitational force softening of 250 pc (or 350 pc spline kernel) and a maximum SPH resolution of 70 pc.Differing from many similar cosmological runs, the dark matter distribution in Romulus25 is oversampled with 3.375 times more dark matter particles than gas particles.This results in a dark matter particle mass of 3.39 × 10 5  ⊙ and a gas particle mass of 2.12 × 10 5  ⊙ .This deviation from the standard approach of simulating equal numbers of gas and dark matter particles minimizes numerical noise and allows for more precise black hole dynamics (Tremmel et al. 2015).
More details about the Romulus25 simulation have been described in a number of published papers.We refer interested readers to Tremmel et al. (2015Tremmel et al. ( , 2017Tremmel et al. ( , 2019Tremmel et al. ( , 2020)); Sanchez et al. (2019); Butsky et al. (2019); Chadayammuri et al. (2021); Jung et al. (2022);and Saeedzadeh et al. (2023).The latter two especially offer a concise yet complete summary.Additionally, for comparisons with galaxy formation models employed in other cosmological simulations, interested readers are referred to Somerville & Davé (2015); Vogelsberger et al. (2020); Oppenheimer et al. (2021) There are, however, a few aspects of the Romulus25 simulation that are important to highlight as these are relevant to the present discussion.These pertain to the treatment of SMBH seeding, growth, and dynamical evolution in the Romulus25 (Tremmel et al. 2017).

SMBH Seeding
Unlike many other cosmological simulations (e.g., Schaye et al. 2015;Weinberger et al. 2017;Pillepich et al. 2018;Nelson et al. 2019;Davé et al. 2019), the Romulus25 SMBH seed model does not depend on a halo or a galaxy to exceed a certain mass threshold for a SMBH to form.Rather, the seeding depends only on the local gas properties (Tremmel et al. 2017).As a result, the SMBHs in Romulus25 can form in low mass halos and tend to form much earlier (z > 5, Tremmel et al. 2017).Additionally, one can also have multiple SMBHs arising in the same halo.
The criteria for converting gas particle into a SMBH seed in Ro-mulus25 are as follows: (i) The gas particle must be both eligible and selected to form a star.The latter is a probabilistic process.(ii) The gas particle must have very low metallicity ( < 3 × 10 −4 ); (iii) its density must be very high; i.e. at least 3   /cc or greater.And, (iv) its temperature is within the range of 9500 -10000 K.This seeding prescription resembles the direct collapse black hole scenario, where high temperatures and low metallicities suppress fragmentation and allow sizeable gas clouds to collapse directly into an SMBH seed (Lodato & Natarajan 2007;Alexander & Natarajan 2014;Natarajan 2021).
In the Romulus25 simulation, SMBHs are seeded with a mass of 10 6 M ⊙ .This seeding mass differs slightly from other simulations like TNG50/100/300, EAGLE, Horizon AGN, and SIMBA-C where the initial SMBH seed masses are set at ∼ 8 × 10 5 M ⊙ , 10 8 M ⊙ , 10 5 M ⊙ and 10 4 M ⊙ respectively (Nelson et al. 2019;Kaviraj et al. 2017;Crain et al. 2015;Hough et al. 2023).The choice of the SMBH seed mass in Romulus25 is constrained primarily by two factors: (i) the resolution of the simulation, with the dark matter and gas mass resolutions in Romulus25 being 3.39 × 10 5  ⊙ and 2.12 × 10 5  ⊙ respectively, and (ii) the necessity to keep SMBHs more massive than dark matter and star particles.This latter requirement is crucial for reducing the occurrence of spurious scattering events (Tremmel et al. 2015).
The resulting SMBH occupation fraction at z=0 is consistent with current observations even on the scale of dwarf galaxies (Ricarte et al. 2019b).Additionally, the SMBH masses correlate with the stellar masses of their host galaxies, following the observed SMBH mass -stellar mass relation (Tremmel et al. 2017;Ricarte et al. 2019c).

SMBH Dynamics and Mergers
The other difference in the SMBH model between Romulus25 and other cosmological simulations is SMBH dynamics.Unlike many other cosmological simulations, where the SMBHs are artificially placed at the gravitational potential minimum of their host galaxies (e.g.Crain et al. 2009;Sĳacki et al. 2015;Davé et al. 2019), the Romulus25 simulation accurately tracks the dynamical evolution of SMBHs down to sub-kpc scales, which is highly advantageous for the present study.To achieve this, a sub-grid correction is employed that accounts for the unresolved dynamical friction from stars and dark matter that the SMBHs ought to be experiencing (Tremmel et al. 2015).For each SMBH in the simulation, this force is estimated by assuming a locally isotropic velocity distribution and integrating Chandrasekhar's equation (Chandrasekhar 1943) from the 90-degree deflection radius (r 90 ) to the SMBH's gravitational softening length (  ).The resulting acceleration is In order for two SMBHs to merge, they must be within a distance of two gravitational softening lengths (0.7 kpc) and possess a low enough relative velocity to be mutually bound; i.e. 1 2 Δv 2 < Δa • Δr, where Δv and Δa are the differences in velocity and acceleration of the two black holes, and Δr is the distance between them (Bellovary et al. 2011;Tremmel et al. 2017) 1 .The separation limit of two gravitational softening lengths is deemed appropriate because once the separation drops below this limit, the simulation's ability to accurately track the SMBH pair's dynamics becomes less reliable.
When a merger takes place, the resulting SMBH is assigned a velocity that conserves momentum, and its mass is the sum of the masses of its progenitors.Mergers are one of the two processes driving the growth of SMBHs.

SMBH Growth and Feedback
The other process by which SMBHs grow is via the accretion of gas.In Romulus25, this accretion rate is estimated via a modified Bondi-Hoyle (Bondi 1952, for modifications see Tremmel et al. 2017) prescription applied to the smoothed properties of the 32 nearest gas particles: where  gas is the ambient gas density,   is the ambient sound speed,   is the local rotational velocity of surrounding gas, and   is the bulk velocity relative to the SMBH.All ambient quantities are calculated using the 32 nearest gas particles.The introduction of   and   terms in the above aims to remedy the neglect of gas bulk motion and angular momentum in the original Bondi-Hoyle formulation.Finally, the coefficient  is introduced to correct for the suppression of the black hole accretion rate due to resolution effects.It is defined as where   ℎ, * is the star formation number density threshold (0.2   /).
Gas accretion onto a SMBH results in energy release into the environment around the black hole.In Romulus25, it is assumed that this energy is electromagnetic and that a fraction of it will couple to the ambient gas and contribute to its internal energy.The thermal energy deposition rate is given by  •, ℎ =      •  2 , where   is the radiative efficiency (assumed to be 10%) and   is gas coupling efficiency (set to 2%).The thermal energy is imparted isotropically to the 32 nearest gas particles, with the energy being distributed among these gas particles according to the smoothing kernel.We refer readers to Tremmel et al. (2017) for further details.
The halos and subhalos exist in a nested hierarchy, where the halos are the primary structures and the subhalos are incorporated within them.To identify these structures, AHF first locates density peaks in an adaptively smoothed density field and identifies all the particles (dark matter, gas, stars, and black holes) that are gravitationally bound to these peak.This process is repeated on successively larger scales until all the structures in the hierarchy have been found.Once the halos are identified, their centers are found by applying the shrinking sphere approach (Power et al. 2003) to the distribution of bound particles associated with each of the halos.
In the case of subhaloes, AHF tracks the local density profile from the peak center outward.At some point, the external gravitational field starts to dominate, altering the shape of the density profile.The distance from the peak to where this happens is taken to be the size of the subhalo, and the mass enclosed is recorded as the subhalo's mass.
We also track all the SMBHs in the Romulus25 simulation volume.We use the resulting information to construct merger trees for all the black holes.At each redshift, we then identify black holes that have experienced a merger during the immediately preceeding output and flag the about-to-merge SMBH pairs as candidate sources of nano-Hertz SGWB.The time resolution (Δ) and redshift resolutions (Δ) for the saved output files within our redshifts of interest (see §5.1) vary in the ranges of 10 Myr < Δt < 400 Myr and 0.002 < Δz < 0.1, getting to smaller values as approaching z = 0.The typical separation of merging SMBH pairs is ∼ 1 kpc and their maximum separation is 2.8 kpc.For completeness, we also identify all black hole pairs separated by ≤ 1.4 kpc and which are not flagged as merging in the next output.We will refer to these as proximate pairs.
We emphasize that the SGWB from flagged SMBBHs with separation scale of ∼ 1 kpc cannot contribute to the nano-Hertz frequency band unless they coalescence down to sub-parsec (10 −5 pc) scale.This journey of the SMBBHs from the scale of ∼ 1 kpc to ≤ 10 −5 pc is governed not only by GW emission but also by environmental effects such as dynamical friction, stellar loss cone and viscous gas drag.These processes are not resolvable in Romulus25 or, for that matter, in any other cosmological simulation.We therefore need to model this coalescence separately.

Modeling SGWB signal from coalescing SMBBHs
In order to calculate the contribution to the SGWB signal from the coalescing SMBBHs, we start with the expression for the characteristic strain of the GW signal ℎ  at frequency  for a source emitting at a rest-frame frequency   = (1 + )  (Rajagopal & Romani 1995;Phinney 2001;Sesana et al. 2008a): where the distribution function, , is the number density of SMBBH GW sources with black hole masses in the range and determines the amplitude and spectral shape of the SGWB signal.The second term, , quantifies the amount of GW energy released per logarithmic rest-frame frequency by a binary of source masses  1 and  2 at redshift .The latter is the product of the GW energy emission rate ( ), and the residence time (i.e. the amount of time a source spends at a frequency:    ln ).Following Kelley et al. (2017a,b), we write the energy released as where for circular orbits emitting signals up to the innermost circular stable orbit (ISCO).Here,   = ( 1  2 ) 3/5 /( 1 +  2 ) 1/5 is the binary's chirp mass, and  is the frequency, which at ISCO is given by   ,ISCO =  3 /(6 3/2   tot ) in terms of total mass of the binary  tot =  1 +  2 .In the presence of higher harmonics, this equation modifies to a sum over the higher harmonics (Enoki & Nagashima 2007).
As for the second term in Eq. ( 5), the ratio  ℎ  GW (  ) captures the residence time of the GW signal at a particular frequency.The numerator ( ℎ ≡  /  ) is the binary hardening time expressed in terms of the semi-major axis of the binary .Initially, this timescale depends on the environmental effects arising due to the interaction between the binaries and their local environment.These effects include (i) dynamical friction, (ii) stellar loss-cone scattering, and (iii) viscous drag.The impact of these environmental effects is among the major sources of uncertainty in the spectral shape of the signal but typically these environmental effects reduce the residence time of the GW signal at a particular frequency and the ratio will be less than one.
As we have noted, the above environmental effects cannot be directly computed from the Romulus25 simulation.Moreover, from an EM observations point of view, resolving galaxies on sub-parsec scales at a cosmological distance is not possible with currently ongoing and upcoming surveys.However, we can determine the average astrophysical properties of a galaxy -like gas density, stellar mass, halo mass, and other properties -on kpc scales from cosmological simulations as well as observations.We therefore model the ratio,  ℎ  GW , in terms of the average astrophysical properties of host SMBBH galaxies: In effect, we want to construct a framework that can relate the nano-Hertz (nHz) GW signal detectable from PTA with the observable quantities of galaxies.

Modelling the environmental effect
In this subsection, we discuss the model for E (  ,  * ,  * ,  ℎ ,  gas , ) in greater detail.But first we note that the impact of the environmental effects is greatest when the binaries are further away from each other and are radiating at lower frequencies of GW (Volonteri et al. 2003(Volonteri et al. , 2020;;Kocsis & Sesana 2011;Sampson et al. 2015;Chen et al. 2017;Kelley et al. 2017a,b).As the SMBBHs inspiral and their separation decreases, they emits GW signals at increasingly higher frequencies.At frequencies of around 1 yr −1 (or about a few ×10 −8 Hz), the environmental effects are no longer dominant.The SMBBHs' evolution is dominantly through GW emission, and the frequency-dependent part of the ratio  ℎ  GW (  ) approaches unity.Impact of these effects on the GW strain are often modelled using parametric forms (Sampson et al. 2015;Chen et al. 2017).
In the present case, Romulus25 allows us to track the evolution of the black holes down to a few hundred parsecs but none of the current generation of cosmological simulations have the resolution to follow their evolution due to the above processes to smaller scales.Moreover, one also needs very high-resolution observations to determine the density profile of stars and gas at these scales.We therefore use a parametric equation to capture the environmental effect E. In effect, E can be thought of as a subgrid model that uses accessible galaxy properties to estimate the overall number of SMBBH sources that will contribute to the nHz signal as well as the amount of their orbital energy that goes into GWs.
As noted in the last section, one of the key aspects of the problem where the astrophysical properties of the galaxies play an important role is in determining the fraction of SMBBHs that can contribute to the SGWB in the nHz range.These SMBBHs can successfully reach from the orbital separation ∼ kpc scales to about 10 −5 pc (i.e. the GW emission-dominated regime) from within the age of the Universe.We model this via a dimensionless parameter , which quantifies how efficiently the SMBBHs identified in the simulations on a ∼ kpc scale will overcome the last parsec problem.We expect that  will depend on the various astrophysical properties of the host galaxy and given the specific nature of the processes involved, we make an ansatz that  will primarily depend on the galaxy's gas density ( gas ), stellar mass ( * ), and star formation rate (  * ), and can be written as  =   log( gas = 10 7  ⊙ /kpc 3 ) log( gas ) +   * log( * = 10 10  ⊙ ) log( * ) where   ,   * , and   * govern how gas density, stellar mass and star formation rate, respectively, are in driving the hardening of the black hole binaries.
The other aspect, which plays a crucial role in controlling the shape of the stochastic GW power spectrum, is the amount of orbital energy that is lost via environmental processes.This energy loss is modeled by the factor (1 + (    ) − ) − .Here the dimensionless factor  captures the frequency-dependent loss of GW signal due to processes like dynamical friction, stellar hardening, and viscous drag, relative to the case where these environmental effects are absent and the hardening of the binary is driven only by GW emission.We model this as The term  controls the spectral behavior of the environmental effects (Sampson et al. 2015).The value for some of the effects such as stellar scattering is 10/3.However, the combination of various effects can lead to a different spectral index.Finally, the parameter  controls the overall tilt of the environmental effects.For a fiducial case with a GW-emission-only scenario, the value of  = 1 can be considered as a fiducial.However, there can be deviations due to astrophysical effects.For the spectral shape of the signal, we use three parameters, namely , , and   where   is the transition frequency at which the GW dominant effects become important over the environmental effects.The transition frequency can be expressed in terms of the stellar density  * (in units of M ⊙ pc −3 ) and velocity dispersion  * (in units of km/s), eccentricity , and chirp mass   (in units of solar mass) as where  () = (1 + (72/24) 2 + (37/96) 4 )/((1 −  2 ) 7/2 ). 0 is a correction factor incorporating any effect that may not be captured by this simplistic approximate formula (such as mass ratio).For  * = 100 ⊙ pc −3 ,  * = 200 km/s,   = 10 9 M ⊙ ,  0 = 1, the value of   is around 0.4 nHz (Chen et al. 2017).Combining all the together, we can write the total contribution from the host galaxy properties as where ,  are defined above in Eq. ( 8) and Eq. ( 9).

SGWB estimation from SMBBHs in Romulus simulation
Having put in place the above elements, we can use Eq. ( 4), Eq. ( 5), and Eq. ( 7) to estimate the SGWB signal from discrete sources in a cosmological simulation as follows ℎ where the sum runs over all the source SMBHHs (or equivalently, host galaxies in which coalescing SMBHs are present) in a simulation box of comoving volume   = (25Mpc) 3 that contributes to GW background.
We use the GW sources identified in the Romulus25 simulation (see §5.1 for details about how these sources are identified) and the above relationships to model the SGWB.As described in §4.2, the environmental effects -the  ℎ  GW (  ) part of Eq.12 -are modelled by E via the  and  parameters defined in Eq. ( 8) and Eq. ( 9).For   ,  * , and  * values that enter into these equations, we use the median values of these properties from the central 1 kpc region of the host galaxies of the Romulus 25 sources.For additional details, see §5.1.As for our procedure, after inputting the median values of   ,  * , and  * into Eq.( 8) and Eq. ( 9), the free parameters of   ,  * ,   * ,   , and   * are chosen such that that the amplitude of the signal matches with the observed nHz signal at frequency  = 1 −1 (Xu et al. 2023b;Antoniadis et al. 2023a;Agazie et al. 2023;Zic et al. 2023).This results in  = 0.2 and  = 1.The results of modeling the SGWB are shown in Fig. 1.For SMBHs with chirp mass   ≥ 10 8 M ⊙ , the results for the case where environmental effects are not considered is shown as solid teal line and denoted as "reference".The dotted and dashed lines show the results for SMBBHs with chirp mass   ≥ 10 8 M ⊙ and   ≤ 10 8 M ⊙ , respectively, when environmental effects are taken into account.The salmon dotted line corresponds to   ≥ 10 8 M ⊙ and parameters  = 0.2,  = 1,   = 5 nHz,  = 10/3 and  = 1.The values of  and  are computed as described above and the values of   ,  and , which control the shape of the signal, are fiducial values from physical models discussed in the previous subsection (cf §4.2).
For the same parameters, the result for SMBBHs with   ≤ 10 8 M ⊙ is shown by a yellow dashed line.An important point to note is that the contribution to the total signal from SMBBHs with chirp mass   < 10 8 M ⊙ is not significant; they only contribute about 10% of the signal arising from sources with   ≥ 10 8 M ⊙ .
To illustrate the impact of variations in environmental properties, we present results for SMBBHs with   ≥ 10 8 M ⊙ where we vary the values of , ,   ,  and .Specifically, we examine the effects of changing   from 5 to 8 nHz (shown in purple),  from 10/3 to 7/3 (in magenta),  from 1.0 to 0.8 (in brown),  from 1. to 0.5 (in cyan), and  from 0.2 to 0.3 (in indigo).As each parameter value is changed, the spectral shape of the signal exhibits only moderate changes.We will explore the parameter estimation using the galaxy catalog in future work (in preparation).If the properties of the underlying host galaxy can be inferred, then the values of the parameters that control the environmental effects can be measured.

Connecting SGWB signal with galaxy properties
In the case of Romulus25 simulation, we have firsthand knowledge of , the number density of GW sources.However, we aim to find a map between the GW sources and the EM observations -specifically, the galaxies in a complete galaxy catalog.
To connect the EM and GW observational sectors, we assert that the total number of GW sources should equal the total number of GW source host galaxies The term , coming from  ℎ  GW (  ) (cf Eq. 11), and the occupation fraction  appear in the above as a multiplicative factor ( × ), which now takes care of the overall occupation of the number of the SMBBHs contributing to the SGWB in terms of both black hole masses and also the astrophysical properties of the galaxies.One can compare this with the observations and make an inference of this quantity.
Finally, we can also write the GW source distribution function, , in terms of the halo mass function where ) is the SMBBHs occupation number density in a halo of mass M h and is the halo mass function, i.e. the number density of halos in halo mass bin M h (We remind the reader that in the present study, we identify M h with M 500 .).
From a simulation, we can estimate the population of SMBBHs which can contribute to the SGWB in the PTA frequency range, and also identify the mass and redshift of the host halo.This gives us a connection between the SMBHs and halo mass  ℎ written in Eq. ( 15).Similarly we can identify the astrophysical properties of the host galaxies of the coalescing binaries from simulations.These would also be accessible from EM observations.This gives us an avenue to connect the astrophysical properties of the host galaxies with the SMBH properties.

PROPERTIES OF THE HOST GALAXIES OF THE NANO-HERTZ GW SOURCES
The dynamics of the SMBBHs and their contribution to the nHz frequency depends on the local astrophysical properties as discussed in the previous section.However, to identify the key astrophysical properties of the host galaxies that can be identified from observations, we need to explore the large-scale properties of the host galaxy.
In this section, we explore both local astrophysical properties in the vicinity of SMBBHs and their host galaxy properties.The properties we consider include gas density, stellar mass, and star formation rate.While discussing the host galaxies, we also comment on the properties of the halos in which these galaxies reside.

Local galactic properties
In this section, we focus on the characteristics of the gas and stars in the vicinity of the SMBBHs, and specifically, within a 1 kpc radius around the more massive black hole in SMBBHs, as the most massive BH drives the chirp mass and hence the strength of the signal.We also restrict ourselves to SMBBHs with a chirp mass   ≥ 10 8 M ⊙ 2 and  ≤ 2. These SMBBHs account for ∼ 90% of the total SGWB power spectrum (see Fig. 1).Our analysis of the Romulus25 simulation reveals six nHz GW sources.Considering the redshift range in which these sources are detected and our simulation volume, we calculate the number density of nHz GW sources to be 7.7 × 10 −6 cMpc −3 .This value closely aligns with the number densities estimated in Mingarelli et al. (2017) and Casey-Clyde et al. (2022) and that derived by Antoniadis et al. (2023c) from PTA data release, which are 1.6 × 10 −6 cMpc −3 , 6.6 × 10 −6 cMpc −3 and 1.5 × 10 −5 cMpc −3 respectively.
In Fig. 1 we have shown the GW background for various different cases based on the median value of the gas density, star formation rate and stellar mass within a 1 kpc radius around the more massive black hole in nHz GW SMBBH sources in the Romulus25 simulation.Fig. 2 shows, from top to bottom, the gas density ( gas ), star formation rate (SFR), and stellar mass (M * ), vs redshift, of the individual host galaxies of the nHz GW SMBBH sources.No significant evolution with redshift is found for the local SFR and the local stellar mass.However, the top panel shows an increasing trend of  gas with redshift.Nonetheless, examining the properties around black holes solely in nHz GW sources doesn't provide a comprehensive view.
2   ≥ 10 8 M ⊙ implies that the most massive SMBH in the pair is at least 8.7 × 10 7 M ⊙ if both SMBHs are equal mass, or more realistically > 2.5 × 10 8 M ⊙ since the BH mass ratios are typically < 0.2.Therefore, we compare these to the same properties surrounding the most massive of the two black holes in our full set of proximate and merging (cf §3.4) pairs of black holes in the Romulus25 simulation.The resulting distributions for gas density, SFR, and stellar mass are shown as salmon histograms in Fig. 3.The corresponding quantities for nHz GW sources are denoted by a vertical line.We find that the local  gas and SFR around our nHz GW sources are typical.Specifically, the increasing gas density with redshift around the nHz GW sources simply reflects the fact that all systems have more gas at earlier epochs.However, the stellar mass in the vicinity of our nHz GW sources falls in the high mass tail of the corresponding distribution, as illustrated in Fig. 3.This indicates that nHz GW sources are more likely to be found in environments with high local stellar mass or equivalently, stellar density.Though it is important to note that this conclusion is based on only 6 events detectable in a simulation box of (25 cMpc) 3 .A study from a bigger box simulation will help in better understanding the statistical properties.

Global galactic and host halo properties
In this section, we consider the global galactic characteristics of the hosts of nHz GW sources, delving into the observable properties of these host galaxies.We also examine the properties of the halos of these host galaxies.
Figure 4 illustrates the relation between galaxies' gas density ( gas ), stellar mass (M * ), star formation rate (SFR), and specific star formation rate (sSFR) vs redshift.Each of these quantities is calculated inside a 25 kpc sphere around the galaxy center.The first panel shows the galaxy gas density.No significant evolution with redshift is observed.On the other hand, the second and third panels clearly show a decrease in the stellar mass and an increase in SFR with redshift, respectively.Not surprisingly, this results in a rising sSFR with redshift, as shown in the lower right panel.Additionally, the stellar mass of the host galaxies suggests that at the redshifts of interest, these galaxies are among the most massive systems in the Romulus volume.Based on analyses of Jung et al. (2022) and Saeedzadeh et al. (2023), we guess that these are the massive central galaxies in group-scale halos.We will discuss this elaborately at a later part in this section.The specific star formation rate (sSFR) is frequently used to classify galaxies as star forming or quenched.We show the evolution of the sSFR with redshift for the host of the SMBBHs with   ≥ 10 8 M ⊙ in the last panel in Fig. 4.
In this paper, we adopt the criteria from Genel et al. (2018), where they label a galaxy as "main sequence" if its sSFR is within ±0.5 dex of the main sequence ridge and as "quenched" if its sSFR is 1 dex below the ridge.Genel et al. (2018) give a relationship between sSFR and  * at few redshifts.We linearly interpolate across these to determine the main sequence ridge sSFR at redshifts and stellar masses of our host galaxies.In Fig. 5, we show Δlog(sSFR) ≡ log(sSFR galaxy ) − log(sSFR ridge ) for our host galaxies.The shaded region shows the main sequence band and the dashed line corresponds to the threshold below which galaxies are classified as quenched.Five of our six hosts either lie in the quenched territory or are on the border.One galaxy, however, the  = 1.061 host, falls within the star-forming main sequence band.
Next, we examine the rest-frame U-V colors of the host galaxies.For completeness, we show (U-V) rest−frame versus rest-frame V-band luminosity (L V ), rest-frame absolute magnitude (M V ), and stellar mass (M * ).In these plots, large circles represent the host galaxies of nHz GW sources with   ≥ 10 8 M ⊙ , while small circles delineate hosts of nHz GW sources with chirp masses ranging between 10 7  ⊙ <   < 10 8  ⊙ .Generally, the latter tend to be bluer and less massive than the hosts of our nHz GW sources.The four  < 1 host galaxies are consistent with the quenched, early-type galaxies in Coma and Virgo clusters (Renzini 2006) as well as the SDSS and the CANDELS Multi-Cycle Treasury Survey (Bell et al. 2012).This aligns with the trends identified by Izquierdo-Villalba et al. (2023), who shows that elliptical galaxies are significant hosts of massive binary black holes at  < 1.The two  > 1 hosts, one of which is star-forming and the other quenched according to the Genel et al. ( 2018) criterion, are both just slightly bluer than the population of quenched early-types at comparable redshifts in the CANDELS Multi-Cycle Treasury Survey (Bell et al. 2012).The comparison to galaxies at the same redshifts is important since the stellar population in higher redshift galaxies is inherently younger and therefore bluer.
To gain better insight into why the two higher redshift exhibit bluer colors and to further clarify the nature of the host galaxies generally, we examine the images of the host galaxies.These are shown in Fig. 7.The top row shows the edge-on view of the galaxies and the bottom row shows the face-on view.The stars are colored based on their magnitudes, determined by their age and metallicity.The magnitudes from the 'i' band influence the red component of the image, 'v' the green, and 'u' the blue.These channels are then   2018).Here we show Δlog(sSFR) = log(sSFR galaxy ) − log(sSFR ridge ) for these host galaxies.The dashed line, which marks 1 dex below the ridge for each redshift, serves as the quenched threshold.The shaded region indicates the "main sequence".See text for detailed definitions.combined to produce a multiband composite image of the galaxy. 3isually, all of the galaxies appear to be early types, and the majority of the stars in these galaxies are old in the sense that 50% had formed within the first 2.5-2.6 Gyrs after the Big Bang.This is illustrated by the histograms in Fig. 8, which show the cosmic age (i.e. the age of the Universe) when the stars in the galaxies formed.Turning back to Fig. 7, we see that nearly all of the galaxies appear to have experienced recent merger(s).They all manifest features like stellar streams and shells (eg Fardal et al. 2007, and references therein).Some of the mergers are gas rich, as in the case of the host galaxy at  = 1.06 where one can clearly see an extended stream with ongoing star formation against a background of older stars.A less prominent stellar stream can also be seen in the host galaxy at z = 1.671.A detailed analysis of the morphology of host galaxy of merging binary black holes by Bardati et al. (2023) shows that the dominant morphological signature of SMBH mergers is the presence of a classical bulge that is also a sign of major mergers of these host galaxies.The panels in Fig. 8 provide additional information, including the fraction of stellar mass at the time of observation that has an age ≤ 1 Gyr.The latter varies from < 1% to as much as 12% in the highest redshift system.This small fraction of very young stars is the likely explanation for the two  > 1 galaxies' bluer colors.As shown by Pipino et al. (2009), even a small fraction (∼ 1%) of young stars (∼ 0.1 Gyr) can have a dramatic impact on UV-optical colors.
For the SMBBHs, we show the time taken by these sources to grow to masses above   = 10 7 M ⊙ from their progenitor mass of 50% of the source masses in Fig. 9 for both the SMBHs in a SMBBH.The distribution of the delay time for the black holes with M  ≥ 10 7 M ⊙ and M  ≥ 10 8 M ⊙ is shown in blue and orange respectively in Figure 9.The shortest delay time observed for these black holes is 1.4 billion years, roughly 10% of the age of the Universe for sources with M  ≥ 10 8 M ⊙ , which contribute significantly to the SGWB.This indicates that the number of mergers of the PTA sources is likely to be more towards low redshift than high redshift, and the  corresponding properties of the host galaxies will be towards older galaxies.
In Fig. 10, left panel, we show the correlation between the chirp mass of SMBBHs with the stellar mass of their host galaxy.The right panel shows the mass of more massive black holes in the SMBBHs vs. the host galaxies' stellar mass.Sources with black holes mass ratio  < 0.14 are shown by stars and ones with  > 0.1 are shown by circles.Large symbols correspond to SMBBHs with chirp mass   ≥ 10 8 M ⊙ and small symbols denote those with 10 7 M ⊙ ≤   ≤ 10 8 M ⊙ .The plot indicates that sources with the highest chirp mass are primarily present at low redshift and are hosted in galaxies with stellar mass greater than 10 11 M ⊙ .Furthermore, the heaviest black holes are found in pairs with chirp masses   ≥ 10 8 M ⊙ , which are hosted by massive galaxies.These observations support the hypothesis we presented at the beginning of this section: that nHz GW sources predominantly inhabit massive, group-central galaxies.We will explore this further below.
In Fig. 11, we compare the global properties (i.e.global stellar mass and sSFR) and the environment (i.e. the halo mass and location therein) of the host galaxies of nHz GW sources with   > 10 8 M ⊙ (vertical lines) against the properties of all galaxies in the simulations within halos with M 200 > 4.5 × 10 9 M ⊙5 (histograms).The plot supports the discussion we have presented above in the context of the colors and quenched/star-forming status of the host galaxies, that the sSFR of these galaxies places them in the low sSFR tail of the sSFR distribution of the galaxies in the Romulus25 simulation volume.From Fig. 11 we deduce that our host galaxies are among the most massive galaxies in the Romulus25 simulation volume.That the host galaxies of the nHz GW sources reside in galaxies with high stellar mass and low sSFR making them a unique class of objects.This is further confirmed by comparing the host halo M 500 to all6 halos M 500 in the Romulus simulation which is shown in the histogram, we note that the host halo masses are in the high-mass tail of the distribution (see the last two panels in Fig. 11).The halos in which the host galaxies reside are group scale systems and based on the results of Jung et al. (2022) and Saeedzadeh et al. (2023), we expected  -and have subsequently confirmed -that these galaxies are massive central group galaxies.Collectively these findings strongly suggest that the nHz GW sources are hosted by massive early-type galaxies at the centers of groups and clusters.However, we assert that the typical hosts of the nHz GW sources will be group-central galaxies.For one, there are many more groups than clusters.Moreover, the lower velocity dispersion of group satellites makes dynamical friction in group halos more efficient, and consequently, group environments are much more conducive to mergers, especially between the satellite and the central galaxies (O'Sullivan et al. 2017;Oppenheimer et al. 2021, and references therein).
For completeness, in Fig. 12, we present the mass of the more massive black hole in nHz GW sources and the host halo mass of these sources in the top and bottom panels, respectively.As expected, these plots indicate that as redshifts increase, the halo mass decreases, and the black hole mass in the pair also reduces.
In this section, we have elucidated the unique astrophysical properties of the host galaxies of the SMBBHs which contribute to the SGWB signal in comparison to all galaxies in the simulation.Our findings highlight that GW source hosts (with chirp mass   ≥ 10 8 M ⊙ ) predominantly reside in galaxies characterized by lower star formation, higher stellar mass, and higher halo masses compared to most counterparts at a given redshift.Specifically, these hosts are located within group-scale halo systems, identifying as massive central group galaxies.These host galaxies are early-type galaxies, displaying a distinct trend in the color-magnitude diagram across redshifts.These astrophysical properties inferred theoretically about the SMBBHs make it possible to correlate electromagnetic observations of the galaxies with the GW sources.Exploring such connections, coupled with comparisons to theoretical models, offers insights into the interplay between galaxy formation and black hole formation.

POSSIBLE TECHNIQUES TO CONNECT OBSERVATIONS WITH THEORETICAL MODELS
In the last two sections, we have discussed a scheme to connect the global astrophysical properties of the galaxies with the spectrum of the SGWB signal in the PTA band and have applied that to the Romulus simulation to understand the underlying theoretical correlation.The next interesting step forward is to connect this with the observations available from currently ongoing/upcoming surveys (see Burke-Spolaor et al. (2019) for review article).The observation of the GW signal from PTA observations in the nHz range can happen as (i) SGWB and (ii) GW signal from individual events.Both of these kind of observations can bring complementary information.SGWB: The measurement using PTA observations provides a measurement of the spectrum of the SGWB signal.However, it is still unclear what are the properties of the host galaxies of the SMBBHs that contribute to the signal.As we have shown in the previous section, the simulations show that the binaries are likely to form in galaxies with high stellar mass, high halo mass, and low SFR, and mostly early-type galaxies, that show signs of mergers in not too distant a past.We also showed that the host galaxies are central group galaxies.The host of the GW sources also shows a trend in the color-magnitude diagram as a function of redshifts.
Based on these understandings, we can classify galaxies from electromagnetic observations based on their color, stellar mass, halo-mass, SFR, and galaxy type and can explore spatial crosscorrelation of the galaxy distribution with the anisotropic SGWB signal (Mingarelli et al. 2013;Hotinli et al. 2019; Sato-Polito & Kamionkowski 2023) and explore cross-correlation between the two quantities (Mukherjee & Silk 2020;Yang et al. 2020;Mukherjee & Silk 2021;Yang et al. 2023).A detailed paper on this formalism will be followed up in a companion paper.The cross-correlation of the SGWB signal with the galaxies of different types will be maximum for types of galaxies that are host of the GW sources.The exploration of the cross-correlation signal will give us an understanding of the population of the GW sources contributing to the background and we can estimate the occupation number of SMBHs.This will be useful in understanding the SGWB measurement in terms of the astrophysical properties of galaxies given in Eq. ( 14) based on observations.In future work, we will explore this aspect from the measurement of the SGWB signal and galaxies detected from optical and infrared surveys.
Signal from individual events: The measurement of the nHz GW signal from individual sources is likely to be possible from the future array of radio antennae such as Square Kilometer Array (SKA) (Ellis 2013;Burke-Spolaor 2013;Kharb et al. 2017).With such observations of individual GW signals, we can fit the astrophysical properties of the galaxies with the frequency dependence of the GW signal and fit the parameters on the occupation number and the signature of the environmental effects on the GW strain by directly comparing the properties of the host galaxy such as the gas density, stellar mass, halo mass, SFR, galaxy morphology, and color.Furthermore, an interesting avenue will be to perform a dedicated study of the hosts of the GW sources with high-resolution spectroscopic surveys to better understand its astrophysical properties.

CONCLUSION AND FUTURE OUTLOOK
In this work, we explore the astrophysical properties of the host galaxies of the SMBBHs which can produce nano-Hertz SGWB using the Romulus25 cosmological simulation.Romulus25 is capable of modeling the astrophysical properties of galaxies and its unique approach to seeding, accretion, and particularly the dynamics of SMBHs makes it especially well-suited for investigating SMBH/SMBBH-galaxy connections.Using this simulation, we have  We found that SMBBHs with chirp mass   > 10 8 M ⊙ are primary source of the SGWB signal.In our simulation up to  = 2, we found six such sources resulting in a number density of 7.7 × 10 −6 cMpc −3 consistent with the results from PTA studies (Mingarelli et al. 2017;Casey-Clyde et al. 2022;Antoniadis et al. 2023c).Although the Romulus25 is a cosmological simulation and not an idealized simulation specifically designed for SMBH-SMBH physics, it still successfully produces a correct number density.This highlights its potential for further studies in this domain.
We then continue by studying the redshift evolution of the astrophysical properties of the host galaxies such as gas density, SFR, stellar mass, and halo mass, across redshifts.These host galaxies are early-type galaxies, characterized predominantly by older star populations.They exhibit a distinct trend in the color-magnitude diagram across redshifts, which could be of particular interest to compare with observations.
Our analysis further reveals that, compared to their counterparts at similar redshifts, the host galaxies of nHz GW sources exhibit lower SFRs, greater stellar masses, and more substantial halo masses.Our findings collectively suggest that nHz GW sources are predominantly hosted by massive early-type galaxies at the centers of groups and clusters.However, we assert that the typical hosts for these GW sources are expected to be group-central galaxies.This is supported by two main factors: (i) groups are more common than clusters, and (ii) the lower velocity dispersion in groups leads to more effective dynamical friction, thereby increasing the likelihood of mergers, especially between satellite galaxies and the central galaxy of the group (O' Sullivan et al. 2017;Oppenheimer et al. 2021).
It is important to note that our conclusions remain robust even if the seed mass in Romulus25 were set lower than 10 6  ⊙ .As discussed in Section 5.1, the SGWB power spectrum signal in the nHz regime is predominantly due to SMBBHs with chirp masses greater than 10 8  ⊙ .These systems are composed of individual SMBHs, each with a mass exceeding 8 × 10 7  ⊙ , which is significantly above our minimum SMBH mass.Thus, our findings regarding the properties of host galaxies of nHz GW sources remain valid.The theoretical connection of the host galaxy properties of the GW sources and the black hole masses indicates which kind of galaxies and their evolution are linked with the black hole merger.This theoretical connection shown in this work will be a guideline for us to explore the connections from GW observations in the nHz band and optical and infrared galaxy observations.By measuring the spatial cross-correlation between the anisotropic SGWB with galaxies as well as a targeted search of individual galaxies for the nHz GW events in the SKA era.
The multi-messenger technique by exploring the connection between the astrophysical properties of host galaxies with the SMBHs and the strain of the GW signal from the coalescing SMBHs will make it possible to establish from observations how the SMBBHs evolution depends on the astrophysical properties of the galaxies.The occupation number of SMBHs in galaxies of different types will make it possible to test theoretical models using observations.In the future with the data from the ongoing International Pulsar Timing Array and upcoming Square Kilometer Array (SKA) (Janssen et al. 2015), we will be able to make high-precision measurements of the nHz GW signal.In synergy with the galaxy surveys up to high redshifts such as Dark Energy Spectroscopic Instrument (DESI Collaboration et al. 2016), Euclid (Laureĳs et al. 2011), Vera Rubin Observatory (LSST Dark Energy Science Collaboration 2012), and Roman Telescope (Akeson et al. 2019) we will make joint estimation of GW and galaxies to unveil the open question of formation of SMBHs and its connection with the galaxy evolution.

Figure 1 .
Figure 1.We show the SGWB strain as a function of the frequency with change in the fiducial values of the parameters  = 0.2,   = 5 nano-Hertz,  = 1,  = 10/3, and  = 1, and also for SMBBHs with chirp mass M c < 10 8 M ⊙ and M c >= 10 8 M ⊙ .

Figure 2 .
Figure 2. The local astrophysical properties of the galaxies: gas density (top), SFR (middle), and stellar mass (bottom) within 1 kpc of the more massive SMBH in merging SMBBHs detected in Romulus25 and contributing most of the SGWB signal (i.e. with chirp masses   ≥ 10 8 M ⊙ ), shown as a function of the host redshift.

Figure 3 .
Figure 3.The vertical line in each panel shows the gas density (first and second row), star formation rate (third and fourth row) stellar mass (fifth and sixth row) within 1 kpc of the more massive SMBH in the merging SMBBH pairs identified as nano-Herz GW source (i.e. with chirp masses   ≥ 10 8 M ⊙ ) at the host galaxy's redshift.The histogram shows the corresponding quantities around the more massive black holes in all SMBH pairs (merging and proximate) in the Romulus25 simulation at the same redshift.

Figure 4 .
Figure 4.The global galactic properties of the galaxies hosting the merging SMBBH pairs identified as nano-Herz GW sources (i.e. with chirp masses   ≥ 10 8 M ⊙ ).The panels show gas density (first panel), stellar mass (second panel), SFR (third panel), and sSFR (fourth panel) as a function of redshift.All quantities are calculated inside a sphere with a 25 kpc radius around the galaxy center.

Figure 5 .
Figure 5. Categorising host galaxies of SMBBHs with   > 10 8  ⊙ as starforming or quenched following definition of Genel et al. (2018).Here we show Δlog(sSFR) = log(sSFR galaxy ) − log(sSFR ridge ) for these host galaxies.The dashed line, which marks 1 dex below the ridge for each redshift, serves as the quenched threshold.The shaded region indicates the "main sequence".See text for detailed definitions.

Figure 6 .
Figure 6.The rest-frame U-V vs rest-frame V-band luminosity (  ) (first panel), vs rest-frame absolute magnitude (  ) (second panel), and vs stellar mass  * (last panel) of the host galaxies of the SMBHs for chirp mass   ≥ 10 8 M ⊙ (represented by big circles) and 10 7 M ⊙ ≤   ≤ 10 8 M ⊙ (represented by small circles).The data points are color-coded as a function of redshift , and the annotations show the exact redshift they are detected.

Figure 7 .
Figure 7. Multi-band composite image of the host galaxies for SMBBHs with chirp mass   ≥ 10 8 M ⊙ is shown with edge-on (top) and face-on (bottom) views at redshifts they are detected at.

Figure 8 .
Figure 8.The age of the Universe when stars in the host galaxies of SMBBHs with a chirp mass   ≥ 10 8 M ⊙ formed.The dashed line represents the age of the Universe at the time the SMBBHs are detected.The values in the top right box, from top to bottom, indicate the stellar mass within 25 kpc from the galaxy center, the median age of the stars, and the fraction of the stellar mass observed that is aged ≤ 1 Gyr.

Figure 9 .
Figure 9.The amount of time taken by the SMBHs to grow the last 50% of its mass it possesses at the time of the merger is shown for sources with chirp mass  c ≥ 10 7 M ⊙ (in blue) and  c ≥ 10 8 M ⊙ (in orange).The distribution indicates that the SMBHs with chirp mass  c ≥ 10 8 M ⊙ need at least 10% of the age of the Universe to grow, indicating these objects are likely to be host in old galaxies.

Figure 10 .
Figure 10.Left panel: The chirp mass of SMBBHs plotted against the stellar mass of their host galaxy.Right panel: Mass of the more massive black hole in SMBBHs versus the stellar mass of their host galaxy.SMBBHs with a mass ratio  < 0.1 are represented by stars, while those with  > 0.1 are represented by circles.Large symbols correspond to SMBBHs with chirp mass   ≥ 10 8 M ⊙ and small symbols denote those with 10 7 M ⊙ ≤   ≤ 10 8 M ⊙ .The data points are colored according to redshift , and the annotations show the exact redshift they are detected at.

Figure 11 .
Figure 11.vertical line in each panel shows properties of host galaxies and host halo of the merging SMBBH pairs identified as nano-Herz GW sources (i.e. with chirp masses   ≥ 10 8 M ⊙ ) at the corresponding redshift.The corresponding distribution for all galaxies or halos with M 200 > 4.5 × 10 9 M ⊙ in the Romulus25 simulation at the same redshift with a background histogram.First and second row: display the stellar mass.Third and fourth row: present the sSFR.Fifth and six row: show M 500 of the host halos.All properties are measured at the host galaxy's redshift.The stellar mass and sSFR are calculated within a 25 kpc sphere around the galaxy center.

Figure 12 .
Figure 12.Top panel: Mass of most massive halo in SMBBH pairs identified as nano-Herz GW sources as a function of redshift.Bottom panel: The M 500 of host halos of nano-Herz GW sources vs. redshift.