The circumgalactic medium in Lyman-alpha: a new constraint on galactic outflow models

TL;DR

This study compares Lyα observations of the circumgalactic medium around Lyman Break Galaxies with simulations to better understand galactic outflows and their impact on galaxy evolution.

Contribution

It introduces a new Monte-Carlo radiative transfer code and demonstrates how Lyα observations can constrain galactic outflow models.

Findings

Fiducial outflow model better matches observed CGM Lyα features.

Outflows influence Lyα signals mainly within 50 kpc of galaxies.

Cold accretion flows dominate Lyα signals beyond 50 kpc.

Abstract

Galactic outflows are critical to our understanding of galaxy formation and evolution. However the details of the underlying feedback process remain unclear. We compare Ly$\alpha$ observations of the circumgalactic medium (CGM) of Lyman Break Galaxies (LBGs) with mock observations of their simulated CGM. We use cosmological hydrodynamical `zoom-in' simulations of an LBG which contains strong, momentum-driven galactic outflows. Simulation snapshots at $z=2.2$ and $z=2.65$ are used, corresponding to the available observational data. The simulation is post-processed with the radiative transfer code \textsc{crash} to account for the impact of ionising photons on hydrogen gas surrounding the simulated LBG. We generate mock absorption line maps for comparison with data derived from observed close galaxy-galaxy pairs. We perform calculations of Ly$\alpha$ photons scattering through the CGM with our newly developed Monte-Carlo code \textsc{slaf}, and compare to observations of diffuse Ly$\alpha$ halos around LBGs. Our fiducial galactic outflow model comes closer to reproducing currently observed characteristics of the CGM in Ly$\alpha$ than a reference inefficient feedback model used for comparison. Nevertheless, our fiducial model still struggles to reproduce the observed data of the inner CGM (at impact parameter $b<30$kpc). Our results suggest that galactic outflows affect Ly$\alpha$ absorption and emission around galaxies mostly at impact parameters $b<50$ kpc, while cold accretion flows dominate at larger distances. We discuss the implications of this result, and underline the potential constraining power of CGM observations - in emission and absorption - on galactic outflow models.

Full text

††pagerange: The circumgalactic medium in Lyman-$\alpha$: a new constraint on galactic outflow models–References††pubyear: 2017

1 Introduction

Understanding processes which govern galaxy formation and evolution is one of the challenges in modern cosmology. In the standard $\Lambda$CDM cosmological scenario strong feedback is needed at both the low and high-mass ends of the galaxy mass function in order for models to match observations (White & Frenk, 1991). ‘Feedback’ typically refers to the complex processes through which star formation and accretion onto black holes deposit energy and momentum back into their surroundings. The details of feedback are not well understood (see Ciardi & Ferrara 2005 and Naab & Ostriker 2016 for reviews on the topic). Because of its importance though, it is fundamental to study as many (potential) observational probes of this process as possible.

There is an increasing amount of data on the so-called ‘circum-galactic’ medium (CGM), which has been defined as the region around galaxies out to a distance of $r\sim 300$ kpc and with a velocity offset from the galaxy’s systemic redshift of up to $\Delta v\sim 300$ km s*-1* (Rudie et al., 2012). Steidel et al. (2010) (hereafter S2010) note that the CGM provides a ‘laboratory in which the effects of galaxy formation and AGN accretion (e.g., radiative and hydrodynamical feedback and its recent history) can be measured on scales that are not accessible using direct observations of galaxies. Indeed analyses of the COS-Halos Survey (Werk et al. 2014; Prochaska et al. 2017) show that the CGM is a dominant baryon reservoir on a galactic scale, and thus reinforces the notion that details of its composition, state, and dynamics may provide strong constraints on galaxy evolution models.

In this paper, we explore whether observations of the CGM in Ly$\alpha$ can constrain galactic outflow models, by comparing simulations of the CGM with observations.

Steidel et al. (2011) (hereafter S2011) found spatially extended Ly$\alpha$ emission from the CGM of Lyman Break Galaxies (LBGs) by stacking Ly$\alpha$ observations of 92 individual galaxies. Similar stacking analyses have revealed (fainter/smaller) Ly$\alpha$ halos around Ly$\alpha$ selected galaxies (i.e. Lyman Alpha Emitters, LAEs) (Matsuda et al., 2012; Momose et al., 2014; Xue et al., 2017)111Feldmeier et al. (2013) did not find Ly$\alpha$ halos around $z\sim 2$ LAEs (also see Smith et al., 2012), and discuss that systematic uncertainties associated with stacking could reduce the statistical significance of previously reported detections. Momose et al. (2014) used larger samples of LAEs to confirm these systematic effects, but still obtained significant detections of Ly$\alpha$ halos.. At low redshift (0.028 $<z<$ 0.18), the Lyman alpha Reference Sample (Hayes et al., 2013) likewise shows extended Ly$\alpha$ emission relative to H$\alpha$ in imaging observations from the Hubble Space Telescope.

S2011 proposed that the diffuse halos arise from Ly$\alpha$ photons, produced in star-forming regions, scattering off outflowing material as they escape the galaxy. The presence of these outflows was inferred from the ubiquitous blueshifted low-ionisation absorption lines, and their interaction with Ly$\alpha$ photons was inferred from the redshifted Ly$\alpha$ emission lines (S2010). Here, star formation is both the source of the Ly$\alpha$ photons (produced by recombination in the HII regions around young stars) and of the stellar feedback which drives the outflowing material. S2011 provided a simple analytic model for the scattering of Ly$\alpha$ photons through the outflow, and showed that it is a good fit to their stacked observations. The radiative transfer of Ly$\alpha$ photons was treated though in an extremely approximate fashion.

Motivated by these models, Dijkstra & Kramer (2012) performed a systematic study of Ly$\alpha$ transfer through phenomenological models of spherically and cylindrically symmetric, large-scale, clumpy outflows. In these models, clumps were exclusively outflowing, with a one-to-one correspondence between outflow velocity ($v$) and distance from the galaxy ($r$). Following S2011, the velocity profile was inspired by ‘momentum-driven’ wind models in which the outflow accelerates as $a\propto r^{-\alpha}$ ($\alpha\sim 1.5$). Dijkstra & Kramer (2012) constrained the HI properties of their clumps by matching the galaxy-galaxy pair absorption line presented by S2010. While these models can simultaneously explain the presence of Ly$\alpha$ halos and the amount of absorption, they also predict that a non-negligible fraction of Ly$\alpha$ photons did not scatter at all. These photons should be visible as a bright point source, which is absent in the data. As Dijkstra & Kramer (2012) pointed out, this problem could be due to the absence of low-column density HI systems, and/or a consequence of the simplified velocity profile of the outflow in their models.

In this work, we take a new complementary approach, and use cosmological hydrodynamical ‘zoom’ simulations from Genel et al. (2012) to generate a model CGM of a LBG. These simulations contain strong galactic outflows which are also momentum-driven, and provide us with a complex CGM that may more closely reflect reality than previous models. Importantly, the simulations contain inflowing ‘cold streams’ (e.g. Kereš et al., 2005; Dekel et al., 2009), which can contribute significantly to the amount of absorption measured in the CGM (van de Voort et al., 2012; Goerdt et al., 2012) and possibly to the emission (e.g. Dijkstra & Loeb, 2009; Rosdahl & Blaizot, 2012). While these simulations do not have the resolution to properly resolve the feedback processes and the kinematics of the cold gas, it is important to check how they compare to the available data on the CGM.

While there are a number of previous works which use simulations to study observational signatures of the CGM, our work distinguishes itself by simultaneously considering the CGM in Ly$\alpha$ emission and absorption. Previous works that used simulations have focussed either on emission (e.g. Laursen et al., 2009; Faucher-Giguère et al., 2010; Barnes et al., 2011; Rosdahl & Blaizot, 2012) or absorption (e.g. Goerdt et al., 2012; Shen et al., 2013). Dijkstra & Kramer (2012) show that joint constraints from Ly$\alpha$ in absorption and emission are much more powerful than either data set individually. Observationally this is demonstrated by Hennawi & Prochaska (2013).

Finally, in contrast to previous studies that modelled the CGM in absorption, our simulations are post-processed with an ionising photon transfer code (crash, Ciardi et al. 2001) and account for local sources of ionising radiation, which can be more important than the overall ionising background (especially at close distances to the galaxy, see e.g. Shen et al., 2013).

This paper is laid out as follows: in § 2 we describe the simulations, § 3 presents the output of our pipeline and compares our results to observations, § 4 discusses the results and our conclusions222Throughout this work, we used the following cosmological parameters based on the 3-year WMAP results (Spergel et al., 2007): $\Omega_{\Lambda}=0.74$, $\Omega_{m}=0.26$, $\Omega_{b}=0.044$, $h=0.72$, $n=0.95$, and $\sigma_{8}=0.77$..

2 Simulations

2.1 Cosmological Hydrodynamical Simulations

It is not currently computationally feasible to perform full hydrodynamic cosmological simulations with sufficient resolution to resolve the detailed gas dynamics of the CGM. Therefore we use the cosmological ‘zoom-in’ simulations from Genel et al. (2012), which start with an N-body dark matter only simulation. A region of space is cut out around a massive dark matter halo and re-simulated, adding baryons and hydrodynamic physics using a modified version of gadget-2 (Springel 2005, Oppenheimer & Davé 2006, and Oppenheimer & Davé 2008).

Genel et al. (2012) use a modified Oppenheimer & Davé (2008) wind model, which implements momentum-driven winds powered by stellar feedback. In brief, gas particles become star-forming particles as per the criteria described in Schaye & Dalla Vecchia (2008), and then subsequently become wind particles which are stochastically kicked perpendicularly to the plane of the galaxy. The strength of the kick given to a particular wind particle is given by $v_{wind}=\sigma(4+4.29\sqrt{f_{L}-1})$, where $\sigma$ is the velocity dispersion of the galaxy, and $f_{L}$ is the luminosity factor stochastically chosen in the range [1.05-2]. With this model $v_{wind}$ is higher than that used in Oppenheimer & Davé (2008), resulting in typical wind velocities of $\approx$ (400-700) km s*-1* for the halos under consideration here. This, combined with the fact that in this model wind particles are temporarily decoupled from hydrodynamics so that they initially encounter no drag, ensures that the wind particles escape the disk. The mass-loading factor, which is the wind mass loss rate divided by the star formation rate, is typically $\eta\approx 4$.

It is worth noting that this wind model was used in a suite of simulations which reproduce the metallicity and ionisation of the intergalactic medium (IGM), the galaxy mass-metallicity relation, the high galactic gas fraction at high redshift ($z\gtrsim 2$), and the fact that galaxies contain a low fraction of cosmic baryons ($\sim 5-10\%$ at $z=0$ Fukugita & Peebles 2004). In other words, despite its simplicity, this model simultaneously reproduces several observational constraints and scaling relations.

The hydrodynamic simulation, covering a region of $\approx$ 5 Mpc comoving, has a mass resolution of $8\times 10^{5}M_{\odot}$ for baryonic particles, and $5\times 10^{6}M_{\odot}$ for dark matter particles. The gravitational softening length of the baryonic particles is 200$h^{-1}$ pc comoving. Since we ultimately want to compare to S2010 and S2011, we select galaxies with a similar stellar mass (i.e. $M_{*}\approx 10^{10.5}M_{\odot}$) and use the gadget snapshots at $z=2.2$ (absorption) and $z=2.65$ (emission), which are equivalent to the mean redshifts of the observations in S2010 and S2011, respectively. Further details of the specific galaxy under consideration in this paper (identified as $s396$) can be found in Table 1 of Genel et al. (2012). In brief, at the $z=2.2$ snapshot, $s396$ resides in a $1.5\times 10^{12}M_{\odot}$ dark matter halo, has stellar mass, $M_{*}=2.5\times 10^{10}M_{\odot}$, a star formation rate of $14M_{\odot}$ yr*-1*, and an intrisic Ly$\alpha$ luminosity of $3.024\times 10^{44}$ erg s*-1*. The intrinsic Ly$\alpha$ luminosity was calculated using starburst99 (Leitherer et al. 1999), summing the contributions from stellar, nebular, and ISM recombination sources and assuming the following ionising and Ly$\alpha$ radiation escape fractions: $f_{esc}^{ion}=0.02$; $f_{esc}^{Ly\alpha}=1$. The star formation rate is comparable to the median star formation rate of the ‘Ly$\alpha$ Em’ sub-sample from S2011, which is $18.6M_{\odot}$ yr*-1*.

The simulations provide, as part of their output, the sites where star formation occurs during the simulation. This is the source of the stellar feedback which drives the galactic wind. In order to calculate the photon budget for the radiative transfer post-processing we model the star forming regions with the starburst99 population synthesis code. Each star-forming particle is treated as a simple stellar population (SSP), where we assume a Chabrier (2003) Initial Mass Function (IMF), to be consistent with the assumptions made for the hydrodynamical simulations. We use the instantaneous star formation mode of starburst99, and integrate in time. Thus for the duration of the burst we ascertain the mean ionising photon count, Ly$\alpha$ luminosity, and the time-averaged ionising photon spectrum of the SSP.

Finally, for comparison to our fiducial model, we also run hydrodynamical simulations with the same initial conditions but with less efficient feedback. For these runs we use the Springel & Hernquist (2002) feedback prescription which does not give an explicit kick to wind particles nor decouple them from hydrodynamics. Throughout this paper we refer to this model as the inefficient feedback model.

2.2 Ionising Photons Radiative Transfer

A smaller box of side $\sim 600$ comoving kpc (corresponding to $\sim$ 328 physical kpc at $z$=2.65 and $\sim$ 187 physical kpc at $z$=2.2) was cut out of the SPH simulation output and gridded onto a discrete, uniform grid with dimension $N_{c}=256$. This ensures a high-enough resolution to give converged results in the radiative transfer calculations. The output provides the temperature and density fields of the gas, as well as the location and luminosity of the stellar populations used as ionising radiation sources. Once augmented with information on the ionisation state of the gas, it can be used as initial conditions for performing the ionising radiation radiative transfer of the galaxy’s local sources.

To obtain the initial conditions for the ionisation state we assume the presence of a Haardt & Madau (2012) uniform UV background (UVB) and photoionisation equilibrium between the gas and the UVB. With the initial conditions defined by the above procedure the effect of the local sources given by the hydro simulations is added with crash (Ciardi et al., 2001; Maselli et al., 2003, 2009; Pierleoni et al., 2009; Graziani et al., 2013), a ray-tracing Monte Carlo 3D radiative transfer code which follows the propagation of the ionising continuum and its effect on the gas it crosses.

The output of crash includes the temperature and ionisation state of each cell in the simulation volume. We use $5\times 10^{4}$ ionising photon packets per source and have checked that the results are converged such that the ionisation state of the simulation box is not significantly affected by using more photon packets. We refer the interested reader to the original papers for more details on the code crash.

The recombination timescale of the ionised gas is long, but nevertheless some recombination should occur. Because the UVB is not explicity included in the radiative transfer calculation, the gas at large galactocentric radius which has been highly ionised by the initial UVB, could artificially recombine during the crash run. To cope with this, we estimate which cells have a UV flux dominated by the UVB and which by local sources. This is done by summing, for each cell, the ionising flux from all sources assuming a $r^{-2}$ falloff and comparing this total to the UVB ionising flux. The results presented in the main body of this paper use a 1:1 domination criterion. That is, a cell is considered to be UVB dominated if the contribution from the UVB to the the cell’s ionising flux is greater than that from local sources (see § A.3 for a discussion of this criteria). In all subsequent post-processing steps we use either the initial ionisation state calculated under UVB photoionisation equilibrium for the UVB dominated cells, or the crash ionisation state for cells dominated by local sources.

2.3 Ly$\alpha$ Photons Radiative Transfer

As we are neither able to resolve nor handle computationally the interstellar medium we simply remove it from the galaxy and allow the Ly$\alpha$ photons to free-stream through the removed cells. The effect of the ISM is then approximated by parameterizing the Ly$\alpha$ and ionising continuum escape fractions ($f_{esc}^{Ly\alpha}$ and $f_{esc}^{ion}$ respectively). This choice is also motivated by the main goal of the paper, i.e. an investigation of the impact of outflowing/inflowing material and not of the ISM. We remove the ISM based on two criteria: a density threshold and a radial distance threshold. That is, grid cells are tagged as belonging to the ISM and subsequently removed if their density is above a density threshold $n_{th}=0.5$ cm*-3* and they also lie within $\sim$10 kpc of the centre of mass of the galaxy. This galactocentric radius threshold serves to ensure that only gas which is part of the galaxy itself is identified as ISM, and avoids removing high-density clumps in the CGM. In practice our results are insensitive to the exact radius threshold used. Likewise we tested density thresholds of $n_{th}=\{0.1,1.0\}$ cm*-3* and found very little variation in our results. Here we assume a value of $f_{esc}^{ion}=0.02$ as in Gnedin et al. (2008). To perform the Ly$\alpha$ radiative transfer we use a total of $\sim 10^{5}$ photon packets and assume $f_{esc}^{Ly\alpha}=1.0$, but as discussed later we renormalise the results to assume a different Ly$\alpha$ escape fraction.

The Ly$\alpha$ radiative transfer is performed with Super Lyman Alpha Fighter (slaf), a new code which we developed during the course of this work. slaf is a Monte-Carlo Ly$\alpha$ radiative transfer code in the vein of many previous works (e.g. Zheng & Miralda-Escudé, 2002; Cantalupo et al., 2005; Dijkstra et al., 2006; Tasitsiomi, 2006; Verhamme et al., 2006; Laursen et al., 2009; Yajima et al., 2012; Gronke & Dijkstra, 2014), and can be applied to arbitrary 3D gas distributions and velocity fields. All Ly$\alpha$ photons are injected into the CGM at line-centre. We have briefly investigated the effect of a Gaussian injection line profile and find that for a Gaussian with a standard deviation of 150 km s*-1* the shape of the surface brightness profile is not altered considerably. However, we leave the discussion of a detailed treatment of the injection profiles for a future paper. We also developed a code, Lyman Alpha Fighter (laf), to calculate absorption line profiles for sightlines through the simulation volume, using the ionisation state and temperature output of the crash radiative transfer. With these two codes we can generate emission maps (using slaf), and absorption equivalent-width maps (using laf). Technical details of slaf and laf are discussed in Appendix A.

We include a very simple dust model in our Ly$\alpha$ radiative transfer calculations. Dust is treated as a grey absorber with a characteristic grain radius and dust-to-hydrogen ratio tuned to satisfy the redshift-dependent dust optical depth, $\tau_{dust}$, relation from Garel et al. (2012). This gives the dust optical depth as a function of HI column density, $\tau_{dust}(\lambda)=\left(A_{\lambda}/A_{V}\right)_{Z_{\odot}}\left(Z/Z_{\odot}\right)^{1.35}\left(N_{HI}/2.1\times 10^{21}cm^{-2}\right)\left(1+z\right)^{-\frac{1}{2}}$. As in Garel et al. (2012), $(A_{\lambda}/A_{V})_{Z_{\odot}}$ is the solar metallicity extinction curve from Mathis et al. (1983), Z is the gas metallicity, $Z_{\odot}$ is the solar metallicity, and $N_{HI}$ is the HI column density. This prescription is applied to each cell in the simulation volume with the assumption of solar metallicity everywhere. Although this assumption is of approximate nature, our testing has shown that the results are insensitive to it.

3 Results

In this section we present the results we obtain by post-processing the galaxy in the manner described above. We break down the results by the two quantities which can be directly compared to observations. For reference we also present the results from the inefficient feedback model, though we leave detailed comparative discussion to § 4.

3.1 Absorption

The upper-left panel of Figure 1 shows a synthetic equivalent width map derived from the absorption line profiles output by laf. The first step to create the map is to choose an orientation from which to view the radiative transfer simulation volume. Then, from each cell on the far face of the box a ray is cast perpendicular to the face, towards the observer, for a total of $N_{c}^{2}$ rays. Each ray starts with a perfectly grey spectrum, which for simplicity we set to 1.0 such that $\{\forall\lambda:I_{0}(\lambda)=1.0\}$. As it traverses the simulation volume the absorption line profile of the intervening HI is imprinted onto the spectrum. This is accomplished by Doppler shifting the existing spectrum into the gas frame of each cell in turn. For the $n^{th}$ cell the optical depth, $\tau_{n}(\lambda)$, is calculated for each sampled wavelength of the spectrum, and the intensity exiting the cell derived as $I_{n+1}(\lambda)=I_{n}(\lambda)e^{-\tau_{n}(\lambda)}$ , where $I_{n}(\lambda)$ is the intensity entering the cell. Thus, when each ray reaches the observer, its spectrum is known. From this, the equivalent width of the total imprinted absorption line for each ray can be calculated as $EW=\int_{line}[1-\frac{I(\lambda)}{I_{0}(\lambda)}]d\lambda$, where $I(\lambda)$ is the final intensity reaching the observer and as noted earlier $I_{0}(\lambda)=1.0$. This results in an $N_{c}\times N_{c}$ equivalent-width map.

Our radiative transfer simulation box is 187 kpc (physical) on a side at redshift 2.2. At this scale the IGM intervening between the edge of the box and the observer could have an impact on the calculation of $EW$, and should in principle be taken into account. As per Laursen et al. (2011), this would primarily affect the spectrum blueward of line-centre. However, since the IGM at z=2.2 is mostly ionised, we do not expect significant interaction with the Ly$\alpha$ line. To explicitly test the effect of the size of our radiative transfer simulation box we simulated a 2x larger box, centred on the same point, and found good agreement with the results presented here.

The equivalent width maps cannot be directly compared to observational data due to the fact that in reality not every line of sight has a sufficiently bright background galaxy or quasar along it. In other words, the relative rarity of galaxy-galaxy pairs and galaxy-quasar pairs prohibits the creation of such a map observationally. Nevertheless it is instructive to see what could be revealed in absorption given large numbers of background sources, especially in light of the other panels of Figure 1. The upper-right panel of Figure 1 shows the HI density averaged along the line of sight of each pixel in the map, the lower-left panel shows the projection of the velocity field of the central slice, and the lower-right panel shows a composite of all the other panels. We can see that generally speaking the morphology of the absorption equivalent width map follows that of the underlying HI distribution. It is clear that there is a correspondence of structures in the absorption map to the structures in the density map, such as the heavy central concentration of HI and the dense filaments. Examination of the velocity field (lower panels) shows material inflowing along the filaments, and being blown out in a biconical outflow. Although the absorption map and velocity fields are somewhat messy, the lower-right panel shows how the velocity field (and hence outflow) affects the absorption map. Unsurprisingly, since the absorption depends on both the density and velocity field, larger velocities tend to coincide with larger absorption equivalent widths. Figure 1, in particular the lower-right panel, demonstrates the complex interplay between these three quantities.

Figure 2 is a plot of absorption equivalent width vs impact parameter, $b$, i.e. an absorption equivalent width radial profile generated from the absorption map in Figure 1 (upper-left). It was created by sampling points in the map radially from the centre of the galaxy (denoted by a ‘+’ in Figure 1). Each pixel from the absorption map in Figure 1 is plotted in Figure 2 as a black point. The solid red line shows the equivalent width profile obtained by azimuthal averaging of the map. The blue points and error bars are the data points from the galaxy-galaxy pairs in S2010. Figure 2 shows that for this viewing angle our simulation provides good agreement with the S2010 observations from $\sim 40$ kpc outwards, mildly underpredicts absorption compared to the data point at $b=31$ kpc, and increasingly underpredicts the absorption towards the centre of the galaxy.

Figure 3 shows the absorption equivalent width profile as viewed from all three axis-aligned orientations as green, red, and blue solid lines. The red line in Figure 3 corresponds to the data in Figure 2. The mean of the three orientations is shown as a black solid line. We note that the mean equivalent width profile from our simulations is in quite good agreement with the observations. Similarly to the single orientation presented in Figure 2, the simulation is within the observational error bars for the outer 2 observational points, $\sim 12\%$ below at $b\sim 31$ kpc and $\sim 43\%$ below at $b\sim 0$. The lower absorption equivalent width predicted by our simulations in the inner $\sim 40$ kpc of the galaxy is indicative of the simulated HI gas density being too low, the simulated velocity dispersion being too low, or a combination of both. We return to discuss this topic later in § 4, taking into account the emission results presented in § 3.2.

For comparison, the inefficient feedback model is also shown in Figure 3 with dashed lines. Since this model lacks strong stellar feedback the galaxy has formed more stars and thus has a stellar ionising photon luminosity around 4 times that of the galaxy formed with the fiducial model. In order to make a meaningful comparison then, we rescale the luminosity to match that of the galaxy formed under the strong feedback model.

Figure 3 shows that the inefficient feedback model also does a good job of reproducing the observations in the outer CGM. However it shows markedly worse agreement with observations in the inner $\sim 40$ kpc, consistently underpredicting the absorption equivalent width.

3.2 Emission

Summing the luminosity of the photon packets exiting the simulation box, with the assumptions that the Ly$\alpha$ escape fraction is 100% and that all Ly$\alpha$ radiation is emitted by the stars in the simulation, gives an observed Ly$\alpha$ luminosity of the simulated galaxy of $3.023\times 10^{44}$ erg s*-1*. This is extremely close to the intrinsic luminosity (see § 2). In the simulations of the Ly$\alpha$ emission of our galaxy, we ignore the component of the signal arising from recombination radiation in the CGM. We justify this by noting that a calculation of the recombination rate, and the resulting Ly$\alpha$ emission, shows that the contribution from recombination within the simulated volume is less than 10% of the total Ly$\alpha$ emission333see § A.4 for further discussion.

Figure 4 shows the Ly$\alpha$ surface brightness profile of the galaxy as viewed from a single side. The data in red assumes $f_{esc}^{Ly\alpha}=0.144$, i.e. 14% of the Ly$\alpha$ photons generated in the galaxy diffuse through the ISM and escape into the CGM without attenuation. This is motivated by S2011, which quotes the Ly$\alpha$ escape fraction as 14.4% for the galaxy sub-sample we compare to. The data in blue show the upper limit of $f_{esc}^{Ly\alpha}=1.0$. Changing $f_{esc}^{Ly\alpha}$ simply modulates the total energy of the Ly$\alpha$ photons injected into the CGM. Since we use enough photon packets to sample the radiative transfer, we similarly modulate the energy assigned to each photon packet to vary $f_{esc}^{Ly\alpha}$, which has the result of shifting the normalisation of the points in Figure 4.

Here, solid lines represent the mean radial profile (which is the observational quantity plotted in S2011) of the simulated surface brightness images, dashed lines represent median radial profiles, and the shaded regions show the lower and upper quartiles. The simulated surface brightness images are degraded to 1” FWHM resolution prior to creating these profiles, in order to match the resolution of the S2011 data, represented by the green dashed line.

Interestingly, at times the mean profile rises above the median profile, and even above the upper quartile of the distribution at some impact parameters. This is indicative of the fact that at these impact parameters the distribution is not Gaussian but skewed, dragging the mean up. This is caused by a small amount of substructure in the surface brightness images - small, bright star-forming clumps. The impact parameter at which these clumps reside depends on projection effects, and thus on the observation angle. In this paper we only simulate one galaxy but given the random projection of substructure we expect that a fairer comparison, stacking many simulated galaxies, would have a similar mean to the median. This is expected because the effect of outlying substructure at a given impact parameter would be diluted when averaged over more galaxies, most of which would not exhibit substructure at the same impact parameter.

We argue then that in the case of Figure 4 where there is little substructure it is perhaps better to compare the median profile to the results from S2011, which are themselves a stack of 52 galaxies. Considering the median profile in Figure 4 it is clear that our simulation is a good fit to observations from $\sim 40$ kpc out to where the observations end at 80 kpc. Below $\sim 40$ kpc the surface brightness profile given by our simulation starts to rise above the observed profile and becomes too peaked. The model proposed in Dijkstra & Kramer (2012) also exhibits a similar behaviour in the inner region.

Figure 5 shows the average surface brightness profile obtained by stacking the six surface brightness maps corresponding to viewing the simulated galaxy from the six sides of the simulation volume. When stacking several viewing angles the shape of the profile remains very similar as there is suprisingly little variation in the profile when the galaxy is viewed from different orientations. The main difference between the profiles seen from different orientations is in the location and size of the bumps and peaks. The large bump at 80 kpc comes from a bright peak visible along a single axis.

Figure 5 also shows the inefficient feedback model in red. As was done for the absorption results, the stellar luminosity of the inefficient feedback model has been scaled down to match that of the fiducial model. The medians of the two models are almost indistinguishable, thus the median surface brightness profile from the inefficient feedback model is also compatible with the data from $\sim 35$ kpc to 80 kpc. However, whereas the fiducial model has good agreement between the mean and median for the vast majority of the profile, the mean of the inefficient feedback model is significantly above the median from $\sim 35$ kpc outwards. This discrepancy arises because in the inefficient feedback model there is much more substructure over a large range of impact parameters compared to the fiducial model, in which there are only a few, isolated clumps. From this single realisation of a galaxy it is not possible to tell if the inefficient feedback model generically predicts an increased presence of substructure, but if this were indeed the case, then we argue that here the median is not a good indicator of what can be expected from mean-stacking many simulated galaxies.

Figure 6 shows the spectrum of all Ly$\alpha$ photons exiting the radiative transfer simulation box. There is a clear distinction between the two models. The fiducial model has a more pronounced redward peak, and a reduced blueward peak relative to the inefficient feedback model. This can be understood as a less extreme example of the effect shown in Figure 8, whereby an expanding medium boosts the redward peak and suppresses the blueward peak. Since the fiducial model has stronger feedback, and thus a stronger outflow, this is exactly the behaviour we expect to see in the spectrum. The larger velocity shift of the red peak has the same origin. As we discuss in §  4 the velocity offset of the peaks in both models is too low compared to observations, but aside from that, the line shape of the fiducial model is generally in better agreement with observations (Steidel et al. 2010).

3.3 Combined Absorption and Emission

Comparing our fiducial model with the inefficient feedback model we can state that judged by absorption alone the former better fits observations, while if we consider emission only the message is less clear. Assuming that for both models the median surface brightness profile is a good proxy for what would result from mean-stacking many more simulated galaxies, then neither model is a better fit to observations than the other. However, this assumption may not be true, and instead the mean profile of the inefficient feedback model could be a better proxy for the mean-stacked profile. If this is the case then the fiducial model again provides a better fit to observations.

Taking both absorption and emission into account it is clear that the fiducial model is favoured, with the caveat that this study is limited to a single simulated galaxy. At this stage, this preference comes primarily from the absorption result. As explained above, whether the emission actively supports this preference or merely does not contradict it hinges on the results of simulating more galaxies.

4 Discussion

We have demonstrated that our fiducial wind model can reproduce some of the observables detailed in S2010 and S2011, although it struggles in the inner regions. At first glance, the partial success of the outflow model is surprising since it is by necessity extremely simplified. In particular, it has the non-physical characteristic of temporary decoupling of the kicked wind particles, which is necessary to produce efficient feedback. In fact, Dalla Vecchia & Schaye (2008) suggest that the density threshold at which hydrodynamical coupling of the kicked particles to their neighbours is reinstated may be the dominant factor in determining the morphology of this class of outflow models.

In more detail, in § 3.1 we noted that the simulated absorption equivalent widths within the inner $\sim 40$ kpc of the galaxy are lower than observations, suggesting that in this region either the simulated HI gas density is too low, the simulated velocity dispersion is too low, or a combination of both. In § 3.2 we saw that the simulated emission profile is too peaked compared to observations, i.e. the Ly$\alpha$ photons are arriving at the observer too directly (to see this intuitively consider the case of a point source in a vacuum - with no scattering the observer would see a 1-pixel wide step function emission profile). The emission surface brightness profile results confirm one of the options from the absorption result: a higher gas density would lead to more Ly$\alpha$ photon scatterings, flattening the emission profile as required to match observations; a higher velocity dispersion would Doppler shift the HI away from Ly$\alpha$ resonance, allowing the Ly$\alpha$ photons a more direct path to the observer and further increasing the peakiness of the emission profile. Thus the emission surface brightness profile results suggest that a more successful wind model should have a higher central HI gas density. However, this alone cannot be the whole explanation. In both cases the emission spectra (Figure 6) show a red peak at a velocity offset less than that observed ($\sim$400 km/s, see Steidel et al. 2010). This is indicative of the outflow velocity not being high enough. It appears then that a combination of both increased HI density in the inner regions and increased outflow velocities is needed to bring simulations into line with all of the available observations.

While we need further simulations to confirm this, we speculate that since both of the tested models show similar absorption and emission profiles at $b>50$ kpc, the CGM signal in Ly$\alpha$ absorption and emission may be largely unaffected by outflows at these impact parameters. Indeed, Shen et al. (2013) show that outflows do not disrupt cold inflows. Thus, if the signal at $b>50$ kpc is dominated by cold accretion flows it is reasonable to expect that our strong feedback model would not significantly affect the absorption and emission profiles in this region. However, we caution against interpretating the fact that the $b>50$ kpc signal does not appear to be coming from outflows to mean that it must necessarily be caused by inflows.

This limited radial influence of the outflow is perhaps a manifestation of the fact that a wind has some finite sphere of influence. Visually, we can see this by referring to the lower-left panel of Figure 1, where we can see that beyond $\sim 50$ kpc the amplitude of the velocity field drops off. The extent of this sphere will surely depend on the velocity of the stochastic kick given to wind particles, but may also be affected by the aforementioned density threshold at which wind particles become re-coupled to hydrodynamics.

So far we only present the results of one simulated galaxy, and it would be premature to judge the fiducial outflow model on this one result. Since S2010 and S2011 both deal with averaged/stacked data the correct comparison to make is to similarly stacked simulated galaxies. Genel et al. (2012) simulate a suite of galaxies and we suggest a stacked analysis of this data would be a significant improvement on the work presented in this paper.

In fact, the need for more simulations is deeper than may be immediately obvious. The two different feedback models we have presented produce galaxies with different properties from the same halo initial conditions. This, of course, is exactly why we are interested in various feedback models, but as we saw in §  3, it makes comparison of the resulting galaxies difficult. Recall that the ‘re-simulated’ halo was chosen such that the galaxy which forms under the fiducial feedback model was similar to the mean of the S2010/S2011 samples. Using the inefficient feedback model the galaxy which formed was brighter and so we needed to rescale the luminosity of the galaxy to compare to the fiducial model and observations. This introduces an inconsistency: the stars which ionise/illuminate the CGM in Ly$\alpha$ do not have the same properties as those which generated the outflow velocity field and CGM properties. If we instead chose to compare two galaxies with the two models, based on the final formed galaxy properties (that is, choose galaxies from both models to match observations) we are faced with the problem that we would be comparing galaxies with different initial conditions. Since what we are really interested in is the impact of the different models on the evolution of the CGM we prefer the approach we have taken where we keep the initial conditions the same across the two models. Nevertheless, we acknowledge that this is not ideal and suggest that the real solution is to simulate a large statistical sample of galaxies for both models, from which a sample with properties matched to observations can be selected and compared.

Other outflow models such as those described in Dalla Vecchia & Schaye (2012) and Dalla Vecchia & Schaye (2008), avoid the unphysical temporary decoupling from hydrodynamics common to Springel & Hernquist (2003)-style winds. The resulting outflows are qualitatively different in spatial and velocity morphology from those produced by Springel & Hernquist (2003) and would thus make for an interesting comparison. Additionally, there are free parameters in the outflow model used by Genel et al. (2012) which presents a parameter space which it may prove fruitful to explore. A comparative test of outflow models can be found in Chung (2016)

Our simulations do not include spatially extended emission from cold accretion streams, either powered by recombination or collisional excitation. Previous works have shown that this may provide a large Ly$\alpha$ luminosity (e.g. Dijkstra & Loeb, 2009).

Accounting for this emission may boost the amount of Ly$\alpha$ emission at large impact parameters, thus flattening the predicted surface brightness profile.

Our simulations show that a stronger outflow appears to diminish the ‘twin peaks’ profile of the Ly$\alpha$ line, emphasising the red peak and diminishing the blue peak. Recent observations by Leclercq et al. (2017) with MUSE444http://www.eso.org/sci/facilities/develop/instruments/muse.html include spectra with similar asymmetric peaks, which may be attributable to the processes discussed in this paper.

This work has possible relevance to Enormous Ly$\alpha$ Nebualae (Cai et al., 2017) and the enigmatic Ly$\alpha$ blobs reported by Steidel et al. (2000) and Matsuda et al. (2004). It is presently unclear what the physical mechanism powering these objects is. One possibility is that Ly$\alpha$ blobs are powered by a central galaxy with the photons scattered off the surrounding medium, in a similar fashion to that considered in this paper. Hayes et al. (2011) presents observations supporting this model of Ly$\alpha$ blobs, but see also Haiman et al. (2000)

4.1 Comparison with Previous Work

As we have already pointed out, the absorption line data is better fit by our fiducial model, but the inefficient feedback model also does a good job of reproducing the observations at higher impact parameters. The similarity of the profiles in Figure 3 goes some way towards explaining why previous analyses of Goerdt et al. (2012) and Shen et al. (2013) can reproduce the absorption line data without strong outflows, instead attributing most of the Ly$\alpha$ absorption to inflowing cold streams.

Goerdt et al. (2012) reproduces the Ly$\alpha$ absorption line data extremely well, even towards the centre of the galaxy. This is in contrast to our inefficient feedback model, which underpredicts the amount of Ly$\alpha$ absorption in the inner $\sim 40$ kpc. One reason for this may be the fact that Goerdt et al. (2012) do not include the effect of local sources, and instead use a simple self-shielding criteria to calculate the ionisation state of the gas assuming a UV background. It is therefore reasonable that our simulations which contain local sources, mostly concentrated in the galaxy, and explicitly compute the ionisation state, should have a higher ionisation fraction in the central region close to the galaxy. This would naturally lead to a lower Ly$\alpha$ absorption equivalent width. In the outer regions, where local sources are less dominant, we expect better agreement with the Goerdt et al. (2012) result and indeed this is exactly what we see.

An alternate explanation for the discrepancy could be the way in which we calculate the initial conditions for our crash runs. Whereas we assume all gas is in photoionisation equilibrium with the UVB (with no shielding), Goerdt et al. (2012) use a self-shielding criteria. In the dense regions towards the centre of the galaxy this may cause us to overestimate the ionisation fraction of the gas. If this is happening, then it is also happening in our fiducial model, which also lies under the observations at very small impact parameter, and a better treatment of the UVB in our calculations may yield better agreement with observations for our strong feedback model.

Shen et al. (2013) does include local sources, albeit in a simplified fashion, placing all sources at the centre of the galaxy. They have problems similar to ours in reproducing the central two absorption data points from S2010.

With respect to emission, S2011 presents a simple analytic model to explain their observations. They consider a spherically symmetric outflowing HI CGM, modulated by a covering fraction which is a power law function of galactocentric radius. Radiative transfer is treated with an extremely basic prescription. Nevertheless this model provides a good fit to their observations. The success of this model provided motivation for us to test whether the underlying assumption of a central source, emitting radiation which scatters in an outflowing CGM, could stand up to the scrutiny applied by our hydrodynamic simulations and a full treatment of radiative transfer.

Laursen et al. (2009) perform a similar radiative transfer treatment to this work. However, their work differs in some important ways. The galaxies they simulate are not targeted to match the sample of galaxies from S2011, as they are at a different redshift ($z$=3.6) and with no strong feedback. Nevertheless, their results are comparable with a very similar surface brightness profile shape.

Finally, on the observation side, it would be disingenuous to omit mention of recent work by Hennawi et al. (2015) and Cantalupo et al. (2014) which cast doubt on the fidelity of current numerical cosmological simulations of the type used in this work. Hennawi et al. (2015) use background quasar absorption lines to study the properties of foreground quasar hosts. They observe much more cool gas than simulations predict, suggesting that ‘essential aspects’ of massive halo hydrodynamics at $z\sim 2$ are not being captured by current cosmological simulations. Cantalupo et al. (2014) use fluorescent Ly$\alpha$ emission to study the cold gas mass of the nebula surrounding quasar UM287, and again find an excess of cold gas relative to simulations. While the case against cosmological simulations brought by these works gives cause for concern, it is consistent with our findings as discussed above - namely that additional HI is required to bring our simulations into agreement with the observations.

5 Summary

We have introduced a new test for galactic outflow models, which combines hydrodynamical simulations and Ly$\alpha$ radiative transfer in a self-consistent way where the stars driving the outflows are also responsible for the ionising and Ly$\alpha$ radiation used in the radiative transfer. Crucially, we use constraints from both Ly$\alpha$ absorption and emission to test our models. The fiducial outflow model which we have presented in this paper can reproduce features of both absorption and emission observations, although the inner region and the Ly$\alpha$ emission spectrum remain problematic.

We also showed that there are differences in these two diagnostics when a different feedback model is used. Furthermore, comparison of our results for the two feedback models hints that galactic outflows may predominantly affect the inner $\sim$50 kpc of the CGM. This suggests that future Ly$\alpha$ observations of the inner CGM may be key to gaining a better understanding of the galactic outflows which appear to be important to galaxy formation and evolution.

Acknowledgements

We thank Luca Graziani for countless fruitful discussions during this work. Many thanks also to Akila Jeeson-Daniel for providing the crash-$\alpha$ pipeline, and all the helpful support along the way. We also thank Kathryn Kreckel, Tyrone Woods, Andressa Jendreieck for constructive comments. Finally, we thank the reviewer Joakim Rosdahl for insightful comments which, when addressed, surely improved the quality and robustness of this paper.

Appendix A code details

In order to perform this work it was necessary to develop some new codes. Here we briefly describe the nature of the two codes, and show the results of some of the verification tests we have applied to them.

A.1 Radiative Transfer (SLAF)

slaf (Super Lyman Alpha Fighter) is a new 3D monte-carlo Ly$\alpha$ radiative transfer code similar to Zheng & Miralda-Escudé (2002); Cantalupo et al. (2005); Dijkstra et al. (2006); Tasitsiomi (2006); Verhamme et al. (2006); Laursen et al. (2009); Yajima et al. (2012). In order to strike a balance between efficiency, extensibility, and maintainability, C++ was chosen as the implementation language, and an object-oriented design paradigm was used.

The papers cited above describe in detail how codes of this nature work. However, here we give some information specific to our code and briefly outline the critical steps involved in Ly$\alpha$ RT. slaf was designed to fit into the existing crash pipeline. Therefore, the input files share the same format as crash.

A number of photon packets are allocated to each source, based on the source luminosity. Each photon packet is then tracked from the source until it exits the simulation box. The first step is to choose a random direction in which to emit the packet from the source. To propagate the photon packet from scattering to scattering we need to randomly choose an optical depth, $\tau$, that the photon will freely stream through from the probability distribution $e^{-\tau}$. We then propagate the photon packet from cell to cell, calculating the optical depth ‘used up’ in traversing each cell, based on the frequency of the photon and the physical properties of the cell. As the photon packet freely streams from cell to cell, a dust attenuation factor is applied, based on the HI density of the cell and a specified dust-to-gas ratio. When the optical depth of the photon packet is exhausted the photon undergoes a scattering event.

For each scattering event we must first generate the thermal motion of the atom which we are scattering off. This is done by picking the thermal velocity from a non-trivial distribution, details of which are given in Zheng & Miralda-Escudé (2002). In the current version of slaf scattering is isotropic. The direction of the scattering together with the thermal motion of the atom determines the gas-frame frequency shift of the scattered photon packet. We are then ready to begin the loop again and propagate to the next scattering event.

To generate surface brightness images we follow the method of Yusef-Zadeh et al. (1984). We set up virtual detectors on each face of the simulation box. Then, for each scattering event we calculate the probability that the photon scatters in the direction of each ‘detector’. In our case this means that we calculate the probability of scattering in the direction of the 6 faces of the simulation box. For each scattering event and for each detector the energy of the photon at the time of the scattering event is weighted by the probability of scattering towards said detector and added to the pixel that the scattering event projects onto at that detector. At the end of the radiative transfer each detector has a relative brightness map. In order to convert these relative brightness maps to absolute values we sum the total energy of all photons which exit the simulation cube and assign it proportionally over the maps.

In order to verify our implementation, we compare the slaf output with the test case in Dijkstra et al. (2006), for which there is an approximate analytic solution (Dijkstra et al., 2006, Equation 9). In brief this test puts a single, monochromatic Ly$\alpha$ source at the centre of a homogeneous static sphere of HI gas at 10 K. In each instance the uniform gas density is chosen so that the line-centre optical depth from the source to the edge of the gas sphere is {$10^{5}$, $10^{6}$, $10^{7}$}.

The outcoming spectrum from slaf is compared to the analytic solution in Figure 7. slaf is shown to be in excellent agreement with the analytic solution with the exception of a systematic amplification of the red peak and suppression of the blue peak compared to the analytic solution. We understand this to be due to the ‘recoil effect’ (Field, 1959; Zheng & Miralda-Escudé, 2002; Tasitsiomi, 2006), the thermalisation of photons which occurs due to the photon scattering in a different direction to the incident photon, which is ignored in the analytic approximation for simplicity. Figure 7 (inset) shows the output of slaf with the ‘recoil effect’ disabled, and demonstrates almost perfect agreement with the analytic solution.

To test that slaf can also correctly handle dynamic scenarios we make a comparison to the results presented in Tasitsiomi (2006) (Figure 3, right panel). Since we want to make a direct comparison, we set up an identical scenario comprising of a single central Ly$\alpha$ source embedded in a spherically symmetric HI cloud.

The HI cloud has a column density of $2\times 10^{18}$ cm*-2* from the centre to the edge of the cloud, and a Hubble-like velocity gradient where the gas velocity scales proportionally to the radius up to a maximum inflow/outflow velocity of 200 km s*-1* at the edge of the cloud.

Figure 8 shows our results overlaid on the results from Tasitsiomi (2006). The results in blue (red) are for the expanding (contracting) cloud case. The solid histograms show our slaf results and the dashed-line histograms show the Tasitsiomi (2006) results. There is excellent agreement between the two codes.

A.2 Absorption along Sightlines (LAF)

laf (Lyman Alpha Fighter) is another new code which calculates Ly$\alpha$ absorption line profiles for sightlines through 3D volumes of arbitrary gas distributions and associated velocity fields. Equivalent widths for sightlines can then be derived from the line profiles.

In order to test laf we construct a simple test case. We create a sphere of HI gas with radius $r=25.6$ kpc, uniform HI number density $n=1\times 10^{-10}$cm*-3*, and a temperature of 5000 K on a 5123 grid.

We then use laf to calculate the absorption line profiles, and hence absorption equivalent width, for lines of sight through the simulation volume. For such a simple geometry the optical depth along lines of sight through the volume can be calculated analytically, and used to derive the equivalent width of the absorption lines along the line of sight.

Figure 9 shows a comparison between the laf output and the analytic solution. The output of laf is in excellent agreement with the analytic solution, with a small deviation arising only from the fact that the sphere has been quantised onto a grid. This can be seen in the fact that the agreement gets worse as the impact parameter approaches the radius of the sphere. This can be understood by realizing that as a sightline moves close to the edge of the sphere the segment intersecting the gas sphere decreases. Thus the relative error due to the discretisation onto a grid grows as the traversed optical depth decreases with higher impact parameter.

A.3 Robustness to the chosen domination criterion

There is scope for concern about the specific criterion chosen to ascertain whether a particular cell is dominated by the ionising radiation from local sources or from the UVB. Our pipeline uses this criterion to decide whether a cell’s ionisation state is calculated via ‘real’ ionising radiative transfer using crash or assumed to be the UVB equilibrium ionisation state. However the issue is not as critical as it may at first appear. This is because the UVB is never fully ignored, even when a cell is determined to be local-source dominated.

To calculate the ionisation state of cells the ionisation equilibrium state with the UVB is first calculated for all grid cells. This is used as the initial conditions for the next step, which is to add in local ionising radiation sources, and perform radiative transfer to calculate the new ionisation state of cells. Because our radiative transfer simulation does not include the UVB, recombinations may occur during the radiative transfer of the local ionising photons which in reality should not (due to the presence of the UVB). This is why we test whether each cell is dominated by the UVB or by local sources. If a cell is UVB dominated, we use the UVB equilibrium state, otherwise we use the result of our local-source radiative transfer. However, this local-source dominated ionisation state already implicitly includes the UVB in the initial conditions.

Robustness to the domination criterion was also tested explicitly. To do this, the criterion used to decide whether a cell is local-source or UVB-dominated was adjusted. In the main paper, a cell is deemed to be dominated by the UVB if the ionising flux from the UVB is greater than the ionising flux from local sources. We modified the criterion so that a cell is considered to be UVB dominated if the ionising flux is 2 or 3 times higher than the local flux. This serves to inform us how sensitive our procedure is to this criterion.

The results (Figure 10) show that altering the domination criterion affects the surface brightness profile by up to $\sim 15\%$ at any given impact parameter (and much less towards the centre of the galaxy). On the logarithmic scale at which the data are presented, and the trends manifest themselves, this $\pm 15\%$ variation is so small as to be irrelevant.

A.4 On ignoring Ly$\alpha$ from recombination

As per § 3.2, Ly$\alpha$ emission resulting from recombinations in the CGM is ignored in our simulations. For brevity this is justified in the main text by noting that recombination in the CGM accounts for less than 10% of the total Ly$\alpha$ emission. However, when deciding whether recombination emission can be ignored, we must also consider the spatial distribution of the emission. If the spatial distributions of the two sources differ significantly, and in such a way that the recombination emission becomes dominant at certain impact parameters, then radiative transfer effects could potentially alter the resultant surface brightness profile.

We directly test for this by estimating the Ly$\alpha$ recombination emission per cell of the simulation. This assumes Case A or Case B recombination, depending on the physical state of each cell. The radial distribution of the recombination emission is then compared to the radial distribution of the stellar emission (Figure 11 and Figure 12).

Figure 11 shows that the distribution of the recombination emission is strongly peaked around the galaxy, which is also the dominant source of Ly$\alpha$ from stellar populations. The peak of the recombination emission is offset from the centre of the galaxy. This is due to the fact that in our simulations the very high density ISM gas is replaced by vacuum and a Ly$\alpha$ escape fraction(see § 2.2). Therefore the highest amount of recombinations occur at the interface of this (artificial) vacuum region. Because of this offset the plot appears (falsely) to show that the recombination emission dominates at a galactocentric radius of $\sim$8 kpc. However, the vacuum which causes the recombination emission to be peaked offset from the centre of the galaxy, also causes the stellar emission to have a similar offset. The Ly$\alpha$ photon packets from the stellar sources in the vacuum region free-stream until they hit this interface, and so the effective position of the sources (and all of the emission in the central $\sim$8kpc in Figure 11 and Figure 12) is also at this interface. This is not manifested in Figure 11 and Figure 12 because they show the spatial distribution of the sources, ignoring radiative transfer effects.

This quirk of our simulations aside, it is clear from Figure 12 that the stellar sources dominate at all galactocentric radii. Due to the vagaries of Ly$\alpha$ radiative transfer, we cannot guarantee that the CGM recombination emission is negligible without explicitly including it in the radiative transfer. Nevertheless, we hope that the above discussion gives some insight into why, on balance, we deem it highly probable that Ly$\alpha$ from CGM recombination can safely be ignored with respect to the main results of this paper.