© V. Mauerhofer et al. 2021

1. Introduction

Decades of observations show that the intergalactic medium (IGM) was completely ionized by redshift z ∼ 6, about one billion years after the Big Bang (e.g., Fan et al. 2006; Schroeder et al. 2013; Sobacchi & Mesinger 2015; Bañados et al. 2018; Inoue et al. 2018; Ouchi et al. 2018; Bosman et al. 2018; Kulkarni et al. 2019). This implies that, before this redshift, sources of ionizing radiation must have emitted enough photons that managed to escape their galactic environment and reach the IGM to photoionize it. The main candidates for the source of reionization are active galactic nuclei (AGNs) and massive stars. Although recent work suggests AGNs could be an important driver of reionization (Madau & Haardt 2015; Giallongo et al. 2019), a consensus seems to be emerging on the idea that their contribution is negligible compared to that of stars (Grissom et al. 2014; Parsa et al. 2018; Trebitsch et al. 2020).

Models based on observations show that stars alone could drive reionization under reasonable assumptions on the Lyman continuum (LyC) production efficiency and on the escape fraction of ionizing photons (f_esc), for example in Atek et al. (2015), Gnedin (2016) and Livermore et al. (2017). The escape fraction is the least constrained quantity, and is often fixed to a value of between 10 and 20% (e.g., Robertson et al. 2015). Finkelstein et al. (2019) argue that reionization can also be explained when assuming lower escape fractions. Additionally, simulations have successfully reionized the Universe with stellar radiation, either at very large scales (e.g., Iliev et al. 2014; Ocvirk et al. 2016), or at smaller scales with high enough resolution to resolve the escape of ionizing radiation from galaxies (e.g., Kimm et al. 2017; Rosdahl et al. 2018).

In order to better understand the process of reionization and to be able to compare observations and simulations, it is essential that we obtain accurate constraints on the LyC escape fractions of observed galaxies. Recently, several low-redshift galaxies were found to be leaking ionizing photons (e.g., Bergvall et al. 2006; Leitet et al. 2013; Leitherer et al. 2016; Puschnig et al. 2017). When this ionizing radiation flux is detected, it is possible to use a direct method to infer the escape fraction (e.g., Borthakur et al. 2014; Izotov et al. 2016a,b). The total LyC production is computed as the addition of the observed LyC photons and the ones absorbed in the galaxy, deduced from the flux of a (dust corrected) Balmer line, assuming that this line is produced by recombinations of photoionized hydrogen ions. The ratio of the observed ionizing flux to the assessed intrinsic LyC production gives the escape fraction.

However, the ionizing radiation escaping galaxies during the epoch of reionization is not directly observable because of IGM absorption, and therefore this method cannot be used. To circumvent this limitation, different ways to infer the escape fractions indirectly have been proposed in the literature. The connections between Lyman α (Lyα) and Lyman continuum emission have been studied (Verhamme et al. 2015; Dijkstra et al. 2016; Kimm et al. 2019) with promising results at z  ∼  3 (Steidel et al. 2018) or in the local Universe (Verhamme et al. 2017; Izotov et al. 2018). Above z  ∼  6, the intergalactic neutral hydrogen may alter Lyα too much, except for galaxies that have such big ionized bubbles that Lyα is redshifted enough before reaching the neutral IGM (Dijkstra 2014; Mason et al. 2018; Gronke et al. 2020). A possible alternative is to use Mg II λλ2796,2803 (Henry et al. 2018), a resonant emission line doublet which is a candidate proxy of Lyα and has the advantage of not being absorbed by the neutral IGM of the epoch of reionization. A second method is to use O32, which is defined as the line ratio [O IIIλ5007]/[O IIλλ3726,3729], and is thought to correlate with the escape fraction (Jaskot & Oey 2013; Nakajima & Ouchi 2014; Izotov et al. 2016b). However, recent studies shed doubt on this method (Izotov et al. 2018; Bassett et al. 2019; Katz et al. 2020). A high O32 ratio seems to be a good indicator of a nonzero escape fraction, but the correlation between the two quantities is weak.

A third way to infer the escape fractions indirectly is to use down-the-barrel absorption lines from low-ionization states (LISs) of metals, such as C II, O I or Si II lines. The shapes of these lines are used as an indicator of the covering fraction of the absorber. Indeed, column densities of metals in galaxies typically lead to saturated absorption lines (i.e., absorption profiles where the minimum flux reaches zero). When nonsaturated lines are observed, the residual flux is interpreted as being due to partial coverage of the sources by the absorbing material, meaning that the ratio of this residual flux to the continuum is used to deduce an escape fraction. As the LISs of metals are thought to be good tracers of neutral hydrogen, this escape fraction is then linked to the escape fraction of ionizing photons (e.g., Erb 2015; Steidel et al. 2018; Gazagnes et al. 2018; Chisholm et al. 2018)

There have not yet been observations of absorption lines at the epoch of reionization because it is challenging to detect the stellar continuum with a sufficiently high signal-to-noise ratio. However, there has been impactful progress in the study of absorption lines for two decades thanks to the observation of galaxies at ever higher redshifts and of local analogs of high-redshift galaxies. Detections of down-the-barrel absorption lines of individual galaxies with the highest redshifts, with z ∼ 3 − 4, are carried out using gravitational lenses (e.g., Smail et al. 2007; Dessauges-Zavadsky et al. 2010; Jones et al. 2013; Patrício et al. 2016). Alternatively, without gravitational lensing, it is possible to study absorption lines at these redshifts by stacking spectra of many galaxies (e.g., Shapley et al. 2003; Steidel et al. 2018; Feltre et al. 2020). The most common lines observed at these high redshifts are Si II (λ1190, λ1193, λ1260, λ1304), O I λ1302 or C II λ1334, which are strong lines whose observed wavelengths fall in the optical range. There are also observations of slightly redder absorption lines, such as Fe II (λ2344, λ2374, λ2382, λ2586, λ2600) or Mg II λλ2796,2803, at intermediate redshifts between 0.7 and 2.3 (e.g., Finley et al. 2017; Feltre et al. 2018). Finally, there are many studies of absorption lines at z <  0.3 (e.g., Rivera-Thorsen et al. 2015; Chisholm et al. 2017, 2018; Jaskot et al. 2019).

Because most atoms and ions in the interstellar medium (ISM) are in their ground state, absorption lines are typically transitions from the ground state to a higher level. The energy of these transitions translates into lines which are mostly in the rest-frame ultraviolet (UV). Some lines are resonant (e.g., Lyα, Mg II λλ2796,2803, Al II λ1670), but most are not completely resonant because the ground state is split into different spin levels, and de-excitation may produce photons with different wavelengths. If the wavelength of the so-called fluorescent channel(s) is different enough from the resonant wavelength the emitted photon will leave the resonance.

Observed absorption lines show a large diversity of spectral profiles, for example in terms of depth, width, and complexity. A common feature is that absorption is often blueward of the systemic velocity, with a minimum of the absorption occurring between around −400 and 100 km s⁻¹. There are also detections of fluorescent emission redward of the absorption (e.g., Shapley et al. 2003; Rivera-Thorsen et al. 2015; Finley et al. 2017; Steidel et al. 2018; Jaskot et al. 2019).

There are several physical processes that challenge the use of absorption lines to measure the covering fraction of neutral gas and the escape fraction of ionizing photons. One of them is the distribution of the velocity of the gas at galaxy scales. The gas in front of the sources moves at different velocities, leading to an absorption that is spread in wavelength. This can lead to an increase of the flux at the line center compared to a case where all the gas has the same velocity (Rivera-Thorsen et al. 2015). Another process is scattering: when photons are absorbed by an ion they are not destroyed but scattered. Thus, the flux in absorption lines can be increased by photons that were initially going away from the observer but that end up in the line of sight because of resonant scattering (Prochaska et al. 2011; Scarlata & Panagia 2015).

This paper investigates the relation between the escape fraction of ionizing photons and LIS absorption lines using mock spectra constructed from a radiation-hydrodynamic (RHD) simulation. There have been studies of various emission lines by post-processing simulations (e.g., Barrow et al. 2017, 2018; Katz et al. 2019, 2020; Pallottini et al. 2019; Corlies et al. 2020), and there are also studies on absorption lines in simulations, but these latter focus on quasar sightline absorption instead of down-the-barrel absorption (e.g., Hummels et al. 2017; Peeples et al. 2019). Concerning down-the-barrel lines, studies have been done using Monte-Carlo methods in idealized geometries (Prochaska et al. 2011) or using semi-analytic computations (Scarlata & Panagia 2015; Carr et al. 2018). Finally, there have been studies of mock down-the-barrel absorption lines in simulations (Kimm et al. 2011), but this present work is to our knowledge the first time that such lines are produced in a high-resolution RHD simulation and including resonant scattering. We focus in particular on the C II λ1334 line, which is among the strongest LIS absorption redward of Lyα, is not contaminated by absorption lines of other abundant elements (unlike O I lines), and is frequently used in the literature (e.g., Heckman et al. 2011; Rivera-Thorsen et al. 2015; Steidel et al. 2018). We also show results using the Si II λ1260 line. In addition, we study Lyβ to see if it is in general a better tracer of the escape fraction, even though it is not observable at the epoch of reionization because its wavelength is 1026 Å.

The paper is structured as follows: in Sect. 2 we provides details of our method, from the simulation that we use to the post-processing that is necessary to make mock observations and to compute escape fractions of ionizing photons. In Sect. 3 we present the resulting spectra and we assess their robustness against changes in modeling parameters. We also study the effect that various physical processes have on the spectra. In Sect. 4 we show the values of the escape fractions, and discuss the effects of helium and dust. In Sect. 5 we compare the line properties with the escape fractions to try to find correlations. We also explain the different sources of complexity of the relations we find. Finally, we summarize our results in Sect. 6.

2. Methods

In this section we provide details of the simulation that we use and the steps that are necessary to build the mock observations. In this paper we focus on the C II λ1334 absorption line but the method can be applied to any absorption line. In Sect. 5.5 we show the main results of the paper for Si II λ1260. Here we use one galaxy from a zoom-in simulation, and in future work we will study a statistical sample of galaxies from various simulations like SPHINX¹ (Rosdahl et al. 2018) or Obelisk (Trebitsch et al. 2020).

2.1. Simulation

We use a zoom-in simulation run with the adaptive mesh refinement (AMR) code RAMSES-RT (Teyssier 2002; Rosdahl et al. 2013; Rosdahl & Teyssier 2015), including a multi-group radiative transfer algorithm using a first-order moment method which employs the M1 closure for the Eddington tensor (Levermore 1984). The collisionless dark matter and stellar particles are evolved using a particle-mesh solver with cloud-in-cell interpolation (Guillet & Teyssier 2011). We compute the gas evolution by solving the Euler equations with a second-order Godunov scheme using the HLLC Riemann solver (Toro et al. 1994). The initial conditions are generated by MUSIC (Hahn & Abel 2013) and chosen such that the resulting galaxy at z = 3 has a stellar mass of around 10⁹ M_⊙. The physics of cooling, supernova feedback, and star-formation are the same as in the SPHINX simulations (Rosdahl et al. 2018). Radiative transfer allows for a self-consistent and on-the-fly propagation of ionizing photons in the simulation, which provides an accurate nonequilibrium ionization state of hydrogen and helium, as well as radiative feedback. We use three energy bins for the radiation: The first bin contains photons with energies between 13.6 eV and 24.59 eV, which ionize H⁰, the second bin between 24.59 eV and 54.42 eV, which ionize H⁰ and He⁰, and the last above 54.42 eV, which ionize H⁰, He⁰ and He⁺. The ionizing radiation is emitted from the stellar particles following version 2.0 of the BPASS² stellar library spectral energy distributions (SEDs; Eldridge et al. 2008; Stanway et al. 2016). In this version of BPASS the maximum mass of a star is 100 M_⊙, and all stars are considered to be in binaries. To account for the ionizing radiation produced by external galaxies that are not present in the zoom-in simulation, we add an ionizing UV background (UVB) to all cells of the simulation with n_H <  10⁻² cm⁻³, following Faucher-Giguère et al. (2009). The outputs of the simulation have a time resolution of 10 Myr, going down to redshift z = 3. The maximum cell resolution around z = 3 is 14 pc. The halo mass, stellar mass, stellar metallicity, and gas metallicity of the galaxy at z = 3 are M_h = 5.6 × 10¹⁰ M_⊙, M_* = 2.3 × 10⁹ M_⊙, Z_⋆ = 0.42 Z_⊙, and Z_gas = 0.43 Z_⊙, respectively.

In this paper we focus on three outputs of the simulation that span a range of ionizing photon escape fractions; we call these outputs A, B, and C. There are roughly 80 Myr of separation between each output. Output C is the last of the simulation, at z = 3. Properties of the outputs, including the star formation rate (SFR), are shown in Table 1. In Fig. 1 we plot the neutral hydrogen column density and the stellar surface brightness at 1500 Å at the virial radius scale and the ISM scale (a tenth of the virial radius), for output C. The column density map in the upper left panel shows gas that is being accreted from the IGM, and regions where the supernova feedback disrupts the gas flows. The upper right panel shows that there are many small satellites orbiting the galaxy. The lower panels highlight the high resolution of the simulation, where small star-forming clusters and dust lanes are well resolved.

Fig. 1. Maps of neutral hydrogen column density and surface brightness at 1500 Å for a snapshot of the simulated galaxy at z = 3. Upper left: column density of neutral hydrogen atoms, with the virial radius highlighted by the white circle. Upper right: intrinsic surface brightness at 1500 Å. Lower left: zoom-in to a tenth of the virial radius. Lower right: surface brightness at 1500 Å after extinction by dust.

2.2. Computing densities of metallic ions

To be able to make mock observations of metallic absorption lines in the simulation, we have to compute the density of the ions of interest in post-processing. The cosmological zoom-in simulation does not trace densities of elements other than hydrogen and helium, but it follows the value of the metal mass fraction Z in every cell. Assuming solar abundance ratios, we obtain the carbon density via the equation , where n_H is the hydrogen density, is the solar ratio of carbon atoms over hydrogen atoms, and Z_⊙ = 0.0134 is the solar metallicity (Grevesse et al. 2010). We do not consider the depletion of carbon into dust grains, which could alter the densities in the gas phase. However, Jenkins (2009) finds that, in the Milky Way (MW), not more than 20 − 50% of the carbon is depleted into dust. Moreover, De Cia (2018) shows that the depletion of metals into dust is less important in high-redshift galaxies than in the MW, which is why we can assume that the depletion does not change the gas-phase carbon density significantly. In order to assess the impact of the uncertainties of abundance ratios and dust depletion on our conclusions, we test the effect of dividing and multiplying the density of C⁺ by two in Sect. 3.2. The effect is visible but mild.

What remains is to compute the ionization fractions of carbon to get the densities of C⁺. To do so we use KROME³ (Grassi et al. 2014), which is a code that computes the evolution of any chemical network implemented by the user. As this computation is done in post-processing, we need to assume equilibrium values for the ionization fractions of carbon; we therefore use a routine of KROME that evolves the chemical network until it reaches equilibrium. For hydrogen and helium we use the (nonequilibrium) ionization fractions from the simulation.

The ionization fractions are determined by three main processes: recombination, collisional ionization, and photoionization. There can also be charge-transfer reactions between carbon ions and hydrogen or helium ions, but we omit them in this work, as we explain below. One of the major advantages of KROME is that the chemical network is completely customizable. The user chooses all the reactions and their rates as a function of temperature. We use the following rates:

• Recombination rates: Badnell (2006).

• Collisional ionization rates: Voronov (1997).

• Photoionization cross-sections: Verner et al. (1996).

We choose recombination and collisional ionization rates to reproduce the carbon ionization fractions as a function of temperature for a collisional ionization equilibrium setup in Cloudy⁴ (Ferland et al. 2017), which is the state-of-the-art tool to compute ionization fractions (among other things), but is too slow to use on a full simulation. There are works that use Cloudy on Ramses simulations, but these adopt strategies to avoid running it on every cell. Katz et al. (2019) use Cloudy on a subset of cells and then extrapolate on the other cells using machine learning algorithms. Pallottini et al. (2019) create Cloudy grids that are then interpolated for each cell of the simulation. We opt to use KROME because it is fast enough to directly compute the ionization fractions in every cell of the simulation, although it does not compute line emissivities of the gas. The collisional ionization rates that we use are slightly different from the ones in Cloudy (Dere 2007). We show in Appendix A that we recover the same ionization fractions of carbon as a function of temperature as Cloudy, despite the fact that we lack charge-transfer reactions and have different collisional ionization rates.

Concerning photoionization, the cross-sections that we use are the same as in Cloudy. To obtain a photoionization rate, those cross-sections have to be multiplied by a flux of photons. Taking advantage of the radiation-hydrodynamics in the simulation, we directly use the inhomogeneous radiation field self-consistently computed by the simulation at energies higher than 13.6 eV. Below 13.6 eV, the simulation does not track the radiation field or its interaction with gas. We nevertheless need to account for this lower energy range as it may ionize metals (e.g., the ionizing potential of C is 11.26 eV). For this, we add a simple model for radiation in the Habing band, from 6 eV to 13.6 eV (Habing 1968), in post-processing, using the UVB as inferred by Haardt & Madau (2012). In Sect. 3.2.2 we look at the effect of dividing the UVB by two or multiplying it by 100, and we see that this does not affect the results of the mock observations. We also use the Haardt & Madau (2012) UVB at all wavelengths instead of using the radiation field of the simulation, for comparison. In this case there can be drastic changes, which demonstrates that it is important to use a simulation that treats radiative transfer on-the-fly. In Appendix B we show in more detail how we compute the photoionization rates using the Verner et al. (1996) photoionization cross-sections and the number density of photons in the simulation.

Another advantage of KROME is that it is possible to fix the state of hydrogen and helium. If we compute the ionization fractions with Cloudy, the convergence to equilibrium of hydrogen and helium ionization fractions modifies the electron density compared to the value of the simulation. This different electron density would lead to different carbon ionization fractions than with the simulation electron density. To summarize, we compute the ionization fractions of carbon at equilibrium, but with the hydrogen and helium ionization fractions fixed at the simulation values. Therefore we set all rates of reactions containing hydrogen and helium to zero in KROME. This is why we do not include charge transfer reactions in KROME, as they would also change the simulation values of hydrogen and helium ionization fractions.

2.3. Fine structure of the ground state

Almost all ions have a ground state that is split into several spin states. The C⁺ ion has one such spin level, which we call level 2, just above the ground state (level 1), as illustrated in Fig. 2. A fraction of C⁺ ions are populating level 2 due to collisions with electrons, and they are absorbing photons at 1335.66 Å or 1335.71 Å instead of 1334.53 Å from level 1. We take this effect into account by computing the percentage of ions that are in level 2 using PyNeb⁵ (Luridiana et al. 2015). For every cell of the simulation we give PyNeb the temperature (assumed to be both the gas temperature and the electron temperature) and the electron density taken from the results of KROME. The fraction of C⁺ ions populating level 3 and 4 is extremely small, due to the short lifetime of these levels, and so we consider that the ions only populate levels 1 and 2.

Fig. 2. Energy levels of the C+ ion. P14 = 100% indicates that a C+ ion in level 1 can be photo-excited only to level 4. The green numbers are the probabilities that a C+ ion in level 4 radiatively de-excites to level 1 or 2. The red numbers are the probabilities that an incoming photon with a wavelength between 1335.66 Å and 1335.71 Å hitting a C+ ion in level 2 excites it to level 3 or 4. P32 = 100% indicates that a C+ ion in level 3 can only radiatively de-excite to level 2.

The three lines at 1334.53, 1335.66, and 1335.71 Å are considered in our radiative transfer post-processing, as explained in Sect. 2.4. The effects of these different channels on the C II λ1334 line are studied in Sect. 3.3.

2.4. Radiative transfer of the lines

The radiative transfer post-processing is computed with the code RAdiation SCattering in Astrophysical Simulations (RASCAS⁶; Michel-Dansac et al. 2020). The stellar luminosity of the simulated galaxy is sampled with one million so-called photon packets, each having a given wavelength and carrying one millionth of the total luminosity in the wavelength range considered (a few angströms on each side of the line). We launch the photon packets from the stellar particles of the simulation with the following initial conditions:

– Direction: The initial direction of every photon packet is randomly drawn from an isotropic distribution.

– Position: Each photon packet is randomly assigned to a stellar particle, with a weight proportional to the stellar luminosity in the wavelength range considered. The initial position of the packet is the position of the stellar particle.

– Frequency: The SED of the stellar particle, given by BPASS using the age and metallicity of the particle, is fitted by a power-law approximation, removing any stellar absorption lines. The frequency is then drawn randomly from this power law.

Photon packets are then propagated through the AMR grid of the simulation. The total optical depth in each cell is the sum of four terms, corresponding to four channels of interaction:

(1)

The optical depth of a line is the product of the cross-section and the column density, τ_line = σ_lineN_C II. The column density is the product of the ion density and the distance to the border of the cell⁷, N_C II = n_C IId. The cross-section for a given line of the ion is given by the formula

(2)

where f_line is the oscillator strength of the line, λ_line is the line wavelength, and b is the Doppler parameter, which is defined below (Eq. (3)). The variable a is defined by a = A_lineλ_line/(4πb), where A_line is the Einstein coefficient of the line. The variable x is defined by x = (c/b)(λ_line − λ_cell)/λ_cell, where λ_cell is the wavelength of the photon packet in the frame of reference of the cell. Finally, the Voigt function is computed with the approximation of Smith et al. (2015).

When an interaction occurs, one of the four channels is randomly chosen with a probability τ_channel/τ_cell. If the photon packet interacts with dust, it can either be scattered or absorbed (see Sect. 2.5). If the photon packet is absorbed by C⁺ in any channel, it will be re-emitted in a direction drawn from an isotropic distribution. If the absorption channel is 1334.53 Å, the photon packet will be re-emitted either via the same channel, which is called “resonant scattering”, or via the 1335.66 Å channel, which is the fluorescent channel. The probabilities of the two channels are determined by their Einstein coefficients, and are shown on Fig. 2. Similarly, if the photon packet is absorbed by the channel 1335.66 Å, it can be re-emitted at 1334.53 Å or at 1335.66 Å. Finally, if the photon packet is absorbed by the channel 1335.71 Å, it will always be re-emitted at 1335.71 Å. The transfer ends when all the photon packets either escape the virial radius or are destroyed by dust. More details, for example on the position of interaction or on the frequency redistribution after scattering, are given in Michel-Dansac et al. (2020).

Mock observations are made from different directions of observation thanks to the peeling-off method (e.g., Dijkstra 2017). This method allows us to make mock images, spectra, or data cubes. We can choose a spatial and spectral resolution and an aperture radius to make a mock observation of the whole galaxy or only of a part of it.

We also make mock Lyβ absorption lines. For this we use the density of neutral hydrogen directly from the simulation. No fine-structure level is taken into account because the resulting differences on the Lyβ wavelength would be of only 0.5 mÅ, which is too small for the photon to leave resonance. Lyβ is also a “semi-resonant” line: photons at the Lyβ wavelength can scatter on H⁰, leading to potential infilling effects, and, like in the case of C II λ1334, there is a way to leave the resonance, namely via an emission of an Hα photon. The probability for Lyβ to be resonantly scattered is 88.2%.

In this paper, the Doppler parameter b, which enters in the computation of the optical depth, depends on the thermal velocity and the turbulent velocity v_turb:

(3)

The turbulent velocity increases the probability of interaction of photons with wavelengths away from line center due to gas moving along or opposite to the direction of propagation, and therefore broadening the effective cross-section. This is important due to two limitations of the simulations that challenge a completely self-consistent treatment of absorption studies. The first limitation is the fact that the simulated galaxy is discretized on a grid, which means there are discrete jumps in the gas velocity from one cell of the grid to the next. Therefore, for some conditions, a photon packet with a wavelength corresponding to the velocity could go through two cells with projected velocities v and v + Δv without being absorbed, while it would be absorbed at intermediate velocities if the grid was more refined. The turbulent velocity is a way to smooth out those discontinuities. It is possible to estimate this quantity in the simulation by comparing the velocity in a cell with the velocities of neighboring cells, as is done for the star formation criteria of Kimm et al. (2017) or Rosdahl et al. (2018). Doing so in the three outputs of Table 1 yields an average n_H I-weighted turbulent velocity of approximately 12 km s⁻¹.

The second limitation of simulations is the fact that there are always gas motions at scales unresolved by the simulation, which would also increase the probability that photons get absorbed. This is called “subgrid turbulent velocities”, which is not taken into account in our simulation or in the previous computation of the turbulent velocity based on neighboring cells. Some other RAMSES simulations quantify those subgrid motions on-the-fly (e.g., Agertz et al. 2015; Kretschmer & Teyssier 2020). In this work we assume a uniform turbulence in every cell of the simulation, accounting for those two kinds of turbulence. We choose a fiducial value of v_turb = 20 km s⁻¹. In Sect. 3.2 we investigate the effects of changing this parameter.

2.5. Dust model

As described in Michel-Dansac et al. (2020), we follow Laursen et al. (2009) to model the effect of dust on radiation. This is a post-processing method, because the dust is not followed directly in the simulation. In this model, the optical depth of dust is the product of a pseudo-density and a cross-section. The pseudo-density is

(4)

where Z₀ is a reference metallicity and f_ion sets the quantity of dust that is put in ionized regions. Throughout this paper we use Z₀ = 0.005 and f_ion = 0.01, corresponding to the calibration on the Small Magellanic Cloud (SMC) of Laursen et al. (2009).

For the cross-section, we use the fits of Gnedin et al. (2008) in the SMC case. This cross-section accounts for both absorption and scattering events. RASCAS has an albedo parameter, representing the probability that a photon-packet that interacts with dust is scattered instead of destroyed. For this albedo we use the value of Li & Draine (2001), which is 0.276 for Lyβ and 0.345 for C II λ1334. Finally, the asymmetry parameter g of the Henyey-Greenstein function (Henyey & Greenstein 1941), which influences the probability distribution of the direction in which a scattered photon goes, is also taken from Li & Draine (2001), and is g = 0.721 for Lyβ and 0.717 for C II λ1334.

2.6. Computation of escape fractions

To compute the escape fractions of ionizing photons, we first create mock spectra of the galaxy between 10 Å and 912 Å using the method described in Sect. 2.4. The only changes are that, first, we use BPASS but without a power-law approximation of the SEDs, and second, the optical depth is computed differently. The optical depth for ionizing photons is a sum of contributions from H, He, He⁺, and dust. The first three are given by Verner et al. (1996) and dust is treated exactly as in Sect. 2.5, including the scattering, with an albedo of 0.24. The mock spectra are then used to compute the total number of ionizing photons going out of the virial radius, which is divided by the intrinsic number to get the escape fraction. We also compute the escape fraction at 900 Å by averaging the spectra between 890 Å and 910 Å. Results are shown in Table 1 and in Sect. 4.

3. Mock spectra

In this section we show a selection of C II λ1334 and Lyβ spectra. From the 5184 mock observations that we made (i.e., 1728 directions of observation for the three outputs A, B, and C), we select seven directions that have absorption profiles that are as diverse as possible, shown in Fig. 3. We highlight that those spectra do not represent the average of the galaxy over time and direction of observation, but they rather depict the large range of possible shapes and properties for one galaxy. Most lines are similar to the one in the top row of Fig. 3, that is to say saturated and wide.

Fig. 3. C II λ1334 and Lyβ absorption lines for seven directions of observation. The vertical dashed lines highlight the wavelength of the fluorescent channels. The spectra are arranged in increasing order of escape fraction of ionizing photons in the direction of observation, fesc (as computed in Sect. 2.6). The noise is Monte-Carlo noise, due to a sampling of the luminosity with a finite number of photon packets. The spectra are smoothed with a Gaussian kernel with FWHM = 20 km s−1. The aperture diameter for the mock observation is a fifth of the virial radius, which corresponds to about 5.6 kpc, or 0.75″ at z = 3. No observational effects such as noise, spectral binning, Earth atmosphere, or MW alterations are considered. We note that the x-axis in the first two rows of the second column differ from the others.

There is a very wide range of shapes of absorption lines for a single galaxy seen from different directions of observation and at time differences of 80 Myr. This shows that the distribution of gas in the galaxy is highly anisotropic. There is however an axis along which the galaxy is more flat, resulting in an irregular disk-like structure. There are slight trends between the properties of the lines and the inclination of this disk: edge-on angles have somewhat stronger absorption than face-on ones. We will explore those relations in a following paper.

The main properties of the lines on which we focus in this paper are the residual flux and the equivalent width (EW). The residual flux is

(5)

where is the flux at the deepest point of the absorption and is the flux of the continuum next to the absorption line. The EW corresponds to the area of the absorption line after the continuum has been normalized to one. The EW of the fluorescence corresponds to the area of the emission part above the continuum. Both EWs are defined as positive. We highlight that in several spectra the EW of fluorescence is larger than that of the absorption. This is because fluorescence is an emission process through which the gas emits in all directions. There are directions of observation with almost no absorption but with a strong fluorescence, which is emitted by gas that is excited by photons traveling in other directions. Lastly, we note that to compensate the movement of the galaxy in the simulation box we add v_CM ⋅ k_obs to the velocity axis of Fig. 3, where v_CM is the velocity of the center of mass of the stars and k_obs is the unit vector pointing toward the observer.

3.1. Comparison with observations

In comparison with observed z ∼ 3 galaxies, our Lyβ and C II λ1334 lines are slightly less wide. In our mock spectra we do not find absorption at very negative velocities. The median velocity at which the absorption reaches 90% of the continuum is −200 km s⁻¹, and the minimum is −450 km s⁻¹. This is to be contrasted with studies such as those of Steidel et al. (2010, 2018) and Heckman et al. (2015) where this velocity is often −800 km s⁻¹. We note however that these observations are averages over many galaxies and typically obtained with lower spectral resolution, meaning that a direct comparison is not really indicative. Furthermore, these absorptions at very negative velocities are produced in outflowing gas, and so an alternative explanation is that our simulation probes galaxies with weaker outflows. This is expected because our galaxy has a much lower mass and SFR than the Lyman Break Galaxies of those studies. However, it is possible that the simulation does not produce enough outflows to be realistic, as is argued in Mitchell et al. (2021), where they use a similar simulation. We will study the gas flows of our galaxy in a forthcoming paper.

In contrast, we find C II λ1334 profiles that are very similar to high-resolution observations of galaxies in the local Universe, as in Rivera-Thorsen et al. (2015) and Jaskot et al. (2019). These latter authors find similar absorption velocities as in our mock spectra, and a large variety of shapes, residual fluxes, and strengths of fluorescence. They even detect fluorescent lines with a P-Cygni profile (i.e., a blue part in absorption and a red part in emission), as in the fifth row of Fig. 3, due to the presence of C⁺ in the excited fine-structure level.

Our Lyβ lines are rather different from current observations of either high or low redshift galaxies. Steidel et al. (2018) and Gazagnes et al. (2020) show Lyβ absorption lines that have a large EW and a significant residual flux. Our lines are either wide with a small residual flux, or have a small EW and a large residual flux. This could be due to low signal-to-noise ratios, spectral resolution, or stacking effects. Alternatively, these effects could come from stellar absorption features or Lyβ emission from the gas, neither of which is included in this work.

We note that some Lyβ spectra display some emission redward of the absorption, especially in the second row of Fig. 3, and very little in the last five rows. This is due to scattering effects, and has not been observed yet because it is rare (the first row, which is the most common, has no emission) and requires an excellent signal-to-noise ratio.

3.2. Robustness of the method

Our procedure to make mock observations involves several free parameters. Here we test the extent to which the spectra are affected if we change their values. We focus in particular on the effect on residual flux, which is the quantity that we use the most in this paper. Figure 4 shows the same spectra as in Fig. 3 but with different alterations, which we explain individually below. Figure 5 shows the statistics of the change of residual flux for the 5184 directions of observation, for every setup that we try.

Fig. 4. Test of the robustness of the spectra against changing free parameters. The seven rows show the same seven directions as in Fig. 3. The vertical dashed lines highlight the wavelength of the fluorescent channels. First column: test of the effect of multiplying and dividing the fiducial carbon density by two. Second column: test of multiplying (dividing) the Habing band UVB by 100 (2) and comparison with the result with only UVB and no radiation field from the simulation. Third and fourth columns: test of changing the turbulent velocity for C II λ1334 and Lyβ respectively. The fiducial turbulent velocity is 20 km s−1.

Fig. 5. Effect of different parameters of the modeling on the residual fluxes. The y-axes display differences of residual fluxes between two setups. The dots are the median of the differences in a given bin of residual flux, and the bars show the 10th and 90th percentiles. The numbers indicate how many spectra are in the corresponding bin of residual flux. Left panel: effects of carbon density, which is divided or multiplied by two. Middle panel: effect of UVB radiation for ionizing photons instead of the radiation field from the simulation. Right panel: effects of turbulent velocity. Turbulent velocities are varied between 0 and 40 km s−1.

3.2.1. Carbon abundance

As discussed in Sect. 2.2, the simulation does not trace carbon directly and we have to assume an abundance ratio for carbon to infer it in post-processing. This setup incurs uncertainties due to the poorly constrained carbon-to-hydrogen ratio and dust depletion factor in z = 3 galaxies, and therefore we vary the density of carbon to see if the result is sensitive to that variation. The first column of Fig. 4 shows the effects of dividing and multiplying the density of carbon by a factor of two. Overall the effect is small. The fluorescence of the second and sixth spectra is affected the most, while the absorption part is never drastically changed.

For a more global view, the left panel of Fig. 5 shows the effects on the 5184 directions. We see that directions with a very small or a very large residual flux are not significantly affected by the change in carbon density. This is because when the residual flux is either very small or very large the line photons have gone through mostly optically thick or thin media, respectively, which are not strongly affected by a change of a factor of two. Alternatively, when the residual flux is between ∼0.2 and 0.8, the gas through which the photons went is neither completely thick nor thin, and therefore a change in density of a factor of two has more of an effect, with a difference of residual flux of up to 0.25. As the carbon-to-hydrogen ratio is poorly constrained, we use solar abundance ratios and test the relations between the absorption lines and the escape fraction of ionizing photons under this assumption.

3.2.2. Implementation of the radiation below 13.6 eV

As the simulation was done with only hydrogen- and helium-ionizing photons, we have to estimate in post-processing the flux of photons between 6 eV and 13.6 eV, which ionize neutral carbon, among other metals. Our fiducial model is to use the Haardt & Madau (2012) UVB in every cell of the simulation, except in cells where n_H I > 10² cm⁻³, because the optical depth of dust in those cells is enough to extinguish the UVB. We note that removing the UVB in those cells does not impact the absorption lines at all because those cells are so dusty that they hide all the stellar continuum. Additionally, as explained in Sect. 2.2, we take the ionization fractions of hydrogen and helium from the simulation, and therefore we remove all reactions including hydrogen or helium from our KROME chemical network.

In the second column of Fig. 4 we compare the fiducial spectra with three other setups:

1. The same as fiducial, but dividing the flux between 6 eV and 13.6 eV by two in every cell. This is done to evaluate the effect of reduction of the UVB due to absorption in the galaxy.

2. The same as fiducial, but multiplying the flux between 6 eV and 13.6 eV by 100 in every cell. This is done to simulate the contribution of a young metal-poor stellar population of 1000 M_⊙ at a distance of 10 pc (which is roughly the resolution in the ISM) at this energy range.

3. Using the Haardt & Madau (2012) UVB at all energies, ignoring the radiation field from the simulation. We let KROME compute the ionization fractions of every species, including hydrogen and helium, unlike in the fiducial setup. We compute the photoionization rate that the UVB imposes on H⁰, He⁰, He⁺, C⁰ and C⁺ at the redshift of the output, and we apply this photoionization rate in every cell with a hydrogen density smaller than a certain threshold. This threshold is 10⁻² cm⁻³ for H⁰ photoionization, 4.3 × 10⁻² cm⁻³ for He⁰, 3.3 × 10⁻¹ cm⁻³ for He⁺, 10² cm⁻³ for C⁰ and 4.3 × 10⁻² cm⁻³ for C⁺. For each ion this threshold is chosen so that a cell with a size of 10 pc having this density would have an optical depth of around one for photons having the smallest ionizing energy of the ion. We test this setup to assess the importance of using an RHD simulation instead of a hydrodynamics simulation where only a UVB can be used for the photoionization.

The two first setups are represented by the dashed lines in the second column of Fig. 4. They are identical to the fiducial spectra. This is because the Haardt & Madau (2012) UVB below 13.6 eV is high enough to completely photoionize C⁰ to C⁺ (in all the regions with n_H I < 10² cm⁻³), and therefore multiplying it by 100 does not change the result, and dividing it by two is not enough to stop photoionizing C⁰. To confirm this we multiplied the UVB by a factor 10 000, because locally the radiation field below 13.6 eV can be even stronger than 100 times the UVB, and the spectra are still not affected. On average, over the 5184 directions the residual flux is changed by only 0.002, both when we divide the UVB by two or multiply it by 100.

In the third setup, the dotted lines in the second column of Fig. 4, the effects are more noticeable. The first, third, and fourth rows are not affected significantly, while the second and the last three are drastically changed. The middle panel of Fig. 5 shows the differences for the 5184 directions. We see that the directions with small residual fluxes are not affected significantly, such as the purple spectrum in the first row of Fig. 4, but that the directions with high fiducial residual fluxes tend to have a much deeper absorption line when we use the UVB instead of the radiation field of the simulation, which is the case for the last three rows of the second column of Fig. 4. This shows that using a UVB at all energies instead of a self-consistent ionizing radiation field simulated on-the-fly does not photoionize C⁺ to C⁺⁺ efficiently enough, and therefore globally overestimates the density of C⁺, which leads to deeper absorption lines. Another factor that explains the change of C⁺ density is the modification of the C⁺ recombination and collisional ionization rates due to the change of electron density when we let KROME compute the equilibrium ionization fractions of hydrogen and helium instead of using the values from the simulation. Those results show that it is important to use an RHD simulation to accurately compute the ionization fractions of metals in post-processing.

3.2.3. Turbulent velocity

Turbulent velocity (with a fiducial value of 20 km s⁻¹) is a parameter we add to smooth out the discontinuities of the velocity field in the simulation and to model the subgrid motion of the gas, as we introduced in Sect. 2.4. The third and fourth columns of Fig. 4 show the effect on C II λ1334 and Lyβ of changing the turbulent velocity to 0 and 40 km s⁻¹. We see that Lyβ is almost unchanged, and for C II λ1334, the absorption parts of our seven spectra are not very sensitive to the turbulent velocity either, except the one in the first row. However, fluorescence grows quickly with turbulent velocity, indicating that there is more absorption overall in the galaxy, leading to a stronger fluorescent emission in every direction of observation.

The right panel of Fig. 5 shows the effects of changing the turbulent velocity on our 5184 mock spectra. The C II λ1334 line shows large differences of residual fluxes, of up to 0.25. An average difference of 0.15 in residual flux is seen for directions with a low residual flux, and a smaller difference is seen for directions with a high residual flux. There are even directions where the residual flux is larger for a turbulent velocity of 40 km s⁻¹, which may seem counter-intuitive. This is because a larger turbulent velocity results in more absorption globally in the galaxy, and so there are also more scattering events and more infilling effects in some directions of observation. Some directions are sensitive to infilling, as we show in Sect. 3.3.2, and this infilling can increase the residual fluxes. Figure 5 confirms that Lyβ is less affected by the turbulent velocity parameter than C II λ1334. The difference of residual fluxes with a turbulent velocity of 0 and 40 km s⁻¹ is on average smaller than 0.1. The turbulence affects Lyβ less than C II λ1334 because the thermal velocity of hydrogen is about 3.5 times larger than that of carbon.

3.3. Effects of various physical processes

In this section we analyze the effects of various physical processes on the spectra to show that the formation of the line is complex. In Figs. 6, 9, and 10 we show the same spectra as in Fig. 3 but altered by removing a number of processes that are taken into account in the fiducial method. Figures 7 and 8 show the statistics of the change of residual flux for the 5184 directions of observation for some of the processes that we analyze. Below we explain these processes one by one.

Fig. 6. Illustration of the effects of different physical processes on the same spectra as Fig. 3. The vertical dashed lines highlight the wavelength of the fluorescent channels. First column: comparison with and without a fraction of C+ ions populating their ground-state fine-structure level. Second column: effect of removing the scattering of photons. Third column: effect of changing the aperture of the observation.

Fig. 7. Same as Fig. 5, to show the effects of different physical processes. Left, middle, and right panels: effects of scattering, dust, and gas velocity, respectively.

Fig. 8. Same as left panel of Fig. 7, in bins of EW of the fluorescent line instead of residual flux. The best fit is 0.06 for EWfluo <  0.6 Å and 0.11 EWfluo − 0.003 otherwise.

Fig. 9. Same as Figs. 6 and 10 but for Lyβ. No fine-structure level is considered here.

Fig. 10. Same as Fig. 6, with other processes. First column: effect of removing dust. We also remove the fine-structure level to simplify the interpretation. Second column: spectra with and without gas motion in the galaxy. Infilling effects give complex shapes to the second and seventh rows. Third column: comparison of fiducial spectra with spectra with neither gas motion, dust, scattering, nor fine-structure level.

3.3.1. Fine-structure level

Around 5% of C⁺ ions in the ISM of our galaxy are excited in the fine-structure level of the ground state, leading to absorption at a different wavelength from the true ground state, as explained in Sects. 2.3 and 2.4. One of the resulting differences is that C⁺ ions in the excited state can resonantly scatter fluorescent photons at 1335.66 Å and 1335.71 Å. The first column of Fig. 6 shows how spectra change if we assume instead that all C⁺ ions are in the true ground state.

The absorption parts of the spectra are almost unaffected, with an average difference of residual flux of less than 0.01. However, there is significant modification of the fluorescent emission. The distribution of EW of fluorescence is redistributed among the different directions of observations when ions in the fine-structure level are considered because of scattering of resonant photons on those ions. This leads to some directions having a larger EW of fluorescence and other directions a smaller one.

Furthermore, the peak of the fluorescence is redshifted when the fine structure is present, once more due to the scattering of photons in the 1335.66 Å and 1335.71 Å channels. Those photons are trapped by C⁺ ions in the fine-structure level; they can only escape if their wavelength is redshifted. The ones that are blueshifted are absorbed by C⁺ ions in their true ground state. This behavior is reflected in the statistics of the number of scattering events of the photons. Without fine-structure level splitting there is up to a maximum of 65 scattering events per photon, while this number becomes ∼27 000 when the fine-structure level is added. Finally, the presence of the fine-structure level can lead to fluorescence with a P-Cygni profile, as illustrated in the first column and fifth row of Fig. 6.

3.3.2. Infilling

The scattering of C II λ1334 or Lyβ photons can lead to an effect called infilling. Without scattering, all the photons that are absorbed are destroyed. However, when scattering is taken into account, the photons are re-emitted, either as an Hα photon for Lyβ or a fluorescent photon for C II λ1334, or at the wavelength of the line. In this case, photons that are absorbed can still contribute to the spectra if they are re-emitted resonantly and manage to escape the galaxy. Furthermore, scattered photons can contribute to the spectrum of a different direction of observation than the initial emission. This is the infilling effect. In order to see how much infilling occurs in our mock spectra, we generate spectra in which any absorption leads to the destruction of the photon. This means that the resulting spectra are the sum of Voigt absorption profiles for every cell in front of every stellar particle in the simulation. In this case, no fluorescent photons are emitted, and so the fluorescent C II⋆ λ1335 line disappears.

The second column of Fig. 6 for C II λ1334 and the first column of Fig. 9 for Lyβ show the fiducial spectra and the spectra without scattering. For C II, we see a second absorption around 250 km s⁻¹ coming from the presence of C⁺ ions in their fine-structure level. The infilling in the absorption is strong for the second row, with a residual flux that is doubled when scattering is included. For the other spectra, the residual flux is only increased by a few percent.

The left panel of Fig. 7 shows the differences of residual flux with and without scattering for the 5184 directions. Both C II λ1334 and Lyβ show the same trends. For residual fluxes below 0.1 the median of the difference is small (less than 0.03). For very high residual fluxes, the difference is larger, with a median of 0.12 for Lyβ and 0.1 for C II λ1334. In between, for residual fluxes from 0.1 to 0.8, the median of the difference varies between 0.05 and 0.1. Overall, the effect of infilling is not very significant on average, but there are extreme directions where the infilling increases the residual flux by more than 0.3.

Figure 8 shows the differences of residual flux with and without scattering, similar as in Fig. 7, but in bins of EW of fluorescence. We choose these bins to show that the effect of infilling is more important for directions with a strong fluorescence. This highlights the fact that there are directions where many scattered photons escape the galaxy. Indeed, both fluorescent photons and photons creating the infilling effect must have scattered at least once. Figure 8 shows that fluorescent photons and resonantly scattered photons both escape the galaxy in preferential directions. This can be used to get a rough estimate of the infilling effect based on the fluorescent line (see Fig. 8).

3.3.3. Aperture size

The effect of infilling depends on the aperture used for the observation, because photons can potentially travel far before being scattered and redirected toward the direction of observation (Scarlata & Panagia 2015). The third column of Fig. 6 for C II λ1334 and the second column of Fig. 9 for Lyβ compare the fiducial spectra with different aperture radii: R_vir/25, R_vir/10 and R_vir, which are approximately 1.1, 2.8, and 28 kpc, or angular diameters of about , , and 7″ at z = 3. Our fiducial spectra are made with the R_vir/10 aperture, because it includes as much as ≈95% of the stellar continuum emission. The R_vir aperture is used to assess the effect of adding the light scattered in the circum-galactic medium (CGM) and the R_vir/25 aperture, which contains ≈40% of the stellar continuum, is useful to show what happens if we do not enclose the whole stellar-continuum -emitting region in the mock observation. We highlight that the aperture is a property of the mock observation, and it does not affect the radius of propagation of the photon packets, which is always one virial radius (Sect. 2.4).

The result is that the residual flux is often barely affected by the aperture of observation, except for the small aperture in the second, third, and fourth rows of Fig. 9. It is noticeable in general that smaller apertures lead to smaller residual fluxes. This is in part because a smaller aperture misses some of the scattered photons, but is mainly because a smaller aperture misses flux from stars that are facing optically thin gas. The light from those stars is not absorbed, and therefore it increases the residual flux if the stars are in the aperture of the observation. Similarly, the fluorescence increases for larger apertures, as does the emission part of Lyβ, because there are photons scattering far away that are missed with smaller apertures. Nevertheless, the differences when taking a larger aperture are small, and the aperture size has little importance as long as it encompass more than ≈90% of the stellar continuum. For example, an aperture corresponding to 1″ gives a result very similar to our fiducial aperture. As discussed in Mitchell et al. (2021), where a simulation similar to ours is used, the small effect of increasing the aperture may be due to a lack of neutral gas in the CGM, possibly due to imperfect supernova feedback modeling.

3.3.4. Dust

Now we study the effects of dust on the spectra. The third column of Fig. 9 for Lyβ and the first column of Fig. 10 for C II λ1334 compare the fiducial spectra with the spectra made without dust. To avoid obtaining a complex shape for the fluorescence, which can be misleading in this case, we omit here the fine-structure of C⁺. Every C⁺ ion is assumed to be in the ground state. All the spectra that we plot are normalized, and so we do not see the difference in the continuum when dust is removed, which would be significantly attenuated. Averaging over the 5184 spectra, the continuum is reduced by 74% at 1330 Å and by 83% at 1020 Å.

One striking difference in the spectra is that the Lyβ lines have especially large EWs when dust is removed. This is because there are stars embedded in optically thick regions, with N_H I > 10²² cm⁻², causing absorption up to high velocities, yet those stars do not contribute to the spectra when dust is present because they are completely attenuated. Another result is that the absorption lines have smaller residual fluxes when dust is absent. This is because the flux at the deepest point of the line does not change significantly with or without dust, whereas the flux of the continuum does. Thus the residual flux, which is the ratio of those two fluxes, is larger with dust than without. The flux of the line center does not change significantly because dust and H⁰ or C⁺ all reside in the neutral regions, and therefore the line center photons have generally a higher probability of being absorbed by the gas than by dust.

The middle panel of Fig. 7 shows the difference in residual flux between the spectra with and without dust for the 5184 directions. For very small fiducial residual fluxes, the difference is almost zero. However, for directions with higher residual fluxes, the spectra without dust have a much smaller residual flux than the ones with dust.

3.3.5. Gas velocity

The gas in front of the stars moves at different velocities, leading to an absorption that is spread in wavelength. This can lead to an increase of the residual flux compared to a case where all the gas has the same velocity. For example, even if each star is facing a column density high enough to saturate the absorption line, the residual flux can be larger than zero if the gas is not going at the same velocity for each star. For the line to be saturated, with zero residual flux, there must be at least one velocity regime for which all the stars are facing a sufficiently high column density of gas going at this (projected) velocity.

To show that the spectra in our simulation are sensitive to the velocity of the gas, the fourth column of Fig. 9 for Lyβ and the second column of Fig. 10 for C II λ1334 show the spectra with and without gas velocity. The EW is generally smaller without gas velocity, because the absorption is not spread in wavelength. Moreover, there are scattering effects that give complex shapes to the absorption profiles, especially in the second row for both lines and in the seventh row for C II λ1334.

The right panel of Fig. 7 shows the difference of residual flux between the spectra with and without gas velocity. As for the case without dust, the difference is negligible for very small fiducial residual fluxes and increases for higher residual fluxes, especially for Lyβ. For very high residual fluxes, above 0.6, the difference becomes smaller again. For C II λ1334, there are even directions where the residual flux is larger when the gas velocity is zero, because of scattering effects.

3.3.6. All processes removed

Finally, we show in the last column of Fig. 9 for Lyβ and of Fig. 10 for C II λ1334 the spectra with all four processes removed: fine-structure level, scattering, dust absorption, and gas velocity. Those spectra are more saturated because all four processes that we remove tend to increase the residual flux. The spectra still do not resemble Voigt profiles because of the fact that there are many sources, each facing a different column density of gas.

To summarize this section, we find that our predictions for the shape of the C II λ1334 and Lyβ lines are rather robust against the free parameters of our modeling, such as metal abundances, nonionizing UVB estimations, and subgrid turbulent velocity. However, there are different processes that are crucial for the formation of the absorption lines. The absorption by dust and the velocity field are the most important, while the absorption from the fine-structure level, the infilling, and the aperture effects alter the lines to a lesser extent.

Additionally, dust has a counter-intuitive impact on the shape of absorption lines: young and luminous stars are often not contributing to the spectrum because of selective suppression of the stellar continuum emerging from the dense and dusty star-forming regions (see also Sect. 5.1), which is why dust increases the residual flux.

4. Escape fractions of LyC

We now focus on the escape fractions of ionizing photons in the simulated galaxy. In order to compare the absorption lines and the escape fractions of ionizing photons, we compute the escape fraction from the galaxy as seen from different directions in post-processing, as explained in Sect. 2.6. The correlations between the mock absorption lines and the escape fractions are studied in Sect. 5.

4.1. Distribution of escape fractions

The resulting escape fractions averaged over all directions are 11% for output A (z = 3.2), 3.8% for output B (z = 3.1) and 1% for output C (z = 3.0). Output A has the largest escape fraction of all outputs with z <  7 in this simulation. Around 10% of the last 35 outputs, at z <  3.5, are similar to output B, with an average f_esc ∼ 3 − 4%. The vast majority, around 90% of outputs, have similar escape fractions to that of output C, with an average f_esc ≤ 1%. Such large variations of the escape fraction with time are expected because the escape of ionizing photons is regulated by small-scale processes, essentially the disruption of star-forming clouds by radiation and supernova explosions, which have typical timescales of the order of a few million years (e.g., Kimm et al. 2017; Trebitsch et al. 2017). Also, at output A, the galaxy has been forming stars relatively intensely over the past ∼10 Myr, at around 5 M_⊙ yr⁻¹, relative to the average of around 3 M_⊙ yr⁻¹ from z = 3.5 to z = 3, and is thus undergoing strong feedback. As can be seen from Sect. 5.3.1 below, this results in lowering the covering fraction of neutral gas in front of stars, and therefore increasing the escape of ionizing radiation. We choose these outputs, in particular output A, to study the largest variety of escape fractions possible, despite the fact that output A is not representative of the typical state of the galaxy.

Figure 11 shows the histogram of escape fractions in the 1728 directions of observation of the galaxy, for the three epochs A, B, and C. The three outputs show very different distributions of escape fractions, even though they are the same galaxy with 80 Myr between the outputs. As for the spectra in Sect. 3, there is a large variety of escape fractions depending on the direction of observation of the galaxy. This large variety implies that a measurement of escape fraction of a galaxy does not necessarily give the correct contribution of the galaxy to the ionizing photon background (e.g., Cen & Kimm 2015). One needs a statistical study of a large sample of galaxies to get information on the potential of galaxies to reionize the Universe.

Fig. 11. Histogram of LyC escape fractions for the three outputs, computed as in Sect. 2.6.

Additionally, we show in Fig. 12 the comparison between the global escape fraction of ionizing photons and the escape fraction at 900 Å, . We find that is mostly smaller than f_esc, from a few percent to around 12% lower. There are a few points with small escape fractions where is slightly larger than f_esc because of the effect of helium absorption.

Fig. 12. Comparison of the escape fraction at 900 Å and the global escape fraction of ionizing photons.

4.2. The effect of helium and dust

As described in Sect. 2.6, the computation of escape fractions includes helium and dust in addition to hydrogen. We want to highlight that in our simulation the role of helium and dust is very subdominant, as also found by Kimm et al. (2019). If we remove helium, the escape fractions smaller than 5% increase on average by only about 0.15% (absolute difference), while the larger escape fractions increase by around 1 − 3%. Those differences are small because helium absorbs only ionizing photons with energies larger than 24.59 eV, which are less numerous than the ones between 13.6 eV and 24.59 eV.

When dust is removed, the change is even smaller, with an average increase of escape fractions of 0.05%, and a maximum increase of 1.1%. This small effect comes from our dust distribution model (Eq. (4)), where the density of dust is proportional to the density of neutral hydrogen. This implies that for ionizing photons the optical depth of dust is always much smaller than the optical depth of H⁰. In other words, dust alone would have a strong impact on the ionizing photons, but neutral hydrogen absorbs them much more efficiently, such that dust becomes ineffective. This is in agreement with other studies that find little or no effect of dust on the escape fraction (Yoo et al. 2020; Ma et al. 2020). However, in Ma et al. (2020) the escape fractions in the most massive galaxies are substantially affected by dust, but these galaxies are more massive than ours. Additionally, their dust modeling puts dust in gas with a temperature of up to 10⁶ K, whereas in our simulation dust is proportional to n_H I + 0.01n_H II, and so the dust is almost completely absent in gas hotter than 3 × 10⁴ K. We note that a model of dust creation and destruction in the simulation could yield a different dust distribution, for example with more dust in ionized regions, which could then affect ionizing photons more, but this is beyond the scope of this paper.

5. The link between escape fractions and absorption lines

Now we turn to the search of correlations between properties of absorption lines and the escape fraction of ionizing photons. Figure 13 shows the relations between the escape fractions and some properties of C II λ1334 on the left and Lyβ on the right, for the three outputs and the 1728 directions of observation of each output. We find that none of the observables give a clear indication of escape fraction in all ranges of f_esc. We can however see some general trends.

Fig. 13. Scatter plots comparing the escape fractions for the three outputs in the 1728 directions of observation to the properties of the C II absorption line on the left and the Lyβ line on the right. The seven colored dots show the position of the spectra displayed in Sect. 3. EW_abs is the EW of the absorption line (negative EW_abs means more emission than absorption), EW_fluo is the EW of the fluorescent emission, vmax is the velocity of the deepest point in the spectrum, vcen is the velocity at which the absorption line has half its EW, and v90 is the velocity where the blue side of the absorption line reaches 90% of the continuum value. The orange vertical lines show the median value of the x-axis variable over the 5184 directions. There is no fluorescent channel for Lyβ, and so one panel is missing.

For C II, the escape fraction is almost always smaller than the residual flux, which may therefore be used as an upper limit for the escape fraction. This means that a zero residual flux implies f_esc = 0. The EW of absorption may also be used to get an upper limit because the following relation holds approximately (for our simulated low-mass galaxy): f_esc <  0.5 − 0.4 × EW_abs [Å]. In particular, directions of observation with EW_abs >  1.25 Å have f_esc <  2%. The EW of the fluorescence is a poorer indicator of f_esc than the EW of the absorption, but a high fluorescent EW also implies small f_esc. More precisely, EW_fluo >  1 Å implies f_esc <  2%. This is compatible with Jaskot & Oey (2014), who find a strong fluorescent line in a galaxy whose Lyα profile hints at a nonzero escape of LyC. Finally, the three velocity indicators, v_max, v_cen, and v₉₀ show the least correlation with f_esc. One potential relation is that directions with v_cen <  −250 km s⁻¹ have f_esc ≳ 10%, which is a similar behavior to what is observed by Chisholm et al. (2017). However we do not have enough data points in this velocity regime to be confident that f_esc could not be smaller.

For Lyβ, the relations are even more scattered. The residual flux does not seem to trace f_esc at all. A high EW again implies a small escape fraction. Indeed, we find that f_esc <  2% when EW_abs >  3 Å. The directions with EW_abs <  −0.6 Å, that is to say with Lyβ dominantly in emission, have f_esc <  2%. However, there are not many points, and so this regime might be undersampled. As in the case of C II, a very small v_cen implies a high f_esc, but again the number of points is too small to be confident. Finally, a very small v₉₀ is a sign of low f_esc: v₉₀ <  −700 km s⁻¹ implies f_esc <  2%, except for two outliers which are above this relation. However, those rare values of v₉₀ are arising only when the EW is very large, and therefore they do not bring new information.

In summary, there is one regime where the absorption lines can provide reliable information about f_esc. When the lines are saturated and wide, with one attribute generally accompanying the other, f_esc is smaller than 2%. For Lyβ, being saturated is not a sufficient condition to have a low escape fraction. In addition, the EW has to be larger than 3 Å. A large EW of fluorescence also implies a low f_esc, and is sometimes not linked with a large EW of absorption, as is the case for the large blue dot in Fig. 13. For the other regimes, there is always significant scatter and no predictive relations. To gain more confidence in the relations we find here, they must be verified in other simulated galaxies with different masses and metallicities, which we postpone to future work.

It might be surprising that there is little correlation between the escape fractions and the properties of the absorption lines, especially for the residual flux. In contrast, it is common in the literature to use the residual flux of LIS absorption lines to evaluate the covering fraction of neutral gas and deduce an estimate of the escape fraction of ionizing photons (e.g., Heckman et al. 2001; Shapley et al. 2003; Grimes et al. 2009; Jones et al. 2013; Borthakur et al. 2014; Alexandroff et al. 2015; Reddy et al. 2016; Vasei et al. 2016; Chisholm et al. 2018; Steidel et al. 2018). This is because there are idealized geometries where the escape fraction of ionizing photons and the residual flux of the LIS absorption lines are equal. If we imagine a plane of stars covered by a slab of gas composed of optically thick clouds and empty holes, the ionizing photons and the line photons can only escape through the holes, and not through the thick clouds. In this case, the escape fraction of ionizing photons and the residual flux of the lines are the same, and are equal to the fraction of the surface of the sources that is behind holes, which is referred to as the covering fraction. This is usually called the picket-fence model, implemented with some variations (e.g., Reddy et al. 2016; Steidel et al. 2018; Gazagnes et al. 2018). Several authors have warned that the residual flux can only give a lower limit on the covering fraction, for example due to the velocity distribution of the gas (e.g., Heckman et al. 2011; Jones et al. 2012; Rivera-Thorsen et al. 2015). We argue that there are other effects that complicate the relations between the lines and the escape fractions of ionizing photons, that we detail in this section.

There are several assumptions linked with the picket-fence model. The sources of ionizing photons and nonionizing UV continuum have to be the same so that they trace the same screens of gas. If the continuum of the lines and the ionizing photons are passing through different media, it is impossible for the absorption lines to be a perfect tracer of the escape fraction of ionizing photons. Concerning dust, there are different treatments in the literature. Some studies do not apply dust corrections, similarly as in Fig. 13. The assumption for the residual flux of absorption lines to be equal to the escape fraction of ionizing photons without dust correction is based on dust being homogeneous, inside and outside of the holes, and the idea that the ionizing photons are affected by dust in the same way as the nonionizing UV continuum. In other studies, the residual flux is corrected under the assumption that the holes have no dust, or to account for the different effects of dust on ionizing and nonionizing photons (e.g., Steidel et al. 2018; Gazagnes et al. 2018; Chisholm et al. 2018). Additionally, one assumption is that the distribution of column densities has to have optically thick parts and holes empty enough to be optically thin for both the ionizing photons and the line photons. In this setup, it is straightforward to define a covering fraction, to measure it with the residual flux of absorption lines, and to deduce the escape fraction as one minus the covering fraction. If there are column densities such that the optical depth is around one, the concept of covering fraction is less well defined, and the escape fraction is no longer one minus the covering fraction. A final assumption is that the infilling effects and the gas velocities do not significantly affect the absorption lines, which would alter the measurement of the covering fraction.

In the following, we compare our simulation with the picket-fence model, addressing these assumptions one by one to decipher the similarities and differences and to explain the large dispersion in the relations between the escape fraction of ionizing photons and the residual flux of the absorption lines. In particular, we apply a dust correction to the residual fluxes using observable quantities, in order to improve the correlations.

5.1. Distribution of the sources

In this section we assess whether the ionizing photons and the Lyβ or C II λ1334 photons are emitted from the same place in order to decipher whether or not they probe the same regions of gas. From stellar models, we know that only stars younger than ∼10 Myr are bright at < 910 Å (i.e., ionizing energies), whereas older stars can still be bright at 1020 Å and 1330 Å, corresponding to the continuum next to the absorption lines. This is the main driver for the difference between sources of ionizing photons and nonionizing photons. The first row of Fig. 14 illustrates those differences. It shows surface brightness maps for one particular direction of observation of the output C, at different wavelengths. Compared to the 910 Å map, we see that at 1020 Å there are more stars shining in an extended disk around the central part of the galaxy, and even more so at 1330 Å.

Fig. 14. Surface brightness maps of output C (z = 3.0). Upper left: intrinsic surface brightness at 910 Å. Upper middle: intrinsic surface brightness at 1020 Å. Upper right: intrinsic surface brightness at 1330 Å. Lower left: surface brightness at 910 Å after H0 and dust absorption. Lower middle: surface brightness at 1020 Å after dust extinction. Lower right: surface brightness at 1330 Å after dust extinction.

To quantify those differences, we note that the brightest 11% of the stellar particles emit 99% of the 910 Å photons, while they emit around 98% of the 1020 Å photons and 85% of the 1330 Å photons. More detailed results for the three outputs are provided in the first two columns of Table 2. While the sources of 1020 Å photons are almost the same as the 910 Å photon sources, the distribution of sources of 1330 Å photons is noticeably different from that of 910 Å photons. Around 15% of the luminosity at 1330 Å is emitted by stars that do not contribute to the ionizing photon budget.

However, it is the dust-attenuated 1020 Å or 1330 Å photons that are directly observable, and not the intrinsic luminosities. Here, we refer to the UV nonionizing photons after dust extinction as the “visible” photons. Once dust is added, the most luminous stars, namely the young stars in dense, neutral regions, are strongly attenuated. As those stars are the main sources of ionizing photons, the distribution of sources of visible photons becomes very different from the distribution of ionizing radiation sources. The middle and right panels of the second row of Fig. 14 illustrate the effect of dust on 1020 Å and 1330 Å photons. Various bright regions are strongly attenuated, leading to a substantially different distribution of sources than for 910 Å photons. The bottom left plot of Fig. 14 shows the 910 Å surface brightness after H⁰ and dust absorption for illustration.

We now compute the fraction of visible 1020 Å and 1330 Å photons emitted by the stars responsible for 99% of ionizing photon emission to assess whether or not the visible photons are tracing the gas screens through which ionizing photons travel. The results are listed in the last two columns of Table 2. When the attenuation by dust is included, the result depends on the direction of observation, and so we present the mean, and the 10th and 90th percentile of the results. In summary, we find that around 15% –but sometimes up to more than 30%– of visible 1020 Å photons are emitted by stars that do not contribute to the budget of ionizing photons, and therefore those photons are affected by screens of gas that do not influence the escape fraction of ionizing photons. For visible 1330 Å photons, the effect is even more pronounced, with around 30% and up to 60% of photons passing through gas that does not screen LyC sources. This is one reason why the residual flux of absorption lines before dust correction of the continuum cannot be a good proxy of the escape of ionizing photons. It is also one motivation to introduce dust correction of the UV nonionizing continuum, as is often done in the literature (e.g., Steidel et al. 2018; Gazagnes et al. 2018; Chisholm et al. 2018).

5.2. Dust correction

In addition to what we see in Sect. 5.1, another motivation to apply a dust correction is the different effects that dust has on both ionizing photons and on the nonionizing UV continuum in our galaxy. We highlight in Sect. 4.2 that dust absorbs very few ionizing photons because the optical depth of neutral hydrogen is always much larger than the optical depth of dust. However, the flux of the continuum next to the absorption line, which enters in the definition of the residual flux of the line, significantly depends on dust. This suggests that a better proxy of the escape fraction of ionizing photons is the ratio of the flux at the bottom of the line to the flux of the intrinsic continuum, rather than the residual flux we plot in Fig. 13, which uses the continuum after dust extinction. To obtain this new ratio, that we refer to here as the dust corrected residual flux, we multiply the residual flux (see Eq. (5)) by the dust extinction factor of the continuum:

(6)

where is the value of the continuum next to the absorption line before dust attenuation. Fortunately this dust extinction factor can be estimated by SED fitting, and therefore the dust correction is feasible in practice. The extinction factor is usually parametrized as 10^{−0.4 E_B − V k_λ}, where E_B − V models the strength of the dust attenuation and k_λ is a law of dust attenuation as a function of wavelength. However, the extinction factor can be obtained directly from the simulation after the transfer of the stellar continuum with RASCAS, which is why we do not use the observational method to apply the dust correction. This dust correction is very similar to the ones in Steidel et al. (2018), Gazagnes et al. (2018) and Chisholm et al. (2018). The only difference is that we use the extinction factor at the wavelength of the continuum next to the absorption lines, while these latter authors use the extinction factor at 912 Å. This is conceptually different but gives similar results.

Figure 15 shows the new relations between the escape fraction of ionizing photons and the dust-corrected residual flux of C II λ1334, on the left, and Lyβ, on the right. The results are more promising than before dust correction, in particular for C II λ1334. The corrected residual flux still shows significant scatter but now follows the one-to-one relation more closely, especially for directions with high escape fraction. Very low corrected residual fluxes still indicate very low escape fractions. For directions with dust-corrected residual fluxes of between 0.02 and 0.3, the prediction of f_esc using the residual flux is unfortunately often too high. Only around 18% of those directions have an escape fraction that is at least equal to 80% of the corrected residual flux. Above a corrected residual flux of 0.3, the predictions are more successful, with around 80% of the directions having an escape fraction at least equal to 80% of the corrected residual flux.

Fig. 15. Scatter plot of the escape fraction of ionizing photons and the dust-corrected residual flux for C II λ1334 on the left and Lyβ on the right. The residual fluxes are multiplied by the ratio between the flux of the continuum with dust extinction and the flux of the continuum without dust extinction.

After dust correction, Lyβ still traces f_esc more poorly than C II λ1334 does, but we see that the dust-corrected residual flux of Lyβ now provides a good lower limit for f_esc, whereas before dust extinction, in the right panel of Fig. 13, the points were extremely dispersed.

5.3. Covering fractions and escape fractions

We now turn to the question of the distribution of the column densities and the concept of covering fractions. At the beginning of Sect. 5 we explain that in order to have equality between the residual flux of absorption lines and the escape fraction of ionizing photons, the screen of gas must be split in two parts: some optically very thick regions and some empty holes. The escape fraction is then equal to one minus the covering fraction f_C, which is the fraction of sources that are covered by optically thick gas. If there is gas with a column density such that τ₉₁₀ ≈ 1, the covering fraction becomes harder to define and it is no longer true that f_esc = 1 − f_C.

In the first part of this section, we study the distribution of neutral hydrogen column densities that the ionizing photon sources are facing. We then define the covering fraction in the simulated galaxy and see how well it traces f_esc in the same three outputs and 1728 directions of observation as in the previous sections.

An additional difficulty arises when using absorption lines to probe this covering fraction: Lyβ does not have the same optical depth as ionizing photons for a given column density of H⁰ and C⁺ is not a perfect tracer of H⁰, which is why C II λ1334 also has a distribution of optical depths that is different from that of ionizing photons. We discuss these issues in Sects. 5.3.2 and 5.3.3.

5.3.1. Covering of ionizing radiation sources

To study the distribution of column densities in a given direction of observation we first compute the column density in front of every stellar particle. As we want to study the coverage of photons rather than the coverage of stars, we give a weight to each stellar particle equal to its luminosity at 910 Å. We then compute the fraction of photons behind a given amount of column density. A distribution perfectly matching the picket-fence model would show two distinct populations (for example, 30% of the ionizing photons facing less than N_H I = 10¹⁵ cm⁻² and 70% facing more than N_H I = 10¹⁹ cm⁻²).

We show in Fig. 16 distributions of column densities from three of the seven directions of observation that we studied in Sect. 3. The pink histogram is typical of the majority of our directions, where the escape fraction is zero. All the photons are emitted behind optically thick gas, with N_H I > 10²⁰ cm⁻². As there are no photons with τ₉₁₀ ≈ 1, the covering fraction is well defined: it is equal to one and the escape fraction is zero. The orange histogram shows a clear bimodality, with a (double)-peak of the curve around N_H I = 10¹⁵ − 10¹⁶ cm⁻², which is the “holes” part, and another peak at high column densities, which is the “cloud” part. There is only a tiny fraction of 910 Å photons that are facing gas with τ₉₁₀ ≈ 1, and so this is also a case where the picket-fence model is a good description. However, the dark-blue histogram is an example of a situation where there is a significant fraction of photons facing gas such that τ₉₁₀ ≈ 1. The photons are partially absorbed in this regime, which makes the notion of the covering fraction unclear. The cyan dashed line corresponds to the average of all 5184 directions, and it shows that the most common configurations have almost exclusively high column densities, and thus very low escape fractions.

Fig. 16. Probability distribution of column densities faced by 910 Å photons. The three filled histograms correspond to the three directions that have the same color in Sect. 3. The cyan dashed line is the average distribution over the 1728 directions of the three outputs. It has a covering fraction of 93%. Molecular hydrogen chemistry is not considered in the simulation, and so the real atomic hydrogen column densities should saturate around 1023 cm−2.

Let us define the covering fraction as the fraction of 910 Å photons facing more than N_H I = 10^17.2 cm⁻², which is the column density for which τ₉₁₀ = 1:

(7)

where is the column density of neutral hydrogen in front of a stellar particle and is the luminosity of the stellar particle at 910 Å. We note that depends on the direction of observation of the galaxy. In an idealized picket-fence geometry with empty holes and optically thick clouds, we have the relation . In the general case, photons that are covered in the sense of Eq. (7) can still have a transmission factor of up to e⁻¹ = 37%. Indeed, if all the gas has exactly N_H I = 10^17.2, the photons are considered as covered, but the optical depth in this case is τ₉₁₀ = 1, and so 37% of photons can still escape. Similarly, photons that are not covered can still be absorbed; we can only say that their transmission factor is larger than 37%. The equation for the allowed region in the relation is:

(8)

Figure 17 displays the relation between f_esc and for our 5184 mock observations and shows that indeed not all the directions are perfectly described by the picket-fence model, in the sense that we do not always have . We note that a few points are below the relation given by Eq. (8) because of absorption by helium and dust that was not considered in this formula. The directions of observations corresponding to the pink and the orange histograms in Fig. 16 are close to the one-to-one relation in Fig. 17, as in the picket-fence geometry. This was expected because those directions have almost no stars facing gas such that τ₉₁₀ ≈ 1. However, the dark-blue dot sits further from the one-to-one relation, which is explained by the significant fraction of photons facing gas with τ₉₁₀ ≈ 1 in the corresponding histogram of Fig. 16. We see that overall the directions are close to the one-to-one relation. Therefore, the picket-fence geometry is on average a reasonable model of the covering of ionizing photons in our galaxy.

Fig. 17. Relation between the covering fraction of ionizing photons, defined by Eq. (7), and the escape fraction of ionizing photons for the 1728 directions of observation of the three outputs, A, B, and C. The seven colored dots are the directions that we study in Sect. 3. The pink area highlights the region of the figure that is allowed by our definition of covering fraction.

Here we considered the description of the escape of ionizing photons in the picket-fence framework and compared to the simulation, but the goal is to use absorption lines to infer this escape fraction. This leads to additional limitations due to the different optical depths of Lyβ and C II λ1334 compared to LyC, and the difference in distribution of C⁺ compared to H⁰. We now explain those new limitations for the two lines.

5.3.2. The case of Lyβ

In order for the residual flux of Lyβ to be an indicator of the covering fraction, the “holes” in the picket-fence have to be optically thin for both LyC and Lyβ. However the column density at which τ_Lyβ = 1, with a turbulent velocity of 20 km s⁻¹, is N_H I = 10^14.22 cm⁻², which is three orders of magnitude below the column density at which τ_LyC = 1. For example, the orange distribution of Fig. 16 represents a case where the picket-fence model works well to describe the escape of ionizing photons, in the sense that , but Lyβ residual flux is not a good indicator of this f_esc because the part that is optically thin for LyC has N_H I of around 10¹⁵ − 10¹⁶ cm⁻², which is already optically thick for Lyβ. This is confirmed in the right panel of Fig. 15, where the orange dot has a residual flux of only 7%, while the escape fraction is as high as 47%. The difficulty of the optical depth of Lyβ being larger than that of LyC cannot be circumvented. A saturated Lyβ line can very well correspond to a direction with a high escape fraction of ionizing photons. This is why Lyβ is a good lower limit of f_esc, but not a direct tracer.

5.3.3. The case of C_IIλ1334

We now analyze the optical depth of C II λ1334 in comparison with the one of 910 Å photons. The column density of C⁺ at which the optical depth of the center of C II λ1334 is one for a turbulent velocity of 20 km s⁻¹ is N_C II = 10^13.89 cm⁻². For 910 Å photons, it is N_H I = 10^17.19 cm⁻², which means that if N_C II/N_H I = 5 × 10⁻⁴, the two optical depths are the same. This ratio is close to the ratio of carbon to hydrogen in the Sun, but metallicity and ionization effects can shift N_C II/N_H I further from this value.

In Fig. 18, we show the ratio of carbon to hydrogen column densities in the upper panel and of C⁺ over H⁰ column densities in the lower panel for one direction of our simulated galaxy. We see that the C⁺ over H⁰ ratio varies by about eight orders of magnitude from region to region. However, the carbon-to-hydrogen ratio is more uniform, with an average of 1.5 × 10⁻⁴, and varying by less than two orders of magnitude. Thus, the variation of the C⁺-to-H⁰ ratio comes primarily from ionization effects. There are regions abundant in H⁰ but lacking C⁺, because most carbon is in the neutral state. In other regions, C⁺ is overabundant compared to the N_C/N_H ratio, which occurs when hydrogen is mostly ionized and carbon is mainly in the C⁺ state, which is possible because the ionization energy of C⁺ is higher than that of H⁰.

Fig. 18. Ratio of column densities in the same output and direction as in Figs. 1 and 14. Upper panel: ratio of carbon to hydrogen column densities. Lower panel: ratio of C+ to H0 column densities. The column densities are computed from the stellar particles to the virial radius. In every pixel, the ratio is averaged over all stellar particles in this pixel and weighted by the luminosity of the stellar particle at 1334 Å. The black pixels have no stellar particles in them.

For the continuum of stars behind screens with N_C II/N_H I < 5 × 10⁻⁴, the C II λ1334 absorption is weaker than the absorption of ionizing photons, which leads to the residual flux of C II λ1334 underestimating f_esc. Conversely, when stars are facing gas such that N_C II/N_H I > 5 × 10⁻⁴, the residual flux of C II λ1334 overestimates f_esc. Therefore, the variation of the value of N_C II/N_H I in the galaxy causes inevitable dispersion in the relation between the residual flux of C II λ1334 and the escape fraction of ionizing photons. This dispersion changes with time and orientation, depending on the alignment between the brightest stars and the screens of gas with different N_C II/N_H I ratios.

5.4. Processes affecting the residual flux

Lastly, we highlight the extent to which processes during the radiative transfer of the line affect its properties; in particular the residual flux. We show in Sect. 3.3 that the absorption lines are sensitive to the velocity distribution of the gas and sometimes to the effect of infilling, especially for directions with strong fluorescence. However, correcting for scattering and gas velocity dispersion does not improve the tightness of the relation between the residual flux of C II λ1334 and the escape fraction of ionizing photons. Figure 19 shows this relation with the scattering ignored and the gas velocity set to zero. We see that the relation has as much dispersion as in Fig. 15. Moreover, the correction of scattering and gas velocity brings the points further away from the one-to-one relation, leading to C II λ1334 being a lower limit of f_esc, similar to Lyβ. In other words, C II λ1334 globally saturates at column densities of gas that are too low for ionizing photons to be significantly absorbed, but the effect of infilling and the dispersion of gas velocity increase the residual flux of C II λ1334 so that it gets closer to the value of f_esc.

Fig. 19. As in the left panel of Fig. 15, but omitting scattering and setting gas velocity to zero.

5.5. Another metallic line: Si_IIλ1260

We find that C II λ1334 traces the escape fractions of ionizing photons better than Lyβ does. We now show that all the results for C II λ1334 apply to Si II λ1260, another LIS absorption line. The modeling of the distribution of Si⁺ is slightly less robust than that of C⁺, first because silicon depletes more onto dust, which we do not model here, and second because of charge transfer reactions (Appendix A). However, the spectra are not very sensitive to a change in density, as shown in Sect. 2.2. First, Fig. 20 shows the seven directions of observations that we selected in this paper, with the spectra of Si II λ1260 on the left and that of C II λ1334 on the right. We see that they have similar properties: the residual flux and EWs are almost the same. This is confirmed in Fig. 21 where we plot the relation between the escape fractions and the different properties of the Si II λ1260 line. All the relations are similar to those seen in the C II plot in the left panel of Fig. 13. Finally, as Fig. 22 shows, the relations after dust correction are also closer to the one-to-one relation. Si II λ1260 is therefore as good a tracer of f_esc as C II λ1334, providing relatively good predictions when the (dust corrected) residual flux is large, but unfortunately predictions that often overestimate reality when the residual flux is smaller than 0.3.

Fig. 20. Same as Fig. 3 but with Si II λ1260 instead of Lyβ. The vertical dashed lines highlight the wavelengths of the fluorescent channels.

6. Summary

In this paper we present a method to post-process radiative hydrodynamics simulations in order to build mock observations of LIS UV absorption lines. We applied this method to a cosmological zoom-in simulation of a galaxy at z = 3 with a stellar mass of 2.3 × 10⁹ M_⊙ in order to study the processes that affect the formation of absorption lines and to find correlations between the escape fraction of ionizing photons and the properties of the absorption lines. Our main results can be summarized as follows:

1. Our procedure reproduces metal absorption lines with a diversity of shape and strength (Fig. 3) that is similar to observations; such as for example those of Jaskot et al. (2019). The Lyβ absorption lines of our galaxy seem to differ from the ones in Steidel et al. (2018) or Gazagnes et al. (2020), probably because our simulated galaxy is less massive and has a lower SFR.

2. Many properties of the observed radiation depend highly on the direction of observation. There are variations of up to more than one magnitude at 1500 Å (Table 1). The LyC escape fractions vary from 0% to 47% (Fig. 11). The residual flux goes from 0% to 95% (Fig. 13).

3. The precise implementation of the radiation between 6 eV and 13.6 eV does not impact the mock absorption lines. There are dense regions where the carbon ionization fraction depends on this implementation, but they are very dusty and do not contribute to our spectra. However, using a UV background at all energies instead of the ionizing radiation from the simulation fails to reproduce the same observed spectra (Fig. 4). This shows that radiative transfer of ionizing photons in the simulation is crucial for this kind of study.

4. The fluorescent lines C II⋆ λ1335 and Si II⋆ λ1265 are clearly visible (Fig. 20), with an EW of the same order as for the absorption. The shape and size of the fluorescent lines is sensitive to the modeling of the fine-level structure of the ground state of their respective ions (Fig. 6). The presence of absorption by ions in the fine-level structure can lead to the formation of a P-Cygni profile in the fluorescence.

5. The infilling effect increases the residual fluxes by around 0.06 on average. There are however directions where the effect is much stronger, with a change of up to more than 0.3 (Fig. 7). A fluorescent line with a large EW is an indicator that the infilling effect is stronger than on average (Fig. 8). Infilling depends very little on the aperture size of the mock observation, as long as it encompasses more than ≈90% of the stellar continuum (Figs. 6 and 9).

6. The large-scale distribution of the gas velocity spreads the absorption over a larger wavelength range than if all the gas had the same velocity (Figs. 9 and 10). This spread leads to an increase in the residual flux of the lines (Fig. 7), which is one factor that challenges the use of absorption lines as tracers of the escape fraction of ionizing photons.

7. The properties of the absorption lines without dust correction correlate poorly with the LyC escape fraction (Figs. 13 and 21). The residual fluxes of C II λ1334 or Si II λ1260 cannot predict the escape fraction, but are good upper limits. A large EW of absorption, either of a metallic line or of Lyβ, or a large EW of fluorescence, indicate an escape fraction smaller than 2%.

8. Multiplying the residual flux of C II λ1334 or Si II λ1260 by the extinction of the continuum by dust makes it a better indicator of escape fraction, following globally the one-to-one relation, although with strong dispersion (Figs. 15 and 22). This dust correction is similar to the ones made in the literature (e.g., Steidel et al. 2018; Gazagnes et al. 2018; Chisholm et al. 2018). For directions with dust-corrected residual flux between ∼2% and 30%, the escape fraction is often overestimated. Concerning Lyβ, multiplying its residual flux by the extinction by dust makes it a good lower limit of the escape fraction, but not a direct proxy (Fig. 15).

9. The covering fraction, defined as the fraction of photons that are facing N_H I ≥ 10^17.2 cm⁻², gives a reasonable estimate for the escape fraction of ionizing radiation for our simulated galaxies, although the scatter in the relation is not negligible (Fig. 17).

10. Many geometrical and physical complexities make it difficult to use absorption lines as reliable tracers of the escape fraction. The different distributions of sources for ionizing photons and nonionizing UV continuum (Fig. 14), the complex distribution of column densities in front of stars (Fig. 16), the different optical depth between the ionizing photons and the lines, the tangled distribution of C⁺ with respect to H⁰ (Fig. 18), and the effects of infilling and gas velocity distribution are not correctable and explain the dispersion in the scatter plots, and the fact that dust-corrected residual fluxes below ∼30% overestimate the escape fraction.

Fig. 21. Same as Fig. 13 but for Si II λ1260.

Fig. 22. Same as Fig. 15 but for Si II λ1260.

In summary, metallic absorption lines can predict the escape fraction of ionizing photons accurately only in two extreme regimes: when f_esc = 0, with saturated absorption lines, or when the escape fraction is very high, with a dust corrected residual flux larger than 30%. As a result, even if it may be difficult to estimate the escape fraction of ionizing radiation from a galaxy, metallic absorption lines have the potential to rule out LyC emission from line-saturated galaxies, and to select promising candidates with high residual fluxes.

We are happy to share the mock observations of our simulated galaxy and encourage those interested to contact us.

Acknowledgments

The authors are grateful to the anonymous referee for her/his very careful and constructive report. We also thank J. Kerutt for her detailed revision of the manuscript. VM, AV, and TG are supported by the ERC Starting Grant 757258 “TRIPLE”. TK was supported in part by the National Research Foundation of Korea (NRF-2019K2A9A1A0609137711 and NRF-2020R1C1C100707911) and in part by the Yonsei University Future-leading Research Initiative (RMS2-2019-22-0216). This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. The simulation could be run thanks to our access to computing resources at the Common Computing Facility (CCF) of the LABEX Lyon Institute of Origins (ANR-10-LABX-0066).

References

Appendix A: Comparison of ionization fractions between KROME and Cloudy

In Fig. A.1 we show a comparison of the ionization fractions of carbon and silicon as a function of temperature. We include here only recombinations and collisional ionizations, and no photo-reactions. The chemical network of KROME is chosen to match the results of Cloudy. We see that the results for carbon are coherent between the two codes. For silicon, the effects of charge transfer, which is included in Cloudy but not in KROME, are more important, and we do not recover the curves of Cloudy as well as with carbon. We still choose to omit those reactions in order to preserve the state of hydrogen and helium, as explained in Sect. 2.2.

Fig. A.1. Ionization fractions of carbon and silicon as a function of temperature. KROME reactions are recombinations from Badnell (2006) and collisional ionization from Voronov (1997).

Appendix B: Computing photoionization rates

In KROME, the usual method to include photo-reactions is in three steps:

• Specify the photo-reaction in the chemical network, for example C⁺ + γ → C⁺⁺ + e.

• Provide a file with a cross-section as a function of photon energy for this reaction.

• Before calling KROME in a cell of a simulation (on-the-fly or in post-processing), specify a discretized spectral energy distribution.

KROME then computes the rate of the photo-reaction based on the cross-section and the radiation field, and uses it to compute the chemical evolution. However, the photo-ionization rates of metallic ions can be easily computed from the RAMSES-RT output (along with the modeling of the radiation below 13.6 eV, see Sect. 2.2) without using KROME, and so we bypass the method above, giving the more simple procedure:

• Specify the photo-reaction in the chemical network.

• For every cell, compute the photoionization rate with the method below.

• Give this photoionization rate to KROME via a routine that we added in the KROME code.

The photoionization rates are products of cross-sections and fluxes of photons integrated over frequency. The flux is given in every cell directly by the simulation. For cross-sections, we have to apply the same method as RAMSES-RT, where the photoionization cross-sections are weighted by the SEDs of all stellar particles in the simulation.

To begin with, we have a stellar SED library, BPASS in our case, with an SED for a grid of stellar ages and metallicities. For all SEDs, at each stellar age and metallicity in the grid, we compute the mean cross-section of photoionization for each bin of radiation of RAMSES-RT, weighted by the number of photons, with this formula:

(B.1)

where bin i is the ith bin of radiation of the simulation, X indicates an ion for which we compute the rate, is the luminosity as a function of frequency for a given stellar age and metallicity in BPASS, in units of erg s⁻¹ Hz⁻¹, hν is the energy of a photon at frequency ν, and is the photoionization cross-section as a function of frequency for species X, given by Verner et al. (1996).

From those mean cross-sections we compute the cross-sections of a given stellar particle, , by interpolating on age and metallicity. In a given domain of the simulation, where we want to compute the photoionization rates in every cell, we average those cross-sections on all stellar particles, weighing by their luminosities:

(B.2)

where is the luminosity of the stellar particle j_star in the energy bin i. Finally, the photoionization rate in a given cell is:

(B.3)

where is the flux of photons in the cell, in units of s⁻¹ cm⁻², for the ith bin of radiation.

All Tables

All Figures

Fig. 1. Maps of neutral hydrogen column density and surface brightness at 1500 Å for a snapshot of the simulated galaxy at z = 3. Upper left: column density of neutral hydrogen atoms, with the virial radius highlighted by the white circle. Upper right: intrinsic surface brightness at 1500 Å. Lower left: zoom-in to a tenth of the virial radius. Lower right: surface brightness at 1500 Å after extinction by dust.

In the text

Fig. 2. Energy levels of the C+ ion. P14 = 100% indicates that a C+ ion in level 1 can be photo-excited only to level 4. The green numbers are the probabilities that a C+ ion in level 4 radiatively de-excites to level 1 or 2. The red numbers are the probabilities that an incoming photon with a wavelength between 1335.66 Å and 1335.71 Å hitting a C+ ion in level 2 excites it to level 3 or 4. P32 = 100% indicates that a C+ ion in level 3 can only radiatively de-excite to level 2.

In the text

Fig. 3. C II λ1334 and Lyβ absorption lines for seven directions of observation. The vertical dashed lines highlight the wavelength of the fluorescent channels. The spectra are arranged in increasing order of escape fraction of ionizing photons in the direction of observation, fesc (as computed in Sect. 2.6). The noise is Monte-Carlo noise, due to a sampling of the luminosity with a finite number of photon packets. The spectra are smoothed with a Gaussian kernel with FWHM = 20 km s−1. The aperture diameter for the mock observation is a fifth of the virial radius, which corresponds to about 5.6 kpc, or 0.75″ at z = 3. No observational effects such as noise, spectral binning, Earth atmosphere, or MW alterations are considered. We note that the x-axis in the first two rows of the second column differ from the others.

In the text

Fig. 4. Test of the robustness of the spectra against changing free parameters. The seven rows show the same seven directions as in Fig. 3. The vertical dashed lines highlight the wavelength of the fluorescent channels. First column: test of the effect of multiplying and dividing the fiducial carbon density by two. Second column: test of multiplying (dividing) the Habing band UVB by 100 (2) and comparison with the result with only UVB and no radiation field from the simulation. Third and fourth columns: test of changing the turbulent velocity for C II λ1334 and Lyβ respectively. The fiducial turbulent velocity is 20 km s−1.

In the text

Fig. 5. Effect of different parameters of the modeling on the residual fluxes. The y-axes display differences of residual fluxes between two setups. The dots are the median of the differences in a given bin of residual flux, and the bars show the 10th and 90th percentiles. The numbers indicate how many spectra are in the corresponding bin of residual flux. Left panel: effects of carbon density, which is divided or multiplied by two. Middle panel: effect of UVB radiation for ionizing photons instead of the radiation field from the simulation. Right panel: effects of turbulent velocity. Turbulent velocities are varied between 0 and 40 km s−1.

In the text

Fig. 6. Illustration of the effects of different physical processes on the same spectra as Fig. 3. The vertical dashed lines highlight the wavelength of the fluorescent channels. First column: comparison with and without a fraction of C+ ions populating their ground-state fine-structure level. Second column: effect of removing the scattering of photons. Third column: effect of changing the aperture of the observation.

In the text

	Fig. 7. Same as Fig. 5, to show the effects of different physical processes. Left, middle, and right panels: effects of scattering, dust, and gas velocity, respectively.
In the text

	Fig. 8. Same as left panel of Fig. 7, in bins of EW of the fluorescent line instead of residual flux. The best fit is 0.06 for EWfluo <  0.6 Å and 0.11 EWfluo − 0.003 otherwise.
In the text

	Fig. 9. Same as Figs. 6 and 10 but for Lyβ. No fine-structure level is considered here.
In the text

Fig. 10. Same as Fig. 6, with other processes. First column: effect of removing dust. We also remove the fine-structure level to simplify the interpretation. Second column: spectra with and without gas motion in the galaxy. Infilling effects give complex shapes to the second and seventh rows. Third column: comparison of fiducial spectra with spectra with neither gas motion, dust, scattering, nor fine-structure level.

In the text

	Fig. 11. Histogram of LyC escape fractions for the three outputs, computed as in Sect. 2.6.
In the text

	Fig. 12. Comparison of the escape fraction at 900 Å and the global escape fraction of ionizing photons.
In the text

Fig. 13. Scatter plots comparing the escape fractions for the three outputs in the 1728 directions of observation to the properties of the C II absorption line on the left and the Lyβ line on the right. The seven colored dots show the position of the spectra displayed in Sect. 3. EW_abs is the EW of the absorption line (negative EW_abs means more emission than absorption), EW_fluo is the EW of the fluorescent emission, vmax is the velocity of the deepest point in the spectrum, vcen is the velocity at which the absorption line has half its EW, and v90 is the velocity where the blue side of the absorption line reaches 90% of the continuum value. The orange vertical lines show the median value of the x-axis variable over the 5184 directions. There is no fluorescent channel for Lyβ, and so one panel is missing.

In the text

Fig. 14. Surface brightness maps of output C (z = 3.0). Upper left: intrinsic surface brightness at 910 Å. Upper middle: intrinsic surface brightness at 1020 Å. Upper right: intrinsic surface brightness at 1330 Å. Lower left: surface brightness at 910 Å after H0 and dust absorption. Lower middle: surface brightness at 1020 Å after dust extinction. Lower right: surface brightness at 1330 Å after dust extinction.

In the text

	Fig. 15. Scatter plot of the escape fraction of ionizing photons and the dust-corrected residual flux for C II λ1334 on the left and Lyβ on the right. The residual fluxes are multiplied by the ratio between the flux of the continuum with dust extinction and the flux of the continuum without dust extinction.
In the text

Fig. 16. Probability distribution of column densities faced by 910 Å photons. The three filled histograms correspond to the three directions that have the same color in Sect. 3. The cyan dashed line is the average distribution over the 1728 directions of the three outputs. It has a covering fraction of 93%. Molecular hydrogen chemistry is not considered in the simulation, and so the real atomic hydrogen column densities should saturate around 1023 cm−2.

In the text

	Fig. 17. Relation between the covering fraction of ionizing photons, defined by Eq. (7), and the escape fraction of ionizing photons for the 1728 directions of observation of the three outputs, A, B, and C. The seven colored dots are the directions that we study in Sect. 3. The pink area highlights the region of the figure that is allowed by our definition of covering fraction.
In the text

Fig. 18. Ratio of column densities in the same output and direction as in Figs. 1 and 14. Upper panel: ratio of carbon to hydrogen column densities. Lower panel: ratio of C+ to H0 column densities. The column densities are computed from the stellar particles to the virial radius. In every pixel, the ratio is averaged over all stellar particles in this pixel and weighted by the luminosity of the stellar particle at 1334 Å. The black pixels have no stellar particles in them.

In the text

	Fig. 19. As in the left panel of Fig. 15, but omitting scattering and setting gas velocity to zero.
In the text

	Fig. 20. Same as Fig. 3 but with Si II λ1260 instead of Lyβ. The vertical dashed lines highlight the wavelengths of the fluorescent channels.
In the text

	Fig. 21. Same as Fig. 13 but for Si II λ1260.
In the text

	Fig. 22. Same as Fig. 15 but for Si II λ1260.
In the text

	Fig. A.1. Ionization fractions of carbon and silicon as a function of temperature. KROME reactions are recombinations from Badnell (2006) and collisional ionization from Voronov (1997).
In the text