Draft version December 16, 2024

Typeset using LATEX twocolumn style in AASTeX631

Jetted Seyfert Galaxies at z = 0:

Simulating Feedback Effects on Galactic Morphology and Beyond
Julianne Goddard,1 Isaac Shlosman,1, 2 and Emilio Romano-Diaz3
1 Department of Physics and Astronomy, University of Kentucky, Lexington KY 40506-0055, USA∗
2 TheoreticalAstrophysics, School of Sciences, Osaka University, Osaka 560-0043, Japan†
3 Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
arXiv:2412.09679v1 [astro-ph.GA] 12 Dec 2024

ABSTRACT
We use high-resolution cosmological zoom-in simulations to model feedback from Seyfert-type super-
massive black hole (SMBH) jets onto galaxies with identical dark matter (DM) halos of log M/M⊙ ∼
11.8. The low mass, ∼ 106 M⊙ , seed SMBHs, have been introduced when the parent DM halos have
reached log M/M⊙ ∼ 11. In a controlled experiment, we vary only the efficiency of the SMBH accretion
and focus on galaxies and their immediate environment properties. Our results show that the AGN jet
feedback has a substantial effect on the basic properties of Seyfert-type galaxies, such as morphology,
gas fraction and distribution, star formation rate and distribution, bulge-to-disk ratio, DM halo baryon
fraction, and properties of circumgalactic medium (CGM) and beyond. These have been compared
to a galaxy with supernovae only feedback. We focus on the energy deposition by the jet in the ISM
and IGM, and follow the expansion of the multiple jet cocoons to ∼ 2 Mpc. We find that the jet-
ISM interaction gradually pushes the star formation to larger radii with increasing accretion efficiency,
which results in increased mass of the outer stellar disk, which is best fit as a double-exponential disk.
Furthermore, we compare our galaxies and their properties with the observed nearby Seyfert galaxies,
including the scaling relations, and find a close agreement, although statistical analysis of observed
Seyferts is currently missing. In a forthcoming paper, we focus on evolution of these objects at z < ∼ 10
and study the effect of the SMBH seeding redshift.

Keywords: AGN host galaxies (2017) — Circumgalactic medium (1879) — Galaxy formation (595) —
Hydrodynamical simulations (767) — Relativistic jets(1390) — Seyfert galaxies (1447)

1. INTRODUCTION magnetically (Blandford & Payne 1982; Emmering et al.

Supermassive black holes (SMBHs) are ubiquitous in 1992; Konigl & Kartje 1994; de Kool & Begelman 1995;
the universe and reside in galactic centers, serving as en- Bottorff et al. 1997; Elitzur & Shlosman 2006; Fukumura
gines of active galactic nuclei (AGN), in the mass range et al. 2014), forming the broad emission and absorption
of ∼ 106 − 1010 M⊙ . About 10% of AGN are radio loud line regions (BEL and BAL). Additional phenomena,
and host relativistic jets (Begelman et al. 1984). The like narrow line regions, warm absorbers, and powerful
AGN generate energy across the electromagnetic spec- molecular outflows are the by-products of AGN activity
trum, as well as mechanically, in the form of collimated which can deposit both energy and momentum in the
(i.e., relativistic jets) and uncollimated outflows in the host galaxy and beyond, in the galactic, circumgalactic
form of accretion disk winds. While jets are collimated (CGM) and intergalactic gas (IGM).
by magnetic fields (Blandford & Rees 1974; Blandford However, the essential details of the AGN feedback,
& Payne 1982; Contopoulos 1995), the disk winds are e.g., how, where and in what way, the AGN energy is de-
either driven radiatively, e.g., by the UV line radiation posited, are still uncertain and subject to intense inves-
(Shlosman et al. 1985; de Kool & Begelman 1995; Mur- tigations. A significant observational evidence points to
ray et al. 1995; Arav et al. 1997, and refs. therein), or the AGN feedback curtailing the high-mass end function
of galaxies, e.g., by quenching the star formation, while
the stellar evolution feedback, e.g., supernovae (SN), de-
∗ E-mail: [email protected] termines the low mass-end function (Kauffmann et al.
† E-mail: [email protected] 2003; McNamara & Nulsen 2012).
2

In this work, we use high-resolution cosmological 0313-192 (Ledlow et al. 2001), 650 kpc in Speca J1409-
zoom-in simulations to investigate the mechanical feed- 030 (Hota et al. 2011), 1.6 Mpc in J2345-0449 (Bagchi
back from relativistic jets in lower luminosity AGN, the et al. 2014), and more. Galaxies with large jets also
Seyfert galaxies. We aim at understanding the effects of display multiple episodes of launching them, e.g., NGC
the energy and momentum deposition by these bipolar 2639, which exhibits relic radio lobes from 3 pc to 9 kpc
jets on their host galaxies, i.e., on the gas distribution in (Rao et al. 2023), typically misaligned1 .
the ISM, CGM and IGM, star formation, stellar masses, About 7% of narrow line Seyfert 1 (NLS1) galaxies
the resulting galactic morphology, and other relevant pa- are radio-loud and, therefore, harbor less powerful rel-
rameters which characterize a galaxy. In the present ativistic jets compared to radio quasars (e.g., Komossa
work, we focus on a MW-type galaxy within a cosmo- et al. 2006). They accrete at a rate close to the Edding-
logical environment, at z = 0. We seed the SMBHs and ton (e.g., Komossa 2018, and refs. therein). But recent
vary only the efficiency of the SMBH accretion, and thus detections of jets in radio-quiet NLS1s have put into
the strength of the AGN feedback. These results have question the usage of radio loudness as an indicator of
been compared to a model without the SMBH. Only SN a relativistic jet whose radio emission could be shielded
feedback and AGN feedback from jets have been used by the ionized gas (e.g., Berton et al. 2020). The radio
here, in order to have a controlled experiment. In a quiet NLS1s can be subject to low efficiency accretion
subsequent work (Goddard et al., in prep., Paper II), we (Heckman & Best 2014).
focus on the evolution of these galaxies at z < ∼ 10 and To summarize, Seyfert host galaxies are typically ra-
vary the seeding time. A number of simulations have dio quiet, but are jetted on all scales, from ∼ 1 pc to
looked at the AGN jet feedback effects in Milky Way- 100s of kpc, including multiple periods of jet activity.
type galaxies over cosmological timescales (e.g., Byrne It is natural to ask the question, what effects jets in
et al. 2024; Irodotou et al. 2022; Wellons et al. 2023), Seyferts have on the ISM of the host, the star forma-
and we aim to expand this body of work. tion, and ultimately on the morphology of these galaxies.
Seyfert galaxies are classified as less luminous than Because accreting SMBHs have their mechanical and ra-
quasars, i.e., with the bolometric luminosity of Lbol < ∼ diation output correlated with the SMBH masses (e.g.,
1012 L⊙ (e.g., Schmidt & Green 1983) and SMBHs which Kormendy & Richstone 1995; Magorrian et al. 1998), it
are less massive, M• ∼ 106 − 108 M⊙ . Long-baseline ar- is natural to expect that the most massive quasar-type
9
rays have resolved compact pc-scale jets in Seyferts and SMBHs, > ∼ 10 M⊙ , have the largest effect on their host
in low-ionization nuclear emission-line region (LINER) galaxies. As Seyferts have lower luminosity than quasars
galaxies (e.g., Falcke et al. 2000). While short-baseline and less massive SMBHs, they host weaker jets, which
interferometry observations, e.g., VLA, have detected are oriented randomly with respect to the midplanes of
relativistic jets in Seyferts on the kpc-scales and larger galactic disks (Osterbrock & Shuder 1982), and so it is
(e.g., Ho & Ulvestad 2001; Rosario et al. 2010). Ex- more difficult to disentangle the influence of jets from
tended radio emission from Seyferts and LINERs has stellar evolution in these objects.
been explained by jet-ISM interaction (e.g., Gallimore Jets interact with the environment by depositing both
et al. 2006). High-velocity winds have been also diag- energy and momentum in the surrounding gas. For a
nosed in UV/optical (e.g., Komossa 2018), similarly to collimated jet, its working surface, i.e., the ‘hot spot,’
quasars (e.g., Emmering et al. 1992). Half of the sam- is relatively small, but the main interaction proceeds
ple of 102 nearby radio-quiet Seyferts and LINERs show through its expanding cocoon. The cocoon consists of
jetted structures, which is a higher fraction than in pow- a shocked jet and ambient matter (Blandford & Rees
erful quasars (e.g., Baldi et al. 2018). Bends and wiggles 1974). Its expansion along the jet, hence, results from
have been frequently observed in Seyfert jets connecting the balance between the jet thrust and the ram pres-
emission at different spatial scales, and jet precession sure of the ambient matter. While in the perpendicular
has been invoked to explain them (e.g., Saikia & Jam- direction, the cocoon expands based on its internal pres-
rozy 2009; Hada et al. 2013). sure. Hence, the cocoon is overpressured with respect to
A number of radio quiet Seyferts display jets on galaxy the ambient gas and drives a shock (Begelman & Cioffi
scale and larger, e.g., 2 kpc jet in NGC4258 (Cecil et al. 1989). Therefore, the cocoon can have a significant effect
2000), 12 kpc jet in NGC 7479 (Laine & Beck 2008), on the ISM, and, if it breaks out of the galaxy, then also
5 kpc in NGC4388 (Damas-Segovia et al. 2016). While
‘normal’ Seyferts in disk galaxies often host kpc-size jets 1 These disk galaxies with such extended radio jets should be clas-
(Keel et al. 2006), some stand out, like spectacular jets sified as radio galaxies rather than radio quiet Seyferts. Never-
of 100 kpc in J1649+263 (Mao et al. 2015), 350 kpc in theless, we emphasize that they are normal disk galaxies.
 3

on the circumgalactic and intergalactic medium, CGM galaxies (Di Matteo et al. 2005), but also in a broader
and IGM, and on the associated processes there. range of galaxy masses. Wellons et al. (2023) looked
In our modeling of the Seyfert jet interaction with the at (∼ 1010−13 M⊙ ) galaxies evolved over cosmological
surroundings, we pay a special attention to the energy timescales with a wide range of AGN feedback mecha-
deposition by the jet, which has the most profound ef- nisms. Of the many models they explored, collimated jet
fect on the star formation. In low luminosity AGN, the feedback was one of only a few that were able to match
feedback from the jets is more difficult to separate from galaxy scaling relations, e.g., M∗ − Mhalo , SFR-M∗ , and
the feedback from stellar evolution, such as stellar winds M• − σ relations, across the mass range.
and SN. Irodotou et al. (2022) studied the effect of the ’radio-
Efficiency of star formation in galaxies which form in- mode’ feedback on MW-mass galaxies by modeling ex-
side DM halos depends on many factors. Observations panding bubbles in the CGM and compared this with
indicate that that stars form most efficiently per unit isotropic mechanical feedback (e.g., winds), and with
DM mass in the range of Mh ∼ X × 1012 M⊙ . For lower a no feedback model. This study found that increas-
or higher DM halo masses, this efficiency decreases (e.g., ing feedback strength decreases galaxy stellar mass, and
Behroozi et al. 2019; Piotrowska et al. 2021). Feedback that varying the type of feedback can have a strong in-
from stellar winds and SN is responsible for this decrease fluence on the galactic morphology. Byrne et al. (2024)
in lower mass halos and correspondingly for lower mass used initially collimated jet feedback in MW-mass and
galaxies embedded in these halos. intermediate mass galaxies along with other feedback
For higher mass halos and associated more massive prescriptions, and found that introducing feedback does
galaxies, the feedback from stellar evolution is proba- reduce the star formation in both mass regimes. Addi-
bly insufficient to decrease the star formation efficiency. tionally, they found that feedback is able to affect mor-
Instead, energy deposition by AGN is suspected to be phology, either increasing the half-mass radius of the
responsible for this effect (Kauffmann et al. 2003; McNa- galaxy or trending toward more spheroidal distributions
mara & Nulsen 2012). From observations we do see that with a stronger feedback.
the prevalence of AGN peaks in a timeline consistent Mukherjee et al. (2018) and Talbot et al. (2022) sim-
with the quenching of many massive galaxies (Kauff- ulated jets over relatively short timescales, <∼ 1 Gyr, in
mann & Haehnelt 2000), which provides further evi- isolated Seyfert-type galaxies, and studied the jet-driven
dence that these two phenomena are correlated. Addi- outflows and their effects on the galaxy and its environ-
tionally, large-scale cosmological numerical simulations ment. Talbot et al. (2022) observed outflows launched
which include AGN feedback have been able to better into the CGM by low power jets and found that the jet
match the observed stellar mass - halo mass relation in direction has a significant influence on the composition
comparison to those without (e.g., Bourne & Yang 2023; and velocity of the outflow. Jets that propagate into the
Rennehan et al. 2024). disk of the galaxy tend to launch slower, colder outflows.
Higher resolution numerical simulations studying the Mukherjee et al. (2018) observed that effects on the SF
influence of AGN feedback in massive galaxies and depend on both the jet angle and its strength, as well as
galaxy clusters have also shown that AGN are excel- on the timeframe. Upon onset, the jets seem to trigger
lent candidates to slowdown and even halt the cooling the SF, but over time they have dissolved cold clouds
flows, and quench the star formation on massive scales through enhanced turbulence, and the overall SF may
(Martı́ 2019; Su et al. 2021; Komissarov & Porth 2021). start to decrease. Additionally, stronger jets and those
Several of these studies have found that feedback in the directed into the disk, are more likely to trigger the SF.
form of jets is especially efficient at heating gas in the Appleby et al. (2021) analyzed the effects of jet feedback
CGM and slowing the cooling flows, and that for galax- on the CGM in the SIMBA simulations. They observed
ies with log(M∗ ) ≥ 10 M⊙ it may be the only way to a reduced baryon fraction and increased metallicity of
fully quench the galaxy (Dutta et al. 2024; Scharré et al. the hot gas in the CGM, due to the presence of jets.
2024). There is a strong degeneracy in the AGN numerical
A number of scaling relations between AGN and their recipes that are able to match observations, as follows
host galaxies have emerged based on observations, e.g., from exploration of the parameter space of the AGN
the M• − σ relation, which relates the mass of the models (e.g., Su et al. 2021; Wellons et al. 2023). A
SMBH with the bulge velocity dispersion (Gebhardt wide variety of AGN numerical recipes replicate ob-
et al. 2000; Ferrarese & Merritt 2000). Promisingly, nu- served galaxy scaling relations and efficacy in quenching
merical simulations have been able to replicate many the star formation.
of these observed scaling relations not only in massive
4

Due to the smaller sample of our models, we look in subset. The final properties of this halo are displayed in
depth into the effects of the simulated jets on various Table 1. Additional results involving other halos will be
parameters of evolving Seyfert-type galaxies over cosmo- published elsewhere.
logical timescales, comparing them with the non-jetted
one. We do verify the galaxy compliance with the major Table 1. DM Halo Properties at z = 0
scaling relations. Our focus is on quantitative morpholo-
gies of modeled galaxies, on the bulge and disk com- Simulation Type log Mvir /M⊙ Rvir [kpc] λ δ
ponents, distribution of star formation, accretion rates DM only 11.8 231 0.03 2.0
on the SMBH, effects on the baryonic components of baryonic 11.8 230 0.02
galaxies, and the CGM. We follow the expansion of the
Note—The columns, from left to right, represent the DM
cocoons generated by the jets and on the distribution
halo virial mass, Mvir , in the DM-only zoom-in and bary-
of the energy deposition by the jets. In this paper we onic zoom-in simulations, the halo virial radius, Rvir , in
present the results at z = 0, and in a forthcoming paper both simulation types, and the DM halo spin, λ.
(Paper II), we shall present the cosmological evolution
at z <
∼ 10. The outline of this work is as following: sec-
tion 2 deals with the numerical methods used, section 3
provides the results, which are discussed and summa-
rized in section 4.
The zoom-in ICs are composed of five nested levels
2. NUMERICS of refinement on top of the base grid, i.e., from 27 to
212 . The DM-only version was first evolved in order to
2.1. Simulation setup check for, and avoid contamination from massive, lower-
We have performed a suite of cosmological zoom-in resolution particles in the highest resolution-level vol-
simulations of a single halo using the N -body/hydro ume. Subsequently, baryons were included at the high-
code gizmo (Hopkins 2015). For better angular mo- est level of refinement in the reconstruction of their re-
mentum conservation and in order to resolve the Kelvin- spective ICs.
Helmholtz instability which is expected to develop when Within this setup, the effective number of particles
the cosmological filaments penetrate the DM halo, the (DM and baryons) for our simulations is 2 × 4, 0963 , i.e.,
MFM hydro solver was employed. The initial condi- this is the number of particles under assumption that
tions (IC) of the simulations used the Planck Collab- the entire simulation box has the resolution of the zoom
oration et al. (2016) ΛCDM concordant model, with region. This leads to a mass resolution per particle of
Ωm = 0.308, ΩΛ = 0.692, Ωb = 0.048, σ8 = 0.82, 3.6 × 104 M⊙ for the gas and stars, and 1.9 × 105 M⊙
and ns = 0.97. We use the Hubble constant h = 0.678 for the DM in the zoom region. The minimal adaptive
in units of 100 km s−1 Mpc−1 . These were generated at gravitational softening in comoving coordinates for the
redshift z = 99 by using the music code (Hahn & Abel gas is 1 pc, for stars 20 pc and 200 pc for DM.
2011) within a box of 50 h−1 Mpc, and were evolved to Galaxies have been identified by the group-finding al-
z = 0. gorithm hop (Eisenstein & Hut 1998), using the outer
From the parent, uni-grid, DM-only simulation, a set boundary threshold of baryonic density of 10−4 nSF crit =
of halos was chosen for re-simulation at a higher resolu- 10−2 cm−3 , which ensured that both the host starform-
tion. DM halos and their properties were identified by ing gas and the lower density non-starforming gas are
the group finder rockstar (Behroozi et al. 2012), with roughly bound to the galaxy (Romano-Dı́az et al. 2014).
a Friends-of-Friends (FoF) linking length of b = 0.28. This assures that identified galaxies are not imposed
The DM halo virial radius and the virial mass, Rvir and with a particular geometry.
Mvir , have been defined by R200 and M200 (e.g., Navarro
et al. 1996). R200 is the radius within which the mean 2.2. Star formation and SN feedback
interior density is 200 times the critical density of the Gas heating and cooling from 1010 K down
universe at that time, and M200 the corresponding en- to 10 K are implemented, including H and He
closed mass. For baryonic+DM halos we use Mhalo and ionization+recombination, collisional, free-free, dust
Rhalo abbreviations. We selected a subset of halos with collisional, cosmic ray, and Compton effects, as well as
a final mass range Mvir =1011.5 − 1012 M⊙ , a local en- metal-line (Wiersma et al. 2009), fine-structure, and
vironment with relative overdensity of δ = 2 ± 1.0, and molecular cooling, as detailed in Hopkins et al. (2018)
halo spin of λ = 0.03 ± 0.02. This paper focuses on and Hopkins et al. (2022). Metal enrichment is in-
the results using just one of the halos selected from this cluded: the metallicity increases in the starforming gas
 5

and scales with the fraction of stars that turn into SN, the spin axis of the SMBH evolves with time. We in-
and the metal yield per SN (see below). A total of 11 clude an artificial velocity-damping term to continuously
metal species were followed in both in gas and stars, move the SMBH toward the most bound particle (e.g.,
including H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe. Wellons et al. 2023). Without this, the SMBHs have a
The H2 abundances used for cooling calculations are well documented tendency to wander from the center,
estimated from the Krumholz & Gnedin (2011) analytic being initially less than one order of magnitude more
fitting function. massive than the high resolution DM particles.
Metal diffusion is not implemented explicitly, but met- The SMBHs are seeded with M• ≃ 106 M⊙ and grow
als can be transported by mechanical feedback from SN by accreting surrounding gas, with an accretion rate
and AGN (see below and section 2.3). Our simulations determined based on gravitational torques (Shlosman
include the redshift-dependent cosmic UV background et al. 1989), calculated using Hopkins & Quataert (2011)
(e.g., Faucher-Giguère 2020; Shen et al. 2020). method,
The density threshold for star formation (SF) was set  1/6
−3
to nSF
crit = 100 cm . Stars form only where gas is self- 3/2 M•
Ṁgrav = αfd ×
gravitating, namely, when 108 M⊙
  −3/2  −1
Md (< R0 ) R0 f0
(∆v)2 + 2c2s < 8πGρ . (1) 1+ ,
109 M⊙ 100 pc fgas (< R0 )
For this, we rely on the virial parameter, αvir , defined (4)
as where α is a normalization factor motivated by differ-
αvir = [(∆v)2 + 2c2s ]/8πGρ , (2) ences in the SF criteria and set to a default value of
α = 5. R0 is the kernel radius of the SMBH, and roughly
such that the self-gravity condition is met when αvir < 1
corresponds to the SMBH radius of influence, limited to
Hopkins et al. (2013). Star formation efficiency (SFE) is
a minimum value of 80 pc. This is the radius inside
calculated from the virial parameter following the model
which the other parameters, e.g., Md (<R0 ) — stellar
by Padoan et al. (2012):
disk mass, fgas (<R0 ) — gas fraction, and fd — disk
1/2 mass fraction of the total mass, are evaluated. Finally,
SF E = exp(−1.4αvir ) . (3)
 1/3
Once the SFE is determined, stars form stochastically Md (< R0 )
f0 = 0.31fd2 . (5)
(Springel & Hernquist 2003). Each star particle repre- 109 M⊙
sents an entire population of stars with a mass distribu-
We multiply the accretion rate by an efficiency pa-
tion following the Chabrier (2003) Initial Mass Function
rameter, ϵ, introduced by Anglés-Alcázar et al. (2017)
(IMF).
to make Ṁgrav more realistic and to capture the effects
Mechanical feedback from SN–type II is implemented
of unresolved processes affecting gas inflow that are not
using the numerical prescription given in Hopkins et al.
addressed in the existing subgrid recipe.
(2018). We follow individual SN events with an assump-
In this work, ϵ is the sole parameter varied between
tion that the SN occur with an IMF-averaged constant
the AGN models. We set it to ϵ = 0, 4.5%, 15%, and
rate of 3 × 10−4 events per Myr and per M⊙ for all stars
50%, and denote these models as ϵ0 , ϵ5 , ϵ15 and ϵ50 ,
less than 30 Myr old. Each SN event injects energy of
respectively. So the final accretion rate is
1 × 1051 erg, mass of 14.8 M⊙ , and the metal mass of
2.6 M⊙ into the surrounding gas within the radius of Ṁ• = ϵṀgrav . (6)
200 pc only.
The efficiency parameter ϵ scales down the accretion
rate onto the SMBH, Ṁgrav , defined by Eq. 4. Smaller
ϵ decreases the jet power, i.e., lowering the feedback
2.3. The SMBH and its mechanical feedback strength. Note however that the jet power depends also
The SMBHs are modeled as sink particles that are on Ṁgrav , which varies from model to model even with
seeded once the DM halo has reached 1011 M⊙ . They identical initial conditions. When we refer to ”AGN
form from gas particles that are located at the minimum strength,” or to ”feedback strength,” we refer solely to
of the gravitational potential of the parent halo and in- the value of ϵ.
herit the initial mass, velocity, angular momentum, and
position of the parent particle. The SMBH also inher-
its angular momentum from the gas it accretes, so that
6

Table 2. Model Galaxies Properties at z = 0

Model log M∗ /M⊙ log Mgas /M⊙ fgas log SFR log Zgas /Z⊙ log Z∗ /Z⊙

ϵ M⊙ yr−1
ϵ0 10.97 9.95 0.09 0.30 -0.77 0.04
ϵ5 10.93 9.82 0.07 -0.44 -0.98 0.05
ϵ15 10.66 9.98 0.17 -0.77 -0.81 0.00
ϵ50 10.58 9.52 0.08 -2.72 -0.84 -0.02
Note—Columns: (1) model name; (2) stellar mass inside HOP-defined galaxies; (3) gas mass inside HOP-defined galaxies; (4)
gas fraction fgas = Mgas /(Mgas + M∗ ); (5) SFR averaged over the final 30 Myr (ϵ50 is averaged over 60 M⊙ to obtain a signal);
(6) gas metallicity; (7) stellar metallicity.

The SMBH accretion rate measured as a fraction of 3. RESULTS

the Eddington accretion rate, ṀEdd , is Here we present the results based on our four model
galaxies in their final state only, at z = 0. All models
fEdd = Ṁ• /ṀEdd , (7)
have been evolved from the same initial conditions, and
where ṀEdd = LEdd /c2 , and LEdd is the Eddington limit differ only in their AGN accretion efficiency, ϵ. Paper II
for M• . will focus on the cosmological evolution of the properties
To model the AGN jets, we use the hyper-refined par- presented here. The SMBH has been seeded when the
ticle spawning (Torrey et al. 2020), modified by Su et al. DM halo mass has reached ∼ 1011 M⊙ at z = 3.7, in
(2021). Particles of mass 103 M⊙ are spawned at a fixed three out of the four models presented here. At this
fraction of the SMBH kernel radius with an initial veloc- redshift, the stellar mass of the HOP-defined galaxy is
ity of 3 × 104 km s−1 and a temperature of 1010 K. This M∗ ∼ 1.3 × 1010 M⊙ , and the gas mass of this galaxy is
gives an initial energy per particle budget of about 22% Mgas ∼ 1.6 × 109 M⊙ , i.e., ∼ 12%. Prior to seeding, the
emitted as thermal energy and 78% as kinetic energy. SFR is ∼ 10 M⊙ yr−1 .
These particles are launched along the spin axis of 3.1. Galaxies: global parameters at z = 0
the SMBH with a mass loading η=0.1, meaning 10% of
In Figure 1, we show the gas and stellar HOP-selected
the accreted gas is returned in the form of jet particles.
galaxies, on a scale 54 kpc × 54 kpc, centered on the
The launching rate is thus Ṁjet = η Ṁ• , giving a jet
potential minimum of the galaxies. All of our mod-
mechanical energy injection rate of
els display a distorted gaseous component at this time.
Ljet = 1/2η Ṁ• v 2 , (8) The source of this distortion and detailed analysis of the
galaxy stellar components are discussed further in sec-
where v is the jet particle velocity. Velocities of the tion 3.2. Only the ϵ50 model has an extended central
spawned particles are initially perfectly collimated, i.e., cavity in the gas — obviously the consequence of the
launched with a zero opening angle along the spin axis AGN feedback, as we discuss here and in Paper II. We
of the SMBH. We refer to this feedback as a ’jet’ feed- also see evidence of recent gravitational interactions in
back. Once created, the jet particles interact hydro- the ϵ15 and ϵ50 models, both of these galaxies feature
dynamically in the same way as any of the other gas tidal stellar streams surrounding the galaxy, and the ϵ50
particles. When they decelerate to at least one fourth galaxy includes a small satellite that does not appear in
of their initial velocity, and enter the kernel radius of the HOP-selected components of the other three models
another gas particle in a head-on trajectory, the jet par- at this redshift.
ticles re-merge and the mass-weighted properties of the The main properties of the simulated galaxies are
two particles are averaged. Note that here and through- listed in Table 2. Those have been calculated within
out the paper, when we refer to the jet luminosity HOP-defined galaxies. The total stellar mass within
we reference only the mechanical luminosity. As men- this radius is ∼ 1011 M⊙ for the ϵ0 model, and a fac-
tioned above, the thermal component of the jet energy tor of ∼ 2.5 smaller for the highest AGN accretion, ϵ50 ,
is non-negligable, and contributes to significant heating model. The stellar–mass decrease with increasing ac-
of the ISM and CGM, particularly when the jet particles cretion efficiency is nearly monotonic. The gas fraction,
merge. This is discussed further in section 3.7. defined as the gas mass in units of total baryonic galaxy
 7

Figure 1. Projected surface density of the HOP-selected galaxies at z = 0 rotated based on the angular momentum, J, of the
central 3 kpc stellar disk. Orange color represents the gaseous component of a galaxy, while blue color represents the stars. The
columns show the ϵ0 , ϵ5 , ϵ15 , and ϵ50 models respectively from left to right. The top row images are the face-on projections,
and the bottom two rows are the two perpendicular edge-on projections. Each frame is 54 kpc × 54 kpc.

mass, is fgas ∼ 0.1. It shows a spread by a factor of plays our models on this diagram adopted from Behroozi
2, but in order to understand the dependency on ϵ, one et al. (2019), Guo et al. (2010) and Hudson et al. (2014).
should follow the galaxy evolution over an extended time All our models lie within 0.5 dex above the median re-
period. The gas distribution within the galaxy will be lation for the Guo et al. (2010) relation. A trend can
discussed later. be observed — the models approach the median with
The star formation rate (SFR), calculated over time increasing efficiency. Positions of the ϵ0 and ϵ5 models
intervals of 30 Myr, exhibits a dramatic decline with ϵ, indicate that their feedback may not be strong enough
by about three orders of magnitude. For ϵ50 , this means to prevent over-cooling, and result in larger M∗ .
essentially quenching the star formation process. The Over-cooling is a long-standing issue in numerical sim-
gas metallicity does not correlate with the accretion ef- ulations (e.g., Shlosman 2013, and refs. therein), still in
ficiency, and is smaller by a factor of a few compared to need of intense investigation, as it is found not only here,
the stellar metallicity, which is about solar. This obvi- but in many other recent and similar bodies of work (i.e.,
ously points to being diluted by the gas accretion — we Byrne et al. 2024; Chen et al. 2023). The addition of a
return to this issue in section 3.5 and Paper II. stronger AGN accretion efficiency does decrease M∗ and
Based on the global parameters of modeled galaxies, brings it closer to the M∗ − Mhalo median, but under-
we have checked their positions on the stellar mass – halo standing of the AGN feedback is far from satisfying.
mass, M∗ − Mhalo , diagram. Figure 2 (left frame) dis-
8

Figure 2. Properties of our HOP-selected galaxies at z = 0, as indicated in the legend. (a) Shows positions of modeled galaxies
on the M∗ − Mhalo scaling relation adopted from Behroozi et al. (2019) (solid blue line), Guo et al. (2010) (dotted red line), and
Hudson et al. (2014) (dashed yellow line). The filled regions around the medians correspond to a scatter of 0.5 dex. (b) Provides
specific SFRs for the HOP-galaxies and are plotted alongside the star-forming galaxy main sequence fits derived in Popesso
et al. (2023) (solid line — their eq.10, dashed line — their eq.14). The colors indicate divisions between the star-forming (blue),
green valley (green), and quiescent (red), as delineated in Salim (2014). (c) The solid red line shows the M• - σ relation fit
from Kormendy & Ho (2013), and the red fill shows the 0.5dex scatter. The 3-D dispersion velocity σ has been corrected for
the rotational velocity (Stewart & Federrath 2022).

The ϵ0 model follows the SF main sequence, as indi- in Figure 15 there. The summary of the fits are given in
cated by the specific SFR (sSFR) provided in Table 2 Table 3.
and Figure 2 (middle frame). The colors in this frame
indicate divisions between star-forming - blue, green val-
ley - green, and quiescent - red, as delineated in Salim The bulge radii have been determined using the inter-
(2014). The final sSFR for ϵ15 and ϵ5 models fall within section of the fitting curves for each component. The ϵ0
the green valley region, while the ϵ50 model is fully and ϵ5 model bulges are similar, 1.58 kpc and 1.81 kpc.
quenched. The values provided for the sSFR have been The models with a higher accretion efficiency, ϵ15 and
averaged over the last 30 Myr of the simulation, except ϵ50 , have smaller bulges, 0.63 kpc and 0.68 kpc, respec-
for the ϵ50 model, which has been averaged over 60 Myr. tively. The resulting bulge masses exhibit a monotonic
The bulge masses and bulge velocity dispersions have decline, by a factor of 5, with increasing accretion effi-
been obtained from Sersic morphological decomposi- ciency. The inner stellar disk mass shows such a decline
tion (section 3.2), to determine the positions on the as well, albeit by a factor of 2 only. The outer stellar
Mbulge − σ relation. The 3-D velocity dispersion has disk mass increases with ϵ by a factor of 5, except for
been corrected for contribution from rotational velocity the ϵ50 model, while the total stellar masses in galaxies
(Stewart & Federrath 2022). The results for our galax- decline by a factor of 3. Note that the total stellar mass
ies are shown over-plotted on data from Kormendy & of the HOP-selected galaxies corresponds closely to the
Ho (2013) in Figure 2 (right frame). All AGN galaxies stellar mass within 0.1Rvir .
lie within the 0.5dex scatter region from the observa- Table 3 shows that lower accretion efficiency models,
tional data. We elaborate on the evolution tracks of ϵ0 and ϵ5 , have similar parameter values among them-
AGN galaxies in this diagram in Paper II. selves, with a dominant bulge, B/D∼ 1.12 − 1.20 (D be-
ing the stellar disk mass and B being the bulge mass.)
3.2. Quantitative morphology The outer stellar disks are less massive than the inner
The bulge and disk components of modeled galaxies ones by a factor of 10. The same trend, B/D∼ 0.28−0.38
have been determined by fitting the 1-D face-on stellar is observed for the higher accretion efficiency, where the
surface density profile of each galaxy within 0.1Rvir , us- total disk dominates the system. The higher ϵ models
ing the Sersic function (e.g., Bi et al. 2022a) for the bulge possess smaller bulge components and disks that domi-
and two exponential disks — the double-exponential nate in mass, and the outer stellar disks which are less
gave the best fits at z = 0. We use the radii 0.1Rvir in massive than the inner ones only by a factor of less than
cylindrical shells instead of the HOP-defined galaxy size, 4.
as they are very similar and easier to calculate. The de- Hence, we conclude that a systematic difference in
tails of the fitting process are shown in appendix A and stellar morphology of modeled galaxies exists along the
 9

Table 3. Quantitative Morphological (Sersic) Decomposition of Modelled Galaxies at z = 0

Model n Rb Re Rdisk1 Rdisk2 B/D B/T log Mbulge /M⊙ log Mdisk1 /M⊙ log Mdisk2 /M⊙ log Mtot /M⊙ Rbar emax

ϵ kpc kpc kpc kpc kpc

ϵ0 1.12 1.58 0.57 1.71 7.99 1.20 0.55 10.67 10.56 9.40 10.93 1.2 0.25
ϵ5 1.28 1.81 0.69 1.68 5.63 1.12 0.53 10.62 10.52 9.60 10.90 2.3 0.50
ϵ15 0.71 0.63 0.37 0.71 3.88 0.28 0.22 9.98 10.37 10.04 10.65 2.0 0.10
ϵ50 1.00 0.68 0.33 0.89 3.84 0.38 0.28 9.97 10.29 9.70 10.53 - -
Note—Fitting the 1-D face-on stellar surface density profile using the Sersic function for the bulge and a double-exponential
disk inside 0.1Rvir . Columns: (1) model name; (2) bulge Sersic index; (3) Rb bulge radius; (4) Re bulge scalelength; (5) Rdisk1
inner disk scalelength; (6) Rdisk2 outer disk scalelength; (7) B/D bulge-to-total stellar disk mass ratio; (8) B/T bulge-to-total
stellar mass ratio; (9) Mbulge bulge mass; (10) Mdisk1 inner stellar disk mass; (11) Mdisk2 outer stellar disk mass; (12) total
stellar disk + bulge mass; (13) Rbar bar or oval distortion radius; (14) emax maximal ellipticity of a stellar bar or an oval
distortion.

sequence of increasing accretion efficiency. Probably the emax and Rbar for all models, if a bar or an oval distor-
most interesting trend is the appearance and increase in tion are present. Here we adopt the definition of a bar
the mass fraction of the outer disk with increasing effi- having emax ≥ 0.4, which approximately corresponds
ciency — the SF is gradually pushed out to larger radii to the ratio of a Fourier amplitude for m = 2 mode,
with the feedback. A2 , normalized by the monopole m = 0 mode, A0 , and
which is A2 /A0 = 0.15 (Martinez-Valpuesta et al. 2006;
Bi et al. 2022b).
The ϵ0 model hosts an oval distortion rather than a
stellar bar, with a maximal ellipticity of emax ∼ 0.25
with a radius of 1.2 kpc. For comparison, A2 /A0 ∼ 0.1,
which justifies our conclusion of an oval distortion. Even
a weaker oval distortion is present in the ϵ15 model, with
emax ∼ 0.1 of a radius of 2 kpc. Not even an oval distor-
tion is present in ϵ50 model. But the ϵ5 model displays
a bar with emax ∼ 0.5 with Rbar ∼ 2.3 kpc.
While only ϵ5 model exhibits a weak bar at z = 0,
other models display such bars at earlier times (Pa-
per II). These bars and oval distortions are typically trig-
gered by prograde tidal interactions, but the retrograde
interactions can also weaken the bars. Additional fac-
tors, such as a rapid influx of gas along the cosmic fila-
ments can have the same effect (Bi et al. 2022b).
Figure 3. Gas surface density of the face-on galaxy at
z=0 measured in cylindrical shells of 1 kpc width and ±5 kpc
The radial 2D distribution of gas in the models, Σgas ,
height. The black arrow indicates 0.1Rvir . The offset gray is more diverse than the stellar distribution, and is
line has been added to show an approximate exponential gas shown in Figure 3. Among all the models, the ϵ0 has the
distribution with the scale length of 3.8 kpc for comparison highest surface density in the central few kpc, and the
only. lowest one in the outer galaxy. In the central ∼ 5 kpc,
the gas surface density forms a plateau, with ϵ0 being
the highest. The ϵ5 model shows a sharp decline in the
We have also analyzed the presence of stellar bars. For
central 3–4 kpc which is the result of the combined ac-
this purpose, we have produced contour maps of stellar
tion of the stellar bar and the AGN feedback. Bars fun-
disk surface density and measured their ellipticity, e(r).
nel the gas inwards, contributing to the gas deficiency
The bar size has been determined by identifying the ra-
within their radii (e.g., Shlosman et al. 1989; Heller &
dius of the maximal ellipticity, emax , and obtaining the
Shlosman 1994). The ϵ50 model displays the strongest
radius of the subsequent 15% decrease from this max-
departure from the stellar surface density, with most of
imum (Martinez-Valpuesta et al. 2006). Table 3 shows
10

Figure 4. 2-D histogram of stellar age vs radius divided into 300 × 300 grid cells with the color representing the stellar mass
in each cell at z = 0. The overplotted solid black line represents the evolution of the gas fraction in galaxies. The stars plotted
here are those residing in the HOP-galaxies at z = 0, and have not necessarily formed within the galaxies.

the gas being pushed to the galaxy edge, as seen in Fig-

ure 1. This effect is also observed in ϵ5 and ϵ15 models,
whose gas surface density outside ∼ 15 kpc exceeds that
of ϵ0 . It affects also the distribution of the SF in the
disk, and we return to this point in the next section.
3.3. Star formation
We can further examine the distribution of stars in
Figure 4. For all the models, about half of the stars have
been formed by z ∼ 2, and another third by z ∼ 1. This Figure 5. The surface density of the SFR vs radius averaged
older population of stars seems to not only populate the in cylindrical shells of height ±5 kpc and 20 logarithmically
bulges, but also the disks, as we observe them at all radii. spaced bins in radius out to 30 kpc. ΣSFR is averaged over
The low AGN efficiency models have stars residing out the last 100 Myr of the simulation to increase the signal and
to ∼ 20 − 25 kpc, while for high AGN efficiency models, minimize the effects of transient behavior in the star forma-
some stars reach to almost 30 kpc. These stars which are tion. The black arrow indicates the radius corresponding to
0.1Rvir .
born at large radii belong to the stellar halo population
and have been accreted during cosmological evolution.
The amount of younger stars decreases with increasing tion and thus the SFR are reduced in all models, for the
effciency, confirming that the ϵ50 model has a quenched stronger efficiency AGN models the feedback works in
SF. tandem to cause a more severe and sometimes a longer-
Figure 4 shows that the distribution of stellar ages cor- lived gas loss. We also observe that the stars formed
relates strongly with the gas fraction in the host galaxy. after z ∼ 1 gradually avoid the central few kpc with in-
This is most visible at z ∼ 1, where we observe a gap in creasing efficiency. They have formed in the outer disk
the stellar age population, which coincides with the min- or been brought in by minor and intermediate mergers
imum in the gas fraction at that time. The prominent (assuming that major mergers would destroy the stellar
gap in stellar age distribution is visible in all models at disk.)
z ∼ 1. The gap becomes more prominent with AGN ef- Interestingly, the gas fraction in galaxies correlates
ficiency, and for ϵ50 , basically determines the quenching with the stellar birth time, only with a large time delay,
of the SF. This reduction in SF seems to be triggered as seen in Figure 4. While the gas fraction declined from
by a series of minor interactions, which substantially re- fgas ∼ 0.25 to < ∼ 0.05 around z ∼ 3, because of numer-
duce the gas fraction in galaxies, and is subsequently en- ous interactions and minor mergers, the SFR declined
hanced by the stellar and AGN feedback. While gas frac- only around z ∼ 1.5.
 11

Figure 6. Face-on gas morphology of the HOP-selected galaxies at z = 0. The white circles display the 0.1Rvir radius for
comparison. The central cavity for the ϵ5 model is the gas response to the stellar bar. The gas in the ϵ50 model has been pushed
away from the stellar disk, which explains the quenching of SF there. Note also the gradual increase in the gaseous disk size
with increasing efficiency.

Note also that the gas fraction, fgas , varies by a factor Table 4. The SMBH-Related Properties at z = 0
of a few over ∼ 3 Gyr time periods after z ∼ 1, strongly
correlating with the stellar age distribution. Finally, fgas Model log M• log Ṁ• log fEdd log Ljet Jet Angle
stabilizes during the last 2 Gyr. ϵ M⊙ M⊙ yr −1
erg s−1
degrees
Figure 5 provides the radial distribution of the SFR,
averaged over the last 100 Myr, and expressed as the ϵ5 7.73 -2.20 -2.28 41.25 20.7
surface density of the SFR. We observe that while the ϵ15 7.86 -2.59 -2.79 40.86 37.0
ϵ0 model SFR peaks close to the center, at ∼ 1 kpc, ϵ50 7.86 -3.26 -3.47 40.19 51.5
models with higher ϵ peak progressively at larger radii. Note—Columns: (1) model name; (2) M• the SMBH final
The ϵ5 model peaks at about 5 kpc, and ϵ15 model at masses; (3) Ṁ• the SMBH accretion rates; (4) the SMBH ac-
∼ 10 kpc. The ϵ50 model shows a very low signal and cretion rates as a fraction of the Eddington rate; (5) Ljet the
only between 20–25 kpc, at the galaxy outskirts. This jet mechanical luminosity; (6) the jet angle with the disk spin
axis. The galaxy spin axis is defined by the stellar angular
distribution of the SFR is well correlated with Figure 4.
momentum vector within the inner 3 kpc. The angle, me-
The gas distribution in modeled galaxies can also chanical luminosity, and accretion rate have been averaged
be visualized in Figure 6. It displays an ϵ-sequence, over the final 30 Myr of simulations.
i.e., based on the AGN accretion efficiency strength.
The most striking observation is the dependency of the
gaseous disk size and its activity of spiral arms on this
efficiency. The ϵ0 model has a gaseous disk located
ther in Paper II. The simplest explanation for this effect
well inside the 0.1Rvir radius. The ϵ5 model exhibits
is that the jet feedback is anisotropic in our models.
a larger disk, supplemented with increased activity of
Mostly, the jet points away from the galactic disk, so
the outer spirals. The ϵ15 model has a disk filling up
the SMBH can continue to accrete material even dur-
the 0.1Rvir radius and extensive spiral arms outside this
ing periods of powerful outflows. The history of the
radius. These models are in sharp contrast with the ϵ50
accretion rate and growth by the SMBH, Ṁ• , differ
model. Here, the gaseous disk is absent and replaced by
as well between the models. The accretion rate, av-
a gaseous ring at 0.1Rvir .
eraged over the last 30 Myr, declines sharply and mono-
tonically from Ṁ• ∼ 6 × 10−3 M⊙ yr−1 for ϵ5 down to
∼ 5 × 10−4 M⊙ yr−1 for ϵ50 , i.e., by more than an order
3.4. The SMBH growth of magnitude.
Table 4 provides the final SMBH masses and other ac- The SMBH jet angle with the disk evolves by inher-
cretion parameters, including the jet mechanical lumi- iting the angular momentum of the accreted gas. The
nosity, Ljet , and the jet angle with respect to the inner angle has been averaged over the last 30 Myr for each
3 kpc stellar disk spin axis. The final SMBH masses model. The jet angle confirms that it points away from
range within M• ∼ 5.4 − 7.2 × 107 M⊙ , so they differ the disk, as the angle ranges between 21◦ − 52◦ . The
by about 30% from each other, increasing with ϵ. The SMBH accretion rate at z = 0, and thus the jet me-
least efficient AGN model has a smaller final SMBH chanical luminosity, appear to scale inversely with the
mass, thus we do observe that accretion rate corre- accretion efficiency. So that the ϵ50 model has Ljet about
lates with the accretion efficiency ϵ, as is explored fur-
12

Figure 7. Radial distribution of gas temperature within 0.1Rvir , at z = 0. The temperature and radius are binned into
100 × 100 cells. The color palette shows the mass within bins.

Figure 8. As Figure 7 but for the ISM metallicity. Black line is the mass-weighted average calculated in 1 kpc wide spherical
shells. The color palette shows the mass within the histogram bins.

an order of magnitude lower than the other two models, While Table 2 shows that the total stellar mass of these
and in accordance with Ṁ• . This result emphasizes that galaxies differs by factor of ∼ 3, their SFRs differ by or-
the current feedback by the SMBH does not reflect the ders of magnitude as a function of the efficiency. The
long time history of the gas supply in the galaxy centers, SF history can be inferred from Figure 4 and differs dra-
as will be shown in Paper II. matically over extended time periods. Why does such a
difference in the SFR result only in a factor 3 difference
3.5. AGN feedback and the ISM in M∗ ? The explanation lies in Figure 4. As we have
Unless otherwise specified, the ISM discussed here estimated in section 3.2, most of the stars in modeled
refers to all of the gas within 0.1Rvir for each of the galaxies have been formed before z = 1. Moreover, half
models. As already seen in the final SFRs of Table 2, of these stars have been formed before the SMBH has
the addition of the AGN feedback decreases the SFR in been seeded in the galaxies, and hence did not experi-
a galaxy with increasing ϵ at z = 0, resulting in a factor ence the AGN feedback.
of ∼ 3 difference in the total stellar mass. Clearly, one should conclude on this point that the
The explanation to this result lies in the history of the SMBH seeds should form in the early stages of the host
SF in the modeled galaxies (Figure 4). While ϵ0 and ϵ5 galaxy evolution in order to have a more profound effect
models exhibit a gradual decline in the SFR over the on the stellar population. This issue is further addressed
cosmological times, the high ϵ models display dramatic in Paper II.
declines after the SMBH has been seeded. The ϵ15 model Table 2 shows the final gas fraction in each of the mod-
nearly recovers its SFR after a couple of Gyrs, the ϵ50 els. At z = 0, all four galaxies show a gas fraction of
model never recovers fully. It takes about 7 Gyr to attain around 7 − 17%. In absolute numbers, some small dis-
∼ 10−3 M⊙ yr−1 after z ∼ 0.3. The dominant effect of persion exists in the gas masses. This dispersion exists
the AGN feedback has been to reduce the SFR in a also for stellar masses in galaxies. In relative numbers,
galaxy.
 13

the gas fraction in ϵ15 stands out, because it is normal- 10−3 − 1. But the radial distribution and the mass of
ized by a smaller stellar mass in this model. the gas which has this metallicity differs between the
Next, we analyze the gas temperature and metallicity models. The higher metallicity in the central regions is
in modeled galaxies. Figure 7 displays the distribution an indication of an inside-out enrichment. The majority
of the gas temperature as a function of distance from of the gas in ϵ0 and ϵ5 models, displays the radial gradi-
the galaxy within 25 kpc, i.e., ∼ 0.1Rvir , coloring it with ent of metallicity in the central few kpc. Beyond 10 kpc,
mass. All the models display a bimodal distribution of most of the gas has Z/Z⊙ ∼ 3 × 10−3 − 10−1 .
gas with temperature. The cold, < 104 K starforming The ϵ15 model is well mixed radially, with most of
gas is present in all galaxies, but its radial distribution the gas having Z/Z⊙ ∼ 3 × 10−2 − 0.3 inside 10 kpc.
differs substantially. In the ϵ0 model, it extends from the Beyond this radius, the distribution is bimodal, with
very center to the outskirts of the galaxy, ∼ 20 kpc, with Z/Z⊙ ∼ 10−2 and 10−3 . Although pushed out, the gas
most of the gas residing at the cooling floor of atomic in the ϵ50 model displays the same bimodality around
hydrogen (although molecular and metal cooling are ac- Z/Z⊙ ∼ 1 and 10−2 .
tive). Already for ϵ5 and ϵ15 , the cold gas does not Note that we do not explicitly model metal diffusion,
extend to the center — the amount of ∼ 104 K gas de- and metals are distributed to the gas through SN feed-
creases at small radii. For the ϵ50 model, the cold gas is back only. So, is the bimodality an artifact? To follow
essentially pushed outside 15 kpc. up, we have analyzed the gas distribution in ϵ15 and ϵ50
models and determined that all models have an offset
outer gaseous ring of a low-metallicity accreting mate-
rial. Recent star formation in this pristine material ac-
creting from the filaments appears to contribute to the
gradient and bimodality observed at large radii for these
models.
To emphasize the presence of an extended gaseous disk
in ϵ50 galaxy, we show the gas distribution within a 240×
240 kpc region in Figure 9. This gas does not appear
in the HOP-selected galaxies shown in Figures 1 and 6,
with the exception of a narrow ring surrounding these
galaxies. This extended disk in tandem with the central
cavity displays the effect of the AGN feedback. While
the stellar component in such galaxies is still present and
displays a normal disk, the SF has been quenched there,
and their colors would represent those of the lenticular
galaxies (This point will be addressed in detail in Pa-
per II).
Figure 9. Projected view of the face-on gas disk in ϵ50 In order to investigate these effects further, we have
galaxy at z = 0. Note the large cavity introduced by the plotted the temperature-density phase diagram for the
AGN feedback. The gas is pushed out of the HOP galaxy,
ISM by coloring it in two-ways (Figure 10). The top
so only a thin ring remains attached, as shown in Figures 1
and 6. The gas outside this ring belongs to the CGM.
row shows the phase diagram with colors representing
the radial distance from the galaxy center (e.g., Bi et al.
6
2024). The bottom row displays the same gas colored
The hot gas with T > ∼ 10 K extends to the central by its metallicity. These frames exhibit a progressively
regions in the ϵ0 model. Not much gas resides in the
reduced higher density gas, especially at >∼ 10 kpc, while
intermediate (warm) temperature between the cold and
lower frames display lower metallicity in the higher den-
hot gas. For ϵ5 model, this trend continues. The hot
sity gas.
gas has a larger distribution in T , which extends even
The top frame of Figure 10 shows that a lower tem-
beyond 107 K. For stronger accretion efficiency, the hot
perature of the ISM is found at higher density, and that
gas is pushed further out in radius, beyond 10 kpc, its
the higher temperature is found at larger radii. The
temperature has decreased below 107 K. Its mass within
bottom frame in this figure offers another view of the
the HOP galaxy is very small, as a result. This gas
gas metallicity, as it relates the temperature and den-
extends well outside the HOP galaxy.
sity distributions in each galaxy colored by metallicity.
We turn to the ISM metallicity (Figure 8). In all mod-
We see that the amount of the highest metallicity gas
els, the dominant gas mass has metallicity of Z/Z⊙ ∼
14

Figure 10. Temperature-density phase diagrams for modeled galaxies within 0.1Rvir at z = 0. The top row is colored by the
particle distance to the galaxy center. The bottom row shows the same ISM colored by its metallicity.

Figure 11. As Figure 10 but showing all the gas in the range of 0.1Rvir − 2Rvir , representing the CGM only. Note, however,
that cold gas residing in the substructure was not removed — it is seen in all models as a distant high metallicity gas.

decreases with ϵ and tends to be ‘pushed’ to lower densi- 3.6. AGN feedback and the CGM
ties. We also observe a decrease in low-metallicity gas at We turn now to the gas properties outside the galax-
high densities (and lower temperature), and a decrease ies, analyzing a good part of the CGM. The top row
of low-metallicity gas at high temperature, which be- of Figure 11 displays the CGM gas diagram within the
comes more diffuse with increasing accretion efficiency. range of of 0.1Rvir − 2Rvir , which is colored similarly
to Figure 10, i.e., by the distance of the individual gas
particles to the center of the main galaxy. Note that the
 15

Figure 12. Gas metallicity in the CGM within 2Rvir . Top row: shows the gaseous component of the CGM with each gas
particle colored by its metallicity. Bottom row: 2D histogram of gas metallicity vs radius. The mass-weighted average metallicity
calculated in 20 kpc wide spherical shells is over-plotted. The color palette shows the mass in each pixel in this phase space.

gas residing in substructures was not removed and can CGM technically. Also, we observe that the hot, low-
be seen in ϵ15 and ϵ50 models at large distance from the density gas displays an intermediate to high metallicity,
central galaxy as a high metallicity gas. and so was likely ejected from the galaxy by the AGN
We can see that the presence of an AGN feedback feedback. Therefore, the presence of an AGN serves to
increases the amount of hot, T > 104 K, low-density gas both heat and enrich the CGM, and to extend the reach
across the shown CGM, and that the density of this gas of the gaseous component in a galaxy.
decreases monotonically with increasing distance. We Figure 12 provides the CGM metallicity distribution
< 104 K, gas just outside
also see that there is more cold, ∼ within 2Rvir . The top row displays snapshots which
of 0.1Rvir . emphasize the increase in reach of the high metallicity
In the bottom row, the gas is colored by metallicity. region with ϵ. At the same time, the ϵ50 model reveals
This distribution displays a progressively larger density the change in shape of this region. It becomes elongated
4
range of the gas with T < ∼ 10 K with accretion effi- horizontally, approximately along the jet orientation.
ciency. The population of the cold gas near the galaxy, The bottom row displays the radial distribution of
∼ 0.1Rvir , is split into two parts. The portion of the gas metallicity, with the mass-weighted average superposed.
that sits near 0.1Rvir , is ∼ 104 K, and has a very low It confirms an increase in the mass-weighted average
metallicity must be either recently accreted or actively metallicity with AGN feedback. It also shows that the
accreting from the cosmological filaments. There is also central metallicity gradient becomes shallower with ϵ, a
< 4
∼ 10 K gas that seems to be quite high in metallicity, sign of increased radial mixing.
with Zgas /Z⊙ > 1, and sits just outside the starforming We summarize the halo baryonic content in Table 5,
disk, outside the HOP galaxy. This high metallicity gas which contains a fraction of the CGM. The stellar frac-
is present in both the ϵ15 and ϵ50 models, and serves as tion shows a decline with increasing efficiency, almost by
further evidence for the extension of the gaseous compo- a factor of 2. The gas fraction remains steady, except
nent of the galaxy being pushed out due to the presence for the ϵ50 model which drops significantly. Overall, the
of a stronger AGN feedback. Note that alongside this baryonic contents decrease with ϵ — a clear signature of
population is the cold gas at large distances, this sits the AGN feedback.
inside the substructure, which does not belong to the
16

Table 5. Halo Baryonic Content at z = 0

Model log Mhalo /M⊙ fh,∗ fh,gas fh,baryons

ϵ
ϵ0 11.92 0.13 0.04 0.17
ϵ5 11.91 0.13 0.04 0.17
ϵ15 11.88 0.08 0.04 0.12
ϵ50 11.87 0.07 0.03 0.10
Note—All quantities have been calculated within
the virial radius of the DM halo. In this case,
Mhalo is the total DM+baryonic mass within Rvir .
Columns: (1) Model name, (2) Mhalo mass, (3) fh,∗
fraction of stellar mass, (4) fh,gas fraction of the gas
mass, (5) fh,baryons fraction of baryons.

3.7. The jet cocoons and the ISM, CGM and IGM
The jet cocoons in our modeled galaxies propagate
well beyond the HOP-galaxies, and, depending on the
AGN feedback, even cover the region beyond the host
halo virial radius. We follow the definition of the cir-
cumgalactic medium (CGM) as contained within 4Rvir .
This distance corresponds roughly to the baryonic back-
splash radius, where some of the inflow misses the cen-
tral galaxy and is diverted out entirely because of gravi-
tational effects (Bi et al. 2024), see also Tumlinson et al.
(2017, and refs. therein). The complicated kinematics
of the CGM, which includes the filamentary and dif-
fuse inflow, as well as outflow generated by the back-
splash, and, in our case, also by the AGN-driven outflow,
provides a challenge to disentangle all these phenom-
ena. Beyond the CGM, the AGN outflow propagates
throughout the IGM.
The AGN jet can deposit its linear momentum as well
as thermal and kinetic energy when interacting with the
ambient gas. Our numerical simulations do not resolve
the jet’s head (the hot spot), where the jet energy is Figure 13. Radial velocities and temperature of expanding
cocoons in AGN and ϵ0 galaxies at z = 0. Top: slice pro-
decollimated. However, we do resolve the backflow of
jections of expanding cocoons into the CGM and IGM, with
the jet shocked particles and their interactions with the boundaries delineated by shocks. Each model is shown on
ambient gas — these result in the formation of an over- a different scale, indicated by the scale bars, to emphasize
pressured and expanding cocoon. In order to investigate the difference in the cocoon sizes. The color palette provides
the interaction between this cocoon and the ambient gas the temperature. The white circles are Rvir . Bottom: Jet
within our modeled Seyfert galaxies and consequently its cocoon expansion provided by radial velocity vs distance di-
effect on the SF, we looked at the dynamic and thermo- agrams for the modeled galaxies. The black arrows mark the
respective Rvir . The color palette provides the gas tempera-
dynamic properties of these cocoons.
ture. Large regions of hot, fast outflowing velocity gas show
In the ϵ0 model, only the SN feedback exists and af- radial propagation of the jets and jet cocoons. For ϵ15 and
fects the ISM and CGM within the virial radius by in- ϵ50 models, previous generations of expanding cocoons have
creasing mixing and temperature, as seen in the top moved beyond the plotted region, affecting the CGM and
frames of Figure 13. The ISM gas is heated, but is able to IGM (not shown here).
cool down, and only a small fraction has T ∼ 106 −107 K.
 17

initial speed exceeded 1000 km s−1 , followed by the next

generation cocoon which changed its direction due to a
change in the jet orientation. The gas temperature be-
hind the front of the cocoon is similar in all models. The
timeline of the breakout, associated shock front veloci-
ties, and temperatures for the cocoons are quoted here
for reference only. We reserve the details of the cocoon
evolution for Paper II.
The lower frames of Figure 13 show the radial veloci-
ties and temperatures of the gas particles in the cocoons
at z = 0. It is evident that all AGN models feature mul-
tiple bursts of hot, fast moving, expanding gas on a scale
of 1.5 Mpc that are absent from the ϵ0 model. The shock
fronts delineate hot, ∼ 106 − 107 K, cocoon fronts, i.e.,
Figure 14. Total (thermal+kinetic) energy deposition the shocks, that propagate away from the galaxy. The
rate into the environment by the jet particles averaged over older cocoons, while hotter than ambient gas, have basi-
the last ∼ 1 Gyr preceding z = 0. The radius is binned in cally stopped their expansion. By comparing the radial
spherical shells with the bin size of 0.3 kpc. Each jet particle
occurrences of the hot outflows with the bubbles, we
has been followed from the point of emission until it’s energy
decreases to 10% of its spawned value or it merges with the can see that they represent the same phenomena, which
surrounding gas, for each model with the SMBH. describe the jet cocoons.
The highest velocities occur during the initial expan-
By z = 0, the majority of the ISM has T < 5 sion of the cocoon ‘bubbles.’ Figure13 provides a snap-
∼ 10 K. The
effect of even weak AGN accretion efficiency, in the ϵ5 shot at z=0 where only the ϵ5 model is undergoing this
model, is more complex. The majority of the ISM still early breakout phase. These outflows have an effect on
remains at t < 5 the CGM and IGM, enrich these regions with metals,
∼ 10 , but the hot component is heated
up to T ∼ 10 − few × 107 K. The hottest gas within
6 and are expected to influence the accretion rate onto
the CGM has T ∼ 105 − 107 K, and has obviously been the central galaxy. Moreover,expanding cocoons should
affected by the interaction with the jet flow. Within the affect the ISM properties and the SF in the substructure
cocoon volume, most of the CGM has been heated up surrounding the parent DM halos.
by a factor of ∼ 2 − 3 compared to the ϵ0 model. While the cocoons reach a large distance from the
The cocoon in ϵ5 model has stagnated at ∼ 750 kpc, galaxy, the majority of the hyper-refined jet particles
so to ∼ 3Rvir . But the breakout has occurred only very typically propagate out to much smaller distances, sub-
late, at z ∼ 0.2. The average velocity of the cocoon ject to hydrodynamical interactions with the surround-
expansion is ∼ 260 km s−1 , and the cocoon extends per- ing gas. The initial expansion of the cocoons is clearly
pendicularly to the major filament where the galaxy re- driven by their over-pressure. Ram pressure can drive
sides, after the breakout. The velocity is clearly super- fast shocks which compress and sweep up the ambient
sonic with respect to the ambient gas, and hence the gas. At large distances, the cocoon fronts appear to
cocoon expansion is delineated by a shock front, which move inertially.
are observable in the top frames of Figure 13. The co- The mechanical energy output by the jets at z = 0
coons are pushed by internal pressure which is supported is typical to that observed in Seyfert jets (e.g., Whittle
by T ∼ 106 − 107 K. et al. 2004), namely, Ljet ∼ 2 × 1041 erg s−1 in ϵ5 , ∼
In the ϵ15 and ϵ50 models, we detect multiple nested 8 × 1040 erg s−1 in ϵ15 , and ∼ 2 × 1040 erg s−1 in ϵ50 .
cocoons — meaning that the outflow follows the same This is about ∼ 10−5−6 of the Eddington mechanical
channel established by the preceding cocoon. In some luminosity of the modeled SMBHs.
cases the subsequent cocoon follows a different path, but We have analyzed the distribution of the deposited en-
after the breakout it always expands nearly perpendicu- ergy by the jet particles as a function of distance from
larly to the main filament. In the ϵ15 model, the cocoon the SMBH, shown in Figure 14. It shows the total (ther-
has reached 2 Mpc at z = 0 in its longest elongation. mal and kinetic) energy deposition rate into the ambient
It has decelerated and stagnated at this distance. The gas by the jet as a function of R averaged over the last
average velocity of this cocoon is about 500 km s−1 . 1 Gyr. The energy deposition rate is found from the
The cocoon in the most energetic AGN ϵ50 model has total change in energy for each jet particle as it propa-
reached ∼ 2.4 Mpc at z ∼ 0.8, where it stagnated. Its gates away from the SMBH until it re-merges with the
18

gas. The changes in energy for all jet particles spawned tum of jet particles emitted anisotropically along the
within the 1 Gyr time period are summed and normal- spin axis of the SMBH. Galaxies at z = 0 differ in the
ized by the radial bin size and time. We find that nearly following ways with increasing accretion efficiency, pre-
all of the thermal energy from the jet is injected when sented here along the ϵ sequence, from ϵ = 0 to ϵ = 0.5:
the jet particles merge with the surrounding gas. The
jet particles barely escape beyond ∼ 0.1Rvir , and only in • Stellar masses in modeled galaxies within DM ha-
ϵ50 model. The ϵ50 model shows that the deposited en- los of ∼ 6.3 × 1011 M⊙ lie within the range of
ergy rate per unit radius is relatively flat out to ∼22 kpc, 3 × 1010 − 1011 M⊙ — the AGN galaxies have
apart from a slight decrease in the inner few kpc. For progressively smaller masses with increased ϵ com-
the smaller efficiency models, dL/dR is a little steeper. pared to the ϵ0 galaxy. Both the stellar and gas
Note that most of the energy has been deposited in fractions within the parent DM halos decrease
the central kpc at low z by the ϵ5 and ϵ15 jets, while not along the ϵ sequence, from 12.9% to 7.2% (stellar)
by the ϵ50 jet, which deposits its energy more uniformly and 4.1% to 2.8% (gas), respectively. Hence, the
with R. This difference is apparently related to the gas baryonic fraction inside the virial radius decreases
distribution in modeled galaxies, where ϵ50 galaxy has from 17% down to 10%.
a large central cavity devoid of gas at this time. The
deceleration of the jet close to the production region has • The radial distribution of the gaseous component
been suggested as the origin of the FR-I morphology in differs from model to model. Some of these differ-
powerful radio sources by De Young (1993). ences can be understood following the cosmologi-
cal evolution, others are a clear consequence of the
AGN feedback. Increase in the size of the gaseous
4. DISCUSSION AND CONCLUSIONS disk with increasing efficiency is the most obvious
We used high-resolution zoom-in cosmological simula- of such trends. The ϵ50 galaxy has nearly all the
tions to model evolution of three Seyfert-type galaxies, gas residing at the galaxy outskirts.
with the SMBHs seeded at z = 3.7 and triggered colli-
mated jets. We find that the SMBH jets in this energy • Sersic decomposition of the stellar distribution
range exert an effect on galactic morphology, their SFRs, shows the presence of a massive bulge in all galax-
their gas fraction, gas distribution, etc. The effects of ies, but B/T decreases ∼1/2 along the ϵ sequence,
this feedback, depending on the associated accretion ef- from ∼ 0.55 to ∼ 0.28. The stellar disk mass and
ficiency, extend outside the galaxy, and influence the its scalelength decrease by a factor of ∼ 1.5 and
CGM properties even beyond Rvir , up to the baryonic ∼ 2, respectively, with ϵ.
backsplash radii, determining the mass accretion onto
• The SFR is inversely proportional to the accretion
the parent DM halo. The evolution of the AGN galax-
efficiency, which decreases the SFR by three orders
ies has been compared to a galaxy without the SMBH.
of magnitude, from 2 M⊙ yr−1 to 2×10−3 M⊙ yr−1 .
We have analyzed the interactions between the SMBH
The ϵ50 model has essentially a quenched SF after
jets, their host galaxies and immediate environment, the
z ∼ 1. In all AGN models, stars forming after
ISM and the CGM. In this work we have only focused
z ∼ 1 have their origin outside the central few
on the objects at z = 0, while cosmological evolution for
kpc. This central gap increases with ϵ. The SFR
z< ∼ 10 will be presented in Paper II. in ϵ0 and ϵ5 models appear comparable, while ϵ15
Feedback from jets launched by the SMBHs deposit
is more similar to the ϵ50 model. Overall, the SF
both energy and momentum in the ambient gas. This
is shifted to larger radii with increasing ϵ.
process can have a potentially far reaching effect on the
gas distribution in the host galaxy, its thermal and kine- • The ISM metallicity is about 5 - 10 times smaller
matic properties, and consequently on the SF. The inte- than the stellar metallicity, meaning that stars
gral effect on the evolution of Seyfert galaxies depends formed from a more metal-rich gas which has
on a number of model parameters characterizing the jets. been subsequently replaced by the inflowing low-
We choose to perform a controlled experiment by vary- metallicity gas. The ϵ0 and ϵ5 models show a radial
ing only a single parameter at a time. Here we present gradient in the average ISM metallicity, while ϵ15
a suite of models which differ by Eddington fraction of and ϵ50 models are radially mixed.
the SMBH accretion rate. The SMBHs of M• ≃ 106 M⊙
have been seeded relatively late, when Mvir ∼ 1011 M⊙ . • The CGM properties differ along the ϵ sequence.
4
The energy output by a seeded SMBH has been chosen The cold, < ∼ 10 K gas displays a progressively
as the thermal energy, kinetic energy, and linear momen- larger density range with ϵ, while hotter gas shows
 19

a lower metallicity and a higher density range with Studies of CO emissions or kinematics of atomic and
increasing ϵ. The mass-weighted average metallic- molecular gas in jetted low-luminosity AGN have shown
ity of the CGM gas is higher in ϵ15 and ϵ50 models that the gas is often displaced from the nuclear region,
by a factor of a few. forming cavities, or rings, that can reach radii of >
∼ 6 kpc
or clouds of displaced gas (Garcı́a-Burillo et al. 2024;
• All AGN models develop AGN-driven outflows, Nesvadba et al. 2021; Nandi et al. 2023; Matsushita
which define expanding cocoons (bubbles) out- et al. 2007). These rings have been shown to be asso-
lined by shocks — extending to ∼ 750 kpc in the ciated with either decreased SF (Nesvadba et al. 2021;
ϵ5 , to ∼ 2 Mpc in ϵ15 , and to ∼ 2.4 Mpc in ϵ50 Rao et al. 2023; Nandi et al. 2023) or localized regions
model. The cocoons are pressure-driven. The of enhanced SF (Duggal et al. 2024; Venturi et al. 2021;
comparison ϵ0 model is subject only to the SN Krause 2023). Gaseous rings can be associated with
feedback which mixes the gas within Rvir , and jets interacting with the MHD winds associated with
no cocoons have developed in this case. The co- the BEL regions and molecular tori (Blandford & Payne
coon expansion in AGN galaxies is supersonic with 1982; Emmering et al. 1992; Konigl & Kartje 1994; Bot-
respect to the ambient gas, with ∼ 260 km s−1 , torff et al. 1997; Elitzur & Shlosman 2006) or with the
∼ 500 km s−1 , and ∼ 1, 000 km s−1 , with increas- gas elevated above the disk plane by SN feedback, and
ing ϵ. The radial velocities of the outflowing gas are not required to be directed into the disk plane. Also,
within the cocoons are much higher of course, the disk plane can be warped, as in our modeled galax-
reaching few×103 km s−1 . ies. Our results confirm these observations as at z=0,
but our jets do not impact the disks directly. We observe
• The cocoons observed in AGN galaxies have oval regions with damped SF around the SMBHs, formation
shapes, with axial ratios of 1:2 – 1:3. The inter- of central cavities and extended gaseous and stellar disks
nal temperatures of the cocoons lie in the range of with increasing accretion efficiency (e.g., Figs. 3 — 6).
few×105 − few × 107 K, where the upper limit de- Looking into the SF distribution, Acharya et al.
lineates the front shocks in the expanding cocoons. (2024) show that disk-dominated and pseudo-bulge
AGN-hosting galaxies tend to have a decrease in SF
4.1. Comparison with observations in the central region and enhanced SF at larger radii,
leading to an overall SF matching that of similar star-
We now turn to comparison with observations and nu-
forming galaxies. This result is further supported by
merical modeling, and quote references in addition to
other observations (Kurian et al. 2024; Lammers et al.
those discussed in Section 1. Observational data on jets
2023; Krause 2023). Bulge dominated systems have a
and their feedback in Seyfert galaxies is scarce, but in-
flat SF distribution across galactic radii. The central
terest is growing. We are still lacking multi-wavelength,
gap which is mostly devoid of SF after z ∼ 1 is observed
high-resolution statistical surveys which can address the
in all of our AGN galaxies, and clearly correlates with
effects that these pc- and kpc-scale jets have on their
the strength of the AGN accretion efficiency.
host galaxies, which are disk-type, unlike their more
A review by Krause (2023) emphasizes how CGM den-
powerful radio galaxies counterparts. Currently, the
sity and virial temperature can strongly influence the
comparison is mostly limited to individual objects (e.g.,
ability of the jet to propagate outward and whether one
Varglund et al. 2022).
finds compression and associated triggered SF, or heat-
Observations of jetted Seyfert galaxies show that they
ing and gas expulsion resulting in quenching of SF.
typically reside in disk galaxies with spiral morphologies
Increased radio emission strength, along the sequence
and host SMBHs in the mass range of 106 − 108 M⊙
from star-forming galaxies to composite (i.e., galaxies
(e.g., Vietri et al. 2024) — our modeled galaxies are
with radio emission from SF and jets) to Seyfert and
bulged disks and their SMBHs fall within this range (see
LINERs, has been claimed in Vitale et al. (2015). This
Table 4).
increase appears to correlate with an analogous decrease
There is abundant evidence that Seyfert jets generate
in Hα emission, which hints that jets are indeed instru-
bow shocks, through interactions with the ISM, that
mental in halting SF not only in large elliptical galaxies,
propagate out with velocities up to 103 km s−1 (e.g., Ce-
but also lower mass disk galaxies.
cil et al. 2000; Tadhunter et al. 2014; Murthy et al. 2022;
Varglund et al. (2023) has studied colors and mor-
Peralta de Arriba et al. 2023; Duggal et al. 2024). Mul-
phologies of a sample of NLS1 galaxies and found disky
tiple generations of misaligned jet radio lobes have been
morphologies and increased dust extinction in the cen-
observed, indicating multiple episodes of jet activity (Se-
ters. The latter feature is often seen in γ-ray-detected
bastian et al. 2019; Rao et al. 2023).
20

jetted NLS1 galaxies. This reddening has been corre- giving an altogether smaller bulge component and a
lated with the jet activity in other observations (e.g., ‘fluffier’ disk. Adding the ‘radio’ mode does not lead to
Olguı́n-Iglesias et al. 2020; Berton et al. 2019). the same outcome (Irodotou et al. 2022). Kinematic de-
Finally, Cunow (2004) observed an increased radial composition finds that the strongest and combined feed-
mixing of metals in stellar populations of Seyfert galax- back modes lead to dissolution of the disk component,
ies in comparison to non-active galaxies, which show en- but that intermediate feedback shows little change from
hanced metallicity in stars at the center. Possibly also the fiducial models (Byrne et al. 2024; Okamoto et al.
a reduced age gradient for Seyferts. 2008). While not quantitatively explored in the paper,
Wellons et al. (2023) mentioned that various AGN mod-
4.2. Comparison with simulations els tested do lead to differences in the host morphology.
Exploration of long-term effects to metallicity and
Numerical work on the effects of AGN jet feedback in
mixing in simulations of Seyfert galaxies is lacking.
Seyfert-type galaxies over cosmological timescales is cur-
While we do not model the metal diffusion explicitly,
rently limited, perhaps in anticipation that AGN feed-
we observe flatter final radial metallicity profiles with in-
back has little effect on MW-mass and smaller galaxies.
creased accretion efficiency, as well as slightly increased
Building up a statistical sample of observational and
metallicity in parent halos. This effect is also noticed
theoretical works is essential to understanding the full
by Appleby et al. (2021), they observed an increase in
impact of low-luminosity AGN jets on their hosts and
the metallicity of hot CGM gas in full-box SIMBA sim-
for lying groundwork for high-z observations of these
ulations which include jet feedback. Talbot et al. (2024)
galaxies.
similarly find increased metallicity in the CGM due to
The aim of our simulation set has been to evaluate
the presence of jets, and so did Qutob et al. (2023), but
the impact of varying efficiency of the SMBH accretion
their jet feedback included cosmic rays and radiation.
in Seyfert galaxies. Isolating the effect of this param-
Modeling the evolution of jet cocoons in Seyfert
eter in a low-luminosity AGN model along with colli-
galaxies under the influence of variation in accretion
mated jet feedback has not been addressed so far in the
efficiency, and specifically the associated timescales,
literature for cosmological simulations. Numerical sim-
shapes, shocks, internal velocities and temperature, has
ulations have demonstrated the ability of AGN jets to
not been addressed in the literature so far. However,
halt the cooling flows in massive galaxies and galaxy
short-term studies about launching and propagation of
clusters. Some works point that jets act to reduce the
cocoons from low-luminosity AGN do exist. Talbot et al.
stellar mass in MW-type and smaller host galaxies as
(2024) found comparable maximal velocities and tem-
well (Wellons et al. 2023; Byrne et al. 2024; Irodotou
peratures inside the cocoon, but for higher luminos-
et al. 2022; Okamoto et al. 2008). However, short-term
ity jets in merging galaxies. Qutob et al. (2023) ob-
studies have noted brief or localized increases in the SFR
served cocoon properties scaling with increasing a multi-
in MW-mass, or Seyfert-type, galaxies (Mukherjee et al.
component feeback. Tanner & Weaver (2022) studied
2018; Feliz & Angles-Alcazar 2023).
the effect of varying the inclination angle of a low-
We observe that SF is pushed to large radii, due to
luminosity jet with respect to the galactic disk.
a more extended gaseous disk with an increased accre-
In summary, similar trends and effects which show
tion efficiency. Irodotou et al. (2022) observed the same
up in our simulations of Seyfert galaxies at z = 0 have
trend but using a ‘radio’ mode feedback, by dumping
been observed in nearby Seyfert galaxies and simulation
thermal energy in the CGM. Byrne et al. (2024) shows
studies. Although statistical analysis of these properties
that central gas surface density decreases and stellar
is currently absent, our results point to the need for a
half-mass radii tend to grow by adding additional feed-
systematic study of feedback from low luminosity AGN,
back in the form of jets, radiation, and cosmic rays.
like Seyfert galaxies. In our study, we have limited the
Formation of central gaseous voids, as observed in our
parameter space to modifying only the efficiency of the
highest efficiency model, ϵ50 , is also reported by Qutob
SMBH accretion which has a profound effect on the ISM
et al. (2023). Short-term (non-cosmological) studies of
of the host galaxy, whose properties will determine the
jets, which vary their orientation with respect to the
morphology of the underlying stellar disk. In the sub-
galaxy, show that the jet can displace gas and dissolve
sequent paper II, we analyze evolution of Seyfert-type
cold clouds in the central regions of galaxies (Mukherjee
galaxies for z < ∼ 10, and focus on effects of the SMBH
et al. 2018; Talbot et al. 2022).
seeding redshift.
Using the Sersic decomposition, we show that increas-
ing SMBH accretion efficiency leads to a decrease in
B/D, B/T, in disk and bulge masses and scale length,
 21

5. ACKNOWLEDGMENTS ter 1601 (SFB 1601 sub-project C5), both funded by

We thank Phil Hopkins for providing us with the lat- the DFG–500700252. This work used Expanse CPU at
est version of the code, and to Alessandro Lupi, Kung-Yi San Diego Supercomputer Center (SDSC) through allo-
Su, and Paul Torrey for their help with GIZMO. Many cation PHY230135 from the ACCESS program, which
thanks to Da Bi and Xingchen Li for invaluable help is supported by NSF grants 2138259, 2138286, 2138307,
with many aspects of the work carried out here. I.S. is 2137603, and 2138296 (Boerner et al. 2023), and by the
grateful for a generous support from the International University of Kentucky Lipscomb Computing Cluster.
Joint Research Promotion Program at Osaka Univer- We thank Vikram Gazula at the Center for Computa-
sity, and acknowledges the hospitality of KITP where tional Studies (UK) for a very generous computational
part of this research has been conducted. The KITP is support throughout this project.
supported by NSF PHY-1748958. E.R.D. acknowledges
DATA AVAILABILITY
support of the Collaborative Research Center 956, sub-
project C4; and from the Collaborative Research Cen- The data used for this paper can be available upon
reasonable request.

REFERENCES
Acharya, N., Bonoli, S., Salvato, M., et al. 2024, A&A, 687, Blandford, R. D., & Rees, M. J. 1974, MNRAS, 169, 395,
A285, doi: 10.1051/0004-6361/202449287 doi: 10.1093/mnras/169.3.395
Anglés-Alcázar, D., Davé, R., Faucher-Giguère, C.-A., Özel, Boerner, T. J., Deems, S., Furlani, T. R., Knuth, S. L., &
F., & Hopkins, P. F. 2017, MNRAS, 464, 2840, Towns, J. 2023, ACCESS: Advancing Innovation: NSF’s
doi: 10.1093/mnras/stw2565 Advanced Cyberinfrastructure Coordination Ecosystem:
Appleby, S., Davé, R., Sorini, D., Storey-Fisher, K., & Services & Support, PEARC ’23 (New York, NY, USA:
Smith, B. 2021, MNRAS, 507, 2383, Association for Computing Machinery), 173,
doi: 10.1093/mnras/stab2310 doi: 10.1145/3569951.3597559
Arav, N., Shlosman, I., & Weymann, R. 1997, in ASP Conf. Bottorff, M., Korista, K. T., Shlosman, I., & Blandford,
Ser., Vol. 128, Mass Ejection from Active Galactic Nuclei R. D. 1997, ApJ, 479, 200, doi: 10.1086/303867
Bagchi, J., Vivek, M., Vikram, V., et al. 2014, ApJ, 788, Bourne, M. A., & Yang, H.-Y. K. 2023, Galaxies, 11, 73,
174, doi: 10.1088/0004-637X/788/2/174 doi: 10.3390/galaxies11030073
Baldi, R. D., Williams, D. R. A., McHardy, I. M., et al. Byrne, L., Faucher-Giguère, C.-A., Wellons, S., et al. 2024,
2018, MNRAS, 476, 3478, doi: 10.1093/mnras/sty342 ApJ, 973, 149, doi: 10.3847/1538-4357/ad67ca
Begelman, M. C., Blandford, R. D., & Rees, M. J. 1984, Cecil, G., Greenhill, L. J., DePree, C. G., et al. 2000, ApJ,
Reviews of Modern Physics, 56, 255, 536, 675, doi: 10.1086/308959
doi: 10.1103/RevModPhys.56.255 Chabrier, G. 2003, ApJL, 586, L133, doi: 10.1086/374879
Begelman, M. C., & Cioffi, D. F. 1989, ApJL, 345, L21, Chen, H.-Z., Kang, X., Macciò, A. V., Buck, T., & Cen, R.
doi: 10.1086/185542 2023, arXiv e-prints, arXiv:2310.13069,
Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. doi: 10.48550/arXiv.2310.13069
2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182 Contopoulos, J. 1995, ApJ, 450, 616, doi: 10.1086/176170
Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2012, ApJ, Cunow, B. 2004, MNRAS, 353, 477,
762, 109, doi: 10.1088/0004-637x/762/2/109 doi: 10.1111/j.1365-2966.2004.08082.x
Berton, M., Congiu, E., Ciroi, S., et al. 2019, AJ, 157, 48, Damas-Segovia, A., Beck, R., Vollmer, B., et al. 2016, ApJ,
doi: 10.3847/1538-3881/aaf5ca 824, 30, doi: 10.3847/0004-637X/824/1/30
Berton, M., Järvelä, E., Crepaldi, L., et al. 2020, A&A, de Kool, M., & Begelman, M. C. 1995, ApJ, 455, 448,
636, A64, doi: 10.1051/0004-6361/202037793 doi: 10.1086/176594
Bi, D., Shlosman, I., & Romano-Dı́az, E. 2022a, MNRAS, De Young, D. S. 1993, ApJL, 405, L13, doi: 10.1086/186754
513, 693, doi: 10.1093/mnras/stac363 Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature,
—. 2022b, ApJ, 934, 52, doi: 10.3847/1538-4357/ac779b 433, 604–607, doi: 10.1038/nature03335
—. 2024, MNRAS, 527, 11095, Duggal, C., O’Dea, C. P., Baum, S. A., et al. 2024, ApJ,
doi: 10.1093/mnras/stad3942 965, 17, doi: 10.3847/1538-4357/ad2513
Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883, Dutta, R., Sharma, P., Sarkar, K. C., & Stone, J. M. 2024,
doi: 10.1093/mnras/199.4.883 ApJ, 973, 148, doi: 10.3847/1538-4357/ad67d7
22

Eisenstein, D. J., & Hut, P. 1998, ApJ, 498, 137, Irodotou, D., Fragkoudi, F., Pakmor, R., et al. 2022,
doi: 10.1086/305535 MNRAS, 513, 3768, doi: 10.1093/mnras/stac1143
Elitzur, M., & Shlosman, I. 2006, ApJL, 648, L101, Kauffmann, G., & Haehnelt, M. 2000, MNRAS, 311, 576,
doi: 10.1086/508158 doi: 10.1046/j.1365-8711.2000.03077.x
Emmering, R. T., Blandford, R. D., & Shlosman, I. 1992, Kauffmann, G., Heckman, T. M., White, S. D. M., et al.
ApJ, 385, 460, doi: 10.1086/170955 2003, MNRAS, 341, 54,
Falcke, H., Nagar, N. M., Wilson, A. S., & Ulvestad, J. S. doi: 10.1046/j.1365-8711.2003.06292.x
2000, ApJ, 542, 197, doi: 10.1086/309543 Keel, W. C., White, Raymond E., I., Owen, F. N., &
Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614, Ledlow, M. J. 2006, AJ, 132, 2233, doi: 10.1086/508340
doi: 10.1093/mnras/staa302 Komissarov, S., & Porth, O. 2021, NewAR, 92, 101610,
Feliz, J. M., & Angles-Alcazar, D. 2023, Bulletin of the doi: 10.1016/j.newar.2021.101610
AAS, 55 Komossa, S. 2018, in Revisiting Narrow-Line Seyfert 1
Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9, Galaxies and their Place in the Universe, 15,
doi: 10.1086/312838 doi: 10.22323/1.328.0015
Fukumura, K., Tombesi, F., Kazanas, D., et al. 2014, ApJ, Komossa, S., Voges, W., Xu, D., et al. 2006, AJ, 132, 531,
780, 120, doi: 10.1088/0004-637X/780/2/120 doi: 10.1086/505043
Gallimore, J. F., Axon, D. J., O’Dea, C. P., Baum, S. A., & Konigl, A., & Kartje, J. F. 1994, ApJ, 434, 446,
Pedlar, A. 2006, AJ, 132, 546, doi: 10.1086/504593 doi: 10.1086/174746
Garcı́a-Burillo, S., Hicks, E. K. S., Alonso-Herrero, A.,
Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511,
et al. 2024, A&A, 689, A347,
doi: 10.1146/annurev-astro-082708-101811
doi: 10.1051/0004-6361/202450268
Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581,
Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJL,
doi: 10.1146/annurev.aa.33.090195.003053
539, L13, doi: 10.1086/312840
Krause, M. G. H. 2023, Galaxies, 11,
Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010,
doi: 10.3390/galaxies11010029
MNRAS, 404, 1111,
Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36,
doi: 10.1111/j.1365-2966.2010.16341.x
doi: 10.1088/0004-637X/729/1/36
Hada, K., Doi, A., Nagai, H., et al. 2013, ApJ, 779, 6,
Kurian, K. S., Stalin, C. S., Rakshit, S., et al. 2024, A&A,
doi: 10.1088/0004-637X/779/1/6
688, A32, doi: 10.1051/0004-6361/202346246
Hahn, O., & Abel, T. 2011, MNRAS, 415, 2101,
Laine, S., & Beck, R. 2008, ApJ, 673, 128,
doi: 10.1111/j.1365-2966.2011.18820.x
doi: 10.1086/523960
Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589,
Lammers, C., Iyer, K. G., Ibarra-Medel, H., et al. 2023,
doi: 10.1146/annurev-astro-081913-035722
ApJ, 953, 26, doi: 10.3847/1538-4357/acdd57
Heller, C. H., & Shlosman, I. 1994, ApJ, 424, 84,
doi: 10.1086/173874 Ledlow, M. J., Owen, F. N., Yun, M. S., & Hill, J. M. 2001,
Ho, L. C., & Ulvestad, J. S. 2001, ApJS, 133, 77, ApJ, 552, 120, doi: 10.1086/320458
doi: 10.1086/319185 Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ,
Hopkins, P. F. 2015, MNRAS, 450, 53, 115, 2285, doi: 10.1086/300353
doi: 10.1093/mnras/stv195 Mao, M. Y., Owen, F., Duffin, R., et al. 2015, MNRAS,
Hopkins, P. F., Narayanan, D., & Murray, N. 2013, 446, 4176, doi: 10.1093/mnras/stu2302
MNRAS, 432, 2647, doi: 10.1093/mnras/stt723 Martı́, J.-M. 2019, Galaxies, 7, 24,
Hopkins, P. F., & Quataert, E. 2011, MNRAS, 415, 1027, doi: 10.3390/galaxies7010024
doi: 10.1111/j.1365-2966.2011.18542.x Martinez-Valpuesta, I., Shlosman, I., & Heller, C. 2006,
Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, ApJ, 637, 214, doi: 10.1086/498338
480, 800, doi: 10.1093/mnras/sty1690 Matsushita, S., Muller, S., & Lim, J. 2007, A&A, 468, L49,
Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2022, doi: 10.1051/0004-6361:20067039
MNRAS, 519, 3154, doi: 10.1093/mnras/stac3489 McNamara, B. R., & Nulsen, P. E. J. 2012, New Journal of
Hota, A., Sirothia, S. K., Ohyama, Y., et al. 2011, MNRAS, Physics, 14, 055023, doi: 10.1088/1367-2630/14/5/055023
417, L36, doi: 10.1111/j.1745-3933.2011.01115.x Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland,
Hudson, M. J., Harris, G. L., & Harris, W. E. 2014, ApJL, R. S., & Silk, J. 2018, MNRAS, 479, 5544,
787, L5, doi: 10.1088/2041-8205/787/1/L5 doi: 10.1093/mnras/sty1776
 23

Murray, N., Chiang, J., Grossman, S. A., & Voit, G. M. Sebastian, B., Kharb, P., O’Dea, C. P., Gallimore, J. F., &
1995, ApJ, 451, 498, doi: 10.1086/176238 Baum, S. A. 2019, MNRAS, 490, L26,
Murthy, S., Morganti, R., Wagner, A. Y., et al. 2022, doi: 10.1093/mnrasl/slz136
Nature Astronomy, 6, 488, Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al.
doi: 10.1038/s41550-021-01596-6 2020, MNRAS, 495, 3252, doi: 10.1093/mnras/staa1381
Nandi, P., Stalin, C. S., Saikia, D. J., et al. 2023, ApJ, 959, Shlosman, I. 2013, in Secular Evolution of Galaxies, ed.
116, doi: 10.3847/1538-4357/ad0c57 J. Falcón-Barroso & J. H. Knapen, 555,
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, doi: 10.48550/arXiv.1212.1463
462, 563, doi: 10.1086/177173 Shlosman, I., Frank, J., & Begelman, M. C. 1989, Nature,
Nesvadba, N. P. H., Wagner, A. Y., Mukherjee, D., et al. 338, 45, doi: 10.1038/338045a0
2021, A&A, 654, A8, doi: 10.1051/0004-6361/202140544 Shlosman, I., Vitello, P. A., & Shaviv, G. 1985, ApJ, 294,
Okamoto, T., Nemmen, R. S., & Bower, R. G. 2008, 96, doi: 10.1086/163278
MNRAS, 385, 161, doi: 10.1111/j.1365-2966.2008.12883.x Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289,
Olguı́n-Iglesias, A., Kotilainen, J., & Chavushyan, V. 2020, doi: 10.1046/j.1365-8711.2003.06206.x
MNRAS, 492, 1450, doi: 10.1093/mnras/stz3549 Stewart, M., & Federrath, C. 2022, MNRAS, 509, 5237,
Osterbrock, D. E., & Shuder, J. M. 1982, ApJS, 49, 149, doi: 10.1093/mnras/stab3313
doi: 10.1086/190793 Su, K.-Y., Hopkins, P. F., Bryan, G. L., et al. 2021,
MNRAS, 507, 175, doi: 10.1093/mnras/stab2021
Padoan, P., Haugbølle, T., & Nordlund, Å. 2012, ApJL,
Tadhunter, C., Morganti, R., Rose, M., Oonk, J. B. R., &
759, L27, doi: 10.1088/2041-8205/759/2/L27
Oosterloo, T. 2014, Nature, 511, 440,
Peralta de Arriba, L., Alonso-Herrero, A., Garcı́a-Burillo,
doi: 10.1038/nature13520
S., et al. 2023, A&A, 675, A58,
Talbot, R. Y., Sijacki, D., & Bourne, M. A. 2022, MNRAS,
doi: 10.1051/0004-6361/202245408
514, 4535, doi: 10.1093/mnras/stac1566
Piotrowska, J. M., Bluck, A. F. L., Maiolino, R., & Peng, Y.
Talbot, R. Y., Sijacki, D., & Bourne, M. A. 2024, Monthly
2021, MNRAS, 512, 1052, doi: 10.1093/mnras/stab3673
Notices of the Royal Astronomical Society, 528, 5432,
Planck Collaboration, Ade, P. A. R., Aghanim, N., et al.
doi: 10.1093/mnras/stae392
2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
Tanner, R., & Weaver, K. A. 2022, The Astronomical
Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS,
Journal, 163, 134, doi: 10.3847/1538-3881/ac4d23
519, 1526, doi: 10.1093/mnras/stac3214
Torrey, P., Hopkins, P. F., Faucher-Giguère, C.-A., et al.
Qutob, N., Emami, R., Su, K.-Y., et al. 2023, arXiv
2020, MNRAS, 497, 5292, doi: 10.1093/mnras/staa2222
e-prints, arXiv:2312.14809,
Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017,
doi: 10.48550/arXiv.2312.14809
ARA&A, 55, 389,
Rao, V. V., Kharb, P., Rubinur, K., et al. 2023, MNRAS,
doi: 10.1146/annurev-astro-091916-055240
524, 1615, doi: 10.1093/mnras/stad1901
Varglund, I., Järvelä, E., Ciroi, S., et al. 2023, A&A, 679,
Rennehan, D., Babul, A., Moa, B., & Davé, R. 2024, A32, doi: 10.1051/0004-6361/202346131
MNRAS, 532, 4793, doi: 10.1093/mnras/stae1785 Varglund, I., Järvelä, E., Lähteenmäki, A., et al. 2022,
Romano-Dı́az, E., Shlosman, I., Choi, J.-H., & Sadoun, R. A&A, 668, A91, doi: 10.1051/0004-6361/202244465
2014, ApJL, 790, L32, doi: 10.1088/2041-8205/790/2/L32 Venturi, G., Cresci, G., Marconi, A., et al. 2021, A&A, 648,
Rosario, D. J., Whittle, M., Nelson, C. H., & Wilson, A. S. A17, doi: 10.1051/0004-6361/202039869
2010, MNRAS, 408, 565, Vietri, A., Berton, M., Järvelä, E., et al. 2024, A&A, 689,
doi: 10.1111/j.1365-2966.2010.17153.x A123, doi: 10.1051/0004-6361/202449247
Saikia, D. J., & Jamrozy, M. 2009, Bulletin of the Vitale, M., Fuhrmann, L., Garcı́a-Marı́n, M., et al. 2015,
Astronomical Society of India, 37, 63, A&A, 573, A93, doi: 10.1051/0004-6361/201423993
doi: 10.48550/arXiv.1002.1841 Wellons, S., Faucher-Giguère, C.-A., Hopkins, P. F., et al.
Salim, S. 2014, Serbian Astronomical Journal, 189, 1, 2023, MNRAS, 520, 5394, doi: 10.1093/mnras/stad511
doi: 10.2298/SAJ1489001S Whittle, M., Silverman, J. D., Rosario, D. J., Wilson, A. S.,
Scharré, L., Sorini, D., & Davé, R. 2024, MNRAS, 534, 361, & Nelson, C. H. 2004, in IAU Symp., Vol. 222, The
doi: 10.1093/mnras/stae2098 Interplay Among BHs, Stars and ISM in Galactic Nuclei,
Schmidt, M., & Green, R. F. 1983, ApJ, 269, 352, ed. T. Storchi-Bergmann, L. C. Ho, & H. R. Schmitt,
doi: 10.1086/161048 299, doi: 10.1017/S1743921304002315
24

Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009,

MNRAS, 393, 99, doi: 10.1111/j.1365-2966.2008.14191.x
 25

APPENDIX

A. SERSIC DECOMPOSITION
The bulge and disk components of the galaxy have been determined by fitting the 1-D face-on stellar surface density
profile of the stellar component of each galaxy within 0.1Rvir , and limited to ±5 kpc in z, to a model that combines
the Sersic function for the bulge, and two exponential disks, namely, the inner and outer disks:
1/n
Σ(R) = Σe e−bn [(R/Re ) −1]
+ Σ01 e−R/Rdisk1 + Σ02 e−R/Rdisk2 , (A1)
where Σe and Re are the effective surface density and radius of the bulge, and n is the bulge Sersic index. The best
fit has been achieved using a double-exponential disk. The fitting parameters are the disk’s central surface densities
Σ01,02 ,and their scalelengths Rdisk1,disk2 .

Figure 15. The Sersic decomposition of the stellar component in modeled galaxies within 0.1Rvir . The total profile fits to the
face-on stellar surface density are shown for each model. The top left is the ϵ0 model, ϵ5 the top right, ϵ15 on the bottom left,
and ϵ50 on the bottom right.

The numerical values of all the parameters of the decomposition are given in Table 3.