The physical origins and dominant emission mechanisms of Lyman-alpha halos: results from the TNG50 simulation in comparison to MUSE observations

Extended Lyman-alpha emission is now commonly detected around high redshift galaxies through stacking and even on individual basis. Despite recent observational advances, the physical origin of these Lyman-alpha halos (LAHs), as well as their relationships to galaxies, quasars, circumgalactic gas, and other environmental factors remains unclear. We present results from our new Lyman-alpha full radiative transfer code voroILTIS which runs directly on the unstructured Voronoi tessellation of cosmological hydrodynamical simulations. We make use of the TNG50 simulation and simulate LAHs from redshift $z=2$ to $z=5$, focusing on star-forming galaxies with $8.0<\log_{10}{(M_\star/\rm{M}_\odot)}<10.5$. While TNG50 does not directly follow ionizing radiation, it includes an on-the-fly treatment for active galactic nuclei and ultraviolet background radiation with self-shielding, which are important processes impacting the cooling and ionization of the gas. Based on this model, we present the predictions for the stacked radial surface brightness profiles of Ly$\alpha$ as a function of galaxy mass and redshift. Comparison with data from the MUSE UDF at $z>3$ reveals a promising level of agreement. We measure the correlations of LAH size and central brightness with galaxy properties, finding that at the masses of $8.5 \leq \log_{10} \left(M_\star/\rm{M}_\odot\right) \leq 9.5$, physical LAH sizes roughly double from $z=2$ to $z=5$. Finally, we decompose the profiles into contributions from diffuse emission and scattered photons from star-forming regions. In our simulations, we find rescattered photons from star-forming regions to be the major source in observed LAHs. Unexpectedly, we find that the flattening of LAH profiles at large radii becomes dominated by photons originating from other nearby halos rather than diffuse emission itself.


INTRODUCTION
The Lyman-(Ly ) line of hydrogen at 121.567 nm is one of the brightest emission lines in the Universe. It allows us to detect, and trace the distribution of, galaxies even out to very high redshifts > 5. These Lyman-emitters (LAEs; Partridge & Peebles 1967) can be used to probe the physics of galaxy formation (Finkelstein et al. 2009;Nagamine et al. 2010;Erb et al. 2014) as well as constrain cosmological parameters and large-scale structure (Hill et al. 2008;Adams et al. 2011). Starting in the 80s, spatially extended Ly emission has been detected, often called Lyman-blobs (LABs) and Lyman-nebulae with large extents between ∼ 10 − 100 pkpc (Mc-Carthy et al. 1987;Heckman et al. 1991;Steidel et al. 2000). ★ E-mail: cbyrohl@mpa-garching-mpg.de More recently, Lyman-halos (LAHs) have been discovered around star-forming galaxies that show Ly emission far beyond the galaxies' optical bodies, tracing the circumgalactic rather than interstellar gas (e.g. Hayes et al. 2013). While LAHs are fainter and smaller than LABs in their Ly extent, they might be a generic feature around Ly emitting galaxies. Narrow-band imaging can efficiently detect LAHs at targeted redshifts through stacking (Hayashino et al. 2004;Steidel et al. 2011;Matsuda et al. 2012;Feldmeier et al. 2013), and narrow-band surveys enable ultra-deep, blind samples of LAHs around distant galaxies (Momose et al. 2014(Momose et al. , 2016Kakuma et al. 2019). In the last years modern surveys performed with integral field unit (IFU) spectrographs on 10m-class telescopes, such as the Multi Unit Spectroscopic Explorer (MUSE) and the Keck Cosmic Web Imager (KCWI), take place. These new instruments allow the study of individual, faint LAHs opposed to previous narrow-band stacks. Along with the IFUs' spectral resolution these recent surveys largely increase the information available from LAH observations. Hundreds of individually extended Lyman-halos at 2 have been revealed since (Wisotzki et al. 2016). Many of these are specifically targeted samples which focus on bright quasars (Borisova et al. 2016;O'Sullivan et al. 2020;Cai et al. 2019;Guo et al. 2020), based on strong earlier evidence of enhanced Ly emission around active galactic nuclei (AGN) (Cantalupo et al. 2014;Arrigoni Battaia et al. 2016;Farina et al. 2019;Arrigoni Battaia et al. 2019). Others exploit the large field of view of MUSE, in particular, to conduct blind surveys for LAHs around more typical, generally star-forming galaxies (Leclercq et al. 2017;Wisotzki et al. 2018;Leclercq et al. 2020). At the same time, follow-up with other instruments such as ALMA reveals complementary views on other gas phases within LAHs including CO (Emonts et al. 2019).
Beyond the circumgalactic medium (CGM), first attempts have been made to detect cosmic web filaments in Ly emission (Gallego et al. 2018;Umehata et al. 2019;Lusso et al. 2019), as well as to study the interface between the IGM and CGM as gas flows feed galaxies (Martin et al. 2019). Finally, Ly emission is also a powerful cosmological tool. The Hobby-Eberly Telescope Dark Energy Experiment (HETDEX) is an IFU survey which will detect up to a million LAEs, as well as many extended LAHs, at more moderate spatial and spectral resolution (Hill et al. 2008).
Despite extensive observational detection and characterization, the physical properties and nature of Ly halos remains an open topic. It is unclear whether the extended profiles are mainly sourced by (i) diffuse emission outside of central galaxies, or whether (ii) scatterings of Ly photons emitted from within central galaxies power observed LAHs. Diffuse emission is commonly considered to be sourced by gravitational cooling (Haiman et al. 2000;Fardal et al. 2001;Faucher-Giguère et al. 2010) and fluorescence (Gould & Weinberg 1996;Cantalupo et al. 2005;Kollmeier et al. 2010;Mas-Ribas & Dĳkstra 2016), while star-formation and quasars can provide significant emission within galaxies (Dĳkstra et al. 2006;Zheng et al. 2011) that can scatter with neutral hydrogen in the CGM. Emission from orbiting satellite galaxies can also lead to extended Ly profiles (Mas-Ribas et al. 2017).
The difficulties in determining the powering source of LAHs are closely linked to the resonant nature of Ly photons that can scatter many times in astrophysical environments before escaping towards the observer. This causes the observed frequency and angular position to significantly change due to radiative transfer (RT). Only in the simple, symmetric geometries, RT can be solved analytically (Harrington 1973;Neufeld 1990;Loeb & Rybicki 1999;Lao & Smith 2020). Moving to more realistic setups and in particular hydrodynamical simulations, RT has to be solved numerically.
Recently our theoretical understanding of the Ly emission around galaxies has been pushed forward with the development of cosmological hydrodynamical simulations (Vogelsberger et al. 2020) which are able to produce broadly realistic galaxy populations (Vogelsberger et al. 2014b,a;Genel et al. 2014;Schaye et al. 2015;Dubois et al. 2016;Pillepich et al. 2018b;Davé et al. 2019). Crucially, these simulations predict the full distribution of gas, including neutral hydrogen, in and around galaxies. With volumes of many Mpc a side, these simulations include thousands of well resolved galaxies across the mass spectrum. This enables explicit Ly RT calculations to solve for the propagation and scattering of Ly photons through the interstellar, circumgalactic, and intergalactic media (Laursen et al. 2009;Behrens et al. 2018). Due to computational expense, Ly RT on these volumes can for now exclusively be done in post-processing, even though Ly pressure can be dynamically important (Smith et al. 2019). Simi-larly, a proper treatment of the temperature and ionization state of the gas requires hydrodynamics to be coupled with the radiative transfer of Lyman Continuum (LyC) photons. While these radiationhydrodynamical (RHD) simulations are becoming increasingly feasible, this remains the case only for individual zoom simulations (e.g. Rosdahl & Blaizot 2012;Mitchell et al. 2020) and for cosmological volumes only down to high redshifts around ∼ 5 (Gnedin et al. 2017;Rosdahl et al. 2018;Ocvirk et al. 2020).
There have been various recent efforts to understand the Ly halos in emission by post-processing hydrodynamical simulations (e.g. Lake et al. 2015;Smith et al. 2019;Kimock et al. 2020) with Ly RT. These works typically consider one, to a few tens, of galaxiesrather than full cosmological volumes -making it difficult to draw conclusions about the environmental dependencies of LAHs. Notable exceptions are Gronke & Bird (2017), who predicted LAH properties from the Illustris simulation, and Zheng et al. (2011) who predicted low surface brightness Ly emission from a cosmological reionization simulation ( = 5.7), although with a limited, 'halo-scale' hydrodynamical resolution ∼ 30 pkpc.
In these theoretical works the most important emission origin(s) and source(s) remain disputed. For example, Lake et al. (2015) find good agreement for their set of 9 LAHs with mass 10 11.5 M contrasted with data from Momose et al. (2014), stressing the importance of gravitational cooling in the outer halo to explain the observed profiles, while Gronke & Bird (2017) simulate Lyman-nebulae with masses 10 11.5 − 10 13.5 M and find the simulation can produce Ly extents as large and luminous as those observed, only using central emission from AGN and star-formation.
Recently, cosmological volumes from modern hydrodynamical simulations have been studied in the context of Ly emission focusing on the detectability of the cosmic web (Elias et al. 2020;Witstok et al. 2021). Latest observations of the Ly cosmic web in Bacon et al. (2021) might point at the importance of emission from (faint) galaxies that has been missing in former theoretical explorations.
In this work we improve on several aspects of previous computational studies, revisiting the nature of Lyman-halos. Specifically, we couple the new, high-resolution cosmological magnetohydrodynamical simulation TNG50 (Pillepich et al. 2019;Nelson et al. 2019b) of the IllustrisTNG project to our new radiative transfer code ILTIS. The former provides a competitive combination of volume (a statistically robust sample of ∼ 6, 800 galaxies with ★ > 10 7 M at = 2) and resolution (∼ 100 parcsecs in the ISM, albeit with a simplified sub-grid treatment of the cold phase, and < 1 kpc in the CGM). The latter includes several Ly emission models and a Monte Carlo treatment of the scattering process directly on the full Voronoi tessellation of the gas distribution of the entire simulation volume, enabling a self-consistent treatment of IGM attenuation (Byrohl & Gronke 2020). We note that the TNG simulations are not RHD, and omit on-the-fly radiation transport. Our setup enables us to statistically contrast the simulation predictions to existing LAH observations, while also probing questions regarding the dominant origins, emission sources and relevance of rescattering for the existence of LAHs, and making future predictions in as of yet unobserved regimes.
The structure of this paper is as follows. In Section 2, we describe our radiative transfer code ILTIS, the Ly emission model, and the analysis details of the underlying IllustrisTNG simulations on which the radiative transfer code is run. In Section 3 we present the results for the radial profiles and related reduced quantities from our simulations and a comparison to observational data. In Section 4 we discuss the radial profile shapes and reduced quantities in more detail. We summarize our findings in Section 5. We also point to an extensive Appendix A that highlights different modeling assumptions in the underlying TNG50 simulations and the Ly radiative transfer.

IllustrisTNG and TNG50
We use the outcome of the IllustrisTNG simulations -both the galaxy properties as well as gas distributions -as the basis for our radiative transfer simulations of Lyman-halos. The IllustrisTNG simulations (hereafter, TNG; Pillepich et al. 2018b;Naiman et al. 2018;Nelson et al. 2018a;Marinacci et al. 2018;Springel et al. 2018) are a series of three large-volume magnetohydrodynamical cosmological simulations of galaxy formation. All are run with the AREPO code (Springel 2010), which solves the coupled equations of selfgravity and ideal, continuum magnetohydrodynamics (Pakmor et al. 2011) with a 'moving mesh' discretization based on an unstructured Voronoi tessellation of space.
The TNG galaxy formation model (Weinberger et al. 2017;Pillepich et al. 2018a) is an evolution of the original Illustris model (Vogelsberger et al. 2013;Torrey et al. 2014) and includes a number of important changes. All TNG simulations include models for the physical processes most important for galaxy formation: primordial and metal-line gas cooling, heating from the metagalactic background radiation field, star formation above a density threshold, stellar population evolution and chemical enrichment following supernovae Ia, II, and AGB stars, and the seeding, merging, and growth via accretion of supermassive black holes (SMBHs). With respect to the earlier Illustris model, two critical updates have been made. First, the galactic-scale winds launched by stellar feedback have been revised (Pillepich et al. 2018a), which impacts the gas (and stellar) contents of low mass galaxies in both their ISM and CGM. Further, TNG includes a two-mode blackhole feedback operating in a thermal 'quasar' state at high accretion rates and a kinetic 'wind' state at low accretion rates. The latter is a new model for low-state SMBH feedback, in the form of a time stochastic, directionally variable, high-velocity kinetic wind (Weinberger et al. 2017).
The temperature and ionization state of hydrogen -crucial to Ly emission and scattering -is computed within TNG incorporating primordial cooling following Katz et al. (1996) with additional metal-line cooling from CLOUDY cooling tables. Both metal and primordial cooling are further modified by the assumption of a uniform, time-varying UV background using the intensities given in Faucher-Giguère et al. (2009) for photoionization and photoheating.
Given the substantial impact of AGN, photoionization and photoheating due to their ionizing flux is incorporated through an AGN radiative feedback channel (for details see Vogelsberger et al. 2013;Weinberger et al. 2017Weinberger et al. , 2018. Intrinsic AGN luminosities are assigned strictly proportionally to the black hole accretion rate above a certain accretion threshold. Furthermore, an obscuration factor derived from observational work in Hopkins et al. (2007) sets the escape of ionizing radiation acting upon the surrounding gas under the assumption of optically thin gas. When an AGN with radiative feedback is present in TNG50, the resulting radiation field dominates over the UVB for the studied 50 pkpc around galaxies.
We discuss the AGN radiation field in TNG along with its impact on Ly emission and reprocessing of Ly photons for the radial profiles in Appendix A5.
The ionizing radiation from the UVB and AGN entering the primordial cooling and CLOUDY metal-line cooling is attenuated by a self-shielding factor based on radiation transfer simulations in Rah-mati et al. (2013). These contributions are included on-the-fly, and not in post-processing.
The TNG model has been shown to produce galaxy and circumgalactic medium properties in broad agreement with available observational constraints. The properties of the CGM have been compared in terms of their OVI column densities (Nelson et al. 2018b), gas contents (Pillepich et al. 2018a;Terrazas et al. 2020), x-ray properties (Barnes et al. 2018;Truong et al. 2020), outflows and dynamics (Nelson et al. 2019b), and MgII covering fractions (Nelson et al. 2020), but this is the first study of Lyman-halos in TNG.
Of the three IllustrisTNG simulations, TNG50, TNG100, and TNG300 (Nelson et al. 2019a) we here exclusively use the highest resolution box TNG50 (Nelson et al. 2019b;Pillepich et al. 2019) with a gas mass resolution of bayron = 8.5×10 4 M and a dark matter mass resolution of DM = 4.5 × 10 5 M which is ∼ 15 (∼ 120) times higher than TNG100 (TNG300). The corresponding spatial resolution of TNG50 is of order ∼ 100 physical parsecs in the dense ISM, and this small-scale structure is useful given the strong resolution dependence of Ly radiative transfer at lower resolutions as demonstrated in Behrens et al. (2018), although a more sophisticated model for the cold phase of the ISM would be needed to fully capture Ly RT effects on small scales.

voroILTIS
The Ly radiative transfer is calculated in post-processing on the TNG simulation output with an updated version of ILTIS 1 , a lightweight line emission transfer code as presented in Behrens et al. (2019). ILTIS implements a Monte Carlo approach, spawning singlewavelength photon 'packages' (representing a large number of actual photons) at emission sites, and following their scattering as they traverse the underlying gas distribution. For efficiency, at each scattering event we output the attenuated luminosity contributions along specified line-of-sights towards assumed observers (the 'peeling-off' algorithm; Whitney 2011). We have developed a new version of the code, ILTIS, which runs directly on the unstructured mesh of a Voronoi tessellation. As IllustrisTNG uses this geometry to represent the gas distribution during the simulation and for the evolution of hydrodynamical quantities, no intermediate interpolation or re-sampling steps are required, and the density field in the RT calculation is self-consistent with the simulation. The mesh is re-created in postprocessing with a parallelized wrapper to the Voronoi tessellation code voro++ (Rycroft 2009). We then spawn photons for each cell in the mesh according to the local emissivity, an update which makes simulating diffuse emission feasible.
After emission, photons propagate through the simulation domain, and across the periodic boundaries of the box as appropriate. Given the amplitude of the Hubble flow (several 100 km/s/Mpc) in the simulated redshift range and typical gas velocities, photons are quickly shifted into the far wings of the line profile where they have negligible cross-section. Therefore, we do not need to construct light-cones and simply propagate all photons for the necessary propagation length of order ∼ 10 Mpc. Here, we chose = 28 cMpc/h. An upcoming methods paper will document the code improvements used in the present work (Behrens & Byrohl, in prep). A predecessor of ILTIS has also been used in Behrens et al. (2018) and Byrohl et al. (2019) where additional details are available.

Lyman-alpha radiative processes
Ly photons are dominantly created by recombination of ionized hydrogen atoms with electrons, and de-excitation of excited neutral hydrogen. Different physical processes can power their creation: For recombination, ionizing radiation originates from star-forming regions as well as from the metagalactic ultraviolet background. For collisional de-excitation, the thermal state of the gas provides the gas heating mechanism. We refer to those processes to emission sources and in figures call them 'rec', and 'coll', respectively. As we discuss in the next paragraphs, we use a special description for the recombinations in star-forming regions, which we abbreviate as 'SF' in figures.
We model recombinations after ionization by a spatially uniform background radiation field and AGN radiative feedback as implemented in TNG50. We assume case-B recombination, resulting in an emissivity rec = rec ( ) e HII ( ) Ly (1) where e and HII are the electron and ionized hydrogen number densities, ( ) (from Draine 2011) is the recombination coefficient, rec is the recombination emission probability (from Cantalupo et al. 2008) assuming case-B, and Ly is the Ly transition energy. Note that ( ) and rec ( ) are both temperature dependent. In addition to recombinations, the de-excitation of excited hydrogen atoms can lead to the emission of Ly sourced by collisional excitations depending on the gas thermal state. The rate is proportional to the colliding species, e and HI . The luminosity density is then coll = 1s2p ( ) e HI Ly . (2) We take the collisional excitation coefficient 1s2p ( ) from tabulated fits (Scholz et al. 1990;Scholz & Walters 1991). For star-forming gas cells, where the TNG model invokes a subgrid effective equation of state model for the two-phase ISM (Springel & Hernquist 2003), the simulation's temperature and hydrogen density do not reflect their physical meaning entering Eqn. (1). Hence, we instead derive the recombination rate from the star-formation rate. The ionization of hydrogen by young and massive stars followed by recombination dominates the local Ly emission and the luminosity within galaxies is thus proportional to the star-formation rate.
Here, we adopt a simple linear model for the Ly luminosity density SF = 10 42 where ★ is the star-formation rate within a star-forming Voronoi cell's volume ★ (see Dĳkstra 2019).
This relation is consistent with the H -SFR relation used in Kennicutt (1998) and an assumed Ly to H recombination ratio in the range of 8-10 (Hummer & Storey 1987). It can be directly calculated by integrating the ionizing flux for a given stellar population but has large modeling uncertainties (Furlanetto et al. 2005). Nevertheless this relation is commonly used in simulations and theory as an estimate for the Ly emission from stellar populations (see e.g. Furlanetto et al. 2005;Zheng et al. 2010;Behrens et al. 2018). Equation (3) does not account for Ly destruction by dust in the complex subgrid ISM structure, and we do not take into consideration the stellar populations' age and metallicity. Incorporating these factors would add suppression and scatter in the relation; we discuss the modeling uncertainties and their implications for this paper in Section 4.2.2. In the discussion, comparing to observations, we find that the modeling uncertainties are more important than the lack of dust attenuation for halos in the mass range 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5. Nevertheless, a treatment of dust (see e.g. the effective treatment in Lake et al. 2015;Inoue et al. 2018) and the introduction of sources of scatter in Equation (3) are desirable for future work, and necessary at the high mass end.
The star-formation rate is taken directly from the TNG output in each cell volume. The TNG star-formation model is described in Pillepich et al. (2018a). In short, a gas cell is star-forming if and only if its physical density exceeds 0.1 hydrogen atoms per cubic centimeter, in which case collisionless star particles are stochastically formed (also see Springel & Hernquist 2003;Vogelsberger et al. 2013).
For recombination and de-excitations, we spawn one weighted photon package per Voronoi cell with the luminosity according to the models above. For star-formation, we only spawn photons in Voronoi cells bound to a halo of mass 200,crit > 10 10 M . In TNG50 this results in a sample of ∼ 5000 to ∼ 13000 LAHs at redshifts ∈ [2, 5].
We found that a constant photon count per cell for Equations (1) and (2) results in good convergence for the radial profiles. In such a scheme, photon packages can carry vastly varying luminosity weights. This however is desired as we need to also trace faint, optically relatively thin regions in the outer CGM where surface brightnesses are orders of magnitudes lower than in the central regions of the LAHs.
In contrast, a Monte Carlo sampling relying on a fixed luminosity weight per photon package appears faster for the emission from starforming cells (see Eqn. (3)). For simplicity, we stick with the fixed Monte Carlo photon count per cell, but confirm that this photon count leads to a converged result (see Appendix A2) for star-forming cells.
We refer to these photon packets as 'photons' for brevity. Photons are always injected at the Ly line-center in the rest-frame of the emitting gas cells. We have also explored different spectral emission distribution and generally find little impact -details are given in Appendix A3.
In our model, we do not account for escaping ionizing radiation from star-forming regions that recombines outside of its emission region ("fluorescent radiation"). Fluorescence in our simulations hence only arises from the uniform UV background and the radiation field from AGN incorporated in TNG. In Appendix A5 we discuss the impact of ionizing sources by means of the AGN radiation implemented within the model underlying TNG50. The radiative feedback from AGN is accounted for in that it affects the ionization state, temperature and cooling rates of the gas and hence can boost recombinations and collisional excitations. We do not consider additional recombinations in gas cells due to the activity of SMBHs.

Post-processing and Observational Realism
The radiative transfer simulations provide us with large outputs of Monte Carlo photons. Each has a corresponding luminosity, fre-Figure 1. Ly surface brightness map for the entire TNG50 cosmological simulation at = 3, highlighting the large-scale structure of the cosmic web as seen in Ly emission. We impose a Gaussian point spread function (PSF) with a FWHM of 0.7 arcsec at a binning resolution Δ res of 8.5 ckpc/h and project through a slice depth of 5.25 cMpc/h. The inset panels (no PSF, Δ res = 0.4 pkpc) show two individual Lyman-halos, on the scale of the halo virial radii, for moderate mass objects: M 200 5 · 10 10 M and M 200 1.2 · 10 11 M (top and bottom, respectively). Ly photons are predominantly emitted in the star-forming regions of the central galaxies, resonantly scatter and illuminate the more extended gaseous halos, including filamentary inflows. The more massive halo (lower right) has a number of star-forming satellites which also contribute Ly emissivity and boost the local surface brightness. quency, emission source (see Section 2.3), and positions of initial emission and last scattering. We save two distinct sets of photons: • "intrinsic": Ly photons as directly emitted from gas cells, neglecting any subsequent interactions with gas.
• "processed": Peeling-off contributions after each scattering of propagated intrinsic Ly photons. These photons include the mod-ifications from scatterings, IGM attenuation and potential dust destruction.
Observed Ly light corresponds to the 'processed' photons only. Comparing results for intrinsic and processed photons allows us to quantify the redistribution of photons since their initial emission and the overall importance of Ly RT in our astrophysical setting.
Both intrinsic and processed photons can be filtered based on the originating Voronoi cell, which is recorded by the intrinsic photons, and inherited by any peeling-off contribution.
Using this, we can classify the emission origin of each photon according to four distinct categories: • "central galaxy": The photon originates in the central subhalo (i.e., galaxy) of the targeted halo.
• "outer halo": The photon originates within the targeted halo, but outside of the central galaxy.
• "IGM": The photon originates in a region not associated with any halo.
• "other halo": The photon originates from a region associated with a halo that is different from the targeted.
The definition of those categories relies on the halo and subhalo catalogues provided by IllustrisTNG, where halos and subhalos are defined via the Friends-of-Friends and Subfind algorithm respectively (see Nelson et al. (2019a)). Each photon falls into exactly one of these categories.
First, we compute 2D surface brightness maps for all galaxies centered on the host halo position, using a pixel size of 0.8 pkpc. We include all scattered (i.e. processed) photons irrespective of their origin, as would be seen observationally. As a result these maps include emission from diffuse gas around the halo and even emission from other nearby galaxies and halos. For the projection depth we include all scattered photons reaching the observer from within ±100 pkpc around the galaxy along the line of sight. By adopting this simple prescription we effectively ignore Ly frequency information (diffusion). Observational studies using integral field spectroscopy commonly adopt an adaptive wavelength window to incorporate all Ly flux of the source based on their varying spectral widths. In our approach, where we can isolate frequency diffusion, we can forgo implementing such an adaptive algorithm. In the Appendix A6, we show the impact when incorporating spectral information and adopting a fixed wavelength window around each emitter. Our quantitative results on LAH sizes and qualitative behaviour at large radii remains unchanged irrespective of the chosen method of depth integration.
Unless otherwise stated, we always apply a Gaussian point spread function (PSF) with a FWHM of 0.7 arcsec. This PSF roughly corresponds to that of MUSE UDF data we compare in the later part of the paper (see Section 3.2).

LAH Sample and Reduced Statistics
In our analysis we focus on halos with galaxy stellar masses of 8.0 < log ( ★ /M ) < 10.5. For those halos, one dimensional radial profiles are computed by averaging the pixel values for a given radial bin. We characterize these radial LAH profiles with a number of 'reduced statistics': two measurements of Lyman-halo size, the half-light radius 1/2 and the exponential scale length 0 , and the 'central' surface brightness value SB 0 , which we take as the value of the surface brightness map pixel(s) closest to the halo center, after smoothing by the PSF.
The Ly half-light radius 1/2 is computed from the onedimensional radial profile as the radius enclosing half of the total surface brightness contained within 50 pkpc. Because SB( ) does not vanish at large distances due to contributions from other halos and diffuse gas, this measure depends on this chosen outer radius.
We also fit the one-dimensional radial profile with a single exponential SB(r) = SB 0,fit exp(−r/r 0 ) (as in Cai et al. 2019). The two parameters are the normalization SB 0,fit and the scale length finite lower limit to exclude the impact of the PSF (as is common observationally; Wu et al. 2020) and the upper limit where the profile commonly transitions from an exponential to more flattened shape. Some observational studies fit a sum of two exponentials (Leclercq et al. 2017;Wu et al. 2020), but this method sets one exponential scale length to that of the UV light, adding additional modeling uncertainties which we avoid with the simpler size measure. For comparison with observations, we impose a simple signal-to-noise criterion, considering only data points with / ≥ 5. We derive the noise from a Gaussian standard deviation of 2 × 10 −19 erg/s/cm 2 /arcsec 2 per pixel of the 2D SB maps. While definitions of LAH sizes vary, we impose the same fitting routine to all compared simulated and observed radial profiles for a fair comparison. Due to the variety of radial shapes for individual LAH, the exponential function can be a bad fit at times. We exclude such cases (<10%) by imposing a maximum relative error of 10% for either fit parameter as given from the estimated covariance matrix of the least square fit. We commonly specify the scale length of a given sample in the form of median high low where "low" and "high" are the 16th and 84th percentile.
In this work we refer to comoving units with a preceding 'c', while physical units are explicitly preceded by 'p'. We always present surface brightness in units of erg/s/cm 2 /arcsec 2 .

RESULTS
We begin in Figure 1 with a visual overview of Ly emission from a large cosmological region encompassing thousands of individual emitters. Here we show the large-scale structure of the entire TNG50 simulation at = 3, with scattered Ly photons illuminating not only LAHs, but also the cosmic web within which they reside. The zoom-in panels show two individual LAHs and their substructure (i.e. satellite galaxies), on the scale of the halo virial radius, for moderate halo masses of 200 5×10 10 M (top) and 200 1.2× 10 11 M (bottom), where Ly photons are predominantly emitted at the sites of star-formation in the central galaxy. These photons then resonantly scatter to illuminate extended gaseous halos of the circumgalactic medium (CGM), where the complex dynamics of high-redshift inflows mix with feedback-driven outflows.
In Figure 2 we show surface brightness maps for a collection of nine LAHs ordered by star-formation rate at = 3, including the two LAHs from Figure 1. The three colored contours all trace a surface brightness of 10 −19 erg/s/cm 2 /arcsec 2 , differentiating between the observable flux (i.e. from processed photons; blue), intrinsic emission (red), and intrinsic emission due to star-formation alone (yellow).
Broadly, the Ly scattering process increases the apparent sizes of LAHs beyond that of intrinsic emission. In some instances, such as the lower right panel, red contoured regions occur without corresponding yellow contours, indicating emission without the presence of star-formation in filamentary structures. We note that the surface brightness distributions of the scattered Ly photons are significantly smoother and isotropic than the more complex structure of the underlying density field. Qualitatively this is compatible with the low eccentricities found for LAHs in observations (Wisotzki et al. 2016). At distances beyond 20 pkpc from the central galaxies, LAHs become increasingly anisotropic, due to the combination of satellite galaxies (for the more massive halos) and anisotropic gas inflows. i.e. accounting for scattering and IGM attenuation) contrasted against intrinsic photons (red) and intrinsic photons due only to star-formation (yellow). We commonly find spatially extended intrinsic emission from star-formation in the most massive halos (log 10 ( ℎ / ) ≥ 11.5), while intrinsic emission from other sources is even more extended. The scattering of Ly photons expands the extent of high surface brightness features, particularly in the more massive halos. Here we adopt a pixel size of Δ res = 0.4 pkpc (no PSF).

Radial Profiles
In Figure 3 we show the radial Ly surface brightness profiles extracted from the same nine intensity maps of Figure 2. We decompose the total profile (including rescattering; black) into its contributions from the three emission sources: star-formation (solid blue), collisional de-excitation (orange), and recombination (green). For star-formation, we additionally show the emission neglecting subhalo/satellite contributions, both intrinsic (dotted blue) and processed (i.e. scattered, dashed blue).
In general, radial profiles are steep near the center of the halo and quickly flatten beyond 15 pkpc. The inner regions are dominated by star-formation sourced recombination, and even at larger radii up to 50 pkpc, star-formation often remains the dominant contribution  Figure 2. The black solid lines show the total, observable radial profiles extracted from the surface brightness maps of processed photons. We decompose these profiles into their three physical emission sources, namely star-formation (SF; blue), collisional excitations (coll; orange) and recombinations (rec; green). Generically, the radial profiles drop steeply within 5 pkpc and flatten towards larger radii. For star-formation, we also split profiles into intrinsic (processed) photons from the central galaxy only as dotted (dashed) blue lines, emphasizing the important and non-negligible role of scattering in outer LAHs. while excitations and recombinations (orange and green lines) reach similar magnitude. Except for the most massive halos in the upper row, intrinsic Ly emission from star-forming regions (blue dotted line) quickly fades within 15 pkpc. At larger radii, the profiles are shaped by scattered photons. However, as expected, there are occasional bumps in the profile from star-formation in satellite galaxies of more massive halos.
At larger radii, profiles are shaped by scattered photons, highlighting two important effects of the radiative transfer. First, the central surface brightness is severely damped by photon rescattering (dashed blue vs. dotted blue curves). Second, large amounts of those Ly photons are scattered further out, which provides an important contribution of star-formation to the extended radial profiles even though little in-situ star-formation may take place at those radii.
Collisional excitations and recombinations become important at radii 20 pkpc for low mass halos, and the typical surface brightness contribution from collisional excitation exceeds that from recombination by a factor of ∼ 2. The most massive galaxies more frequently host nearby satellites, which result in the occasional bumps in the profile from star-formation in these subhalos.
We move from the case study of individual profiles to a quantitative exploration of the average predicted LAH profiles. In Figure  4 we stack galaxies based on stellar mass (main panel) and halo mass (lower left panel) at = 3. We also show the evolution with redshift from = 2 to = 5 in a fixed stellar mass bin (lower right panel). Profiles always show the median stacked profile after radial binning, which we note is not the same as first median stacking the two-dimensional surface brightness maps. For the radial binning we calculate the mean surface brightness for a given annulus. Shaded regions show the central 68 percentiles. The dashed line in the main panel indicates the 0.7 arcsec FWHM Gaussian PSF we adopt, which dominates the smoothing of the radial profiles at small distances < 10 pkpc. At larger radii, the surface brightness rapidly flattens, as we explore below.
In the top panel, stacking surface brightness in stellar mass bins from low-mass galaxies with ★ = 10 8 M to Milky Way mass systems with ★ = 10 10.5 M , we see that Ly surface brightness increases monotonically, at all radii, with increasing stellar mass. Despite this strong correlation between peak surface brightness and stellar mass, the overall shape of the median radial profiles is largely independent of the stellar mass. We note that a mass-dependent Ly photon escape probability, e.g. due to dust physics, would impact the trend of overall luminosity and stellar mass, which we explore further in Section 3.5.
The lower left panel of Figure 4 shows the radial profiles in three halo mass bins. The central surface brightness of the Ly profiles  Overview of the predicted median Ly surface brightness profiles around TNG50 galaxies, as a function of mass and redshift. Main panel: Stacked radial profiles for five different stellar mass ranges from ★ = 10 8 M to 10 10.5 M , at = 3. The respective count of contributing halos to the stack is given in parentheses in the legend. At fixed cosmic time, the central Ly surface brightness increases monotonically with stellar mass, roughly as SB(r = 0) ∝ M 0.9 ★ . We overplot the median radial profile of 58 observed LAHs at redshifts between 2.9 and 3.5 presented in Leclercq et al. (2017) from the MUSE UDF (see text for details). Lower left: Stacked radial profiles as a function of dark matter halo mass, also at = 3. The solid lines show the observed radial profiles, while the dashed lines only consider those photon contributions originating from the central halo. Lower right: Redshift evolution of stacked radial profiles for a fixed galaxy stellar mass bin of 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5. At fixed stellar mass, Ly halos are more luminous towards lower redshift. We contrast the full radiative transfer result (RT; solid lines) with the intrinsic emission profiles (dashed lines). The scattering which occurs during the RT lowers the Ly surface brightness at halo center ( 5 pkpc) while increasing it at larger radii, producing an overall flatter profile.
rises as a function of halo mass. In addition, we more clearly find a change in shape of the radial profiles as a function of halo mass, whereby flattening begins at smaller radii for lower mass halos. Considering central galaxy emission only (dashed lines), we observe that as star-formation decreases towards lower mass, the luminosity budget available for rescattering in the halo does likewise. Hence, diffuse emission outside of the central galaxy becomes dominant at smaller radii. Equivalently, since the central (or total) luminosities are lower at lower mass, external emission (i.e. from other halos) can take over more quickly.
In the lower right panel of Figure 4 we explore the redshift evolution of radial profiles in a fixed stellar mass bin of 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5. Towards higher redshift, the central surface brightness drops significantly while changes to the overall shape are minor. Here we also show the radial profiles of the intrinsic photons (dashed lines), where Ly photons are allowed to escape directly to the observer without scattering, contrasting against the processed emission (solid lines). The intrinsic profiles' lower brightnesses at higher redshifts are driven by surface brightness dimming, which is however largely countered by the higher specific star formation rates at fixed stellar mass. The intrinsic central SB luminosities decrease much more slowly towards higher redshift, implying that the photon redistribution due to resonant scattering is significantly more important at higher redshifts.

Comparison of TNG50 and MUSE data
In the top panel of Figure 4 we also overplot the median radial profile from the observational LAH dataset of Leclercq et al. (2017) for objects with redshifts between 2.9 and 3.5. The data set is based on the MUSE Ultra Deep Field (UDF), which finds extended Ly emission around 145 of 184 star-forming galaxies at 3 6 with a median of ∼ 3.7. The observed galaxies extend in stellar mass down to ∼ 10 7 M with an average stellar mass of ★ ∼ 10 8.5 M (Boogaard et al. 2018). As a result, the most appropriate comparison is against the orange line, where we find a excellent agreement of the normalization and radial shape of the surface brightness profile between MUSE and TNG50. At larger radii > 20 pkpc the observed stacked profile becomes uncertain given the large errorbars, so we cannot ascertain whether or not the strong flattening we observe in the TNG50 LAHs is also seen in the data. Additionally, this flattening is affected by MUSE's more extended Moffat PSF (Bacon et al. 2017). The agreement of stacked radial profiles between observed and simulated samples degrades toward higher redshifts, particularly = 5.0. This discrepancy is largely driven by the flattening of the simulated profiles at high redshifts (see bottom right panel) that is not found in the observed sample.
We note that quantitatively comparing TNG50 to Lyman-halo observations has many challenges and subtleties. With respect to the MUSE data of Leclercq et al. (2017) in particular, although we address the issue of spatial smoothing and the PSF, remaining systematics could include details of (i) spectral smoothing, resolution, and the construction of the surface brightness map with a frequencyspace integration; (ii) noise, including relevant surface brightness limits and sky backgrounds; (iii) the position chosen as the center of the halo; (iv) the derivation of galaxy properties, including stellar mass, where SED fitting based on HST broad-band photometry lacks near-IR points, and neglecting emission lines including Ly which can fall into the F606W filter (Feltre et al. 2020); and (v) sample selection and selection biases, i.e. choosing appropriate analogs for comparison to the observed halos and/or matching the observed galaxy population in general.
Beyond the stacked profile comparison, we also contrast individual LAHs as observed in the MUSE UDF to those from TNG50. Figure  5 shows the 24 radial profiles from Leclercq et al. (2017) closest in redshift to = 3 (blue data points with errorbars). For each, we search for the best 'match' from among our catalog of simulated LAHs at that redshift, and select the single LAH with the minimum least-squared difference 2 which are overplotted (orange).
In general, we are able to find excellent matches to the observed data, demonstrating that TNG50 can reproduce the diversity and variety of observed LAH profile shapes. For instance, we show good 2 Note that the low surface brightness measurements at large radii > 2 ", have only small weights in the linear least-square fit. matches for more compact and more extended objects (e.g. id 2178, id 6521), and although noise starts dominate at larger radii, we also find good matches for candidates with very flattened profiles (e.g. id 1226). The stellar masses of the observed galaxies, and the stellar mass of the selected TNG50 analog, are both shown in the legend for reference. For the subsample of observed LAHs for which we have a stellar mass estimate (from Feltre et al. 2020), we find a mean difference in ★,MUSE / ★,TNG of 0.11 dex with a standard deviation of 0.43 dex. This indicates that the simulated LAHs selected as good matches surround galaxies with comparable stellar masses as the observed systems. Figure 6 shows a similar comparison of matched LAH profiles between the MUSE data and TNG50. However, we now restrict the least-square fit to those values at distances below ≤ 2 arcsec, and include the five matches in each case to highlight the range of predicted large radii behaviour from simulated profiles. Most of those fits are compatible with the observed radial profiles and their errorbars at large radii (see also Song et al. 2020, for fits to L17 SB(r) profiles using 3D Ly RT coupled to idealized shell models rather than cosmological simulations). A notable exception is the profile of MUSE-id 1281 which has an excess at ∼ 2 arcsec, possibly due to the existence of a satellite galaxy.

The Origin and Source of Lyman-alpha Halo Photons
Although LAHs are observed localized around galaxies and their dark matter halos, the photons which contribute to that emission can arise from a number of different origins. In Figure 7 we show the relative contribution to a stacked Ly surface brightness profile, depending on the origin of emission, for galaxies in the stellar mass bin 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5. We categorize the emission origin of each photon as one of the four categories previously introduced in Section 2.4.
We find that emission from the central galaxy (blue) clearly dominates the radial profile below 20 pkpc after which emission originating from other halos but scattering onto the targeted halo becomes increasingly important. Above 40 pkpc this 'other halo' origin (red) even dominates the radial profiles. The radial profile from rescattered photons originating in other halos has a very shallow slope, thus leading in large part to the flattening of the overall profiles at larger radii. In particular, we find that the contribution from other halos is significantly boosted if more massive halos are nearby, an effect we explore more in Section 4.1.2.
Emission originating in the IGM (green) and outside of the central subhalo (particularly from satellites; orange) is generally negligible, and never contributes more than a few percent to the total stacked profile. However, there can be infrequent radial profiles of individual halos with larger contributions from IGM and satellites than for the shown average. We note that the IGM contribution in particular will depend upon the line-of-sight integration depth.
For higher mass halos (not shown) we find that the emission from the central subhalo grows more rapidly than any of the other origins, pushing the observed flattening to larger radii. Analogously, lower mass halos flatten at smaller radii. Other halos start to significantly contribute (≥ 10% to the total stacked profile) at 7, 22 and ≥ 50 pkpc 3 for the respective 1 dex stellar mass binned halos starting at 7.5, 8.5 and 9.5 log ( ★ /M ). There is very little redshift evolution for these radii from = 2 to = 5.  = 3. Every observed LAH is matched to a simulated halo from TNG50 by choosing the best least-squares fit profile. All simulated radial profiles are smoothed with a Gaussian 0.7 arcsec FWHM PSF, and overplotted (orange lines). This demonstrates that the simulation has the diversity and sample statistics to recover excellent matches to all observed halos, and that TNG50 can successfully reproduce every observed profile, in both normalization and shape, with at least some simulated halo. For each halo, we include the ID and redshift of the observed MUSE object, as well as its stellar mass estimate if available (Feltre et al. 2020). We also show the exponential scale length 0 (in pkpc) fitted to the MUSE data. For the simulated profile we also include its galaxy stellar mass: observed LAHs and their simulated matches from TNG50 surround galaxies of similar mass.  [8.5, 9.0] TNG L17 Figure 6. Gallery of observed Lyman-halos from the MUSE UDF (blue points with uncertainties; Leclercq et al. 2017), where here we show the first 8 LAHs shown from Figure 5. As before, we search for the best matching LAH from TNG50, but now only fit data points for small radii < 2 arcsec, and show the best five fits for each profile. As before, we contrast the MUSE stellar mass (if available) with the stellar mass range of the five simulated profiles, which are shaded from light to dark orange with increasing stellar mass. We find that although the simulated profiles show object-to-object variances at large distance, in most cases they are consistent with the observed profile within its errorbars. Median stacked radial Ly profile for galaxies with stellar masses 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 at = 3 in TNG50. We decompose this profile into the photons with differing emission origins: from the central galaxy (blue), outer halo (orange), the IGM (green), and other halos (red). At < 20 pkpc emission from the central subhalo dominates, beyond which contributions from other halos start contributing to the overall shape. Beyond 40 pkpc the 'other halo' origin is critical and produces the flattening of the profiles towards large radii. The contributions from emission originating in the outer parts of the halo and the IGM are negligible.
In Figure 8 we similarly decompose the stacked profile into the relative contributions of different emission sources: star-formation sourced rescattered photons (blue), collisional excitation (orange), and recombination (green). In addition to the 'processed' signal (solid lines) we also show the intrinsic emission signal (i.e. ignoring scattering effects; dashed lines). In both cases we find that star-formation makes up the bulk of the SB within the central 10 pkpc. At larger radii however rescattered photons from star-forming regions drop to a ∼ 50% relative contribution as diffuse collisional excitations and recombinations rise to ∼ 30% and ∼ 20% respectively. These relative fractions, shown here at = 3, are similar at other redshifts (not shown).

Lyman-alpha Halo Sizes
To study the dependence of LAH profiles on galaxy/halo properties and redshift, we summarize each profile by a characteristic size and surface brightness, as defined in Section 2.5. Figure 9 shows the Ly half-light radius 1/2 as a function of stellar mass across the studied redshift range. We find that the half-light radii of our LAHs is typically between 5 and 15 pkpc. There is little correlation with stellar mass compared with the scatter. Half-light radii are systematically larger in physical kpc towards higher redshift, roughly doubling in size from = 2 to = 5 across the redshift range. This implies a stronger redistribution through scatterings into the outskirts of halos at higher redshifts, which is largely driven by a higher neutral hydrogen density at those redshifts.
In Figure 10 we show the fitted single exponential scale length 0 . In all panels we impose the noise modeling described in Section 2.5, which effectively imposes a sensitivity limit and induces a lower cutoff mass below which modelled LAHs show no observable extent. This limit starts to affect the median at a stellar mass below ∼ 2 × 10 8 M (∼ 10 9 M ) at = 2 ( = 5), causing the artificial drop-offs at low ★ .
Except for this sensitivity limit, there are no clear trends of 0 with stellar mass. Similarly, no clear redshift evolution is evident (upper panel). However, when looking at our fiducial stellar mass range of 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 only, we do find that 0 increases with redshift from 3.4 +1.3 −0.7 pkpc at = 2 to 3.7 +1.7 −0.9 pkpc (4.0 +1.7 −1.0 pkpc, 5.5 +2.4 −1.6 pkpc) at = 3 ( = 4, = 5). The middle panel of Figure 10 shows the scale radius as a function of stellar mass at = 3. Here we decompose the contribution to LAH size by emission source, by determining 0 based on each of three respective emission sources alone (colored lines). The black solid line shows the median 0 relation combining all three emission sources, as would be observable. We find that the overall scale radius 0 is largely determined by the emission from star-formation due to its high surface brightness. Scale lengths from star-formation show no correlation with stellar mass, while both collisions and recombinations do show a strong positive correlation. Thus, 0 follows the lack of evolution with mass seen in the star-formation source. In contrast, both diffuse collisions and recombinations show a strong positive correlation with mass. For halos between 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 we find 0 to be 3.7 +1.7 −0.9 pkpc for the overall profile and 3.1 +1.4 −0.7 pkpc, 8.0 +2.8 −1.7 pkpc or 9.2 +3.2 −2.0 pkpc for SFR, excitations, recombinations only respectively.
Solid lines indicate processed photons, while the dashed lines show 0 based on the intrinsic photons only. Not surprisingly, intrinsic photons typically give rise to a smaller scale length than processed photons. Most importantly, intrinsic photons from star-formation give rise to a scale length close to a point-like source convolved with the PSF, while this value doubles from rescattering in the CGM.
In the bottom panel of Figure 10 we show the explicit redshift evolution of LAH sizes, from = 2 to = 5, and compare the TNG50 result with that of Leclercq et al. (2017). For the comparison with MUSE, we show simulated halos in stellar mass bins log 10 ( ★ /M ) from 8.0 to 9.5, which overlap the observed distribution. At fixed stellar mass, we find that the median 0 increases towards higher redshift and the magnitude of this trend is similar to the scatter of 0 at each redshift. For the MUSE UDF observations, we find 3.1 +1.0 −0.6 pkpc, 2.9 +1.3 −0.5 pkpc and 2.4 +1.2 −0.3 pkpc for the redshift bins centered at = 3, = 4 and = 5, respectively.
Our face value comparison with the MUSE data implies that LAHs in TNG50 are slightly more spatially extended than in observations, by roughly a factor of ∼ 1.2 at redshift three, increasing to ∼ 2.3 by redshift five with a simulated 0 ∼ 5.5 pkpc, the two size distributions being roughly one standard deviation apart. We caution, however, that the stellar mass distribution of the observed galaxies is not fixed as a function of redshift. More massive galaxies are more easily observed at higher redshift, and ideally we would match the joint ( ★ , ) distribution to make this comparison. However, remaining methodological differences likely still dominate the uncertainty in this comparison, as sizes are not measured in exactly the same way in both the observational data and TNG50 simulation. We note that observed sources of Ly emission have been characterized by many different sizes ranging from roughly 2 pkpc to 9 pkpc at ∼ 3 (Bond et al. 2009;Momose et al. 2014). This large range of radii hints at the large diversity of different galaxy selection (LBGs/LAEs), functional fitting forms (single/double exponential) and methodologies (individual/stacking). We leave a quantitative comparison of Ly sizes for future work, and do not try to explicitly compare LAE sizes to the observational literature here.

LAH Central Brightness
In addition to the two LAH size measures 1/2 and 0 we also calculate the central surface brightness SB 0 . Figure 11 shows the trend between SB 0 and galaxy stellar mass, from = 2 to = 5 (different line colors). We derive this value for both the intrinsic (dashed lines) and processed photons (solid lines). For the latter, we shade the central 68% of outcomes in a given mass bin. Additionally, we show the results from the MUSE UDF data set (circles), colored to match the nearest simulated redshift for comparison.
First, we see a strong correlation of the peak LAH brightness with galaxy stellar mass. This trend was previously noted in Figure 4 for the simulated halos. Over four orders of magnitude in stellar mass, the central surface brightness value increases by roughly two orders of magnitude. The evolution of SB 0 as a function of redshift is minimal for intrinsic photons. While the intrinsic surface brightnesses are subject to cosmological dimming, this seems to be countered by the increased specific star formation rate for halos at higher redshifts. For Decomposition of 0 at = 3.0 (black) into its three emission origins (SF, coll, rec), for both intrinsic and processed photons. As the luminosity budget in the proximity of the halos' center is dominated by star-formation, the latter effectively sets the scale length. Bottom: 0 versus redshift compared to observational data from the MUSE UDF (Leclercq et al. 2017). For simulations and observations, we fit the scale length 0 using the same procedure. We show three different bins of fixed stellar mass for TNG50 (blue lines). At fixed stellar mass, LAH sizes are overall larger towards higher redshift. No clear redshift trend is evident in the observations which are consistent with no size evolution, although the galaxy stellar mass distribution as a function of redshift in the data is uncontrolled (see text for details). There also is significant redshift evolution of SB 0 for the processed photons, roughly scaling with (1 + ) 4 . In the MUSE data, we find no significant correlation, neither with stellar mass nor redshift.
processed photons, CGM diffusion and IGM attenuation suppress SB 0 with a scaling of roughly (1 + ) 4 . In comparison to the clear correlation in the simulations, the MUSE data does not show such a strong relationship between central surface brightness and either stellar mass or redshift. Although the SB 0 values are in reasonable agreement where the bulk of the observed systems reside, 10 8 < ★ /M < 10 9 , the flat trend in the data leads to lower inferred values at higher stellar masses, when compared to those in TNG50. This is certainly caused in part by systematic uncertainties in the observational determinations of stellar mass, which are acknowledged for the observational estimates in Feltre et al. (2020). These uncertainties in part wash out the trend with stellar mass found here. Even more importantly, we speculate that this difference arises due to our omission of a model for unresolved dust attenuation and stochasticity on the smallest scales, as we discuss below.

The Shape and Nature of Lyman-alpha Halos
Based on our results we now discuss the implications for the shape and nature of Lyman-halos. In particular, is there a typical shape for LAH radial profiles, is there a common cause for this shape, and what is the resulting interpretation of the observations? In Section 4.1.1, we focus on 'small' scales, of order of ∼ 10 pkpc, where LAHs are detected around star-forming galaxies. In Section 4.1.2 we discuss larger scales, and profile flattening, as accessible in current and future intensity mapping studies.

Ly Profiles at Small Scales
We have shown that the radial profiles of our simulated LAHs are primarily shaped through rescattered photons which originate from star-forming regions in the central galaxy of a halo on scales around ∼ 10 pkpc. This rescattering gives rise to a smoothing of the surface brightness maps and radial profiles that is larger than typical PSFs (0.7 arcsec adopted herein). Beyond this distance median radial profiles tend to steepen rapidly before flattening on scales above ∼ 20 pkpc.
We find typical exponential scale lengths 0 of ∼ 4 pkpc with little to no correlation with stellar mass. As radial profiles are largely dominated by emission from star-formation on these scales, the typical shape of rescattered photons from the central galaxy sets this typical 0 and leads to the lack of correlation with stellar mass. If diffuse emission through collisions and recombinations were the dominant source of LAH photons, we would infer much larger scale radii 0 and a strong correlation with stellar mass with scale lengths starting at ∼ 5 pkpc growing to ∼ 15 pkpc between 10 8 and 10 10 M . Thus, 0 and its mass dependency can serve as a discriminator between rescattered photons from star-forming regions and diffuse emission in observations. This holds even if the relative importance of the different emission sources is not precisely correct in our simulations.
An important modification of the Ly emission stems from local ionizing sources. TNG's simplified on-the-fly treatment of radiation stemming from AGN is already reflected in our results. In particular, in our fiducial sample of LAH candidates with stellar masses of 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 1030 out of 1189 halos host a SMBH and incorporate their ionizing flux. While we expect that additional ionizing flux from stellar populations could additionally boost diffuse emission, we only find a small impact of local ionizing radiation from AGN on the radial profiles on small scales measured by the exponential scale radius 0 . As the AGN UV radiation (when present) dominates over that of the stellar population at lower redshifts considered, we similarly do not expect our 0 results to significantly change if stellar radiation was incorporated.
Interestingly, as the scale length 0 keeps growing toward higher masses for emission sourced by cooling and the UV background, high mass halos in TNG50 reach the lower end of extents observed in LABs. While a future dedicated study is required, this might strengthen the case of such diffuse emission sources (without local sources ionizing flux) causing observed LABs.
In comparison to data, our results for the stacked profiles are consistent with MUSE UDF observations presented in Leclercq et al. (2017). At the level of individual LAH profiles, we similarly find good agreement, such that there are numerous analogs in TNG50 which have compatible Ly radial profiles. This agreement is noteworthy as we have no tuned or calibrated parameters in our Ly modeling.
In the quantitative comparison of 0 for individual LAHs we find up to a 20% mismatch at = 3, which grows to a factor of two at = 5. Similarly, the central surface brightness values from Leclercq et al. (2017) show significant scatter and tend to be below those obtained from TNG50 at the highest stellar masses. Despite these regimes of tension, the observations show no clear mass or redshift evolution in either 0 or SB 0 values, which together with the relatively small sizes implies that observed LAHs in the MUSE UDF sample are sourced by star-formation.
Although a rigorous comparison between simulated and observed LAHs is complicated by several systematic uncertainties, the SB 0 and 0 tensions hint at possible shortcomings of our Ly modeling. In Figure 10 we found that considering star-formation alone gives the smallest exponential scale lengths. At = 3, we found 0 = 3.1 +1.4 −0.7 pkpc in this case, which is very similar to the MUSE UDF estimate of 3.1 +1.0 −0.6 pkpc. The simplest explanation is that we have underestimated the Ly luminosity from star-forming regions, or overestimated the Ly luminosity from other sources. Further, the strong correlation of SB 0 with stellar mass found in our models is not clearly present in the data. As we assume a fixed relation between starformation and Ly luminosity this outcome is not surprising, and can be alleviated by developing a more realistic model for the underlying relation between Ly and SFR, as discussed in Section 4.2.1.

Ly Profiles at Large Scales
Although much of our analysis has focused on the inner regions of LAH profiles, 20 pkpc, we here consider the general shape of extended Ly profiles where the impact of environment reveals itself.
The most prominent feature in the median stacked profiles of Figure 4, regardless of stellar mass, is the significant 'flattening' beyond the inner few 10s of pkpc. In Figure 12 we take advantage of our global, large-volume RT calculation and present the radial profiles of halos with 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 out to 3000 pkpc. Here we neglect the PSF, and integrate through ±600 pkpc from the emitter's position along the line of sight, which roughly reflects the spectral resolution of the HETDEX survey.
We decompose this median surface brightness profile based on photon origin (different colored lines). For reference, the expected surface brightness given the average luminosity density in TNG50 is shown with the black dotted line. We find that the contribution from other halos (solid red line) is nearly constant with distance and reaches down to roughly match the mean global. This 'other halo' contribution starts to dominate the surface brightness profile beyond 100 pkpc by a factor of more than 10. To emphasize this environmental effect, we also show a stack of the subsample of halos that have a close (≤ 0.5 pMpc) massive neighboring halo, by a factor of ten or more, as the red dashed line. For those galaxies there is an elevated plateau above the global mean that only starts to drop around 100 pkpc.
This behavior is a proximity effect whereby nearby halos give rise to an effectively elevated background. As a result, emission from star-formation from nearby neighbors dominates the observed radial profiles over diffuse emission from the outer halo and the IGM. Zheng et al. (2011) find a similar flattening effect from rescattered photons from other halos in their simulations, but are unable to assess the relative importance of diffuse emission, which remains untreated. In contrast, Lake et al. (2015) do incorporate both central sources and diffuse emission, reporting a flattening up to scales of 80 pkpc that in equal parts stems from the two emission sources, in contrast to our results where emission from other halos dominates over IGM emission.
We note that modeling of this proximity effect remains tricky: off-scattered photons in a targeted halo originating in a neighboring halo could potentially not be modeled by a classic 2-halo term as this would merely capture the overlap between different halos' profiles.
Here -in addition to the possible overlap of profiles -we have Ly photons freely traveling through the IGM to a neighboring halo at which point scatterings will trace out part of the targeted halo.
Our findings are qualitatively compatible with stacking results by Matsuda et al. (2012) showing a strong correlation of the flattening in the overdense regions. Wisotzki et al. (2018) find significant flattening up to scales of around ∼ 50 pkpc, for which they consider UVB fluorescence as a potential source. In contrast, we find that towards large radii diffuse emission outside of halos contributes less than 10% to the surface brightness (see Figure 12), while we find that scattering of Ly photons from galaxies dominate the flattening. More recent Ly intensity mapping results at = 5.7 and 6.6 by Kakuma et al. (2019) found reasonable agreement with the MUSE stacking (see also Matthee et al. 2020) and the simulated profiles of Zheng et al. (2011) on small scales ( < 150 ckpc). However, they cannot confirm the proximity effect and flattening found here for larger scales. We note a further complication, that this flattening could be removed in large part or entirely in narrow-band surveys due to the required background subtraction. However, a careful examination of upcoming data, such as with HETDEX, could reveal this proximity effect, particularly by stacking based on the presence of massive neighbor number or environmental density.

Model Limitations: TNG and its gas state
Our results depend inseparably on the outcome of the underlying cosmological simulation as well as our Ly radiative transfer postprocessing, inheriting limitations from each. The basis of our Ly calculations is the TNG galaxy formation model (Weinberger et al. 2017;Pillepich et al. 2018a) and the TNG50 simulation in particular. Although the TNG model has been shown to reproduce a large variety of galaxy and CGM properties in broad agreement with observations (see Section 2.1), it does not treat local ionizing radiation self-consistently as done in radiation-hydrodynamical (RHD) simulations, that are now becoming increasingly feasible for individual zoom simulations (e.g. Rosdahl & Blaizot 2012;Mitchell et al. 2020) as well as cosmological volumes down to redshifts = 5 − 6 (Gnedin et al. 2017;Rosdahl et al. 2018;Ocvirk et al. 2020). As a result, our knowledge of the gas ionization and temperature state is limited to the outcome of the physical modeling in our simulation, as is true in all galaxy formation simulations. As TNG50 remains unfeasible as RHD simulation, it substitutes such approach with an on-the-fly treat-ment of AGN as the dominant ionizing radiation source as studied in Appendix A5. Towards higher redshifts, the relative importance of stellar ionizing sources increases, thus potentially hinting at larger uncertainties of TNG's model only incorporating AGN's radiative feedback.
While particularly emission from collisional excitation strongly depends on the temperature (Furlanetto et al. 2005;Faucher-Giguère et al. 2010) and subsequently can boost emission, we found that the diffuse emission remains subdominant to scattering from the central sources even if radiation from AGNs is present. If central Ly emission from SMBHs themselves is incorporated, this finding will be further strengthened. Predicted central surface brightnesses when only considering collisional excitations from the CGM in TNG, even with AGN radiative feedback, are commonly an order of magnitude too small compared to the MUSE observations. Apart from this normalization problem, collisional excitations from TNG's predicted temperature and ionization state show scale lengths significantly too extended for compatibility with observations (see Figure 10). Reconciliation with observations in case collisional excitations are the dominant source of LAHs thus requires significantly different ionization and temperature radial profiles from those found in TNG. Independent constraints on the density, ionization state and temperature from other observational probes would be useful to assess TNG's model and check resulting conclusions for predicted LAHs. In addition, these findings should be revisited with future RHD simulations that come closer TNG's sample size and studied redshift range.
In addition, the TNG50 simulation has a resolution of ∼ 100 pc in the dense interstellar medium, which is furthermore modeled with an effective two-phase model. As above, this implies that a number of approximations and simplifications must be adopted in both the emission and transport of Ly (see Section 2.3). Importantly, in the present work, we have adopted a simple mapping from the starformation rate of a gas cell to its Ly emissivity. We have also neglected dust in its entirety, as well as sub-resolution (unresolved) density structure during the radiative transfer step.

Model Limitations: ISM emissivity
As discussed, we do not include the destruction of Ly photons by dust in our Ly radiative transfer (although we have investigated its impact, see Appendix A4). Such destruction would primarily take place in the ISM, where resolution is marginal for capturing the complex density and ionization structure relevant for the Ly radiative transfer. If however we do not explicitly treat dust, we potentially overestimate the SF luminosities in Eqn. (3), where we do not explicitly model dust destruction. We review the emission model for stellar populations here, as such an overestimation would strongly affect our conclusions of central Ly emission dominating observed LAH profiles through CGM scatterings.
Equation (3) is derived by integrating the ionizing flux from stellar population synthesis models, and a conversion of ∼ 2/3 of ionizing photons into Ly (Dĳkstra 2019). The proportionality factor between SFR and Ly emission depends on a range of factors such as the stellar population synthesis model, the escape of ionizing flux, the stellar population's age and metallicity, and the initial mass function (IMF) (Furlanetto et al. 2005). Overall, this relation should only be seen as a rough estimate with an uncertainty of a factor of a few. Furthermore, observations, such as from H , Ly / and Ly /H show a large scatter in the relation of hydrogen line emission to SFR (Kennicutt 1998;Blanc et al. 2011;Weiss et al. 2021).
Given the large modeling uncertainties shown, we adopted this SFR-Ly relation as a rough estimate, which could potentially be rescaled after comparison with observations. However, as we show in the following, such a calibration is not necessary for the mass range 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5, where we already obtain reasonable agreement with observations. This also implies that, while not used in the derivation of the relation, Ly dust extinction is nevertheless effectively captured for the 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5. We note however, that at the high mass end luminosities can be significantly overestimated, and the role of dust becomes critical.
In Figure 11, we show the central luminosities as predicted in our simulations compared to MUSE observations in Leclercq et al. (2017). We find that the simulations show similar or slightly lower central surface brightnesses when compared with the observational data.
We have also compared the LAH luminosity function (LF) of our samples with observational LAE luminosity functions at = 2.2, 3.1 and 5.7 (respectively based on Konno et al. 2016;Ouchi et al. 2008;Konno et al. 2018) (not shown). We calculate LAH luminosities from the photon contributions of the targeted halos falling within a 3 arcsecond aperture. We find very reasonable agreement of the observations with our intrinsic luminosity functions in the luminosity range ∼ 3 × 10 41 erg/s to ∼ 2 × 10 42 erg/s, which includes at least 68% of all halos with stellar masses 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 for each redshift = 2 to = 5 in TNG50. Quantitatively comparing to Konno et al. (2016), we find that at = 1.4 × 10 42 erg/s, our intrinsic emission in fact underestimates the LF by just 0.16 (0.32, 0.17, 0.06) dex for = 2 (3,4,5). For the processed photons, our LFs underestimate the observed LFs by 0.19 (0.10, -0.11, -0.63) dex at = 2 (3,4,5). That is, our preliminary look already shows promising agreement, while a detailed analysis of the LFs remains a topic for future work.
We are thus confident that, on average, for the fiducial stellar mass range of 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5, the luminosities reasonably match observations. We thus adopt the simple SFR-Ly relation without readjustment for dust attenuation. However, we overestimate the high luminosity end of the LF, and thus inaccurately capture the high-mass LAHs' shape. This might also affect the overall amplitude of the flattening from neighboring galaxies for the fiducial stellar mass range.
In observations there is a large scatter between the star-formation rate inferred from Ly luminosity compared to UV based estimates (Santos et al. 2020;Runnholm et al. 2020). Typically, the median SFR estimate from Ly exceeds that from UV measurements below ∼ 10 M /yr, but decreases above this value. This systematic trend implies a larger suppression of Ly escape for emitters with higher SFR and stellar mass. In addition, there is non-negligible scatter between individual objects with similar properties. It is thought that Ly emission traces the most recent star-formation (Santos et al. 2020), which could be highly time variable, and is modified by the complex small-scale neutral hydrogen distribution and kinematics (Blanc et al. 2011).
In our simulations we have instead assumed a strict proportionality between Ly luminosity and SFR at the resolution scale, as given by Equation (3). We find that subsequent scattering only adds minor scatter to this correlation, implying that sub-resolution stochasticity may be required. As discussed in Section 4.1.1, we tend to slightly overestimate the 0 metric of LAH sizes in comparison to the MUSE UDF sample, and generally find strong trends of both 0 and SB 0 with stellar mass, which are less clear in the data. We suspect that adopting an observationally motivated scatter between Ly luminosity and UV inferred SFR would also alleviate the tensions noted in the central surface brightness SB 0 trends between our simulations and the data. We ran a simple test to test the impact of including a naive dust treatment (Appendix A4), which demonstrates that the central brightness SB 0 is increasingly attenuated towards higher stellar masses, as dust counters the generally higher SFRs and thus Ly luminosities of those halos.
While subgrid modeling and dust, along with intrinsic scatter between Ly luminosity and SFR, should help reconcile the central surface brightness SB 0 comparison in Figure 11, we would additionally require the mean Ly luminosity from star-formation to increase in order to obtain smaller LAH sizes. For the typical objects considered as LAH candidates with 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 here, Santos et al. (2020) find SFR Ly /SFR UV ∼ 2 for SFR < 10 M /yr with significant scatter. Beside the simplistic model assumptions with significant uncertainties as previously discussed, we would need to incorporate Ly emission from AGN in future work, which would boost emission and naturally induce scatter in the SFR-Ly relation.
In conclusion, we demonstrated that the current ISM emissivities appear reasonable given the simplistic model, but future work on the effective modeling of dust and a modified SFR-Ly relation (reflecting varying ISM environments and potential AGN presence) in both proportionality factor and its scatter might help to simultaneously diminish tensions of SB 0 and 0 in future work.

Future Directions
An important improvement for future Ly radiative transfer modeling in large-scale cosmological simulations such as TNG50 will be a treatment of complexities below the resolution scale (∼ 100 pc). Such a model could either implement explicit numerical subgrid models (Hansen & Oh 2006;) and retain our current principle of a parameter-untuned model. Alternatively, we could adopt an effective parameter-based model, incorporating observational findings to motivate Ly production and escape fractions from star-forming regions (Weiss et al. 2021).
In future work we will also use our coupling of ILTIS and TNG50 to examine two key areas: environmental imprints on largescale observations, and the information content of the spectral dimension. We found a significant redistribution of photons from starforming regions to large scales ( 100 kpc) which might affect the interpretation of Ly intensity mapping experiments. Furthermore, spectral modeling of Ly emission in cosmological volumes remains challenging given its multi-scale nature (Byrohl & Gronke 2020;Song et al. 2020). However, detailed and spatially resolved spectral information on Ly -halos is increasingly becoming available (Claeyssens et al. 2019;Leclercq et al. 2020), and promises to offer significant insight into the kinematics and small-scale structure of hydrogen gas in the CGM of dark matter halos across cosmic time.
Finally, we have here focused on lower mass star-forming galaxies with halo < 10 12 M . Above this threshold, the AGN in more massive galaxies are known to have a significant impact on the ionization state of the CGM and the Ly scattering processes occurring therein. Our current radiative transfer methodology does not account for local radiation fields from AGN, but this is a natural extension of ILTIS which will allow us to compare to the many rich observational data sets of Ly emission around quasars.

CONCLUSIONS
In this paper we develop a technique to perform full radiative transfer calculations to trace resonantly scattered Ly emission, and the Ly halos (LAHs) around galaxies, at 2 < < 5. We do so by post-processing the TNG50 cosmological magnetohydrodynamical simulation (Nelson et al. 2019b;Pillepich et al. 2019) of the Illus-trisTNG project.
This large volume offers a powerful statistical sample of thousands of LAHs across an unprecedented halo mass range and across diverse environments, together with a high resolution of order 100 physical parcsecs in the dense interstellar medium. Furthermore, our new radiative transfer code ILTIS (Byrohl et al., in prep) operates natively on the global Voronoi tessellation of the TNG simulation volume, incorporates both diffuse and galaxy emission, and selfconsistently accounts for attenuation within the intergalactic medium.
This allows us to carry out a detailed investigation into the origins, physical properties, and emission sources that shape Lyman-halos. At the same time, the realism of the underlying TNG50 simulation enables us to make quantitative connections between LAH and galaxy properties. Our main findings are: • Star-forming galaxies with 10 7 < ★ /M < 10 10.5 at 2 < < 5 emitting Ly photons are surrounded by extended Lymanalpha halos (LAHs). We present the stacked, median predictions for TNG50 LAHs as a function of galaxy mass and redshift. The radial surface brightness profiles of LAHs have a characteristic shape comprised of a rapid, exponential decline followed by a large-distance flattening. This flattening arises from the density structure probed by rescattering photons as well as in-situ diffuse emission.
• Scattered photons from star-forming regions are the dominant contributor to LAH profiles on typically observed scales 20 pkpc. Given the importance of scattered photons, we stress the need to use radiative transfer simulations or semi-analytic expressions capturing such behavior. At larger distances, contributions from diffuse emission via recombinations and de-excitations become equally important.
• On larger scales 30 pkpc we find that the flattening of LAH profiles is actually dominated by rescattered photons that originate from other nearby halos, rather than the primary halo itself. This proximity effect is boosted in high density environments, and should be observable. A careful reproduction of a survey's background subtraction and wavelength window are needed for a detailed comparison of this flattening to observations.
• Characterizing LAH sizes, we find that their half-light radii 1/2 are of order 5 pkpc at = 2, increasing to ∼ 15 pkpc at = 5. This signposts a significant redistribution of photons due to higher neutral hydrogen densities at higher redshifts. The exponential scale lengths 0 also increase slightly towards higher redshift for galaxy stellar masses of 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5. Neither 1/2 and 0 show clear trends with mass. In contrast, our fiducial model shows a strong positive correlation of the central surface brightness SB 0 with both stellar mass and redshift.
• While AGN radiative feedback adds significant heating and ionization to the surrounding CGM, subsequently boosting intrinsic emission from both collisional excitations and recombinations, we only find a marginal impact on the emission mechanisms' relative importance and the overall scale lengths 0 for LAHs with AGN activity in TNG.
• We compare the stacked, median LAH radial profile between TNG50 and data for galaxies with 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 from the MUSE UDF and find good qualitative agreement. We also demonstrate, by finding statistically consistent simulated analogs for individual observed profiles, that the simulation successfully reproduces the diversity of observed LAHs.
• For the quantitative comparison to observational results of LAH sizes as measured by half-light radii 1/2 and exponential scale lengths 0 , we find agreement at the level of ∼ 20% at = 3, with the 0 ∼ 4 pkpc of our simulated profiles slightly above those of the MUSE UDF data set ( 0 ∼ 3 pkpc). This difference increases with redshift, to a factor of two at = 5. Similarly, we find that the central surface brightness SB 0 of LAHs is in good agreement at low stellar mass, but such agreement becomes progressively worse towards higher ★ . Both tensions arise because LAH properties tend to correlate strongly with both galaxy mass and redshift in the simulations, but less so in the data. Barring important differences in the property distributions of the selected galaxies between the simulated and observed samples, we attribute these correlations to a number of simplifying assumptions in our Ly modeling, and discuss future improvements.
The extended Ly emission around galaxies and quasars at > 2 provides an insightful window into many aspects of galaxy formation and evolution. Here we have demonstrated the power of coupling a large-volume cosmological hydrodynamical simulation with a global Ly radiative transfer modeling. Future improvements in our treatment of unresolved small-scale gas structure and local radiation fields from AGN and stellar populations will enable interpretation of additional datasets and upcoming surveys, from large-scale intensity mapping experiments to highly detailed, spectral data from targeted IFU surveys.

DATA AVAILABILITY
Data directly related to this publication and its figures is available on request from the corresponding author. The IllustrisTNG simulations, including TNG50, are publicly available and accessible at www.tng-project.org/data (Nelson et al. 2019a Figure A1. Median stacked Ly surface brightness radial profiles at = 3 for a fixed stellar mass 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 at different resolution runs of TNG50. TNG50-1 has a mass resolution of baryon = 8.5×10 4 M , which increases by a factor of eight for each lower resolution run. Only star-formation is considered as an emission source, and scattered light contributions from other halos are ignored. The dashed lines correspond to the intrinsic emission radial profiles with a corresponding solid line for the profiles after radiative transfer. The shaded regions show the 16th to 84th percentiles. Differences between resolution levels are generally minor.

APPENDIX A: MODELING ASSUMPTIONS
In this appendix, we assess various modeling assumptions along with numerical convergence of our results. Unless stated otherwise we focus on the star-formation emission source only, and derive radial profiles neglecting rescattered light from other halos, as we limit ourselves to rerunning the radiative transfer on halos with stellar masses 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5.  Figure A2. Median stacked Ly surface brightness radial profiles (solid lines) at = 3 in TNG50 for a fixed stellar mass 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 for 1 (fiducial, blue) and 10 Monte Carlo photon packages per star-forming cell (orange). To show the largest possible discrepancy, we also include the leastconverged individual radial profile (dashed lines), which also demonstrates near perfect convergence in count . Our results are clearly insensitive to this parameter.

A1 Hydrodynamic Resolution
In Figure A1 we show the Ly surface brightness radial profiles at different hydrodynamical resolutions of TNG50 runs stacked at fixed stellar mass 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5. Overall, we find the stacked profiles to be very robust over the different resolution levels shown. For fixed stellar bins, we find the intrinsic emission to slightly expand with increasing resolution. For the processed photons, there is a slight flattening of the stacked profiles between 10 and 20 pkpc.
The resolution dependency appears to be significantly smaller than that in Zheng et al. (2011). We note that even the lowest resolution run presented here has a higher resolution than that study in the proximity of halos.

A2 Photon Package Count
In our fiducial runs, we spawn one photon package per Voronoi cell and emission source (star-forming region). Particularly for photons originating from central star-forming regions, there might be a convergence issue given the large volumes that are to be traced out at larger radii. We thus focus on convergence checks for the photon package count from star-forming regions.
In Figure A2 we show the radial profiles with varying count of initial photon packages. The median stacked profiles (solid) are already fully converged at our fiducial choice of photons = 1. As the most stringent test, we also plot the radial profile of the single halo with the largest sum of squared errors between both these two test runs -the least converged individual profile. In this case we also find only minuscule deviations confirming that our choice of the photon count leads to a well converged result.
We note that for halos with stellar masses between 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5, our fiducial setup spawns ∼ 35000 intrin- Here we contrast four different input photon wavelength offsets. In our fiducial model, photons are injected at the line-center within the rest-frame of the respective star-forming cell, indicated here as an offset of Δ = 0 km/s from the line center (blue). Three non-zero offsets of Δ = 100 km/s (orange), Δ = 200 km/s (green) and Δ = 500 km/s (red) are shown, where the shaded regions enclose the central 68 % of SB(r) at fixed r. We find little difference in the resultant LAH profiles between the different input spectra choices (Δ), except for the unrealistically large case of 500 km/s. sic photons: the intrinsically high resolution of TNG50-1 gives us already sufficient sampling.

A3 Input Spectrum
In our fiducial model, photons are injected at the Ly line-center in the rest-frame of the respective hydrodynamic cell. For diffuse emission, Equations (1)/(2), this is a fair assumption given that optical depths are moderate, although subgrid clumpiness would add complexity. However, in dense star-forming cells radiative transfer is significantly more complex due to small-scale dust, clumpiness, and ionization. Our radiative transfer simulations do not self-consistently capture these details. While attenuation on these scales only changes the overall radial profile normalization of the attenuated component, the spectral shape of emitted photons might have an effect on the radial profile shape itself.
As shown in similar cosmological simulations by Byrohl et al. (2019), Ly spectra often appear to have too much flux at wavelength larger than the line-center (too 'blue') compared to observations. We thus consider an additional redshift for the injected photons, which naturally leads to a more realistic 'red' spectrum. In Figure A3 we show the results of radiative transfer simulations where the redshift Δ from the line center has been varied between 0 km/s (fiducial model) and 500 km/s, considering emission from star-formation only. The corresponding intrinsic (dashed) and emergent (solid) Ly spectra are shown in Figure A4.
the frequency distribution at injection will be important for future studies, it is not crucial for our study of LAH profiles.

A4 Dust
In our fiducial model we have neglected dust. Comparison of luminosities with observations in Section 4.2.2 indicate that they are reasonable for the fiducial stellar mass range 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 without further explicit dust modeling, although this will be necessary for the high mass end. Dust modeling can either take the form of an additional attenuation factor of Ly emission, particularly of star-forming gas cells, or an explicit dust treatment in the Ly RT. In this section, we explore the option of including an explicit dust treatment in the Ly RT, as a preliminary study. We neglect clumpiness, such that dust is smooth at the resolution scale of the simulation. We use the model for Milky Way like dust (Weingartner & Draine 2001) as implemented in Behrens et al. (2019), which relies on the metallicity field in TNG50.
We expect that dust will primarily modulate the overall normalization of the radial profiles for the emission from the ISM given in Equation (3). Nevertheless, this can boost the relative importance of other emission mechanisms, and change the overall radial profiles. As dust content is related to gas-phase metallicity and galaxy mass, the impact of dust could alter the trends of LAH properties with mass, and be particularly important at higher galaxy masses.
In Figure A5 we show the impact of dust in our simulations at = 3 for emission from star-formation only. We contrast LAH radial profiles between the fiducial dust-free case (solid lines) and the dust included model (dashed lines). Our findings on the impact of dust are similar to those in Laursen et al. (2009). In particular, surface brightness is increasingly suppressed in overdense, dusty regions. In addition, regions of lower density are uniformly suppressed as dust limits the escaping (and then rescattering) contributions from star-forming regions.
Dust attenuation strongly scales with the stellar mass. For example, the median attenuation for the central surface brightnesses is less than  Figure A5. Median stacked Ly surface brightness radial profiles in various stellar mass bins at = 3 in TNG50. In this plot, we only consider contributions from star-forming regions, where the impact from dust is most severe. We show the results with (dashed lines) and without (solid lines) dust modeling. We generally find a significant suppression of flux that significantly increases for higher stellar mass halos. one order of magnitude for 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5, but roughly two orders of magnitude for 10.5 ≤ log 10 ( ★ /M ) ≤ 11.5.
The inner 10 pkpc are particularly suppressed, decreasing the central surface brightness values and thus increasing the exponential scale lengths. For the scale lengths of halos with 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5, we find an average 0 3.1 +1.4 −0.7 pkpc without dust and 5.1 +2.4 −1.2 pkpc with dust. Both the median scale length and its variance within the sample increase with dust. The increase in scale length with this dust modeling is even larger when incorporating other emission sources, as they become relatively more important (see Section 3). It is clear that future models will need to include at least a basic dust model.

A5 Impact of Local Ionizing Sources
Properly accounting for the impact of local ionizing sources on the temperature and ionization state of the CGM requires, ultimately, full radiation-hydrodynamical simulations. However, this remains computationally infeasible for the present combination of galaxy sample size, resolution and redshift.
In this section, we assess the impact of local ionizing radiation from AGN and stellar populations. An on-the-fly treatment (i.e. selfconsistently considering local sources' impact on cooling in each time step) is implemented for AGN radiative feedback in TNG. Hence, we can measure such ionizing sources' impact on the LAHs by comparing halo samples with and without SMBHs, as in Section A5.1. In the subsequent Section A5.2 we then discuss what impact ionizing radiation from stellar populations would have, given that this is not included in our models.

A5.1 Active Galactic Nuclei
Ionizing radiation from SMBH is modelled in TNG, which we will describe in the following. Note that we use the terms AGN and SMBH interchangeably here.
The intrinsic luminosity from black hole mass accretion at a rate BH in TNG is taken to be with 1 = 0.3, 2 = 0.07 (Hopkins et al. 2007) that results in a median obscuration of AGN esc ∼ 0.19 +0.10 −0.08 for all AGN in TNG50 at = 3. The sub-and superscript denote the interval of the central 99.73% of AGN. Higher obscuration factors typically occur for AGN in higher mass halos.
For AGN with low accretion / Eddington < 0.002, the radiative feedback is set to zero in the TNG model. The radiative output of AGN differs substantially with redshift. While at = 0.0 the vast majority of halos are inactive given this threshold, only 5% are radiatively inactive for the studied redshift range ≥ 2.
More details of the implemented radiative feedback model in TNG can be found in Vogelsberger et al. (2013);Weinberger et al. (2017Weinberger et al. ( , 2018. TNG incorporates "on-the-fly" ionizing radiation from SMBH in each time step as given in Equation (A2), which impacts gas cooling, temperature, and ionization state. For the gas cooling, these local sources' photoionization and photoheating are incorporated. The model assumes the bolometric intensity as SMBH UV,esc ∝ SMBH UV,esc / 2 at a distance of each nearby AGN under the assumption that the gas is optically thin to SMBH radiation. Photoionization and photoheating are calculated as a superposition of the local ionizing radiation from AGN and the uniform metagalactic background. For halos with a radiating black hole in TNG, the AGN field dominates over the UV background within the halo.
We will now assess the impact of the SMBHs' ionizing radiation on the halos' surrounding gas and the modifications to their Ly halos. To do this, we split the sample in halos with and without SMBHs. In TNG50, this is equal to splitting the samples by AGN radiative feedback activity.
In Figure A6 we show the median radial profiles for the neutral hydrogen fraction HI and temperature for the sample of halos in the stellar mass range 8.45 ≤ log 10 ( ★ /M ) ≤ 8.55 at = 3.0 in TNG50. For this mass range, we have a roughly equal sample of halos with (148) and without (122) SMBH activity (also see cutoff in Figure A9). In addition to the central 68 percentiles as shaded area, we also show the 2 median arcsecond radius given the virial radii, along the x-axis (vertical dotted line).
Note that only 1 of the 148 AGN is radiatively inactive. Hence, for our sample, AGN presence nearly always implies a significant ionizing flux in the surrounding CGM according to Equation (A2). As a consequency of the SMBH (radiative) feedback, we find a strongly suppressed neutral hydrogen content in affected halos. At the same time, the temperature is significantly higher outside the starforming regions at > 0.2 vir . Given the decreased neutral hydrogen fraction, we expect a boosted Ly emission from recombinations, and in particular a boosted Ly emission from collisional excitations given the strong temperature dependency of the latter.
In Figure A7 we radial profiles for = 3.0 and 8.45 ≤ log 10 ( ★ /M ) ≤ 8.55 split into their respective emission mechanisms and divided into the sample with and without SMBH activity similar to the dashed lines in Figure 8. Without scattering, we generally find collisional excitations to be the dominant emission mechanisms for the majority of halos above 7 pkpc irrespective of SMBH activity. Particularly recombinations are strongly boosted in the presence of a SMBH for 0 20 pkpc. Interestingly, while the overall emission from collisional excitations is boosted, the relative fraction decreases not just because of the higher recombination rates but also due to a larger fraction of emission from star-forming regions.
To understand latter point, we stress that the underlying samples  Figure A7, but for the processed emission (i.e. with scatterings). We observe a similar behaviour for both samples as in Figure 8. Slight differences in shape seem to arise from boosted collisional excitations and recombinations for 0 < 20 pkpc relative to the contribution from starforming cells, leading to a slightly flatter slope. Star-formation remains the dominant contributor to the radial profiles across the shown radii. differ. Particularly, the sample hosting SMBHs, even though with a similar stellar mass, are typically more massive with a median total halo mass log 10 ( ℎ /M ) higher by about 0.2 dex and a significant scatter towards higher masses. The high mass objects show a larger amount of satellites contributing to the additional emission from star-formation in the halos' outskirts. Figure A8 shows the processed median radial profiles, i.e. after running the Ly radiative transfer for the intrinsic emission shown in A7. While all emission mechanisms remain boosted in absolute terms when SMBH activity is present, collisional excitations and particularly recombinations become more important in relative terms. However, irrespective of SMBH activity the radial profiles remain dominantly sourced by emission from star-forming cells. At 0 20 pkpc, contributions from collisional excitations become close to equally important to photons from star-forming regions. However, this trend is the same for both the AGN and no-AGN sample.
From the median profiles it appears that there is a slight additional flattening for the AGN hosting subsample compared to the no-AGN sample. To quantify the difference of the LAHs' shape, we calculate the individual exponential scale lengths 0 and find an overall increase of roughly 11% from 3.6 +1.2 −0.9 pkpc (no AGN) to 4.0 +1.4 −1.0 pkpc (AGN). This increase is however not solely driven by the ionizing AGN radiation but also by the different halo population, given the fixed stellar mass ranges of the sample. If we instead constrain the total halo mass to be within a fixed range of 10.6 ≤ log 10 ( ℎ /M ) ≤ 10.8, which corresponds to the central ∼ 68% of the non-AGN 8.45 ≤ log 10 ( ★ /M ) ≤ 8.55 sample, we obtain a reversed trend with 0 being 3.8 +1.4 −0.9 pkpc (3.6 +1.7 −0.8 pkpc) for the sample without (with) AGN.
In conclusion, the impact of ionizing radiation from SMBHs in TNG appears to be significant in terms of ionization state, temperature and thus intrinsic Ly emission of the gas, but the findings on the LAH shape, particularly through the scale radius 0 , remain unchanged irrespective of modelled SMBH activity. We therefore do not split our sample within the main body with respect to the presence of AGN.
Note that here, we only discuss the AGN's impact through its ionizing budget onto the surrounding. In particular, we do not discuss nor implement a description for the unresolved Ly emission from AGN itself at this point. We would expect that this emission would scatter outwards and produce a contribution to the LAH of similar shape to that of the star-formation, given the concentrated emission source scattering into the surrounding CGM. As differences in the intrinsic emission in terms of spatial and spectral distribution exist, this would need to be explored in future investigations.

A5.2 Stellar Populations
TNG does not incorporate the ionizing flux from local stellar populations as is done for SMBHs. Hence, the possible impact of those sources on our predicted LAH profiles cannot be assessed within the existing simulations. However, we can derive an estimate of the upper limit of this effect, by comparing to the local ionization from AGN.
The ionizing luminosity escaping from the stellar populations, adopting a formulation consistent with our model, is given as: As in Equation (3), we assume no dust and Case B recombination with B = 0.68 being the conversion factor from ionizing to Ly photons.
,UV / ,Ly is the ratio of the average ionizing photon energy in the population and the Ly line transition energy. We assume UV,esc = 0.1 in Equation (17) of Dĳkstra (2019).
In Figure A9 we show the escaping UV radiation from stellar sources and AGN as a function of the stellar mass of the host halo at = 3.0 in TNG50. We find AGN activity for the majority of halos with a stellar mass above log 10 ( ★ /M ) 8.5. Above this mass, we find an approximate power law UV,esc = 0 · ★ 10 8.5 M with 0 = 10 41.7 and = 1.31 ( 0 = 10 41.1 , = 0.99) for AGN (SF). We thus find that AGN UV,esc 6 · SF UV,esc at = 3.0 for the bulk of halos with AGN activity and a growing disparity between the luminosities with larger mass given the larger slope for AGN. The discrepancy  Figure A9. We show the median of the UV luminosity escaping the galaxy at a given stellar mass of the host halo at = 3.0 in TNG50. The UV luminosities for star-formation (SF) and SMBH are derived from Equations (A2)/(A4) respectively. grows (shrinks) at lower (higher) redshifts across the stellar mass range.
Changing UV,esc or adopting different assumptions concerning e.g. metallicity, binary fraction or IMF can boost the stellar populations' escaping UV luminosity. Given the large margin at = 3.0 between AGN and stellar luminosities, qualitative findings should be robust. However, such adjustments might lower the redshift at which stellar populations' UV luminosity becomes dominant into the upper studied redshift range. For the given UV luminosities, the AGN ionizing flux integrated over the halo population dominates over the ionizing flux from stellar populations. This ratio peaks at = 2 in TNG50, where AGN are most relevant, while stellar populations start to dominate the overall ionizing flux at higher redshifts 6. We found that the ionizing flux of SMBHs, if present, dominates over the ionizing flux from stellar populations in the same host halo. A consideration of the impact photoionization and photoheating thus needs to primarily address the impact of AGN, for which a simplified description is indeed implemented in TNG50.
In our fiducial sample with stellar masses of 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5, 87% of halos host at least one SMBH with a radiation field according to Eqn. A2. Therefore, our analysis contains the primary ionization source. As demonstrated in Appendix A5.1, while the gas becomes significantly hotter and ionized, the LAH size measurements and our conclusions remain largely unaffected. As a consequence, we also find the inclusion of stellar ionizing sources, with overall weaker UV luminosities, would have a small impact on our overall findings.

A6 Impact of a Wavelength Window
For simplicity, we ignored spectral information throughout this paper and instead sum all photons scattering last within ±100 pkpc depth around the halos' center. Here, we show the difference when instead using a simple spectral bandwidth from the photons' spectral information incorporating Hubble flow, peculiar velocity and spectral diffusion. Just as in Appendix A5, we consider contributions from all three implemented emission sources and do not neglect the rescattered light from other halos as done in the other appendices. We note that a fair comparison will need to carefully match the various  Figure A10. Median stacked radial Ly profile for galaxies with stellar masses 8.5 ≤ log 10 ( ★ /M ) ≤ 9.5 at = 3 in TNG50. We decompose the radial profiles into their two dominant origins from central galaxies (blue) and other halos (orange). For the solid line, we show the fiducial method throughout this paper that ignores spectral information and integrates all photon contributions that scatter last in the ±100 pkpc around the halo's position. The dashed and dotted lines show the contributions when incorporating spectral information in a ∼ 12 Å (dashed) and ∼ 16 Å (dotted) observed wavelength window around the Ly line center. The central surface brightness is slightly suppressed when considering spectral information due to some emitters with a spectral diffusion from the central galaxies exceeding the imposed wavelength window. At large radii, the wavelength windows lead to a stronger flattening due to the larger physical depth.
methodologies used in observational studies. For example Leclercq et al. (2017) uses an adaptive spectral bandwidth per emitter. In Figure A10 we show the median stacked radial profiles using the fiducial depth integration based on the last scattering of photons (solid line) opposed to a fixed spectral bandwidth of ±6.1 Å (dashed line) and ±8.1 Å (solid line) for the overall radial profile and split into the dominant emission origins (central galaxies in blue, other halos in orange). Ignoring spectral diffusion, we effectively model an adaptive spectral bandwidth to capture potentially large spectral diffusion. Thus, a fixed spectral bandwidth leads to a slightly suppressed median radial profile in Figure A10.
The scale lengths 0 = 4.2 +2.0 −1.1 pkpc (±6.1 Å) and 4.0 +1.9 −1.0 pkpc (±8.1 Å) for the fixed bandwidth window slightly increase compared to the fiducial setup due to the exclusion of Ly contributions that diffused outside of these windows.
In Figure 7 and 12 we demonstrate a flattening of the radial profiles that is dominated by contributions from other halos. This flattening will therefore be heavily influenced by the environment in a chosen field of view and the methodology for its detection and stacking. A fair comparison will thus require a thorough reproduction of factors such as field overdensity, source masking and chosen wavelength depth along the line of sight. Generally, a fixed bandwidth as in Figure A10 causes a larger flattening due to the larger relative background contribution from other halos which is expected to scale roughly linear with the integration depth. The integration depth from the differential Hubble flow corresponds to 2.4 pMpc and 3.2 pMpc respectively compared with the 0.2 pMpc in the fiducial setup. Future observational comparisons will benefit from close attention to the chosen wavelength depth along with a careful account of the background subtraction.