<figure><img src="image_1.png" /></figure>

# Constraining Lyman-alpha spatial offsets at $3 < z < 5.5$ from VANDELS slit spectroscopy

A. Hoag<sup>1</sup>,<sup>2</sup> T. Treu,<sup>1</sup> L. Pentericci,<sup>2</sup> R. Amorin,<sup>3,4</sup> M. Bolzonella,<sup>5</sup> M. Bradač<sup>1,6</sup>,<sup>6</sup>

M. Castellano,<sup>2</sup> F. Cullen,<sup>7</sup> J. P. U. Fynbo,<sup>8</sup> B. Garilli,<sup>9</sup> L. Guaita,<sup>2,10</sup> N. Hathi<sup>1,11</sup>

A. Henry,<sup>11</sup> T. Jones<sup>1,6</sup> C. Mason<sup>1,12</sup>† D. McLeod,<sup>7</sup> R. McLure,<sup>7</sup> T. Morishita,<sup>11</sup>

L. Pozzetti,<sup>5</sup> D. Schaerer,<sup>13,14</sup> K. B. Schmidt,<sup>15</sup> M. Talia<sup>5,16</sup> and R. Thomas<sup>17</sup>

*Affiliations are listed at the end of the paper*

Accepted 2019 June 24. Received 2019 June 16; in original form 2019 May 23

## ABSTRACT

We constrain the distribution of spatially offset Lyman-alpha emission (Ly α) relative to rest-frame ultraviolet emission in ~300 high redshift ($3 < z < 5.5$) Lyman-break galaxies (LBGs) exhibiting Ly α emission from VANDELS, a VLT/VIMOS slit-spectroscopic survey of the CANDELS Ultra Deep Survey and Chandra Deep Field South fields (≈0.2 deg² total). Because slit spectroscopy only provides one spatial dimension, we use Bayesian inference to recover the underlying two-dimensional Ly α spatial offset distribution. We model the distribution using a two-dimensional circular Gaussian, defined by a single parameter $\sigma_{r,\text{Ly}\alpha}$, the standard deviation expressed in polar coordinates. Over the entire redshift range of our sample ($3 < z < 5.5$), we find $\sigma_{r,\text{Ly}\alpha} = 1.70_{-0.08}^{+0.09}$ kpc (68 per cent conf.), corresponding to ~0.25 at $\langle z \rangle = 4.5$. We also find that $\sigma_{r,\text{Ly}\alpha}$ decreases significantly with redshift. Because Ly α spatial offsets can cause slit losses, the decrease in $\sigma_{r,\text{Ly}\alpha}$ with redshift can partially explain the increase in the fraction of Ly α emitters observed in the literature over this same interval, although uncertainties are still too large to reach a strong conclusion. If $\sigma_{r,\text{Ly}\alpha}$ continues to decrease into the reionization epoch, then the decrease in Ly α transmission from galaxies observed during this epoch might require an even higher neutral hydrogen fraction than what is currently inferred. Conversely, if spatial offsets increase with the increasing opacity of the intergalactic medium, slit losses may explain some of the drop in Ly α transmission observed at $z > 6$. Spatially resolved observations of Ly α and UV continuum at $6 < z < 8$ are needed to settle the issue.

**Key words:** galaxies: high-redshift.

## 1 INTRODUCTION

The Lyman-alpha (Ly α; rest frame 1215.7 Å) emission line has been used as a beacon to spectroscopically confirm the redshifts to the most distant galaxies for decades now (e.g. Steidel et al. 1996; Shapley et al. 2003; Stark et al. 2010; Finkelstein et al. 2013; Schenker et al. 2014; Zitrin et al. 2015). It is well suited for this task for two main reasons. First, it is typically the strongest emission line in the rest-frame UV spectra of galaxies. Second, at $z > 2$ it is redshifted in the optical/near-infrared where detector quantum efficiency is high, sky backgrounds are relatively low, and it is often the only line accessible for spectroscopic confirmation.

In addition to its detectability out to high redshift, Ly α experiences resonant scattering with neutral hydrogen gas. As a result, it has been used as a diagnostic to probe the state of cosmic reionization (e.g. Haiman & Spaans 1999; Malhotra & Rhoads 2004; Fontana et al. 2010; Stark et al. 2010; Treu et al. 2012). Evidence from the detection of Gunn & Peterson (1965) troughs in quasars (Becker et al. 2001; Fan et al. 2006) suggests that reionization is completed by $z \sim 6$, and ongoing at $z > 6$. If this is the case, then it is expected that the fraction of Lyman-break galaxies (LBGs) exhibiting strong Ly α should decline during reionization (e.g. Fontana et al. 2010; Stark et al. 2010). This was in fact found to be observationally true (Pentericci et al. 2011; Ono et al. 2012; Schenker et al. 2012; Treu et al. 2012, 2013; Pentericci et al. 2014; Schenker et al. 2014; Tilvi et al. 2014; Mason et al. 2018a; Hoag et al. 2019), providing additional evidence for the onset of reionization at $z > 6$. The decline in Ly α fraction is especially

<footnote>* E-mail: athoag@astro.ucla.edu</footnote>

<footnote>† Hubble Fellow.</footnote>
significant for fainter galaxies ($M_{UV} \gtrsim -20.25$) (e.g. Pentericci et al. 2014).

While the drop in Ly α fraction at $z > 6$ is potentially compelling evidence for reionization, especially given the complementary evidence from quasars, many argue that reionization is not the only possible explanation. For example, evolution in Lyman-continuum escape fraction (Dijkstra et al. 2014) or increasing numbers of absorption systems present at the end of reionization (Bolton & Haehnelt 2013) provide potential alternative explanations. Regardless of the reason(s) for the decline, Haiman (2002), Santos (2004), Dijkstra, Mesinger & Wyithe (2011), and Mesinger et al. (2015) found that the velocity offset imparted on Ly α in the interstellar medium (ISM) and circumgalactic medium (CGM) strongly affects the transmission of Ly α through the intergalactic medium (IGM). For example, Ly α lines with large velocity offsets ($\gtrsim 200$ km s$^{-1}$) are less attenuated by the IGM. While this is not necessarily evidence against reionization as an explanation of a declining Ly α, it suggests that the interpretation is complicated due to ISM and CGM physics.

A less explored explanation for the drop in Ly α fraction at $z > 6$ is differential slit losses from spatial variations in the distribution of Ly α emission relative to the UV continuum. This is relevant because most of the Ly α fraction measurements are made using slit spectroscopy, as this is currently the most efficient way to probe the $z > 6$ Universe spectroscopically. Ly α radiative transfer in the ISM and CGM is known to affect both the spectral and spatial distribution of the line. In particular, Laursen & Sommer-Larsen (2007), Zheng et al. (2011), and Dijkstra & Kramer (2012) showed using theoretical models that scattering in the ISM and CGM can produce Ly α haloes an order of magnitude larger in size than the rest-frame UV. Observational evidence of extended Ly α haloes around galaxies first came from narrow-band imaging (e.g. Møller & Warren 1998; Swinbank et al. 2007; Nilsson et al. 2009; Finkelstein et al. 2011). The ubiquity of Ly α haloes was later convincingly shown first from stacks of LBGs (Hayashino et al. 2004; Steidel et al. 2011) and later in individual galaxies (Rauch et al. 2008; Wisotzki et al. 2016; Leclercq et al. 2017; Wisotzki et al. 2018).

Spatially offset Ly α is less well understood. 3D models of Ly α radiative transfer indicate that Ly α escape is strongly dependent on the inclination angle (e.g. Laursen & Sommer-Larsen 2007; Verhamme et al. 2012; Behrens & Braun 2014) in systems with discs. One explanation for this is that Ly α escapes more readily perpendicular to discs as opposed to through them. Another scenario that could lead to spatially offset Ly α is the accretion of less enriched satellites undergoing bursts of star formation (Ly α) on to a more enriched host (UV). Ly α spatial offsets have been sporadically reported in the literature. For example, Bunker, Moustakas & Davis (2000) demonstrated a convincing Ly α spatial offset relative to the rest-frame UV continuum in a longslit observation of a bright lensed galaxy at $z = 4$. After correcting for lensing, the spatial offset is $\sim 1$ kpc. Similarly, Fynbo, Møller & Thomsen (2001) found a $\sim 4$ kpc spatial offset at $10\sigma$ significance using narrow-band imaging. In a more comprehensive narrow-band imaging survey, Shibuya et al. (2014) studied a large sample of $z \sim 2.2$ Ly α emitters, finding statistically significant offsets as large as $\sim 4$ kpc ($\sim 0.5$ at $z = 2.2$). The authors do not quantify the frequency or size distribution of such offsets, and they also provide the caveat that some of their offsets are most likely due to mergers.

The advent of the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2014) on the Very Large Telescope (VLT) has made detailed spatially resolved Ly α spectroscopy of individual galaxies

at $z > 3$ possible, without the stringent constraints on the redshift from narrow-band imaging. Using MUSE, Wisotzki et al. (2016) found that nearly all of the 26 Ly α emitting galaxies in their sample have an extended Ly α halo which is $\gtrsim 5–15$ times larger than their rest-frame UV continuum size, and $\sim 5$ times larger than Ly α haloes measured in the local Universe (Hayes et al. 2013; Guaita et al. 2015). However, Wisotzki et al. (2016) were unable to reliably measure spatial offsets between Ly α and the rest-frame UV due to the large astrometric uncertainty in their HST-MUSE registration, even with an ultradeep exposure. Leclercq et al. (2017) did observe Ly α spatial offsets up to $\sim 0.3$ in MUSE observations, but they were more interested in constraining the extent of the Ly α haloes. Astrometry issues aside, MUSE has a limited FOV ($\sim 1$ arcmin$^2$), so a statistical measurement of Ly α spatial offsets would require many MUSE pointings. The MUSE-Wide survey (Urrutia et al. 2019) is currently in progress to obtain 100 MUSE pointings at 1 h per pointing. Assuming the astrometric issues mentioned by Wisotzki et al. (2016) can be overcome, this may be a promising avenue to constrain Ly α spatial offsets.

In order to determine the relative importance of the different phenomena that could be causing the drop in Ly α fraction, it is crucial to first establish the evolution of LBG/Ly α properties at $z < 6$, i.e. after reionization is complete. At redshifts $3 < z < 6$, Stark et al. (2010), Stark, Ellis & Ouchi (2011), Hayes et al. (2011), Curtis-Lake et al. (2012), and Cassata et al. (2015) found that the Ly α fraction actually increases, possibly due to decreasing fractions of dust and neutral absorbing gas in the ISM (e.g. Finkelstein et al. 2012; Jones et al. 2013). A recent study with MUSE, however, found no evidence for an increase in Ly α fraction over $3 < z < 6$ (Caruana et al. 2018). The discrepancy in these results may be due to the fact that MUSE integrates the light over the true Ly α spatial profile, while slits may miss spatially offset or extended emission. Perhaps the evolution in the Ly α fraction measured by previous authors with slit spectroscopy is due to an evolution in the morphology of Ly α emission which would make it more easily observed in slits at higher redshifts. If this were true, then narrow-band imaging would find a flat Ly α fraction. Ouchi et al. (2008) investigated the evolution of the Ly α fraction above an equivalent width of 240 Å from narrow-band imaging at $z = 3.1, z = 3.7$, and $z = 5.7$. The authors found a tentative increase in Ly α fraction, but the results were statistically consistent with no evolution.

In this work, we use slit spectroscopy of a large sample ($\sim 300$) of LBGs exhibiting Ly α in emission to constrain the distribution of Ly α spatial offsets at high redshift ($3 < z < 5.5$). We aim to understand the impact of slit losses due to these spatial offsets on the interpretation of current and future Ly α fraction observations. In Section 2, we summarize the data sets that are used in this work. We describe our method for measuring Ly α spatial offsets in Section 3. In Section 4, we describe our Bayesian inference methodology used to recover the physical Ly α offset distribution. We apply our inference to constrain the offset distribution in Section 5 and discuss our findings and their context in Section 6. We summarize in Section 7.

We adopt a concordance cosmology with $\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$, and $h = 0.7$. All magnitudes are reported in the AB system, and all physical distance measurements are in proper kpc, unless otherwise specified.

## 2 DATA

Here we describe the spectroscopic and imaging data that we used in this work.
## 2.1 Spectroscopic data

The primary data used in this work are spectra from the VANDELS survey (McLure et al. 2018; Pentericci et al. 2018b). VANDELS is a deep optical spectroscopic survey with the VIMOS spectrograph on the VLT. VANDELS targeted two fields: one centred on the UKIDSS Ultra Deep Survey (Almaini et al. in preparation; UDS: 02:17:38, −05:11:55) and the other centred on the Chandra Deep Field South (CDFS: 03:32:30, −27:48:28). Both fields have high-quality ancillary multiwavelength data. All spectra were obtained with the VIMOS medium resolution grating. Details on the overview of the survey strategy and target selection can be found in McLure et al. (2018) while observations and data reduction can be found in Pentericci et al. (2018b).

The VANDELS team produces 2D (1 spatial axis, 1 dispersion axis) and 1D spectra (dispersion axis only) along with a catalogue of spectroscopic redshifts ($z_{\text{spec}}$) of varying quality. The team determined redshifts via the EZ software<fnref n="1" /> (Garilli et al. 2010). The software allows one to look simultaneously at the 1D, 2D sky, and signal-to-noise ratio (S/N) spectra, and also at the HST imaging thumbnails. In the case of Ly $\alpha$ emission that is spatially offset from the UV continuum, this line will appear in the 2D spectra and will still be identified by the inspectors. Qualities were assigned to the spectra based on the criteria outlined by Le Fèvre et al. (2005), using two independent human inspectors. In brief, each inspector assigned a quality of 0, 1, 2, 3, 4, or 9, where 0 means redshift was not able to be measured, 1–4 means the confidence in the redshift measurement was 50 per cent, 75 per cent, 95 per cent, or 100 per cent, and 9 means the redshift was assigned based on a single emission line (Pentericci et al. 2018b). The two inspectors were required to come to an agreement on the quality of each redshift after independently grading each spectrum. A final double check on the flag was done independently by the two co-PIs (see McLure et al. 2018).

The VANDELS team provide spectra of reference point sources on each mask, but no estimate of the seeing. Because the full-depth spectra were often observed on multiple masks from different nights, the seeing for each object must be calculated individually as a combination of the seeing from the different masks. The seeing we used for each target was the median of the seeing calculated on all of the masks on which it was observed. For each individual mask, which often contained multiple reference point sources, we calculated the seeing by taking the median of the seeing calculated from each reference point source spectrum. For each reference source, we fit a 1D Gaussian to the spatial profile of the spectral continuum. We found that the centroid of the continuum varied significantly with wavelength, so we fit 1D Gaussians in 28 bins of 75 pixels (~190 Å) each, calculated the standard deviation and then took the median of all standard deviations. The bin size was chosen to balance sufficient S/N for centroiding with the ability to measure wavelength dependence.

## 2.2 Imaging data and VANDELS target selection

The construction of the photometric catalogues used for target selection are described in detail by McLure et al. (2018). To briefly summarize, the VANDELS footprint within the UDS and CDFS fields covers both the central areas which have deep HST imaging as well as the wider areas where only shallower ground-based imaging

are available. As a result, a distinct photometric catalogue is used for each of the four regions: UDS-HST, UDS-GROUND, CDFS-HST, and CDFS-GROUND. The two regions with HST coverage employ the H-band selected catalogues provided by the CANDELS collaboration (Galametz et al. 2013; Guo et al. 2013). The VANDELS team produced photometric multiwavelength catalogues for the UDS-GROUND and CDFS-GROUND regions as there were no publicly available multiwavelength catalogues for these regions. The CDFS-GROUND images had variable seeing and were PSF-homogenized to a seeing of 1.0 arcsec full width at half-maximum (FWHM). The UDS-GROUND images had stable seeing so PSF-homogenization was unnecessary. H-band selected catalogues were produced using 2 arcsec diameter circular apertures for photometry. UDS-GROUND (CDFS-GROUND) spans 12 (17) filters from the U band to the K band. For full details on the photometry see Mortlock et al. (2017) and McLure et al. (2018).

For target selection, the VANDELS team made use of the photometric redshifts provided by the CANDELS team for the UDS-HST and CDFS-HST regions (Galametz et al. 2013; Santini et al. 2015). For the UDS-GROUND and CDFS-GROUND regions, the VANDELS team generated their own photometric redshifts. The photometric redshifts were derived by taking the median best-fitting value from 14 different redshift codes using a broad range of spectral energy distribution (SED) templates, star formation histories, metallicities, and emission-line prescriptions. The photometric redshifts from each code were tested and validated against previous spectroscopic redshift data sets from the 2 GROUND field regions, e.g. 3D-HST (Brammer et al. 2012; Momcheva et al. 2016), UDSz (Almaini et al., in preparation) in UDS and Le Fèvre et al. (2005), Vanzella et al. (2008), and Momcheva et al. (2016).

After removing potential stellar sources from the catalogues, the VANDELS team performed SED fitting on all sources in the four fields using the Bruzual & Charlot (2003) templates with solar metallicity, no nebular emission, exponentially declining star formation histories, and the Calzetti et al. (2000) dust attenuation law. For more details see McLure et al. (2018).

From this photometric sample in the four regions, potential spectroscopic targets were selected by the VANDELS team to be in these main categories (McLure et al. 2018):

(i) Bright star-forming galaxies in the range $2.4 \le z \le 5.5$

(ii) LBGs in the range $3.0 \le z \le 7.0$

(iii) Passive galaxies in the range $1.0 \le z \le 2.5$

Because Ly $\alpha$ is only observable at 2.99 < $z_{\text{spec}}$ < 7.38 with VIMOS, our final spectroscopic sample (see Section 3 for a description of the spectroscopic sample selection) almost entirely consists of objects in categories (i; 38/305) and (ii; 266/305). The one remaining object in our sample came from a sample of Herschel detected objects supplied by D. Elbaz, which was provided after the initial object selection.

The apparent H-band magnitude distribution of the sources in our final spectroscopic sample is well described by a Gaussian distribution with $\mu_H = 25.2$ mag and $\sigma_H = 0.6$ mag. We determined the absolute magnitude of the sources in our final sample using the H-band magnitude and the $z_{\text{spec}}$ recorded by the VANDELS team. The absolute magnitude distribution computed in this way is similarly well described by a Gaussian with $\mu_{M_{\text{UV}}} = -20.6$ mag and $\sigma_{M_{\text{UV}}} = 0.6$ mag.

<footnote n="1">¹http://pandora.lambrate.inaf.it/docs/ez/</footnote>
## 3 LY $\alpha$ SPATIAL OFFSET MEASUREMENTS

To constrain the Ly $\alpha$ spatial offset distribution, we first assembled a catalogue of all galaxies showing Ly $\alpha$ in emission in the VANDELS data base.<fnref n="2" /> We queried the data base for the UDS BEST SPECTRA and CDFS BEST SPECTRA, which represent the full-depth co-added spectra from the UDS and CDFS on 2019 January 8. We note that the slit orientation was always the same for each science target when observed on different masks, such that co-adding the masks does not affect the measured Ly $\alpha$ offset, if present. Care was taken during the reduction stage to ensure miscentring did not occur when stacking spectra from multiple masks. We filtered the output by: $2.95 < z_{\text{spec}} < 8$, which resulted in 611 UDS spectra and 658 CDFS spectra. This filtering step ensured that all of the spectra we downloaded had an assigned spectroscopic redshift, $z_{\text{spec}}$. The medium resolution VIMOS wavelength coverage is sensitive to Ly $\alpha$ at $2.99 < z_{\text{spec}} < 7.38$, but we include a wider range in our filtering step in case the redshift assignment was slightly incorrect. Including this extra range will not affect our results as explained below. We do not perform a filter on redshift quality assigned by the VANDELS team because of the case where a large Ly $\alpha$ spatial offset could have been misinterpreted as coming from another source. As we point out below, the redshift qualities of the spectra in our final sample turn out to be primarily (~95 per cent of them) $Q=3$ and $Q=4$, i.e. $\gtrsim$ 95 per cent confidence in the $z_{\text{spec}}$.

Using the $z_{\text{spec}}$ for each of our downloaded spectra, we visually searched the 2D spectrum for Ly $\alpha$ emission within a ~30 Å window of the predicted Ly $\alpha$ wavelength. We chose a large spectral window to ensure that we did not miss Ly $\alpha$ that was offset in velocity from the $z_{\text{spec}}$ reported in the catalogue. Within this spectral window, we searched the entire spatial extent of the spectrum for Ly $\alpha$ so that we could detect spatially offset Ly $\alpha$ emission. We also produced a collapsed spatial profile within this window to visually inspect. After inspecting the 2D spectrum and the collapsed spatial line profile, we flagged each spectrum as either having Ly $\alpha$ emission or not. We flagged 426 (194 in UDS, 232 in CDFS) targets as having Ly $\alpha$. While we were very inclusive in our visual inspection, in the following steps we removed low S/N and spurious features from our selection.

To obtain the Ly $\alpha$ spatial centroid, we first found the optimal spectral line centre by collapsing the spectra along the spatial axis and fitting the resulting spectrum to a 1D Gaussian. We used the Ly $\alpha$ wavelength inferred from the catalogue $z_{\text{spec}}$ as the wavelength prior for this step. We then produced a Ly $\alpha$ spatial profile by collapsing the 2D spectrum along the spectral axis in a 30 Å window centred on the optimal line centre we found in the previous step. For this step, we used the imaging catalogue position as the spatial prior.<fnref n="3" />

While the imaging catalogues provide a spatial continuum centroid in the spectrum, we cannot use this centroid to compare to the Ly $\alpha$ spatial centroid when calculating the offset between Ly $\alpha$ and the UV continuum. This is because there is a noticeable drift in the continuum spatial centroid in many of the spectra. Examples of the drift are shown in Fig. 1. The drift is likely due to atmospheric refraction, which varies with the airmass of the observations. Because this effect was not corrected for during the reduction, the full-depth reduced spectra, which consist of many sets of exposures taken at varying airmass values, have a blend of

spectra with varying degrees of the distortion. As a result, in any given spectrum the spectral continuum centroid may be spatially shifted relative to the expected imaging catalogue position, and this shift varies with wavelength. Because the drift is wavelength dependent, it affects the continuum centroid and the Ly $\alpha$ centroid, and it could result in an artificial spatial offset between the two.

In order to reliably measure true Ly $\alpha$-UV spatial offsets in our spectra, we need to be able to calculate and remove the drift. The only constraint we have on the drift is the continuum spatial centroid, which is independent of wavelength in the absence of the drift. We specifically do not use the Ly $\alpha$ spatial position to estimate the drift because if Ly $\alpha$ is truly spatially offset from the continuum then using it could bias our drift measurement and therefore our spatial offset measurement. We fit a second-order polynomial to the continuum centroid measured in bins of 50 Å, over a bandpass of 1500 Å, starting 35 Å redward of the optimized Ly $\alpha$ wavelength. We found that $S/N_{1D} \ge 2$ per pixel was required for a sufficiently accurate fit to the continuum. We also required an integrated Ly $\alpha$ S/N ratio, $S/N_{\text{int, Ly }\alpha} \ge 5$ to accurately measure the Ly $\alpha$ spatial centroid. We used the 1D signal and noise extensions of the VANDELS data products to calculate $S/N_{1D}$ and $S/N_{\text{int, Ly }\alpha}$. These two S/N requirements resulted in 320 objects from the 426 which we flagged as having Ly $\alpha$. We visually inspected these objects to assess the polynomial fitting. 15 of the 320 objects had poor fits to their continuum centroid, either due to the presence of continuum from other objects falling serendipitously in the slits or from spectral artefacts. We removed those 15 objects from our sample, resulting in a final spectroscopic sample of 305 objects (131 in UDS and 174 in CDFS).

For the 305 objects in our final sample, we calculated the 1D spatial offset as $y_{\text{Ly }\alpha, \text{opt}} - y_{\text{cont, opt}}$, where $y_{\text{Ly }\alpha, \text{opt}}$ is the mean of the Gaussian fit to the spatial profile of the line and $y_{\text{cont, opt}}$ is the value of the second-order polynomial fit to the UV-continuum, evaluated at the best-fitting Ly $\alpha$ wavelength. We extrapolated the fit to the Ly $\alpha$ wavelength because there could still be a small drift between the bluest wavelengths used in our continuum bandpass and the Ly $\alpha$ wavelength. This polynomial and its extrapolated spatial position are shown for a few example spectra in Fig. 1. Typical drifts are small (≤1 pixel per 500 Å) but in some cases they exceed this. The polynomial fit performs well for the range of drifts we observe in the spectra. By using $y_{\text{cont, opt}}$, we correct for the spectral drift, which allows us to reliably calculate an accurate Ly $\alpha$-UV spatial offset.

The distribution of 1D spatial offsets measured from the 305 spectra is shown in Fig. 2. We note that the maximum spatial offset we can measure in the spectra is set by the slit length and the position of the target within the slit. VANDELS used variable slit lengths, but with a minimum of 7 arcsec, or $\gtrsim$40 kpc over the entire redshift range probed in this work, and objects are positioned away from the edges of the slit. Given that the maximum spatial offset that we observe is <10 kpc in magnitude, it is extremely unlikely that we missed spatial offsets due to Ly $\alpha$ falling outside of the slit along the length of the slit.

Some objects show significantly offset Ly $\alpha$ emission. For objects with the largest spatial offsets (>2 pixels ≃ 0.4 ≃ 2.8 kpc in magnitude at $z=4$), we inspected their continuum image thumbnail (from the 'THUMB' extension of the object fits file) with the slit overlaid to ensure that the emission line was not coming from another object at a different redshift or perhaps an interacting galaxy at the same redshift. The CDFS thumbnails are in the *r* band, while the UDS thumbnails are in the *i* band. None of these offsets appeared to be coming from a nearby source. We also checked whether the objects with large offsets came primarily from the ground-based

<footnote n="2">²http://vandels.inaf.it/dr2.html</footnote>

<footnote n="3">³This imaging catalogue position is saved under the 'HIERARCH PND WIN_OBJ_POS' keyword in the header of the image extension of each downloaded spectrum fits file.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of how we compute Ly α spatial offsets from the VANDELS spectra. Top row: A bonafide spatially offset Ly α-emitting galaxy. Left: r-CDFS (or *i* (UDS)-band Image of the galaxy showing a portion of the slit overlaid (red vertical lines). Horizontal dotted lines are: (black) continuum position estimated from the polynomial fit (see right-hand panel), (blue) the best-fitting Ly α centroid, and (red) the catalogue spatial position which is shown at y = 0. Right: 2D spectrum from the slit to the left showing the strong Ly α emission and rest-frame UV continuum, with a clear spatial offset between the two. The continuum emission centroid drifts, so we fit it with a second-order polynomial (black line), which we extrapolate to estimate the continuum position at the Ly α wavelength (dotted black line in left-hand panel). The top two rows are examples of large Ly α-UV spatial offsets, whereas the bottom two rows show examples of coincident Ly α and rest-frame UV continuum.</figcaption></figure>

imaging part of the VANDELS footprint, but in fact there is an even split (10 in HST coverage, 11 in ground-based coverage) for objects with absolute spatial offsets of >1.5 pixels. The choice of a 1.5 pixel threshold here versus the 2 pixel threshold above is simply to obtain a larger sample to mitigate small number statistics when comparing the number of objects in space-based and ground-based footprints.

While we cannot completely rule out the possibility that the emission line originates from a source too faint to detect in the continuum image, we stress that it is very unlikely for two main reasons: First, the sources would have to have very high equivalent widths in order for them to be detected in the spectra but not in the continuum images. For the 21 objects with offsets >1.5 pixels in magnitude, we estimated the minimum rest-frame Ly α equivalent width (*W*<sub>0,LYα</sub>) that would be required to produce the line flux observed in the spectrum. The limiting magnitudes in the four footprints (UDS-*HST*, UDS-GROUND, CDFS-*HST*, CDFS-GROUND) differ, so we use the limiting magnitude in the corresponding region for each object in that region to estimate the minimum *W*<sub>0,LYα</sub>. The median minimum *W*<sub>0,LYα</sub> is >156 Å in UDS-*HST*, >115 Å in UDS-GROUND, >71 Å in CDFS-*HST*, and >57 Å in CDFS-GROUND. Most of these equivalent widths, particularly in both UDS regions, are significantly larger than typical values of *W*<sub>0,LYα</sub> for star-forming galaxies at 3 < *z* < 6 (e.g. Shapley et al. 2003; Stark et al. 2011). The second reason we find it unlikely that the spatially offset emission lines come from undetected sources in

the continuum images is that the emission line would have to appear at exactly the right wavelength to mimic a spatially offset Ly α. Such bright serendipitous emission lines with no image counterpart are not detected at other wavelengths in the VANDELS spectra.

In Fig. 1, we show two examples of large spatial offsets. The top two rows show cases where the Ly α emission and UV continuum are likely coming from different regions of the same galaxy. The bottom two rows show cases where the Ly α emission is spatially coincident with the rest-frame UV continuum.

We show the redshift distribution, *W*<sub>0,LYα</sub> distribution, and rest-frame UV absolute magnitude (*M*<sub>UV</sub>) distribution of the sample in Fig. 3. All but one galaxy in our sample is at *z* < 5.5. so this is where we have the statistical power to constrain the Ly α spatial offset distribution. The equivalent widths were calculated using the Ly α flux and UV continuum from the 1D spectra provided in the VANDELS data products. The spatial apertures on the Ly α and the continuum are spatially centred on their respective peak emission, as to account for the potential Ly α spatial offsets as well as the centroid drift discussed above. These spatial apertures generally overlap. To obtain the Ly α flux, we fit a Gaussian to the 1D Ly α flux density and sum the Gaussian. A bandpass of 150 Å was used to calculate the continuum flux density from the 1D spectra.

We also investigated the distribution of redshift quality flags of our final sample. There were 0 (Q=0), 4 (Q=1), 4 (Q=2), 86 (Q=3), 203 (Q=4), and 8 (Q=9) spectra with the various quality flags. The vast majority (~95 per cent) of our sample has quality
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Distribution of physical 1D Ly α spatial offsets measured in this work. The top five panels show the distributions in distinct redshift bins, while the bottom panel shows the distribution for all objects. Because there is only one object in our sample at $z > 5.5$ (with offset = 0.3 kpc), we do not show that bin.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top: Redshift distribution of the galaxies showing Ly α emission used to constrain the intrinsic Ly α spatial offset distribution. The majority of objects have $z_{spec} < 4$, and there is only one object with $z_{spec} > 5.5$. Middle: The Ly α rest-frame equivalent width ($W_{0, Lyα}$) distribution of the sample. Bottom: The rest-frame UV absolute magnitude distribution of the sample.</figcaption></figure>

flags $Q \ge 3$, i.e. confidence of $\ge 95$ per cent in the $z_{spec}$ assigned to the galaxy by the visual inspectors. We also note that for the 8 ($\sim 3$ per cent) spectra with $Q = 9$ (i.e. a definite single emission line), the photometric redshift derived from the deep CANDELS imaging was also used to infer the redshift based on the single line. As a result, we expect that the contamination fraction in our final spectroscopic sample, i.e. the fraction of galaxies that are at a different $z_{spec}$ than the one listed in the catalogue, is very small.

The spatial offsets calculated from the 2D spectra are in units of pixels. However, we wish to constrain the physical Ly α spatial offset, so we convert the pixel offsets to proper kpc, using the $z_{spec}$ to do the angular diameter distance correction. While the redshift inferred from Ly α may be shifted by up to $\sim 500$ km s$^{-1}$ from systemic (e.g. Shapley et al. 2003; Song et al. 2014; Verhamme et al. 2018), the difference introduced by this shift into the angular diameter distance is insignificant. We show the distribution of physical 1D Ly α spatial offsets in proper kpc for the entire sample in Fig. 2. The one object in our sample at $z_{spec} > 5.5$ has $z_{spec} = 5.784$ and a spatial offset of $\sim 0.3$ pix $\simeq 0.06 \simeq 0.3$ kpc, which is not a significant offset given the uncertainty of this measurement (see Section 4).

## 4 LY α SPATIAL OFFSET MODELLING

Here we describe the Bayesian inference method we employ to constrain the intrinsic distribution of 2D Ly α spatial offsets from our 1D projected spatial offset measurements described in Section 3. We begin with Bayes’ theorem:

$$p(m|\{x_i\}) = \frac{p(\{x_i\}|m)p(m)}{E}, \quad (1)$$

where $p(m|\{x_i\})$ is the posterior for model parameter $m$ given the data set $\{x_i\}$, $p(\{x_i\}|m)$ is the total likelihood, $p(m)$ is the prior on the model parameter, and $E$ is the evidence.

We choose a simple model to represent the distribution of Ly α spatial offsets from rest-frame UV centre: a circular 2D Gaussian with $\mu_x = \mu_y = 0$. This choice is motivated by our ignorance of the shape of the offset distribution with the added bonus that it will allow us to write down the likelihood analytically. The choice to use $\mu_x = \mu_y = 0$ is motivated by the fact that there should be no preferred orientation to the offset distribution. Fig. 2 also provides evidence that the spatial offset distribution is centred at 0. The single parameter we want to constrain is the radial standard deviation of this symmetric Gaussian, $\sigma_{r, Lyα}$, which we define formally below. We opt not to use the morphology observed in the imaging data as a prior on the Ly α morphology because we do not want to bias the inference on the Ly α spatial distribution.

Our data consist of the spatial offset measurements made in the slits. These are 1D projected spatial offsets and therefore do not constitute true two-dimensional offsets. To account for this, we write the Gaussian distribution not in 2D but projected along one spatial dimension.

The 2D likelihood is

$$p(x, y|\sigma_x, \sigma_y) = \frac{1}{2\pi\sigma_x\sigma_y} \exp\left(-\left(\frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2}\right)\right), \quad (2)$$

where $\sigma_x$ and $\sigma_y$ are the standard deviations of the Gaussian in the $x$ and $y$ dimensions. For a Gaussian symmetric in $x$ and $y$, $\sigma_x^2 = \sigma_y^2$. And given that $\sigma_r^2 = \sigma_x^2 + \sigma_y^2 = 2\sigma_x^2$, we can write the above more
simply as

$$p(x, y|\sigma_r) = \frac{1}{\pi\sigma_r^2} \exp\left(-\frac{1}{\sigma_r^2}(x^2 + y^2)\right). \quad (3)$$

As we will see, $\sigma_r$ is related to the model parameter we want to constrain, $\sigma_{r, Ly\alpha}$. The two are not identical because there is measurement uncertainty that we must include in our likelihood function. This will not alter the overall form of the likelihood, so we continue with this expression to derive the 1D likelihood function.

When we take a spectrum of a source, we only obtain spatial information along the major axis of the slit, i.e. the axis perpendicular to the dispersion axis. This is illustrated in Fig. 4. The y-axis is the slit major axis, and therefore the y-component of the spatial offset is imprinted into the two-dimensional spectrum. As a result, only spatial offsets with a non-zero y-component (bottom panel of Fig. 4) will have a spatial offset in the slit.

We want to find the likelihood $p(y|\sigma_{r, Ly\alpha})$, where y is the spatial offset we measure in the 2D spectrum. We obtain this by integrating the 2D expression (equation 3) over x:

$$p(y|\sigma_r) = \int dx p(x, y|\sigma_r) = \frac{1}{\sqrt{\pi}\sigma_r} \exp\left(-\frac{y^2}{\sigma_r^2}\right). \quad (4)$$

This is simply the equation for a 1D Gaussian centred at y = 0 with standard deviation $\sigma = \sigma_r/\sqrt{2}$. We note that while the slit width could affect our 2D likelihood, it does not enter the final 1D likelihood $p(y|\sigma_r)$ because it is only a function of x. $p(x, y)$ is separable in x and y, so no matter what form the slit width enters into the 2D likelihood, it will always integrate to a constant when doing the x integral in equation (4).

If we were able to perfectly measure the spatial offsets in our slit, then we could just substitute $\sigma_{r, Ly\alpha}$ for $\sigma_r$ in equation (4) and then apply Bayes' theorem to obtain the posterior on $\sigma_{r, Ly\alpha}$. However, there is uncertainty in our measurement of y, the projected spatial offset in the slit. This uncertainty depends on a few factors, but most importantly the seeing and the S/N of the Ly $\alpha$ emission line. If we assume that this uncertainty manifests as Gaussian noise, we can include it in our likelihood via

$$\sigma_r^2 = \sigma_{r, Ly\alpha}^2 + \sigma_{r, msr}^2, \quad (5)$$

where $\sigma_{r, msr}$ is the measurement uncertainty on the radial offset. The radial offset uncertainty is related to the uncertainty in the offset in y, $\sigma_{y, msr}$, via $\sigma_{r, msr} = \sqrt{2}\sigma_{y, msr}$.

We estimate $\sigma_{y, msr}$ by performing simulations over a grid of FWHM of the seeing values (FWHM$_{seeing}$) and integrated Ly $\alpha$ S/N (S/N$_{int, Ly\alpha}$) values. We use a range of FWHM$_{seeing}$ from 0.1 to 1.0 arcsec in steps of 0.1 arcsec. For each step in FWHM$_{seeing}$, we simulate 1000 Ly $\alpha$ spectra with S/N$_{int, Ly\alpha}$ drawn from the uniform distribution in the range S/N$_{int, Ly\alpha}$ = 1–20, resulting in a total of 10000 simulations. The simulated spectra have the same dimensions as the real VANDELS spectra, and we place the simulated Ly $\alpha$ line at the same location in the spectrum in each simulation. For each simulation, we measure the centroid of the Ly $\alpha$ emission line using the same process we use to find the line centroid on the real data, as described in Section 3. We calculate the offset between the correct centroid and the recovered centroid for each simulation. To obtain the measurement uncertainty in each FWHM$_{seeing}$, S/N$_{int, Ly\alpha}$ cell, we take the standard deviation of these offsets in 2D bins given by $\Delta$FWHM$_{seeing} = 0.1$, $\Delta$S/N$_{int, Ly\alpha}$ = 2. The resulting grid of simulated measurement uncertainties is shown in Fig. 5. As expected, the measurement uncertainty decreases with S/N$_{int, Ly\alpha}$ and increases with FWHM$_{seeing}$.

Given that each spectrum has an arbitrary value of FWHM$_{seeing}$ and S/N$_{int, Ly\alpha}$, we want to be able to estimate the measurement uncertainty based on these parameters. To achieve this, we fit a function to the 2D histogram shown in Fig. 5:

$$\sigma_{y, msr}(m, \alpha, c) = m \left( \frac{\text{FWHM}}{\text{S/N}^{\alpha}} \right) + c, \quad (6)$$

finding best-fitting values of $m = 13.3$, $\alpha = 1.4$, $c = 0$. This is similar to the standard assumption of $\sigma_{y, msr} \sim \frac{\text{FWHM}}{\text{S/N}}$, but provides a better fit to the simulated data. $\sigma_{y, msr}$ depends much more strongly on S/N$_{int, Ly\alpha}$ than FWHM$_{seeing}$, especially at low S/N$_{int, Ly\alpha}$. Above S/N$_{int, Ly\alpha}$ = 10, $\sigma_{y, msr}$ flattens and higher S/N$_{int, Ly\alpha}$ do not yield significantly better measurement uncertainty. For the spectra with higher S/N than the range that we simulated, we extrapolate the function to estimate the measurement uncertainty for those objects. For an emission line with the median FWHM$_{seeing}$ (0.8 arcsec) and S/N$_{int, Ly\alpha}$ (10.9) of our sample, the measurement uncertainty is $\sim 0.4$ pixels $\simeq 0.08$ arcsec, corresponding to $\sim 0.5$ kpc at $z = 4$.

With this function in hand, we can rewrite the likelihood for a single spectrum in terms of our model parameter, $\sigma_{r, Ly\alpha}$ and the measurement uncertainty, $\sigma_{r, msr}$:

$$p(y|\sigma_{r, Ly\alpha}) = \frac{1}{\sqrt{\pi} \sqrt{\sigma_{r, Ly\alpha}^2 + \sigma_{r, msr}^2}} \exp\left(-\frac{y^2}{\sigma_{r, Ly\alpha}^2 + \sigma_{r, msr}^2}\right). \quad (7)$$

# 5 INFERENCE ON THE INTRINSIC LY $\alpha$ DISTRIBUTION

We wish to find the combined posterior for $\sigma_{r, Ly\alpha}$, $p(\sigma_{r, Ly\alpha} | \{y_i\})$, where $\{y_i\}$ are the set of all measured spatial offsets. The combined posterior is simply the product of the individual posteriors, $p(\sigma_{r, Ly\alpha} | y_i)$. We evaluate the combined posterior using the MCMC sampler from the python package EMCEE<fnref n="4" /> (Foreman-Mackey et al. 2013). The inputs to EMCEE are a likelihood, prior, and two parameters specific to the MCMC sampler. We used the likelihood in equation (7) and a flat prior on $\sigma_{r, Ly\alpha}$ over the interval 0–4 kpc. We originally explored a flat prior extending to larger values of $\sigma_{r, Ly\alpha}$, but the resulting posterior was zero-valued at larger values. The two MCMC parameters are the number of walkers and number of steps per walker. The number of walkers represents the number of independent paths through the parameter space that are taken by the sampler. We use 100 walkers and 250 steps per walker, chosen so that convergence is achieved. We discard the first 30 steps for each walker as these represent the burn-in steps when plotting our posterior or sampling from it to obtain derived quantities.

We show the final posterior using our entire data set (305 spectra) in Fig. 6. $\sigma_{r, Ly\alpha}$ is well constrained by our data, with a 68 per cent credible interval of $\sigma_{r, Ly\alpha} = 1.70^{+0.09}_{-0.08}$ kpc. We also compute the posterior after separating our data into five redshift bins: $3 \le z < 3.5$, $3.5 \le z < 4$, $4 \le z < 4.5$, $4.5 \le z < 5$, and $5 \le z < 5.5$. We show these five posteriors on Fig. 6. Interestingly, the posteriors suggest an evolution to smaller $\sigma_{r, Ly\alpha}$ with increasing redshift. In the bottom two panels of Fig. 6, we show this evolution more clearly. These values, along with the number of objects in each redshift bin, are listed in Table 1. The bottom left panel of Fig. 6 shows physical Ly $\alpha$ offset, whereas the bottom right panel shows apparent Ly $\alpha$ offset. $\sigma_{r, Ly\alpha}$ declines with redshift over the interval spanned by these data $3 < z < 4.5$ both in physical and apparent units. At $z >$

<footnote n="4">⁴http://dfm.io/emcee/current/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cartoon showing a spectroscopic slit for a simulated Gaussian source at $z = 4$ with FWHM$_{seeing} = 0.7$. **Top:** Left-hand panel shows the simulated Gaussian source perfectly centred in the slit. Right-hand panel shows a simulated 2D spectrum with the emission line arising from the simulated observation in the left-hand panel. The emission line is not spatially offset in the spectrum. **Bottom:** Same format as top panel, but now the same Gaussian source is offset by $r = 5$ kpc (~0.7) at $\theta = 45^\circ$ in the left-hand panel. The right-hand panel shows the resulting offset in the spectrum. In this configuration, the x-axis is the dispersion axis, so only the y-component of the spatial offset in the left-hand panel produces an offset in the spectrum. The x-component of the spatial offset results in slit loss, producing a fainter emission line in the bottom spectrum than in the top spectrum.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Simulation to estimate the measurement uncertainty on the spatial offset of Ly α in the 2D spectrum as a function of integrated Ly α S/N (S/N$_{int, Lyα}$) and FWHM$_{seeing}$. **Left:** Simulated measurement uncertainty on a grid of S/N$_{int, Lyα}$ and FWHM$_{seeing}$, **Middle:** Best-fitting model to the simulated measurement uncertainty, **Right:** Fractional residual between the simulated data and model.</figcaption></figure>

4.5 uncertainties are too large (due to small numbers) to conclude whether the trend continues.

To ensure that the decrease in $\sigma_{r, Ly\alpha}$ was not related to the fact that the higher redshift bins have fewer numbers, we reran the MCMC in redshift quartiles, i.e. four bins of increasing redshift, each containing an equal number of objects. We also find a decreasing trend in $\sigma_{r, Ly\alpha}$ in the four increasing quartiles: Q1: $2.36_{-0.21}^{+0.23}$, Q2: $1.63_{-0.15}^{+0.17}$, Q3: $1.39_{-0.12}^{+0.16}$, Q4: $1.35_{-0.12}^{+0.15}$.

# 6 DISCUSSION

## 6.1 Possible causes for Ly α spatial offsets

Our inference on $\sigma_{r, Ly\alpha}$ shows that over the redshift interval $3 < z < 5.5$ Ly α emission in galaxies can be spatially offset from the rest-frame UV continuum emission. We also found that the scale of the physical offsets decreases with redshift, at least up to $z = 4.5$. While a rigorous theoretical investigation into the origin of these
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Top: Posteriores on the standard deviation of the physical 2D Ly α offset distribution, σ<sub>r, Ly α</sub>, from the VANDELS spectra. We show the posterior for all (black) objects in the sample, as well as posteriors when the data are separated into multiple redshift bins. Bottom left: The median and 68 per cent credible intervals on σ<sub>r, Ly α</sub> (kpc) from the posteriors in each redshift bin in the top panel. σ<sub>r, Ly α</sub> shows a decreasing evolution with redshift. Bottom right: Same as bottom left but here we show σ<sub>r, Ly α</sub> in arcsec. For reference, the VIMOS slit width is 1 arcsec. In both bottom panels, the horizontal error bars represent the bin width of Δz = 0.5.</figcaption></figure>

<table><caption>**Table 1.** Constraints on σ<sub>r, Ly α</sub>.</caption><thead><tr><th>Redshift bin</th><th>σ<sub>r, Ly α</sub><br>(kpc)</th><th>N<sub>obj</sub></th></tr></thead><tbody><tr><td>3.0 ≤ z < 3.5</td><td>2.17<sup>+0.19</sup><sub>-0.14</sub></td><td>118</td></tr><tr><td>3.5 ≤ z < 4.0</td><td>1.46<sup>+0.14</sup><sub>-0.11</sub></td><td>99</td></tr><tr><td>4.0 ≤ z < 4.5</td><td>1.20<sup>+0.17</sup><sub>-0.13</sub></td><td>49</td></tr><tr><td>4.5 ≤ z < 5.0</td><td>1.38<sup>+0.29</sup><sub>-0.13</sub></td><td>35</td></tr><tr><td>5.0 ≤ z < 5.5</td><td>1.19<sup>+1.29</sup><sub>-0.33</sub></td><td>3</td></tr><tr><td>All</td><td>1.70<sup>+0.09</sup><sub>-0.08</sub></td><td>305</td></tr></tbody></table>

Note. N<sub>obj</sub> lists the number of objects in each redshift bin.

offsets is beyond the scope of this work, we explore some simple hypotheses using our data.

If the spatial offsets are mostly due to Ly α scattering, we might expect larger offsets in systems with more scattering, such as more luminous systems and systems with lower W<sub>0, Ly α</sub> (e.g. Shibuya et al. 2014). To test this hypothesis, we divided the sample into a ‘faint’ bin and ‘bright’ bin using the median absolute magnitude, M<sub>UV</sub> = -20.67 mag, to divide them so the two bins had the same number of objects. We then constrained σ<sub>r, Ly α</sub> in both bins using the MCMC approach described in Section 5. We found σ<sub>r, Ly α</sub> (faint) = 1.66<sup>+0.12</sup><sub>-0.10</sub> kpc and σ<sub>r, Ly α</sub> (bright) = 1.76<sup>+0.13</sup><sub>-0.12</sub> kpc. The two results are statistically consistent. This could simply be due to our sample size; a larger sample size would allow for a higher precision test which may indicate that either brighter or fainter galaxies have larger

σ<sub>r, Ly α</sub>. We also divided the sample into a low equivalent width bin and a high equivalent width bin (both with equal numbers of objects) and inferred σ<sub>r, Ly α</sub> in both bins. We found that σ<sub>r, Ly α</sub> was larger in the low equivalent width bin: σ<sub>r, Ly α</sub> (low W<sub>0, Ly α</sub>) = 1.92<sup>+0.14</sup><sub>-0.12</sub> kpc compared to the high equivalent width bin: σ<sub>r, Ly α</sub> (high W<sub>0, Ly α</sub>) = 1.51<sup>+0.11</sup><sub>-0.10</sub> kpc. The difference is statistically significant at ~3.5σ. This is evidence in favour of the scattering hypothesis.

A second potential cause for spatially offset Ly α is dust screening. In this scenario, dustier regions in the ISM could preferentially absorb Ly α photons over non-ionizing UV photons. If this were the case, then we might expect galaxies with more dust in general to exhibit a larger σ<sub>r, Ly α</sub>. Using the visual extinction magnitudes (A<sub>V</sub>) derived from SED fits to the galaxies in our sample (McLure et al. 2018), we find the opposite to be true. Splitting our sample into two equal-sized bins of A<sub>V</sub>, we find: σ<sub>r, Ly α</sub> (A<sub>V</sub> < A<sub>V, med</sub>) = 1.89<sup>+0.22</sup><sub>-0.16</sub> kpc for the less dusty half of the sample and σ<sub>r, Ly α</sub> (A<sub>V</sub> > A<sub>V, med</sub>) = 1.38<sup>+0.16</sup><sub>-0.13</sub> kpc for the dustier half, where A<sub>V, med</sub> = 0.25 mag. The difference is statistically significant at ~3σ. Because less dusty galaxies tend to have larger Ly α spatial offsets, the dust-screening hypothesis is unlikely, at least with the interpretation we put forward.

Finally, we consider a hypothesis where spatial offsets are driven by mergers, wherein one component preferentially emits the Ly α. In this scenario, we would expect the spatial offset to correlate with ellipticity in the UV image because mergers would have more elliptical isophotes. To estimate the UV ellipticity, we make use of the catalogue of structural parameters by van der Wel et al. (2012). This catalogue only overlaps with the 162/305 VANDELS objects
with space-based imaging, i.e. those in the UDS-HST and CDFS-HST footprints. The VANDELS team did not fit for the structural parameters of objects in the ground-based footprints due to the lower angular resolution in those regions.

We split the 162 spectra with derived ellipticities in half by the axis ratio, *q*, defined as the ratio of the semiminor axis to the semimajor axis, and infer $\sigma_{r,Ly\alpha}$ from both halves of the data. From the more elliptical half of the data we find: $\sigma_{r,Ly\alpha}$ ($q < q_{med}$) = 1.62$^{+0.14}_{-0.16}$ kpc, and from the less elliptical half of the data we find: $\sigma_{r,Ly\alpha}$ ($q < q_{med}$) = 1.83$^{+0.15}_{-0.18}$ kpc. The objects that are less elliptical in the UV tend to have larger spatial offsets, although the difference is not very significant ($<2\sigma$). As a result, our data do not show strong evidence that mergers play a significant role in driving the offsets, although we cannot rule out the possibility.

The results of these tests indicate that the origin of Ly α spatial offsets may be in part due to Ly α scattering. We point out that the scattering hypothesis and dust-screening hypothesis may not be independent because the more scattering events the Ly α photons experience the greater chance that the photons encounter dust. In reality, Ly α radiative transfer is more complicated than the simple scattering or dust-screening hypotheses we put forward due to the complex structure and kinematics of stars and gas in the ISM. Furthermore, if the Ly α offsets are caused by the accretion of less enriched satellites on to a more evolved host, interpreting the dust screening and merger hypotheses could be misleading. For example, if dust attenuation is measured in the host while Ly α arises from the satellite(s), then the dust screening hypothesis we put forward is not meaningful. Similarly, in the merger scenario it is possible that Ly α could arise from one merger component that is not detected in the UV. While it is unlikely this would occur over the entire population of mergers, if it did we would not expect spatial offset to correlate with ellipticity, invalidating the test we performed. Detailed observations of local starbust galaxies, such as those performed by the Lyman Alpha Reference Sample (LARS; Guaita et al. 2015) may yield further insight into the origin of Ly α spatial offsets.

## 6.2 Implications for slit-spectroscopic surveys

The fact that the Ly α spatial offsets we observe tend to decline in size with redshift is potentially interesting in the context of slit spectroscopy of Ly α. For example, several authors have found an evolving Ly α fraction as a function of redshift over the redshift range probed in this work (Stark et al. 2010, 2011; Hayes et al. 2011; Schenker et al. 2012; Cassata et al. 2015). We consider whether the evolution in the scale of the Ly α spatial offset distribution we infer could partially account for these trends simply due to differential slit losses.

In order to test this, we first must convert the physical offset scales, $\sigma_{r,Ly\alpha}$, into apparent offsets. We show the evolution of the apparent size of the Ly α offset scale in the bottom right panel of Fig. 6. Like the physical offsets, the apparent offsets inferred from our data exhibit evolution with redshift. The evolution is similar because the apparent size of objects at fixed physical size is relatively flat over the redshift range in our sample: 1 proper kpc = 0.13 ($z = 3$) = 0.14 ($z = 4$) = 0.16 ($z = 5$).

The magnitudes of the largest apparent radial offset scales in our sample are $\sigma_{r,Ly\alpha} \sim 0.2-0.3$ arcsec at $z \sim 3-4$. As illustrated in Fig. 4, the x-component of the radial offset, $r \cos \theta$, is what determines the slit loss. Here we are assuming the slit length is much larger than the seeing FWHM, which is a good approximation for our VIMOS observations, which use a minimum slit length of 7 arcsec.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Cumulative distribution functions (CDFs) for the slit loss due to spatial Ly α offsets. The CDFs are obtained by sampling Ly α radial offsets from Gaussian distributions with standard deviation determined from our inferred $\sigma_{r,Ly\alpha}$ values and uniform distribution in $\theta$. Solid (dashed) lines show the CDFs when using a slit width of 1 (0.5) arcsec and a seeing of 0.8, the median value from our observations. The VIMOS slit width is 1 arcsec. The slit loss is non-zero even for perfectly centred Ly α due to the seeing.</figcaption></figure>

We estimate the distribution of slit losses in each redshift bin by drawing radial offsets from a Gaussian distribution with a standard deviation given by our inferred $\sigma_{r,Ly\alpha}$ in that bin and drawing theta from a uniform distribution. We show the cumulative distribution function (CDF) of slit losses for the five redshift bins in Fig. 7. As expected, at the redshifts where we inferred the largest $\sigma_{r,Ly\alpha}$ ($z \sim 3-4$), the slit losses are largest. At $3 \le z < 3.5$, 50 per cent (20 per cent) of Ly α-emitting galaxies will be observed with slit losses $\gtrsim 0.17$ ($\gtrsim 0.24$). We note that there is a floor in the slit loss at $\sim 0.04$ for the VIMOS slit/seeing configuration due to the seeing blurring some flux outside of the slit even when Ly α is perfectly centred in the slit.

We also explored the slit losses for a slit width of 0.5 to bracket the slit widths used in ground-based spectroscopy. Because we model the seeing as a Gaussian, this simply shifts the CDFs to along the horizontal axis, i.e. to higher slit losses. We discuss the implications for the slit losses on future surveys with smaller slits in Section 6.

In calculating the slit loss, we assume that before convolution with the seeing, the Ly α emission is point-like. While Ly α has been shown to have significant spatial extent, often much larger than the UV continuum (Steidel et al. 2011; Wisotzki et al. 2016; Leclercq et al. 2017), what matters for the following analysis is the differential slit losses. In this work, we assume that the spatial extent of Ly α haloes is constant on the redshift range we probe. Future spatially resolved Ly α surveys with VLT/MUSE, for example, may be able to test this assumption at these redshifts.

To assess the impact of these slit losses on the evolution of the Ly α fraction, we consider an intrinsic rest-frame Ly α equivalent width distribution, $p(W_{0,Ly\alpha})$ before slit losses and then calculate the differential fraction of Ly α emitters we would measure if those emitters suffered the slit losses we found in each redshift bin. We use the $p_{z \sim 6}(W_{0,Ly\alpha})$ distribution at $z \sim 6$ compiled by De Barros et al. (2017) with the parametrization in terms of absolute magnitude ($M_{UV}$) by Mason et al. (2018a, see their equation 4). Their compilation is the largest sample with a well-defined selection function and homogenous observations available at $z \sim 6$, i.e. before Ly α is attenuated by the IGM neutral hydrogen due to reionization. In each redshift bin, we take the product $p_{z \sim 6}(W_{0,Ly\alpha}) \times (1 - \text{SL})$,
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Contribution of Ly α spatial offsets to the differential Ly α fraction (dx$_{Lyα}$/dz) over the interval $3 < z < 5.5$. The red shaded region shows the 68 per cent confidence interval on dx$_{Lyα}$/dz from Stark et al. (2011). The grey shaded region shows the 68 per cent confidence interval on dx$_{Lyα}$/dz from the slit losses induced by Ly α spatial offsets found in this work. The dotted lines show the best-fitting value of dx$_{Lyα}$/dz in each case. Both lines are parametrized to go through the same point in our lowest redshift bin, $3 < z < 3.5$, such that the evolution in the Ly α fraction can be easily compared. The two slopes are inconsistent at $2.4\sigma$, indicating that it is likely that Ly α spatial offsets are not entirely responsible for the increase in Ly α fraction observed at $3 < z < 6$. Given the large uncertainty in the Ly α fraction evolution, we cannot completely rule out the case where spatial offsets cause the evolution, however.</figcaption></figure>

where SL is the slit-loss distribution. We note that in each redshift bin we assume that the intrinsic EW distribution is $p_z \sim 6(W_{0,LY\alpha})$ and that only the slit losses are evolving with redshift. The result is an adjusted equivalent width distribution in each redshift bin which has suffered slit losses. The Ly α fraction is by definition the integrated probability of the distribution above $W_{0,LY\alpha} > X$, where X is a threshold value often chosen to be 25 or 50 Å for observational convenience. Stark et al. (2011) used $X = 25$ Å, so we adopt this threshold when calculating our Ly α fractions for the purposes of comparison. We generate errors on our Ly α fractions by resampling from our $\sigma_{r,LY\alpha}$ posteriors to generate resampled slit-loss CDFs. We then recalculate the Ly α fraction for each resampled CDF and the underlying distribution, $p(W_{0,LY\alpha})$, allowing us to construct Ly α fraction probability distributions in each redshift bin.

We show the differential Ly α fraction induced by slit losses that we obtain, $dx_{Lyα}/dz = 0.014 \pm 0.002$, in Fig. 8. This is an order of magnitude smaller than the evolution Stark et al. (2011) found from their sample of $3 < z < 6$ LBGs: $dx_{Lyα}/dz = 0.11 \pm 0.04$, although the difference is only $\sim 2.4\sigma$. We note that we are only interested in comparing the differential Ly α fraction, $dx_{Lyα}/dz$, so we scaled the Stark et al. (2011) Ly α fraction so that it equals the Ly α fraction measured in our lowest redshift bin. The normalization of the Ly α fraction is irrelevant in this comparison. Stark et al. (2010) computed $dx_{Lyα}/dz$ over their fainter luminosity range: $-20.25 < M_{UV} < -18.75$, but found that $dx_{Lyα}/dz$ was similar if they split their sample into two luminosity bins. We adopted $M_{UV} = -20.25$, the mean absolute magnitude from their entire sample, when evaluating the underlying distribution $p(W_{0,LYα}|M_{UV})$. We found that if we instead adopted $M_{UV} = -19.5$, the mean of their faint sample, we obtain a slope of: $dx_{Lyα}/dz = 0.016 \pm 0.004$, consistent with our fiducial result. We also note that the Stark et al. (2011) Ly α spectroscopy was performed with the same size slit widths (1.0 arcsec) as the

VIMOS slit widths, which are the slit widths we assumed in our slit-loss calculations.

The fact that the slit losses due to Ly α spatial offsets at $3 < z < 5.5$ are small is reassuring. Studies inferring the neutral hydrogen fraction during reionization from Ly α spectroscopy at $z > 6$ (e.g. Stark et al. 2011; Treu et al. 2012, 2013; Schenker et al. 2014; Tilvi et al. 2014; Mesinger et al. 2015; Mason et al. 2018a, 2019; Hoag et al. 2019) typically anchor to the Ly α equivalent width distribution at $z \lesssim 6$. If Ly α spatial offsets produced a large differential evolution in slit losses at these redshifts, and it was not accounted for, the inferred neutral fractions (and hence the reionization timeline) would be biased. What ultimately matters for these studies is whether slit losses at $z > 6$ are significantly different than what we have measured at $z < 6$. This could arise if the IGM or CGM is patchy on galaxy scales during reionization. UV-bright galaxies may clear channels in their CGM and the IGM through which Ly α can escape (e.g. Zitrin et al. 2015; Stark et al. 2017; Mason et al. 2018b), potentially resulting in an apparent spatial offset between Ly α emission and the non-ionizing continuum.

Finally, we note that in some cases, Ly α is bright enough to influence the rest-frame UV continuum images. In these cases, the spatial peak of the continuum image may be near the Ly α even if it is slightly offset. Depending on how common this is, it could mean that Ly α spatial offsets are already somewhat accounted for in the Ly α fraction studies. This is easily avoided by using a longer wavelength continuum image whose passband excludes Ly α during target selection and slit-mask design (e.g. Pentericci et al. 2018a; Hoag et al. 2019).

## 6.3 Implications for higher redshift Ly α surveys

The measurements in this work establish a baseline for Ly α spatial offsets at $z < 6$. The slit losses we have inferred will impact measured Ly α fractions at $z < 6$ and should be taken into account when Ly α fractions are reported from slit-spectroscopic observations. If spatial offsets measured at $z > 6$ are significantly different than our measurements at $z < 6$, it can be concluded that the difference is likely due to the neutral IGM during reionization. If the decreasing trend in Ly α spatial offsets we observed continues out to $z > 6$, then Ly α spatial offsets are not responsible for the decreasing Ly α fractions noted widely in the literature. In fact, the neutral hydrogen fractions inferred from current studies would need to be higher than reported in order to account for the fact that slit losses are larger at $z < 6$. Given the large uncertainty in $\sigma_{r,LY\alpha}$ in our highest redshift bin, $5 \le z < 5.5$, and the fact that the trend seems to be changing we cannot meaningfully constrain the impact that anchoring to the $z \sim 6$ rest-frame EW distribution would have on reionization studies at $z > 6$.

An assumption we made when calculating the differential Ly α fraction due to slit losses was that the spatial extent of Ly α haloes does not vary over $3 < z < 5.5$. While there is evidence that Ly α haloes are larger at $z \sim 3$ than in the local Universe (Wisotzki et al. 2016), it is not clear whether this trend continues to higher redshift. If it does, then slit losses will become more severe at higher redshift and will result in a steeper slope, $dx_{Lyα}/dz$, than we measured purely from spatial offsets alone.

While we do not constrain $\sigma_{r,LY\alpha}$ at $z > 5.5$, we can explore what spatial offsets might mean for future spectroscopic surveys at these redshifts. The Near InfraRed Spectrograph (NIRSpec) on JWST will have the capability to perform highly multiplexed multi-object spectroscopy of Ly α at $z > 7$ over a large $(3.6 \times 3.4 \text{ arcmin}^2)$ FOV. However, the effective slit size is $0.2$ (width) by $0.46$ (height),
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Projected slit losses for JWST/NIRSpec multi-object slit spectroscopy at $z \sim 7$ if Ly$\alpha$ spatial offsets endure to this redshift. The NIRSpec microshutters which form the slits have dimensions $0.′2$ wide by $0.′46$ tall, much smaller than the VIMOS slits used in this work. We explore two scenarios: *Top*: Expected slit-loss distribution assuming $\sigma_{r,\text{Ly}\alpha}$ is the same at $z=7$ as what we measured over our entire data set ($3 < z < 5.5$; Table 1). The solid black is for a single microshutter, while the dotted black is for multiple adjacent microshutters which effectively form a longer slit. *Bottom*: Expected slit-loss distribution if $\sigma_{r,\text{Ly}\alpha}$ is extrapolated from a power-law fit to our binned $\sigma_{r,\text{Ly}\alpha}$ measurements. If $\sigma_{r,\text{Ly}\alpha}$ does not decrease with redshift out to $z \sim 7$, then Ly$\alpha$ slit losses with JWST/NIRSpec will be catastrophically large. If $\sigma_{r,\text{Ly}\alpha}$ evolves according to our prescription, then slit losses will still be significant, but less severe.</figcaption></figure>

much smaller than typical slit spectrographs. If Ly$\alpha$ spatial offsets are not negligible at $z \gtrsim 7$, then slit losses may severely impact the detectability of Ly$\alpha$ with NIRSpec.

We forecast the slit-loss distribution at $z = 7$ from JWST/NIRSpec observations, using a similar method to the forecast shown in Fig. 7. We explore two scenarios: (i) no evolution in $\sigma_{r,\text{Ly}\alpha}$ versus (ii) evolution in $\sigma_{r,\text{Ly}\alpha}$ relative to $3 < z < 5.5$. For scenario (i), we set $\sigma_{r,\text{Ly}\alpha}$ at $z = 7$ equal to what we measured over our entire VIMOS data set ($3 < z < 5.5$), i.e. $\sigma_{r,\text{Ly}\alpha} = 1.70^{+0.09}_{-0.08}$ kpc. In scenario (ii) we use $\sigma_{r,\text{Ly}\alpha}$ at $z = 7$ from an extrapolation of a power-law fit to our constraints on $\sigma_{r,\text{Ly}\alpha}$ in the five redshift bins. Using this method we find $\sigma_{r,\text{Ly}\alpha} = 0.93^{+0.01}_{-0.15}$ at $z = 7$. In the absence of seeing, the angular size of the emission is governed by the size of the object and its distance. We used the size–luminosity relation compiled by Kawamata et al. (2015) to estimate the UV-continuum size at $z = 7$, and then assumed the Ly$\alpha$ emission is the same size. In reality, Ly$\alpha$ is generally more spatially extended than the UV (often significantly so up to $z \sim 3$; Wisotzki et al. 2016), so our projected slit losses are likely underestimated. For luminosities comparable to those studied in this work ($M_{\text{UV}} \sim -20.5$), this results in effective radii of $\sim 0.8 \pm 0.2$ proper kpc, or $\sim 0.1–0.2$ arcsec at $z = 7$. When simulating the slit-loss distributions we draw Ly$\alpha$ sizes from a Gaussian distribution with $\mu = 0.8$ kpc and $\sigma = 0.2$ kpc.

The slit-loss distributions for both scenarios are shown in Fig. 9. Slit losses will be significant in either case, but they are largest if $\sigma_{r,\text{Ly}\alpha}$ is comparable at $z = 7$ to what we measured over the entire

redshift range in this work. This is primarily due to the narrow width of the slits. With NIRSpec, one can open adjacent microshutters (assuming the shutter is not disabled), effectively increasing the slit height (cross-dispersion axis). We show how this would affect the slit-loss distributions in both scenarios in Fig. 9. As expected, it decreases the slit losses, yet the losses are still significant in either case. While there is a small gap between adjacent microshutters, we ignored these gaps when simulating the slit losses. As a result, the slit losses in this scenario are slightly underestimated. Our assumption that the Ly$\alpha$ size is comparable to the rest-frame UV size likely results in a much larger underestimate of the slit loss.

Measuring the spatial distribution of Ly$\alpha$ emission relative to the rest-frame UV light at $z > 6$ will be challenging, in part due to the sensitivity requirements but also because the Ly$\alpha$ fraction plummets at these redshifts, regardless of the mechanism. Surveys in lensed fields will show Ly$\alpha$ spatial offsets because spatial offsets in the source plane increase by a factor of $\mu$, the magnification factor, in the image plane. For the same reason, slit spectroscopic surveys in lensed fields will suffer from more severe slit losses.

It may already be possible to constrain $\sigma_{r,\text{Ly}\alpha}$ at $z \gtrsim 6$ with an instrument like Keck/DEIMOS. However, this would likely require a large dedicated effort. The first generation of instruments on upcoming 30-m class telescopes will very likely include wide-field optical spectrographs (i.e. TMT/WFOS, GMT/GMACS, ELT/MOSAIC), which will be capable of tackling this problem. VLT/MUSE also holds promise for constraining Ly$\alpha$ spatial offsets (e.g. Urrutiia et al. 2019), as long as sufficient astrometric precision can be achieved.

The natural multiplexing of space-based grism spectroscopy is also a potential avenue for constraining $\sigma_{r,\text{Ly}\alpha}$. While the HST WFC3/IR grisms have the wavelength coverage and spatial resolution to detect spatially resolved Ly$\alpha$ at $z > 5.5$, the sensitivity requirements are prohibitive (cf. Schmidt et al. 2016). With the grisms on JWST/NIRISS and JWST/NIRCAM, it may be possible to constrain $\sigma_{r,\text{Ly}\alpha}$ well into the reionization epoch ($z \gtrsim 7$).

## 7 SUMMARY

Using a large sample ($\sim 300$) of galaxies showing Ly$\alpha$ in emission from the VANDELS spectroscopic survey, we constrained the distribution of spatial offsets of Ly$\alpha$ emission relative to the rest-frame UV continuum. While we used slit spectroscopy which contains less spatial information than, e.g. IFU spectroscopy or narrow-band imaging to constrain the offsets, we employed a large sample which enabled us to make a statistically powerful measurement.

We parametrized the Ly$\alpha$ spatial offset distribution with a 2D circular Gaussian with zero mean and a single free parameter, $\sigma_{r,\text{Ly}\alpha}$, the standard deviation of this Gaussian expressed in polar coordinates. We constrained $\sigma_{r,\text{Ly}\alpha}$ using Bayesian inference by constructing a likelihood in terms of the measured spatial offset in the 2D VANDELS spectra. Using spectra over the entire redshift range ($3 < z < 5.5$) in our sample, we inferred a value of $\sigma_{r,\text{Ly}\alpha} = 1.70^{+0.09}_{-0.08}$ kpc. $\sigma_{r,\text{Ly}\alpha}$ declines from $2.17^{+0.19}_{-0.14}$ kpc at $3 \le z < 3.5$ to $1.19^{+1.29}_{-0.33}$ kpc at $5 \le z < 5.5$, or $\sim 0.′3$ to $\lesssim 0.′2$ in terms of apparent size.

We proposed three simple explanations for the origin of Ly$\alpha$ spatial offsets: scattering, dust screening, and mergers. The fact that $\sigma_{r,\text{Ly}\alpha}$ is higher in our lower rest-frame Ly$\alpha$ equivalent width bin supports the scattering explanation. We found that systems with lower dust content experienced significantly larger spatial offsets, contrary to the dust-screening explanation we put forward. There
was no conclusive correlation between ellipticity and Ly α spatial offset, though this may be due to the smaller sample size available for this test. We plan to investigate the origin of spatial offsets further in future work.

We examined the effect that the decreasing spatial offsets would have on slit losses for Ly α spectroscopic surveys. Slit losses alone could result in increasing Ly α fractions over the range 3 < z < 5.5: $dx_{Lyα}/dz = 0.014 \pm 0.002$. The effect is smaller than, but not entirely inconsistent (2.4σ) with the Ly α fraction evolution in the literature, from e.g. Stark et al. (2011): $dx_{Lyα}/dz = 0.11 \pm 0.04$. If Ly α spatial offsets continue to decline with redshift, they are not responsible for the decreased Ly α transmission measured at z > 6, typically attributed to reionization. Conversely, if spatial Ly α offsets become larger as the covering fraction of primeval galaxies increases, then they may represent a significant effect. Future Ly α surveys with JWST/NIRSpec may experience significant (> 50 per cent) slit losses if Ly α spatial offsets do not decline more rapidly out to z = 7 than the redshift evolution in our work suggests. The methodology we developed to infer $\sigma_{r, Lyα}$ at z < 6 can be readily applied to slit-spectroscopic surveys at z ≳ 6.

## ACKNOWLEDGEMENTS

Support for this work was provided by NASA through grant number JWST-ERS-01324 and HST-GO-15212 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. Support was also provided through NSF grant: COLLABORATIVE RESEARCH: The Final Frontier: Spectroscopic Probes of Galaxies at the Epoch of Reionization (AST 1815458, AST 1810822). The Cosmic Dawn Center is funded by the DNRF. MB acknowledges support by NSF through award NSF-AAG-1815458 and through NASA ADAP grant 80NSSC18K0945. The authors wish to thank Simon Birrer for helpful discussions during the development of the Bayesian inference analysis.

## REFERENCES

* Bacon R. et al., 2014, Messenger, 157, 13

* Becker R. H. et al., 2001, AJ, 122, 2850

* Behrens C., Braun H., 2014, A&A, 572, A74

* Bolton J. S., Haehnelt M. G., 2013, MNRAS, 429, 1695

* Brammer G. B. et al., 2012, ApJS, 200, 13

* Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

* Bunker A. J., Moustakas L. A., Davis M., 2000, ApJ, 531, 95

* Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storch-Bergmann T., 2000, ApJ, 533, 682

* Caruana J. et al., 2018, MNRAS, 473, 30

* Cassata P. et al., 2015, A&A, 573, A24

* Curtis-Lake E. et al., 2012, MNRAS, 422, 1425

* De Barros S. et al., 2017, A&A, 608, A123

* Dijkstra M., Kramer R., 2012, MNRAS, 424, 1672

* Dijkstra M., Mesinger A., Wyithe J. S. B., 2011, MNRAS, 414, 2139

* Dijkstra M., Wyithe S., Haiman Z., Mesinger A., Pentericci L., 2014, MNRAS, 440, 3309

* Fan X. et al., 2006, AJ, 132, 117

* Finkelstein S. L. et al., 2011, ApJ, 735, 5

* Finkelstein S. L. et al., 2012, ApJ, 756, 164

* Finkelstein S. L. et al., 2013, Nature, 502, 524

* Fontana A. et al., 2010, ApJ, 725, L205

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Fynbo J. U., Møller P., Thomsen B., 2001, A&A, 374, 443

* Galametz A. et al., 2013, ApJS, 206, 10

* Garilli B., Fumana M., Franzetti P., Paioro L., Scodeggio M., Le Fèvre O., Paltani S., Scaramella R., 2010, PASP, 122, 827

* Guaita L. et al., 2015, A&A, 576, A51

* Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633

* Guo Y. et al., 2013, ApJS, 207, 24

* Haiman Z., 2002, ApJ, 576, L1

* Haiman Z., Spaans M., 1999, in Holt S., Smith E., eds, AIP Conf. Proc. Vol. 470, After the Dark Ages: When Galaxies were Young (the Universe at 2 < z < 5). Am. Inst. Phys., New York, p. 63

* Hayashino T. et al., 2004, AJ, 128, 2073

* Hayes M., Scherer D., Östlin G., Mas-Hesse J. M., Atek H., Kunth D., 2011, ApJ, 730, 8

* Hayes M. et al., 2013, ApJ, 765, L27

* Hoag A. et al., 2019, ApJ, 878, 12

* Jones T. A., Ellis R. S., Schenker M. A., Stark D. P., 2013, ApJ, 779, 52

* Kawamata R., Ishigaki M., Shimasaku K., Oguri M., Ouchi M., 2015, ApJ, 804, 103

* Laursen P., Sommer-Larsen J., 2007, ApJ, 657, L69

* Le Fèvre O. et al., 2005, A&A, 439, 877

* Leclercq F. et al., 2017, A&A, 608, A8

* Malhotra S., Rhoads J. E., 2004, ApJ, 617, L5

* Mason C. A., Treu T., Dijkstra M., Mesinger A., Trenti M., Pentericci L., de Barros S., Vanzella E., 2018a, ApJ, 856, 2

* Mason C. A. et al., 2018b, ApJ, 857, L11

* Mason C. A. et al., 2019, MNRAS, 485, 3947

* McLure R. J. et al., 2018, MNRAS, 479, 25

* Mesinger A., Aykutalp A., Vanzella E., Pentericci L., Ferrara A., Dijkstra M., 2015, MNRAS, 446, 566

* Momcheva I. G. et al., 2016, ApJS, 225, 27

* Mortlock A., McLure R. J., Bowler R. A. A., McLeod D. J., Mármol-Queraltó E., Parsa S., Dunlop J. S., Bruce V. A., 2017, MNRAS, 465, 672

* Møller P., Warren S. J., 1998, MNRAS, 299, 661

* Nilsson K. K., Tapken C., Møller P., Freudling W., Fynbo J. P. U., Meisenheimer K., Laursen P., Östlin G., 2009, A&A, 498, 13

* Ono Y. et al., 2012, ApJ, 744, 83

* Ouchi M. et al., 2008, ApJS, 176, 301

* Pentericci L. et al., 2011, ApJ, 743, 132

* Pentericci L. et al., 2014, ApJ, 793, 113

* Pentericci L. et al., 2018a, A&A, 619, A147

* Pentericci L. et al., 2018b, A&A, 616, A174

* Rauch M. et al., 2008, ApJ, 681, 856

* Santini P. et al., 2015, ApJ, 801, 97

* Santos M. R., 2004, MNRAS, 349, 1137

* Schenker M. A., Stark D. P., Ellis R. S., Robertson B. E., Dunlop J. S., McLure R. J., Kneib J.-P., Richard J., 2012, ApJ, 744, 179

* Schenker M. A., Ellis R. S., Konidaris N. P., Stark D. P., 2014, ApJ, 795, 20

* Schmidt K. B. et al., 2016, ApJ, 818, 38

* Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., 2003, ApJ, 588, 65

* Shibuya T. et al., 2014, ApJ, 788, 74

* Song M. et al., 2014, ApJ, 791, 3

* Stark D. P., Ellis R. S., Chiu K., Ouchi M., Bunker A., 2010, MNRAS, 408, 1628

* Stark D. P., Ellis R. S., Ouchi M., 2011, ApJ, 728, L2

* Stark D. P. et al., 2017, MNRAS, 464, 469

* Steidel C. C., Giavalisco M., Pettini M., Dickinson M., Adelberger K. L., 1996, ApJ, 462, L17

* Steidel C. C., Bogosavljević M., Shapley A. E., Kollmeier J. A., Reddy N. A., Erb D. K., Pettini M., 2011, ApJ, 736, 160

* Swinbank A. M., Bower R. G., Smith G. P., Wilman R. J., Smail I., Ellis R. S., Morris S. L., Kneib J.-P., 2007, MNRAS, 376, 479

* Tilvi V. et al., 2014, ApJ, 794, 5

* Treu T., Trenti M., Stiavelli M., Auger M. W., Bradley L. D., 2012, ApJ, 747, 27

* Treu T., Schmidt K. B., Trenti M., Bradley L. D., Stiavelli M., 2013, ApJ, 775, L29

* Urrutia T. et al., 2019, A&A, 624, A141
van der Wel A. et al., 2012, ApJS, 203, 24

Vanzella E. et al., 2008, A&A, 478, 83

Verhamme A., Dubois Y., Blaizot J., Garel T., Bacon R., Devriendt J., Guiderdoni B., Slyz A., 2012, A&A, 546, A111

Verhamme A. et al., 2018, MNRAS, 478, L60

Wisotzki L. et al., 2016, A&A, 587, A98

Wisotzki L. et al., 2018, Nature, 562, 229

Zheng Z., Cen R., Weinberg D., Trac H., Miralda-Escudé J., 2011, ApJ, 739, 62

Zitrin A. et al., 2015, ApJ, 810, L12

¹Department of Physics and Astronomy, University of California, Los Angeles, CA 90095-1547, USA

²INAF, Osservatorio Astronomico di Roma, via Frascati 33, I-00078 Monteporzio Catone, Italy

³Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile

⁴Departamento de Física y Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena, Chile

⁵INAF, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti 93/3, I-40129 Bologna, Italy

⁶Department of Physics, University of California, Davis, 1 Shields Ave, Davis, CA 95616, USA

⁷Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK

⁸The Cosmic Dawn Center, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen Ø, Denmark

⁹INAF, Istituto di Astrofisica Spaziale e Fisica Cosmica Milano, via Bassini 15, I-20133 Milano, Italy

¹⁰Núcleo de Astronomia, Facultad de Ingeniería, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

¹¹Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹²Harvard-Smithsonian Center for Astrophysics, 60 Garden St, Cambridge, MA 02138, USA

¹³Geneva Observatory, University of Geneva, 51 ch. des Maillettes, CH-1290 Versoix, Switzerland

¹⁴Institut de Recherche en Astrophysique et Planétologie CNRS, IRAP, 14 Avenue E. Belin, F-31400 Toulouse, France

¹⁵Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany

¹⁶Department of Physics and Astronomy (DIFA), University of Bologna, Via Gobetti 93/2, I-40129 Bologna, Italy

¹⁷European Southern Observatory, Avenida Alonso de Cordova 3107, Vitacura, 19001 Casilla, Santiago de Chile, Chile

This paper has been typeset from a TeX/LA TeX file prepared by the author.