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

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 with our current model that local stellar radiation can modify gas cooling in the CGM and thus suppress star formation and the surface densities of young stars and H i gas, while having little impact on the total gas content. 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 ultraviolet background and dominates the radiation field within the virial halo at all redshifts. Nevertheless, we find that the local stellar radiation, as implemented in the current study, has overall little impact on the radial density and temperature profile of the CGM gas. However, for the ion species H i and Mg ii, the column densities within |$\sim 0.3\, R_{\rm vir}$| are reduced, while the O vi column density is hardly impacted by the radiation field due to a lack of soft X-ray components in our current model. Additional effects can be expected from the radiation of the central active galactic nucleus during phases of quasar activity and from soft X-ray sources, which have not yet been included in the simulations of the present study.

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 (Rees & Ostriker 1977; Silk 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 centre 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 haloes.

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, 2012; Faucher-Giguère et al. 2009; 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 haloes (Babul & Rees 1992; Cen & Ostriker 1992; Efstathiou 1992; Haiman, Thoul & Loeb 1996; Quinn, Katz & Efstathiou 1996). As a result, star formation in these small haloes is suppressed, as the gas is no longer able to efficiently cool in them (Wang et al. 2009; Finlator, Davé & Özel 2011; Borrow et al. 2023).

Current cosmological simulations commonly implement the UVB as a spatially uniform, time-dependent radiation field (e.g. Katz, Weinberg & Hernquist 1996; Quinn et al. 1996; Thoul & Weinberg 1996; Weinberg, Hernquist & Katz 1997; Vogelsberger et al. 2014; Crain et al. 2015; Pillepich et al. 2018), with a time-dependent spectrum as modelled in several works (Haardt & Madau 1996, 2012; Faucher-Giguère et al. 2009; 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 photionization 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 DM haloes (Cen & Ostriker 1993; Benson et al. 2002; Okamoto, Gao & Theuns 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 (Quinn et al. 1996; Thoul & Weinberg 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 haloes of MW-like galaxies. However, at the centre 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 centre 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, Quataert & Murray (2012) performed several galaxy formation simulations with different halo masses, considering both the short- and long-range radiation pressures as well as the photoionization heating from H ii 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 disc. However, these works did not take the photoionization heating caused by the certain stellar populations into account. Ceverino et al. (2014) have performed a cosmological zoom-in simulation of the progenitor of an MW-sized galaxy with radiation pressure, photoionization, and photoheating calculated by a detailed stellar population description and found that radiation feedback is helpful for reducing and flattening the rotation curve. However, long-range effects of the local stellar radiation were not included in the study.

Cantalupo (2010) investigated the impact of X-ray radiation emitted by hot winds from young stars, supernova (SN) 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, Stadel & Quinn 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 per cent 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; Nelson et al. 2018; Oppenheimer et al. 2018; Qu & Bregman 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 O vi column density with and without the galaxies’ central stellar radiation and found that it helps to increase the O vi column density in the inner region. Oppenheimer et al. (2018) consider the flickering AGN (active galactic nucleus) radiation and non-equilibrium ionization of the CGM gas in Evolution and Assembly of GaLaxies and their Environments zoom simulations (EAGLE, Crain et al. 2015), finding that the AGN radiation can enhance the O vi column density allowing the COS (Cosmic Origins Spectrograph)-haloes observations (Tumlinson et al. 2011; Werk et al. 2014) to be matched. Finally, Nelson et al. (2018) studied the O vi column density in the IllustrisTNG simulation (TNG, Pillepich et al. 2018; Weinberger et al. 2018) with several variations of the feedback prescriptions. They found that the O vi column density becomes lower when turning off the BH (black hole) radiative feedback compared to a simulation with the fiducial TNG model, for halo masses in the range |$10^{11.5}\lt M_{\rm halo}/{\rm {\rm M}_{\odot }}\lt 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 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 disc 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, discuss the potential effects of some radiation sources that are not considered in this work, and test the non-linear coupling of the interaction between AGN feedback and local stellar radiation. 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 were 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 TNG (Pillepich et al. 2018; Weinberger et al. 2018). This in turn employs a star formation and stellar feedback model which is based on the two-phase interstellar medium (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 equal or less than |$\sim 10^4~{\rm K}$|⁠. 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 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 multiphase 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).

BHs are modelled as sink particles in the centres 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 BH 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 colour distribution at |$z = 0$| consistent with observations.

The Auriga project selected 30 haloes with mass |$M_{\rm halo}=1-2\times 10^{12}\, {\rm {\rm M}_{\odot }}$| from the DM-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 XVI (2014) analysis.

In this work, we select a subset of five typical haloes from the Auriga galaxy set to study the impacts of the local stellar radiation on the galaxy and CGM properties. The selected haloes correspond to Au3, Au5, Au6, Au15, and 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 haloes 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 haloes at |$z=0$|⁠, with and without the local stellar radiation, are given in Table 1.

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 two 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 SFR.

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

The escape fraction at the Lyman limit frequency |$f^{\rm LL}_{\rm esc}$| is set to 5 per cent. Here, only absorption of the neutral hydrogen is taken into account, thus |$\tau _{\nu }=\sigma _{\nu } N{\rm (H^{0}})$|⁠. |$\sigma _{\nu }$| is the cross-section of neutral hydrogen and |$N{\rm (H^{0}})$| is the column density of neutral 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 follow the description presented in Cerviño, Mas-Hesse & Kunth (2002). They reported that the total energy released in the EINSTEIN band (⁠|$0.1-3.4 \, {\rm keV}$|⁠) of an individual SNe is |$\sim 1.07\times 10^{40} {\rm erg}$|⁠. Based on this value, we calculate the SNe event rate from a young stellar population with a Chabrier (2003) initial mass function and constant |$1~{\rm {\rm M}_{\odot }~yr^{-1}}$| SFR and normalize the SED accordingly. The shape of the SED is as follows:

where |$f_{\rm RS}(kT)$| is the Raymond–Smith model (Raymond & Smith 1977).

The HMXB emission in this work is assumed to follow a |$\Gamma =2$| power law. Based on Anderson, Bregman & Dai (2013), the total SED of HMXB is in the range from 0.5 to 8.0 keV band, and the scaling relation of SFR and total X-ray luminosity is:

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 BHs 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. Following the description in Anderson et al. (2013), the SED is in the range from 0.3 to 8 keV, and the LMXB luminosity as a function of total stellar mass scales as

The intensities of the different components mentioned above are shown in Fig. 1. The radiation intensity is calculated by assuming an SFR or |$1~{\rm {\rm M}_{\odot }~yr^{-1}}$| for young stellar populations, |$M_{\mathrm{star}}=10^{11}~{\rm {\rm M}_{\odot }}$| for the old stellar population, and adopting a distance of 200 kpc away from the galaxy.

The hard truncation of the SED in the LMXB and HMXB regime is not a particularly accurate physical model since there are certainly some photons emitted below 0.3 keV. For example, Upton Sanderbeck et al. (2018) use an empirical extrapolation of the HMXB SED below 0.3 keV based on a power-law index of |$-0.7$|⁠. The corresponding emission is uncertain and missing in our present analysis. We will examine the potential impact of these lacking photons in future work.

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 EoR.

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

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

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 (5)–(9) 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

where

|$\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.

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

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 Fig. 2.

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 6D 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\gt 10$|⁠.

For the table of the cooling/heating function, the gas number density ranges from |$10^{-8}\, \mathrm{ to}\, 10^{2}~{\rm cm^{-3}}$| in spacing of 1 dex in log space. The temperature ranges from |$10^{1}\, \mathrm{ to}\, 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}\, \mathrm{ to}\, 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  per cent in 99  per cent.

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, that is, 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 {\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

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.

3 THE IMPACTS OF LOCAL IONIZATION FEEDBACK ON GALAXY PROPERTIES

In this section, we analyse and discuss the differences in the time evolution, in the properties of gas and stellar discs, 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 analysing 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 Fig. 2 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}$|⁠, and |$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}\gtrsim 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

Fig. 3 shows the time evolution of the total baryonic, the BH, the stellar, and the gas mass in five selected haloes from |$z=5.5$| to 0. The first panel gives the evolution of the baryonic mass and the relative difference of five simulated haloes with the NoRad and StarRad galaxies. In all but the Au3 galaxies, the relative differences between the NoRad and StarRad galaxies are |$\lt 20~{{\ \rm per\ cent}}$| and do not have a clear trend, implying that the current implementation of local photoionization sources has little impact on the gas accretion of the haloes. 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 Fig. 3 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 Fig. 3, 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 non-linearity of the fluid behaviour, combined with the local stellar radiation, the subsequent evolution may experience strong differences. A detailed analysis of the non-linearly combined influence of AGN feedback and local stellar radiation during and after the major merger of the Au3 galaxy is presented in Section 5.3.

To further confirm this conclusion, we present the time evolution of the 100 Myr-averaged SFR in Fig. 4. We focus on the Au3 and Au16 galaxies. The first panel of Fig. 3 shows that the time evolution of the baryonic mass is similar in Au3 and Au16 at |$z\gtrsim 1$|⁠. However, at |$z\sim 1$| the growth rate of the baryonic mass rapidly decreases in the Au3 StarRad model. However, we can observe in Fig. 4 that the SFR in Au3 at |$z\sim 1$| does not show significant differences compared to Au16. This confirms our conclusions mentioned above that our model of local stellar radiation does not have a significant impact on the gas accrection of the haloes. The differences of the evolution of the baryonic mass in Au3 at |$z\lesssim 1$| is due to the AGN feedback during the major merger events.

The second panel of Fig. 3 shows the growth of the BH mass and the relative differences between the NoRad and StarRad galaxies in all simulated haloes. 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 per cent in the first 2 Gyr, but then becomes lower at later times. This is also reflected in the right panel of Fig. 5, 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 enhance 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~{{\ \rm per\ cent}}$| 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  per cent, 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.1.

Note, however, that these trends are not always very clear at |$z\ge 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 disc has formed. The results in Kannan et al. (2014), where star formation was found to be suppressed by |$\sim 30~{{\ \rm per\ cent}}$| in a disc 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 discs 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.

We further investigated the stellar age distribution of simulated galaxies at |$z=0$|⁠. The left panel of Fig. 5 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 |$\sim 0.7$| is systematically lower in the StarRad galaxies than in the NoRad galaxies. At |$z\gtrsim 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 Fig. 3. At |$z\lesssim 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 SFR 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 Fig. 3 shows the time evolution of the total gas mass and the relative differences between the NoRad and StarRad galaxies in all simulated haloes. 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 Fig. 6. 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 Fig. 6 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 Fig. 6 shows the time evolution of the star-forming gas mass within |$0.1 R_{\rm vir}$|⁠. Unlike the time evolution of the total gas mass shown in the first panel of Fig. 6, 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 Fig. 5, 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 Fig. 6 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.

Fig. 7 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 Disc 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 analyse the impact of local ionization feedback on both stellar and gas disc properties |$z=0$|⁠.

3.2.1 Stellar disc

Fig. 8 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 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 centre. Except for the surface brightness, both the red and blue lights of NoRad galaxies are also more extended. These results imply that local stellar radiation can also affect the stellar disc 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 1D surface density of total, old, and young stars in Fig. 9. Young 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 0. The solid line represents the median value of the surface density, while the shaded area in the figure represents the 10–90 percentile range of the distribution of measurements.

The first panel of Fig. 9 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 centre, but becomes larger outside 10 kpc. For determining the 10–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 Fig. 9, 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 centre 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 Fig. 9. 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 Fig. 8. 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 disc-like, these results imply that local stellar radiation is comparatively effective in suppressing the star formation in disc-like configurations.

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 Sérsic 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 Sérsic 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.

Fig. 10 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\lt 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 Fig. 5 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\lt 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 disc galaxy formation.

The second panel of Fig. 10 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–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 disc

Fig. 11 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 disc in the NoRad and StarRad galaxies highlight that they have similar morphology and size. The disc 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 Fig. 6, 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 discs 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 discs in the NoRad and StarRad galaxies, we compare the radial profile of the total gas surface density in the first panel of Fig. 12. The profiles are stacked and averaged from |$z=0.1$| to 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.

Fig. 13 shows the H i gas surface density of different simulated galaxies at |$z=0$|⁠. Similarly to Fig. 11, 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 H i gas surface density exhibits significant differences between the NoRad and StarRad cases. The face-on view indicates that the H i disc of the NoRad galaxies features a relatively smooth distribution. However, the H i disc in the StarRad galaxies appears clumpy. In addition, the H i disc in the StarRad galaxies has a smaller size. From the edge-on view, we see that the H i disc in the StarRad galaxies has also a lower disc height. These results are consistent with the findings in Obreja et al. (2019).

These differences are quantified in the stacked radial profiles of the H i gas surface density shown in the second panel of Fig. 12. From the figure, we see that the H i 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 {\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 H i disc surface density drops rapidly at |$\sim 25-35\, {\rm kpc}$| in the StarRad galaxies. These differences in H i morphology and surface density between NoRad and StarRad galaxies indicate that local stellar radiation can effectively ionize the H i gas in the disc, making the H i disc smaller, clumpier, and giving it a lower surface density.

The final bottom panel of Fig. 12 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 H i gas disc. This lower star-forming gas surface density of the StarRad galaxies can also explain the lower disc surface density of young stars seen in Fig. 9.

3.3 Rotation curve

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 haloes 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.

Fig. 14 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  per cent. 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 equation of state (EOS), suppressing resolved multiphase 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 centre 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 CIRCUMGALACTIC MEDIUM

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; Nelson et al. 2018; Oppenheimer et al. 2018). However, the radiation from nearby stars may play a significant role for the CGM region.

To illustrate this point, Fig. 15 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\lt 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 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 per cent, 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 analyse the influence of the local stellar radiation on the CGM properties.

4.1 Radial profiles of gas properties

We begin by studying the impact of local stellar radiation on the global structure of the gaseous haloes. Fig. 16 gives the stacked radial profiles of the gas number density, temperature, and metallicity, in both the NoRad and StarRad galaxies. Only non-star-forming gas is included, that is, 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 0 to eliminate temporal fluctuations. The solid line represents the median value of the individual measurements, while the shaded area indicates the 10  per cent–90  per cent percentiles.

The number density profile, shown in the first panel of Fig. 16, 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 disc. 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 Fig. 12, 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 Fig. 16 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. Since the virial temperature of our simulated haloes is |$\sim 10^6\, {\rm K}$|⁠, the coolants are dominated by |${\rm Ne^{\rm 5+}}$| and |${\rm Fe^{\rm 8+}}$| (Wiersma, Schaye & Smith 2009). Note that to photoionize these ions, the energy of the photons should be |$\sim 100-300\, {\rm eV}$|⁠, which are lacking in our model of the local stellar radiation. Nevertheless, there are still some indications in our 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 Fig. 16 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 centre of the galaxy. Recall that Fig. 7 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

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\lt r/R_{\rm vir}\lt 0.3$|⁠) and the outer halo (⁠|$0.3\lt r/R_{\rm vir}\lt 1$|⁠). Fig. 17 shows the gas number density versus 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 1 to obtain statistically more robust results. The star-forming gas of the ISM is excluded.

In the top row of Fig. 17, 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 centre. The similar gas mass in this region for the NoRad and StarRad galaxies indicates that our model of 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\lesssim 0.1\, {\rm cm^{-3}}$| is missing in the NoRad galaxies (because this is all star-forming), whereas non-star-forming 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 rows show the gas n–T diagram in the inner and outer haloes, 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 behaviour 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 Fig. 18. The left panel of Fig. 18 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 Fig. 18 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 our model of local 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. H i, Mg ii, and O vi 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

Fig. 19 shows stacked images of the H i, O vi, and Mg ii column density in the galaxies with the same physics model, averaged in time over the outputs from |$z=0.1$| to 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 H i gas is systematically lower in the StarRad galaxies. Since local stellar radiation can easily photoionize the H i gas, this result is to be expected. For the O vi 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 Mg ii 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 Fig. 19. The first panel gives the radial profile of the H i column density. From the figure, we can see that the H i column density of the NoRad galaxies is higher than that of the StarRad galaxies, which confirms the conclusion we obtained by eye from Fig. 20. Inside |$\sim$| 0.3 |$R_{\rm vir}$|⁠, the simulation matches the observations well. However, recent 21-cm observations (Das et al. 2024) show that the H i gas column density within 0.5 |$R_{\rm vir}$| declines from |$\sim 5\times 10^{18}~{\rm cm^{-2}}$| at |$r\sim 0.1R_{\rm vir}$| to |$\sim 10^{17}~{\rm cm^{-2}}$| at |$r\sim 0.5R_{\rm vir}$|⁠, which is higher than in our simulations. In addition, outside |$\sim$| 0.3 |$R_{\rm vir}$|⁠, 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. Both lower estimates of the H i column density may be due to selection effects, or it may reflect resolution limitations of the simulations (Hummels et al. 2019; Peeples et al. 2019; van de Voort et al. 2019; Nelson et al. 2020).

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

Finally, the third panel gives the Mg ii column density in both the NoRad and StarRad galaxies. As for H i, the Mg ii column density in the StarRad galaxies is lower than in the NoRad galaxies. Mg ii 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 reside 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 Mg ii 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 (Hummels et al. 2019; Peeples et al. 2019; van de Voort 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 Mg ii column density substantially, which may be another way to reproduce the observed Mg ii column densities.

4.3.2 Ion mass phase diagram

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 H i, O vi, and Mg ii gas mass in the whole, inner, and outer haloes in Fig. 21. The first panel of the figure gives the |$n-T$| phase diagram of the H i 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 H i gas occupies and the total H i gas with the same density and temperature, are reduced in the StarRad galaxies. The region where |$T\lesssim 10^{4}\, {\rm K}$| disappears and the amount of H i gas mass at |$T\sim 10^4\, {\rm K}$| and |$10^{-3}~{\rm cm^{-3}}\lesssim n\lesssim 10^{-3}~{\rm cm^{-1}}$| also becomes smaller. This result indicates that local stellar radiation affects the H i gas through both photoionizing and heating the gas. The second and third columns show the H i gas phase diagram in the inner and outer halo. The H i gas in both the inner and outer haloes is also reduced, indicating once more that the influence of local stellar radiation on the H i gas spreads throughout the whole halo.

The second panel of Fig. 21 shows the |$n-T$| phase diagram of the O vi gas mass in the NoRad (left) and StarRad (right) galaxies. The O vi gas mainly resides in the low-density, high-temperature region in both the NoRad and StarRad galaxies. The difference of the O vi gas mass in this region is negligible in both types of simulated galaxies, indicating that O vi in this region is dominated by the collisional ionization process. However, compared to the NoRad galaxies, the StarRad galaxies have the O vi gas in the region where |$10^{-3}\, {\rm cm^{-3}}\lesssim n\lesssim 10^{-1}\, {\rm cm^{-3}}$| and |$T\lt 10^{5}\, {\rm K}$|⁠. Since the typical ionization temperature of the O vi is |$10^{5-6}\, {\rm K}$| if only collisional ionization is considered, the O vi 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 O vi gas mass. However, we do not see a corresponding increase in the radial profile of the O vi column density. This is because the O vi gas is still dominated by hot gas. Since the hot gas has a large covering fraction, the column density of the O vi should be dominated by the hot gas. The high-density, low-temperature O vi gas on the other hand is distributed in both the inner and outer haloes. Since high-density, low-temperature gas exists preferentially in the inner halo, O vi is also more prevalent in the inner halo.

The third panel of Fig. 21 shows the |$n-T$| phase diagram of the Mg ii gas mass in the NoRad (left) and StarRad (right) galaxies. The Mg ii gas resides in a similar region in the phase diagram as the H i gas, but it occupies a smaller region in the diagram than the H i gas. This is because both the Mg ii and H i gas trace the cold gas, but the H i gas can also trace the cold gas accreted from the IGM, which is not yet polluted with metals. Similar to the H i gas, Mg ii 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 Mg ii gas in the StarRad galaxies is more significant than that of the H i gas, because the Mg ii fraction is more sensitive to temperature than the H i gas. The second and third panels show the Mg ii gas phase diagram in the inner and outer halo. Here, the Mg ii gas mass is reduced in both the inner and outer haloes, similar to the H i gas.

5 DISCUSSION

5.1 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 (Locally Extincted Background Radiationin Optically ThinNetworks) 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.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 disc 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 centre 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 this 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.1.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, 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  per cent at all frequencies, which is different from Kannan et al. (2014) and this 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.1.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 {\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.

5.2 Other possible radiation sources

There are various sources of extreme ultraviolet (EUV) and soft X-ray photons that could potentially impact the properties of galaxies and the CGM, as noted in recent studies (Upton Sanderbeck et al. 2018). Here, we focus on two potentially significant contributors: super soft sources (SSS) and the warm-hot interstellar medium (ISM), which emits EUV and soft X-ray photons.

SSS arise from radiation emitted during accretion onto white dwarfs, while the warm-hot ISM emission stems from radiation generated by shocked ISM due to SNe or stellar winds from OB stars. We utilize the SEDs presented in Upton Sanderbeck et al. (2018) and compare them with those used in our present work. Given the uncertainty in these SEDs, we incorporate two different models for SSS (with temperatures of |$3\times 10^5$| and |$5\times 10^5$| K) and consider the upper and lower bounds of radiation intensity from the warm-hot ISM as reported in Upton Sanderbeck et al. (2018).

Fig. 22 displays the intensities of SSS and warm-hot ISM (with |$f_{\rm esc}f_{\rm rad}=0.01$| and 0.1) at a distance of 200 kpc from the radiation sources, assuming an SFR of 1 |${\rm {\rm M}_{\odot }~yr^{-1}}$|⁠. These intensities are compared with those of the young stellar population (SFR = 1 |${\rm {\rm M}_{\odot }yr^{-1}}$|⁠), the old stellar population (⁠|${M_{\mathrm{star}}=10^{11}\mathrm{{\rm M}_{\odot }}}$|⁠), and the UVB at |$z=0$| for context.

Regarding the intensity from the warm-hot ISM compared to the young stellar population, we find that the former is lower, even at the upper limit of the estimated radiation from the warm-hot ISM. This suggests that the warm-hot ISM does not significantly add to the radiation observed from the young stellar population.

For the SSS intensity, it is comparable to that of the young stellar population. However, since SSS radiation peaks at approximately 1 Gyr after star formation, there is minimal overlap in time with the radiation from the young stellar population, which typically dominates during the first 100 Myr after star formation. None the less, SSS could potentially extend the period of the local radiation impact beyond this initial phase. Investigating the effects of SSS will thus be an important task of our future work.

5.3 Possible non-linear coupling of AGN feedback and local stellar radiation

In comparison with previous studies, our investigation also incorporates BH growth and consequently AGN feedback. The interaction between AGN feedback and local stellar radiation is complex: an increase in AGN feedback typically leads to a decrease in the SFR, thereby reducing local stellar radiation. Conversely, lower local stellar radiation limits the photoionization of gas, potentially suppressing BH accretion. Thus, local stellar radiation may indirectly influence star formation by regulating BH accretion and subsequent AGN feedback.

To directly assess the impact of local stellar radiation on star formation, we conducted re-simulations of the Au3 galaxy with and without local stellar radiation, while maintaining fixed BH accretion rates at specific redshifts. The BH accretion rate function we prescribed for this experiment is:

We refer to the simulations with self-regulated BH growth as Fiducial simulations, and to the ones with the fixed BH growth as FixBH simulations. Fig. 23 compares the baryonic, the BH, the stellar, and the gas mass evolution from |$z=5.5$| to 0. Inspection of the time evolution of the stellar mass in the lower left panel shows that the stellar mass in the StarRad-FixBH simulation is |$\sim$| 10  per cent lower than in the NoRad-FixBH simulation, implying that the local stellar radiation has a direct effect on the star formation independent of the AGN feedback.

We next consider the evolution of other baryonic components. The upper left panel of Fig. 23 displays the evolution of the baryonic mass and the relative difference of the Fiducial and FixBH simulations of the Au3 galaxy. The large difference of the baryonic mass after |$z\lesssim 1$| disappears. Since the galaxy experiences a merger event at |$z\sim 1.5$|⁠, this result indicates that this huge variation of the baryonic mass within the DM halo is due to the AGN feedback during the merger. The upper right panel shows the evolution of the BH mass. The BH mass in the Fidicial simulations is |$\sim$|2–4 times larger than in the FixBH simulations. However, we can see that there is no significant difference of the baryonic and the stellar mass in the NoRad simulation between the Fiducial and FixBH models. Since higher BH mass indicates a higher AGN output over the evolutionary history, this result indicates that the AGN feedback has only a mild affect on the star formation of this galaxy.

Interestingly, in the StarRad-Fiducial simulation the evolution of the baryonic mass and the gas mass is lower than in the other three simulations after the merger event at |$z\lesssim 1$|⁠. Comparing also the stellar mass in the StarRad-Fiducial and StarRad-FixBH simulations, we can infer that the further suppression of star formation in the StarRad-Fiducial simulation is due to the lower baryonic mass/gas mass compared to the one in the StarRad-FixBH simulation, which is due to the AGN feedback during the merger event. Overall we conclude that the AGN feedback alone hardly affects the star formation in this |$\sim 10^{12}~{\rm {\rm M}_{\odot }}$| halo despite the merger event, which is consistent with previous work about the quenching of the galaxies with the TNG model (Zinger et al. 2020). However, when taking the local stellar radiation into account, the AGN feedback during the merger event becomes amplified since the local stellar radiation can prevent the gas from forming star, making it instead available for fuelling the BH. The AGN feedback will then become stronger, allowing it to reduce the gas fraction in the simulation more effectively, and ultimately suppressing the star formation more strongly.

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~{{\ \rm per\ cent}}$| higher than in NoRad galaxies, whereas the stellar mass is |$\sim 25~{{\ \rm per\ cent}}$| 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 centre 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 H i 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 disc, but it can heat and/or ionize the cold gas in the disc.

• By comparing the rotation curve of five NoRad 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\lt 10^4\, {\rm K}$|⁠) and cool (⁠|$10^4\, {\rm K}\lt T\lt 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 H i and Mg ii column density of CGM gas within |$0.4\, R_{\rm vir}$|⁠, while it does not have a significant impact on the O vi column density in our simulations. This result arises from a combination of heating and photoionization effects due to local stellar radiation. For H i and Mg ii 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, H i and Mg ii experience a more significant reduction in the inner region. While the O vi gas traces the hot gas, however our model of local stellar radiation is too soft to heat and photoionize the hot gas. Thus, the O vi 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) have 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 O vi in CGM. Upton Sanderbeck et al. (2018) also proposed many components of the local stellar radiation like SSS and warm-hot ISM that can possibly affect the properties of the ISM and CGM. These 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 would like to thank the anonymous referee for detailed, insightful and constructive feedback that improved the quality of the manuscript. 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