© ESO 2019

1. Introduction

The Lyα emission line plays a fundamental role in the investigation of several astrophysical problems, from detecting very distant sources to inferring physical parameters of galaxies. Understanding the mechanisms that regulate this strong emission is critical for interpreting the galaxy populations during the reionization epoch because the most luminous Lyα emitters (LAE) are frequent sources in the population of high-z galaxies (e.g., Stark et al. 2010, 2017; Pentericci et al. 2009). In addition, it is believed that the mechanisms that regulate the escape of Lyα photons are also responsible for the escape of LyC photons (Verhamme et al. 2017; Steidel et al. 2018; Marchi et al. 2018). Constraining these processes is therefore fundamental to investigate the ionizing radiation from high-redshift star-forming galaxies.

Because Lyα photons are mainly produced in HII regions by recombination processes, the line strength is related to the UV radiation that is emitted by young stars in galaxies, and so to the ongoing star formation. The escape of Lyα photons depends on a variety of physical processes, however. The interpretation from theoretical studies is that these photons are scattered by the neutral gas and absorbed by dust in the interstellar medium (ISM) (e.g., Verhamme et al. 2006; Gronke & Dijkstra 2016) and that these processes have a fundamental role in shaping the Lyα line profile, which is usually characterized by a single peak that is redshifted with respect to the systemic redshift and has an asymmetric shape. However, multiple-peaked profiles are observed in about 30% of the objects with Lyα emission at a redshift of about 2–3 in a sample of UV-selected star-forming galaxies (Kulas et al. 2012), and in most of these cases, the line is characterized by a red and a blue peak. The red peak frequently dominates the blue peak (e.g., Yamada et al. 2012; Verhamme et al. 2017). Double-peak profiles are now observed up to z >  6 (Matthee et al. 2018).

According to models, the shift of the red peak, the intensity of the blue peak, and the distance of the blue from the red peak critically depend on the neutral hydrogen column density (N(HI)) and the dust content of galaxies (e.g., Verhamme et al. 2015). The Lyα line shape is also affected by the temperature of the gas in the ISM and its clumping factor (e.g., Neufeld 1990, 1991).

The physical properties of Lyα emitting galaxies and the relation between the Lyα emission and galaxy physical properties have been extensively investigated by several groups in the past decade, using observations of large samples of galaxies at different redshifts. In particular, a bright Lyα emission appears to be associated with galaxies showing lower UV luminosities, lower metallicities (e.g., Shapley et al. 2003; Reddy et al. 2006), and lower dust content. It is also believed that lower stellar masses and lower star formation rates (SFRs) characterize the brightest LAEs, although a general consensus for these two quantities is still lacking (see, e.g., Pentericci et al. 2010; Kornei et al. 2010; Hathi et al. 2016; Du et al. 2018). While some authors advocate a scenario in which Lyα emitting galaxies represent early stages in galaxy formation (e.g., Finkelstein et al. 2007; Cowie et al. 2011), others find that these galaxies span a wide range of physical properties (e.g., Pentericci et al. 2009; Finkelstein et al. 2009, 2015; Kornei et al. 2010).

The recently completed spectroscopic survey VANDELS (a survey conducted with the Visible Multi Object Spectrograph, VIMOS, of the fields surveyed by the Cosmic Assembly Near-infrared Deep Extragalacic Legacy Survey, CANDELS, Ultra Deep Survey, UDS, and the Chandra Deep Field South, CDFS; Pentericci et al. 2018; McLure et al. 2018) has obtained spectra of star-forming galaxies spanning a wide range of redshifts and with unprecedented depth. With this new dataset we can analyze the Lyα line properties in extremely deep and moderate-resolution spectra for a large sample of z  ∼  3.5 galaxies, and infer the galaxy properties that correlate best with the intensity and kinematics of this line. In particular, we aim to assess the role of ISM kinematics and neutral hydrogen column density in the escape and distribution of Lyα photons.

The paper is organized as follows. In Sect. 2 we describe the criteria that we used to select the sample. In Sect. 3 we describe the methods we used to measure the spectral properties from the deep VANDELS spectra and the physical properties from SED fitting of the deep multiband imaging. In Sect. 4 we describe the relation of our sample to a more generic sample of galaxies with Lyα emission also from the VANDELS survey. In Sect. 5 we present the results of this analysis and show the correlations that we found between Lyα and galaxy properties. In Sect. 6 we finally discuss the main result, that is, the correlation between the Lyα peak shift and the velocity of the gas inside the ISM, and interpret this relation in the context of the shell model by assessing the role of the gas kinematics and HI column density in the escape of Lyα photons. Throughout the paper we adopt the Λ cold dark matter (Λ-CDM) cosmological model (H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3 and Ω_Λ = 0.7) and express all magnitudes in the AB system. Because the VANDELS spectra are calibrated in air and not in vacuum, we use air wavelengths to fit the position of the different lines. We also use positive values for the equivalent widths (EWs) of emission lines and negative for absorption lines. Our EWs are always rest-frame values.

2. Sample selection

To study the Lyα emission properties of high-redshift star-forming galaxies, we exploited the extremely rich and deep VANDELS dataset. In this section we first describe the VANDELS survey and then the criteria that we applied to select our sample.

2.1. VANDELS survey

VANDELS is an ESO public spectroscopic survey mainly targeting high-redshift galaxies at z >  2.5. It is specifically designed to be the deepest-ever spectroscopic survey of the high-redshift Universe. The survey is presented in two companion papers: McLure et al. (2018), who provide an overview of the survey strategy and target selection, and Pentericci et al. (2018), who focus on the observations and the first data release. VANDELS exploited the VIMOS multi-object spectrograph to obtain ultra-deep optical spectra of about 2100 high-redshift galaxies within the redshift interval 1 ≲ z ≲ 7, covering a total area of 0.2 deg² centered on the CANDELS-UDS and CANDELS-CDFS fields (Grogin et al. 2011; Koekemoer et al. 2011). Through the ultra-long exposure-time strategy, each target has been observed from at least 20 h to 80 h at most, using the red medium-resolution grism of VIMOS, which covers a wavelength range of 4800 Å ≤ λ ≤ 10 000 Å with an average resolution of 580. To select VANDELS sources, we adopted in the CANDELS/Hubble Space Telescope (HST) areas (CDFS-HST and UDS-HST) the photometric catalogs based on H₁₆₀-band detections that are provided by the CANDELS team (Galametz et al. 2013; Guo et al. 2013). Specifically, the CDFS-HST catalog is based on 17 photometric broadband filters (Guo et al. 2013), while the UDS-HST catalog includes photometry in 19 broadband filters (Galametz et al. 2013). However, because the VIMOS footprint is wider than that of the HST, many VANDELS sources fall outside the CANDELS footprints but are still covered by ground-based imaging and imaging obtained with the Infrared Array Camera (IRAC) on board the Spitzer Space Telescope. The two photometric catalogs we used for the sources in the wide-field regions were generated by McLure et al. (2018) using the publicly available imaging in 12 filters for UDS-GROUND and 17 filters for CDFS-GROUND, spanning the range from U to K band in both cases. These photometric catalogs will not only be combined with the spectroscopic data to create this unique dataset to study the high-redshift galaxy properties and evolution, but have also allowed the VANDELS team to carefully preselect the VANDELS targets on the basis of their accurate photometric redshift measures (McLure et al. 2018). About 85% of the galaxies were selected to be at z ≥ 3, the majority of which are Lyman break galaxies (LBG), which are star-forming galaxies selected on the basis of their colors redward and blueward of the observed-frame Lyman limit. The masks were prepared using the VIMOS mask preparation software (Bottini et al. 2005) that is distributed by ESO. The VANDELS spectroscopic data were reduced using the fully automated pipeline easy-life (Garilli et al. 2012), an updated version of the algorithms and dataflow of the original VIPGI system (Scodeggio et al. 2005), which generated wavelength- and flux-calibrated spectra directly from the raw data. Finally, the spectroscopic redshifts were estimated using the Pandora software package, within the EZ environment (Garilli et al. 2010). The same procedure as was adopted in VUDS was applied (see Le Fèvre et al. 2015), that is, each spectrum was analyzed by two different team members. In addition, a final independent check by the two principal investigators was carried out. We currently use the same quality flags as for VUDS (see Le Fèvre et al. 2015) to asses the reliability of the spectroscopic redshifts.

The second VANDELS data release became public at the beginning of October 2018 and contains the spectra and redshifts of 1339 objects (557 in CDFS and 762 in UDS) with more than 200 spectra with ultra-deep 80-h exposures¹.

2.2. Sample

We selected all the objects in the VANDELS data release 2 that showed Lyα in emission and had a secure spectroscopic redshift with VANDELS reliability flags 3 or 4 (corresponding to a probability greater than 95% for the spectroscopic redshift to be correct, see Pentericci et al. 2018). This sample, which we will call the parent sample, contains 305 star-forming galaxies in the redshift range 3 ≲ z ≲ 4.5. One hundred and forty-six of these galaxies are in the CDFS field and 159 are in UDS. One hundred and fifty-four galaxies (78 in the CDFS and 76 in UDS) fall within the CANDELS regions (Grogin et al. 2011; Koekemoer et al. 2011) and are therefore covered by very deep multiwavelength HST imaging, while the remaining 151 (68 and 83 in the CDFS and UDS, respectively) are in the wide-field areas of the VANDELS fields and mainly benefit from ground-based imaging.

We wish to study the kinematics of the Lyα emission, which gives us important information about the ISM and the HI column density as shown by Verhamme et al. (2015). To do this, we need a precise measure of the velocity of the Lyα line with respect to the systemic redshift of the galaxy. We therefore restricted our sample to the galaxies whose systemic redshift we were able to estimate in a reliable way. The systemic redshift (z_sys) is ideally defined by the positions of photospheric absorption lines that are generated by the absorption of the stellar radiation by the most external layers of the star. The most important photospheric absorption lines in the UV regime for highly star-forming galaxies would be OIV λ1343 Å, SiII λ1417 Å, and SV λ1500 Å (Talia et al. 2012). However, these lines are either too weak or even not detected in our ultra-deep spectra. As an alternative, we can use nonresonant nebular emission lines that originate from the photoionization of nebular regions by the radiation from very massive O and B stars and can therefore be used to trace the systemic redshift. The only nebular emission line that enters the wavelength range covered by the VANDELS spectra at redshift ∼4 is the CIII] λ1909 Å emission. This is a semi-forbidden line at λ = 1907.7 Å and it is known to be a doublet at a rest-frame vacuum wavelength of 1906.68 Å and 1908.68 Å (because the VANDELS spectra are air-calibrated, we here use the corresponding air wavelengths 1906.05 Å and 1908.05 Å). Another emission line often appears in our spectra, however, that can trace the systemic redshift. This is the HeII line at rest-frame vacuum wavelength 1640 Å (1639.83 Å in air). This feature can be produced either by Wolf-Rayet stars (WR), which are thought to represent a stage in the evolution of O-type stars with masses > 20 M_⊙, or by young star clusters with intense star formation (e.g., Talia et al. 2012; Cassata et al. 2013). Recently, Schaerer et al. (2019) also proposed that high-mass X-ray binaries could be the main source of nebular HeII emission in low-metallicity star-forming galaxies.

The line can show a broad profile because of the strong winds powered by the intense star formation, or can appear as a nebular narrow line as the strong UV emission from the stars photoionizes the surrounding medium. Because of the possible effects of the winds, this line is not always a good tracer of the systemic redshift. However, the results presented in this paper do not change when we exclude the sources whose systemic redshift was derived using HeII λ1640 Å, suggesting that this line provides a good estimate of the systemic redshift at least for our galaxies.

We therefore selected all galaxies from the parent sample that showed a CIII] λ1909 Å and/or HeII λ1640 Å emission lines to derive the systemic redshift. We only included objects where the CIII] λ1909 Å (or HeII λ1640 Å) were reasonably clear from skyline residuals and had a signal-to-noise ratio (S/N) ≥ 3 to allow a precise measurement of the line peak (see below). With these criteria we selected 60 star-forming galaxies in a redshift range of 3 ≲ z ≲ 4.5.

From this sample we further excluded six objects whose Lyα was too noisy to give a precise estimate of the peak shift and the EW, mainly because they were affected by skyline residuals, which prevented us from carrying out the measures in a reliable way. We also excluded two objects with intense CIV λ1549 Å emission that could indicate the presence of an AGN.

The final selected sample contains 52 sources. Twenty-nine (25) of these are in the CDFS (UDS) field. Thirty sources (18 in the CDFS and 12 in UDS) are in the CANDELS regions (Grogin et al. 2011; Koekemoer et al. 2011) and are therefore covered by very deep multiwavelength HST imaging, while the remaining 24 (11 and 13 in the CDFS and UDS, respectively) are in the wide-field areas of the VANDELS fields. Throughout the paper, we refer to this group of galaxies as the selected sample.

The distributions in redshift and z-band magnitude (which probes the rest-frame UV continuum at ∼1800 Å for all sources) of the parent and the selected sample are shown in Fig. 1 in yellow and green, respectively. The median redshift and the median z-band magnitude are z_med = 3.47 and zmag_med = 25.4 for the selected sample and z_med = 3.47 and zmag_med = 25.5 for the parent sample. The UV continuum luminosity distributions of the two samples are very similar. However, because the selected sample mostly contains galaxies that also show CIII] λ 1909 Å emission, we expect this sample to be in part biased to brigher Lyα luminosities because these two emissions appear to be correlated according to several studies (e.g., Shapley et al. 2003; Stark et al. 2014; Maseda et al. 2017).

Fig. 1. Redshift distribution (left panel) and z magnitude distributions (right panel) of the sources in the selected sample (green empty histogram) and in the parent sample (yellow histogram).

Finally, we remark that although our samples were selected for their Lyα emission, only a fraction might strictly be considered LAEs if the traditional threshold of EW(Lyα)> 25 Å were applied. In our parent sample, 26% of all galaxies have EW(Lyα) ≥ 25 Å, in agreement with the statistics at this redshift (Shapley et al. 2003; Stark et al. 2010).

3. Method

For each galaxy in our samples, the VANDELS dataset provides not only a deep spectrum but also deep optical, near-IR, and Spitzer photometry. To best exploit this unique legacy dataset, we therefore used the optical spectra to characterize the Lyα emission, both spectrally and spatially, the interstellar absorption lines and the CIII] λ1909 Å properties, while we used the VANDELS photometric catalogs (McLure et al. 2018) to derive the physical parameters. We note that for the parent sample we only measured the Lyα EW from the one-dimensional spectra and the physical properties from the available photometry.

3.1. Spectroscopic measures

From the VANDELS spectra we measured the properties listed below.

The systemic redshift (z_sys). We evaluated the systemic redshift from the CIII] λ1909 Å line or from the HeII λ1640 Å line. In the majority of galaxies with CIII] λ1909 Å emission, the two carbon lines appear blended; for these and for the objects with HeII λ1640 Å emission, we derived the systemic redshift from the observed-frame spectra by fitting a single Gaussian profile assuming that the continuum has a well-defined slope (linear fit). In the minimization procedure we took the VANDELS noise spectrum into account, which contains the observational and instrumental uncertainties. While fitting the continuum, we excluded the wavelength regions with absorption lines or where strong residuals from sky lines were present because these regions might alter the slope of the fit. This procedure was done by visual inspection in each spectrum. We finally estimated the systemic redshift from the best-fit mean of the Gaussian profile. We show an example of a fit of the CIII] λ1909 Å in the bottom right panel of Fig. 2. The blue line is the VANDELS spectrum, the magenta range is the part of the spectrum that we fit, and the dark blue line is the best fit. For the few objects that clearly show the two peaks of the CIII] λ1909 Å line due to the better S/N, we used the same procedure, but included two Gaussian profiles, and we evaluated the systemic redshift directly from the fit by imposing that the position of the two mean peaks lay at 1906.68 (1 + z) Å and 1908.68 (1 + z) Å. We did not apply any restriction to the CIII] doublet ratio during the fitting procedure, but the ratios that we obtain have values in the range ∼1−1.5. These are typical values for the electronic densities of local star-forming regions (see, e.g., Osterbrock & Ferland 2006; Quider et al. 2009). To derive the errors on these measures, we used the equation given by Lenz & Ayres (1992) for a Gaussian fit.

Fig. 2. Different regions of the spectrum of object CDFS122764 at z = 3.49 are shown for reference. Upper panel: UV continuum between about 1250–1750 Å rest frame, where the most important absorption features are clearly visible, lower left panel: Lyα line profile, and lower right panel: CIII] λ1909 Å emission. We used the magenta parts of the spectrum to perform the fits. The dark blue lines in the upper panel and in the lower right panel are the best fits, and the shaded region in the lower left panel represents the normalized integrated flux of the Lyα that we used to evaluate the EW of the line. The vertical green lines show the positions of the lines at rest-frame air wavelength. The blue vertical dashed line in the lower left panel indicates the position of the Lyα peak. The cyan line is the VANDELS noise spectrum.

The Lyα redshift and Lyα velocity shift (Lyα_shift). The Lyα redshift was evaluated as the peak of the Lyα line, that is, at the position where the flux has a maximum. In contrast to the CIII] emission, we did not attempt a Gaussian fit of the Lyα line because its shape is known to be asymmetric (e.g., Shimasaku et al. 2006). For the few spectra (5 out of 52) where we were able to clearly identify a double-peak profile of the Lyα line, we took the peak position of the red peak, which is also the highest peak. We then evaluated the Lyα velocity shift as the velocity difference between the Lyα redshift and the systemic redshift. The error was evaluated from the usual propagation formula using the error on z_sys and the error for the Lyα wavelength position was evaluated in a similar way.

In some spectra the wavelength range around the Lyα was affected by noise, for example, caused by residuals from skylines, and therefore it was not possible to derive an accurate peak position of the Lyα line. We decided to keep these spectra to derive the other spectral parameters, such as EW(Lyα) (see below). In the end we were able to measure the Lyα velocity shift for only 46 of 52 galaxies. We show the Lyα line profile of the galaxy CDFS122764 in the bottom left panel of Fig. 2 for reference. Here, the Lyα velocity shift is the velocity difference between the peak of the line (dotted line) and the green line that represents the position of the line at systemic redshift.

The Lyα FWHM. For the same 46 objects for which we measured the Lyα peak position, we also evaluated the full width at half-maximum (FWHM) of the line, directly taking the wavelength interval at half of the flux of the peak of the line. Again we preferred a direct measure of the FWHM rather than a Gaussian fit because the line is clearly asymmetric in many cases and the Gaussian does not represent its shape well. We estimated the error using the formula in Lenz & Ayres (1992).

The Lyα and CIII] λ1909 Å EWs. To estimate the Lyα EW, we first normalized the line flux. We used the mean flux in a spectral range between the Lyα and the SiII line at 1260 Å as continuum. The SiII line is the first strong absorption line redward of the Lyα line that we can clearly see in the VANDELS spectra. We point out that we did not include possible effects from HI absorption around the Lyα line in our continuum estimate. This might result in a slightly underestimated EW. We then evaluated the EW(Lyα) from the integrated normalized flux. To estimate the error on the EW(Lyα), we used Eq. (13) in Ebbets (1995). We show the Lyα line profile of the galaxy CDFS122764 in the bottom left panel of Fig. 2 for reference. The shaded region represents the integrated flux of the line.

The EW(CIII]) was evaluated directly from the fit we used to derived the systemic redshift by normalizing the line flux with the best-fit line of the continuum. For the EW(Lyα) we used the Ebbets (1995) formula to derive the errors. As continuum uncertainty we used the standard deviation in a wavelength range of 100 Å observed-frame redward of the CIII] λ1909 Å line.

The LIS redshift (z_IS) and velocity shift (IS_shift). The low-ionization absorption lines (LIS) are among the strongest lines that can be detected in UV rest-frame spectra of star-forming galaxies. In particular, the strongest lines that we can observe in the majority of our spectra are the SiII λ1260 Å, CII λ1334 Å, SiII λ1526 Å, FeII λ1608 Å, and the AlIII λ1679 Å lines. In high-redshift star-forming galaxies whose UV continuum is dominated by young stars, these lines are mainly produced by the absorption of the stellar radiation by the ISM and are a fundamental tool to study its properties.

We used the velocity shift of the LIS to probe the velocity of the ISM. We chose to use only the lines with the highest S/N in our spectra: SiII λ1260 Å CII λ1334 Å, and SiII λ1526 Å. We evaluated the redshift of these lines using a combined fit of all the three features, even though in some cases we excluded the lines that did not have a sufficient S/N (lower than 3). We fit the continuum with a fifth-order polynomial function plus three Gaussian profiles to fit the lines (or fewer when we were unable to fit all the lines) to derive z_IS, imposing the line centers at 1259.33 (1 + z) Å, 1333.95 (1 + z) Å, and 1526.13 (1 + z) Å, respectively. We note that even if in some cases the interstellar (IS) absorption lines can be asymmetric, at our moderate resolution, a Gaussian profile is a good approximation of the line shape. We again excluded the wavelength regions where other absorption lines were visible or where strong residuals from sky lines were present.

We then evaluated the LIS velocity shifts as the velocity difference between the LIS redshift and the systemic redshift. The error was evaluated from the error propagation formula using the error on z_sys and the error for the LIS redshift.

We were able to estimate the IS_shift for only 26 galaxies because the S/N of the LIS for the remaining sources was not sufficient to evaluate the IS velocity. We show the spectral continuum of the galaxy CDFS122764 in the top panel of Fig. 2 for reference. The blue line is the VANDELS spectrum and the magenta line is the part of the spectrum that we used for fitting the LIS and the continuum. For this object we were able to fit all the three LIS, and we excluded the strong absorption lines from the fit (for instance OIdb λ1303 Å, the SiIV+OIV λ1397 Å and the CIV λ1549 Å) that would have altered the shape of the continuum.

Spatial extensions in Lyα and UV. The spatial extent in Lyα spatial extent was directly evaluated from the two-dimensional spectrum of the sources. To estimate the background, we selected four spatial windows on the 2D spectrum, as shown in Fig. 3, which do not contain the continuum or the Lyα residuals that can affect the regions below and above the Lyα emission in the 2D spectrum. We then separately collapsed the two bottom and upper windows and the region that contains the Lyα line and evaluated the mean value for each y-pixel. We then fit it with a Gaussian profile using a straight line to fit the background and finally evaluated the Lyα_ext as the FWHM of the best fit. To determine the uncertainties in these measurements, we evaluated not only the error on the Gaussian fit as in Lenz & Ayres (1992), but also included the error on the background, which we evaluated through Monte Carlo simulations by varying the value of each y-pixel within the standard deviation of the background and fitting each time with a Gaussian profile. The final uncertainties were computed by combining the two errors in quadrature, assuming that they are independent. We note that the two different errors are comparable on average. We evaluated the UV spatial extension in the same way as we measured the Lyα extent, but we used the spectral region between the two spatial windows redward of the Lyα line for continuum, as shown with a red rectangle in Fig. 3. The typical spatial windows are ∼25 Å observed-frame × 30 pixels for the Lyα region and ∼150 Å observed-frame × 10 pixels for the UV continuum. We further checked these regions objects by objects visually inspecting the 2-dimensional spectra. We also made sure that the Lyα emission does not contaminate the UV measurements by allowing in most cases for a gap between the two spatial windows.

Fig. 3. Two-dimensional spectrum of object CDFS122764 that we show to visualize the method that we used to evaluate the Lyα and the UV extents. The large green rectangles show the spatial regions of the two-dimensional spectrum that we used to estimate the background, while the yellow (red) square around the Lyα (UV continuum) is the region where we collapsed the spectrum.

Our measurements of the Lyα and UV sizes cannot be taken as absolute estimates of these quantities because we did not deconvolve them by the point spread function (PSF) of the observations. Deriving the precise PSF of each spectrum is particularly difficult because observations were spread over a long time. In most but not all masks, point-source reference objects were inserted, but almost all of our targets were observed in two or four different masks. However, this is probably not a crucial issue for our analysis because we are only interested in relative values. First, the Lyα and UV sizes are evaluated in a very close wavelength range, so that we can safely assume that the PSF affects both measurements in the same way. In addition, all observations were taken under very similar average conditions over many separate nights (see Pentericci et al. 2018). We monitored the PSF of the reference objects in each 20-h mask, and they are spread over a very small interval with an average value of 0.8 arcsec and a standard deviation of 0.07. This implies that the final PSF of all sources should be very similar. We also particularly monitored the objects whose size determination was equal to or slightly below the average PSF and conservatively assumed that their size was equal to this average PSF.

3.2. Physical parameters

3.2.1. M^*, SFR, E(B–V), and mass-weighted galaxy age

We performed SED fitting to derive the physical parameters of the sources in our sample. For the objects that fall within the CANDELS CDFS and UDS footprints, we used the CANDELS official photometric catalogs (Galametz et al. 2013; Guo et al. 2013), while for the galaxies in the wide-field areas, we used the photometric catalogs described in McLure et al. (2018). Both catalogs span a wavelength range that extends from the U band to the K band observed frame.

To perform the SED fitting we used BEAGLE (which stands for Bayesian analysis of galaxy SEDs, Chevallard & Charlot 2016), a new-generation tool that allows modeling spectroscopic and photometric galaxy observations with a self-consistent physical model. We used BEAGLE with the most recent prescriptions of Charlot & Bruzual (in prep.) to describe the emission from stars and the models from Gutkin et al. (2016) to account for emission from photoionized gas. The redshift was fixed at the systemic value that we evaluated from the VANDELS spectra for the objects in the selected sample, and we used the official VANDELS spectroscopic redshift for the galaxies in the parent sample. We used the Chabrier (2003) initial mass function and a Calzetti et al. (2000) attenuation law because Cullen et al. (2018) recently showed that VANDELS galaxies at 3 <  z <  4 are consistent with this law. We used a constant star formation history (SFH) and fit the metallicity using a Gaussian prior centered at 0.14 Z_⊙ ² with a σ of 0.07 Z_⊙. The value of 0.14 Z_⊙ corresponds to the average metallicity of our sample evaluated applying the method described in Sommariva et al. (2012) on the stacked spectrum of all the galaxies in the selected sample. Briefly, this method exploits spectral indices that are present in the UV spectra of star-forming galaxies in the range between ∼1300 Å and 1600 Å rest frame, whose EWs only depend on the metallicity of the objects. Of the indicators described in Sommariva et al. (2012), we used the spectral indices F1370, F1425, F1460, and F1501 because they are not affected by the interstellar absorption lines at our resolution. The value of the standard deviation of the Gaussian prior is the dispersion of the values derived by the four indices. We remark that the aim of this measurement is only to determine an average metallicity to lower the number of free parameters in the SED fitting.

We derived the stellar masses, SFR, mass-weighted ages, and the color excess E(B − V), which measures the dust content. We show in Fig. 4 an example of an SED fit of a galaxy in our sample. In the top panel we show the probability distribution functions of the derived physical parameters. We also show the best-fit values obtained by allowing the metallicity to vary in a uniform range [0.006 Z_⊙; 1.74 Z_⊙]. For most of the galaxies the results are fully consistent within 1σ, with the age being the most scattered parameter. Finally, in the bottom panel of Fig. 4, we compare the observed photometry with the photometry predicted by BEAGLE.

Fig. 4. Top panel: posterior probability distributions of the physical parameters derived for object CDFS020583 at z = 3.5. From left to the right we show the logarithmic mass in units of solar masses, the logarithmic metallicity in units of solar metallicity, the dust attenuation coefficient, the age, and the logarithmic specific SFR. The gray bands are the 1σ confidence intervals. The red crosses are the 1σ confidence intervals obtained without the prior on the metallicity. The different blue contours indicate the 1σ, 2σ, and 3σ confidence levels. Bottom panel: spectral energy distribution of the same galaxy. The blue points are the photometric observations, and the orange shaded regions are the 68% confidence intervals derived by BEAGLE. The black solid lines are ten of the theoretical SEDs with the highest likelihood values.

We show in Fig. 5 the SFR as a function of the stellar mass for the sources in our selected sample in green and the parent sample in yellow, along with the star-formation main sequence relation evaluated by Santini et al. (2017) for star-forming galaxies in the first four HST Frontier Fields in the redshift range 3 ≲ z ≲ 4 (blue line). We also show in the same figure the distributions of SFRs and stellar masses for the sources in our sample (green histogram) and in the parent sample (yellow histogram). Both the samples lie on average on the star formation main sequence, as do the majority of high-redshift star-forming galaxies. This relation is believed to characterize galaxies that have grown on long timescales as a consequence of smooth cold gas accretion from the intergalactic medium (IGM) (e.g., Dekel et al. 2009). We note that there is some scatter around the main sequence in Fig. 5, which is consistent with that of the Santini et al. (2017) sample, however. The reason for this scatter could be related to the SFH we used in the SED fitting procedure (e.g., Cassarà et al. 2016). To determine this, we also performed the same SED fitting using free metallicity in the range [0.006 Z_⊙; 1.74 Z_⊙] and two different types of SFH: an exponentially declining SFH and a delayed exponentially declining SFH. In both cases the scatter around the main sequence is consistent with that obtained with a constant SFH. Moreover, the individual values are also consistent. In particular, while the stellar masses are usually not affected by the choice of SFHs, the SFRs can instead vary significantly. In this case, however, when we compare the results obtained with constant SFH and exponential and delayed SFHs, they are still consistent within the 1σ errors. In a future work we will use BEAGLE to assess the implications of this choice with a larger sample of galaxies in the VANDELS survey.

Fig. 5. Star formation rate as a function of the stellar mass for the objects in our sample (green points) and the parent sample (yellow points). The blue line is one of the latest star formation main sequence relation evaluated by Santini et al. (2017) for galaxies in the first four HST Frontier Fields in the redshift range 3 ≲ z ≲ 4.

3.2.2. β slope

We also determined the slope β of the UV continuum. We adopted the common power-law approximation for the UV spectral range F_λ ∝ λ^β (see, e.g., Castellano et al. 2012) and estimated the slope β by fitting a linear relation through the observed AB magnitudes: M_i = −2.5(β + 2.0)logλ_i + c, where M_i is the magnitude in the ith filter at effective wavelength λ_i. The filters we used span the rest-frame wavelength range λ ∼ 1250−2500 Å, such that the Lyα emission line is not included in any of the filters. We divided our sources into two different redshift bins: 3 ≲ z ≲ 3.5 and 3.5 ≲ z ≲ 4.5, resulting in the observed AB magnitudes I, Z, Y, and J for the first and Z, Y, J, and H for the second. For the SED fitting procedure, we used the CANDELS official photometric catalogs (Galametz et al. 2013; Guo et al. 2013) for the sources in CANDELS and the photometric catalogs described in McLure et al. (2018) for the galaxies in the wide-field areas.

4. Parent sample versus selected sample

To determine to what extent our sample (which is mostly based on the presence of the CIII] emission line) is representative of the parent sample of Lyα emitting galaxies, we performed a Kolmogorov-Smirnov test. Because we were unable to determine residuals of skylines or other extraction problems around the position of the Lyα emission line for all the > 300 individual objects in the parent sample, we restricted the parent sample to galaxies where the fundamental parameter, the Lyα EW, was determined with an accuracy of at least 50%. We find that in terms of stellar masses, SFRs, dust content, ages, and β, we cannot reject the hypothesis that the two samples are drawn from the same distribution. For the Lyα EW distribution, we find instead a very low p-value (p = 0.002), meaning that our galaxies are not drawn from the same distribution of the parent sample. The difference can be clearly appreciated in Fig. 6, where we show the Lyα EW distribution of the selected sample in green and that of the parent sample in yellow. The median Lyα EWs are 21.4 Å and 13.4 Å, respectively. This difference was expected because we selected our sample mostly on the basis of the presence of the CIII] λ1909 Å, which was necessary to evaluate the systemic redshift, and given the apparent correlation between the strength of Lyα and that of CIII] (see below). However, the Lyα EW distribution spans the same range of values as the parent sample and only changes in shape.

Fig. 6. Lyα EW distributions for the sources in the parent sample (yellow histogram) and those in the selected sample (green empty histogram).

Along with the fact that the physical properties of the two samples are drawn from the same statistical distribution, we can therefore reasonably conclude that the results that we find in this analysis for the selected sample are also valid, at least in qualitative terms, for the general population of galaxies with Lyα emission. As we show in the next section, the Lyα EW and the physical properties in both samples show the same correlations. This corroborates our results as well.

5. Results

To study how the kinematics and appearance of Lyα emission depend on the physical properties of the galaxies and on the kinematics of the ISM, we investigated the relations between the spectral and physical properties of the galaxies we evaluated in the previous section. In particular, to give a non-parametric statistical measure of the correlations between the different parameters, we evaluated the Spearman rank correlation coefficients (ρ). We show in Table 1 the Spearman rank correlation coefficients for the variables that we find (anti)correlated, and the significance of these relationships. We trusted the correlation only when the significance was higher than 2σ.

5.1. Correlations between Lyα EW and physical properties

In this subsection we discuss the main correlations that we find between the Lyα EW and the physical properties of the galaxies. We first discuss the results obtained for the selected sample, then we compare them with the corresponding relations for the parent sample.

EW(Lyα) versus M^*. We show in Fig. 7 (blue points) the relation between EW(Lyα) and stellar mass M^*: the two parameters are moderately anticorrelated with a Spearman rank correlation coefficient of −0.34. While the lack of objects with high masses and large LyαEW is surely real, we might be losing objects with low masses and small EWs (lower left part of the figure) because of the flux limit of the survey. We verified that this not the case because of the sensitivity of our observations to emission lines and the median R-band magnitude of our targets.

Fig. 7. Correlations between the Lyα EW and the stellar mass (top panel) and between the Lyα EW and the color excess E(B − V) (bottom panel) for the selected sample (blue) and the parent sample (green).

The existence of this relation is quite controversial because it has previously been observed in some works (e.g., Finkelstein et al. 2007; Jones et al. 2012; Hathi et al. 2016) but not by others (e.g., Kornei et al. 2010). In particular, Pentericci et al. (2007, 2010) found that galaxies without Lyα in emission tend to be more massive and dustier than the rest of the LBG sample and found a significant lack of massive galaxies with high EW(Lyα), which could be explained if the most massive galaxies were either dustier and/or if they contained more neutral gas than less massive objects. However, the authors did not find a precise one-to-one correlation. Finally, Du et al. (2018) found a clear relation for stacks of spectroscopically confirmed LBGs at z ∼ 4, but they showed that this becomes weaker as the redshift decreases with a almost flat trend at z ∼ 3.

EW(Lyα) versus E(B − V). As shown in the bottom panel of Fig. 7, we also observe a moderate anticorrelation between dust extinction, indicated by the color excess E(B − V), and EW(Lyα) with a Spearman rank correlation coefficient of −0.43, in agreement with several previous studies (Shapley et al. 2003; Kornei et al. 2010; Pentericci et al. 2010; Jones et al. 2012; Erb et al. 2016; Hathi et al. 2016). This is also consistent with theoretical expectations: Lyα emission is easily quenched by dust because Lyα photons can be absorbed by dust grains that then reemit them as thermal emission.

EW(Lyα) versus galaxy age. We observe a correlation between Lyα EW and the mass-weighted galaxy age, in the sense that galaxies with a higher Lyα EW are older than galaxies with lower Lyα EW. The Spearman rank correlation coefficient is 0.40. While the relation between Lyα and E(B − V) is now quite well established, the relation between Lyα and galaxy age is still highly debated. Some studies claim that Lyα emitting galaxies are mostly young galaxies in the very early stages of formation (e.g., Pentericci et al. 2007; Du et al. 2018, while others claim that Lyα emission is typically found in older galaxies (e.g., Shapley et al. 2003; Kornei et al. 2010. Finally, other studies find no apparent relation between Lyα EW and age (e.g., Pentericci et al. 2009, 2010.

According to our observations, the brightest Lyα emitters are also the oldest. We know that Lyα emission typically comes from populations of stars with ages < 10 Myr. A Lyα luminous phase in a galaxy can be a long-lasting or recurring phenomenon, however. In the case of galaxies with estimated ages log(Age) >8.5 that also have relatively strong Lyα(EW), we can therefore invoke a scenario in which the oldest galaxies have experienced a strong burst of star formation in the past that caused massive stellar winds, with the consequent expulsion of neutral gas and dust. This resulted in a reduced covering fraction and therefore a smaller attenuation for the observed Lyα photons, which are produced by a second or recurring episode of star formation. On the other hand, the very young galaxies still have more dust, and their Lyα is therefore still more attenuated. A similar scenario is supported by various works of IRAC-detected versus non-detected LAEs (e.g. Lai 2008; Guaita et al. 2010).

We remark that the ages are one of the most uncertain parameters in SED fitting also because they are somewhat dependent on the SFH adopted (see, e.g., Carnall et al. 2019; Leja et al. 2019): we tested the exponentially declining SFH and the delayed exponentially declining SFH instead of a constant one, but did not find significant differences in the observed correlation. Using a more complex SFH might be more physically motivated, but this is beyond the scope of this paper.

EW(Lyα) versus SFR. We do not find any correlation between the Lyα EW and the SFR, measuring a Spearman rank correlation coefficient of 0.14, with very high confidence. This might contradict the results from several previous studies (e.g., Kornei et al. 2010; Hathi et al. 2016) that instead found that these two quantities are anticorrelated, with the brightest LAEs having lower SFRs. However, all these studies analyzed galaxies with Lyα both in absorption and in emission, while we only selected galaxies with Lyα emission. When we restrict previous results to only those sources with Lyα in emission, no clear correlations are observed. This suggests that galaxies with Lyα emission have on average lower SFRs when they are compared to the non-emitting galaxies, but that there is no particular trend with the actual strength of Lyα. Moreover, while the sources presented in Hathi et al. (2016) are in the same mass range as those analyzed in this paper, the sources in Kornei et al. (2010) have much higher average stellar masses (log(M/M_⊙)∼9.9).

EW(Lyα) versus β slope. We reproduce the well-known anticorrelation between the UV spectral slope β and the Lyα EW, with a Spearman rank correlation coefficient of −0.44 and high significance (> 3σ). This relation is shown in Fig. 8 for the sources in the selected sample in blue and for those in the parent sample in green.

Fig. 8. Correlation between the β slope and the Lyα EW.

The UV β slope depends on several parameters including dust, age, and metallicity, although it is thought to be primarily driven by the dust content (e.g., Reddy et al. 2010; Hathi et al. 2016). The anticorrelation between Lyα EW and the β slope is therefore a direct consequence of the fact that less dusty galaxies tend to show brighter Lyα emission, as also probed by the anticorrelations that we observe between the Lyα EW with E(B − V).

We finally tested whether the relations between EW(Lyα) and the physical properties are also valid for the parent sample, which contains galaxies with similar physical properties but different distributions of Lyα EWs and β slopes. The results of this analysis are also shown in Fig. 7. The correlation between EW(Lyα) and stellar mass and E(B − V) is comparable to the correlations found for the selected sample, and similarly, there is no relation between SFR and EW(Lyα). We do not observe any correlation between the Lyα EW and the galaxy mass-weighted age for the galaxies in the parent sample, however, which was instead present for the sources in the selected sample. We remark, however, that as we mentioned above, the galaxy age is the most uncertain parameter from SED fitting and that the Spearman rank correlation coefficients and their significance are evaluated without taking the uncertainty in the individual parameters into account. Finally, a correlation between the Lyα EW and β slope is also observed for the parent sample, although it is weaker than the correlation for the selected sample (ρ = 0.2 with > 2σ significance).

5.2. Correlations between Lyα EW and spectral properties

We discuss in this subsection the relations obtained between the Lyα EW and the other galaxy spectral properties.

EW(Lyα) versus EW(CIII]). We observe a strong correlation between EW(Lyα) and EW(CIII]) with a Spearman rank correlation coefficient of 0.62, which we show in the top panel of Fig. 9. This relation has previously been investigated by other studies using both stacks and individual spectra with somewhat contradictory results. Shapley et al. (2003) found a positive trend between the two quantities on spectral stacks of a sample of LBGs at z ∼ 3. More recently, Stark et al. (2014) found a tighter relation between EW(Lyα) and EW(CIII]) in individual galaxies in a sample of 17 strongly lensed low-mass galaxies at z ∼ 2.

Fig. 9. Correlations between the Lyα EW and the CIII] EW (top panel) and between the Lyα EW and the Lyα spatial extent (bottom panel).

In contrast to these results, other studies just show tentative correlations with a very large scatter. For instance, Rigby et al. (2015) studied a sample of galaxies at 1.6 <  z <  7 and found an apparent correlation, although with appreciable scatter especially for the weakest Lyα and CIII] emitters (EW(Lyα)≲50 Å and EW(C III]) ≲5 Å). A very similar large scatter was found by Le Fèvre et al. (2019) for individual star-forming galaxies at 2 <  z <  3.8 with CIII] emission in VUDS, with a clear trend that is only observed using stacks as in Shapley et al. (2003). Finally, we mention that Du et al. (2018), using stacks of z ∼ 2 star-forming galaxies, observed a nearly null correlation between EW(CIII]) and EW(Lyα), with the exception of the strongest Lyα bin, where they observe a higher EW(CIII]) than in the other subsamples. All these works suggest that on average, the presence of Lyα (CIII]) emission increases the probability of detecting the CIII] (Lyα) emission as well, but that there is not a clear one-to-one relation. There are indeed numerous examples of galaxies with CIII] emission that show only a weak Lyα emission or no emission at all, and several LAEs that do not show any sign of CIII] λ1909 Å emission line (e.g., Le Fèvre et al. 2019). This behavior is expected. The production of both Lyα and CIII] photons indeed strongly depends on the ongoing star formation, but there are significant differences in the production and radiative transfer processes of the two lines (see, e.g., Nakajima et al. 2018, for the CIII production). In our data we find a rather strong correlation between the Lyα EW and CIII] EW in individual galaxies. The correlation is unchanged when we exclude the point in the upper right part of Fig. 9 that might appear to drive it. We note that our sample contains galaxies with higher stellar masses than were reported for instance by Stark et al. (2014), who probed a mass range of 9 × 10⁶ − 1.3 × 10⁹ M_⊙. Clearly, the tight relation that we observe is somewhat biased because we preselected the sample to have both Lyα and CIII] emission. We therefore miss objects with CIII] (Lyα) emission and no Lyα (CIII]) emission. A similar result was also found by Amorín et al. (2017) for a sample of ten low-mass galaxies at z ∼ 3, selected for their strong Lyα emission (EW(Lyα)> 45 Å) in VUDS. In a future work we will exploit the whole VANDELS dataset to evaluate the unbiased relation between Lyα and CIII].

EW(Lyα) versus Lyα_ext and UV_ext. We observe a mild anticorrelation between EW(Lyα) and the UV and Lyα spatial extents, in the sense that galaxies with higher EW appear to be more compact in size both in the UV continuum and in Lyα. The Spearman rank test gives correlation coefficients of −0.37 and −0.39, respectively. The relation between Lyα EW and UV size has previously been observed in LBGs (e.g., Jones et al. 2012). The relation with Lyα_ext has previously been observed by Guaita et al. (2017) using stacks of VUDS galaxies at z ∼ 3, although their definition of Lyα extension was slightly different than ours because they subtracted the UV continuum in quadrature from the Lyα_ext to implicitly consider the deconvolution with the PSF of the observations. We now find a similar correlation also using individual objects.

We tested these correlations also on a small subsample of 17 galaxies for which we were able to measure all the quantities listed in Sect. 3. We find the same correlations and the same significance as observed before for this sample.

5.3. Correlations with IS and Lyα shifts

We tested whether the velocity shifts of the ISM and the Lyα emission line are related to other galaxy parameters. While we do not find any correlation between these shifts and other physical or spectral properties, we find a correlation between these two parameters that we show in Fig. 10. The Spearman rank correlation coefficient for this relation is 0.44, with a significance that is > 2σ. This relation has previously been observed by Guaita et al. (2017) using stacks of z ∼ 3 star-forming galaxies and by Steidel et al. (2010) on a sample of ∼90 galaxies at z ∼ 2. We discuss this trend and its interpretation in the next section. We also show in Fig. 10 the distributions of the IS and Lyα velocity shifts in yellow and green, respectively. The mean values are IS_shift = −100 km s⁻¹ and Lyα_shift = 358 km s⁻¹. These values are very similar to those found by Steidel et al. (2010), although their sample covers a slightly broader range of values than those of our sample (⟨IS_shift ⟩ = − 164 km s⁻¹ and ⟨Lyα_shift ⟩ =400 km s⁻¹).

Fig. 10. Left panel: anticorrelation between the velocity of the main peak of the Lyα profile with respect to the systemic redshift and the IS velocity shift. Right panel: histogram of the measured (centroid) velocities of interstellar absorption lines in yellow and Lyα emission in green with respect to the galaxy systemic redshift derived from the centroid of either the CIII]λ1909 Å or the HeIIλ1640 Å emission lines. The yellow sample contains only the 26 objects for which we were able to give a precise measure of the IS velocity shifts, and the green sample contains the 46 spectra where we were able to estimate the Lyα peak shift.

As in Steidel et al. (2010), we do not observe any correlation between IS and Lyα shifts and other physical parameters. In particular, we expected some trend between IS shift and SFR, but the Spearman rank correlation coefficient is only −0.15, with a significance that is ≪2σ. This trend is well established in low-redshift samples (e.g., Weiner et al. 2009; Heckman et al. 2015) and a relation would also be expected at high redshift because outflows are believed to be driven by pressure that is created by supernovae, stellar winds, and the radiation field caused by star formation. However, our SFRs are derived from SED fitting and have uncertainties: to properly test this relation, an SFR derived from Hα line flux would be needed. In addition to the SFR, we also specifically tested the correlation between outflow velocity and SFR surface density as defined in Heckman et al. (2011) and between outflow velocity and specific SFR defined as SFR/M^*, but again we find no significant trend in our sample.

Finally, we remark that the only tentative correlation (with significance < 3σ) found by Steidel et al. (2010) for the IS velocity offset was with the total baryonic mass, which contains the contribution of gas and stars in galaxies. However, we cannot test this relation on our sample because our spectra do not cover the wavelength range of the Hα emission line, which is needed to estimate the gas mass.

6. Discussion: origin of the relation between ISM and Lyα velocity shifts

The correlation observed in Fig. 10 is the most interesting result that we find in this analysis because it relates the Lyα line kinematics to that of the ISM. According to theoretical predictions, the Lyα kinematics depends on the velocity of the ISM, but also on the HI column density, which gives a measure of the amount of neutral gas along the line of sight. In this section we interpret the observed correlation between the Lyα velocity shift and the IS outflow velocity.

Several detailed studies of the resonant radiative transfer of Lyα photons have been conducted in the past decades. In most of these studies, the source of UV and Lyα photons is assumed to be surrounded by an envelope of neutral gas and dust that can be modeled as a simplified homogeneous shell or as a clumpy medium. The impact of the geometry, the clumpiness, and the multi-phase nature of the ISM in affecting the observed Lyα line has been emphasized since the very early studies (Neufeld 1991; Charlot & Fall 1993; Giavalisco et al. 1996). Hydrodynamical simulations suggest a fragmented HI morphology (e.g., McCourt et al. 2018), coupled with the now well-established clumpiness of the ISM in our galaxy (Stutzki & Guesten 1990; Marscher et al. 1993). This seems to indicate that a clumpy ISM is the best representation of the real ISM in galaxies. Many studies of Lyα radiative transfer have included clumpy media in their prescriptions (e.g., Dijkstra & Kramer 2012; Laursen et al. 2013; Duval et al. 2014; Gronke & Dijkstra 2014). A clumpy ISM has also been invoked to explain the anomalously high EWs that are sometimes seen in high-redshift LAEs because Lyα can preferentially escape from a dusty multi-phase interstellar medium if most of the dust lies in cold neutral clouds (Hansen & Oh 2006). On the other hand, models with very simple and homogeneous geometries, such as an expanding and spherical shell, can in some cases reproduce the observed Lyα spectra better (Verhamme et al. 2008; Yang et al. 2017). The reason for this seems to be that most the models have focused on moderately clumpy media, while simulations indicate larger numbers of smaller clumps (McCourt et al. 2018). In particular, Gronke et al. (2016) showed through hydrodynamical simulations that the predicted Lyα spectra emerging from clumpy geometries with many clouds along the line of sight are very similar to the spectra emerging from galaxies with simplified homogeneous media.

For these reasons, we decided to use the model that was first described in Verhamme et al. (2006). It represents the galaxy as a central region of young stars that are surrounded by an expanding, spherical, homogeneous, and isothermal shell of neutral hydrogen. This model has been extensively used to interpret the kinematics of LAEs (e.g., Verhamme et al. 2008, 2015; Guaita et al. 2016), and it is also particularly well suited for our analysis because it accurately predicts the quantities considered in this paper, such as Lyα shift, Lyα EW, and spatial extension as functions of the ISM kinematics and HI column density.

In this scenario, small Lyα shifts and high (but lower than ∼300 km s⁻¹) outflow velocities are predicted only for low HI column densities (∼10¹⁹ cm⁻²), while large Lyα shifts with low or null outflow velocities can only be observed in case of high HI column densities (∼10²⁰ − 10²¹ cm⁻²) because the Lyα photons must have undergone many scatterings before they can escape the galaxy, resulting in higher Lyα velocity offsets with respect to the systemic redshift. According to this, we expect the objects in the lower left part of Fig. 10 to be characterized by lower HI column densities and the objects in the upper right part of the plot by higher HI column densities.

This scenario is also consistent with the distribution of the Lyα spatial extension with respect to the UV continuum of our sources. We show in Fig. 11 the same plot as in Fig. 10, but color-coded for the Lyα spatial extent with respect to the UV extension (Lyα_ext/UV_ext). The objects that appear to have lower HI column densities (from the Lyα versus IS kinematics) also show a more compact Lyα spatial profile on average, while the objects that should have a higher HI column density (from the Lyα kinematics) also have a larger Lyα extension than the UV extensions. These observations are perfectly in agreement with the shell-model predictions because broader Lyα spatial profiles with respect to the continuum are expected when the Lyα photons undergo many scatterings before they escape the galaxy (caused by higher HI column densities), while more compact spatial profile are predicted for galaxies with lower HI column densities because the Lyα photons can escape galaxies in a more direct way (Verhamme et al., in prep.). Assuming that this model works for our sample of galaxies and that the trend observed in Fig. 10 is due to a different HI column density, we then try to assess whether N(HI) and other galaxy physical properties might be related. In particular, we show in Fig. 12 the same relation between IS_shift and Lyα_shift color-coded for the stellar mass, the attenuation of the dust, indicated by the color excess E(B − V), the SFR, and the Lyα EW. None of these physical parameters appears to correlate with the HI column density, that is, we do not find any significant change in the physical properties in the galaxies as we move from the lower left part of the figure (low HI column densities) to the upper right part (high HI column densities). A possible limitation of this analysis is that we probe only a relatively small range for these quantities. For example, our galaxies span only a factor of 10 in mass: it is possible that if a weak correlation between N(HI) and mass exist, it would emerge only by probing much larger value ranges.

Fig. 11. Lyα velocity offset versus IS velocity offset with respect to the systemic redshift. The plot is color-coded for the Lyα spatial extension with respect to the UV continuum, as shown in the color bar at the right side. The star in the plot is the stack of the sources with no individual ISshift measure. Its color is given by the median of the measured Lyαext/UVext of the IS undetected sample.

Fig. 12. ISshift versus Lyαshift color-coded for different galaxy properties. From left to the right and top to bottom we show the stellar mass, E(B − V) color excess, SFR, and EW(Lyα). The star in each plot is the stack of the sources with no individual ISshift measure. The color is given by the median values of the measured quantity that is indicated in the color bar for the galaxies in the IS undetected sample.

We do not observe any trend between N(HI) and Lyα EW either. A tentative anticorrelation between these quantities has been observed in local galaxies (e.g., Pardy et al. 2014; Yang et al. 2017). From a theoretical point of view, this would only be expected if the intrinsic EW or the gas-to-dust ratio were approximately constant. In the first case, with the increase of N(HI), the path length of Lyα photons would increase and so would their chance of being absorbed by dust, while in the second case, an increase in N(HI) would also lead to an increase in the dust optical depth, suppressing the Lyα EW. We do not observe this relation. This does not contradict our results, however, because we do not have reasons to believe that the conditions to observe this relation are met in our sample of galaxies. From Fig. 12 we do not observe any relation between dust attenuation and the HI column density, suggesting a varying gas-to-dust ratio. Moreover, we do not observe any correlation between the Lyα peak shift and EW(Lyα) relation that would be expected only if the Lyα EW were correlated with N(HI) because higher HI column densities cause more scattering of the Lyα photons, which shifts the line peak to redder wavelengths.

Figure 12 shows that galaxies with brighter Lyα emission on average show no outflows and Lyα peak shifts around 400 km s⁻¹. It is not yet clear if there is a real physical motivation for these observations. This is still under investigation.

7. Summary and conclusions

We have investigated the correlations between the physical properties, the Lyα properties, and the kinematics of the ISM in a sample of 52 star-forming galaxies selected for their Lyα emission and systemic redshift from the VANDELS survey. We find that the Lyα EW and the CIII] λ1909 Å EW are strongly related, in agreement with previous observational studies at similar redshifts (e.g., Stark et al. 2014), suggesting an intrinsic connection between the production of the two lines. We note, however, that for this correlation the sample is biased because we only selected galaxies with both Lyα and CIII] λ1909 Å. This relation could therefore change if the entire population of galaxies with CIII] λ1909 Å emission were considered. We also observe a good anticorrelation, although with some scatter, between the Lyα EW and the UV and Lyα spatial extents, with compact sources being the brightest Lyα emitters. The Lyα line also appears to be brighter in galaxies with lower mass that are less dusty, while no correlation appears to exist with the SFR. We finally observe a very interesting correlation between the Lyα velocity shift and IS shift in the sense that galaxies with high ISM outflow velocities also show small Lyα velocity shifts.

According to the shell model of Verhamme et al. (2006), this relation can be explained if the galaxies with small Lyα shifts also have low HI column densities: the morphology of these galaxies indeed shows compact Lyα_ext/UV_ext, in agreement with expectations. Compact Lyα spatial profiles with respect to the continuum are indeed expected for galaxies with small N(HI) because the Lyα photons did not undergo many scatters with the neutral gas while escaping the galaxies (Verhamme et al., in prep.). On the other hand, objects without ISM outflows show large Lyα shifts (∼500 km s⁻¹). In the context of the shell model, this is only expected to be observed for galaxies with high HI column densities (∼10²⁰ − 10²¹ cm⁻²). These galaxies also show large Lyα_ext/UV_ext, in agreement with expectations because in these cases the Lyα photons scatter many times before escaping, which results in broader spatial profiles.

We finally do not observe any correlation between the trend observed for the IS and Lyα shifts (and therefore N(HI)) and any other physical properties such as EW(Lyα), SFR, Mass, sSFR, Σ_SFR, and dust content. To better understand these correlations and assess their validity, we plan to enlarge the analysis we carried out so far to the entire VANDELS dataset, which will be available in the near future, because our work here was only applied to half of the final VANDELS sample.

We remark in conclusion that the sources we presented here were selected as LBGs and not as LAEs. They are therefore on average brighter in the UV continuum than the objects selected with a narrow-band technique that are tuned to the Lyα line. This means that we did not probe the least massive sources, which are thought to be the most common galaxies in the very early Universe, and which may also be those that contributed most to the reionization photon budget. The correlations that we find might therefore be stronger (or different) if they were investigated in a sample that also included these lower mass galaxies. We plan to test this hypothesis by collecting a sample of strong LAEs that are selected based on narrow-band studies to compare them with the VANDELS galaxies in our sample here and to determine whether there are significant differences in the correlations that we observe based on the sample selection.

Acknowledgments

We thank the ESO staff for their continuous support for the VANDELS survey, particularly the Paranal staff, who helped us to conduct the observations, and the ESO user support group in Garching. FM and LP acknowledge financial support from Premiale 2015 MITiC. We would like to thank Alice Shapley and Max Pettini for their valuable comments on an earlier version of this manuscript.

References

All Tables

All Figures

	Fig. 1. Redshift distribution (left panel) and z magnitude distributions (right panel) of the sources in the selected sample (green empty histogram) and in the parent sample (yellow histogram).
In the text

Fig. 2. Different regions of the spectrum of object CDFS122764 at z = 3.49 are shown for reference. Upper panel: UV continuum between about 1250–1750 Å rest frame, where the most important absorption features are clearly visible, lower left panel: Lyα line profile, and lower right panel: CIII] λ1909 Å emission. We used the magenta parts of the spectrum to perform the fits. The dark blue lines in the upper panel and in the lower right panel are the best fits, and the shaded region in the lower left panel represents the normalized integrated flux of the Lyα that we used to evaluate the EW of the line. The vertical green lines show the positions of the lines at rest-frame air wavelength. The blue vertical dashed line in the lower left panel indicates the position of the Lyα peak. The cyan line is the VANDELS noise spectrum.

In the text

	Fig. 3. Two-dimensional spectrum of object CDFS122764 that we show to visualize the method that we used to evaluate the Lyα and the UV extents. The large green rectangles show the spatial regions of the two-dimensional spectrum that we used to estimate the background, while the yellow (red) square around the Lyα (UV continuum) is the region where we collapsed the spectrum.
In the text

Fig. 4. Top panel: posterior probability distributions of the physical parameters derived for object CDFS020583 at z = 3.5. From left to the right we show the logarithmic mass in units of solar masses, the logarithmic metallicity in units of solar metallicity, the dust attenuation coefficient, the age, and the logarithmic specific SFR. The gray bands are the 1σ confidence intervals. The red crosses are the 1σ confidence intervals obtained without the prior on the metallicity. The different blue contours indicate the 1σ, 2σ, and 3σ confidence levels. Bottom panel: spectral energy distribution of the same galaxy. The blue points are the photometric observations, and the orange shaded regions are the 68% confidence intervals derived by BEAGLE. The black solid lines are ten of the theoretical SEDs with the highest likelihood values.

In the text

	Fig. 5. Star formation rate as a function of the stellar mass for the objects in our sample (green points) and the parent sample (yellow points). The blue line is one of the latest star formation main sequence relation evaluated by Santini et al. (2017) for galaxies in the first four HST Frontier Fields in the redshift range 3 ≲ z ≲ 4.
In the text

	Fig. 6. Lyα EW distributions for the sources in the parent sample (yellow histogram) and those in the selected sample (green empty histogram).
In the text

	Fig. 7. Correlations between the Lyα EW and the stellar mass (top panel) and between the Lyα EW and the color excess E(B − V) (bottom panel) for the selected sample (blue) and the parent sample (green).
In the text

	Fig. 8. Correlation between the β slope and the Lyα EW.
In the text

	Fig. 9. Correlations between the Lyα EW and the CIII] EW (top panel) and between the Lyα EW and the Lyα spatial extent (bottom panel).
In the text

Fig. 10. Left panel: anticorrelation between the velocity of the main peak of the Lyα profile with respect to the systemic redshift and the IS velocity shift. Right panel: histogram of the measured (centroid) velocities of interstellar absorption lines in yellow and Lyα emission in green with respect to the galaxy systemic redshift derived from the centroid of either the CIII]λ1909 Å or the HeIIλ1640 Å emission lines. The yellow sample contains only the 26 objects for which we were able to give a precise measure of the IS velocity shifts, and the green sample contains the 46 spectra where we were able to estimate the Lyα peak shift.

In the text

Fig. 11. Lyα velocity offset versus IS velocity offset with respect to the systemic redshift. The plot is color-coded for the Lyα spatial extension with respect to the UV continuum, as shown in the color bar at the right side. The star in the plot is the stack of the sources with no individual ISshift measure. Its color is given by the median of the measured Lyαext/UVext of the IS undetected sample.

In the text

Fig. 12. ISshift versus Lyαshift color-coded for different galaxy properties. From left to the right and top to bottom we show the stellar mass, E(B − V) color excess, SFR, and EW(Lyα). The star in each plot is the stack of the sources with no individual ISshift measure. The color is given by the median values of the measured quantity that is indicated in the color bar for the galaxies in the IS undetected sample.

In the text