text:
\section{Comparing to stacked \mbox{Ly$\alpha$}\xspace surface brightness profiles from MUSE observations}
\label{ap_11_comp}

\begin{figure}
\includegraphics[width=20pc]{Figures/Lya_projected_SB_stacked_mocks.pdf}
\caption{A comparison of the predicted \mbox{Ly$\alpha$}\xspace surface brightness profile from our simulation
(for $3<z<4$) compared to the stack of LAEs produced using deep MUSE data by \protect \cite{Wisotzki18}.
Solid (dashed) lines show profiles without (with) convolution with the appropriate MUSE point spread function (PSF).
In this case we more closely reproduce the observational selection and methodology by median stacking (after azimuthal averaging),
by using mock data cubes orientated along random lines of sight (rather than angle averaging)
and by applying a flux cut to mimic the sensitivity of MUSE in detecting LAEs.
We choose a flux cut of $\log_{10}(F_{\mathrm{\mbox{Ly$\alpha$}\xspace}} \, / \mathrm{erg s^{-1} cm^{-2}}) > -17.6$
in order to match the peak surface brightness of the observed stack, which is comparable
to the approximate flux limit for which the MUSE observations are complete \protect \cite[see, e.g., figure 1 in][]{Leclercq17}.
With this flux cut, we reproduce the observed profile to within at least $0.2 \, \mathrm{dex}$
over the plotted range.
}
\label{lya_sb_mock}
\end{figure}

In Section~\ref{MUSE_comp_sec}, we compared a mean-stacked and angle-averaged surface brightness
profile of our simulated galaxy with a stack of observed LAEs from \cite{Wisotzki18}.
In Fig.~\ref{lya_sb_mock}, we perform a more realistic comparison with the observations by median stacking
(after first azimuthal averaging), by using mock data cubes orientated along $12$ random lines of
sight for each simulation output (rather than angle averaging), and by applying a flux cut
to mimic the sensitivity of MUSE in detecting LAEs.
We choose a flux cut of $\log_{10}(F_{\mathrm{\mbox{Ly$\alpha$}\xspace}} \, / \mathrm{erg s^{-1} cm^{-2}}) > -17.6$
in order to match the peak surface brightness of the observed stack, which is comparable
to the approximate flux limit for which the MUSE observations are complete \protect \cite[see, e.g., figure 1 in][]{Leclercq17}.
With this flux cut, we reproduce the observed profile to within at least $0.2 \, \mathrm{dex}$
over the plotted range.

The mock data cubes used here make use of the peeling off technique in order to obtain
sufficient signal along specific lines of sight (which in turn is needed to apply an observational
flux cut). By using peeling off, we lose the information about the origin of each individual photon,
which is why we do not use the mock cubes for the results presented in \protect Sec~\ref{flows_sec}.
With the adopted flux cut, the simulated surface brightness profile in any case very closely resembles
the mean, angle-averaged stack shown in \protect Fig.~\ref{SB_profile}.

\section{Introduction}

Spatially extended ``haloes'' of \mbox{Ly$\alpha$}\xspace emission are observed around
galaxies across a range of scales, ranging from
massive radio galaxies and quasars \cite[e.g.][]{Heckman91,Reuland03,Cantalupo14},
to UV bright Lyman Break Galaxies \cite[LBGs][]{Steidel11},
down to comparatively UV-faint \mbox{Ly$\alpha$}\xspace-emitters (LAEs)
\cite[e.g.][]{Rauch08,Matsuda12,Hayes13,Wisotzki16}.
These observations demonstrate that there must be substantial amounts
of dense hydrogen gas in the circum-galactic
medium (CGM) around high-redshift galaxies. Comparisons to the incidence
rates of strong hydrogen absorbers (specifically Lyman limit systems)
seen in quasar spectra support the idea that this extended emission is tracing
similar gas \cite[][]{Wisotzki18}, and direct evidence for connection between faint
high-redshift LAEs and strong hydrogen absorbers is now starting to appear
\cite[][]{Lofthouse19,Mackenzie19,Muzahid19,Diaz20}.

The detection of spatially extended \mbox{Ly$\alpha$}\xspace emission around high-redshift
galaxies poses a number of questions. First, what component(s) of the
CGM is being traced? Cosmological simulations predict for example the presence
of inflowing filamentary streams of dense gas \cite[e.g.][]{Keres05,Dekel06}, and
in addition dense gas could be either entrained or formed in situ within galactic-scale
outflows \cite[e.g.][]{Wang95,Heckman00,Scannapieco17}.
Second, a related but distinct question is
what is the production mechanism for the extended \mbox{Ly$\alpha$}\xspace emission?
An obvious mechanism is that \mbox{Ly$\alpha$}\xspace is emitted within ionized HII
regions inside the interstellar medium (ISM), which then resonantly
scatter from neutral hydrogen in the CGM \cite[e.g.][]{Laursen07,Steidel10,Zheng11}. Alternatively,
the observed photons could be emitted within the HII regions of satellite
galaxies, or in situ within the CGM due to photo-ionization by
escaping ionizing photons (fluorescence) originating from the host galaxy or quasar
\cite[e.g.][]{Haiman01,MasRibas16,Gallego18}, from
satellites \cite[e.g.][]{MasRibas16,MasRibas17}, or from the wider ultraviolet background \cite[UVB, e.g.][]{Furlanetto05}, or
because of collisional excitations and recombinations that dissipate
the energy gained by compressive heating, which can be triggered
by cosmological gas infall \cite[dubbed ``cooling radiation'', e.g.][]{Fardal01,Dijkstra06,FaucherGiguere10,Rosdahl12},
or by supernova shock waves that break out of the ISM \cite[][]{Mori04}.

\cite{Steidel10} propose that spatially extended cool gas within galactic
outflows can explain simultaneously the kinematics of metal absorption
lines detected ``down the barrel'' in the spectra of high-redshift LBGs, the \mbox{Ly$\alpha$}\xspace
spectra of LBGs (which are red-shifted from systemic), the equivalent width
of absorbing cool gas as a function of transverse distance (detected via
background sources), and \cite{Steidel11} argue that this picture can
also explain the spatially extended \mbox{Ly$\alpha$}\xspace emission from the CGM
\cite[see also][]{Dijkstra12}.
The question of extended \mbox{Ly$\alpha$}\xspace emission aside, simple outflowing neutral gas
shell models are generally used to explain the spectral morphology of observed
\mbox{Ly$\alpha$}\xspace spectra \cite[e.g.][]{Ahn04,Verhamme06}. In this scenario, \mbox{Ly$\alpha$}\xspace photons escaping a central source
backscatter from side to side of the inner surface of a circum-galactic outflowing shell, until
they shift in frequency sufficiently to be out of resonance with
the blue-shifted side of the shell moving towards the observer.
Outflowing shell models can explain the asymmetry and often
complex spectral morphology of observed \mbox{Ly$\alpha$}\xspace lines, the commonly
observed single peaked \mbox{Ly$\alpha$}\xspace line that is red-shifted from the systemic
velocity, and why there is often less (or no) emission observed at (and blueward of)
the systemic velocity \cite[e.g.][]{Verhamme08,Wofford13,Hashimoto15,Gronke17b}.
Similar models with a more continuous expanding medium can also explain
observed surface brightness profiles \cite[e.g.][]{Song20}.

These observations do not rule out the other scenarios however.
The spectral morphology of observed \mbox{Ly$\alpha$}\xspace lines could be imprinted
on small scales within the ISM, rather than on CGM scales.
Similarly, the spatial location of strongly outflowing gas seen down-the-barrel
in galaxy spectra is highly uncertain, and could be tracing be sub-ISM
scales of order $\sim 100 \, \mathrm{pc}$ \cite[][]{Chisholm16},
and so not connected to the absorbing gas seen in transverse sight-lines
(where the kinematics relative to the galaxy are poorly constrained).
As such, the door remains open for other CGM components to provide
the main explanation for the observed spatial distribution of \mbox{Ly$\alpha$}\xspace haloes,
without the CGM necessarily being responsible for the observed \mbox{Ly$\alpha$}\xspace
spectral morphology.

A natural way to explore these questions is with cosmological simulations,
which in principle can simultaneously capture the complex mix of
circum-galactic gas flows that surround galaxies, sourced by the
filamentary infall of matter from larger scales. It must
be acknowledged that the finite resolution of simulations
\cite[particularly in the CGM, see][]{Hummels19,Peeples19,Suresh19,vandevoort19} limits their ability to accurately represent
multi-phase gaseous media; the resolution required to resolve
the formation of dense clumps in the CGM or in winds has been
estimated to be as small as $0.1 \, \mathrm{pc}$ \cite[e.g.][]{Fujita09,Mccourt18}.
Nonetheless, many studies
of \mbox{Ly$\alpha$}\xspace emission in connection to high-redshift galaxies
have been performed with cosmological simulations
\cite[e.g.][]{Cantalupo05,Tasitsiomi06,Laursen07,Zheng10,Barnes11},
with the simulations progressively increasing in both complexity and resolution.
Simulations have been used to explore the effect of scattering
of centrally emitted \mbox{Ly$\alpha$}\xspace photons from the CGM to create extended emission
\cite[e.g.][]{Barnes11,Zheng11,Lake15}, to study in-situ gravitational ``cooling'' radiation
\cite[e.g.][]{FaucherGiguere10,Goerdt10,Rosdahl12},
as well as fluorescent emission powered by the ultraviolet background (UVB), and by proximate
quasars \cite[e.g.][]{Cantalupo05,Kollmeier10}.

\cite{Verhamme12} use non-cosmological high-resolution simulations of
idealised galaxies and emphasise the importance of allowing a cold ISM
phase to form for studying \mbox{Ly$\alpha$}\xspace escape from the ISM.
Studies of \mbox{Ly$\alpha$}\xspace transfer in cosmological simulations of high-redshift galaxies
with the resolution to start addressing this challenge are now starting to appear.
\cite{Behrens19} study
\mbox{Ly$\alpha$}\xspace escape from a $z \sim 8$ galaxy with maximum spatial resolution
of $25 \, \mathrm{pc}$ and gas mass resolution of $2 \times 10^4 \, \mathrm{M_\odot}$,
employing full \mbox{Ly$\alpha$}\xspace radiative transfer within the ISM and CGM, but using
simplified modelling to account for photo-ionization by local sources.
They find generally low \mbox{Ly$\alpha$}\xspace escape fractions as a result of strong local dust attenuation
around young stellar clusters.
Notably, \cite{Smith19} analyse a galaxy at $z=5-7$ from the FIRE-2 suite of
cosmological zoom-in simulations \cite[][]{Hopkins18}, with
a gas mass resolution of $7 \times 10^3 \, \mathrm{M_\odot}$.
They account for photo-ionization from local
sources by post-processing their simulation with ionizing UV photons,
and also perform full \mbox{Ly$\alpha$}\xspace radiative transfer through both the ISM and CGM.
They study the angular and temporal variations of escaping \mbox{Ly$\alpha$}\xspace emission,
finding that temporal variations are primarily driven simply by the
star formation history of the primary galaxy. They also explore in post-processing
the possible impact of \mbox{Ly$\alpha$}\xspace radiation pressure in and around the simulated
galaxy, finding that it should be dynamically important.

Here, we present results from a full cosmological radiation hydrodynamics (RHD) simulation
of a high-redshift galaxy for $z \geq 3$, taken from the parent sample of
(non-RHD) simulations presented in \cite{Mitchell18}. The simulation self-consistently
includes the full range of expected \mbox{Ly$\alpha$}\xspace powering mechanisms, including
photo-ionization and photo-heating from local sources and the UVB.
The simulation achieves a maximum spatial resolution of $14 \, \mathrm{pc}$
within the ISM at $z=3$,
with a characteristic gas cell (and star particle) mass of $\sim 10^{3} \, \mathrm{M_\odot}$.
This simulation represents a powerful tool to study the radiative
transfer problem for \mbox{Ly$\alpha$}\xspace photons in both the ISM and the CGM simultaneously.
Our simulated galaxy (with a halo mass of $M_{200}= 10^{11.1} \, \mathrm{M_\odot}$ at $z=3$)
is deliberately chosen to be representative of the UV-faint LAEs detected
by the Multi Unit Spectroscopic Explorer (MUSE) in deep MUSE fields ($\geq 10$ hour exposure)
for $3<z<6$ \cite[e.g.][]{Wisotzki16,Hashimoto17,Inami17,Leclercq17}. The substantial
increase in sensitivity afforded by MUSE is significant in this sense,
as this greatly facilitates the computational challenge of producing observable
galaxies in cosmological simulations with high resolution and full radiation
hydrodynamics.

Building on our previous work presented in \cite{Mitchell18}, we set out
in this study to first understand what shapes the dense phases
of hydrogen in the high-redshift CGM. We then use this information
to interpret predictions from our simulation for spatially extended
\mbox{Ly$\alpha$}\xspace emission, and for the spectral morphology of the
\mbox{Ly$\alpha$}\xspace line. We address the questions of both how \mbox{Ly$\alpha$}\xspace traces the CGM,
and of the origin of the extended emission (in terms of in-situ
versus scattered emission that originates from the host galaxy or from satellites),
albeit we do not undertake here the non-trivial task of
separating the \mbox{Ly$\alpha$}\xspace emitted in situ from the CGM between fluorescence
and non-radiative heating sources.
Complementing our study of the CGM, a forthcoming paper (Blaizot et al., in preparation) will
extend the sample of radiation hydrodynamics simulations, and
will focus on the observability of \mbox{Ly$\alpha$}\xspace emission from
the host galaxy, the angular and temporal dependence of the
\mbox{Ly$\alpha$}\xspace escape fraction, equivalent width, and spectral morphology,
and the relative role of collisional excitations and recombinations to
the escaping \mbox{Ly$\alpha$}\xspace signal.

The layout of this article is as follows: we describe our simulation setup,
subgrid models, and post-processing strategy in Section~\ref{methods},
we present our results in Section~\ref{flows_sec}, we discuss various caveats and
implications of these results in Section~\ref{discussion_section}, and summarise
in Section~\ref{summary_sec}.

\section{Methods}
\label{methods}

\subsection{Simulation setup}
\label{simul_description}

The results presented in this paper are taken from a single cosmological zoom-in simulation of a high-redshift
galaxy, simulated down to $z=3$ using the {\sc{ramses-rt}}\xspace code \cite[][]{Rosdahl13,Rosdahl15}, which is a radiation hydrodynamics extension
of the {\sc{ramses}}\xspace code \cite[][]{Teyssier02}. {\sc{ramses}}\xspace is an Eulerian code for hydrodynamics that employs adaptive
mesh refinement. The zoom-in simulation presented uses the same subgrid physics models and simulation parameters
as the fiducial {\sc{sphinx}}\xspace simulations presented in \cite{Rosdahl18}, and effectively extends the project
from $z=6$ (the final redshift currently achieved in {\sc{sphinx}}\xspace) to $z=3$. The initial conditions for the simulation
are taken from the parent sample described in \cite{Mitchell18}, with this simulation representing a resimulation
of the most massive halo from that study, with a halo mass of $M_{200} = 10^{11.1} \, \mathrm{M_\odot}$ and
virial radius of $R_{\mathrm{vir}} = 37 \, \mathrm{kpc}$ at $z=3$.

\defcitealias{Planck13}{Planck Collaboration, 2014}
We assume an underlying $\Lambda$ Cold Dark Matter cosmological model, with $\Omega_{\mathrm{M}} = 0.3175$, $\Omega_{\mathrm{\Lambda}} = 0.6825$,
$\Omega_{\mathrm{B}} = 0.049$, $H_{0} = 67.11 \, \mathrm{km s^{-1} Mpc^{-1}}$, $n_{\mathrm{s}} = 0.962$, and $\sigma_8 = 0.83$
\citepalias{Planck13}. We use the {\sc{music}}\xspace code to generate initial conditions \cite[][]{Hahn11}. The zoom-in simulation
is performed within a low-resolution box of volume $(20 \, \mathrm{cMpc}/h)^3$, with the high-resolution
region selected as a sphere of $150 \, \mathrm{pkpc}$ radius around the target halo at $z=3$.
We assume hydrogen and helium mass fractions of $X=0.76$ and $Y=0.24$ respectively, and we initialise cells
with a starting metallicity of $3.2 \times 10^{-4} \, Z_{\odot}$, assuming the Solar metal mass fraction $Z_{\odot}=0.02$.
This higher-than-pristine metallicity value is chosen to account for the enrichment from stars forming
in unresolved haloes with masses below the atomic cooling limit.

For the hydrodynamics, we solve the Euler equations using a second order Gudunov scheme, using the Harten-Lax-van Leer-Contact (HLLC)
Riemann solver with a MinMod slope limiter to construct fluid variables at cell interfaces. We employ a Courant factor
of $0.8$ to regulate the simulation timestep, and we assume an adiabatic index of $\gamma=5/3$ to close the
relationship between gas pressure and internal energy, appropriate for an ideal monoatomic gas. Gravitational dynamics are solved using a particle-mesh solver
and cloud-in-cell interpolation \cite[][]{Guillet11}. Advection of radiation between cells is solved
using a first-order moment method, using the M1 closure method \cite[][]{Levermore84} and the Global-Lax-Friedrich
flux function to construct the radiation field at cell interfaces.

\subsubsection{Resolution and refinement strategy}

The simulation uses a dark matter particle mass of $m_{\mathrm{DM}} = 10^{4} \, \mathrm{M_\odot}$ within the zoom-in region.
The base refinement level of the mesh in the zoom-in region is $12$ (proper cell width of $1.8 \, \mathrm{pkpc}$ at $z=3$),
and it is allowed to refine up to a maximum refinement level of $19$ by $z=3$, corresponding to a minimum proper cell
size of $14 \, \mathrm{pkpc}$ at $z=3$. At $z=3$ there are $1.2 \times 10^7$ leaf cells within $R_{\mathrm{vir}}$,
the median cell width in the diffuse CGM ($0.2 < r / R_{\mathrm{vir}} < 1$, excluding satellite galaxies) is
$227 \, \mathrm{pc}$ (weighted by mass) and $454 \, \mathrm{pc}$ (weighted by volume), and the median cell
size in the ISM ($r < 0.2 R_{\mathrm{vir}}$) is $14 \, \mathrm{pc}$ (weighted by mass) and $113 \, \mathrm{pc}$
(weighted by volume).

Cells are refined either if: \textit{i)} $M_{\mathrm{DM,cell}} + \frac{\Omega_{\mathrm{M}}}{\Omega_{\mathrm{B}}} M_{\mathrm{baryon,cell}} > 8 \, m_{\mathrm{DM}}$,
where $M_{\mathrm{DM,cell}}$ and $M_{\mathrm{baryon,cell}}$ are respectively the total dark matter and baryonic mass enclosed
within a cell, \textit{ii)} the cell width is larger than one quarter of the Jeans length $\lambda_{\mathrm{J}} \equiv \sqrt{\frac{\pi c_{\mathrm{s}}^2}{G \rho_{\mathrm{gas}}}}$,
where $c_{\mathrm{s}}$ is the gas sound speed, $\rho_{\mathrm{gas}}$ is the gas mass density, and $G$ is the gravitational constant.

\subsubsection{Radiation and thermochemistry}

To represent the locally produced ultraviolet (UV) radiation field, we employ three discrete radiation bins, with frequencies bracketed by the
ionization energies of hydrogen and helium, including both the first and second ionization state. Full non-equilibrium
interactions between UV photons, hydrogen and helium are modelled as described in \cite{Rosdahl15}, including the effects
of photo-ionization, heating and momentum transfer. Radiation advection is subcycled relative to the hydrodynamical timestep,
and we use a reduced speed of light approximation with an effective advection velocity of $1/80$ times the true speed of light.
As described in \cite{Rosdahl18}, we use {B\sc{pass}}\xspace stellar evolution models \cite[][]{Elridge08} to set the input spectral
energy distribution of star particles. We use {B\sc{pass}}\xspace version 2.0 \cite[][]{Elridge16,Stanway16} which assumes
a $100 \, \%$ binary fraction. Since our zoom-in simulation targets only a single galaxy, we include the effects of a homogeneous
evolving ultraviolet background (UVB) radiation field, in addition to the radiation provided by local sources. We adopt the
model for the UVB from \cite{FaucherGiguere09}, but apply a self-shielding correction that exponentially damps the UVB intensity
for gas cells with $n_{\mathrm{H}} > 10^{-2} \, \mathrm{cm^{-3}}$ \cite[][]{FaucherGiguere10,Rosdahl12}.

In addition to the full non-equilibrium thermochemistry solved for hydrogen and helium, metal line cooling is modelled
for gas with temperature $T>10^4 \, \mathrm{K}$ using {\sc{cloudy}}\xspace \cite[][]{Ferland98}, assuming ionization equilibrium
with a \cite{Haardt96} UVB (and so is not modelled self-consistently with the non-equilibrium hydrogen and helium thermochemistry).
For $T < 10^4 \, \mathrm{K}$ we apply fine structure cooling rates from \cite{rosen95}.

\subsubsection{Star formation}

We use a model for star formation inspired by \cite{Federrath12}, as described in \cite{Kimm17}, \cite{Trebitsch17}, and
\cite{Rosdahl18}. The model accounts for the stabilising effects of both thermal pressure and turbulent motions against
gravitational collapse. The precise criteria for stars to form in a cell are that: \textit{(i)}, the local hydrogen number
density $n_{\mathrm{H}} > 10 \, \mathrm{cm^{-3}}$ (and that the local overdensity is greater than $200$ times the cosmic mean),
\textit{(ii)}, the gas is locally convergent and the host cell represents a local density maxima compared to the six nearest cells,
\textit{(iii)}, the cell width $\Delta x$ is smaller than the turbulent Jeans length given by

\begin{equation}
\lambda_{\mathrm{J,turb}} = \frac{\pi \sigma_{\mathrm{gas}}^2 + \sqrt{36 \pi c_{\mathrm{s}}^2 G \Delta x^2 \rho + \pi^2 \sigma_{\mathrm{gas}}^4 }}{6 G \rho \Delta x},
\end{equation}

\noindent where $\sigma_{\mathrm{gas}}$ is the local gas velocity dispersion,
computed using the velocity gradients to cells that share vertices with the host cell. As described in \cite{Rosdahl18},
the computation of $\sigma_{\mathrm{gas}}$ is modified with respect to the model described in \cite{Kimm17} by subtracting
rotational and symmetric divergent components of the velocity field, such that the turbulent support against a converging flow
is only provided by the relevant anisotropic motions.

For cells that meet these criteria, star particles are allowed to form stochastically with (on average) a rate given by
$\dot{\rho}_\star = \epsilon_\star \rho_{\mathrm{gas}} / t_{\mathrm{ff}}$, where
$t_{\mathrm{ff}} = [3 \pi / (32 G \rho_{\mathrm{gas}})]^{1/2}$ is the local freefall time, and $\epsilon_\star$ is
a variable efficiency that depends on the local turbulent properties of the gas and on the sound speed
\cite[see][for details]{Kimm17,Trebitsch17}. This efficiency can often reach or exceed values of unity, leading
to a bursty scenario for star formation.

\subsubsection{Stellar evolution and supernova feedback}
\label{sec_sf_fb}

We model the effects of supernova explosions using the mechanical feedback scheme described in \cite{Kimm14}
and \cite{Kimm15}; specifically we use the mechanical feedback scheme with multiple explosions
from the list of variant models considered in \cite{Kimm15}. For each star particle with age $<50 \, \mathrm{Myr}$ and $> 3 \, \mathrm{Myr}$, a stochastic sampling
is performed over the supernova delay time distribution, with individual supernova explosions
of $10^{51} \, \mathrm{erg}$ that inject energy/momentum into the host cell and neighbouring cells around the particle.
We assume a \cite{Kroupa02} stellar initial mass function, with $20 \%$ of the star particle mass
returned to the neighbouring cells across all supernova events, and in addition metals are injected with a metal
yield $y=0.075$. As discussed in \cite{Rosdahl18} and \cite{Mitchell18}, we break consistency with the Kroupa stellar initial mass function (IMF) to decide how
many supernova explosions occur for each star particle, increasing the number by a factor four such that there
are four supernova explosions (each of $10^{51} \, \mathrm{erg}$) per $100$ Solar masses of stars formed. This choice was
made by \cite{Rosdahl18} to calibrate the {\sc{sphinx}}\xspace simulations against observational constraints for $z \geq 6$
and was made by \cite{Mitchell18} to improve agreement with constraints for $z \geq 3$.

Depending on the cell metallicities and densities, the mechanical scheme determines for each cell whether the simulation is able
to resolve the adiabatic phase (in which case the injected energy is $10^{51} \, \mathrm{erg}$ and the injected
momentum scales $\propto \sqrt{m}$, where $m$ is the mass of the neighbouring cells that is considered to be
entrained into the outflow), and otherwise injects energy and the maximum momentum achieved in the
momentum-conserving snowplow phase, with a maximum injected momentum of $3 \times 10^5 \, \mathrm{M_\odot kms^{-1}}$
for gas with $n_{\mathrm{H}} = 1 \, \mathrm{cm^{-3}}$ and $Z=Z_\odot$ \cite[][]{Blondin98,Thornton98}.
Following \cite{Kimm17}, we then boost the maximum momentum to $4.2 \times 10^5 \, \mathrm{M_\odot kms^{-1}}$
if the local Stromgren radius is unresolved, accounting for unresolved photo-ionizing heating from the star particle,
which can reduce the densities around the star and so reduce the radiative losses \cite[][]{Geen15}.

\subsection{Halo finder and satellite identification}

We identify both central and satellite dark matter subhaloes using the AdaptaHOP algorithm \cite[][]{Aubert04}, and construct
merger trees following \cite{Tweed09}. We store simulation outputs with a temporal cadence of $10 \, \mathrm{Myr}$.
Halo centres are defined based on the dark matter density
maxima. All halo masses ($M_{200}$) and virial radii ($R_{\mathrm{vir}}$) quoted are measured by computing the spherical radius within
which the mean enclosed density is $200$ times the critical density of the Universe.

For each satellite identified by the halo finder, we assign particles and gas cells to that satellite if
they fall within the tidal radius, which we approximate following \cite{Binney08} as

\begin{equation}
r_{\mathrm{t}} = R_0 \left(\frac{m}{M(R_0) \left(3 - \frac{\mathrm{d}\ln{M}}{\mathrm{d}\ln(R)}|_{R=R_0}\right)}\right)^{1/3},
\end{equation}

\noindent where $R_0$ is the distance from the satellite centre to the host
centre, $M(R_0)$ is the mass of the host enclosed within $R_0$ and $m$ is the mass of
the satellite (the mass within $r_{\mathrm{t}}$). This relation is appropriate
for a satellite on a circular orbit within a spherically symmetric host potential.

\subsection{Tracer particles and classifying circum-galactic gas}
\label{tracer_methods}

The results presented in this paper make extensive use of tracer particles to provide Lagrangian information
about the trajectories of gas in our simulation. We use the Monte Carlo implementation of tracer particles
introduced by \cite{Genel13}, and adapted for use in {\sc{ramses}}\xspace simulations by \cite{Cadiou19}. \cite{Cadiou19}
show that these tracers accurately reproduce the underlying gas density field, and self-consistently follow
our subgrid models for star formation and feedback. We use a relatively high tracer sampling, with ten tracers
initialised per gas cell in the zoom-in region of the initial conditions.

One application of tracer particles is to classify the nature of the gas contained within a given cell in the
simulation. We can then use these components to (for example) decompose \mbox{Ly$\alpha$}\xspace emission from the CGM into
different components, thereby gaining insight into the CGM-\mbox{Ly$\alpha$}\xspace connection.
Unless otherwise stated, we use the following definitions to define different (mutually exclusive)
gaseous components within the virial radius of the primary host halo:

\begin{itemize}
\item \textit{ISM}: any gas with radius $r < 0.2 \, R_{\mathrm{vir}}$.
In practice this also includes material that could be argued to belong to an ISM-CGM interface, or to the inner CGM, but this is qualitatively unimportant for our results.
\item \textit{Satellite gas}: any gas that is within the tidal radius of a satellite galaxy.
\item \textit{Stripped CGM}: circum-galactic gas that was within the tidal radius of a satellite (for at least $30 \, \mathrm{Myr}$)
as it entered the virial radius of the host halo, but is subsequently located beyond the tidal radius.
This could be related to ram pressure stripping, gravitational tidal forces, or stellar feedback from the satellite; we do not attempt to distinguish between these mechanisms.
\item \textit{First-infall CGM}: circum-galactic gas that is radially infalling from the inter-galactic medium (IGM) for the first time.
We account for occasional fluctuations in the radial velocity of tracers by requiring that a tracer be radially outflowing
consecutively for $30 \, \mathrm{Myr}$ after the current snapshot in order to have left the ``first-infall'' component at the current snapshot.
\item \textit{Post-pericentric CGM}: circum-galactic gas that has moved past pericenter of the orbit associated with first infall.
This material can be outflowing (after pericenter) or inflowing (after apocenter).
Diffuse gas that has been inside the CGM for more than a halo dynamical time will almost always belong to this component,
unless it is identified as having been influenced by feedback, or was stripped from a satellite.
\item \textit{Feedback-influenced CGM}: circum-galactic gas that has been positively identified as having been likely influenced by stellar feedback in the past.
This category overrides all of the other CGM classifications.
\end{itemize}

Gas outside the virial radius of the primary host halo is always labelled as ``first-infall'', unless
it is within the virial radius of another halo (in which case it is then labelled ``satellite gas''), or was ejected from the primary host.

We estimate which tracers have been influenced by feedback (within the main halo) at a given snapshot by requiring that all of the following conditions are met:
\begin{itemize}
\item the gas is outflowing with a radial velocity (relative to the primary host halo center) of $v_{\mathrm{r}} > 50 \, \mathrm{km s^{-1}}$,
\item the gas has accelerated radially by $\Delta v_{\mathrm{r}} > 30 \, \mathrm{km s^{-1}}$ over the last $10 \, \mathrm{Myr}$,
and the gas has moved outwards by $\Delta r > 0.1 \, R_{\mathrm{vir}}$ over the same interval,
\item the radial kinetic plus thermal specific energy of the gas has increased by at least $25 \, \%$ over the same interval,
\item the tracer is not going to enter within the tidal radius of a satellite over the next $30 \, \mathrm{Myr}$.
This can be important because the halo finder often loses track of satellites within the inner regions of the host halo,
\item the tracer has \textit{not} accelerated radially outwards simply due to moving past pericenter of its (gravity-driven) orbit.
\end{itemize}

We estimate if tracers are undergoing such an orbital transition at a given snapshot by first computing the
minimum radial velocity (where inflowing is negative) of the tracer at previous snapshots, searching back
in time for as long as the tracer velocity is monotonically decreasing. For a particle that has just moved past pericenter,
this procedure identifies the maximum past infall velocity (and will return a smaller absolute value for particles that fluctuate
in velocity while settled within the ISM, corresponding roughly to the ISM velocity dispersion).
We then search for the corresponding maximum radial
velocity of the tracer for future snapshots, searching forwards for as long as the radial velocity is monotonically increasing
with time. This then returns the maximum future outflowing velocity of the particle before it starts to slow
as it approaches apocenter. For particles that are on pure infalling orbits (independent of feedback or other
energy injection mechanisms), the maximum future outflow velocity will always be similar or less than the
maximum past inflow velocity, and so for the tracer to be considered influenced by feedback we require that
the maximum future outflow velocity be greater than $3/2$ times the magnitude of the maximum past infall velocity.

This scheme is of course approximate, and was designed and tested primarily by manual inspection of many individual tracer trajectories.
In practice, feedback events will affect many more of the tracers in the simulation at some level, and there
can be other physical mechanisms that might cause a tracer to suddenly accelerate outwards.

\subsection{\mbox{Ly$\alpha$}\xspace emission and radiative transfer}

With a full radiation hydrodynamics simulation, we self-consistently account for both the photo-ionization
and photo-heating of hydrogen by local sources (plus a uniform UVB), and so we can compute the
\mbox{Ly$\alpha$}\xspace luminosity of all gas cells within the zoom-in region based on the relevant collisional excitation
and recombination rates.

For recombinations, we assume case B (see Blaizot et al., in prep for a full justification of this choice)
and follow \cite{Cantalupo08}, computing the rate of \mbox{Ly$\alpha$}\xspace photons produced by recombinations from a single
gas cell as

\begin{equation}
\varepsilon_{\mathrm{rec}} = n_{\mathrm{e}} \, n_{\mathrm{p}} \, \epsilon_{\mathrm{Ly\alpha}}^{B}(T) \, \alpha_{B}(T) \times (\Delta x)^3,
\end{equation}

\noindent where $n_{\mathrm{e}}$ and $n_{\mathrm{p}}$ are the electron and proton number densities,
$\alpha_{B}(T)$ is the case B recombination rate, $\epsilon_{\mathrm{Ly\alpha}}^{B}(T)$
is the fraction of recombinations that result in the production of a \mbox{Ly$\alpha$}\xspace photon, and $(\Delta x)^3$
is the volume of the cell. We evaluate $\epsilon_{\mathrm{Ly\alpha}}^{B}(T)$ using the fit from
\cite[][given by their equation 2]{Cantalupo08}, and $\alpha_{B}(T)$ with the fit from
\citet[][their appendix A]{Hui97}.

For collisional excitations, we compute the rate of resulting \mbox{Ly$\alpha$}\xspace photons produced from a single cell
as

\begin{equation}
\varepsilon_{\mathrm{col}} = n_{\mathrm{e}} \, n_{\mathrm{HI}} \, C_{\mathrm{Ly\alpha}}(T) \times (\Delta x)^3,
\end{equation}

\noindent where $n_{\mathrm{HI}}$ is the number density of neutral hydrogen atoms, and $C_{\mathrm{Ly\alpha}}(T)$
is the rate of collisional excitations from level $1s$ to $2p$, which we evaluate using the fit from
\citet[][their equation 10]{Goerdt10}. In addition, we set the collisional excitation rate to zero
for cells within which the net cooling timescale is less than ten times the simulation timestep,
defining the net cooling timescale as $t_{\mathrm{cool,net}} \equiv (\mathcal{C}-\mathcal{H})^{-1}$, where $\mathcal{C}$ and $\mathcal{H}$
are the cooling and heating rates. In such cells, we do not adequately resolve the balance between
heating and cooling, which may lead to a significant overestimate or underestimate of the true
collisional excitation rate. Setting the collisional excitation rate to zero means that we are computing
a conservative lower limit for these cells. We do however resolve the net cooling timescale in at least
$90 \, \%$ of the cells within the virial radius, and so we do not believe that this introduces a large uncertainty to the
collisional excitation contribution to \mbox{Ly$\alpha$}\xspace (Blaizot et al., in prep).

With the \mbox{Ly$\alpha$}\xspace emission rates from recombinations and collisional excitations, we then perform Monte Carlo
radiative transfer in post-processing to compute the \mbox{Ly$\alpha$}\xspace signal that escapes the halo, using the
{\sc{rascas}}\xspace code \cite[][]{MichelDansac20}. We sample each emission component
($\varepsilon_{\mathrm{rec}}$ and $\varepsilon_{\mathrm{col}}$) with $4\times 10^5$ Monte Carlo photon
packets cast from cells with a probability proportional to their luminosity. Each photon packet is emitted
with a frequency in the frame of the emitting cell which is drawn from a Gaussian distribution of width
fixed by the cell's temperature \citep[see Sec. 2.2 of][]{MichelDansac20}. In order to measure the
extended emission out to large distances, all cells within a radius of $11 \, \mathrm{arcsec}$ around the main galaxy are
allowed to emit. In turn, photon packets are propagated until they either pass that radius or are
absorbed by a dust grain. We have included dust following the model of \citet{Laursen09b}, assuming SMC
dust properties. In this model, the amount of dust scales primarily with the amount of neutral hydrogen
and of metals, although traces of dust will be present in the ionised phase as well. We refer to
Blaizot et al. (in prep.) for a detailed discussion of the effect of dust on the \mbox{Ly$\alpha$}\xspace properties of our
simulated galaxy. This model results in UV attenuations
(typically $\sim 1-2$ mag at $M_{1500} \sim -20$ AB) and \mbox{Ly$\alpha$}\xspace escape fractions (typically a few percent)
which are broadly consistent with what is observed for LAEs. In the scope of the present paper, the main
effect of dust is to limit the \mbox{Ly$\alpha$}\xspace flux emerging from the ISM of the central plus satellite galaxies.

The main output of RASCAS is the state of photon packets as they escape the computational domain. With this
data, we can perform part of our analysis and for example understand how different gas phases contribute to
the observed extended emission. However, this limits our analysis to angle-averaged properties of the
\mbox{Ly$\alpha$}\xspace emission and does not allow us to compare directly our simulations to MUSE observations. In order to
do so, we also use the peeling algorithm \citep{YusefZadeh84,Zheng02} as described
in \citet[][Sect. 9.3]{Dijkstra17} to collect flux in 12 mock MUSE data-cubes at each snapshot.
These 12 mocks correspond to observations of the simulated galaxy from 12 different directions (the same at
all outputs). We use them in Sec. \ref{MUSE_comp_sec} to reproduce the observational selection and analysis
of \cite{Wisotzki18}.

In order to accelerate the computation, we use the core-skipping algorithm described
in \citet[][Sect. 3.2.4]{Smith15}, which very slightly underestimates the effect of dust (introducing
a relative error of order 1\% on the escape fraction of \mbox{Ly$\alpha$}\xspace radiation) but produces a speedup of the
computation of a factor $\sim 1000$.
We do not apply any attenuation correction for the intergalactic medium, but this is generally expected
to have at most a modest $20-30 \, \%$ effect for the redshifts ($z \sim 3$) studied here \cite[e.g.][]{Inoue14,Hayes20}.

\section{Results}
\label{flows_sec}

Our results are presented as follows: we start by using tracer particles to explore the origin of the
dense gas phases of the CGM in Section~\ref{lag_orig}, we connect this to the properties of \mbox{Ly$\alpha$}\xspace emission
related to the CGM in Section~\ref{composition_subsec}, we decompose the different components
that comprise observed \mbox{Ly$\alpha$}\xspace haloes in Section~\ref{nature_obs_haloes}, we compare \mbox{Ly$\alpha$}\xspace
surface brightness profiles with stacked MUSE observations in Section~\ref{MUSE_comp_sec}, and
we finish by presenting the predicted \mbox{Ly$\alpha$}\xspace spectral morphology in Section~\ref{spectral_morphology_sec}.
Since our focus is primarily directed at the connection between \mbox{Ly$\alpha$}\xspace emission and the circum-galactic
medium, we generally exclude gas (and photons) associated with satellite galaxies for the analysis presented in
Sections~\ref{lag_orig} and \ref{composition_subsec}. We then introduce the contribution of satellites
when considering the observed properties of our simulated galaxy, in Sections~\ref{nature_obs_haloes},~\ref{MUSE_comp_sec} and ~\ref{spectral_morphology_sec}.

\subsection{Lagrangian origin of the neutral hydrogen in the CGM}
\label{lag_orig}

\begin{figure}
\includegraphics[width=20pc]{Figures/mp_radial_profiles_RhdRun-tracer_radfix.pdf}
\caption{Mass and radial momentum profiles of circum-galactic gas (excluding satellite galaxies),
plotted as a function of radius.
Data is mean-stacked over all available outputs between $z=4$ and $z=3$.
Solid lines show outflowing gas, and dashed lines show inflowing gas.
Blue shows neutral hydrogen (i.e. the mass/momentum of hydrogen in each cell
multiplied by the neutral hydrogen fraction), green shows warm-ionized hydrogen ($T<10^{4.5} \,\mathrm{K}$),
and magenta lines show hot-ionized hydrogen ($T>10^{4.5} \,\mathrm{K}$).
Most of the neutral hydrogen content of the CGM is
located within the inner halo, whereas the hot-ionized profile rises with radius,
comprising the majority of circum-galactic gas in the outer halo.
The sharp drop in neutral hydrogen mass/momentum with radius could be associated with a galactic fountain
process, or with a change in ionization state of gas as it flows inwards/outwards.
}
\label{mp_profiles}
\end{figure}

In \cite{Mitchell18}, we analysed the properties of circum-galactic gas from a sample of zoom-in simulations
at $z \approx 3$, including the halo studied here. Focussing on the neutral phase of hydrogen (which we expect
to be the closely related to the gas phases associated with extended \mbox{Ly$\alpha$}\xspace emission), we found that the mass per
unit radius in neutral hydrogen\footnote{By neutral neutral hydrogen mass we mean $0.76 \, x_{\mathrm{HI}} \, \rho_{\mathrm{gas}} \, \Delta x^3$,
where $0.76$ is the assumed primordial hydrogen mass fraction, and $x_{\mathrm{HI}}$ is the hydrogen neutral fraction
of a given cell.}
drops exponentially as a function of radius. In contrast, the mass profile of warm and hot ionized CGM components
is generally flat or rising with radius. This is shown for the simulation studied here in the top panel
of Fig.~\ref{mp_profiles}. Most of the neutral hydrogen is radially inflowing (dashed blue line), but a non-negligible
fraction ($\approx 1/3$) is outflowing (solid blue line) in the inner regions ($r < 0.5 \, R_{\mathrm{vir}}$).
Under the approximation of constant radial flow velocity, a drop in the neutral hydrogen mass with radius
implies that mass is being lost from the phase as gas flows outward, and correspondingly that the inflowing
neutral component is gaining mass as gas flows inwards. This is confirmed by the lower panel of Fig.~\ref{mp_profiles},
which shows similar gradients for the radial momentum profiles.

Understanding this behaviour is an important prerequisite for understanding in turn the predicted spatially extended \mbox{Ly$\alpha$}\xspace
signal. Our simulation predicts copious amounts of outflowing neutral hydrogen in the
inner CGM, which may be consistent with the picture of outflowing neutral gas shells that is commonly invoked
to explain the \mbox{Ly$\alpha$}\xspace spectral morphology of observed LAEs \cite[e.g.][]{Ahn04,Verhamme06}. At larger radii
($r>0.5 \, R_{\mathrm{vir}}$), there is comparatively very little neutral hydrogen in our simulation,
which poses the question of whether the outflowing neutral gas recorded in the inner CGM is being
photo-ionized as it moves outwards (joining the warm-ionized component, solid green line, which is more
radially extended than the neutral outflowing component), or alternatively is falling back towards the ISM in a fountain flow (before
reaching the outer CGM). Accordingly, the increase in inflowing neutral gas in the inner CGM could be attributed
to a change in ionization state of gas that is infalling for the first time, or could be connected
to the afore-mentioned fountain flow of neutral gas.
Another question that arises is whether the neutral hydrogen
content of the CGM in the simulation is being strongly affected by feedback.
In \cite{Mitchell18}, we show that not including stellar feedback in our simulations approximately leaves
the neutral content of the CGM unchanged,  implying that neutral circum-galactic gas flows in the simulation
are driven primarily by other processes.

\begin{figure*}
\begin{center}
\includegraphics[width=40pc]{Figures/fate_distribution_tstacks_Zoom-7-5092_RhdRun-tracer_radfix.pdf}
\caption{An overview of the Lagrangian history of circum-galactic gas in the simulation.
Each panel shows weighted distributions of tracer trajectories, for tracers that meet sets of selection criteria.
Each column corresponds to a different set of selection criteria, as indicated by the column titles.
Selections are performed either at $0.25 < r/R_{\mathrm{vir}} < 0.3$, or at $0.95 < r/R_{\mathrm{vir}} < 1$,
and include either inflowing or outflowing gas.
Individual tracers that fulfil these criteria are shifted along the x-axis such that they pass
the selection at $t = t_{\mathrm{cross}}$.
Tracer distributions are flux-weighted either by the total radial momentum at the time of selection (right-side columns),
or by the radial momentum of neutral hydrogen at the time of selection (left-side columns).
Solid black lines indicate the $16^{\mathrm{th}}$, $50^{\mathrm{th}}$, and $84^{\mathrm{th}}$ percentiles of the distributions.
Each row corresponds to a different gas quantity, including radius, neutral hydrogen fraction, temperature, density, and
the mass fraction, $f_{\mathrm{mass}}$,
belonging to the discrete and mutually exclusive Lagrangian components described in \protect{Section~\ref{tracer_methods}}.
Lagrangian components include gas that is infalling for the first time (blue),
that has been stripped from a satellite (cyan), that is within the tidal radius of a satellite (cyan),
that has been influenced by feedback (red),
that has moved past the first pericentric passage (green),
that belongs to the central ISM (orange),
that belongs to a satellite (khaki),
and tracers that are locked inside stars (black).
Distributions of gas temperature, density, and ionization fraction are truncated once $20 \, \%$ of the selected gas
has either been locked into stars (first, second columns),
or is located beyond $2 \, R_{\mathrm{vir}}$ (third, fourth columns),
at which we point we stop tracking trajectories.
}
\label{fate_neut_outflows}
\end{center}
\end{figure*}

We can now address these questions directly by using the Lagrangian information afforded by tracer particles.
Fig.~\ref{fate_neut_outflows} shows the Lagrangian history of gas selected as it crosses either
the virial radius, or a surface at $0.3 \, R_{\mathrm{vir}}$. We temporally stack simulation outputs around
the crossing redshift, such that the trajectory distributions are representative for this halo over $3.5<z<4$.
Focussing first on inflowing neutral hydrogen crossing $0.3 \, R_{\mathrm{vir}}$ (far-left column), the distribution
of tracer trajectories shows that most of the selected gas is flowing inwards for the first time from the wider
environment, as opposed to being recycled in a fountain flow. The temperature of the
inflowing gas ($T \sim 10^4 \, \mathrm{K}$, middle row) is fairly steady as the gas flows inwards down to
$r \sim 0.5 \, R_{\mathrm{vir}}$, consistent with it being maintained in thermal equilibrium by
photoheating. At the same time, the density of the gas gradually increases (fourth row) until the gas reaches
the threshold to be self-shielded from ionizing radiation at $n_{\mathrm{H}} \sim 10^{-2} \, \mathrm{cm}^{-3}$, at which
point the neutral hydrogen fraction rapidly increases (second row).
From this, we conclude that the exponential increase in the mass/momentum of inflowing neutral hydrogen with decreasing
radius seen in Fig.~\ref{mp_profiles} is caused primarily by compression of infalling gas leading to a change of ionization
state\footnote{Note however that the distributions shown in Fig.~\ref{fate_neut_outflows} are weighted by flux, which
will somewhat down-weight (relative to weighting by mass) the importance of a neutral fountain flow that is turning around at a radius comparable to $0.3 \, R_{\mathrm{vir}}$.
Fig.~\ref{mp_profiles} does indeed show a jump in both the inflowing and (more prominently) the outflowing
neutral hydrogen radial mass profiles at $r \approx \, 0.5 R_{\mathrm{vir}}$, implying that fountain flows
do indeed also play a role.}.

Shifting our attention to the neutral \emph{outflowing} hydrogen crossing $0.3 \, R_{\mathrm{vir}}$ (second-left column
in Fig.~\ref{fate_neut_outflows}), the past distribution of trajectories shows that this gas is comprised primarily
of infalling gas that has moved past first pericenter. The majority of this outflowing gas will soon move past apocentre of its
first orbit and will subsequently settle down towards the central ISM. The median average
ionization fraction stays at $x_{\mathrm{HI}} \sim 0.8$ after the selection redshift, meaning
that the sharp drop in outflowing hydrogen mass with increasing radius is primarily connected to
gas falling back down onto the ISM, rather than because of a change in ionization state.
Employing the Lagrangian tracer classifications described in Section~\ref{tracer_methods}, this component of the
gas is, as expected, associated with the ``post-pericenter'' component (green line in the bottom row), which comprises the majority
of the neutral outflow at $0.3 R_{\mathrm{vir}}$.
At the same time, a subset of the gas does move out into the outer CGM, reflected by the upper ($84^{\mathrm{th}}$) percentile in
the radius distribution (top-row).  This is associated with a growth in the
``feedback-influenced'' component (red line in the bottom row), implying that feedback is responsible for pushing
this fraction of the gas out into the outer halo.

The third and fourth columns of Fig.~\ref{fate_neut_outflows} give a broader perspective by looking at the trajectories of gas
before/after crossing the halo virial radius (in this case without selecting any specific gas phase).
Starting with inflowing gas at the virial radius (third column), the future trajectories of the gas are diverse,
with $\approx 40 \, \%$ of the gas settling into the central ISM (orange line, bottom row), $\approx 20 \, \%$ moving past pericenter
but remaining in the inner CGM (green line, bottom row), and $\approx 20 \, \%$ moving back out of the halo, which
is connected with the component that is estimated to have been influenced by feedback (red line, bottom row).
This diversity is reflected by a broad range in temperature and density, spanning three orders of magnitude
in temperature and seven orders of magnitude in density (at $300 \, \mathrm{Myr}$ after entering the halo).
It is worth noting that the inflowing gas is almost entirely ionized at the virial radius, and on average remains mostly
ionized for at least $300 \, \mathrm{Myr}$ after entering the halo (so for at least one to two halo dynamical times).

\begin{figure*}
\begin{center}
\includegraphics[width=40pc]{Figures/cgm_profiles_lagrangian_Zoom-7-5092_RhdRun-tracer_radfix.pdf}
\caption{Radial profiles of total gas mass (upper-left), neutral hydrogen mass (upper-right),
and both intrinsic (lower-left) and escaping \mbox{Ly$\alpha$}\xspace luminosity (lower-right) for circum-galactic gas
(excluding satellite galaxies).
Unlike in \protect Fig.~\ref{mp_profiles}, here the radial profiles are split between different Lagrangian CGM components,
as indicated by the legend.
Data is mean-stacked over all available outputs between $z=4$ and $z=3$.
The intrinsic \mbox{Ly$\alpha$}\xspace luminosity profile (lower-left) shows the luminosity of emitting gas, irrespective of whether photons escape the halo.
The escaping \mbox{Ly$\alpha$}\xspace luminosity profile (lower-right) shows how escaping \mbox{Ly$\alpha$}\xspace photons last trace the CGM before escape, meaning that photons are
binned at the radial position of last-scattering (or at the emission location if no scattering occurs), and are assigned to the CGM component
associated with that position.
The dashed black line in the lower-right panel shows the total intrinsic \mbox{Ly$\alpha$}\xspace luminosity profile from the lower-left panel,
for comparison.
Note that these the profiles are binned as a function of three-dimensional radius, and so the \mbox{Ly$\alpha$}\xspace profiles are not directly comparable to observed surface brightness profiles.
\mbox{Ly$\alpha$}\xspace emission is dominated by first-infalling gas in the outer CGM, both in terms of intrinsic luminosity and as the component being
traced prior to escape. \mbox{Ly$\alpha$}\xspace emission in the inner CGM traces gas associated with a mix of components.
Radiative transfer effects (i.e. scattering) enhances the escaping \mbox{Ly$\alpha$}\xspace profile relative to intrinsic by around a factor two in the CGM.
}
\label{cgm_profiles_lagrangian}
\end{center}
\end{figure*}

Considering instead outflowing gas at the virial radius (far-right column), this selection is dominated by the feedback
influenced component (red line, bottom row). The majority of this gas did not come from the ISM, but was instead accelerated
outwards while infalling (for the first time) within the CGM, with an associated increase in temperature from $T \sim 10^{4} \, \mathrm{K}$
to $10^6 \, \mathrm{K}$, consistent with being heated and entrained by hot and diffuse supernova-driven outflows.
Half of this component did not even reach $r<0.5 \, R_{\mathrm{vir}}$ before being accelerated outwards.

\subsection{The composition of the high-redshift \mbox{Ly$\alpha$}\xspace-emitting CGM}
\label{composition_subsec}

Combined together, Fig.~\ref{mp_profiles} and Fig.~\ref{fate_neut_outflows} present a picture in which most of the
the neutral (and so higher-density) phase of the CGM is undergoing gravitational infall, with the majority of
neutral inflows comprised of gas infalling before first-pericentric passage, and the majority of neutral outflowing
gas comprised of gas that has moved past first pericenter, and is in the process of settling into an equilibrium
configuration either within the inner CGM, or down in the dense ISM at the halo centre.

\subsubsection{Radial profiles}

\begin{figure*}
\begin{center}
\includegraphics[width=40pc]{Figures/phase_radial_profiles.pdf}
\caption{Distributions of neutral hydrogen fraction ($x_{\mathrm{HI}}$, top-left),
radial velocity ($v_{\mathrm{rad}}$, top-right),
hydrogen number density ($n_{\mathrm{H}}$, bottom-left),
and gas temperature ($T$, bottom-right), all plotted as a function of radius.
Gas within satellite galaxies is excluded.
Points indicate the $16$, $50$, and $84^{\mathrm{th}}$ percentiles of distributions weighted
either by mass (black), volume (red), neutral hydrogen mass (blue),
intrinsic \mbox{Ly$\alpha$}\xspace luminosity (green),
or the luminosity of \mbox{Ly$\alpha$}\xspace photons that last interact with a given gas cell before escaping (magenta).
The vast majority of escaping photons scatter at least once before escaping.
Different colour points are horizontally offset from each other for plotting purposes only, and are
all computed from a common set of five radial bins, ranging from $r=0$ to $r=R_{\mathrm{vir}}$.
Generally speaking, \mbox{Ly$\alpha$}\xspace photons are emitted and (last) scatter from partially (or near fully)
ionized and inflowing gas that is heated to slightly above $10^4 \, \mathrm{K}$.
In the CGM, \mbox{Ly$\alpha$}\xspace generally traces gas that is denser than the average (defined either by mass or volume),
but that is less dense than the neutral hydrogen component of the CGM.
}
\label{phase_profiles}
\end{center}
\end{figure*}

This picture is shown explicitly in Fig.~\ref{cgm_profiles_lagrangian}, which splits the total CGM (top-left panel)
and neutral hydrogen content of the CGM (top-right) between the various Lagrangian components described in
Section~\ref{tracer_methods}. The total gas content (top-left panel) of the outer CGM is dominated by the inflowing
gas that is infalling for the first time (``first infall'', blue line), with the contribution from gas that has moved past first
pericentric passage (``post pericenter'', green line) gradually increasing until it becomes comparable to the first-infall
component at $r \sim 0.2 \, R_{\mathrm{vir}}$. Gas that has been stripped/removed from satellites (``stripped'', cyan line)
is also an important contributor in the inner CGM (and in practice is comprised of both gas that is infalling for the
first time, and gas that has moved its first pericentric passage). Gas that has been estimated to have been clearly
influenced by stellar feedback (``feedback influenced'', red line) is subdominant at all radii, but with a flat radial
mass profile, corresponding to the approximately flat mass and momentum profiles seen for the hot-ionized outflowing phase
in Fig.~\ref{mp_profiles}.

Shifting attention to the neutral hydrogen content of the CGM (top-right panel in Fig.~\ref{cgm_profiles_lagrangian}),
the drop in neutral hydrogen mass with radius is evident for each individual Lagrangian component.
The relative contribution of stripped gas is higher for the neutral phase, and there is a greater level of parity between
the first infall, post pericenter, and stripped components, particularly within $r < 0.6 \, R_{\mathrm{vir}}$.
At larger radii, the first-infalling component increasingly dominates the neutral gas profile, although the
vast majority of the hydrogen at $R_{\mathrm{vir}}$ is in an ionized phase.

The bottom panels of Fig.~\ref{cgm_profiles_lagrangian} show the connection with \mbox{Ly$\alpha$}\xspace emission,
splitting the luminosity profile between the same Lagrangian components.
Interestingly, the intrinsic \mbox{Ly$\alpha$}\xspace luminosity emitted in-situ within the CGM (lower-left panel) appears to scale with radius
in a manner that is between the total and neutral hydrogen mass profiles.
As with the neutral hydrogen mass profile,
the intrinsic \mbox{Ly$\alpha$}\xspace luminosity decreases with radius for each component. The \mbox{Ly$\alpha$}\xspace luminosity declines less strongly
with radius however, and the relative contribution from first-infalling gas is more reminiscent
of the total mass profile.  This implies that the intrinsic \mbox{Ly$\alpha$}\xspace emission is not being
dominated by the neutral phase (despite the higher associated densities), and a significant contribution
to the total intrinsic \mbox{Ly$\alpha$}\xspace luminosity of the CGM comes from warmer ionized gas, particularly in the outer CGM.
Note that satellite galaxies are excluded from our analysis at this stage, which would otherwise dominate
the intrinsic \mbox{Ly$\alpha$}\xspace luminosity outside of the central ISM (but not the escaping luminosity).

The lower-right panel of Fig.~\ref{cgm_profiles_lagrangian} focuses on escaping \mbox{Ly$\alpha$}\xspace emission as a
tracer of the CGM. The solid black line shows the total \mbox{Ly$\alpha$}\xspace luminosity that escapes the halo after last
scattering from a given radius within the CGM (or from the radius of emission if the escaping photons do not scatter).
This (unlike the intrinsic CGM emission profile) therefore does include a contribution from photons that are emitted within the central ISM or within
satellite galaxies, but only if those photons subsequently scatter off of circum-galactic gas.
As with the other panels, the coloured lines then indicate the contribution of photons that last scatter
off of different CGM components. In general, the escaping \mbox{Ly$\alpha$}\xspace photons trace the CGM in much the same
way as the intrinsically emitted photons shown in the lower-left panel.
As a reference, the dashed black line shows the intrinsic \mbox{Ly$\alpha$}\xspace luminosity of the CGM, showing
that radiative transfer effects (i.e. scattering of photons emitted in the central ISM or
within satellites) boost the escaping \mbox{Ly$\alpha$}\xspace signal tracing the CGM out to the virial radius ($\sim 30 \, \mathrm{kpc}$)
by roughly a factor two relative to the intrinsic profile.

The gas phases traced by \mbox{Ly$\alpha$}\xspace emission are shown directly in Fig.~\ref{phase_profiles}, in which
distributions of gas properties are plotted, weighted by mass, volume, or associated \mbox{Ly$\alpha$}\xspace luminosity.
As was indicated by the shape of the radial mass/luminosity profiles shown in
Fig.~\ref{cgm_profiles_lagrangian}, \mbox{Ly$\alpha$}\xspace generally traces partially (or fully) ionized gas that
is heated to slightly above $10^4 \, \mathrm{K}$, rather than tracing the fully neutral phase
directly. \mbox{Ly$\alpha$}\xspace does still trace gas that is generally denser than the average within the CGM, defined
either by total mass or volume weighting. Escaping \mbox{Ly$\alpha$}\xspace photons (magenta points) last-scatter off gas
that is warmer and lower in density than the \mbox{Ly$\alpha$}\xspace emitting gas (green points) in the CGM.
The vast majority of escaping photons do scatter at least once.
Notably, within the ISM (far-left radial bin) escaping \mbox{Ly$\alpha$}\xspace photons last-scatter from warm, partially ionized gas (median
$x_{\mathrm{HI}} \approx 0.2$) that sits at densities above the imposed threshold for self-shielding
from the UVB ($n_{\mathrm{H}} = 0.01 \, \mathrm{cm^{-3}}$), highlighting the impact of photo-heating and photo-ionization from local radiation sources.

Putting the information in Fig.~\ref{cgm_profiles_lagrangian} and Fig.~\ref{phase_profiles} together,
we conclude that \mbox{Ly$\alpha$}\xspace is a fairly faithful differential tracer of mass in the different Lagrangian components that we
consider, albeit with declining luminosity as a function of radius due to the associated drop in average gas densities
at larger radii. That said, our analysis so far has only considered both time and angle-averaged quantities,
and in principle \mbox{Ly$\alpha$}\xspace may trace specific components of the CGM preferentially if viewed from a specific direction
or at a given time. This is pertinent, given that only a subset of high-redshift galaxies are observed to be LAEs,
and that at least for now there is an open question in observations as to whether all rest-frame UV-visible high-redshift
galaxies have associated spatially extended \mbox{Ly$\alpha$}\xspace haloes.

\subsubsection{Temporal variations}

\begin{figure*}
\begin{center}
\includegraphics[width=40pc]{Figures/cgm_mass_lum_evo_rsplit_Zoom-7-5092_RhdRun-tracer_radfix.pdf}
\caption{Time evolution in mass/luminosity of different CGM components (excluding satellite galaxies).
Different rows show total gas mass (top row), neutral hydrogen mass (second row),
intrinsic \mbox{Ly$\alpha$}\xspace luminosity (third row), and escaping \mbox{Ly$\alpha$}\xspace luminosity (bottom).
The left column shows the inner CGM over $ 0.2<r \, / R_{\mathrm{vir}} < 0.6 $,
and the right column shows the outer CGM over $ 0.6<r \, / R_{\mathrm{vir}} < 1 $.
Generally speaking, temporal fluctuations are relatively modest after $z=5$; with
escaping \mbox{Ly$\alpha$}\xspace luminosity tracing the CGM only varying by factors of a few over short timescales.
}
\label{cgm_origin}
\end{center}
\end{figure*}

Fig.~\ref{cgm_origin} shows the time dependence of different
spatially integrated CGM properties for both the inner CGM ($0.2<r\,/R_{\mathrm{vir}}<0.6$, left column) and the outer
CGM ($0.6<r\,/R_{\mathrm{vir}}<1$, right column). The total mass in the CGM (black line, top row) is
approximately constant for $z<4.5$, and generally speaking this applies to the individual Lagrangian components, with
only factors of a few variations over timescales of around $100 \, \mathrm{Myr}$.
Neutral hydrogen in the CGM varies more strongly in time (second row), particularly in the outer CGM where order
of magnitude variations are seen over $100 \, \mathrm{Myr}$ timescales. The contribution from gas stripped (or
otherwise removed) from satellites varies particularly strongly, as individual massive satellites can dominate
the signal.

Strong time variations in the intrinsic \mbox{Ly$\alpha$}\xspace luminosity of the inner CGM (third row, left column) are seen at high
redshift ($z>5$), but subsequently variations are fairly modest, with the exception of a strong feature at $z \sim 3.6$.
From visual inspection this appears to be caused by dense star-forming ISM gas briefly moving outside of $0.2 \, R_{\mathrm{vir}}$
during a satellite merger. Considering also the outer CGM, the temporal variations in intrinsic \mbox{Ly$\alpha$}\xspace luminosity are
modest at later times, particularly when compared to the mass in neutral hydrogen.

The bottom row of Fig.~\ref{cgm_origin} shows the last CGM component that \mbox{Ly$\alpha$}\xspace photons interacted with
before escaping the halo. The total escaping \mbox{Ly$\alpha$}\xspace luminosity that traces the CGM varies only modestly in time. Short term
fluctuations of about a factor two are apparent, but are partly driven by the Poisson noise in the finite distribution
of \mbox{Ly$\alpha$}\xspace photons used to perform Monte Carlo radiative transfer.
We conclude that temporal variations are not a dominant leading-order effect for \mbox{Ly$\alpha$}\xspace emission from the CGM in our simulation,
at least when averaged over all angles. This has the implication that high-redshift galaxies should generally
be surrounded by a diffuse \mbox{Ly$\alpha$}\xspace halo, irrespective of the brightness of the central \mbox{Ly$\alpha$}\xspace component (which
does vary strongly in time, Blaizot et al., in preparation).

\subsection{The nature of observed \mbox{Ly$\alpha$}\xspace haloes}
\label{nature_obs_haloes}

\begin{figure}
\includegraphics[width=20pc]{Figures/Lya_projected_SB_stacked_lagrangian_components_Zoom-7-5092_RhdRun-tracer_radfix.pdf}
\caption{Mean stacked ($3<z<4$) and angle-averaged \mbox{Ly$\alpha$}\xspace surface brightness profiles,
plotted as a function of projected separation from the halo center.
The thick solid black line in both panels shows the total escaping \mbox{Ly$\alpha$}\xspace surface brightness profile from the simulation.
The thick dashed black line shows the same profile after convolving with the typical
point spread function (PSF) inferred for deep-field MUSE observations.
The virial radius of the simulated halo at $z=3.5$ is $R_{\mathrm{vir}} = 29 \, \mathrm{kpc}$,
as marked by a small grey indicator at the bottom of each panel.
In the top panel, thin coloured lines show how escaping \mbox{Ly$\alpha$}\xspace emission traces different components of the CGM,
with each photon associated with the component it last interacted with before escape.
For the projected profiles shown here, we now also include photons that escape
from the central ISM (orange line), and from satellite galaxies
(khaki line), without interacting with the intervening CGM.
In the bottom panel, thin lines show the origin of escaping \mbox{Ly$\alpha$}\xspace emission, indicating whether
escaping photons at a given projected radius were originally emitted within
the central ISM (dashed line), in situ within the CGM (dash-dotted),
or within a satellite galaxy (dotted).
For comparison, the stacked \mbox{Ly$\alpha$}\xspace surface brightness profile of $92$ detected \mbox{Ly$\alpha$}\xspace emitters from the MUSE HDFS and UDF fields
is shown by blue crosses, taken from \protect \cite{Wisotzki18}.
}
\label{SB_profile}
\end{figure}

\begin{figure}
\includegraphics[width=20pc]{Figures/sb_profile_collrec_origin.pdf}
\caption{The fraction of the projected \mbox{Ly$\alpha$}\xspace surface brightness contributed by photons produced in recombinations
(as opposed to collisional excitations), plotted as a function of projected separation from the halo center.
As in \protect Fig.~\ref{SB_profile}, data are taken from a mean, angle-averaged stack of our simulated galaxy over $3<z<4$
(though note that no PSF convolution is applied here).
Lines show the fractional contribution of recombinations for the total surface brightness profile (solid line),
and for photons that were originally emitted inside the ISM (dashed), inside satellite galaxies (dotted), and in-situ within
the CGM (dash-dotted). Collisional excitations provide a fractionally larger contribution to the escaping
\mbox{Ly$\alpha$}\xspace flux produced in-situ within the CGM ($\approx 50 \%$), whereas recombinations contribute the majority
of the escaping flux for photons emitted within the ISM ($\approx 65 \%$).
Recombinations generally contribute $\approx 60 \%$ of the total escaping flux.
}
\label{rec_col_split}
\end{figure}

Fig.~\ref{SB_profile} shows radial surface brightness profiles of escaping \mbox{Ly$\alpha$}\xspace photons,
now plotted as a function of projected circular radius (relative to the direction of each
escaping photon). The profile is mean stacked over redshifts $3<z<4$, and over all
possible viewing angles.
In the top panel, we split the escaping \mbox{Ly$\alpha$}\xspace photons based on the Lagrangian CGM
component that each photon last interacts with before escape, using the Lagrangian components
described in Section~\ref{tracer_methods}. In addition to the CGM components discussed
previously, we now also show the contribution from photons that
last scattered on the central ISM (orange lines, defined here as
$r < 0.2 \, R_{\mathrm{vir}}$, where $r$ in this case is the three-dimensional radius), or
last-scattered within satellites (khaki lines), without interacting with the intervening CGM.

While it is expected that the central ISM dominates the production of photons observed near
the centre of the projected profile, it is notable that the majority of the photons observed
in the central region of the profile ($r_{\mathrm{projected}} < 7 \, \mathrm{kpc}$) do not
scatter in the CGM (orange line). At larger separations, the two main contributors to the signal
are the ``first infall'' gaseous component (blue line) and the component of photons
that escape from satellite galaxies (khaki line).
If the simulation picture is correct, this implies that observed \mbox{Ly$\alpha$}\xspace haloes are
tracing primarily inflowing circum-galactic gas along with emission escaping
from satellites (at least when stacked).

The lower panel of Fig.~\ref{SB_profile} splits the escaping \mbox{Ly$\alpha$}\xspace emission based on the
original source of the emission (i.e. before scattering), splitting between the ISM of the central galaxy (dashed lines), the ISM of
satellites (dotted), and in-situ photons emitted within the CGM (dash-dotted). The
largest contributor at projected separations $< 10 \, \mathrm{kpc}$ is from photons that
were emitted within the ISM of the central galaxy (note that this is before any point spread function
convolution is applied). At larger separations the escaping radiation is comprised of
a roughly equal mix of scattered photons from the central ISM, in-situ CGM photons, and photons emitted
within satellites. Over $10 < r_{\mathrm{projected}} \, / \mathrm{kpc} < 20$ the in-situ
component provides marginally the largest contribution, with satellites increasingly powering the
profile at larger separations \cite[forming an effective ``two-halo term'' outside the halo
virial radius, see also][]{Zheng11, MasRibas16, MasRibas17}, though at these radii
the stacked surface brightness has dropped significantly below the $\sim 10^{-19} \, \mathrm{erg s^{-1} cm^{-2} arcsec^{-2}}$
limit of non-stacked MUSE deep-field observations of individual LAEs \cite[e.g.][]{Wisotzki16,Leclercq17}.
Note also that after convolving with the MUSE PSF, the central ISM component is pushed outwards,
and contributes comparably to the other components out to $\sim 25 \, \mathrm{kpc}$ in projection.

An important outstanding question for the interpretation of observed \mbox{Ly$\alpha$}\xspace haloes is the nature of the
powering mechanism for the emission, which can be scattering of \mbox{Ly$\alpha$}\xspace photons escaping the ISM, emission associated with
satellite galaxies (which are likely unresolved even in deep observations of LAEs), in-situ
emission powered by ionizing UV radiation escaping the galaxy or from the UVB, or alternatively
in-situ emission powered by compressive heating associated with gravitational infall or feedback processes. While
we have addressed the role of scattering and satellite galaxies here, the powering of the in-situ
CGM emission is not trivially known even in a simulation. It is sometimes
assumed that collisional excitations in the CGM are predominantly powered by gravitational heating,
and that recombinations in the CGM are powered by fluorescence, but in practice compressive heating
of dense gas will lead to recombinations (though admittedly we expect photo-ionization to be the dominant
ionization mechanism) and photo-heating from UV photons will lead to collisional
excitations.

We leave the task of separating these effects for future work, but we do consider the
relative contribution of recombinations versus collisional excitations in Fig.~\ref{rec_col_split}.
Recombinations provide the majority of escaping photons that were emitted within the ISM
(about $65 \, \%$), but collisional excitations provide around $50 \, \%$ for the in-situ
CGM component. This implies that different physical processes are likely powering
the in-situ CGM emission relative to the ISM emission.

\subsection{Comparison to surface brightness profiles from stacked MUSE observations}
\label{MUSE_comp_sec}

Fig.~\ref{SB_profile} also compares our stacked simulation profile with the stacked
\mbox{Ly$\alpha$}\xspace surface brightness profiles measured using MUSE observations of $92$ $3<z<4$ LAEs from the Hubble Deep
field South, and the Hubble Ultra Deep Field \cite[][]{Wisotzki18}. The thick solid-dashed
line in Fig.~\ref{SB_profile} shows our simulation prediction after convolving with
the appropriate point spread function (PSF). We adopt a Moffat profile with
$\beta = 2.8$, and $\mathrm{FWHM} = 0.875 - \frac{1}{3}(\lambda \, / 10^4 \,\mathrm{\AA})$,
as inferred for the relevant deep-field MUSE observations \cite[][]{Bacon17}.

After PSF convolution the agreement between the observed and simulated stacked profiles
is generally excellent, though perhaps partially fortuitous given that we have
only simulated a single galaxy. The simulated profile does under-predict the observed
profile by up to $50 \, \%$ over the range $10 < r_{\mathrm{projected}} \, / \mathrm{kpc} < 20$,
which may be significant given that this is the range where we predict in-situ CGM emission
to play the largest role. This could be interpreted as evidence that we have insufficient
dense and/or neutral gas in the CGM of our simulated halo, possibly because our simulated galaxy is
too low in halo mass (we believe our simulated galaxy has a stellar mass which is
too large relative to its halo mass, see discussion in Section~\ref{mhalo_problem}),
or possibly because we are missing a dense/neutral phase associated with supernova-driven outflows.
Overall however, the level of agreement with observations inspires a level of confidence
in our results regarding the origin of observed \mbox{Ly$\alpha$}\xspace haloes (though see the forthcoming
discussing on \mbox{Ly$\alpha$}\xspace spectral morphologies).

Note that the simulated surface brightness profile does not exactly follow the
observational selection and stacking procedure. Specifically, the
simulated profile shown here is mean stacked and uses angle-averaged \mbox{Ly$\alpha$}\xspace
outputs. This is necessary in order to split photons into different components\footnote{
Median stacking requires the production of direction-dependent mock data cubes,
which in turn requires the use of the peeling-off technique to obtain adequate
sample statistics. With peeling-off we lose information about the origin of
individual photons however.}.
The observational profile is instead median stacked, with the motivation being
to prevent faint undetected companions from enhancing the signal in the profile outskirts.
In our simulation, we resolve $\approx 100$ satellites per projected annular bin of width
$0.5 \, \mathrm{arcsec}$, of which there are generally $\approx 20$ satellites that contain
stars \cite[see, e.g., figure 2 in][for a visual impression]{Mitchell18}. When looking at radial
profiles (not in projection) of escaping \mbox{Ly$\alpha$}\xspace luminosity, we therefore find that
while median (as opposed to mean) stacking does reduce the contribution from satellites
in the average profile (by effectively excluding the brightest satellites), the reduction is modest
(roughly a factor two in the number of the escaping photons that escape from satellites,
and a $25 \%$ reduction in the number of escaping photons that were originally emitted
within satellites before scattering). The contribution from diffuse inflowing gas
is insensitive to mean versus median stacking, and is therefore slightly up-weighted
in relative importance for a median stack (compared to the mean stack shown in Fig.~\ref{SB_profile}).

A more realistic comparison to the observed profile is shown in Appendix~\ref{ap_11_comp},
using mock data cubes orientated along specific lines of sight (rather than angle-averaging),
median stacking (matching the observed procedure), and an explicit flux selection
(rather than simply averaging all simulations outputs over $3<z<4$).
The more realistically produced profile almost exactly resembles the simulation profile shown
in Fig.~\ref{SB_profile}.

\subsection{\mbox{Ly$\alpha$}\xspace spectral morphology}
\label{spectral_morphology_sec}

\begin{figure}
\includegraphics[width=20pc]{Figures/spectra_stacked_lagrangian_components_Zoom-7-5092_RhdRun-tracer_radfix.pdf}
\caption{Stacked ($3<z<4$), angle-averaged spectra of escaping \mbox{Ly$\alpha$}\xspace photons.
The top panel shows the relative spectral morphology of the inner versus outer region, split at $5 \, \mathrm{kpc}$ in projection.
The bottom panel shows the contribution of different components, splitting photons by which galaxy/CGM component each photon last interacted
with before escaping, as indicated by the legend.
The dashed grey vertical line shown in the top panel indicates the peak shift of the \mbox{Ly$\alpha$}\xspace line predicted by the empirical relation from \protect \cite{Verhamme18}, given the full width at half maximum of our simulated spectrum for $r_{\mathrm{projected}} < 5 \, \mathrm{kpc}$.
Contrary to observed LAEs at $z \sim 3$ (for which IGM attenuation is expected to be modest), our simulated
spectrum exhibits flux bluewards of the systemic velocity that is comparable to the flux redwards of systemic.
As a consequence, the red peak of the simulated (inner region) spectrum is insufficiently redshifted relative to systemic,
given the FWHM and the \protect \cite{Verhamme18} relation.
}
\label{spectra}
\end{figure}

Fig.~\ref{spectra} shows the angle-averaged \mbox{Ly$\alpha$}\xspace spectrum of escaping photons for our simulation, mean stacked over $3<z<4$.
Since we are considering angle-averaged spectra, each \mbox{Ly$\alpha$}\xspace photon is adjusted in frequency to account for the peculiar
velocity of the halo within the cosmological box (otherwise the stacked spectrum would be artificially broadened).
The top panel of Fig.~\ref{spectra} shows the dichotomy between the spectrum of the spatial centre in projection
($r_{\mathrm{projected}} < 5 \, \mathrm{kpc}$),
and the spectrum of the spatially extended component. Most of the flux escapes from the central component, exhibiting
a double peaked spectral morphology, with a dominant blue peak slightly bluewards of the systemic velocity,
and a more offset and a subdominant red peak. In contrast, the spatially extended component is asymmetric, single peaked, and peaking
slightly bluewards of systemic.

The lower panel of Fig.~\ref{spectra} shows the relative contribution of different Lagrangian components to
the overall profile (divided on the basis of which component each escaping photon last interacts with before
escaping the system). The spectral morphology of the flux tracing the central ISM (orange line) is the broadest
component, reflecting presumably the higher neutral HI optical depth of the dense ISM, thus requiring
larger frequency excursions before escape (though there is still significant flux escaping at the
systemic velocity, implying the existence of copious low-column density sight lines). The central ISM component is
double-peaked, and near-symmetric around
systemic. \mbox{Ly$\alpha$}\xspace escaping from satellites (khaki line) is also double-peaked, but is less
symmetric. Satellites are responsible for the slight red ``bump'' feature that is apparent
for the spatially extended ($r_{\mathrm{projected}} > 5 \, \mathrm{kpc}$) component in the top panel.
All of the other CGM components in the lower
panel are slightly blue-shifted and single peaked, explaining the overall elevation of blue
relative to red flux in the total spectrum shown in the top panel.

For the CGM-tracing components, an overall blueshift is expected for an inflowing medium
\cite[e.g.][]{Dijkstra06,Verhamme06}, albeit with a fairly modest kinematic
shift due to the smaller expected absolute radial velocities of inflows compared to outflows. Typical radial inflow
velocities in our simulation are of order $50$ to $100 \, \mathrm{kms^{-1}}$ (see Fig.~\ref{phase_profiles}).
It is notable that the component tracing the gas affected by stellar feedback (red line) is
not clearly redshifted, given that feedback-driven outflows are commonly invoked to explain the
systematically redshifted spectral morphology of observed LAEs \cite[e.g.][]{Ahn04, Verhamme06}.
Note that our definition of the ``feedback-influenced'' component includes gas that was
affected by stellar feedback at any time in the past however, and is therefore comprised of
a mix of both inflowing and outflowing gas. Nonetheless, the lack of a strong redshifted
component reflects that fast-moving outflows in our simulation are generally fully ionized (see,
e.g., Fig.~\ref{fate_neut_outflows} and Fig.~\ref{phase_profiles}),
and are therefore transparent to \mbox{Ly$\alpha$}\xspace radiation.

Our simulated spectrum does not resemble the typical \mbox{Ly$\alpha$}\xspace spectrum of observed high-redshift galaxies.
Observed \mbox{Ly$\alpha$}\xspace spectra (usually measured in an aperture that is comparable in size to the central ISM)
are nearly always observed with a red peak that is substantially redshifted
relative to systemic, are asymmetric with an extended tail towards the red, occasionally show a secondary blue peak,
and often show a strong deficit of photons at the systemic redshift if the secondary peak is present
\cite[e.g.][]{Shapley03,Tapken07,Kulas12,Trainor15}.
There are observed counter-examples: for example a LAE with a double peak, a dominant blue peak and significant
emission at the systemic velocity was reported recently by \cite{Erb19}, and objects similar to our stacked spectrum
(i.e. are more symmetric around the systemic velocity) are occasionally seen
\cite[e.g.][]{Tapken07,Erb16}, but are not representative.
Recently, \cite{Hayes20} presented mean and median stacks of LAEs from MUSE observations (albeit for LAEs that are generally
brighter than the object studied here), demonstrating that any flux bluewards of systemic is completely
subdominant to a dominant redshifted component. Furthermore, using observed \mbox{Ly$\alpha$}\xspace absorption statistics,
they find that the intervening IGM has a modest $\approx 25 \%$ attenuation of emission bluewards
of the systemic velocity, and furthermore demonstrate that this is consistent with a comparison
with low-z LAEs (for which IGM attenuation is thought to be negligible).

\cite{Verhamme18} show that a positive correlation exists between the \mbox{Ly$\alpha$}\xspace FWHM, and the velocity shift of
peak from the systemic redshift of the galaxy. Given the FWHM of our stacked spectrum, we plot the predicted
\mbox{Ly$\alpha$}\xspace peak velocity ($V_{\mathrm{peak}} \approx 400 \, \mathrm{km s^{-1}}$) as the dashed grey vertical line in the top
panel of Fig.~\ref{spectra}.
While there is systematic scatter in the \cite{Verhamme18} empirical relation, this serves to underline
that even for the ISM-tracing component, we likely significantly underestimate the velocity shift of the red peak
in the spectrum.

Another consideration is the comparison of the \mbox{Ly$\alpha$}\xspace spectral morphology between the spatial center
and the spatially extended halo.
Recent spatially resolved observations of high-redshift \mbox{Ly$\alpha$}\xspace spectra
reveal significant diversity in \mbox{Ly$\alpha$}\xspace halo spectra (the halo can be both red or blue-shifted with
respect to the central region), but that the halo spectra are still significantly
redshifted from the systemic velocity of the system \cite[][]{Swinbank07,Smit17,Claeyssens19,Leclercq20},
although \cite{Erb18} present a more complex case where a double-peaked morphology is
seen in the outer halo.
These observations indicate that the unrealistic spectral morphology of our simulated galaxy
is therefore not likely confined to a problem on scales within the central ISM alone, though it is still
conceivable that the ISM is where large velocity shifts are imparted if most of the extended
emission is scattered from the CGM after being emitted within the ISM.

While it is possible that our simulated galaxy is an outlier in some way relative to the broader population, the differences
between (typical) observed spectra with our simulation more likely imply that we do not have sufficient
HI column densities, and that we are missing a strongly outflowing and volume-filling neutral gas component in our simulation.
We note that if we inspect viewing angle-specific spectra (using the peeling off technique)
instead of angle-averaged spectra, we do see a very significant diversity in the spectral
morphologies (even if viewed from different angles at a single redshift), including
rare cases where there is a dominant redshifted component. We do not however find examples of spectra with
negligible escaping emission bluewards of the systemic velocity, reinforcing the tension
between our stacked angle-averaged spectrum and observations.

\section{Discussion}
\label{discussion_section}

\subsection{\mbox{Ly$\alpha$}\xspace spectra from cosmological simulations}

Given that our simulation fails on average to produce the characteristic single peaked, redshifted from systemic, and
asymmetric spectral profile characteristic of observed LAEs at $z \sim 3$, we here review results presented
from other cosmological simulations in the literature, to assess the generality of the problem with simulations.
We focus on analyses of simulations that attempt \mbox{Ly$\alpha$}\xspace radiative transfer within the ISM, and that account (with
varying degrees of sophistication) for ionizing UV radiation from local sources.

\cite{Smith19} present an analysis of a single cosmological zoom-in simulation from the FIRE simulation suite
\cite[][]{Hopkins14} for redshifts $z>5$, with comparable resolution to the simulation presented here.
Their figure 11 shows that before IGM attenuation is applied,
their spectra also show strong blue peaks, and also high flux levels at the systemic velocity between the two
peaks, similar to the ISM-tracing component for our galaxy shown in Fig.~\ref{spectra}. As they focus
on higher redshifts, they are able to achieve reasonable looking single-peaked redshifted \mbox{Ly$\alpha$}\xspace spectra
after IGM attenuation is applied, as the IGM attenuates almost all of the flux at line centre and on
the blue side. IGM attenuation is not likely to be a sufficient explanation at lower redshifts however \cite[e.g.][]{Dijkstra07,Laursen11, Inoue14, Hayes20},
and it would be interesting to see analysis of the \mbox{Ly$\alpha$}\xspace spectra from the FIRE simulations presented
at lower redshifts.

\cite{Behrens19} analyse a single cosmological zoom-in simulation run with the {\sc{ramses}}\xspace code (without
radiation hydrodynamics). Again, as their analysis is presented at $z=7$ the IGM completely attenuates
the \mbox{Ly$\alpha$}\xspace spectrum both at, and bluewards of the systemic velocity. Before their IGM correction their
spectra appear dominated by emission from a blue peak.
\cite{Laursen07}, \cite{Laursen09} and \cite{Laursen09b} present analyses of an older and lower resolution (gas particle mass $\sim 10^5 \, \mathrm{M_\odot}$)
cosmological zoom-in simulation of a set of galaxies at $z=3.6$, and find similar results to \cite{Smith19},
in that the escaping \mbox{Ly$\alpha$}\xspace spectrum is double-peaked, with (in most cases) a stronger blue peak,
and with substantial flux emerging at the systemic velocity.

As a counter-example, while not a live cosmological simulation, \cite{Verhamme12} analysed the emerging \mbox{Ly$\alpha$}\xspace
spectra from two idealised low-mass disk galaxies. When the ISM was allowed to cool radiatively below $10^4 \, \mathrm{K}$
(their ``G2'' simulation), the emerging \mbox{Ly$\alpha$}\xspace profile is double-peaked when viewed edge-on, but is
asymmetric when viewed face-on, and with a single peak that is located redwards of the systemic velocity,
as is typical in real observations. This is attributed to outflowing gas in the simulation, but
it is not clear how the outflows in their simulation compare to the full cosmological simulations
discussed previously.
A second counter-example is provided by the recent analysis of a comparatively low-resolution cosmological
simulation by \cite{Chung19}, who do find redshifted spectra. In this case however, they utilise a subgrid
feedback model that temporarily decouples ``wind'' particles from the hydrodynamical scheme,
and ejects these particles from the ISM at high velocity. In such a scheme, supernova-driven
outflows can remain in a neutral phase even while outflowing at many hundreds of $\mathrm{km s^{-1}}$,
presumably explaining the difference in the predicted spectral morphology with the other
cosmological simulations discussed.

It can be argued that cosmological simulations produce unrealistic \mbox{Ly$\alpha$}\xspace spectral morphologies
because they do not sufficiently resolve the internal structure of star-forming clouds within the ISM,
where many of the photons are produced. Using idealised radiation hydrodynamics simulations, \cite{Kimm19} study the escape of \mbox{Ly$\alpha$}\xspace photons
from individual star-forming clouds, isolated from the effects of the surrounding diffuse ISM and
CGM. They find that radiation feedback from young stars can drive cloud expansion and
so create red-peak dominated \mbox{Ly$\alpha$}\xspace spectra, though they also show that if star formation efficiencies are
high ($10 \%$) clouds are disrupted early, leading to luminosity-weighted, time integrated
\mbox{Ly$\alpha$}\xspace spectra that are narrow and not kinematically shifted. \cite{Kimm19} also find that
they not succeed in reproducing the large velocity shifts of the red peak that are often
observed, concluding that larger-scale radiation
transfer effects in the ISM/CGM are required.

Putting all of these examples together, we conclude that we have yet to see clear evidence
that current state-of-the-art cosmological simulations are capable of producing
realistic \mbox{Ly$\alpha$}\xspace spectra, unless wind particles are explicitly decoupled from
the hydrodynamical scheme. This point has also been discussed by \cite{Gronke17},
who interpret this problem as evidence that simulations indeed are not resolving
an important population of tiny optically thick clumps, concluding that sub-resolution
models for unresolved clumping in the ISM/CGM should be implemented when post-processing cosmological
simulations with \mbox{Ly$\alpha$}\xspace radiative transfer.
The relative lack of spatial resolution in the CGM has in particular received
a great deal of attention recently \cite[][]{Hummels19,Peeples19,Suresh19,vandevoort19},
and may be particularly pertinent given that we find that most of the escaping \mbox{Ly$\alpha$}\xspace
radiation from our galaxy does not scatter from the CGM (orange line, top panel, Fig.~\ref{SB_profile}).
Another aspect that has received attention is the potential for cosmic rays to play
an important role in regulating escaping \mbox{Ly$\alpha$}\xspace spectra \cite[][]{Gronke18}, as cosmic rays can be
more resistant to dissipative losses and so more smoothly distribute the injected
feedback energy and momentum (both spatially and temporally), leading to less impulsive
but more persistent galactic-scale outflows that are better able to entrain and
maintain neutral hydrogen clouds.

At the same time, the unrealistic spectral morphology of our simulated galaxy
needs to be reconciled against the realistic \mbox{Ly$\alpha$}\xspace surface brightness profile,
which closely reproduces stacked MUSE observations
(Fig.~\ref{SB_profile} and Fig.~\ref{lya_sb_mock}). At face value, adding
copious amounts of outflowing neutral gas to the CGM of our simulation
could lead to an over-production of \mbox{Ly$\alpha$}\xspace surface brightness at larger
projected separations. We do however under-predict the surface
brightness at projected separations of $\sim 15 \, \mathrm{kpc}$
by about $50 \, \%$, leaving room to add outflowing neutral gas
at these scales to act as a scattering medium. Recently, \cite{Song20}
have demonstrated that it is possible to reproduce both the surface brightness
profiles and spectral morphology of MUSE-observed LAEs, using idealised
spherically-symmetric models that feature a continuous expanding medium extending
outwards to scales of $\sim 20 \mathrm{kpc}$ with outflowing velocities
$\sim 200 \, \mathrm{kms^{-1}}$ (their table 3). In future work we
plan to see if similar results can be obtained empirically by modifying our simulations
in post-processing such that they they contain trace amounts of neutral
hydrogen in diffuse outflowing winds.

\subsection{Producing realistic stellar distributions in simulations of high-redshift galaxies}
\label{mhalo_problem}

\begin{figure}
\includegraphics[width=20pc]{Figures/mstar_mhalo_vrot.pdf}
\caption{Basic properties of our simulated galaxy.
\emph{\textbf{Top}:} stellar mass as a function of dark matter halo mass, plotted at $z=3$.
Red points show galaxies from the fiducial RHD simulation used in this study (which uses a realistic time delay distribution for supernova explosions).
Magenta points show galaxies from a hydrodynamics-only simulation, performed without stellar feedback, using the same initial conditions \protect \cite[][]{Mitchell18}.
Black points show galaxies from a corresponding hydrodynamics-only simulation that includes stellar feedback, using a fixed time delay of $10 \, \mathrm{Myr}$ for supernova explosions.
The main target galaxy is the most massive galaxy in each simulation with $M_{\mathrm{200}} \sim 10^{11.1} \, \mathrm{M_\odot}$
(other points are from neighbouring lower-mass galaxies within the zoom-in region).
Extrapolated constraints from abundance matching are shown by the grey lines bands, and are taken from \protect \cite{Behroozi13} and \protect \cite{Moster13}.
The blue diagonal line indicates the point where stellar mass would be equivalent to halo mass multiplied by the cosmic baryon fraction.
Our simulated galaxy is uncomfortably high in stellar mass given the halo mass, with any reduction in star formation from including RHD seemingly more than offset by using
a more realistic time delay distribution function for supernova explosions (less temporally clustered supernovae apparently result in higher radiative losses).
\emph{\textbf{Bottom}:} Rotation curve ($V_{\mathrm{rot}} \equiv G \, M(r) / r$) of the main galaxy from our fiducial RHD simulation at $z=3$,
decomposed into the contribution from stars, gas and dark matter.
The rotation curve shows that the stellar distribution is highly centrally concentrated, indicative of likely radiative overcooling.}
\label{mstar_mhalo_vrot}
\end{figure}

An important caveats in our comparison to observed \mbox{Ly$\alpha$}\xspace surface brightness profiles
is the question of whether we are simulating a LAE with the correct halo mass
\cite[see for example][who predict that the in-situ \mbox{Ly$\alpha$}\xspace luminosity of the CGM strongly correlates with
halo mass]{Rosdahl12}.
The top panel of Fig.~\ref{mstar_mhalo_vrot} shows the stellar mass of our simulated galaxy
(red point at $M_{200} = 10^{11.1} \, \mathrm{M_\odot}$), compared
to extrapolations of empirical constraints from \cite{Behroozi13} and \cite{Moster13}.
With a stellar mass of $M_\star = 10^{9.5} \, \mathrm{M_\odot}$, our simulated
galaxy has an uncomfortably high stellar mass given its halo mass, sitting
an order of magnitude above the extrapolated best-fit empirical relations.

We also show the same galaxy but run without radiation hydrodynamics, both for the case of no stellar
feedback (magenta points) and with supernova feedback using a fixed time delay of $10 \, \mathrm{Myr}$ for
supernova explosions (black points). As a reminder, the fiducial RHD simulation used in this
study includes radiation feedback from local sources, boosts the injected supernova momentum
to account for unresolved photoheating \cite[][]{Geen15}, and uses a more realistic time
delay distribution function for supernova explosions. Excluding the boost associated with unresolved
photoheating, both the RhD and non-RhD simulations with feedback inject four times the supernova
energy expected for a Kroupa IMF (see Section~\ref{sec_sf_fb}).
While the RHD simulation forms about half the
mass in stars compared to the no-feedback case, it forms twice as many stars as the pure hydrodynamics
simulation with a fixed time delay. Evidently the expected reduction in stars formed due
to increased momentum input plus radiation feedback is more than offset by reducing
the temporal clustering of supernova explosions, which presumably increases the radiative cooling losses
\cite[somewhat contrary to expectations from recent work on idealised simulations, e.g., ][]{Gentry20,Keller20}.

In addition, the lower panel of Fig.~\ref{mstar_mhalo_vrot} shows the galaxy rotation curve (assuming
spherical symmetry), decomposed into contributions from stars, gas, and dark matter. The stellar component
is highly centrally concentrated, such that the resulting rotation curve is far from flat. While such
rotation curves cannot be currently ruled out for faint galaxies at $z>3$, we nonetheless interpret
this as an additional sign that the simulation likely suffers from overly high radiative losses relative to the energy
injected by feedback. We note that at this halo mass we are on the cusp of the regime where
feedback from a supermassive black hole may start to become relevant \cite[][]{Dubois15}, which may be an important
missing component needed to prevent the formation of such a compact stellar core (and could
in principle also act as a mechanism to accelerate dense gas to the velocities required to explain observed \mbox{Ly$\alpha$}\xspace spectral morphologies).

\section{Summary}
\label{summary_sec}

In this study we set out to explore the origin of dense gas in the
circum-galactic medium (CGM) at high-redshift ($z \sim 3$), and the
implications for spatially extended \mbox{Ly$\alpha$}\xspace emission from the CGM.
We have presented results taken from a radiation hydrodynamics simulation of a
high-redshift galaxy with a halo mass of $\sim 10^{11} \, \mathrm{M_\odot}$ at $z=3$.

We find that the spherically averaged radial mass profile of
circum-galactic neutral hydrogen in the simulation drops strongly with radius
(Fig.~\ref{mp_profiles}), and that
this profile is shaped primarily by the compression of cosmologically
infalling gas pushing the typical hydrogen number density above
the threshold for efficient self-shielding from ionizing radiation
(Fig.~\ref{fate_neut_outflows} and Fig.~\ref{cgm_profiles_lagrangian}).
We also demonstrated that the
(subdominant) component of neutral outflowing circum-galactic hydrogen
is primarily comprised of gas that is still undergoing gravitational
infall, but has moved past first-pericentric passage, and will subsequently
fall back towards the ISM or settle into a approximate rotational equilibrium
within the inner CGM (Fig.~\ref{fate_neut_outflows} and Fig.~\ref{cgm_profiles_lagrangian}).

We then explore the implications of this scenario for the extended \mbox{Ly$\alpha$}\xspace
haloes recently detected around faint $z \sim 3$ \mbox{Ly$\alpha$}\xspace emitting galaxies
(LAEs) with the MUSE
instrument \cite[][]{Wisotzki16,Leclercq17}. We find that \mbox{Ly$\alpha$}\xspace emitted
in-situ within the CGM roughly (but not exactly) follows the radial
profile of neutral hydrogen in the simulation (Fig.~\ref{cgm_profiles_lagrangian}),
and traces primarily dense and partially ionized infalling gas,
spanning a range of densities and ionization conditions (Fig.~\ref{phase_profiles}).
We find that angle-averaged emission from the CGM does
not vary significantly in time for our simulated galaxy, with variations
of a factor of a few at most (Fig.~\ref{cgm_origin}). This implies that high-redshift galaxies
may be generally surrounded by a diffuse \mbox{Ly$\alpha$}\xspace halo, irrespective
of the \mbox{Ly$\alpha$}\xspace escaping from the central region.

Combining in-situ emission from the CGM with emission emitted within
the central ISM and satellites, and accounting for the MUSE PSF,
we find that scattered photons from the central ISM, in-situ emission
and emission from satellites contribute comparably to the spatially
extended signal (Fig.~\ref{SB_profile}). Satellites provide much
of the very extended signal (beyond the spatial scales probed by observations
of individual LAEs, but which are accessibly via stacking).
Extended emission generally
last traces inflowing gas before escaping the halo, or otherwise escapes
from satellite galaxies without scattering from the CGM.
About $60 \%$ of the escaping photons are generated in recombinations,
with the remainder generated by collisional excitations (Fig.~\ref{rec_col_split}).
Compared to a stack of MUSE-detected LAEs from \cite{Wisotzki18},
the time and angle-averaged surface brightness profile from
our simulation reproduces the observations to at least $0.2 \, \mathrm{dex}$
over two orders of magnitude in surface brightness, out to beyond
the halo virial radius (Fig.~\ref{SB_profile}, but see also Appendix~\ref{ap_11_comp}).

Finally, we find that the average \mbox{Ly$\alpha$}\xspace spectrum produced by our galaxy
is very different from the typical spectrum of an observed LAE (Fig.~\ref{spectra}).
We find a complex spectral morphology with a dominant peak slightly bluewards
of systemic. The total spectrum is composed of a double-peaked spectral
component contributed by photons that escape from the ISM without scattering from
the CGM, and of a blue single peaked spectral component that traces inflowing
circum-galactic gas. This contrasts strongly with the spectral morphology
of typical high-redshift LAEs in observations; the spectra of observed LAEs
at $z \sim 3$ are
usually significantly redshifted from the systemic velocity,
are asymmetric with a broad redward tail, and have minimal observed
flux bluewards of the systemic velocity (even accounting for IGM attenuation).

This problem seems to apply to other cosmological simulations
before IGM attenuation is applied \cite[e.g.][]{Behrens19,Smith19},
but since in our case we analyse
emission at redshifts $3<z<4$, IGM attenuation is unlikely to be a sufficient
explanation. We speculate that this is indeed a general problem
with current simulations
\cite[see Section~\ref{discussion_section} as well as the discussion in][]{Gronke17},
and that the problems could be related to insufficient resolution in the CGM,
and/or because cosmological simulations do not (usually) include cosmic
ray creation and transport.

We conclude that observations of spatially resolved \mbox{Ly$\alpha$}\xspace emission
from high-redshift galaxies are providing highly challenging
constraints for state-of-the-art cosmological simulations.
These constraints are particularly timely given the current
debate regarding the sensitivity of the phase distribution of
the simulated CGM to numerical resolution, and the question
of whether cosmic rays play a significant role in shaping
the ISM and CGM of galaxies. With the upcoming launch of JWST
set to constrain whether high-redshift LAEs are surrounded by extended
haloes of Balmer-line emission, it will soon be possible to
establish whether in-situ emission from the CGM provides
a significant contribution to the extended signal.
If confirmed, and given that the extended emission is generally observed to be red-shifted
from systemic \cite[][]{Swinbank07,Smit17,Claeyssens19,Leclercq20}, the implication would then be that the red-shifted
spectral morphology of high-redshift \mbox{Ly$\alpha$}\xspace lines is indeed being shaped by
a fast-moving and neutral phase of CGM-scale (as opposed to only
ISM-scale) galactic outflows
with a high covering fraction, a component that is seemingly
missing from current state-of-the-art cosmological simulations.

\section*{Acknowledgements}

We would like to thank Thibault Garel, Taysun Kimm, Joakim Rosdahl, and Anne Verhamme for providing comments and suggestions
that improved the quality of this manuscript.
We are grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its financial support within the program ``Investissements d'Avenir'' (ANR-11-IDEX-0007) of the French government operated by the National Research Agency (ANR).
This work was supported by Vici grant 639.043.409 from
the Netherlands Organisation for Scientific Research (NWO).
Simulations were run at the Common Computing Facility (CCF) of LABEX LIO.
CC acknowledges support from the Institut Lagrange de Paris and from
the European Union Horizon 2020 research and innovation programme under
grant agreement No. 818085 GMGalaxies.

\section*{Data availability}

The data underlying this article will be shared on reasonable request to the corresponding author.

\bibliographystyle{mn2e}

redpajama_set_name: RedPajamaArXiv