The cosmic web in Lyman-alpha emission

ABSTRACT

We develop a comprehensive theoretical model for Lyman-alpha (Ly α) emission, from the scale of individual Ly α emitters (LAEs) to Ly α haloes (LAHs), Ly α blobs (LABs), and Ly α filaments (LAFs) of the diffuse cosmic web itself. To do so, we post-process the high-resolution TNG50 cosmological magnetohydrodynamical simulation with a Monte Carlo radiative transfer method to capture the resonant scattering process of Ly α photons. We build an emission model incorporating recombinations and collisions in diffuse gas, including radiative effects from nearby AGN, as well as emission sourced by stellar populations. Our treatment includes a physically motivated dust model, which we empirically calibrate to the observed LAE luminosity function. We then focus on the observability and physical origin of the z = 2 Ly α cosmic web, studying the dominant emission mechanisms and spatial origins. We find that diffuse Ly α filaments are, in fact, illuminated by photons that originate not only from the intergalactic medium itself but also from within galaxies and their gaseous haloes. In our model, this emission is primarily sourced by intermediate mass haloes (10¹⁰–10¹¹ M_⊙), principally due to collisional excitations in their circumgalactic media as well as central, young stellar populations. Observationally, we make predictions for the abundance, area, linear size, and embedded halo/emitter populations within filaments. Adopting an isophotal surface brightness threshold of 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻², we predict a volume abundance of Ly α filaments of ∼10⁻³ cMpc⁻³ for lengths above 400 pkpc. Given sufficiently large survey footprints, detection of the Ly α cosmic web is within reach of modern integral field spectrographs, including MUSE, VIRUS, and KCWI.

1 INTRODUCTION

Within the ΛCDM cosmological paradigm, gravitationally unstable initial matter density fluctuations evolve into a filament-dominated structure on large scales: the cosmic web (Bond, Kofman & Pogosyan 1996). At late times, the majority of dark matter haloes, as well as galaxies, reside in the filaments and nodes of this cosmic web (Meiksin 2009). The same is true for the majority of dark matter and baryons, including diffuse gas. As a result, the formation and evolution of galaxies is mediated in large part by their gaseous environments, including gas gravitationally bound within dark matter haloes – the circumgalactic medium (CGM; Tumlinson, Peeples & Werk 2017).

The large-scale filaments of the cosmic web can indirectly be observed through galaxy clustering in galaxy redshift surveys (e.g. Colless et al. 2001; Abazajian et al. 2009). At high redshift, a direct detection of these filaments is possible via absorption by the Lyman-alpha (Ly α) line of neutral hydrogen. In this case, spectra of background sources, mainly quasars, probe the hydrogen distribution in the intergalactic medium (IGM) along the line of sight (Gunn & Peterson 1965; Meiksin 2009). In recent years, high sampling density of such quasar spectra has enabled the reconstruction of the 3D density field of neutral hydrogen (Lee et al. 2014, 2018; Newman et al. 2020). However, the coarse resolution of the order of megaparsecs makes it difficult to resolve the filamentary structure of the cosmic web, a limitation inherited from the sparseness of background quasars on the sky.

In contrast to absorption, the Ly α emitting cosmic web offers a complementary approach. However, direct imaging of large-scale Ly α filaments (LAFs) remains challenging given the low emissivities of the diffuse gas (Gallego et al. 2018). For denser environments, Ly α emission is already a frequently used tracer of cold gas. For example, Ly α emission is commonly used to identify high-redshift galaxies in blind surveys (Cowie & Hu 1998). In targeted surveys, extended emission with sizes of ∼10–100 pkpc around massive galaxies has been detected for decades (McCarthy et al. 1987; Heckman et al. 1991; Steidel et al. 2000). More recently, the Ly α emission around smaller star-forming galaxies, tracing the CGM around these objects, has been revealed on scales of ∼10 pkpc – first through narrowband stacking (Hayashino et al. 2004; Steidel et al. 2011; Matsuda et al. 2012; Momose et al. 2014; Kakuma et al. 2021) and then through integral field spectroscopy (Wisotzki et al. 2016; Leclercq et al. 2017; Lujan Niemeyer et al. 2022b).

The latest observations of Ly α emission around star-forming galaxies show flattened extended radial profiles (Wisotzki et al. 2018; Kakuma et al. 2021; Kikuchihara et al. 2022; Lujan Niemeyer et al. 2022a), potentially hinting at the faint Ly α glow of the cosmic web. Large filamentary structures with extents of ≳1 pMpc have been detected by targeting known overdense fields: the SSA22 protocluster (Umehata et al. 2019; Herenz, Hayes & Scarlata 2020) and the Hyperion proto-supercluster (Huang et al. 2022). Aligned stacking of galaxy pairs at z ∼ 2 shows extended Ly α emission from the CGM, but no signal from intergalactic scales (Gallego et al. 2018). However, the Ly α cosmic web at z ≳ 3 has potentially recently been detected in a blind survey using integral field spectroscopy (Bacon et al. 2021).

The Ly α emission line of the neutral hydrogen atom is a promising tool to study large-scale structure and can be probed at redshifts z ≥ 2 with current ground-based instruments on ∼8–10 m telescopes such as MUSE, KCWI, and VIRUS (Bacon et al. 2010; Morrissey et al. 2018; Hill et al. 2021) on the VLT, Keck, HET telescopes, respectively. Upcoming ∼30 m class telescopes will further enable the detection of LAFs. The lower end of the redshift range at z ∼ 2 is particularly promising given its favourable cosmological surface brightness (SB) dimming scaling as (1 + z)⁻⁴. For example, the majority of filament candidates in Bacon et al. (2021) are located towards the lowest accessible redshifts, around z ∼ 3. However, the overall redshift trend of filament detectability depends on a complex evolution of the physical properties of filaments, including their density, temperature, and ionization state, as well as galaxy clustering, global star-formation, dust content, and IGM opacity.

Predictions for the observability of cosmic web filaments have been made for intensity mapping (Silva et al. 2013; Silva, Kooistra & Zaroubi 2016; Heneka, Cooray & Feng 2017) as well as direct observation (Elias et al. 2020; Witstok et al. 2021). However, this requires comprehensive and accurate emission models for Ly α photons. The physical processes involved include emission from excitations and recombinations in the diffuse gas, and the effective emission that arises due to the radiative output of young stars during the process of star formation, as well as due to radiation from AGN.

The uncertainties in the modelling of emission mechanisms are further complicated by the complex radiation transfer that Ly α photons experience (Neufeld 1990, 1991; Hansen & Oh 2006). The Ly α emission line is resonant and optically thick in astrophysical environments, leading to numerous scatterings before photons eventually escape, or are destroyed. This causes substantial spatial and spectral redistribution of photons, making this forward modelling step imperative from the simulation point of view, as well as complicating the interpretation of any observed emission.

A common approach is to post-process cosmological (radiation-) hydrodynamical simulations with Monte Carlo-based Ly α radiative transfer codes (e.g. Cantalupo et al. 2005; Laursen & Sommer-Larsen 2007; Goerdt et al. 2010; Kollmeier et al. 2010) to study different Ly α observables such as Ly α emitter (LAE) clustering (Zheng et al. 2011; Behrens et al. 2018; Byrohl, Saito & Behrens 2019), spectral signatures from the IGM (Byrohl & Gronke 2020; Park et al. 2022), and extended emission.

Many recent theoretical studies of extended Ly α emission have included scattering effects and focused on CGM scales (Lake et al. 2015; Gronke & Bird 2017; Behrens et al. 2019; Smith et al. 2019; Byrohl et al. 2021; Mitchell et al. 2021). However, investigations dedicated to the cosmic web in Ly α emission have universally neglected the impact of radiative transfer (Elias et al. 2020; Witstok et al. 2021). In that context, Witstok et al. (2021) conclude that observations of the cosmic web at lower redshifts are most promising in overdense regions, where emission is dominated by the haloes and galaxies within filaments. In this regime, collisional excitations produce more Ly α photons than recombinations. Simultaneously, Elias et al. (2020) suggest that Ly α SB predictions can be used to constrain the underlying galaxy formation model physics. However, the lack of a quantitative assessment to date for the occurrence of these filaments hinders observational constraints on Ly α radiative transfer simulations of the cosmic web.

In this study, we model and characterize the Lyman-α cosmic web in emission. To do so, we adopt the high-resolution, large-volume TNG50 cosmological galaxy formation simulation. We focus on redshift z = 2 and post-process the original simulation output with our sophisticated Monte Carlo radiative transfer method. We furthermore introduce a physically motivated dust rescaling model calibrated against the observed Ly α luminosity function (LF).

This paper is organized as follows: In Section 2, we introduce the TNG50 simulations, the Ly α radiative transfer code, the underlying emission model, and our analysis pipeline. In Section 3, we present results regarding global Ly α-related properties and a study of filamentary Ly α structures. In Section 4, we discuss our results for the dominant physical mechanisms, which light up the Ly α cosmic web, the origin of Ly α photons from filaments, and the detectability of the cosmic web with current and upcoming Ly α emission surveys.

2 METHODOLOGY

2.1 TNG50

The TNG50 simulation (Pillepich et al. 2019; Nelson et al. 2019b) is the highest-resolution simulation of the IllustrisTNG suite, a series of three large-volume magnetohydrodynamical cosmological simulations (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018a; Pillepich et al. 2018b; Springel et al. 2018) with a volume of ∼50³ cMpc³. The simulations were run with the AREPO code (Springel 2010), which solves the coupled equations for self-gravity and ideal, continuum magnetohydrodynamics (Pakmor, Bauer & Springel 2011) discretizing space using an unstructured Voronoi tessellation. The TNG galaxy formation model (Weinberger et al. 2017; Pillepich et al. 2018a) includes a treatment for the majority of physical processes shaping galaxy formation: primordial and metal-line cooling, heating from ultraviolet background (UVB) radiation, star formation above a density threshold, stellar feedback-driven galactic winds, stellar population evolution and chemical enrichment from supernovae Ia, II and AGB stars, and the seeding, merging, and growth via accretion of supermassive black holes (SMBHs).

The temperature and ionization state of gas, which is relevant for Ly α emission and scattering, is computed following the primordial cooling network of Katz, Weinberg & Hernquist (1996) with additional metal-line cooling from cloudy cooling tables. In addition, a heating and ionization term arises from the assumption of a uniform, time-varying UVB using the intensities given in Faucher-Giguère et al. (2009, FG11 update). An additional local ionization field is introduced for active galactic nuclei (AGNs) up to three times the hosting haloes’ virial radius. The underlying AGN luminosities are proportional to the accretion rates above a certain accretion threshold, modulated by an obscuration factor based on Hopkins, Richards & Hernquist (2007), assuming optically thin gas (Vogelsberger et al. 2013). Radiation from the UVB and AGN is attenuated on the fly, according to Rahmati et al. (2013), to account for self-shielding. The additional AGN radiation field is important for the gas state and subsequent Ly α emission. Particularly, at z ≥ 2, substantial changes arise given the high accretion rates of SMBHs (see Byrohl et al. 2021). Note that ionizing radiation escaping from local stellar sources is not included in the model.

TNG50 has a gas mass resolution of m_bayron = 8.5 × 10⁴ M_⊙ and a dark matter mass resolution of m_DM = 4.5 × 10⁵ M_⊙. This roughly corresponds to a spatial resolution of ∼100 physical pc in the ISM. The simulations use a set of cosmological parameters consistent with recent results by the Planck Collaboration (2016), namely Ω_{Λ, 0} = 0.6911, Ω_{m, 0} = 0.3089, Ω_{b, 0} = 0.0486, σ₈ = 0.8159, n_s = 0.9667, and h = 0.6774.

2.2 Lyman-alpha emission and radiative transfer

We calculate the Ly α radiative transfer in post-processing using our updated, light-weight line emission radiative transfer code (originally introduced in Behrens et al. 2019). It propagates Monte Carlo photon packages according to the given gas structure, accounting for scattering and destruction. Upon each scattering, we calculate the luminosity contribution escaping towards predefined observers (the ‘peeling-off’ algorithm; Whitney 2011). Our radiative transfer method supports a range of geometries, including the underlying Voronoi tessellation of TNG50. Using the Voronoi tessellation of the TNG50 snapshot, we capture the native resolution of the simulation. For this work we compute and ray-trace through the global (entire snapshot) mesh at once, in order to self-consistently capture environmental and large-scale IGM effects in the radiative transfer (as introduced in Byrohl et al. 2021).

The Monte Carlo implementation closely resembles that of other codes used for Ly α radiative transfer (e.g. Ahn, Lee & Lee 2001; Zheng & Miralda-Escudé 2002; Cantalupo et al. 2005; Semelin, Combes & Baek 2007; Gronke & Dijkstra 2014; Smith et al. 2015; Michel-Dansac et al. 2020) and many more. Differences to other codes primarily arise from the technical implementation, focusing on a highly performant, distributed, geometry-agnostic code, being able to run on both structured and unstructured grids, particularly the Voronoi tessellation underlying the TNG50 simulation here. Our code has been verified against well-studied setups such as a static or expanding homogeneous setups (Harrington 1973; Neufeld 1990; Dijkstra, Haiman & Spaans 2006).

We follow the Ly α emission model as introduced in Byrohl et al. (2021), where we include the emission of diffuse gas by recombinations and collisional excitations and emission from star-forming regions. The emission model for the diffuse gas is unchanged, with luminosity densities given by

and

which scale with the number density of electrons (n_e), neutral (n_HI), and ionized hydrogen (n_HII). The temperature-dependent recombination coefficient α(T) and collisional excitation coefficient γ_1s2p(T) are taken from (Scholz et al. 1990; Scholz & Walters 1991; Draine 2011). Recombinations are calculated assuming case-B, i.e. under the assumption of optical thickness of Lyman-series photons.

For the emission from dense gas around star-forming regions, we update the previous model and do not emit Ly α radiation based on the instantaneous star formation rate of gas cells. Instead, we model the Ly α emission from stellar populations as follows. We calculate the ionization rate |$\dot{N}_{\mathrm{ion}}$| of all stars according to their age, mass, and metallicity with bpass (Stanway & Eldridge 2018), assuming a Chabrier initial mass function (Chabrier 2003). From this, we derive the nebular emission of Ly α under the case-B assumption via

where the conversion factor f_B was fixed to a reasonable value of 0.68 while the UV escape fraction is set to f_{UV, esc} = 0.1, physically due to the gas surrounding each stellar population. This approach becomes feasible given the high resolution of TNG50, which allows a reasonable sampling of young <10 Myr stellar populations, a common problem for cosmological simulations at lower resolution (Trayford et al. 2017; Nelson et al. 2018a).

In addition to the ionization rate |$\dot{N}_{\mathrm{ion}}$|⁠, we calculate the luminosity per wavelength of the stellar continuum l_cont at the Lyman-α resonance wavelength with a Gaussian smoothing kernel of σ = 50 Å. We then compute the intrinsic rest-frame equivalent width (REW) as

For the diffuse emission mechanisms, we spawn one photon per gas cell. For stellar emission, we spawn 11 000 photons per 10⁴² erg s⁻¹ in luminosity, with a minimum of three photons per stellar particle.¹ This is a significant boost in photon count from previous work in Byrohl et al. (2021), which already demonstrated convergence on the Ly α halo (LAH) scale. The luminosity-dependent sampling was found to further reduce noise.

In the radiative transfer code, we use the temperature and neutral hydrogen density as directly inferred from the TNG50 snapshot data. The effective temperature and average hydrogen density cannot be used in a straightforward way in star-forming cells invoking a sub-grid effective equation of state for a two-phase ISM (Springel & Hernquist 2003). For these cells, we adopt the temperature and density values from the snapshot but do not generate any emission from recombinations and excitations.

The role of dust for Ly α radiative transfer is important (Laursen, Sommer-Larsen & Andersen 2009; Byrohl et al. 2021). Its impact is particularly susceptible to the unresolved small-scale structure (Gronke et al. 2017) due to the high-optical depths within individual cells. We have developed new models and strategies for incorporating dust, and for the present work, we have decided on an empirical calibration strategy. In particular, we introduce an effective dust attenuation by rescaling the stellar luminosity contributions as described in Section 2.4.

All photons are emitted at the Ly α line-centre frequency with a random initial direction. Frequency shifts with ≤200 km s⁻¹ injected on the ISM scale have little impact on the radiative transfer for the described setup (see appendix of Byrohl et al. 2021).² Larger wavelength shifts could, however, largely change the outcome. Generally speaking, a large spectral redshift (blueshift) decreases (increases) the redistribution of photons into their surroundings.

2.3 Scattered Lyman-alpha photon properties

For each photon contribution, we keep track of the following details: the spatial location of last scattering, luminosity, and wavelength. In addition, we save the global Voronoi cell indices at the points of initial emission and last scattering. This allows us to trace and analyse the photon contributions with respect to the underlying simulation, including its halo and galaxy populations and their properties.

We distinguish between two different sets of photons:

• intrinsic: Ly α photons as emitted from gas cells and stellar particles, intentionally neglecting further gas interaction, i.e. scatterings and destruction.

• scattered: Ly α photons that escape towards the observer after each scattering of a previously emitted Ly α photon, including attenuation and scattering.

Only the scattered photons represent observable Ly α signatures. However, the intrinsic photons, which are only accessible theoretically, enable us to study the origin of Ly α emission and the impact of the Ly α radiative transfer.

By identifying the initial and final Voronoi gas cell for each photon, we classify each photon into exactly one of the following five spatial categories, at the time of emission as well as last scattering:

• IGM: IGM gas, i.e. does not belong to any collapsed halo.

• outer halo: gas that is part of a dark matter halo, but gravitationally unbound, i.e. on the outskirts.

• CGM: gas in the halo, gravitationally bound to the central galaxy, and outside 10 per cent of the halo virial radius.

• central: gas in the halo, gravitationally bound to the central galaxy, and inside 10 per cent of the halo virial radius.

• satellite: gas gravitationally bound to a satellite galaxy, which is within a larger host halo.

These categories rely on the friends-of-friends halo and subfind subhalo identification algorithms (see Nelson et al. 2019a).

Furthermore, we can study the physical gas state, as well as galaxy and halo properties, at two distinct times:

• at origin, using the gas cell, where the Ly α photon was initially emitted, or

• at last scattering, using the gas cell within which the Ly α photon finally escaped to the observer.

We create 2D SB projections with a depth of 5.7 Å in the observed frame (R ∼ 650), similar to the HET-VIRUS resolution (Hill et al. 2021) and 50 per cent higher than the spectral binning used in Bacon et al. (2021) to study the Ly α cosmic web at higher redshift.

For SB maps covering the entire extent of the simulation box, we use a map with 4096² pixels corresponding to resolution elements of ∼0.5² arcsec². This resolution is also used for all reduced statistics, such as filament sizes and shapes. In such cases, we create and use as many projections of the given depth as possible, equally spaced and non-overlapping along the line of sight. When zooming into sub-regions, we use resolution elements of 0.2² arcsec².

2.4 Observational calibration

The escape or destruction of Ly α emission in the ISM strongly depends on small-scale gas structure. The relevant scales are at least partially unresolved at the resolution of TNG50, motivating us to develop a sub-grid attenuation model that is empirically calibrated. Instead of explicitly modelling the abundance, distribution, and physics of dust, we instead rescale the Ly α luminosities in the ISM. Given the significantly lower optical depths of dust in more diffuse gas, we do not rescale the contributions from recombinations and collisional excitations, which occur only in non star-forming gas.

The rescaling of emission arising from nebular emission around stellar populations is done as follows. We first calculate a Ly α luminosity for each galaxy by summing the luminosities of all gas and stars bound to the halo, restricted to a circular aperture with a diameter of 3 arcsec. This procedure is done for the scattered photons of a full radiative transfer run, i.e. fully observable Ly α photons. Next, all photon contributions from the ‘stellar origin’ are rescaled downward to match observational constraints. In particular, we use the z = 2 Ly α LF as our only calibration.

We propose a coarse but physically motivated attenuation model intended to represent Ly α destruction by dust. We rescale the original luminosity of each stellar population L_{orig, i} as

where f_j is the per-halo rescaling factor for halo j,

that is set by the attenuating optical depth τ_j for that particular halo. These optical depths are assumed to follow a simple relation with host halo properties, including scatter. In particular, they are drawn from a Gaussian with mean μ_τ and standard deviation σ_τ, similar in spirit to Inoue et al. (2018). We parametrize the mean value using the mass-weighted average gas metallicity Z_j of the host halo:

and assume the scatter increases with the mean σ_τ = sμ_τ. The three free fit parameters of our model are therefore: Z_⋆, β, and s. Here, τ_j can be interpreted as the effective optical depth experienced by the Ly α photons. In contrast, the three fit parameters do not map to physical properties, while they can be interpreted in terms of their effective impact on the Ly α luminosities: Z_⋆ characterizes the overall scaling factor of metallicities given an optical depth τ, β the power-law coefficient of the τ(Z) relation, and an increasing s reflects the scatter in the τ(Z) relation relative to its mean. However, it can also capture potential disagreements of TNG50’s star formation rates with observations, compensate for modelling deficiencies of the diffuse emission, and encapsulate the impact of different |$L_{{\mathrm{Ly}\alpha }}\left(\dot{N}_{\mathrm{ion}}\right)$| relations.

We minimize the mean-squared error between the mocked TNG50 Ly α LF and the observational data by Konno et al. (2016). We impose a REW cut of 20 Å, which is most appropriate for comparison with Konno et al. (2016) and only fit observational data points between 10⁴² and 10⁴³ erg s⁻¹. At z = 2, the global best-fitting model occupies a well-defined minimum with Z_⋆ = 10^−4.78, β = 0.47, and s = 0.32. Note that other local minima exist, although they provide worse fits, and we expect that the parameter values can vary substantially with larger galaxy samples, simulations, and redshift.

The intrinsic equivalent widths are modelled using the stellar continuum estimate as described in Section 2.2 attenuated by dust. For the dust attenuation of the continuum, we compute the optical depth

for dust similar to Nelson et al. (2018a), taking attenuation strictly proportional to the gas metallicity Z_g and hydrogen column density N_H with Z_⊙ = 0.0127 and N_{H, 0} = 2.1 × 10²¹ cm⁻². We adopt the attenuation curve |$\left(\frac{A_\lambda }{A_V}\right)$| from Calzetti et al. (2000). The optical depth is calculated along each line of sight by ray-tracing through the metallicity and hydrogen density in each Voronoi cell.

In all cases, we use the rescaled luminosities throughout the following analysis by rescaling all photons, whether scattered or intrinsic, originating from stellar populations in halo j by f_j provided by our best-fitting model. Overall, this calibration ensures that our Ly α emission model is reasonable, and so fulfills a necessary condition to study the observability of the cosmic web hosting the Ly α emitting objects contained in the LF.

3 RESULTS

We introduce the outcome of our theoretical modelling with Fig. 1, which shows the Ly α SB across cosmological scales for TNG50 at z = 2. We project through a relatively narrow slice depth of 5.7 Å and adopt our fiducial emission model. The volume is suffused with Ly α light across a range of scales, from compact emission sources to elongated filamentary structures spanning megaparsecs to tens of megaparsecs in extent. SB levels within these filaments vary significantly, with brighter structures exceeding 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻².

3.1 The Lyman-alpha luminosity function and relation to galaxy and halo properties

The realism of our Ly α modelling results bears close inspection. We first assess the outcome by considering the LF, the observable against which we calibrate the emission model. In Fig. 2, we show the LF for the calibrated luminosities from scattered photons (blue) and for the intrinsic uncalibrated luminosities (orange). The calibrated (i.e. ‘rescaled’, two terms we use interchangeably) LF is in good agreement with the observational LF of Konno et al. (2016). At higher luminosities L > 10⁴³ erg s⁻¹, where AGN contamination in the observed sample increases (de La Vieuville et al. 2019) and the volume of TNG50 is too small to include these rare systems, we no longer (aim to) match observational data points. At luminosities below 10⁴² erg s⁻¹, the calibrated TNG50 LF gradually flattens with a plateau at ∼10⁴¹ erg s⁻¹ before once more steepening, and then finally turning over and decreasing below ∼10³⁷ erg s⁻¹ (not shown). The plateau roughly coincides with star formation ceasing to be the dominant emission mechanism in this luminosity range.

The global Ly α luminosity density inferred from the LF integrated for all L > 10^41.75 erg s⁻¹, i.e. where constrained by the data points, is 5.2 · 10³⁹ erg s⁻¹ cMpc⁻³, which is roughly a factor of two smaller than the Konno et al. (2016) data itself. This is due to the simulation’s drop-off at the high-luminosity end compared to the distinct AGN induced bump in observations. We also point out that the majority of the global Ly α luminosity density is below Konno et al. (2016) lower limit of L_min = 10^41.75 erg s⁻¹ with a total luminosity density from the LF of 1.33 × 10⁴⁰ erg s⁻¹ cMpc⁻³, when we also neglect any REW threshold.

In Fig. 3, we show a fundamental outcome of the model: the Ly α mass-observable relations. Specifically, median Ly α luminosity as a function of halo mass (top panel) and stellar mass (bottom panel). The upper grey line shows the median of the total, intrinsic, uncalibrated luminosities. All other lines show the luminosity of our fiducial model after the calibration method and radiative transfer calculation. In particular, the black line shows the median of the total, scattered calibrated luminosities. Shaded regions show the central 68 per cent and 95 per cent of luminosities for scattered photons. Coloured lines show the mean luminosities, for scattered photons, separating into the three physical origins: star formation i.e. young stellar population sourced (blue), collisions (orange), and recombinations (green).

All solid lines adopt our fiducial aperture and definition: summing photons from bound gas within an aperture radius of 1.5 arcsec. The dashed lines include contributions outside this aperture radius. The virial mass corresponding to this aperture radius is 5 × 10⁹ M_⊙, above which the dashed lines rise above the solid lines in the upper panel. We find that for the most massive haloes, the majority of escaping Ly α photons originate from radii beyond this aperture.

Overall, we see that Ly α luminosity monotonically increases with mass. The steepest scaling between halo mass and Ly α luminosity occurs between 10⁹ and 10¹⁰ M_⊙. At lower masses, the relation flattens out as scattered photons from other haloes start to dominate. At higher masses, the relation flattens significantly, in part due to the halo extent exceeding the fiducial aperture radius (compare to dashed lines), but primarily due to the impact of dust. Ly α luminosity after radiative transfer and calibration decreases substantially compared to intrinsic uncalibrated values, up to 2 dex for the most massive haloes.

Our Ly α luminosity versus mass relations from Fig. 3 can be contrasted against other simulations, providing a benchmark comparison. As with other mass–observable relations, they can in the future also be compared against observational data, when large non-Ly α-selected surveys are available.

To provide a first comparison, red markers show results from selected computational models at z = 3 (Faucher-Giguère et al. 2010; Goerdt et al. 2010; Rosdahl & Blaizot 2012). All data points include the luminosity for all bound halo gas, and thus should be compared against the black dashed lines. Given substantial physical model and simulation differences, including the details of stellar and AGN feedback, cooling, and on-the-fly self-shielding in TNG with respect to these previous simulations, differences are to be expected. We discuss this in more detail in Section 4.3.1. Nevertheless, our TNG50 results are roughly bracketed from below by the Faucher-Giguère et al. (2010)³ result and Goerdt et al. (2010) and Rosdahl & Blaizot (2012) from above. These studies focus on emission from collisional excitations and recombinations⁴ without an explicit treatment for star formation and with exception of Faucher-Giguère et al. (2010), without considering Ly α scattering.

In Fig. 4, we show four observable properties and relations for Ly α emitting galaxies. None were used for calibration, enabling us to use them to assess the performance and realism of the emission and rescaling model. In the first panel, we show the UV LF with (orange) and without (blue) accounting for dust, as treated in equation (8). The dust-free UV LF has a nearly constant slope down to magnitudes of M_UV = −22 after which a substantial decline sets in (not shown). With dust attenuation, this decline sharpens and already sets in at M_UV = −19. Red error bars show observational data points from the photometric redshift sample in Mehta et al. (2017) at z = 2. Generally, the simulated UV LF is in reasonable agreement with the observational data, except for the high luminosity end beyond M_UV = −20 where the simulated LF drops off faster than observed.

The second panel shows Ly α luminosity as a function of the UV magnitude. The median (orange line) indicates a positive correlation between UV and Ly α luminosity, which begins to flatten towards higher UV luminosities. This is a consequence of our dust model and empirical calibration, which increasingly suppress Ly α emission from massive, dust-rich galaxies.

This suppression can be seen clearly in the third panel, showing the Ly α escape fraction as a function of stellar mass. Here, we include only galaxies with L_Ly α > 3 × 10⁴¹ erg s⁻¹ to be roughly consistent with the data selection function. We calculate the escape fraction as the ratio between the scattered Ly α luminosity after rescaling and the intrinsic Ly α luminosity before rescaling ignoring contributions from excitations, which is closest to the methodology used to infer the Ly α escape fraction in observations. The Ly α escape fraction is high for low-mass galaxies, approaching unity below 10⁸ M_⊙, before rapidly dropping to only ∼10 per cent by 10⁹ M_⊙ and reaching ∼2 per cent for 10¹⁰ M_⊙. In red, we show observations for the escape fraction from Weiss et al. (2021) based on emission-line galaxies at z ∼ 2 (Bowman et al. 2019), broadly consistent with other observations (Hayes et al. 2010; Ciardullo et al. 2014; Sobral et al. 2017; Snapp-Kolas et al. 2022).

In Weiss et al. (2021), the escape fraction is estimated using the flux ratio of Ly α to dust-corrected Hβ, adopting |$f_{\mathrm{esc}}^{\mathrm{Ly}\alpha }=1/23 \frac{{\mathrm{Ly}\alpha }}{\mathrm{H}\beta }$|⁠. This relation is derived assuming all emission arises from recombinations and neglecting minor changes due to the temperature dependence of the ratio between the Ly α and Hα recombination coefficients. The reasonable agreement in comparison to observations suggests that our model primarily captures dust attenuation (as intended) rather than, e.g. modifying effective galaxy star formation rates. When including emission from excitations, we obtain a slightly steeper relation, which is still broadly consistent with observations.

The last panel shows the distribution of REWs, for all galaxies with L_Ly α > 3 × 10⁴¹ erg s⁻¹, to best match the selection in the observational data. The distribution is positively skewed and unimodal, peaking around 25 Å. Sixty-eight per cent of emitters have REW ≥20 Å, while the median value is 36 Å. The tail follows an exponential decay similar to the form found in observations by Oyarzún et al. (2017).

In all four panels, we account for only the most basic and zeroth-order selection effects. These comparisons will in the future benefit from more sophisticated mocks. Nevertheless, our overall result is that the calibrated Ly α emission model is reasonably realistic and broadly consistent with available z = 2 data.

3.2 The physical origin of Lyman-alpha emission

In Fig. 5, we show density-temperature phase space diagrams weighted by the total Ly α emissivity. Different rows represent the various emission mechanisms: nebular emission from stellar populations (top), collisional excitation (middle), and recombination (bottom). In the left-hand panels, we show the emissivities for the intrinsic photons, i.e. the gas state at the site of emission. In the right-hand panels, the gas state at last scattering is instead shown, i.e. the gas state at the location where the Ly α photon last interacts before escaping towards the observer. The contours show the brightest regions of phase space containing 99 per cent, 95 per cent, and 50 per cent of the total luminosity.

The first row shows the emissivities from nebular emission sourced by stellar populations. As expected, most of the emission originates in the star-forming gas above n_H ≥ 0.13 cm⁻³. However, nearly the entire phase space contains stellar emission due to a small fraction of stars that have migrated out of star-forming regions (left-hand panel). At last scattering (right-hand panel), a significant fraction of photons have been redistributed out of star-forming regions. We quantify this by integrating the luminosity in the cold-dense region of the phase diagram, defined as the area below the black dashed line in Fig. 5. We find that 99.3 per cent of the intrinsic emission stems from this region, dropping to 69.9 per cent after radiative transfer, indicating the large spatial redistribution of photons originating around stellar populations.

For diffuse emission from excitations and recombinations, emission primarily stems from a similar region of the phase diagram: dense, cold gas. Emission for recombinations is also significant from low-density regions characteristic of the IGM, where ionization is powered by the UVB radiation field. On the other hand, emission from collisional excitations quickly drops faster towards low densities. Comparing the two columns, we see a significant redistribution of photons in phase space when comparing intrinsic and scattered photons. After accounting for radiative transfer, photons more homogeneously sample larger regions of phase space. A significant fraction of the total Ly α emission for both recombinations and collisional excitation arises, after scattering, from the cool, low-density IGM.

In Fig. 6, we give an overview of the global Ly α luminosity budget at z = 2, considering the different emission mechanisms and how photons are redistributed by scattering between spatial components. In the bar charts on the right, we show the relative importance of emission mechanisms and spatial components reaching the observer. The relative fractions of intrinsic emission (‘I’) stemming from stars, excitations, and recombinations are shown in blue, red, and green, respectively, with corresponding total contributions of 29 per cent, 48 per cent, and 23 per cent, making collisions the most important emission mechanism for the overall volume. Note that if we exclusively consider the denser environments of filaments, with dark matter overdensities from 3 to 30, emission from stars starts to dominate the luminosity budget (not shown). Each coloured area is subdivided into five regions, corresponding to the five distinct spatial components: CGM, central, satellites, outer halo, and IGM, with increasing transparency. In addition, the two grey pie charts show the global fraction of each spatial component, regardless of emission mechanism, by mass and volume.

Focusing on observable Ly α photons (at last scattering; ‘LS’) while the IGM encompasses 99.8 per cent of the entire volume of the Universe and contains 83 per cent of all matter, only |${\sim }16~{{\ \rm per\ cent}}$| of emission originates here. On the other hand, the CGM is the single largest spatial component for all three emission mechanisms, contributing more than half of all emission, with only 10 per cent of matter in the Universe. Central galaxies contribute a similar fraction of emission due to stars but add little via recombinations and collisions, with a combined |${\sim }14~{{\ \rm per\ cent}}$| of all emission. The outskirts of haloes and satellites combined contribute |${\sim }18~{{\ \rm per\ cent}}$| to the overall emission, signifying sub-dominant but non-negligible spatial components.

On the left of Fig. 6, we show four matrices for the emission mechanisms: all combined (i.e. the total, in grey), nebular emission due to star formation (blue), excitation (red), and recombination (green). The matrices reveal the degree to which photons emitted in a certain component are redistributed to another due to radiative transfer effects. Each column is normalized to unity. As a result, the numbers give the fraction of Ly α luminosity emerging (i.e. last scattering) from a given spatial component, which originated (i.e. intrinsically) elsewhere. For satellites, centrals, and the CGM, most of the observed emission does originate from those components themselves. However, this does not hold for the outer halo and the IGM, where most of the emission comes from the CGM of galaxies (for recombinations and excitations) and central galaxies themselves (for stellar populations). Specifically, 94 per cent of all Ly α radiation reaching us from the IGM – that is, the large-scale cosmic web – actually does not originate there, stressing the importance of the Ly α radiative transfer modelling in the IGM.

The total luminosity density of Ly α, based on our fiducial emission and radiative transfer models applied to TNG50 at z = 2, is |$\dot{\rho }_{\mathrm{Ly}\alpha }=4.4\cdot 10^{40}$| erg s⁻¹ cMpc⁻³. Central galaxies emit a quarter of the total luminosity density with 1.1 × 10⁴⁰ erg s⁻¹ cMpc⁻³. Moreover, 98 per cent of all emission originates within haloes, with only 9.1 × 10³⁸ erg s⁻¹ cMpc⁻³ originating in the IGM. When considering the luminosity density at last scattering, this picture changes significantly: central galaxies retain only half of their emission upon reaching the observer with 5.8 × 10³⁹ erg s⁻¹ cMpc⁻³ and haloes only contain 73 per cent of the luminosity budget scattering there last. The remainder is redistributed into the IGM, boosting its luminosity density by an order of magnitude to 1.2 × 10⁴⁰ erg s⁻¹ cMpc⁻³.

3.3 Visual inspection of the large-scale topology, origin, and physics of the Lyman-alpha cosmic web

Fig. 7 visualizes a large 3 × 10 pMpc region selected to include both large-scale and small-scale filamentary gas structures (white rectangle in Fig. 1). From left to right, we show: the ionizing photon rate, temperature, line-of-sight velocity, neutral hydrogen column density for a slice depth of 5.7 Å, the dark matter density distribution, and the resulting Ly α SB map, based on our fiducial model, including radiative transfer effects.

The ionizing photon rate projection is sparsely populated, showing that the stars within galaxies trace only the highest overdensities. These sites of ionizing photon production are correlated with the filaments visible in neutral hydrogen column density, for example. In the temperature, line-of-sight velocity and neutral hydrogen density panels, we can clearly recognize filamentary structures across a variety of length scales from shorter ∼100 kpc filaments to larger ∼ Mpc long filaments. In general, filaments have a velocity relative to their background, are relatively cold and optically thick. They are often surrounded by hot regions powered by feedback and shocks. The dark matter overdensity field traces the neutral hydrogen column density. However, while the dark matter has a more pronounced and clumpy substructure, it does not show the distinct string-like filaments visible for neutral hydrogen. Most importantly, the Ly α SB clearly traces the filamentary structure, suggesting that the cosmic web in Ly α emission is closely linked to the underlying gas and dark matter distributions on large scales. However, Ly α does not have a simple one-to-one mapping with gas density, as visible for instance in the more diffuse and less pronounced edges and boundaries in comparison to the neutral hydrogen distribution.

In Fig. 8, we study the Ly α SB projection for this same region of space in more detail. We decompose the total emission into the contributions from the three emission mechanisms (different columns). The top panels show intrinsic emission, i.e. Ly α photons where they are initially emitted. The top left panel shows nebular emission sourced by star formation, where we see sparsely distributed emission stemming from massive star-forming galaxies, similar to the projection of ionizing photon rates in Fig. 7. There is no strict proportionality of ionizing photon rates and Ly α emission [see equation (3)] due to our application of the empirically calibrated rescaling model. Collisions and recombinations (top middle and top right panels) both roughly follow the distribution of neutral hydrogen and dark matter. For collisions, we see a higher contrast with strong emission in the filamentary structure that quickly drops below 10⁻²⁴ erg s⁻¹ cm⁻² arcsec⁻², compared to recombinations that maintain a higher SB, even in voids.

The bottom row of panels in Fig. 8 shows the SB maps after accounting for the radiative transfer, i.e. after processing these photons to account for the resonant scattering process. The most striking difference arises in star formation sourced emission. While intrinsically confined to within galaxies, after scattering, these same photons illuminate the filaments of the cosmic web. While photon redistribution for collisions and recombinations is less pronounced, the spatial diffusion of photons increases the SB of filament outskirts. In terms of peak SB and contrast, filaments are most pronounced due to Ly α emission from collisional excitation. In this case, radiative transfer effects lead to somewhat fuzzier filaments due to spatial diffusion.

All mechanisms independently give rise to SBes above 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻² around massive objects, but filamentary structures above 3 × 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻² are only visible in the proximity of massive galaxies and cosmic web nodes. In this case, AGN provide additional photoheating and -ionization. With decreasing SB, larger connected filamentary structures appear. Thin filaments extending more than a physical Megaparsec are easily identifiable by eye at ∼10⁻²² erg s⁻¹ cm⁻² arcsec⁻² due to collisional excitations and recombinations. We expand on the quantitative brightness of LAFs, and their observational detectability in the following sections.

In Fig. 9, we visually decompose the total Ly α SB of the same region of space into the dominant emission mechanism (left-hand panels), spatial component, i.e. origin (middle panels), and contributing mean halo mass (right-hand panels). In all cases, pixels of the images are coloured by the mechanism/component/halo mass that contributes most of the Ly α luminosity to that pixel. The left-hand panel of each pair shows the intrinsic emission, while the right-hand panel of each pair instead shows the result after scattering, i.e. the observable Ly α emission after accounting for the radiative transfer process. With respect to the emission mechanism in the intrinsic case (left-most panel), we find excitations within filaments and recombinations outside of them to dominate, and the separation is clear. There exist only small regions within filaments situated within the most massive haloes where star formation dominates. After radiative transfer, most filaments remain dominated by excitations, and the regions dominated by star formation sourced emission enlarge. Most importantly, emission in filament outskirts and in voids is now dominated by either recombinations or stellar emission rather than recombinations.

The middle panels of Fig. 9 show the dominant component, i.e. spatial origin, before and after radiative transfer. Before radiative transfer, we see that most of the area around the filaments is dominated by the CGM of galaxies, but there remain large unbound areas that are dominated by the IGM. Occasionally, regions within filaments occur where satellites, centrals, or outer halo origins dominate. After radiative transfer, most of the SB of cosmic web filaments is dominated by emission originating from the CGM of galaxies, residing both within and outside of filaments. This is a key result of our study. Only small regions remain where emission originating from the IGM dominates.

The fifth and sixth panels of Fig. 9 show the mean dark matter halo mass, from which photons dominate, weighted by Ly α luminosity. Photons that originate outside of all haloes are not considered in this case. In general, significant fractions of the filaments are filled by gas from massive haloes ≳10¹⁰ M_⊙, and the intrinsic emission from these relatively high-mass haloes is directly responsible for emission from the filament regions. The role of emission from haloes as a function of mass changes after we account for radiative transfer effects (right-most panel). In this case, we find that the emission from massive galaxies often overshadows that of the smaller haloes. This is facilitated by a significant photon flux of massive haloes reaching and scattering off the CGM and IGM surrounding smaller haloes and scattering into, and then out from, smaller haloes often hundreds to thousands of physical kiloparsecs away (Byrohl et al. 2021).

3.4 Surface brightness distributions

We now quantify the sky area covered by Ly α emission at different SB levels. In Fig. 10, we show the fraction of pixels with a certain SB for the respective emission mechanism before (i.e. intrinsic; dotted lines) and after (i.e. scattered; solid lines) radiative transfer. For collisions and recombinations (orange and green), the intrinsic SB distribution drops sharply at around 10^−17.5 erg s⁻¹ cm⁻² arcsec⁻², with collisions extending to SB values higher by roughly a factor of 3. If we further decomposing the distribution based on the existence of an AGN radiation field at a given location, we find the high SB tail for collisions and recombinations is dominated by emission from gas experiencing additional photoheating and -ionization (not shown). At lower SB, recombinations and collisions show a similar qualitative trend. For star formation sourced emission (blue), a substantial number of pixels in the intrinsic SB distribution extend to significantly higher values compared to collisions and recombinations.

Comparing dashed (for intrinsic photons) and solid (for scattered photons) lines, we find the largest impact of radiative transfer at the high luminosity end where the occurrence of peak SB values decreases by a factor of a few, and at SB values around 10⁻²¹ erg s⁻¹ cm⁻² arcsec⁻² where it is enhanced by an order of magnitude when radiative transfer is applied. That is, the importance of Ly α resonant scattering is actually maximal at the SB values, which correspond to the cosmic web filaments, making a full radiative transfer treatment essential. Regarding the three emission mechanisms, the impact of radiative transfer is largest for emission from stellar populations. Before scattering, the distribution is relatively flat, but after applying radiative transfer, large fractions of the previously Ly α dim sky are boosted in brightness, such that the SB distribution roughly resembles the distribution of the other two mechanisms.

In Fig. 11, we quantify the relative luminosity contributions for different emission mechanisms (top), components (i.e. spatial origin, middle), and contributing dark matter halo mass ranges (bottom) as a function of SB. The coloured area at a given SB reflects the relative fraction of each legend item, stacked vertically. The results are shown for the scattered photons, i.e. after accounting for radiative transfer effects.

The upper panel shows the contribution by emission mechanism. We find star formation (blue) sourced Ly α radiation dominates SB values above 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻². For smaller values, star-formation remains a relevant emission channel, contributing 10 per cent to 35 per cent. Between 10⁻²¹ and 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻², collisional excitation (orange) sources most of the observed SB. For values below 10⁻²¹ erg s⁻¹ cm⁻² arcsec⁻², recombination (blue) starts dominating.

The middle panel shows the fraction of Ly α luminosity, as a function of SB, originating from each spatial component. Luminosity originating from central galaxies (orange) dominates at high SB values, and down to 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻², coinciding with the trend of star formation as in the upper panel. Between 10⁻²⁰ and 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻², luminosity from the CGM (yellow) contributes more than 60 per cent. Below 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻², the contribution from the CGM declines, and IGM (black) contributions gradually grow until the IGM eventually becomes the most important origin below 10⁻²² erg s⁻¹ cm⁻² arcsec⁻². Satellites and outer haloes (magenta and purple) contribute |${\sim }20-30~{{\ \rm per\ cent}}$| to the budget throughout the full dynamic range of SB.

The lower panel shows the luminosity fraction contributed by emission from haloes as a function of halo mass. After we account for radiative transfer effects and consider Ly α photons as they would be actually observable, we find that intermediate-mass haloes are crucial Ly α sources. Specifically, between 10¹⁰ and 10¹¹ M_⊙, they contribute more than 50 per cent to Ly α filament luminosity, down to 10⁻²¹ erg s⁻¹ cm⁻² arcsec⁻².

We can evaluate these same contributions for intrinsic photons, i.e. neglecting radiative transfer (not explicitly shown). If we do so, this same halo mass range dominates at high SB above 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻², while the other three halo mass bins have roughly equal contributions below, and the IGM becomes the major contributor at the lowest SB levels, below 10⁻²¹ erg s⁻¹ cm⁻² arcsec⁻². Similarly, stellar contributions quickly approach zero, reaching <10 per cent below 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻², while the CGM contributes significantly less to the luminosity budget, and emission from the IGM dominates below 10⁻²⁰.

Finally, when including radiative transfer, we can also consider the spatial component at the point of last scattering (also not explicitly shown) rather than at the origin. When doing so, we find that there is a sharp transition at ∼10^−19.5 erg s⁻¹ cm⁻² arcsec⁻² below which nearly all photons are scattered by the IGM.

That is, we find that the Ly α cosmic web in emission is illuminated predominantly by photons that (i) originate within central galaxies and the CGM of intermediate-mass haloes and (ii) scatter into the IGM before reaching the observer.

In Fig. 12, we show the Ly α SB as a function of dark matter overdensity δ_DM. We construct the dark matter overdensity field and impose a Gaussian filter with σ = 100 pkpc. We then project the maximum overdensity in each pixel of the slice and create a 2D histogram together with the Ly α SB for the pixels in the slice. We show the median (luminosity weighted mean) in green (blue), separately for scattered (intrinsic/unscattered) photons as solid (dashed) lines. In the background, we show the probability density of logarithmic SB at given dark matter overdensity.

We find that the observed mean SB is a strictly monotonic function of overdensity starting at ∼10⁻²¹ erg s⁻¹ cm⁻² arcsec⁻² for overdensities of a few characteristic of the cosmic web filaments. This increases to 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² at overdensities around 200 characteristic of collapsed haloes. The mean SB is boosted by around an order of magnitude over the full overdensity range by redistributing Ly α emission from the densest regions. For scattered photons, the central 68 per cent of values around the median show a scatter of roughly 1.5 dex across the range of overdensities shown. In comparison, for intrinsic photons, this spread grows from 0.5 to 3 dex with overdensity. The tight correlation for intrinsic photons at low SB values is set by UVB photoionized hydrogen outside of filaments. This correlation broadens at larger values significantly due to the complex temperature and density structure within haloes and due to the lack of strong correlations with the smoothed dark matter field. For scattered photons, the scatter at low SB increases as a fraction of underdense regions are significantly boosted in SB when in proximity to Ly α bright haloes.

At high SB, the scatter decreases as photon contributions are smoothed out in hydrogen rich, overdense regions. There is a rapid increase between 10⁻²³ and 10⁻²¹ erg s⁻¹ cm⁻² arcsec⁻² where scattered contributions start dominating over otherwise lower SB values from intrinsic contributions. The SB mean (median) value grows from 10^−20.8 (10^−21.2) to 10^−19.2 (10^−19.9) erg s⁻¹ cm⁻² arcsec⁻² from an overdensity of 3–30.

3.5 Lyman-alpha filament identification and detectability

Next, we aim to identify and characterize LAFs. To do so, we adopt a relatively simple approach. We search for LAFs as connected regions above a certain SB threshold. We label this SB threshold SB_{fil, 10}, defined as log₁₀ of the SB in erg s⁻¹ cm⁻² arcsec⁻². The SB maps are evaluated on a fixed smoothing scale, for which we use a Gaussian filter with a full width at half-maximum (FWHM) of σ_{FWHM, fil}. Such a smoothing scale would also be imposed on data in order to identify large-scale structure rather than small-scale details by maximizing available the signal-to-noise ratio. We then measure the filament size L as the maximal distance between any points of a connected region with area A. From this region, we also determine the circularity c = 4A/(πL²) with a value of 0 indicating a 1D line and a value of 1 indicating a perfect circle.

Considering a region as elongated for circularity values below c < 0.5, we find that Ly α structures with lengths of L ≳ 400 pkpc are typically elongated. In the following, we thus commonly adopt a length minimum of 400 pkpc for LAFs. For reference, this length threshold roughly corresponds to an area of ∼0.25 arcmin for an elongated structure with c = 0.5 at z = 2.0. We use a fiducial value of σ_{FWHM, fil} = 3.5 arcsec. We find the number of filaments above >400 pkpc remains nearly constant for smoothing scales ≤10 arcsec, irrespective of the fixed SB threshold considered. Finally, we adopt a fiducial SB threshold of SB_fil = 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻² as a reasonable value in reach of future surveys. While this results in a robust set of identified filaments, we note that smaller structures with lengths less than 400 pkpc are more sensitive to the analysis choices.

In Fig. 13, we show Ly α SB maps of 12 selected zoom-in regions from Fig. 1. We select each region by hand to contain filamentary structures spanning a range of filament sizes and luminosities. We quantify the frequency of such structures below. We include contour levels at 10⁻²⁰, 10⁻¹⁹ and 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² in light yellow, orange, and red, respectively. At 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻², a range of large filamentary structures can be seen, while at higher SB thresholds structures resemble LAHs (mostly orange), i.e. more circular emission centred on haloes, or point sources, as they appear given the smoothing scale (mostly red). On the left side within each panel, we show information for the emission mechanism, spatial component, halo masses (pie charts and bar charts, colour coding consistent with other plots) as well as the length, area, and dark matter overdensity. These quantities are only evaluated for the diffuse part of filaments with SBes between 10⁻²⁰ and 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻² and are calculated for the largest filament overlapping the field of view.

While there are variations, particularly for the fraction of Ly α photons sourced by star-formation, we find the collisional excitations and the CGM to be the dominant emission mechanisms contributing to the luminosity of filaments. The bar charts indicate contributing halo masses at emission for intrinsic emission (‘I’, left), after radiative transfer by origin (‘P’, centre) and last scattering (‘LS’, right). We find that the majority of Ly α photons originate from haloes between 10¹⁰ and 10¹¹ M_⊙, except in cases where more massive emitters exist in the filament. While most of these photons originate within haloes, they scatter off the diffuse IGM before reaching the observer.

In Fig. 14, we quantify filaments by their linear extent (size), area, and SB. The upper panel shows a histogram (i.e. space volume density) of filaments as a function of size for different SB thresholds. Overall, small structures are much more common than larger ones. Similarly, dim structures are much more common than brighter ones. The abundance of filaments with lengths below ∼400 pkpc grows by a factor of a few when decreasing the SB threshold by a factor of 10. For larger filaments, however, the number count is independent of SB threshold, at and below our fiducial value of 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻². For 10⁻²² erg s⁻¹ cm⁻² arcsec⁻², the distribution drastically changes its behaviour as many of the larger structures merge. Finally, for smaller filament sizes from ∼100 to 1000 pkpc, the decrease in number counts roughly follows a power-law with a slope parameter α > −2.0 indicating larger filaments cover a comparable or larger sky fraction relative to smaller filaments.

The middle panel of Fig. 14 shows a scatter plot of average filament SB versus area. Each filament is represented by a single circle, where colour indicates its circularity parameter. For detected structures above 100 arcsec², the number of objects drops significantly above a few times 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻², while the largest structures occur around this value. The number of filaments quickly drops towards the detection threshold of 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻², implying that low SB filaments exist around brighter objects rather than on their own. With increasing area, the circularity decreases, demonstrating that the largest Ly α structures are elongated and filamentary in nature. On average, objects have increasingly circular shapes at larger SBes, indicative of halo-centred emission such as LAHs rather than intergalactic filaments.

To offer a face-value comparison, we also show the observational data points of the filament detections from Bacon et al. (2021), including only confident detections. We emphasize that the filament identification and measurement methods in that work differ from ours, and we do not intend a quantitative comparison at this stage. Each data point is plotted twice: once with the SB and area at the redshift of its actual detection, and again with a second point representing the same filament at redshift z = 2. To do so, we assume that the physical size and luminosity are both constant, which increases the SB while decreasing its angular area.

Qualitatively, the confident detections in Bacon et al. (2021) scaled to z = 2 represent the brightest filaments in our simulations. These observed filaments have typical sizes around ∼1000 arcsec² at the upper end of sizes that can be expected given the limited field of view. In the simulations, the majority of objects are smaller and therefore harder to detect. In addition, observational searches within known LAE overdensity regions will bias the filament size distribution towards higher values, as we find to be the case in our models.

The lower panel of Fig. 14 shows a scatter plot of filament luminosity versus area. While this plot is a simple mapping from the data shown in the middle panel, it offers complementary insights. First, we note that a relation between Ly α luminosity and area is unavoidable. Shaded regions show the allowed space for objects extracted for a given SB threshold and smoothing scale. The upper limit is set by the minimum area for a smoothed point source at the SB threshold. The lower limit is set by the minimum luminosity covered by an area above the SB threshold.

We show to shaded regions for LABs and LAFs, defined by different SB thresholds and smoothing scales. In each region, circles show detected Ly α structures at the respective SB threshold and smoothing scale. Just as in previous panel, we colour the filaments by their circularity. Additionally, we show circles with a black edge colour for LAB detections with ≥100 pkpc and for LAF detections with ≥400 pkpc.

The upper shaded region uses a SB threshold of (1.0 + 2.0)/(1.0 + 3.0)⁻⁴1.4 × 10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² and a Gaussian smoothing with FWHM of 1.4 arcsec as in the observational LAB data sets in Matsuda et al. (2004, 2011) at z = 3. In this case, extended Ly α structures in TNG50 have roughly similar characteristics, although we leave a detailed exploration of the abundance and properties of Ly α blobs (LABs) for future work.

The lower shaded region is based on the values of our fiducial definition for extended LAFs. We find a power-law scaling of slope 1.1 for the luminosity as a function of area for LAFs, which is in excess of the expected 1.0 by construction. The typical SB values of LAFs are roughly an order of magnitude above the SB threshold for filaments with large areas. The largest filaments in TNG50 at z = 2 have an area of ∼10⁷ pkpc² and luminosities of a few times |$10^{44}\, \rm {erg s^{-1}}$|⁠.

In Fig. 15, we show the mean SB distributions of detected filaments with L > 100 pkpc as a number of violin plots. We split by emission mechanism (top panel), spatial component (middle panel), and halo mass (bottom panel). The left half of each coloured regions shows the probability density function of the mean SB within the defining contour at a SB threshold of 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻². The right half of each coloured region shows the mean SB of the diffuse parts of the filaments, defined as contributions at SB values below 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻². The circles indicate the respective median values. In the panel for the spatial origin, we also show the median for the contributions at last scattering as squares.

Overall, we find that emission from stars and collisions contribute equally, with a substantial contribution from recombinations. For the diffuse regions (right-hand side of violins), star formation plays a subdominant role and collisions continue to dominate. We note that there is a large filament-to-filament variation in the relative contribution from stars, depending on the existence of nearby dust-poor massive galaxies. Emission from the CGM and central galaxies dominate in filaments (middle panel). Even in diffuse low SB regions, the emission originating in the CGM remains the dominant contributor. Finally, halo masses between 10¹⁰ and 10¹¹ are the most important in terms of contributing to filaments as well as in their more diffuse regions (bottom panel).

While not explicitly shown, we also find that the relative importance of different contributions does not significantly depend on the filament area or length. In addition, if we decompose the total Ly α luminosity of detected filaments rather than SB, we obtain consistent results with those described above. Namely, for both small and large filaments – independent of angular area – these structures are dominated by Ly α photons sourced by collisional excitations (the other two mechanisms not far behind), which originate predominantly (>50 per cent) in the CGM of intermediate-mass (⁠|$11 \lt \log {(M_{\rm halo}/\rm {M}_\odot)} \lt 12$|⁠) haloes before scattering into filaments and towards the observer.

In Fig. 16, we show the dark matter overdensity traced by LAFs, as a function of their area. In blue, we consider the entirety of each filament, while in orange, we restrict to the diffuse, i.e. low SB <10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻² regions of filaments, thereby avoiding bright, halo-centric emission, which may be embedded within them. We also show intrinsic (dashed) and scattered (solid) photons separately to assess radiative transfer effects.

Our main finding is that LAFs trace out increasingly overdense regions as their sizes increase. This relation increases monotonically, with overdensity rising from ∼20 to ∼150 as LAF area increases from 10² to |$10^{5}\,$|arcsec². In contrast, the overdensity traced by the diffuse component is comparably constant between ∼20 and ∼50. For filaments in their entirety, the corresponding overdensity is higher for intrinsic emission, as emission from the most massive objects escapes. For the diffuse component, the overdensity traced by intrinsic photons is lower by up to a factor of 2, as photons from compact sources scatter into more diffuse regions.

Diffuse LAFs across a range of size scales are an excellent tracer of the underlying matter of the cosmic web.

4 DISCUSSION

4.1 The Lyman-alpha cosmic web and its physical origin

Our synthesized Ly α SB maps enable us to assess the origin and illumination of Ly α emitting gas across a large cosmological volume at z = 2. Our emission model self-consistently captures all major scenarios invoked to explain extended Ly α emission: scattered light from central sources, gravitational cooling, satellite galaxies, fluorescence from AGN, and the recent suggestion of contributions from unresolved LAEs (Bacon et al. 2021).

Measuring isophotal areas above a fiducial SB threshold |$\rm {SB} \gt 10^{-20}$| erg s⁻¹ cm⁻² arcsec⁻² and an imposed Gaussian PSF with a FWHM of 3.5 arcsec, we identify a range of Ly α emitting structures. At these fiducial values, large structures with L ≥ 400 pkpc, which we call LAFs, are increasingly more elongated and less round in shape. Our careful tracking of Ly α photons in terms of their emission mechanism, spatial origin, and originating dark matter halo mass lets us evaluate which scenarios contribute to the emergence of a ‘Ly α cosmic web’.

4.1.1 What powers Lyman-alpha filaments?

On average, LAFs source roughly half of their global emission from collisional excitation, while 30 per cent is sourced by radiation from stars, leaving 20 per cent from recombinations. Our results suggest that emission originating from cold gas outside of haloes, i.e. from cold filaments of the IGM, does not significantly contribute to observable LAFs. Instead, most of the LAF emission arises from the CGM of haloes, together with central galaxies and satellites. Even when considering the low SB component (<10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻²) of our filaments, the vast majority of emission stems from within gravitationally collapsed haloes. Satellites provide between 10 and 20 per cent of the luminosity and are a visible though minor contributor. Fluorescence from AGN, which is captured (in a simplified manner) within our model due to the inclusion of AGN radiative effects in the underlying TNG model, is likewise only a minor player. Instead, the CGM surrounding central galaxies is the most important spatial origin of Ly α photons for most filaments. In these denser environments, emission from collisions in the cold gas dominates. These photons, however, do not reach us directly – 50–80 per cent of this emission reaches us only after scattering in the IGM, i.e. in filaments, causing them to shine in Ly α.

We can quantify the importance of these radiative transfer effects by defining a ‘boost factor’, as the ratio between scattered versus intrinsic photon luminosity, for a given spatial component (not explicitly shown). If we do so, the IGM has a boost factor of 1 dex, and the CGM itself is boosted below 10¹⁰ M_⊙ by a factor of a few due to contributions from higher mass objects (see Byrohl et al. 2021). On the other hand, satellites and centrals have a net loss of luminosity with boost factors less than unity (∼0.5–0.7 on average), regardless of halo mass, representing intrinsic emission that escapes into the CGM and beyond. We now explore this scenario in more detail.

Fig. 15 shows that haloes with total mass between 10¹⁰ and 10¹¹ M_⊙ source most filament luminosity. Contributions from fainter emitters with L < 10⁴¹ erg s⁻¹, even in filament regions of low SB, only play a minor role. This relative unimportance of low luminosity sources is due to three reasons.

First, our emission model has a steep L_Ly α(M_halo) relation at the faint end, which reduces the slope of the LF (where no observational constraints are available). Secondly, regardless of emission model, low luminosity emitters simply cannot contribute significantly due to the low clustering of faint galaxies within LAFs. To quantify this, we can model the SB stemming from the LF only, i.e. ignoring environmental effects, as

Here, d_L and d_A are the luminosity and angular diameter distance, Δz is the slice depth in terms of physical length, and δ_fil is the overdensity of emitters in a given filament compared to their cosmic mean. Using the Schechter form with parameters Φ_⋆ = 6.32 × 10⁻⁴ cMpc⁻³, L_⋆ = 5.29 × 10⁴² erg s⁻¹, and α = −1.8 from Konno et al. (2016), faint emitters below 10^41.75 erg s⁻¹ contribute up to 68 per cent of the total luminosity budget.

However, this assumes that the clustering, and thus the overdensity δ_fil, is constant across the luminosity range. This is not correct, as the clustering of haloes is strongly mass-dependent. LAFs most commonly occur around more massive objects, and we find a median overdensity of 30.7 for dark matter haloes with total mass 10¹¹ M_⊙ ≤M_halo ≤ 10¹² M_⊙ in filaments with lengths above 400 pkpc. The overdensity drops to 8.5 and 3.7 for mass ranges lower by 1 and 2 dex, respectively. Incorporating clustering, we find that haloes with luminosities below 10⁴¹ erg s⁻¹ (10⁴⁰ erg s⁻¹, 10³⁹ erg s⁻¹) contribute at most 22 per cent (7 per cent, 3 per cent), rather than 49 per cent (31 per cent, 19 per cent) when assuming equal clustering for all masses. Simply put, the abundance and distribution of low-mass haloes (i.e. LAEs) prevents them from playing an important role in large-scale LAFs.⁵

Thirdly and finally, the LF captures only a fraction of the total Ly α emission in the Universe. While emission from the IGM itself is indeed negligible (see Fig. 6), emission from outside of observationally accessible apertures, i.e. below observational SB limits, together with scattering from within these same apertures, out into the surroundings, contributes most of the global Ly α emission.

We turn to the recent observational detection of LAFs of Bacon et al. (2021). That work identifies five filamentary Ly α structures at z ∼ 3 and z ∼ 4.5 at high significance. The authors conclude that most of the emission cannot be accounted for by the UV background, nor by detected LAEs. Instead, they infer that very faint and unresolved LAEs, down to 10³⁸–10⁴⁰ erg s⁻¹, are the major source of photons in the observed LAFs.

A full comparison with the findings of Bacon et al. (2021) is not feasible due to differences in methodology as well as redshift. Principally, this is due to the complexity of the filament identification algorithm in comparison to our more simple SB threshold approach. With that caveat, our modelling does not appear to require any faint LAE population, i.e. emitters which are not simply already present in TNG50, to explain the existence and abundance of LAFs. In contrast to their interpretation, our work instead suggests that rescattered contributions from resolved haloes, particularly in the mass range of 10¹⁰ and 10¹¹ M_⊙, above L = 10⁴⁰ erg s⁻¹ produce significant numbers of LAFs. However, we do agree with the conclusion that most of the filament emission does not stem from unbound IGM gas but instead from the CGM and, to a lesser degree, other halo-bound gas such as the satellites and from central galaxies themselves.

4.1.2 How can we constrain different contributions in the Lyman-alpha cosmic web?

Our analysis is based on the combination of the underlying TNG50 hydrodynamical simulation, the Ly α emission model, and the radiative transfer process. To disentangle these potentially degenerate effects, e.g. photon emission mechanisms, contributions by halo mass, the impact of radiative transfer, and the role of supernova and SMBH feedback, we need observational measures that are separably sensitive to each of these aspects.

We have undertaken a number of variation runs, including switching each individual emission mechanism on and off, and turning the radiative transfer on and off. In each case, we run the entire calibration procedure from scratch in order to match the observed LAE LF at z = 2. Overall, we find that these experiments all produce a similar, and thus robust, space number density of L > 400 pkpc filaments. This suggests that our primary predictions, related to the properties and abundance of LAFs, are fairly robust. As a downside, however, degeneracies seem to be present, which preclude the ability to clearly differentiate between many of the underlying physical processes. In the future, we plan to therefore explore:

• Alternative filament identification and characterization methods, such as multiscale filtering and two-point statistics, or dedicated filament detection algorithms such as disperse (Sousbie 2011).

• Correlations with complementary Ly α data in filaments, particularly embedded LAE populations, including emitter luminosities, radial profiles, equivalent widths, and spectra.

• Correlations with complementary galaxy properties, particularly regarding AGN activity, stellar populations, and dust content.

As a first step, we study the relation between LAEs and LAFs in Fig. 17 for our fiducial model. In the top panel, we show the fraction of filaments with a LAE above a given luminosity L, denoting |$\log _{10}(\mathrm{L}/(\mathrm{erg}\, \mathrm{s}^{-1}))$|⁠. Ly α luminosities are calculated as in our fiducial model with an aperture radius of 1.5 arcsec, see Section 2. We find that all filaments have at least one LAE above 10³⁹ erg s⁻¹, and the vast majority also have at least one LAE above L > 10⁴⁰ erg s⁻¹. However, at a luminosity threshold of 10⁴¹ erg s⁻¹ (10⁴² erg s⁻¹), only half of filaments with length above 100 pkpc (300 pkpc) host such an LAE. With increasing length, filaments are more and more likely to contain such a bright emitter, while smaller filaments rarely do so.

In the middle and bottom panels, we investigate the degree to which filaments are actually made up of emitters. The middle panel shows the fraction of filament area that remains after masking all LAEs above a given luminosity threshold, while the bottom panel similarly shows the fraction of filament luminosity. Solid lines give the median for filaments with at least one emitter above the luminosity threshold, while dashed lines show the median for filaments without any such emitter. Shaded regions show the associated central 68 per cent scatter.⁶ Considering bright emitters of L > 10⁴² erg s⁻¹, we find that they account for the bulk of the Ly α luminosity for smaller filaments, and up to 40 per cent for the largest. In terms of area, these emitters are responsible for up to half of the area for smaller filaments, but only 15 per cent for the largest. If we instead consider LAEs down to a luminosity threshold of 10⁴¹ erg s⁻¹, they represent a similar ∼50 per cent of the area of small filaments, and ∼30 per cent of the area of the largest filaments, as well as the majority (≳80 per cent) of filament luminosity.

Going to even lower LAE luminosities does not appreciably change these results. That is, fainter emitters with L < 10⁴⁰ erg s⁻¹ do not contribute to LAFs in terms of either area or luminosity. For large filaments, more than half of the area and 15 per cent of the luminosity remains unassociated with emitters L > 10³⁹ erg s⁻¹. This is consistent with the true fraction of 15 per cent of Ly α luminosity within filaments reaching us from unbound gas. Further, only 10 per cent of this initial 15 per cent luminosity fraction originates within the IGM. As previously found, the SB level of the diffuse IGM within filaments is thus boosted by an order of magnitude through scattered photons from nearby haloes compared to its intrinsic emission.

In addition to our fiducial virial radii mask, we have also considered a fixed 3 arcsec aperture for all emitters (not shown), which is more accessible observationally. We find that more than 90 per cent of filament area and 60 per cent of filament luminosity remain after masking above a luminosity threshold of 10⁴¹ erg s⁻¹. Decreasing the emitter luminosity threshold removes an increasing fraction of area and luminosity, as opposed to the plateau effect seen for the virial radii apertures. This is, however, driven by the substantial amount of IGM masked. For L > 10³⁹ erg s⁻¹, the luminosity and area fractions remain above those of the virial radius aperture masking.

4.1.3 Connection with Lyman-alpha halo radial profiles

While the detection of diffuse LAFs remains challenging, observations can increasingly identify brighter emission centred on galaxies in the form of LAHs and LABs (Steidel et al. 2011; Borisova et al. 2016; Momose et al. 2016; Leclercq et al. 2017; Cai et al. 2019). Their stacked radial profiles are observed to flatten at large radii in HETDEX survey data (Gebhardt et al. 2021), far beyond the halo boundary (Lujan Niemeyer et al. 2022a), a phenomenon identified in our previous analysis of TNG50 (Byrohl et al. 2021). We now connect these halo-centric profiles to our current cosmic web study.

Fig. 18 shows the median stacked radial profile of REW>20 Å LAEs with luminosities between 10⁴² and 10⁴³ erg s⁻¹. We show the contributions split by the spatial component they last scatter from, i.e. observable Ly α photons (left) and intrinsically originate from (right). We see that the large distance flattening primarily arises due to contributions last scattered by the IGM. However, these photons mainly originate in the CGM, with smaller contributions by central galaxies and centrals. This boosts the IGM luminosity by roughly an order of magnitude over its intrinsic emission, consistent with our findings for the diffuse parts of LAFs.

The original source of these photons in the outskirts of LAH profiles, as shown in the bottom panels, are intermediate mass haloes with 10¹⁰ M_⊙ ≤M_halo ≤ 10¹¹ M_⊙. Even more massive haloes start to non-negligibly contribute at large, ≳100s of kpc, distances (see also Byrohl et al. 2021). Our current work demonstrates that at these SB values, LAH radial profiles in fact trace filaments of the cosmic web itself. In particular, they begin to represent contributions from the fuzzy outskirts of filaments illuminated by scattered photons.

4.1.4 Lyman-alpha luminosity global budget

The luminosity density of Ly α photons gives us a handle on the global emissivity across a large-scale, cosmological volume. The LF of Fig. 2 has a total luminosity density of |$\dot{\rho }_{{\mathrm{Ly}\alpha }}=4.3 \cdot 10^{39}$| erg s⁻¹ cMpc⁻³ when integrated across the observed luminosity range L > 10^41.75 erg s⁻¹. By construction, this is close to the integrated Schechter form from Konno et al. (2016) within |${\sim }10~{{\ \rm per\ cent}}$|⁠. However, this is an order of magnitude lower than the true global luminosity density of |$\dot{\rho }_{{\mathrm{Ly}\alpha }}=4.4 \cdot 10^{40}$| erg s⁻¹ cMpc⁻³. For Ly α emission included within the LF, the majority of emission stems from faint sources. In particular, LAEs below 10^41.75 erg s⁻¹ contribute |$\dot{\rho }_{{\mathrm{Ly}\alpha }}=6.1 \cdot 10^{39}$| erg s⁻¹ cMpc⁻³. The majority, 5.5 × 10³⁹ erg s⁻¹ cMpc⁻³, of this stems from emitters above 10⁴⁰ erg s⁻¹. Furthermore, the observational Ly α LF is constructed by including only emitters with an equivalent width cut of REW>20 Å. Our fiducial model shows that this REW threshold removes LAEs with a total luminosity density of |$\dot{\rho }_{{\mathrm{Ly}\alpha }}=2.8 \cdot 10^{39}$| erg s⁻¹ cMpc⁻³. If we instead consider all LAEs, keeping fixed the 1.5 arcsec aperture, they account for |${\sim }30~{{\ \rm per\ cent}}$| of the global Ly α luminosity.

Of the ∼70 per cent of the global Ly α budget not captured within this aperture⁷, 60 per cent reach us from bound halo gas outside of the aperture and 40 per cent from the IGM. Again, we stress the importance of Ly α radiative transfer for this result. Intrinsically, without radiative transfer, only 2 per cent of the global budget originates in the IGM, but this value is boosted to 27 per cent after scattering. Similarly, intrinsic emission from within the aperture would account for |${\sim }50~{{\ \rm per\ cent}}$| of the total budget, but a substantial fraction of these photons are redistributed outside the aperture after radiative transfer.

In agreement with this finding, Ly α intensity mapping results at z ∼ 2 have found a significantly larger Ly α photon budget than accounted for by detected LAEs (Croft et al. 2016, 2018). However, we caution that since these Ly α intensity mapping experiments are constructed using quasars for their cross-correlation analysis, their large Ly α luminosity density estimate might result from the biased quasar environment and proximity effects. None the less, a recent analysis by Lin, Zheng & Cai (2022) yields |$\dot{\rho }_{{\mathrm{Ly}\alpha }}=6.6^{+3.3}_{-3.1} \cdot 10^{40}$| erg s⁻¹ cMpc⁻³, and the authors argue that Ly α emission from star-forming galaxies not associated with resolved LAEs could be responsible for such a high luminosity density. This observational estimate of the luminosity density is 50 per cent higher than our fiducial result but consistent to better than 1σ.

LAFs with linear extents above 100 pkpc (400 pkpc, 1000 pkpc) contain 76 per cent (59 per cent, 47 per cent) of the global Ly α photon budget, a substantial amount compared with the Ly α luminosity captured by the LAE LF.

With respect to the brightest luminosity Ly α sources, the observed LF above |$L \gt 10^{43} \rm {erg s^{-1}}$| flattens considerably, and these rare, bright emitters are directly related to AGN (Konno et al. 2016). Our model does not include the intrinsic Ly α emission from AGN, which are regardless quite rare in the relatively small volume of TNG50. As a result, we underestimate the abundance of LAEs above this threshold. However, only a small fraction of the total luminosity density arises from these sources.

4.2 Prospects of observing the Lyman-alpha cosmic web

Our modelling provides quantitative predictions for the expected occurrence of observable LAFs. For filaments of a linear size L > 400 pkpc, we predict a moderate density of 2 × 10⁻³ cMpc⁻³ above an SB threshold of 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻² for a line-of-sight slice depth of 5.7 Å and a smoothing with a FWHM of 3.5 arcsec (see Fig. 14).

MUSE has a 1 arcmin² footprint. This covers a volume of 1800 cMpc³ (3600 cMpc³, 5300 cMpc³) between z = 2.0 and 2.5 (3.0, 3.5). Neglecting SB dimming, redshift evolution, and noise complexities beyond an SB threshold, this would imply an average of 2.8 (5.5, 8.2) filaments above an SB threshold of 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻² for Ly α structures with L > 400 pkpc. At a higher SB threshold of 3 × 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻² (10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻²), this number reduces to 1.4 (0.4) between z = 2.0 and 2.5.

Three other relevant IFU instruments are KCWI, BlueMUSE, and HET-VIRUS, with field-of-view footprint sizes 0.2, 2.0, and 54.2 times the footprint of VLT-MUSE, respectively.⁸ The expected detection counts can be scaled linearly for a first estimate. Clearly, statistically significant samples of ≳10 filaments can only be achieved for large survey volumes. With current instruments, this will require many pointings and mosaicing, with the exception of HET-VIRUS.

We further quantify the observability of LAFs in Fig. 19, which shows the expected number of detected LAFs with L ≥ 400 pkpc for a hypothetical survey, as a function of survey sky footprint (linear length, assumed to be square) and an SB threshold. We again adopt a fiducial redshift slice depth of 5.7 Å. Lines show contours for detection probability in a single redshift slice, between 1 per cent and 80 per cent. We consider two cases: random pointings (left-hand panel) and targeted pointings centring on a previously known LAE with L > 10⁴² erg s⁻¹ (right-hand panel).

For reference, a single pointing of VLT-MUSE corresponds to 60 arcsec, which is the left edge of this figure, while a single pointing of HET-VIRUS corresponds to 951 arcsec (neglecting the non-unity filling factor), at the right edge of this figure. Regarding the sensitivity, for example, the 10-h MUSE UDF data, with a 3 × 3 tiling (∼180 arcsec total linear extent), reaches a 1σ SB sensitivity of |$5.5 \cdot 10^{-20} \, \rm {erg\, s^{-1}\, cm^{-2}\, arcsec^{-2}}$| Å⁻¹ in a 1 arcsec aperture (Leclercq et al. 2017).

For large footprints, the probabilities for random and targeted pointing strategies at z = 2 are similar.⁹ We predict a 50 per cent detection probability at 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻² for a 1000-arcsec survey footprint. As these are achievable SB levels, our modelling suggests that a sufficiently large survey at such a sensitivity will easily detect filaments of the Ly α cosmic web. The situation is more difficult for small survey footprints. At 60 arcsec and the same SB threshold of 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻², the targeted pointing has only a 10 per cent detection probability, while the random pointing has only a 1 per cent chance. To achieve the same probability, the random pointing would need to reach a SB threshold roughly an order of magnitude lower, ∼2 × 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻².

However, one promising avenue for the observational detection of diffuse LAFs is with integral field spectrograph/unit (IFS/IFU) instruments. In this case, a significant wavelength coverage along the spectral dimension implies that the observed volume along the line-of-sight direction is significantly larger than our fiducial, narrow slice width of only 5.7 Å. A rough estimate for the detection of filamentary structures in an extended survey volume along the line of sight can be estimated from the number densities given at the beginning of this section and Fig. 14.

Our use of a single SB threshold number to characterize the sensitivity of an observation is only meant to provide intuition. The real power of the modelLing approach presented here is that, for any proposed observational strategy, detection probabilities can be quantitatively assessed. Specific targeting and/or stacking concepts can be evaluated a priori (for example, oriented pair experiments; Gallego et al. 2018). The ideal observational experiment, given instrumental and practical constraints, can be designed to detect the cosmic web in Ly α emission. We can also test new data analysis/statistical methods to better identify and characterize such structures in difficult regimes. For example, in low signal-to-noise data, in the presence of substantial foreground contamination, and/or when extended low SB emission could be lost due to overaggressive background subtraction on a limited field of view.

4.3 Modelling uncertainties and future directions

The emission of Ly α photons and their radiative transfer are sensitive to the physical state of gas, from galactic, to circumgalactic, to intergalactic scales. We therefore discuss current limitations of our Ly α emission model, the Ly α radiative transfer method, and the underlying hydrodynamical simulation.

4.3.1 Lyman-alpha emission modelling

Building a physically motivated and realistic emission model for Ly α on top of a hydrodynamical simulation result is challenging. Various Ly α emission mechanisms have been considered in the context of numerical studies (Kollmeier et al. 2010; Cen & Zheng 2013; Lake et al. 2015; Smith et al. 2019; Byrohl et al. 2021; Mitchell et al. 2021), although consensus on which mechanisms dominate for different observed classes of Ly α emitting objects has not been reached in the past. Cold dense gas accounts for the majority of the Ly α emission in our model, before scattering outwards into lower density surroundings. This implies that modelling the emission from overdense regions is always critically important. Because of radiative transfer effects, this is true even if one is primarily interested in lower density regions, such as cosmic web filaments, as we are here.

Characterizing the mean diffuse emission for collisions and recombinations within the virial radius as a function of halo mass (Fig. 3, calibrated model), we find a steep relation with a power law slope α > −2 below 10¹⁰ M_⊙. There is a mild discontinuity around 10¹¹ M_⊙ where haloes with radiative feedback from AGN start dominating the sample. Compared to other models, our fiducial, calibrated model shows a similar mass dependence as the lower-bound model by Faucher-Giguère et al. (2010). Our L_Lyα(M_halo) relation predicts moderate luminosities between the conservative model of Faucher-Giguère et al. (2010) and brighter alternatives, which have been previously applied to cosmological simulations (Yang et al. 2006; Goerdt et al. 2010; Rosdahl & Blaizot 2012).¹⁰

Ionizing radiation from stellar populations dominates Ly α emission in star-forming regions. In our previous work (Byrohl et al. 2021), we adopted a simple model where Ly α emission was proportional to the local star-formation rate of gas. Instead, the new model presented in this work is fundamentally different: we model emission for stellar populations according to their age and metallicity of stars, and apply the calibration described in Section 2.4. The empirically calibrated model accounts for dust, including the destruction of Ly α photons by dust, which would otherwise not be captured. This small-scale physics constitute a sub-resolution model, as it occurs below the resolution scale of the simulation. The previous lack of a dust model led to an overestimate for emission from star-forming regions in massive galaxies. Now, by calibrating the stellar luminosities against observed LFs, we substantially limit the stellar Ly α photon budget. Fig. 4 indicates that our rescaled model performs decently against other observables, and sufficiently well for our purposes. Future advancements could incorporate additional observables into the fitting procedure. This would help to break remaining degeneracies, and improve the overall physical fidelity of the emission model.

One such observable is the observed extent of LAHs. For example, Leclercq et al. (2017) find an exponential scale length r₀ of 4.4 ± 1.5 pkpc at z ∼ 3. With our fiducial model, for LAHs with 41.5 ≤ log (L/L_⋆) ≤ 42.5 and REW ≥ 20 Å, we find an exponential scale length of |$r_{0}=7.0_{-1.9}^{+3.7}$| pkpc (see Fig. 18), which is ∼60 per cent larger than the observational profiles. However, as previously demonstrated in Byrohl et al. (2021), the exponential scale lengths in TNG50 LAHs match observations when dust attenuation is negligible. At the same time, Figs 2 and 4 make clear that significant dust attenuation is needed to match the observed Ly α LF at z = 2. It is possible that a more sophisticated dust model could reconcile these differences, i.e. simultaneously matching the exponential scale lengths and the LF. At present, this tension tentatively suggests that our model produces too much diffuse emission from collisions and recombinations and has too little escaping Ly α luminosity from dense star-forming regions.

In the current approach, dust impacts only the stellar population sourced component and does not modify any of the diffuse emission mechanisms. The physical motivation is that optical depths and dust destruction is substantially smaller away from the ISM. In our best-fitting model, diffuse emission (i.e. recombinations and collisions) dominates the luminosity budget for haloes with observed luminosities below 10⁴⁰ erg s⁻¹ as star formation ceases in the TNG50 simulation. However, there is an additional peak around 10⁴² erg s⁻¹ where diffuse emission dominates over stellar emission, potentially hinting at an overestimate in the diffuse emission.

As our dust model is based on the rescaling of intrinsic luminosities, it will naturally overestimate spectral as well as spatial redistribution, when compared to runs self-consistently incorporating dust (see e.g. Laursen et al. 2009). To explore this scenario, we test an alternative model that explicitly includes a dust description for diffuse emission. This leads to a reduction of 15 per cent and 18 per cent for recombinations and collisions, respectively, which would not change our overall conclusions.¹¹ Even if star-forming regions contributed more than our fiducial model suggests, making them the dominant emission mechanism, our finding that compact spatial components and star-forming galaxies illuminate LAFs through resonant scatterings would remain robust.

To bracket the uncertainties in the modelling of collisional excitations, we ran a simulation without any collisions. In this case, we can still find a rescaling model to match the observed LF in agreement with the Ly α escape fractions and equivalent width distribution. We find that the majority of the global luminosity budget due to collisions in the fiducial model (2.1 × 10⁴⁰ erg s⁻¹ cMpc⁻³) is countered by an increased Ly α emission around stars in this variation run (1.3 × 10⁴⁰ erg s⁻¹ cMpc⁻³ → 2.4 × 10⁴⁰ erg s⁻¹ cMpc⁻³). This leads to an increased fraction of photons scattered into filaments originating from galaxies. Interestingly, the CGM still contributes 34 per cent of the global luminosity budget, primarily through recombinations but also stellar populations. We expect similarly robust conclusions when ignoring recombinations rather than collisions. The filament density for L > 400 pkpc decreases by ∼20 per cent (40 per cent, 80 per cent) at an SB threshold of 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻² (3 × 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻², 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻²) compared to our fiducial model.

Using a multiscale approach, such as the starlet wavelength transform (Starck, Murtagh & Bijaoui 1998) used in Bacon et al. (2021), could offer complementary information to differentiate underlying contributions. In the future, purely Ly α observation-based statistics, such as filament number densities, shapes, and presence of luminous LAEs, in addition to the UV luminosities of galaxies within filaments, can be explored to differentiate these scenarios (see Figs 14 and 17).

4.3.2 Lyman-alpha radiative transfer

LAFs have recently been studied with hydrodynamic simulations in two past works, which are particularly relevant to our study (Elias et al. 2020; Witstok et al. 2021). Neither of these studies incorporates radiative transfer or scattering, i.e. their results would correspond to our intrinsic emission. However, our analysis shows that the SB of the Ly α cosmic web is considerably boosted by radiative transfer. This finding directly invalidates the assumption of previous work, hypothesizing that scattering does not substantially affect, or even decreases, Ly α filament detectability.

For example, a consistency check of an individual filament in Elias et al. (2020) shows that Ly α radiative transfer redistributes emission to lower SB, boosting the occurrence of spaxels with SB values below 10⁻²¹ erg s⁻¹ cm⁻² arcsec⁻² (see their figure A1). Such low SB values are observationally unreachable, leading to the idea that radiative transfer decreases detectability. Their model, however, only includes emission from gas outside of haloes. When incorporating emission from within haloes, which contributes 98 per cent of our global emission, we find the occurrence of SB values below 10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻², typical for LAFs, to be boosted by radiative transfer (see our Figs 10 and 12). Similarly, Witstok et al. (2021) argue that the impact of Ly α radiative transfer is negligible for the observability of LAFs; however, this assumption is not actually tested with radiative transfer simulations. Based on the lack of a Ly α radiative transfer treatment and the lack of emission from within haloes, we consider the expectations from Elias et al. (2020) and Witstok et al. (2021) to be lower limits for the purposes of Ly α filament detectability.

In particular, Witstok et al. (2021) mitigate modelling uncertainties at high densities by imposing a overdensity threshold above which Ly α luminosity contributions are ignored. This limits the SB from recombinations to the ‘mirror limit’ at which an optically thick cloud reflects most of the UVB as Ly α emission. This would correspond to a value of ∼2.5 × 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻² for the UVB of TNG50 at z = 2.0. However, this assumption yields a very conservative lower limit for the recombinational (as well as collisional) Ly α emission in filaments, as it neglects any recombinations that can still be sourced in collisional equilibrium, as well as local sources and subsequent enhancements through Ly α scatterings. In our model, we do not adopt such an assumption.

On the other hand, Elias et al. (2020) find that LAFs are, in general, brighter and more detectable in simulations that incorporate the TNG galaxy formation model, in comparison to simpler ‘no feedback’ simulations. The authors suggest that the ejection of gas from lower mass galaxies, due to stellar feedback-driven winds, can distribute more gas into filaments. They then speculate that the lack of detection of LAFs to date suggests that this gas redistribution may be too strong in the TNG model. However, given the caveats discussed above, namely the importance of radiative transfer together with the lack of quantitative and/or statistical comparisons with observational data in Elias et al. (2020), we do not find this suggestion compelling.

In fact, we conclude that our TNG50+RT predictions for Ly α filament abundance and detectability are roughly consistent with the recent available observations of Bacon et al. (2021). While we cannot directly compare results given the different filament identification method, our model suggests a low single digit count of extended LAFs similar to the largest structure in (‘group2’ of Bacon et al. 2021). The observed structure has an linear extent above 450 pkpc and a circularity below 0.2, consistent with our expectations. Our simulated sample contains numerous smaller filaments, which are however not identified in Bacon et al. (2021), by construction, as they focus on filaments around LAE overdensities.

Also of interest is thatLABs represent another class of extended Ly α structures beside LAHs and LAFs. The ∼100 pkpc sizes of LABs fall between those of LAHs and LAFs, but their SB threshold is the highest at ∼10⁻¹⁸ erg s⁻¹ cm⁻² arcsec⁻² (Matsuda et al. 2004, 2011). LAHs, LABs and LAFs are primarily defined by selection in SB and typical size. Nevertheless, they stem from the same underlying gas, and the same underlying emission mechanisms, albeit potentially in different configurations. When we consider the same SB threshold and point-spread function of Matsuda et al. (2011), we find one LAB with an area >100 arcsec² and circularity of 0.33 (see bottom panel of Fig. 14). This is broadly consistent with the observational occurrence rate and observed circularities (Yang et al. 2010; Matsuda et al. 2011), although the small volume of TNG50 makes this comparison challenging.

4.3.3 The TNG model and physical gas state

With respect to underlying cosmological simulation, the main limitations of interest based on our use of TNG50 of the IllustrisTNG project are discussed in Byrohl et al. (2021), to which we refer the reader for more discussion.

The TNG simulations do not solve the equations of radiation-hydrodynamics (RHD), i.e. they do not directly incorporate radiation on-the-fly during the simulation. In the current context, this could be used to include currently missing local sources of ionizing radiation, namely stellar populations, as photoionization and photo-heating terms. While such simulations are increasingly feasible for in cosmological volumes down to z ∼ 5 in the study of cosmic reionization (Rosdahl et al. 2018; Ocvirk et al. 2020; Kannan, Garaldi & Smith 2022), as well as zoom simulations of individual galaxies (Rosdahl & Blaizot 2012; Mitchell et al. 2021; Costa et al. 2022), they are currently not feasible for large cosmological volumes to low redshift. This is an absolute requirement for predicting the observability of the Ly α cosmic web at z ∼ 2–3, which is why TNG50 and its combination of resolution and volume is our preferred tool.

While TNG50 is not an RHD simulation, it does incorporate simplified models of radiation in the most important regimes. Namely, the TNG model has an on-the-fly treatment for ionizing radiation from AGN, approximated as a spherically symmetric radiation field in the optically thin limit. Its impact on Ly α emission around galaxies is critically important (as discussed in Byrohl et al. 2021). Through photoheating and photoionization, AGN radiation significantly boosts emission from both recombinations and collisions (see Section 3.4). As a result, for haloes above 3 × 10¹⁰ M_⊙, Ly α luminosities increase by a factor of 10 for recombinations, and a factor of 2 for collisions. For reference, ∼40 per cent (∼80 per cent) of the gas mass (Ly α emission) in TNG50 at z = 2 arises from gas experiencing a local AGN radiation field in addition to the UVB. Beside the impact on the temperature and ionization state of the surrounding gas, AGN feedback drives large-scale outflows in massive TNG50 galaxies (Nelson et al. 2019b), transporting large amounts of gas out to large radii, reshaping the site and intensity of emission, and scattering of Ly α, which we will explore in future work.

In addition to impacting their surrounding gas, AGN can also be a powerful Ly α source. One could model Ly α emission from AGN with an effective model similar to our model for stellar populations by converting a fraction of the AGN bolometric luminosity into Ly α (see e.g. Costa et al. 2022). For the current study, we have neglected such emission model for three reasons: First, given the sparsity of broad-line Ly α emitting AGN (Spinoso et al. 2020), we only expect a low abundance in the TNG50 simulation volume, making their luminosity calibration difficult. Secondly, LAEs with broad-line Ly α features are more than an order of magnitude sparser than the LAFs studied here, hence demonstrating limited impact on our results. Thirdly, the overall luminosity contribution of fainter, more abundant AGN is effectively captured by the calibration to the observed LF. As we, however, ascribe all unresolved small-scale Ly α luminosity from AGN and stars to the latter, we introduce additional inaccuracies in the clustering and spectral shapes of these LAEs, which we hope to address in future work.

We note that TNG does not include a radiation field for stellar populations similar to that of AGN, such that the ionizing radiation from local stars is not included in our modelling. However, TNG does incorporate a spatially uniform, time variable metagalactic background radiation field (UVB; Faucher-Giguère et al. 2009), which is critically important for setting the physical state of low-density gas. In denser regions, gas can self-shield from this external radiation, and TNG includes a self-shielding model based on Rahmati et al. (2013) to reproduce the average ionizing field suppression. Collisional excitations are particularly sensitive to the gas state in the simulation due to the exponential temperature dependence around |${\sim }10^{4}\,$|K (Furlanetto et al. 2005; Faucher-Giguère et al. 2010), as we discuss below. Future variation runs of underlying galaxy formation simulations with respect to the AGN radiation field and the time-varying UVB intensity will allow to determine the robustness/sensitivity of Ly α observables to these processes.

While TNG50 is the highest resolution cosmological hydrodynamical simulation that exists at its volume (Nelson et al. 2019b), the halo mass function is only resolved, if we require a minimum dark matter particle count of 1000, down to halo masses of 4 × 10⁸ M_⊙, which corresponds to the median halo mass at an intrinsic Ly α luminosity of 2 × 10³⁹ erg s⁻¹ for our fiducial (rescaled) model. In our model, Ly α luminosities are drawn from the averaged initial mass function, as star particles in the TNG model represent entire stellar populations. Particularly, at the faint end, we thus do not incorporate the scatter introduced by sampling individual, massive stars that could potentially alter the Ly α emission in low-mass galaxies.

Despite the promising level of consistency between our model results and available Ly α observations, as discussed above, the feedback model of the TNG simulations, both in terms of expelled gas and as a photoionization source in the TNG model, is effective in nature and its details have significant theoretical uncertainty. Our methodology to predict Ly α emission across scales enables new empirical constraints on the underlying galaxy formation physics, through future, more rigorous comparisons with observations.

These constraints on the Ly α emitting gas probe scales and phases, distinct to other observable tracers, thereby providing orthogonal tests of gas in the TNG simulations, beyond existing explorations, e.g. total gas fraction (Pillepich et al. 2018b; Terrazas et al. 2020; Ramesh, Nelson & Pillepich 2023), absorption as well as emission from cool MgII (Nelson et al. 2019b, 2021), warm-hot OVI abundance Nelson et al. (2018b), hot X-ray emission (Truong et al. 2020; Truong, Pillepich & Werner 2021), the Sunyaev–Zeldovich signal (Pop et al. 2022), and the large-scale (re)distribution of gas by feedback (Ayromlou, Nelson & Pillepich 2022).

The observables of Ly α, from interstellar medium to IGM scales, are fundamentally linked to the underlying hydrodynamical simulation and its baryonic feedback models.

5 CONCLUSIONS

In this paper, we develop a comprehensive theoretical model for Ly α emission across scales, from individual LAEs to extended structures such as LAHs, blobs, and filaments of the cosmic web. To do so, we combine two ingredients: (i) TNG50 of the IllustrisTNG project, a high-resolution magnetohydrodynamical galaxy formation simulation in a cosmological volume of ∼50³ cMpc³ and (ii) a full Monte Carlo Ly α radiative transfer method to capture the impact of resonant scattering. We use the resulting Ly α predictions to explore the properties and detection rates of filamentary structures of the large-scale cosmic web in Ly α emission. Our main results are:

• Our Ly α emission model includes collisional excitations and recombinations in diffuse gas, as well as emission sourced by star-forming regions. We couple this to an empirically motivated dust treatment, which we calibrate against the observed LAE LF at z = 2 (Section 2). Our calibrated model reasonably reproduces other related observables: the galaxy UV LF, the Ly α escape fraction versus stellar mass, and the rest equivalent width distribution. It predicts a flattening Ly α luminosity versus galaxy/halo mass observable relation (Section 3.1).

• Within the fiducial model, radiative transfer of Ly α photons substantially alters the budget, spatial distribution, and observability within haloes, filaments, and voids. Collisional excitations in cold gas in the CGM of galaxies dominate the global emission budget. Star formation from within galaxies plays a significant role, and recombinations can be substantial around local ionizing radiation from nearby AGN (Section 3.2).

• Our main theoretical result is that diffuse Ly α cosmic web filaments are dominated by emission that originates from within galaxies and their gaseous haloes and not from the IGM itself. LAFs are therefore illuminated by scattered photons, which are intrinsically emitted within nearby haloes. This emission arises predominantly from ‘intermediate mass’ haloes between 10¹⁰ and 10¹¹ M_⊙. The ultimate physical origin is collisional excitations in their circumgalactic media as well as due to central, young stellar populations (Sections 3.3–3.4).

Considering the prospects for observing LAFs, we find that:

• Identifying filaments as extended regions (i.e. isophotal areas) above a SB threshold of 10⁻²⁰ erg s⁻¹ cm⁻² arcsec⁻², we predict that ‘large’ filaments with extent L ≥ 400 pkpc have a space number density of 2 × 10⁻³ cMpc⁻³. These Ly α structures directly trace large-scale dark matter filaments with characteristic overdensities of ∼10 and are responsible for |$\sim 60~{{\ \rm per\ cent}}$| of the global Ly α density |$\dot{\rho }_{\mathrm{Ly}\alpha }=4.4\cdot 10^{40}$| erg s⁻¹ cMpc⁻³ (Section 3.5).

• With respect to the relationship between LAFs and LAEs, we find that essentially all filaments host one or more low luminosity LAEs with L > 10³⁹ erg s⁻¹. However, only large filaments tend to contain a bright L > 10⁴¹ erg s⁻¹ LAE. In all cases, the diffuse regions (<10⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻²) of LAFs are powered by faint LAEs due to scattered photons. The emission arising within the IGM itself is negligible. (Section 4.1).

• We quantify the probability of detecting LAFs for both blind and targeted survey strategies, covering instrument footprint sizes ranging from KCWI, to MUSE, BlueMUSE, and HET-VIRUS. The SB levels of LAFs are within reach of current IFU instruments. However, at such moderate SB thresholds, filaments are sufficiently rare that large survey areas are needed (Section 4.2).

In the near future, new telescopes and instruments will continue to push to higher sensitivities and larger survey sizes, ultimately with the hope of imaging the large-scale cosmic web directly in Ly α emission. From the theoretical side, our overall model approach and its associated Ly α photon-level output will enable us to explore several promising directions. First, we can study the efficacy of more advanced filament detection and characterization algorithms. Secondly, we have not yet considered the spectral (i.e. wavelength) dimension of these data sets, the information content therein, including comparisons of emergent Ly α spectra with data.

Simultaneously, two key areas of improvement are the inclusion of intrinsic Ly α emission from AGN, and a more physically motivated dust treatment. To this end, galaxy formation simulations beyond TNG50 will enable us to push towards smaller spatial scales, permitting self-consistent dust models and explicit treatments of the escape of Ly α radiation from the star-forming ISM. Complexity is both a curse and a blessing, as the rich multiscale nature of Ly α emission simultaneously informs the physics of emission, radiative transfer, and the underlying galaxy formation model.

ACKNOWLEDGEMENTS

CB and DN acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Group (grant number NE 2441/1-1). We also thank the Hector Fellow Academy for their funding support. This work was further supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1–390900948 (the Heidelberg STRUCTURES Excellence Cluster). The TNG50 simulation was run with compute time granted by the Gauss Centre for Supercomputing (GCS) under Large-Scale Projects GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS). GCS is the alliance of the three national supercomputing centres HLRS (Universität Stuttgart), JSC (Forschungszentrum Jülich), and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF). Additional simulations and analysis were carried out on the Vera machine of the Max Planck Institute for Astronomy (MPIA) and systems at the Max Planck Computing and Data Facility (MPCDF).

DATA AVAILABILITY

Data directly related to this publication and its figures are 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).

Footnotes

References