# The effect of local photoionization on the galaxy properties and the
circumgalactic medium in simulations of Milky Way-sized galaxies

Bocheng Zhu1,2 and Volker Springel2
1Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical
Observatory, Chinese Academy of Sciences,
80 Nandan Road, Shanghai 200030, People’s Republic of China
2Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741
Garching, Germany<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this study, we investigate the impact of local stellar radiation in
cosmological zoom simulations of the formation of Milky Way-sized galaxies. We
include the radiation field as an additional feedback component that is
computed alongside gravity with a tree code in an optically thin
approximation. We resimulate the initial conditions of five Milk Way-like
systems taken from the Auriga project with and without stellar radiation, and
study the effects of local stellar radiation on several properties of the
galaxies and the circumgalactic medium (CGM). Similar to previous findings, we
observe that local stellar radiation can modify gas cooling in the
circumgalactic medium and thus suppress star formation and the surface
densities of young stars and HI gas, while having little impact on the total
gas content. Even though the magnitude of the effect we find is smaller than
reported in earlier work, the radiation field is thus clearly an important
component in simulations of late time galaxies. In particular, it also
suppresses the peak of the rotation curve and reduces the mass of the stellar
bulge. In the CGM region, the young stellar radiation exceeds the external UVB
and dominates the radiation field within the virial halo at all redshifts.
Nevertheless, we find that the local stellar radiation has overall little
impact on the radial density and temperature profile of the CGM gas. The
metallicity profile is slightly reduced, however, as well as the HI and MgII
column densities within $\sim 0.3\,R_{\rm vir}$. In contrast, we find that the
OVI column density is hardly impacted by the radiation field. Additional
effects can be expected from the radiation of the central AGN during phases of
quasar activity, which has not yet been included in the simulations of the
present study.

###### keywords:

galaxies: formation – methods: numerical – galaxies: evolution – hydrodynamics

††pubyear: 2024††pagerange: The effect of local photoionization on the galaxy
properties and the circumgalactic medium in simulations of Milky Way-sized
galaxies–The effect of local photoionization on the galaxy properties and the
circumgalactic medium in simulations of Milky Way-sized galaxies

##  1 Introduction

In the classical theory of galaxy formation, gas follows the collapse of dark
matter (DM) due to the pull of gravity and reaches virial equilibrium through
shock-heating (Silk, 1977; Rees & Ostriker, 1977; White & Rees, 1978; White &
Frenk, 1991). However, unlike DM, gas is also subject to radiative cooling.
Upon dissipating its thermal support, it can further collapse into the center
of a DM halo, eventually becoming dense enough for star formation, or becoming
available to feed the growth of a supermassive black hole (SMBH). On the other
hand, SMBHs and stars will also inject energy due to the release of
gravitational potential energy and as a result of nuclear reactions,
respectively. The released photons can heat and photoionize the gas, thereby
in principal counteracting the cooling processes in the halos.

The injected radiation energy plays a significant role in galaxy formation and
evolution. The associated photons give rise to a background radiation field in
the universe, an important part of which comprises the so-called ultraviolet
background (UVB, Haardt & Madau, 1996; Faucher-Giguère et al., 2009; Haardt &
Madau, 2012; Faucher-Giguère, 2020). In the early universe, the UVB is
important for photoionizing the neutral hydrogen that fills the universe after
recombination, converting the universe from a neutral state back to a highly
ionized state. This transition is known as the epoch of reionization (EoR,
Shapiro & Giroux, 1987; Loeb & Barkana, 2001; Haiman, 2016). Importantly, the
UVB can also suppress the formation of small structures through the heating
associated with ionizing the intergalactic medium (IGM), making the gas
difficult to capture by low mass halos (Efstathiou, 1992; Cen & Ostriker,
1992; Babul & Rees, 1992; Haiman et al., 1996; Quinn et al., 1996). As a
result, star formation in these small halos is suppressed, as the gas is no
longer able to efficiently cool in them (Wang et al., 2009; Finlator et al.,
2011; Borrow et al., 2023).

Current cosmological simulations commonly implement the UVB as a spatially
uniform, time-dependent radiation field (e.g. Thoul & Weinberg, 1996; Quinn et
al., 1996; Katz et al., 1996; Weinberg et al., 1997; Vogelsberger et al.,
2014; Crain et al., 2015; Pillepich et al., 2018), with a time-dependent
spectrum as modeled in several works (Haardt & Madau, 1996; Faucher-Giguère et
al., 2009; Haardt & Madau, 2012; Faucher-Giguère, 2020). The photoionization
rate and the radiative heating and cooling rates associated with the UVB are
usually calculated with spectral synthesis codes like Cloudy (Ferland et al.,
1998), assuming that the gas is in collisional and photionized ionization
equilibrium. The inclusion of an externally imposed homogeneous UVB in
cosmological simulations makes them able to reproduce the average thermal
history of the IGM, and helps to match the observed galaxy population by
suppressing galaxy formation in very small dark matter halos (Cen & Ostriker,
1993; Benson et al., 2002; Okamoto et al., 2008).

On the other hand, in cosmological simulations of Milky Way (MW)-like galaxies
the UVB has been found to not significantly affect their evolution (Thoul &
Weinberg, 1996; Quinn et al., 1996; Weinberg et al., 1997). Since the virial
temperature of MW-like galaxies is about $10^{6}~{}{\rm K}$, the effective UVB
temperature is too low to significantly affect the gas in the halos of MW-like
galaxies. However, at the center of such galaxies, the gas temperature is
relatively low, suggesting that the gas may here be more prone to be affected.
On the other hand, the gas density is so high there that the intensity of the
UVB is likely not sufficient to significantly affect the high-density gas.

But the UVB is in reality not a spatially uniform radiation field, rather it
arises from the collective light emitted by a discrete population of galaxies.
It should therefore have a higher intensity when approaching the center of a
DM halo, where the photon sources – the stars and SMBHs – are located. Some
cosmological simulation projects try to account for this effect by considering
the radiation emitted from each discrete SMBH via a $r^{-2}$-law (Vogelsberger
et al., 2014; Pillepich et al., 2018; Davé et al., 2019). These studies found
that including such a local radiation contribution from SMBHs has non-
negligible effects on the galaxy evolution. For example, Davé et al. (2019)
reported that massive galaxies at $z=0$ become redder when accounting for the
local radiation from SMBHs, which implies that this local radiation field can
help quench massive galaxies.

There are also some theoretical works and simulation studies that considered
the radiation feedback from stars. For example, Hopkins et al. (2012)
performed several galaxy formation simulations with different halo masses,
considering both the short-range and long-range radiation pressures as well as
the photoionization heating from HII region by using cosmological zoom-in
simulation. They found that radiation pressure plays a key role in generating
galactic winds in massive galaxies. Roškar et al. (2014) simulated MW-sized
galaxies with short-range radiation pressure on dust, and found that such
radiative feedback can suppress star formation and the clumpiness of the cold
gas in the disk. However, these works did not take the photoionization heating
caused by the certain stellar populations into account.

Cantalupo (2010) investigated the impact of X-ray radiation emitted by hot
winds from young stars, supernova remnants, and X-ray binaries on the gas
cooling time. This differed from previous works which had only considered the
radiation from young stars themselves. They found that this additionally
included X-ray radiation significantly affects the gas cooling time, and that
the cooling time increases with increasing star formation rate (SFR). This
“stars quench stars” paradigm constitutes a form of preventive feedback, as
the X-ray radiation prevents the gas from cooling and thus ceases star
formation.

Inspired by Cantalupo (2010), Kannan et al. (2014) studied the impacts of
local stellar radiation on the evolution of MW-like galaxies, which is the
first work to consider local stellar radiation in cosmological simulations.
They separately considered the radiation of young and old stars, as well as
the X-ray radiation from the hot gas components considered by Cantalupo
(2010), then calculated the gas cooling and heating rate using Cloudy. They
incorporated the local stellar radiation in the McMaster Unbiased Galaxy
Simulations suite (MUGS, Stinson et al., 2010), which is a cosmological zoom-
in simulation project performed with the GASOLINE2 (Wadsley et al., 2004)
code. The radiation field was calculated by attaching the radiation
information to the gravitational tree structure in an optically thin
approximation. Kannan et al. (2014) found that local stellar radiation can
significantly affect the evolution of MW-like galaxies, suppressing star
formation by about 30% and lowering the central peak of the rotation curve.
Based on these findings, they introduced the term local photoionization
feedback (LPF) for the locally varying stellar radiation field.

Obreja et al. (2019) further improved the LPF implementation of Kannan et al.
(2014) and applied it in the project Numerical Investigation of a Hundred
Astrophysical Objects suite (NIHAO, Wang et al., 2015), which is another
cosmological simulation project performed by GASOLINE2. Obreja et al. (2019)
simulated a sub-$L^{*}$, an $L^{*}$, and a super-$L^{*}$ galaxies with and
without the LPF, and confirmed the basic conclusions obtained by Kannan et al.
(2014) for MW-like galaxies. Further, they found that the effects of the LPF
become weaker in more massive galaxies, and that the LPF also affects the CGM
properties. In independent work, Hopkins et al. (2018, 2020) implemented two
different numerical schemes as part of the FIRE project to study the impact of
stellar radiation feedback on galaxy evolution, both for dwarf- and MW-sized
galaxies. Differently from the conclusions reached in Obreja et al. (2019),
they found that the local stellar radiation has no impact on the evolution of
MW-like galaxies.

Besides galaxy evolution, the local radiation field can potentially also
affect the CGM properties (Stern et al., 2016; Suresh et al., 2017;
Oppenheimer et al., 2018; Qu & Bregman, 2018; Nelson et al., 2018; Fielding et
al., 2020; Strawn et al., 2021). However, most current CGM simulation studies
consider a uniform UVB when studying the ionization state of the CGM, and
there is so far comparatively little work including the local radiation.
Notably, Suresh et al. (2017) compare the OVI column density with and without
the galaxies’ central stellar radiation and found that it helps to increase
the OVI column density in the inner region. Oppenheimer et al. (2018) consider
the flickering AGN radiation and non-equilibrium ionization of the CGM gas in
EAGLE zoom simulations, finding that the AGN radiation can enhance the OVI
column density allowing the COS-halos observations (Tumlinson et al., 2011;
Werk et al., 2014) to be matched. Finally, Nelson et al. (2018) studied the
OVI column density in the TNG-100 simulation with several variations of the
feedback prescriptions. They found that the OVI column density becomes lower
when turning off the BH radiative feedback compared to a simulation with the
fiducial TNG model, for halo masses in the range $10^{11.5}<M_{\rm halo}/{\rm
M_{\odot}}<10^{12.5}$.

As mentioned above, there is no consensus in the field about how strongly
local stellar radiation can affect the evolution of MW-like galaxies, and the
impact of the local stellar radiation on the CGM is also not well understood.
To further investigate these two questions, we here revisit simulations of MW-
like galaxies that include local stellar radiation. We use the successful
model of the Auriga project (Grand et al., 2017) as the basis for our
simulations, and implement the local stellar radiation with a scheme
technically similar to the LPF treatment introduced by Kannan et al. (2014).

This paper is organized as follows. In Section 2, we briefly introduce the
Auriga project, describe the spectral energy distribution (SED) of the local
stellar sources we used, and the numerical implementation of the local
photoionization feedback (LPF) in the code. In Section 3, we compare the
simulated Auriga galaxies with and without LPF to study the impacts of LPF on
the evolution, the gas as well as stellar disk properties, and the rotation
curve of the simulated galaxies. In Section 4, we investigate the impacts of
LPF on the CGM gas properties, including the radial profile, the gas phase
distribution, and the properties of some selected species. In Section 5, we
compare our results with previous studies. Finally, we summarize our findings
and conclude in Section 6.

##  2 Methodology

###  2.1 Simulations

This work is an extension of the Auriga project (Grand et al., 2017), which
consists of a large suite of high-resolution magnetohydrodynamical simulations
of the formation of MW-like galaxies, and whose data was recently been made
publicly available (Grand et al., 2024). The computations utilize the ‘zoom-
in’ technique to follow galaxies within a full cosmological context. The
project has been carried out using the AREPO code (Springel, 2010), which is a
gravitational N-body and moving-mesh magnetohydrodynamical code. Due to its
quasi-Lagrangian nature, the moving-mesh approach can provide an automatically
adaptive resolution and manifest Galilean invariance at the discretized level
of the fluid equations. These properties are particularly useful for the often
highly supersonic flows occurring in galaxy formation.

The galaxy formation physics model we employ in this work is based on that of
IllustrisTNG (Weinberger et al., 2018; Pillepich et al., 2018). This in turn
employs a star formation and stellar feedback model which is based on the two-
phase ISM and non-local stellar wind subgrid prescription described in
Springel & Hernquist (2003). The gas in the simulation is eligible for star-
formation gas when its density is higher than $\sim 0.1\,{\rm cm^{3}}$ and the
temperature is sufficiently low. This specific density threshold value is
calculated in the SN feedback-regulated two-phase ISM model as a function of
the overall gas consumption timescale, with the latter matched to
observations. When a gas cell becomes star-forming, its thermal state is
assumed to follow an effective equation of state, and the star formation rate
(SFR) is taken to be proportional to the dynamical time of the gas. As a
central feedback channel, stellar winds from the star-forming region are
invoked. Since the multi-phase structure of the ISM is not spatially resolved,
the numerical implementation of this wind feedback is realized in terms of
temporary wind particles that are decoupled from the gas until they leave the
star-forming region, thereby mediating a non-local kinetic feedback just
outside the star forming region. For a more detailed description of the star
formation, stellar wind, and stellar yield subgrid model and its parameters,
we refer to Pillepich et al. (2018).

Black holes (BH) are modelled as sink particles in the centers of galaxies.
They can accrete gas with a rate estimated with the Bondi-Hoyle formula using
the local gas properties at the BH’s position. We employ the two-mode BH
feedback model introduced in the TNG galaxy formation model (Weinberger et
al., 2018), which includes both a thermal and a kinetic mode. When the black
hole accretion rate (BHAR) is higher than a critical value (i.e. for high
Eddington ratio), the BH feedback is assumed to be in a purely thermal form
(‘quasar mode’), whereas for low accretion rate it occurs in a kinetic form
(sometimes called ‘radio mode’). The critical BHAR where the transition occurs
is scaled with the BH mass itself, making it easier for massive BHs to
transition into the kinetic mode. This two-mode AGN feedback model can quite
successfully produce a galaxy stellar mass function and a bimodality in the
galaxy color distribution at $z=0$ consistent with observations.

The Auriga project selected 30 halos with mass $M_{\rm halo}=1-2\times
10^{12}\,{\rm M_{\odot}}$ from the dark matter only simulation of the EAGLE
project (Crain et al., 2015), and resimulated them with the the zoom-in
technique keeping the original box length of 100 cMpc. The adapted
cosmological parameters for the simulations were taken from the Planck
Collaboration et al. (2014) analysis.

In this work, we select a subset of five typical halos from the Auriga galaxy
set to study the impacts of the local stellar radiation on the galaxy and CGM
properties. The selected halos correspond to Au3, Au5, Au6, Au15, Au16 from
Auriga, as presented in Grand et al. (2017). However, note that we do not
employ the original Auriga galaxy formation model, but rather we simulate the
halos with a fiducial, slightly adjusted TNG galaxy formation model in which
we omit AGN radiative feedback. In addition, we resimulate them again by
accounting for a local stellar radiation field as described below. For
reference, the primary galaxy properties of these five halos at $z=0$, with
and without the local stellar radiation, are given in Table 1.

Table 1: The properties of the simulated galaxies at $z=0$. The columns are (1) model name; (2) halo mass defined by the FoF algorithm; (3) halo virial radius; (4) total stellar mass; (5) BH mass; (6) half-mass radius of young stars; (7) half-mass radius of total stars; (8) gas fraction within the halo. Run | $\frac{M_{\rm halo}}{10^{10}~{}{\rm M_{\odot}}}$ | $\frac{R_{200}}{\rm kpc}$ | $\frac{M_{\star}}{\rm 10^{10}~{}M_{\odot}}$ | $\frac{M_{\rm BH}}{\rm 10^{7}~{}M_{\odot}}$ | $\frac{R_{\rm young~{}star,~{}1/2}}{\rm kpc}$ | $\frac{R_{\rm total~{}star,~{}1/2}}{\rm kpc}$ | $\frac{f_{\rm gas}}{\%}$
---|---|---|---|---|---|---|---
Au3NoRad | 135.77 | 185.69 | 5.96 | 4.78 | 10.06 | 2.48 | 4.40
Au3StarRad | 132.54 | 184.25 | 4.49 | 9.13 | 14.71 | 1.60 | 3.53
Au6NoRad | 98.79 | 167.06 | 3.14 | 3.41 | 13.43 | 4.53 | 6.78
Au6StarRad | 97.50 | 166.33 | 2.48 | 2.78 | 7.75 | 4.01 | 7.10
Au16NoRad | 133.70 | 184.79 | 5.16 | 2.71 | 14.61 | 5.98 | 7.21
Au16StarRad | 136.29 | 185.97 | 4.19 | 4.41 | 18.50 | 6.15 | 8.51
Au5NoRad | 110.63 | 173.48 | 4.10 | 3.57 | 2.65 | 1.50 | 4.16
Au5StarRad | 111.54 | 173.96 | 3.17 | 4.68 | 2.84 | 1.37 | 5.67
Au15NoRad | 93.18 | 163.84 | 2.21 | 2.07 | 6.86 | 1.60 | 7.20
Au15StarRad | 91.97 | 163.12 | 1.87 | 2.73 | 4.12 | 1.79 | 6.48

###  2.2 Radiation and photoionization from stellar populations

In this section, we specify the photoionization feedback model that we
introduced as a physics extension in our simulations.

In general, the photoionizing radiation in the universe comes from AGN and
stars. For a single galaxy, the radiation sources can be further divided into
local radiation sources and radiation from the other galaxies. The cumulative
radiation of AGN and stars from other galaxies results in the well-known
cosmic UV background (UVB). This UVB radiation has been extensively studied
(Faucher-Giguère et al., 2009; Haardt & Madau, 2012; Faucher-Giguère, 2020)
and has been incorporated into many cosmological simulations. Here, we follow
the previous work but use the updated UVB spectrum produced by Faucher-Giguère
(2020).

For representing local radiation sources, one should aim to include both AGN
and stellar radiation. A simple model for AGN radiation has been incorporated
in the Auriga simulations following Vogelsberger et al. (2014). However, our
technical implementation of stellar radiation, as described in the following
section, is not readily compatible with this AGN radiation module. More
importantly, since the aim of this paper is to study the impact of local
stellar radiation, we prefer to consider this first in isolation for clarity.
We thus turn off the AGN radiation in the following and focus on the stellar
radiation alone, deferring an investigation of the joint impact of both AGN
and stellar radiation to future work.

####  2.2.1 Stellar radiation

Our setup for the radiation treatment from stars follows the approach of
Kannan et al. (2014) but with several modifications. In the following, we
concisely summarize our approach and refer readers to Kannan et al. (2014) for
additional details.

The stellar radiation is divided into too parts – radiation from star forming
regions and from old stellar populations. For radiation from ongoing star
formation, Kannan et al. (2014) consider the young stellar population and the
X-ray emission from SNe remnants. Kannan et al. (2016) further account for the
X-ray emission from high-mass X-ray binaries (HMXB), whose companions are O or
B type stars. We combine these three kinds of radiation in our model in a
single channel. For the radiation from the young stellar component, following
again Kannan et al. (2014), we do not consider the time evolution of the SED.
Instead, the SED of the young stellar radiation is fixed to the SED of a 10
Myr old stellar population for simplicity. We use the Binary Population and
Spectral Synthesis v2.3 (BPASS, Eldridge et al., 2017) model to generate the
spectrum of a stellar population at 10 Myr for a given constant star formation
rate.

The young stellar radiation will be partly absorbed by hydrogen and dust when
it leaves the star-forming region. Here, we follow the approach of Kannan et
al. (2014), which is similar to Cantalupo (2010), and describe the frequency-
dependent escape fraction $f^{\nu}_{\rm esc}$ as

$f^{\nu}_{\rm esc}=f^{\rm LL}_{\rm esc}+(1-f^{\rm LL}_{\rm esc})\,{\rm
e}^{-\tau_{\nu}}.$ (1)

The escape fraction at the Lyman limit frequency $f^{\nu}_{\rm esc}$ is set to
5%. Here, only absorption of the neutron hydrogen is taken into account, thus
$\tau_{\nu}=\sigma_{\nu}N{\rm(H^{0}})$. $\sigma_{\nu}$ is the cross-section of
neutron hydrogen and $N{\rm(H^{0}})$ is the column density of neutron hydrogen
in the star-forming region. The value of $N{\rm(H^{0}})$ is set to
$10^{20}~{}{\rm cm^{2}}$, as in Kannan et al. (2014).

For the X-ray emission from SNe remnants, we assume that it follows a Raymond-
Smith law (Raymond et al., 1975). Cerviño et al. (2002) reported that the
total energy released in the EINSTEIN band ($0.1-3.4\,{\rm keV}$) of an
individual SNe is $\sim 10^{36}\,{\rm erg~{}s^{-1}}$. Based on this value, we
calculate the SNe event rate from a young stellar population with a Chabrier
(2003) IMF and constant $1~{}{\rm M_{\odot}~{}yr^{-1}}$ SFR and normalize the
X-ray spectrum accordingly.

The HMXB emission in this work is assumed to follow a $\Gamma=2$ power law in
the $0.5-2.0\,{\rm keV}$ band. The scaling relation of SFR and total X-ray
luminosity comes from Anderson et al. (2013):

$L_{\rm HMXB}\approx 1.4\times 10^{39}\frac{\rm SFR}{\rm
M_{\odot}~{}yr^{-1}}~{}{\rm erg~{}s^{-1}}.$ (2)

For the radiation from the evolved stellar population, both radiation from old
stars and low-mass X-ray binaries (LMXB) is included. Similar to the radiation
from young stellar populations, we do not consider the time evolution of the
radiation of old stellar populations, for simplicity. Instead we fix the SED
of an old stellar population at an age of 2 Gyr under the simple stellar
population (SSP) assumption, i.e., the stars are assumed to have formed at the
same time. We also use the BPASS v2.3 model and assume that the escape
fraction of the old stellar radiation is unity.

The energy of the LMXBs comes from black holes or neutron stars accreting
material from their low-mass companion stars. We assume that the shape of the
SED from LMXBs is a $\Gamma=1.5$ power law in the $0.5-2.0\,{\rm keV}$ band,
and that the LMXB luminosity as a function of total stellar mass scales as

$L_{\rm LMXB}=10^{40}\frac{M_{\star}}{10^{11}{\rm M_{\odot}}}~{}{\rm
erg~{}s^{-1}}.$ (3)

####  2.2.2 Gas heating and cooling with irradiation

Having specified the SEDs of the radiation components as described above, we
can calculate the gas heating and cooling rates in the presence of the
radiation field. For the heating and cooling of the primordial gas, Compton
cooling due to the interaction between the cosmic microwave background (CMB)
and the gas, and other cooling/heating processes are considered separately in
the simulation. The treatment of the heating and cooling of the primordial
gas, and the gas cooling due to CMB, follows the calculation in the
IllustrisTNG galaxy formation model with small modifications. The primordial
gas is assumed to be always in an ionization equilibrium state since the
simulation time span (Hubble time) is much larger than the relaxation time
($\sim 10^{7}\,{\rm yr}$). This assumption is valid for most of the time
except at the epoch of reionization.

The number density of different species of the primordial gas can be
calculated by the following equilibrium rate equations (Katz et al., 1996):

$\alpha_{\rm H~{}II}~{}n_{\rm H~{}II}~{}n_{e}-\Gamma_{e{\rm
H~{}I}}~{}n_{e}~{}n_{\rm H~{}I}-\Gamma_{\gamma{\rm H~{}I}}~{}n_{\rm H~{}I}=0,$
(4) $(\alpha_{\rm He~{}II}+\alpha_{\rm d})~{}n_{\rm
He~{}II}~{}n_{e}-\Gamma_{e{\rm He~{}I}}~{}n_{e}~{}n_{\rm
He~{}I}-\Gamma_{\gamma{\rm He~{}I}}~{}n_{\rm He~{}I}=0,$ (5)
$\displaystyle\alpha_{\rm He~{}III}~{}n_{\rm He~{}III}~{}n_{e}+\Gamma_{e{\rm
He~{}I}}~{}n_{e}~{}n_{\rm He~{}I}+\Gamma_{\gamma{\rm He~{}I}}~{}n_{\rm
He~{}I}$ (6) $\displaystyle-(\alpha_{\rm He~{}II}+\alpha_{\rm d})~{}n_{\rm
He~{}II}~{}n_{e}-\Gamma_{e{\rm He~{}II}}~{}n_{e}~{}n_{\rm
He~{}II}-\Gamma_{\gamma{\rm He~{}II}}~{}n_{\rm He~{}II}=0,$ $\alpha_{\rm
He~{}III}~{}n_{\rm He~{}III}~{}n_{e}-\Gamma_{e{\rm He~{}II}}~{}n_{e}~{}n_{\rm
He~{}II}-\Gamma_{\gamma{\rm He~{}II}}~{}n_{\rm He~{}II}=0,$ (7)

where $\alpha_{i}$ is the radiative recombination coefficient, $\alpha_{\rm
d}$ is the dielectric recombination coefficient, $\Gamma_{ei}$ is the
collisional ionization rate for different species, and $\Gamma_{\gamma i}$ is
the photonionization rate for different species. The values of the radiative
recombination coefficient, the dielectric recombination coefficient, and the
collisional ionization rate are given in Katz et al. (1996). With the UVB and
stellar radiation field, $\Gamma_{\gamma i}$ can be calculated by

$\Gamma_{\gamma i}=\int^{\infty}_{\nu_{Ti}}\frac{4\pi}{h\nu}\sigma_{\nu
i}(J_{\nu,~{}{\rm UVB}}(z)+J_{\nu,~{}\star}){\rm d}\nu,$ (8)

where $J_{\nu,~{}{\rm UVB}}(z)$ is the redshift-dependent, frequency-
dependent, and solid-angle-averaged intensity of UVB, and $J_{\nu,~{}\star}$
is the frequency-dependent, solid-angle-averaged intensity of the total
radiation from the stellar population. As described in Section 2.2, the UVB
radiation is described by the tabulated results in Faucher-Giguère (2020), and
the radiation from the stellar populations is obtained by summing up the total
radiation of the stellar populations in the galaxies.

Based on equations (4)-(8) and the number conservation equations, the
abundance of different species can be solved. Gas cooling and heating can then
be calculated. The gas cooling due to collisional excitation, collisional
ionization, recombination, dielectric recombination and free-free collision
can be calculated using Table 1 in Katz et al. (1996). The photoionization
heating is calculated by

$H=n_{\rm H~{}I}\varepsilon_{\rm H~{}I}+n_{\rm He~{}I}\varepsilon_{\rm
He~{}I}+n_{\rm He~{}II}\varepsilon_{\rm He~{}II}$ (9)

where

$\varepsilon_{i}=\int_{v_{T_{i}}}^{\infty}\frac{4\pi}{h\nu}\sigma_{\nu
i}(h\nu-h\nu_{T_{i}})(J_{\nu,~{}{\rm UVB}}(z)+J_{\nu,~{}\star}){\rm
d}\nu~{}{\rm erg~{}s^{-1}}$ (10)

$\varepsilon_{i}$ due to UVB is based on the tabular results from Faucher-
Giguère (2020), while $\varepsilon_{i}$ due to the stellar radiation is
obtained by using CLOUDY.

Figure 1: Upper panels: Cooling curves of the gas for different ionization
parameters. The left panel shows the cooling curve when the ionization
parameter $U_{\rm YS}$ of young stellar radiation (as labeled, in units of
${\rm erg\,s^{-1}\,cm^{3}}$) is varied, while the right panel shows the
cooling curve for different ionization parameters $U_{\rm OS}$ due to old
stellar radiation. Lower panels: The time evolution of the temperature for a
gas cell with different initial temperature and ionization parameters, and a
prescribed gas density corresponding to a number density of $1\,{\rm
cm^{-3}}$. The left panel shows the time evolution of the temperature for
different $U_{\rm YS}$, while the right panel gives the temperature decline
for different $U_{\rm OS}$. The lines with different colors represent the time
evolution of the temperature of the gas cells with ionization parameters equal
to 0, $10^{-3}$, $10^{-1}$, $10^{1}$, and $10^{3}$ ${\rm
erg~{}s^{-1}~{}cm^{3}}$. The different lines with the same color represent the
time evolution of the temperature of gas cells with initial temperatures equal
to $10^{3}$, $10^{4}$, $10^{5}$, $10^{6}$, and $10^{7}\,{\rm K}$.

Other cooling/heating processes, including metal cooling/heating and Compton
cooling due to the stellar radiation field, are calculated in the form of an
input table for the simulations with CLOUDY under the ionization equilibrium
assumption. Combined with all the processes described above, the total
cooling/heating rate $\Lambda$ in the simulation can then be obtained as

$\Lambda=\Lambda_{\rm p}+\Lambda_{\rm C}+\Lambda_{\rm o},$ (11)

where $\Lambda_{\rm p}$ is the cooling/heating rate of the primordial gas,
$\Lambda_{\rm C}$ is the Compton cooling due to CMB, and $\Lambda_{\rm o}$ is
the cooling/heating rate due to other cooling/heating processes. The resulting
cooling/heating curves for different ionization parameters of young and old
stellar radiation $z=0$ are shown in the first row of Figure 1.

Since the cooling/heating rate due to other processes, $\Lambda_{\rm
o}=\Lambda(n,T,J_{\rm ys},J_{\rm os},Z,z)$, is a function of gas density $n$,
temperature $T$, radiation field from young stellar population $J_{\rm ys}$,
radiation field from old stellar population $J_{\rm os}$, metallicity $Z$ and
redshift $z$, we generate a 6-dimensional cooling/heating table and store it
in a HDF5 file. The redshift dependence of the function is due to the
variation of the UVB spectrum. Because metal cooling and heating due to the
stellar radiation only occur after stars have formed, we ignore the metal
cooling/heating from the UVB at $z>10$.

For the table of the cooling/heating function, the gas number density ranges
from $10^{-8}-10^{2}~{}{\rm cm^{-3}}$ in spacing of 1 dex in log space. The
temperature ranges from $10^{1}-10^{9}~{}{\rm K}$ with 0.5 dex steps in log
space. The mean intensity of the radiation from young and old stellar
radiation ranges from $10^{-7}$ to $10^{7}~{}{\rm erg~{}s^{-1}~{}cm^{-2}}$
with 1 dex in log space. The metallicity ranges from $10^{-3}-10^{1}~{}{\rm
Z_{\odot}}$ with steps of 1 dex in log space. The redshift covers $0-10$,
again with 1 dex in log space. These bins are relatively coarse due to the
high memory cost of more densely sampled parameter space. To investigate the
validity of our produce table, we have randomly chosen $10^{6}$ particles with
different states $(n,T,J_{\rm ys},J_{\rm os},Z,z)$ and evolved them with
different resolutions of the cooling/heating table and calculated the relative
difference of the results. Reassuringly, we found relative difference less
than 10% in $99\%$.

###  2.3 Numerical implementation of the radiation field

As described in the previous subsection, the gas cooling/heating rate is a
function of the local radiation field. Thus, to implement the photoionization
feedback into the simulation, we need to obtain the radiation intensity where
the gas grid is located.

For the ambient UVB radiation, we follow the treatment in previous work, i.e.
the UVB radiation is approximated as a redshift-dependent homogeneous
radiation field. The UVB intensity value we adopt follows the results of
Faucher-Giguère (2020), as mentioned above.

In contrast, we calculate the radiation from the local stellar population
similar to the method described in Kannan et al. (2014). First of all, to get
the strength of the radiation emitted from a young stellar population, the
luminosity is scaled with the total mass in young stars. We identify star
particles borne more recently than 100 Myr as young stars. The total
luminosity and spectrum of young stars are obtained using the method described
in Section 2.2.1. We assume that star particles whose age is older than 100
Myr belong to the old stellar population. Again, the luminosity of the old
stellar population is scaled in proportion with the stellar mass and its
luminosity can be calculated easily. We note that our criterion for the young
stellar population is less strict than (Kannan et al., 2014), who only
identified star particles borne more recently than 10 Myr as young stars.
However, as the typical age of stars on the main sequence with mass $\sim
8~{}{\rm M_{\odot}}$ is $\sim 30-50$ Myr, we adopt 100 Myr as an upper bound
in the current work to safely include them. We also have run various test
simulations with different criteria for the age cut, finding that the value of
the age criterion for young stars does affect the results at a quantitative
level, for example, smaller value of the age criterion for young stars lead to
weaker effects of the local star radiation. However, our conclusions are not
affected qualitatively by the value for the age cut. Hence, for the rest of
this study we adopte the 100 Myr value for our study of the effects of the
local stellar radiation.

Our approximate radiation transfer method for the local stellar emission
follows the methods described in Kannan et al. (2014) and Woods (2015). The
radiation intensity impinging on a gas cell from local radiation sources is
calculated using an octree by assuming that the gas outside the star-forming
region is optically thin. Therefore, the radiation intensity from a single
source at any location is calculated by

$j=\frac{L}{4\pi r^{\prime 2}},$ (12)

where $L$ is the luminosity of the single source and $r^{\prime}$ is the
distance between a given location and the radiation source. The luminosity of
a single source in turn can be obtained by the methods described in the
previous subsections. The total radiation intensity at any location can then
be simply calculated by adding the radiation from all local sources, with the
gravitational tree algorithm being repurposed for this task, allowing the
corresponding summation to be computed efficiently.

Figure 2: Upper panels: The time evolution of the total baryonic mass, the BH
mass, the stellar mass, and the gas mass from $z=5.5$ to $z=0$, as labelled,
in our five selected galaxies. The solid lines give the results from the NoRad
simulations, while the dashed lines represent the results from the StarRad
simulations. Bottom panels: The relative difference between the simulations
with and without stellar radiation. The NoRad simulations are treated as here
as fiducial reference simulations. Figure 3: Left panel: The stellar mass
growth for different subsequent time bins in five selected Auriga galaxies at
$z=0$. The time interval for each bin is 1 Gyr, and this is used here to
obtain a time-averaged version of the star formation rate history in order to
smooth out temporal fluctuations in the SFR. Right panel: The BH mass growth
for different time bins in five selected Auriga galaxies. Again, we use time
intervals equal to 1 Gyr in order to smooth the BHAR. Solid lines in both
subfigures represent the results from the NoRad simulations, while the dashed
lines give the results from the StarRad simulations. The relative difference
of the NoRad and StarRad simulations is shown in the bottom panels of each
panel. Figure 4: The time evolution of the gas mass selected according to
different criteria. From the upper to the bottom panels we show: all gas
within $0.1\,R_{\rm vir}$, cold gas ($T<10^{5.5}\,{\rm K}$) in the same
region, star-forming gas in the same region, and finally the gas mass within a
1 kpc aperture around the central BH. The solid lines represent the median
values of the results of our five simulated galaxies, while the shaded regions
give the corresponding maximum and minimum values over this set of galaxies.
Both simulations with and without local radiation are shown, as labelled.

##  3 The impacts of local ionization feedback on galaxy properties

In this section, we analyze and discuss the differences in the time evolution,
in the properties of gas and stellar disks, and in the rotation curve between
simulated galaxies with and without local stellar radiation.

Before we consider our results for the simulated galaxies, we investigate the
time evolution of the temperature for several fiducial gas states with
different initial temperatures and different ionization parameters of young
stellar radiation $U_{\rm YS}\equiv J_{\rm YS}/n$ and old stellar radiation
${U_{\rm OS}}\equiv J_{\rm OS}/n$. By analyzing the temperature evolution of
these illustrative test situations, we can develop a first rough understanding
of the heating effects of the radiation on the expected evolution of gas
thermal states.

The second row of Figure 1 shows the time evolution of the temperature of test
gas cells with initial temperatures $T=10^{3}$, $10^{4}$, $10^{5}$, $10^{6}$,
and $10^{7}\,{\rm K}$, and $U_{\rm YS}$ or $U_{\rm OS}$ equal to $0$,
$10^{-3}$, $10^{-1}$, $10^{1}$, $10^{3}\,{\rm erg\,s^{-1}\,cm^{-2}}$,
respectively. At first glance, we can infer from the figure that the
temperature will eventually reach an equilibrium temperature when radiation is
included, which corresponds to the point where the cooling rate equals the
heating rate. This equilibrium temperature is sensitive to the strength of the
young stellar radiation. It can increase from $\sim 10^{3.5}\,{\rm K}$ at
$U_{\rm YS}=10^{-3}\,{\rm erg\,s^{-1}\,cm^{-2}}$ to $\sim 10^{5}\,{\rm K}$ at
$U_{\rm YS}=10^{-3}\,{\rm erg\,s^{-1}\,cm^{-2}}$. Note that the equilibrium
temperature due to the old stellar radiation is always close to or even lower
than $10^{4}\,{\rm K}$, and in general, it is lower than the equilibrium
temperature due to the young stellar radiation.

The stellar radiation can also delay the cooling when the gas temperature is
mildly higher than the equilibrium temperature, and this effect is more
pronounced for the old stellar radiation. For example, for both the young and
old stellar radiation, the gas temperature of the test gas configuration with
an initial temperature equal to $10^{5}~{}{\rm K}$ stays longer higher at any
given time for a higher ionization parameter, before the equilibrium
temperature is reached. This conclusion still holds in the presence of young
stellar radiation for gas with an initial temperature equal to $10^{6}$ or
$10^{7}\,{\rm K}$, but for high ionization parameter. For the old stellar
radiation, this effect becomes weaker, but it can still be found for $U_{\rm
OS}\ga 10^{3}\,{\rm erg~{}s^{-1}\,cm^{-2}}$ in gas with initial temperature
equal to $10^{6}$ or $10^{7}\,{\rm K}$.

###  3.1 Global time evolution of different baryonic components

Figure 5: The distribution of the radial velocity of gas within $0.1\,R_{\rm
vir}$ for times $z<0.1$. The data is averaged over five selected Auriga
galaxies. The black line represents the simulation results from the NoRad
runs, while the red line is for the corresponding StarRad runs.

Figure 2 shows the time evolution of the total baryonic mass, the BH mass, the
stellar mass, and the gas mass in five selected halos from $z=5.5$ to $z=0$.
The first panel gives the evolution of the baryonic mass and the relative
difference of five simulated halos with the NoRad and StarRad galaxies. In all
but the Au3 galaxies, the relative differences between the NoRad and StarRad
galaxies are $<20\%$ and do not have a clear trend, implying that the local
photoionization source has little impact on the gas accretion of the halos.
The only exception is Au3, where the total baryonic mass decreases in the
StarRad galaxies compared to the NoRad galaxies after 6 Gyr. The third panel
of Figure 2 also shows the signal of the large variation in the time evolution
of the stellar mass of Au3. The stellar mass of Au3 doubles around 6 Gyr,
which implies that this system encounters a major merger. Combined with the
third panel of Figure 2, the large difference in the Au3 galaxy may be caused
by the AGN feedback during the major merger at 6 Gyr. Since AGN feedback will
cause strong disturbances of the gas and nonlinearity of the fluid behavior,
the subsequent evolution may experience strong differences.

The second panel of Figure 2 shows the growth of the BH mass and the relative
differences between the NoRad and StarRad galaxies in all simulated halos. In
the figure, most of the simulated galaxies have a positive value for the
relative differences, which indicates that the BH is larger in the StarRad
galaxies than in the NoRad galaxies. The difference grows quickly and can be
greater than $50\%$ in the first 2 Gyr, but then becomes lower at later times.
This is also reflected in the right panel of Figure 3, which shows the growth
of the BH mass over time intervals of 1 Gyr. In the first 6 Gyr, the relative
difference of the BH mass is in general larger in the StarRad model than in
the NoRad model. The BH mass growth can even become six times higher in some
simulated galaxies. After 6 Gyr, the difference becomes small, however, and
the relative difference can even become negative, thus reversing the trend.

These results indicate that local stellar radiation can strongly suppress the
growth of BH at high redshift, while having only a minor impact at later time.
This may be caused by the suppression of star formation in dense gas when
local stellar radiation is considered. Since our model prevents gas cells from
becoming star-forming when the gas ionization parameter is larger than 1,
these gas cells can stay at relatively low temperature because their
equilibrium temperature is lower than the effective equation of state
temperature calculated by the star formation model in TNG.

We similarly find that the relative differences in the growth of the stellar
mass between the NoRad and StarRad galaxies are higher at earlier times. In
general, the relative differences increase with time in the first 6 Gyr and
then reduce to $\sim 25\%$ after 6 Gyr. This result implies that the local
stellar radiation can suppress star formation preferentially at high redshift,
which is a robust conclusion for the Auriga galaxies with our photoionization
feedback described above. At late times, the stellar mass can be suppressed by
around $25\%$, which is still stronger than the results reported in Obreja et
al. (2019) and Hopkins et al. (2020). We will discuss the comparison with
these two works in Section 5.

Note, however, that these trends are not always very clear at $z\geq 2$. In
fact, models Au5 and Au15 have actually larger stellar mass with the StarRad
model than with the NoRad model at high redshift. Given the more robust trend
at low redshift, this may be because local stellar radiation is only effective
in suppressing star formation once a disk has formed. The results in Kannan et
al. (2014), where star formation was found to be suppressed by $\sim 30\%$ in
a disk galaxy, are consistent with this conjecture. In contrast, at high
redshift, the galaxies are still in an intense formation phase, where their
virial temperatures are lower, the gas infall rates are very high, and stable
disks have not yet formed. The impact of local stellar radiation in this
environment is apparently weaker and likely more complicated to understand. In
addition, we also expect a larger system to system variation at these early
times. A further discussion of the high-redshift evolution is beyond the scope
of this work, but we intend to investigate it in more detail in future work.

Figure 6: Stellar images of different Auriga galaxies, stacked from $z=0.1$ to
$z=0$, a timespan still short compared to the bar rotation period. From the
left to the right columns we show Au3, Au6, Au16, Au5 and Au15. The top row
shows the stellar distributions of the NoRad simulations, while the bottom row
gives the stellar images of the corresponding StarRad simulations. Figure 7:
Stacked radial stellar surface density profiles of different Auriga galaxies,
averaged from $z=0.1$ to $z=0$. We show results for simulations both with and
without local stellar radiation, as labelled. The shaded areas illustrate the
10 to 90 percentiles of the five galaxies in our simulation set, while the
solid lines are the medians.

We further investigated the stellar age distribution of simulated galaxies at
$z=0$. The left panel of Figure 3 shows the stellar age distribution of the
NoRad and StarRad galaxies at $z=0$ and their relative differences. The
differences indicate that the stellar mass formed at $z\sim 2$ to $z\sim 0.7$
is systematically lower in the StarRad galaxies than in the NoRad galaxies. At
$z\ga 2$, the value of the relative differences is slightly higher than zero
in Au5 and Au15, which is consistent with the results in the third panel of
Figure 2. At $z\la 0.7$, the relative differences exhibit large scatter, which
may be related to the growing importance of AGN feedback at late times. Since
the BH mass becomes large there, the cumulative energy released by the AGN
grows as well, and the impact of the associated energy feedback becomes
significant (Zinger et al., 2020). Note that although the stellar mass also
increases with time, the impact of the radiation from old stars on the gas is
much weaker than that of radiation from young stars. Given also the declining
star formation rate at low redshift, it is then not surprising that the impact
of the AGN feedback becomes more important than the local stellar radiation.
Thus, the impact caused by local stellar radiation becomes subdominant to the
impact of AGN feedback and cannot be distinguished separately any more.

The final panel of Figure 2 shows the time evolution of the total gas mass and
the relative differences between the NoRad and StarRad galaxies in all
simulated halos. Similar to the time evolution of the baryonic mass, the
figure shows good convergence in the first 8 Gyr but highlights larger
relative differences between the NoRad and StarRad galaxies at later time.
Since the stellar mass exhibits a systematic decrease with radiation included,
the gas mass shows a corresponding increase in its proportion relative to the
stellar mass and the the total baryon mass later than $\sim 6$ Gyr. This
result again confirms that less gas mass is converted to stars in the StarRad
galaxies.

The time evolution of the total baryon mass in both the NoRad and StarRad
galaxies indicates that the local stellar radiation will not affect the total
gas accretion into the halo. To further investigate whether gas accretion into
galaxies remains equally unaffected, we show the time evolution of the gas
mass within $0.1\,R_{\rm vir}$ in both the NoRad and StarRad galaxies in the
first panel of Figure 4. We directly put the five NoRad/StarRad galaxies
together for a statistical robust mean trend. Similarly to the total baryon
mass, the difference in the gas within $0.1\,R_{\rm vir}$ in the first 8 Gyr
is small, indicating that the feedback processes are still weak and do not
provide enough energy to push the gas out. After 8 Gyr, the median difference
between the two is still small. However, the range between the 10 and 90
percentile values in the StarRad galaxies becomes lower and smaller than in
the NoRad galaxies. This is because local stellar radiation can heat the gas,
either preventing it from flowing in, or even pushing it to larger radii
through thermal pressure.

The second panel of Figure 4 shows the time evolution of the cold gas. The
cold gas is here defined as gas with a temperature lower than $10^{5.5}\,{\rm
K}$. The result is almost the same as seen for the time evolution of the total
gas within $0.1\,R_{\rm vir}$, indicating that the dominant gas component
within $0.1\,R_{\rm vir}$ is cold gas. The reduction of the cold gas at low
redshift can be ascribed to the heating effect of local stellar radiation. The
third panel of Figure 4 shows the time evolution of the star-forming gas mass
within $0.1R_{\rm vir}$. Unlike the time evolution of the total gas mass shown
in the first panel of Figure 4, the total amount of star-forming gas in the
StarRad galaxies is systematically lower than in the NoRad galaxies from high
redshift. This follows the trend of the star formation history shown in the
first panel of Figure 3, again implying that local stellar radiation affects
star formation by reducing the amount of star-forming gas, but not the total
gas mass.

The final panel of Figure 4 gives the time evolution of the gas mass within 1
kpc around the central BH. From the figure, we can see that the median value
of the gas mass evolves similarly in both the NoRad and StarRad galaxies.
However, the variation of the gas mass in StarRad galaxies is much lager than
in NoRad galaxies, especially in the first 5 Gyr. A higher gas mass in
principle implies that the BHs have more gas to accrete from, so this result
can partly explain why the BH mass can grow faster in the StarRad galaxies at
high redshift.

Figure 5 provides another view on the impact of the stellar radiation by
showing the distribution of radial velocities of the gas, in units of the
virial velocity. From the figure we can infer that the amount of gas with a
negative inflow velocity as well as with a positive outflow velocity (wind
feedback) is reduced in the StarRad galaxies compared to the NoRad galaxies.
This implies that the local stellar radiation effectively reduces the amount
of material available for star formation, and thus indirectly the feedback
strength. Note that the local stellar radiation thus acts more as a preventive
process that keeps gas from flowing in in the first place.

###  3.2 Disk properties at $z=0$

The previous section examined the impact of local ionization feedback on the
formation history of the simulated galaxies. In this section, we analyze the
impact of local ionization feedback on both stellar and gas disk properties
$z=0$.

####  3.2.1 Stellar disk

Figure 6 shows stacked stellar images of five simulated galaxies with and
without local stellar radiation, where the surface density of each simulation
is averaged from $z=0.1$ to $z=0$. The upper row of panels gives the simulated
galaxies without local stellar radiation, while the bottom row of panels shows
the corresponding simulated galaxies with local stellar radiation. At first
glance, the NoRad galaxies have slightly higher surface brightness in the
face-on view. The NoRad galaxies display more blue light outside 5 kpc, and
more red light at the center. Except for the surface brightness, both the red-
light and blue-light of NoRad galaxies are also more extended. These results
imply that local stellar radiation can also affect the stellar disk properties
by influencing the spatial location of star formation. Similar conclusions can
also be drawn from the edge-on views.

To compare the stellar distribution quantitatively in both the NoRad and
StarRad galaxies, we show the one-dimensional surface density of total stars,
old stars, and young stars in Figure 7. Young stars and old stars are
separated by their age. When the stellar age is younger than 3 Gyr, the stars
are identified as young stars, otherwise as old stars. Again, we stack all of
the NoRad/StarRad galaxies from $z=0.1$ to $z=0$. The solid line represents
the median value of the surface density, while the shaded area in the figure
represents the 10 to 90 percentile range of the distribution of measurements.

The first panel of Figure 7 gives the surface density of the total stars of
the stacked NoRad and StarRad galaxies. We find that the total stellar surface
density in the NoRad galaxies is slightly higher than in the StarRad galaxies.
The difference of the median values between the StarRad and NoRad galaxies is
relatively small at the center, but becomes larger outside 10 kpc. For
determining the 10 to 90 percentiles, note that the coverage within 10 kpc is
generally lower for StarRad galaxies than for NoRad galaxies. Outside 10 kpc,
the coverage is roughly the same except that the lower limit in the StarRad
galaxies is slightly smaller.

To investigate the relative contributions of old and young stars to the
surface density, we also show the surface density of old and young stars in
the second and third panels of Figure 7, respectively. As for the total, the
solid lines in the second panel represent the median values of the surface
density of old stars, while the shaded area indicates the 10 and 90 percentile
values. For the median value, we find that the surface density of the NoRad
and StarRad galaxies is similar to each other within 50 kpc. The median and 90
percentile values of the surface density in the StarRad galaxies are slightly
lower than the NoRad galaxies within $\sim 15$ kpc, which is similar to the
surface density of total stars. This is consistent with the finding that the
NoRad galaxies are brighter in the red band at the center of the stellar
images. However, in general, the differences in the surface density of the old
stars between the NoRad and StarRad galaxies are rather small.

This is different for the surface density of young stars, shown in the third
panel of Figure 7. Here the median value for the StarRad galaxies is lower
than for the NoRad galaxies over the entire 50 kpc range. This relative
difference becomes larger when the radius increases. This result can explain
the higher surface brightness in the blue band in the NoRad galaxies outside
of 5 kpc in Figure 6. The large spread between the 10 and 90 percentile values
is here primarily induced by the large size variations of the five simulated
galaxies. Since the stellar mass distribution in the outer regions is disk-
like, these results imply that local stellar radiation is comparatively
effective in suppressing the star formation in disk-like configurations.

Figure 8: Time evolution of the half-mass radius of total stars (left panel)
and young stars only (right panel). The solid lines represent the simulation
results for our NoRad runs, while the dashed lines give the results for
StarRad simulations. The relative differences of the corresponding NoRad and
StarRad runs is shown in the bottom panels.

The above results confirm that local stellar radiation can efficiently
suppress the stellar surface density. We further investigate its impacts on
the concentration of the stellar distribution. To this end we have fitted the
stellar surface density with Sersic profiles, using separate fits for the
surface density of old and young stars. We note, however, that the fitting
parameters can sometimes be sensitive to details of the fitting procedure,
calling for caution in interpreting their precise values. This is also
reflected in relatively noisy time evolutions of the Sersic index and the
effective fit radius, which obscures clear trends. Instead of using a
parametric fit, we therefore in the following characterize the concentration
of the stellar mass through the stellar half-mass radius. This non-parametric
measure avoids a sensitivity to details of the fitting procedure and thus
tends to be more robust.

Figure 8 shows the time evolution of the half-mass radius of all stars and of
the young stars, and of the relative difference in our five NoRad and five
StarRad galaxies at $z<1$. The left panel of the figure gives the time
evolution of the half-mass radius of all stars. Each of the simulated galaxies
has a smaller half-mass radius of all stars in the StarRad variant compared
with the NoRad variant, except for Au15. Recall that in the first panel of
Figure 3 we found that the formation history of Au15 is different from those
of the other galaxies, making this galaxy a bit of an outlier. In particular,
the stellar age of Au15 is younger than that of the other galaxies, and its
star formation reaches its peak at $z<1$. So aside from Au15, the smaller
half-mass radius of all stars in the other StarRad galaxies implies that local
stellar radiation can increase the concentration of the galaxies. On other
hand, through its different formation history Au15 corroborates again the
result that the impact of local stellar radiation differs in the early- and
late-stages of disk galaxy formation.

The second panel of Figure 8 shows the time evolution of the half-mass radius
of young stars. Unlike the half-mass radius of the total stars, the half-mass
radius of young stars does not yield a consistent trend over the entire
period. At 6 to 8 Gyr, it is smaller in the StarRad than in the NoRad galaxies
– except for Au15. After 8 Gyr, it does however not show an obvious and clear
systematic trend any more for the full sample of simulated galaxies. Here the
half-mass radius of young stars in the StarRad galaxies can be found to be
both larger or smaller than in the NoRad galaxies. This may be due to the AGN
feedback becoming important after 8 Gyr and beginning to affect the gas,
especially in the inner regions.

####  3.2.2 Gas disk

Figure 9: Maps of the total gas surface density of different Auriga galaxies
at $z=0$. From the left to the right columns, we display Au3, Au6, Au16, Au5,
and Au15. The top row gives results for the NoRad simulations, both in face-on
and edge-on projections (large and small subpanels, respectively), while the
bottom row shows the corresponding StarRad simulations. Figure 10: Radial
surface density profiles of total gas (top panel), HI gas (middle panel) and
just star-forming gas (bottom panels) of our sample of five Auriga galaxies,
averaged for the cases with and without local stellar radiation (as labelled),
and over time from $z=0.1$ to $z=0$. The shaded area shows the 10 to 90
percentiles of the measured values for the stacked simulation set, while the
solid lines are the medians. Figure 11: Same as Figure 9, but for the HI gas
surface density.

Figure 9 shows the total gas surface density in different simulated galaxies
at $z=0$. The first and second panels display the face-on and edge-on views of
the total gas surface density in our five NoRad galaxies, while the third and
fourth panels give the corresponding information for our five StarRad
galaxies. The face-on views of the gas disk in the NoRad and StarRad galaxies
highlight that they have similar morphology and size. The disk heights also do
not indicate significant differences, based on comparing the edge-on views of
both NoRad and StarRad galaxies. This result is consistent with the finding in
the first panel of Figure 4, which shows that the median value of the total
gas within $0.1\,R_{\rm vir}$ is similar in the NoRad and StarRad galaxies.
Since the gas disks in both models have similar morphology and size, the
angular momentum of the inflowing gas should also be similar, again providing
indirect evidence that the local ionization feedback does only weakly affect
the inflow and is not an overly effective preventive feedback process.

To further investigate the differences of the gas disks in the NoRad and
StarRad galaxies, we compare the radial profile of the total gas surface
density in the first panel of Figure 10. The profiles are stacked and averaged
from $z=0.1$ to $z=0$. In general, the total gas surface density in the NoRad
galaxies is higher than in the StarRad galaxies, although the difference is
very small. This confirms the conclusion obtained by eye from the 2D
projections of the total gas surface density. The slightly lower median value
in the StarRad runs within $10\,{\rm kpc}$ occurs because local stellar
radiation can heat the cold gas. As the gas temperature increases, the gas
density needs to reduce to satisfy the pressure equilibrium.

Figure 11 shows the HI gas surface density of different simulated galaxies at
$z=0$. Similarly to Figure 9, the top row of panels gives the face-on and
edge-on views of five NoRad galaxies, and the bottom row of panels is for our
five StarRad galaxies. Unlike for the total gas surface density, the HI gas
surface density exhibits significant differences between the NoRad and StarRad
cases. The face-on view indicates that the HI disk of the NoRad galaxies
features a relatively smooth distribution. However, the HI disk in the StarRad
galaxies appears clumpy. In addition, the HI disk in the StarRad galaxies has
a smaller size. From the edge-on view, we see that the HI disk in the StarRad
galaxies has also a lower disk height. These results are consistent with the
findings in Obreja et al. (2019).

These differences are quantified in the stacked radial profiles of the HI gas
surface density shown in the second panel of Figure 10. From the figure we see
that the HI gas surface density in the StarRad galaxies is systematically
lower than in the NoRad galaxies. The median value at 5-25 kpc in the StarRad
galaxies is around $10^{2}\,{\rm M_{\odot}\,pc^{-2}}\sim 10^{21}\,{\rm
cm^{-2}}$, which is one order of magnitude lower than in the NoRad galaxies.
Another difference is that the HI disk surface density drops rapidly at $\sim
25-35\,{\rm kpc}$ in the StarRad galaxies. These differences in HI morphology
and surface density between NoRad and StarRad galaxies indicate that local
stellar radiation can effectively ionize the HI gas in the disk, making the HI
disk smaller, clumpier, and giving it a lower surface density.

The final bottom panel of Figure 10 shows the stacked surface density of star-
forming gas in the NoRad and StarRad galaxies. We find that the surface
density of star-forming gas is systematically lower in the StarRad galaxies
compared with the NoRad galaxies, reflecting the results for the HI gas disk.
This lower star-forming gas surface density of the StarRad galaxies can also
explain the lower disk surface density of young stars seen in Figure 7.

###  3.3 Rotation curve

Figure 12: Comparison of the $z=0$ rotation curve of our five simulated Auriga
galaxies. Solid lines give the results for the NoRad runs, while dashed lines
are for the corresponding StarRad simulations where local stellar radiation is
included. The bottom panel shows their difference relative to the NoRad
simulations.

Arguably one of the most notorious difficulties in simulating MW-like galaxies
is the overcooling problem (Balogh et al., 2001). The low angular-momentum gas
accreted into the DM halos at the early stages of galaxy formation needs to be
removed to avoid forming too many stars in a central bulge (Navarro & Benz,
1991; Governato et al., 2004; Scannapieco et al., 2012). One of the symptoms
of the overcooling problem, if not sufficiently prevented, is the existence of
a central peak in the rotation curve. Kannan et al. (2014) found in their
simulations that local stellar radiation can effectively suppress this peak,
indicating that local stellar radiation is a process that can contribute to
solving the overcooling problem. Furthermore, Obreja et al. (2019) implemented
local ionization feedback scheme in a comparable way to Kannan et al. (2014)
in the NIHAO simulations, finding a similar effect on the rotation curve as in
Kannan et al. (2014). We here therefore also check the effect of the local
stellar radiation on the rotation curve in our simulations.

Figure 12 shows the averaged $z=0$ rotation curve of our five NoRad and five
StarRad galaxies, respectively. Compared to the NoRad case, the peak of the
rotation curve shows some suppression. However, this suppression is relatively
small compared to the results in Kannan et al. (2014) and Obreja et al.
(2019). The relative difference in the peak of the rotation curve between the
NoRad and StarRad galaxies is less than $10\%$. The reason why the local
stellar radiation in our work shows an overall weaker effect may be related to
the star formation model. In Kannan et al. (2014) and Obreja et al. (2019),
the multiphase gas is in part directly resolved even at high density. Their
temperature and density thresholds for the star formation model are
$15000~{}{\rm K}$ and $10~{}{\rm cm^{-3}}$, and the gas can still be affected
by radiation even if it is star-forming. The density threshold for the star
formation model in our work is instead $\sim 0.1~{}{\rm cm^{-3}}$, and the
thermal state of the denser gas is set by an effective EOS, suppressing
resolved multi-phase structure in the star-forming gas. While gas with
temperature lower than $\sim 10^{4.5}~{}{\rm K}$ is more vulnerable to local
stellar radiation, once the gas becomes star-forming, its effective
temperature is no longer affected by local stellar radiation in our model,
except for our added criterion that gas with high ionization parameter will
not be allowed to form stars. Especially at the center of galaxies, the gas
density is relatively high, and we expect only mild effects from this
ionization parameter criterion. This may be the main reason why the
suppression effects of local stellar radiation on the peak of the rotation
curve are relatively weak in our simulations.

##  4 Structure of the circum-galactic medium

Figure 13: Time evolution of the physical radius in units of the virial radius
where the old and young stellar radiation intensity equals the ambient UVB
intensity, separately shown for our five simulated Auriga galaxies. The bottom
panel is for the old stellar component, whereas the bottom panel is for the
young stars. Note the different vertical scale in the two panels. Unlike the
radiation field from the old stars, the young stellar component can dominate
the UVB throughout the halo, and even beyond.

Previous work studying the impact of radiation on CGM gas usually considered
the influence of a UV radiation background, whereas only few studies
considered the radiation from the galaxies themselves (e.g., Suresh et al.,
2017; Oppenheimer et al., 2018; Nelson et al., 2018). However, the radiation
from nearby stars may play a significant role for the CGM region.

To illustrate this point, Figure 13 shows the time evolution of the radius
where the young/old stellar radiation is equal to the UVB intensity, in units
of the virial radius. For old stars as sources, this radius gradually
increases with time, due to the monotonic growth of the stellar mass and the
decrease of the UVB intensity at $z<3$. The radius varies between
$0.1-0.4\,R_{\rm vir}$ over time, which means that the stellar radiation of
old stars is comparable of higher than the UVB only in the innermost regions
of the CGM. Furthermore, since the older stellar populations emit few ionizing
photons, this significantly limits the overall impact this radiation component
can have on the CGM gas.

For younger stars, on the other hand, the radius where their influence
dominates gradually decreases with time. It evolves from $\sim 10$ times the
virial radius to $\sim 2$ times the virial radius from $z=3$ to $z=0$, which
reflects that the star formation history reaches a peak at $z\sim 3$ and then
decreases with time. This result implies that within the whole CGM gas the
radiation field is actually dominated by the young stellar radiation from the
central galaxies instead of the UVB. Unlike for the radiation from old stars,
the fraction of ionizing photons in the young stellar radiation can reach 23%,
indicating that this radiation component can potentially have a significant
impact on the CGM gas, even outside the halo itself. In this section, we will
thus analyze the influence of the local stellar radiation on the CGM
properties.

###  4.1 Radial profiles of gas properties

Figure 14: Radial profiles of number density (top panel), temperature (middle
panel) and metallicity (bottom panel) of our Auriga galaxies, stacked for the
different models simulated with and without radiation (as labelled) and
averaged from $z=0.1$ to $z=0$ to reduce the influence of temporal and
stochastic fluctuations. The shaded areas show the 10 to 90 percentile over
the simulation set, while the solid lines give median values.

We begin by studying the impact of local stellar radiation on the global
structure of the gaseous halos. Figure 14 gives the stacked radial profiles of
the gas number density, temperature, and metallicity, in both the NoRad and
StarRad galaxies. Only non-starforming gas is included, i.e. the ISM component
is excluded. The data are stacked for all galaxies with the same physical
model, and are averaged over $z=0.1$ to $z=0$ to eliminate temporal
fluctuations. The solid line represents the median value of the individual
measurements, while the shaded area indicates the 10% to 90% percentiles.

The number density profile, shown in the first panel of Figure 14, can
evidently be separated into two characteristic regions. Within $0.3\,R_{\rm
vir}$, the number density is still comparatively high, reflecting
contributions from the extended gas disk. Outside $\sim 0.3\,R_{\rm vir}$, the
number density then quickly drops to $10^{-4}-10^{-5}\,{\rm cm^{-3}}$. Outside
this radius of $\sim 0.3\,R_{\rm vir}$, the differences in the number density
profiles between the NoRad and StarRad galaxies appear negligible. Some of the
discrete peaks in the shaded region are due to cold gas debris or satellites,
and these features would be expected to disappear when averaging over a still
larger galaxy sample. In the radial range $0.1-0.3\,R_{\rm vir}$, the StarRad
galaxies have a slightly higher number density than the NoRad galaxies. Note
that this corresponds to radii of $\sim 20-60\,{\rm kpc}$ in our simulated
galaxies. For the total gas surface density shown in the first panel of Figure
10, the StarRad galaxies have a slightly lower value than the NoRad galaxies
in the corresponding region. This difference arises from the omission of the
star-forming gas in this profile. Since the local stellar radiation prevents
some gas from becoming eligible for star-formation, the star-forming gas
itself will have a lower density, and the median value of the gas number
density in the StarRad galaxies is expected to be lower than in the NoRad
galaxies.

The second panel of Figure 14 shows the radial profile of the temperature.
Unlike the number density, the median values of the temperature profiles show
good consistency for the entire CGM region. While this result suggests that
the local stellar radiation does not have a strong heating effect in the CGM
region, it should be noted that the temperature tends to closely track the
virial temperature of the halo even in the presence of strong heating, and it
is set ultimately by gravity. However, there are still some indications even
in this result that local stellar radiation influences some aspects of the CGM
gas. In particular, the lower limit of the shaded region is slightly higher
the StarRad galaxies than for the NoRad galaxies. This can be understood due
to the ability of the local stellar radiation to heat gas with temperature
lower than $\sim 10^{4.5}~{}{\rm K}$.

Finally, the third panel of Figure 14 shows the radial profile of the
metallicity. From the figure we can see that the metallicity in the StarRad
galaxies is lower than that in the NoRad galaxies over the entire CGM region.
The difference decreases with increasing radius and disappears at $\sim R_{\rm
vir}$. This radially dependent difference suggests that it originates in the
center of the galaxy. Recall that Figure 5 has shown that local stellar
radiation damps the gas outflow, which is due to a suppression of star
formation and thus a reduction of feedback energy injection. Because the gas
outflows become weaker, they carry fewer metals into the CGM region.

###  4.2 CGM phases by mass

Figure 15: The temperature–number density distributions of CGM gas
($0.1<r/R_{\rm vir}<1$) in the upper row, inner halo gas ($0.1<r/R_{\rm
vir}<0.3$) in the middle row, and outer halo gas ($0.3<r/R_{\rm vir}<1$) in
the bottom row. The left column represents the results from the NoRad
simulations, while the right column gives the results of the StarRad
simulations that include local stellar radiation. The data is stacked for all
five simulated Auriga halos, and over outputs from $z=0.1$ to $z=0$. Figure
16: Left panel: Distribution of the ionization parameter due to the young
stellar radiation as a function of radius in the StarRad simulations. Right
panel: The distribution of the ionization parameter of the CGM gas due to
young stars as a function of temperature. The color key indicates the amount
of gas mass in the corresponding region of the distribution.

In the following we discuss the amount of gas in different CGM phases. To this
end we divide the CGM into two regions: the inner halo ($0.1<r/R_{\rm
vir}<0.3$) and the outer halo ($0.3<r/R_{\rm vir}<1$). Figure 15 shows the gas
number density vs temperature distribution of the whole halo, and separately
for the inner and outer halo, both in the NoRad and StarRad galaxies. As done
earlier, we again combine the data of all galaxies with the same physics
model, and average over the time span $z=0.1$ to $z=1$ to obtain statistically
more robust results. The star-forming gas of the ISM is excluded.

In the top row of Figure 15, we show the gas number density-temperature
distribution of the full CGM gas mass. Two distinct regions are readily
apparent. In the low-density, high-temperature phase, the gas mass
distribution is similar in both the NoRad and StarRad galaxies. The number
density $n$ and temperature $T$ are related as $n\propto T^{3/2}$, which is an
isentropic line. This happens because gas, after passing through the virial
shock, is adiabatically compressed when flowing to the halo center. The
similar gas mass in this region for the NoRad and StarRad galaxies indicates
that the local stellar radiation does not have a significant impact on this
virialized gas component. On the other hand, the high-density, high-
temperature region shows a more obvious difference between the NoRad and
StarRad galaxies. The gas in the NoRad galaxies can reach temperatures of
$\sim 10^{3.5}\,{\rm K}$, while this phase disappears in the StarRad galaxies.
Note that gas with number density $n\la 0.1\,{\rm cm^{-3}}$ is missing in the
NoRad galaxies (because this is all star-forming), whereas non-starforming gas
denser than this threshold value can still be present in the StarRad galaxies,
since gas with ionization parameter higher than unity is not eligible for
star-formation in our model.

We now divide the gas into inner and outer halo regions and investigate the
gas distribution in the $n-T$ plane separately. The panels of the second and
third row show the gas $n$-$T$ diagram in the inner halo and outer halo,
respectively. Both phases exhibit a wide temperature range. The density in the
inner halo is larger than $10^{-4}\,{\rm cm}^{-3}$, and most of the gas mass
actually resides in the low temperature regime. For the outer halo gas, the
density is significantly lower, and a significant amount of the gas resides at
the virial temperature. Although the total amount of cold gas in this region
is low, it still shows a similar behavior as in the inner region, implying
that local stellar radiation can still affect the gas even outside
$0.3\,R_{\rm vir}$.

The ionization parameter $U$ is a physical quantity that can measure the
strength of the radiation effects. To understand the impact of local stellar
radiation on different phases of the gas, we show the ionization parameter of
the gas at different radii and for different temperatures in Figure 16. The
left panel of Figure 16 gives the gas mass distribution for different
ionization parameters and radii in the StarRad galaxies. In the inner region,
the gas mass experiences mainly $U=0.1-1$, even though the gas is close to the
radiation sources, but this effect is counter acted by the high gas density in
this region. For radiation from the young stars, its heating effect becomes
important when the gas ionization parameter is greater than 1, and the
temperature is lower than $10^{4.5}\,{\rm K}$. This result implies that local
stellar radiation will have a mild heating effect on the inner halo gas. For
the outer halo, it is interesting that the ionization parameter of the gas is
around 10 and is approximately independent of radius. This is because the
radiation flux decreases with $r^{-2}$ in the outer region, and the density
happens to decrease with a similar trend, making the ionization parameter
roughly constant.

The right panel of Figure 16 shows the gas in the $T-U$ diagram. In the
figure, most of the hot gas has a typical ionization parameter $U\sim 1-10$,
while most of the cold gas resides at $U\sim 0.1-10$. Although the ionization
parameter of the hot gas is higher, local stellar radiation does not have a
significant heating effect on it. Thus, the high ionization parameter does not
modify the thermal state of the hot gas. The cold gas, on the other hand, is
much more vulnerable to radiation. The presence of cold gas with high
ionization parameter indicates that there is a significant amount of cold gas
that is affected by local stellar radiation. Compared to the left panel, the
total cold gas mass with a high ionization parameter is higher than that of
the inner halo, implying that a certain amount of radiation-dominated cold gas
resides throughout the entire CGM region.

###  4.3 The abundance of different ion species

Although stellar radiation does not have a significant heating effect on the
hot gas, it may still affect its ionization state. In this section we
therefore investigate the impact of local stellar radiation on different ionic
species. HI, MgII, and OVI are three particularly important ions in CGM
observations. Since they trace the cold gas within the CGM, they probe the
outflow of cold gas from galaxies, as well as the associated metal enrichment.
We thus focus on them in the following.

####  4.3.1 Column density

Figure 17: Stacked column density maps of HI (top row), MgII (middle row) and
OVI (bottom row) gas of differently simulated Auriga galaxies. The left column
gives results from our NoRad simulations, while the right column shows the
corresponding results for our StarRad simulations which included local stellar
radiation. Figure 18: The radial profile of stacked HI, MgII and OVI column
density distributions of our different Auriga galaxies, averaged from $z=0.1$
to $z=0$. The blue star symbols represent the observational results, coming
from Johnson et al. (2015) for HI and OVI column densities, and from Werk et
al. (2013) for observed MgII column densities.

Figure 18 shows stacked images of the HI, OVI and MgII column density in the
galaxies with the same physics model, averaged in time over the outputs from
$z=0.1$ to $z=0$. The three dashed lines in the figure indicate radii equal to
0.1, 0.3 and $1.0\,R_{\rm vir}$. The left column gives the result for the
NoRad galaxies, while the right column is for the StarRad galaxies. Compared
to the NoRad galaxies, the column density of HI gas is systematically lower in
the StarRad galaxies. Since local stellar radiation can easily photoionize the
HI gas, this result is to be expected. For the OVI gas, the column density in
both the NoRad and StarRad galaxies is however quite similar, apart from a
difference in the central column density, which is slightly higher in the
StarRad galaxies than in the NoRad galaxies. For the MgII gas, the NoRad
galaxies feature a uniform distribution within $0.3\,R_{\rm vir}$, whereas the
StarRad galaxies exhibit a more fragmented distribution. This can be explained
by the patchy photoionization effects from the local stellar radiation.

To compare the differences in the column density in the NoRad and StarRad
galaxies more quantitatively, we show the radial profile of the column density
of different ions in Figure 18. The first panel gives the radial profile of
the HI column density. From the figure, we can see that the HI column density
of the NoRad galaxies is higher than that of the StarRad galaxies, which
confirms the conclusion we obtained by eye from Figure 17. The observational
data lies within the 10 and 90 percentile regions for both the NoRad and
StarRad galaxies, but their median value is lower than the observations. This
may be due to selection effects, or it may reflect resolution limitations of
the simulations (Hummels et al., 2019; Peeples et al., 2019; Nelson et al.,
2020; van de Voort et al., 2019).

The second panel shows the OVI column density in both the NoRad and StarRad
galaxies. Similar to our observation in Figure 17, the OVI column densities in
both the NoRad and StarRad galaxies do not display a significant difference.
The OVI column density in the StarRad galaxies is slightly lower than that of
the NoRad galaxies, but the relative difference is only around 20%. The OVI
shows an overall good consistency with observations.

Finally, the third panel gives the MgII column density in both the NoRad and
StarRad galaxies. As for HI, the MgII column density in the StarRad galaxies
is lower than in the NoRad galaxies. MgII also traces the cold gas, and since
local stellar radiation can heat or photoionize the cold gas, it can modify
the abundance of this ion. The observational data resides in the upper range
of the shaded region, indicating a tension between the observations and the
simulations. However, since current instrumentation can only observe luminous
MgII sources, the observational data may be affected by a significant
selection biases. From this point of view, it appears reasonable that the
observational data reside in the upper parts of the shaded region. But it is
also quite plausible that out results are affected by limited resolution,
preventing us from resolving cool outflows accurately (van de Voort et al.,
2019; Hummels et al., 2019; Peeples et al., 2019; Nelson et al., 2020). Also
note that Ji et al. (2020) found that including cosmic rays in the CGM can
increase the MgII column density substantially, which may be another way to
reproduce the observed MgII column densities.

####  4.3.2 Ion mass phase diagram

Figure 19: The temperature–number density distribution for OVI gas (left) and
MgII gas mass (right). In both cases, the upper, middle and bottom rows show
the phase distribution of the CGM gas ($0.1<r/R_{\rm vir}<1$), inner halo gas
($0.1<r/R_{\rm vir}<0.3$) and outer halo gas ($0.3<r/R_{\rm vir}<1$),
respectively, while the left column is for our NoRad simulations and the right
column for the StarRad simulations that include local stellar radiation. The
corresponding distribution for HI gas shows only minor differences, and is
therefore left out for clarity.

Ions can be produced or destroyed by photoionization or collisional
ionization. The collisional ionization process becomes more important as the
temperature increases, while the strength of the photoionization process
depends on the strength of the radiation field, which is characterized by the
ionization parameter. Local stellar radiation can affect the ions both by
photoionizing and by heating the gas to enhance collisional ionization. To
further investigate the effects of local stellar radiation on different phases
of the gas, we show the $n-T$ phase diagram of the HI, OVI, and MgII gas mass
in the whole halo, inner halo, and outer halo in Figure 19. The first panel of
the figure gives the $n-T$ phase diagram of the HI gas mass. The left column
shows the results for the NoRad galaxies, while the right column is for the
StarRad galaxies. Compared to the NoRad galaxies, it can be seen that both,
the area that the HI gas occupies and the total HI gas with the same density
and temperature, are reduced in the StarRad galaxies. The region where $T\la
10^{4}\,{\rm K}$ disappears and the amount of HI gas mass at $T\sim
10^{4}\,{\rm K}$ and $10^{-3}~{}{\rm cm^{-3}}\la n\la 10^{-3}~{}{\rm cm^{-1}}$
also becomes smaller. This result indicates that local stellar radiation
affects the HI gas through both photoionizing and heating the gas. The second
and third columns show the HI gas phase diagram in the inner and outer halo.
The HI gas in both the inner and outer halos is also reduced, indicating once
more that the influence of local stellar radiation on the HI gas spreads
throughout the whole halo.

The second panel of Figure 19 shows the $n-T$ phase diagram of the OVI gas
mass in the NoRad (left) and StarRad (right) galaxies. The OVI gas mainly
resides in the low-density, high-temperature region in both the NoRad and
StarRad galaxies. The difference of the OVI gas mass in this region is
negligible in both types of simulated galaxies, indicating that OVI in this
region is dominated by the collisional ionization process. However, compared
to the NoRad galaxies, the StarRad galaxies have the OVI gas in the region
where $10^{-3}\,{\rm cm^{-3}}\la n\la 10^{-1}\,{\rm cm^{-3}}$ and
$T<10^{5}\,{\rm K}$. Since the typical ionization temperature of the OVI is
$10^{5-6}\,{\rm K}$ if only collisional ionization is considered, the OVI gas
that appears in the low-temperature region in the StarRad galaxies should be
photoionized. This result implies that local stellar radiation can increase
the OVI gas mass. However, we do not see a corresponding increase in the
radial profile of the OVI column density. This is because the OVI gas is still
dominated by hot gas. Since the hot gas has a large covering fraction, the
column density of the OVI should be dominated by the hot gas. The high-
density, low-temperature OVI gas on the other hand is distributed in both the
inner and outer halos. Since high-density, low-temperature gas exists
preferentially in the inner halo, OVI is also more prevalent in the inner
halo.

The third panel of Figure 19 shows the $n-T$ phase diagram of the MgII gas
mass in the NoRad (left) and StarRad (right) galaxies. The MgII gas resides in
a similar region in the phase diagram as the HI gas, but it occupies a smaller
region in the diagram than the HI gas. This is because both the MgII and HI
gas trace the cold gas, but the HI gas can also trace the cold gas accreted
from the IGM, which is not yet polluted with metals. Similar to the HI gas,
MgII is also reduced in the StarRad galaxies compared to the NoRad galaxies.
The reduction is both due to heating and photoionizing effects. The reduction
of the MgII gas in the StarRad galaxies is more significant than that of the
HI gas, because the MgII fraction is more sensitive to temperature than the HI
gas. The second and third panels show the MgII gas phase diagram in the inner
and outer halo. Here the MgII gas mass is reduced in both the inner and outer
halos, similar to the HI gas.

##  5 Comparison with previous work

In the following we briefly discuss three other literature studies that have
investigated the “long-range” impact of local stellar radiation on galaxy
formation and evolution in a comparable fashion as done here (Kannan et al.,
2014; Obreja et al., 2019; Hopkins et al., 2020). In fact, Kannan et al.
(2014) and Obreja et al. (2019) have used a numerical approach similar to the
one we employ in this work, whereas Hopkins et al. (2020) used both the LEBRON
method and the M1 method for radiative transfer. LEBRON is a ray-based scheme
for radiative transfer, which is conceptually still quite close to the method
used in this work, whereas the M1 method refers to the low-order closure used
in a moment-based radiative transfer approach. We now discuss the similarities
and differences between the findings of these three studies and the current
work.

###  5.1 Comparison with Kannan et al. (2014)

Kannan et al. (2014) studied the impact of local stellar radiation on MW-like
galaxies in the MUGS suite, which is a cosmological zoom-in simulation
project. Kannan et al. (2014) employed the same numerical approach as our
current work to calculate the radiation field. However, they identified star
particles younger than 10 Myr as young stars and used the SED of the stellar
population with constant SFR at 10 Myr. In the current work, we identify star
particles borne more recently than 100 Myr as young stars, which implies a
higher radiation flux from young stars.

In Kannan et al. (2014), it was found that both the SFR and the disk gas mass
of the galaxy with local stellar radiation were reduced by a factor of 2. The
amount of gas with temperature $T\sim 10^{4}~{}{\rm K}$ was strongly
suppressed. The peak of the rotation curve at the center was also strongly
suppressed. In general, the effect of local stellar radiation in Kannan et al.
(2014) is similar to but stronger than that found in the present work.

Since we identify more stars as young stars, the radiation flux from young
stars should be larger than in Kannan et al. (2014), and the effects should in
principle also be stronger. However, our results show weaker effects than
those in Kannan et al. (2014). This may be due to the star formation model. In
Kannan et al. (2014), instead of using an effective EOS to determine the
thermal state of the star-forming gas, they directly resolve the cold gas, and
this gas can also be affected by local stellar radiation, potentially
explaining that star formation was more significantly affected by radiation in
their model.

###  5.2 Comparison with Obreja et al. (2019)

Obreja et al. (2019) have implemented a similar scheme of local stellar
radiation as Kannan et al. (2014), but in the NIHAO suite, which is another
cosmological zoom-in simulation project. Their implementation is based on
Kannan et al. (2014) and Kannan et al. (2016). In comparison to our current
work, they also include radiation from the hot gas. The escape fraction in
their work is set to 5% at all frequencies, which is different from Kannan et
al. (2014) and the present work. Such a setup reduces the total amount of
ionizing photons.

In Obreja et al. (2019), very similar conclusions as in Kannan et al. (2014)
and the current work were reached for MW-like galaxies, although the effects
were weaker. Since ionizing photons are more abundant in the current work
because a larger fraction of stars are considered young, and because the
escape fraction is assumed to be higher at certain frequencies, the effects of
local stellar radiation would have been expected to be more significant,
contrary to what we have found.

###  5.3 Comparison with Hopkins et al. (2020)

Hopkins et al. (2020) studied the impacts of different components of radiation
and different methods of radiative transfer on the evolution of galaxies with
stellar masses ranging from $10^{4}$ to $10^{11}{\rm M_{\odot}}$. Unlike in
our work, the radiation pressure of the continuous spectrum was also included.
Their SED of stellar radiation was calculated by STARBURST99, and the SED of
the star particles was evolved with stellar age. However, since their SED of
young stellar radiation was calculated by STARBURST99, they lacked a certain
amount of ionized photons coming from the SNRs and HXMBs. They used a ray-
based LEBRON method and a momentum-based M1 method to calculate radiative
transfer. The LEBRON method is similar to the method we implement in this
work, but they also included local extinction calculated with the Sobolev
approximation as well as radiation pressure.

We focus on a comparison with results for their m12i and m12m galaxies, which
are MW-like galaxies. For the time evolution of the stellar mass of m12i and
m12m, different from the result presented in this work and Obreja et al.
(2019), the local stellar radiation did not have an impact on the stellar
mass, both for the LERBON and M1 methods. For the gas in the $n-T$ phase
diagram of m12i, the total amount of gas mass lower than $10^{4}\,{\rm K}$ was
reduced, which is qualitatively consistent with the results of Obreja et al.
(2019) and our work. However, the degree of gas reduction in Hopkins et al.
(2020) is nevertheless weaker than found by Obreja et al. (2019) and in our
current work.

In Hopkins et al. (2020), there may be two reasons for the quite weak impact
of the stellar radiation feedback. First, since Hopkins et al. (2020) resolve
GMCs directly, they do not need to assume an escape fraction of the radiation
from the star-forming region. The real escape fraction of radiation from the
star-forming region may be lower than what we have used in the current work.
Second, as mentioned above, the radiation from star-forming regions lacks some
ionized photons from SNRs and HXMBs, which may make the effects of local
stellar radiation become weaker.

##  6 Conclusions

In this work, we have studied the impacts of local stellar radiation on MW-
like galaxies taken from the Auriga project, which are cosmological zoom-in
simulations performed with the moving-mesh code AREPO. Our numerical
implementation for calculating the local stellar radiation field closely
follows the method presented in Kannan et al. (2014), where the local stellar
radiation field is calculated in an optically thing approximation by using the
octree method. In particular, we assume that a certain fraction of the
radiation emanating from a star-forming region will be absorbed internally,
but once the photons have escaped to outside the star forming region the
radiation from old stars is not absorbed when it propagates. The spectrum of
the local stellar radiation is obtained using a BPASS v2.3 model, which is a
stellar population synthesis code that accounts also for the evolution of
binary stars. Since the simulation results exhibit stochastic fluctuations and
significant system-to-system variations, we compute five Auriga galaxies and
simulate them each with (‘StarRad’ simulations) and without stellar radiation
(‘NoRad’ models), allowing us to obtain robust average trends for their
differences. Our main findings can be summarized as follows:

  * •

By comparing the differences in the time evolution of different baryon
components in the NoRad and StarRad galaxies, we find that local stellar
radiation does hardly affect the total baryonic mass in the five simulated
systems, implying that it cannot prevent or modify the gas accretion into the
halo. However, the growth of the central BH mass is enhanced, while the growth
of the stellar mass is suppressed. We find that the BH mass in our StarRad
galaxies is $\sim 40\%$ higher than in NoRad galaxies, whereas the stellar
mass is $\sim 25\%$ lower. The associated suppression of star formation occurs
at all redshifts, but the enhancement of the BH mass growth happens
predominantly at high redshift. In fact, at low redshift, the inclusion of
local stellar radiation leads to a reduction of BH growth.

  * •

The suppression of star formation by local stellar radiation happens because
some gas is prevented from forming stars by being heated slightly above $\sim
10^{4}\,{\rm K}$, which is the equilibrium temperature in the presence of the
local stellar radiation field, but it is not because gas is ejected outside of
the galaxy. The enhancement of the BH growth can be partially understood as a
consequence of converting less gas to stars, leaving more gaseous fuel for
flowing into the galaxy center at high redshift. However, at low redshift,
this channel is becoming less effective, which may be in part related to the
higher specific angular momentum of the corresponding gas. In addition, the
heating effect of the local stellar radiation trends to increase the gas
temperature and decreases the density near the BH, slowing BH accretion.

  * •

As a consequence of the modified star formation, maps of the stellar light are
redder for our StarRad galaxies than for the NoRad galaxies at low redshift.
By comparing the average radial profiles of the young and old stellar surface
density, we find that the surface density of old stars in the StarRad galaxies
is only slightly lower than in the StarRad galaxies, whereas the surface
density of young stars in the StarRad galaxies is significantly lower than in
NoRad galaxies. In general, the half-mass radius of the total stars in our
StarRad galaxies is slightly larger than that in NoRad galaxies, while there
is no obvious trend for the young stars.

  * •

The surface density distribution of the total gas in both NoRad and StarRad
galaxies is similar, while the surface density distribution of the HI gas in
StarRad galaxies is significantly lower than the NoRad galaxies, which
confirms that local stellar radiation does not significantly affect gas inflow
onto the disk, but it can heat and/or ionize the cold gas in the disk.

  * •

By comparing the rotation curve of five NoRad galaxies and StarRad galaxies,
we find that local stellar radiation can suppress the peak of the rotation
curve, implying that this feedback can contribute to alleviating the
overcooling problem, albeit only weakly. This influence on the rotation curve
is due to the reduction of old stellar mass within 5 kpc, and originates in
the suppression effect of star formation at early times.

  * •

The young stellar component is the dominated factor for the feedback influence
of the local stellar radiation. The ionization parameter of the young stellar
radiation for the cold CGM gas is $\sim 0.01-10$, whereas it is $\sim 10$ for
the hot CGM gas. Local stellar radiation can significantly affect the thermal
state of cold ($T<10^{4}\,{\rm K}$) and cool ($10^{4}\,{\rm K}<T<10^{5}~{}{\rm
K}$) CGM gas. As a result, the cold CGM gas is significantly reduced even
outside $0.3\,R_{\rm vir}$ due to the heating effect of local stellar
radiation.

  * •

Local stellar radiation can also reduce the HI and MgII column density of CGM
gas within $0.4\,R_{\rm vir}$, while it does not have a significant impact on
the OVI column density. This result arises from a combination of heating and
photoionization effects due to the local stellar radiation. For HI and MgII
that trace the cold gas, the local stellar radiation can reduce the column
density both by reducing the total amount of cold gas by heating, and by
reducing the ion fraction for the gas at certain temperatures by
photoionization. Since the cold gas mainly resides in the inner region, HI and
MgII experience a more significant reduction in the inner region. While the
OVI gas traces the hot gas, local stellar radiation is too soft to heat and
photoionize the hot gas. Thus, the OVI column density does not show
significant differences between the NoRad and StarRad galaxies.

The findings presented above demonstrate that local stellar radiation acts as
a mild preventive feedback and clearly has important quantitative impact on
the formation and evolution of MW-like galaxies, consistent with earlier
findings. In particular, this radiative feedback needs to be considered when
calculating the physical properties of CGM gas. This form of radiative
feedback should therefore be included by default in future cosmological
simulations of galaxy formation, something that has largely not been done thus
far.

We note, however, that the model for local radiation considered here is not
complete and still calls for many improvements in future work. With respect to
the radiation sources, Kannan et al. (2016) has shown that the radiation from
the hot gas can suppress the star formation of satellites in the massive halo,
while Oppenheimer et al. (2018) and Nelson et al. (2018) show that the flicker
of AGN radiation can help photoionize the OVI in CGM. Both radiation sources
have not been included in the current work. Also, in our numerical
implementation of local radiation, the self-shielding of high-density clouds
is not yet included, and the critical value of the ionization parameter we
used for stopping star formation is quite ad-hoc and poorly constrained. The
quantitative impacts of the local radiation field we have found here provide a
strong motivation to address these points in future work.

## Acknowledgements

The authors acknowledge support by a Max Planck Partner group between the
Shanghai Astronomical Observatory and the Max Planck Institute for
Astrophysics (MPA). Computations were performed on the Freya compute cluster
of MPA, operated by the Max Planck Computing and Data Facility.

## Data Availability

The data of this study are available from the corresponding author upon
reasonable request.

## References

  * Anderson et al. (2013) Anderson M. E., Bregman J. N., Dai X., 2013, ApJ, 762, 106
  * Babul & Rees (1992) Babul A., Rees M. J., 1992, MNRAS, 255, 346
  * Balogh et al. (2001) Balogh M. L., Pearce F. R., Bower R. G., Kay S. T., 2001, MNRAS, 326, 1228
  * Benson et al. (2002) Benson A. J., Lacey C. G., Baugh C. M., Cole S., Frenk C. S., 2002, MNRAS, 333, 156
  * Borrow et al. (2023) Borrow J., Kannan R., Garaldi E., Smith A., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2023, MNRAS, 525, 5932
  * Cantalupo (2010) Cantalupo S., 2010, MNRAS, 403, L16
  * Cen & Ostriker (1992) Cen R., Ostriker J., 1992, ApJ, 393, 22
  * Cen & Ostriker (1993) Cen R., Ostriker J. P., 1993, ApJ, 417, 404
  * Cerviño et al. (2002) Cerviño M., Mas-Hesse J. M., Kunth D., 2002, A&A, 392, 19
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Efstathiou (1992) Efstathiou G., 1992, MNRAS, 256, 43P
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Australia, 34, e058
  * Faucher-Giguère (2020) Faucher-Giguère C.-A., 2020, MNRAS, 493, 1614
  * Faucher-Giguère et al. (2009) Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416
  * Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
  * Fielding et al. (2020) Fielding D. B., et al., 2020, ApJ, 903, 32
  * Finlator et al. (2011) Finlator K., Davé R., Özel F., 2011, ApJ, 743, 169
  * Governato et al. (2004) Governato F., et al., 2004, ApJ, 607, 688
  * Grand et al. (2017) Grand R. J. J., et al., 2017, MNRAS, 467, 179
  * Grand et al. (2024) Grand R. J. J., Fragkoudi F., Gómez F. A., Jenkins A., Marinacci F., Pakmor R., Springel V., 2024, arXiv e-prints, p. arXiv:2401.08750
  * Haardt & Madau (1996) Haardt F., Madau P., 1996, ApJ, 461, 20
  * Haardt & Madau (2012) Haardt F., Madau P., 2012, ApJ, 746, 125
  * Haiman (2016) Haiman Z., 2016, in Mesinger A., ed., Astrophysics and Space Science Library Vol. 423, Understanding the Epoch of Cosmic Reionization: Challenges and Progress. p. 1 (arXiv:1511.01125)
  * Haiman et al. (1996) Haiman Z., Thoul A. A., Loeb A., 1996, ApJ, 464, 523
  * Hopkins et al. (2012) Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3522
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Hopkins et al. (2020) Hopkins P. F., Grudić M. Y., Wetzel A., Kereš D., Faucher-Giguère C.-A., Ma X., Murray N., Butcher N., 2020, MNRAS, 491, 3702
  * Hummels et al. (2019) Hummels C. B., et al., 2019, ApJ, 882, 156
  * Ji et al. (2020) Ji S., et al., 2020, MNRAS, 496, 4221
  * Johnson et al. (2015) Johnson S. D., Chen H.-W., Mulchaey J. S., 2015, MNRAS, 449, 3263
  * Kannan et al. (2014) Kannan R., et al., 2014, MNRAS, 437, 2882
  * Kannan et al. (2016) Kannan R., Vogelsberger M., Stinson G. S., Hennawi J. F., Marinacci F., Springel V., Macciò A. V., 2016, MNRAS, 458, 2516
  * Katz et al. (1996) Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
  * Loeb & Barkana (2001) Loeb A., Barkana R., 2001, ARA&A, 39, 19
  * Navarro & Benz (1991) Navarro J. F., Benz W., 1991, ApJ, 380, 320
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 477, 450
  * Nelson et al. (2020) Nelson D., et al., 2020, MNRAS, 498, 2391
  * Obreja et al. (2019) Obreja A., Macciò A. V., Moster B., Udrescu S. M., Buck T., Kannan R., Dutton A. A., Blank M., 2019, MNRAS, 490, 1518
  * Okamoto et al. (2008) Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920
  * Oppenheimer et al. (2018) Oppenheimer B. D., Segers M., Schaye J., Richings A. J., Crain R. A., 2018, MNRAS, 474, 4740
  * Peeples et al. (2019) Peeples M. S., et al., 2019, ApJ, 873, 129
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Qu & Bregman (2018) Qu Z., Bregman J. N., 2018, ApJ, 856, 5
  * Quinn et al. (1996) Quinn T., Katz N., Efstathiou G., 1996, MNRAS, 278, L49
  * Raymond et al. (1975) Raymond J. C., Smith B. W., Cox D. P., 1975, in Bulletin of the American Astronomical Society. p. 538
  * Rees & Ostriker (1977) Rees M. J., Ostriker J. P., 1977, MNRAS, 179, 541
  * Roškar et al. (2014) Roškar R., Teyssier R., Agertz O., Wetzstein M., Moore B., 2014, MNRAS, 444, 2837
  * Scannapieco et al. (2012) Scannapieco C., et al., 2012, MNRAS, 423, 1726
  * Shapiro & Giroux (1987) Shapiro P. R., Giroux M. L., 1987, ApJ, 321, L107
  * Silk (1977) Silk J., 1977, ApJ, 211, 638
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Stern et al. (2016) Stern J., Hennawi J. F., Prochaska J. X., Werk J. K., 2016, ApJ, 830, 87
  * Stinson et al. (2010) Stinson G. S., Bailin J., Couchman H., Wadsley J., Shen S., Nickerson S., Brook C., Quinn T., 2010, MNRAS, 408, 812
  * Strawn et al. (2021) Strawn C., et al., 2021, MNRAS, 501, 4948
  * Suresh et al. (2017) Suresh J., Rubin K. H. R., Kannan R., Werk J. K., Hernquist L., Vogelsberger M., 2017, MNRAS, 465, 2966
  * Thoul & Weinberg (1996) Thoul A. A., Weinberg D. H., 1996, ApJ, 465, 608
  * Tumlinson et al. (2011) Tumlinson J., et al., 2011, Science, 334, 948
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Wadsley et al. (2004) Wadsley J. W., Stadel J., Quinn T., 2004, New Astron., 9, 137
  * Wang et al. (2009) Wang L., Mao J., Xiang S., Yuan Y. F., 2009, A&A, 494, 817
  * Wang et al. (2015) Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83
  * Weinberg et al. (1997) Weinberg D. H., Hernquist L., Katz N., 1997, ApJ, 477, 8
  * Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056
  * Werk et al. (2013) Werk J. K., Prochaska J. X., Thom C., Tumlinson J., Tripp T. M., O’Meara J. M., Peeples M. S., 2013, ApJS, 204, 17
  * Werk et al. (2014) Werk J. K., et al., 2014, ApJ, 792, 8
  * White & Frenk (1991) White S. D. M., Frenk C. S., 1991, ApJ, 379, 52
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Woods (2015) Woods R., 2015, PhD thesis, McMaster University, Canada
  * Zinger et al. (2020) Zinger E., et al., 2020, MNRAS, 499, 768
  * van de Voort et al. (2019) van de Voort F., Springel V., Mandelker N., van den Bosch F. C., Pakmor R., 2019, MNRAS, 482, L85