{'fraction_non_alphanumeric': 0.03388508290427358, 'fraction_numerical': 0.0279137057997532, 'mean_word_length': 4.410482695111657, 'pattern_counts': {'":': 0, '<': 84, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Keck Spectroscopy of Faint 3<z<7 Lyman Break Galaxies: III. The Mean Ultraviolet Spectrum at z=4

We present and discuss the mean rest-frame ultraviolet spectrum for a sample of 81 Lyman Break Galaxies (LBGs) selected to be B-band dropouts with a mean redshift of z=3.9 and apparent magnitudes z_AB<26. Most of the individual spectra are drawn from our ongoing survey in the GOODS fields with the Keck DEIMOS spectrograph, and we have augmented our sample with published data taken with FORS2 on the VLT. In general we find similar trends in the spectral diagnostics to those found in the earlier, more extensive survey of LBGs at z=3 undertaken by Shapley et al (2003). Specifically, we find low-ionization absorption lines which trace the presence of neutral outflowing gas are weaker in galaxies with stronger Lyman-alpha emission, bluer UV spectral slopes, lower stellar masses, lower UV luminosities and smaller half-light radii. This is consistent with a physical picture whereby star formation drives outflows of neutral gas which scatters Lyman-alpha and gives rise to strong low-ionization absorption lines, while increasing the stellar mass, size, metallicity, and dust content of galaxies. Typical galaxies are thus expected to have stronger Lyman-alpha emission and weaker low-ionization absorption at earlier times (higher redshifts). Indeed, our mean spectrum at z=4 shows somewhat weaker low-ionization absorption lines than at z=3 and available data at higher redshift indicates a rapid decrease in low-ionization absorption strength with redshift. We argue that the reduced low-ionization absorption is likely caused by a decrease in the covering fraction and/or velocity range of outflowing neutral gas at earlier epochs. Our continuing survey will enable us to extend these diagnostics more reliably to higher redshift and determine the implications for the escape fraction of ionizing photons which governs the role of early galaxies in cosmic reionization. [Abridged]


INTRODUCTION
Considerable progress has been made over the past decade in charting the demographics of high redshift galaxies. Multi-wavelength surveys have defined the luminosity functions of UV and sub-mm selected starforming sources (Reddy & Steidel 2009;Wardlow et al. 2011) as well as the coeval population of quiescent massive red galaxies (Brammer et al. 2011). Spitzer data has revealed the time-dependent stellar mass density -a complementary quantity which represents the integral of the past star formation activity (e.g. Stark et al. 2009). Through these surveys a well-defined picture of the history of star formation and mass assembly over 0 < z < 6 has been empirically determined (Hopkins & Beacom 2006;Ellis 2008;Robertson et al. 2010). The redshift range 2 < z < 3 corresponds to the peak of star formation activity where the Hubble sequence starts to emerge, and the earlier era corresponding to 3.5 < z < 5 is an even more formative one where mass assembly was par-ticularly rapid.
Intermediate dispersion spectroscopy of carefullyselected Lyman break galaxies (LBGs) has been particularly important in defining population trends that cannot be identified from photometric data alone. A very influential study at z ≃ 3 was undertaken by Shapley et al. (2003) who used composite Keck LRIS spectra of various subsets of nearly 1000 LBGs to examine the role of hot stars, Hii regions and dust obscuration, as well as to measure the outflow kinematics and absorption line properties of neutral and ionized gas. Composite spectra are particularly useful for measuring weak lines which cannot be studied in detail for individual objects. Through these careful studies, a detailed picture of the mass-dependent evolution of LBGs has emerged (see Shapley 2011 for a recent review).
In earlier papers in this series , hereafter Paper I; Stark et al. 2011, hereafter Paper II), we introduced an equivalent spectroscopic survey of LBGs selected from a photometric catalog of more distant LBGs with 3 < z < 7 in the Great Observatories Origins Deep Survey (GOODS) fields Stark et al. 2009). Whereas the Shapley et al. (2003) study targetted the study of LBGs close to the peak of activity in the overall cosmic star formation his-tory, this earlier period corresponds to a less well-studied era when the rate of mass assembly is particularly rapid. From photometric data alone, Stark et al. (2009) deduced some significant changes in the characteristics of star formation at z ≃ 4 − 6 compared to later times, for example a shorter timescale of activity (≃ 300 Myr). We considered it crucial to understand these changes in LBG properties if these galaxies are to be used as probes of cosmic reionization at higher redshifts.
At the time of writing, our Keck survey is continuing with increasing emphasis at high redshift. Paper I presented the first substantial results from a survey of LBGs at 3.5 < z < 6 observed with the Keck/DEIMOS spectrograph. Paper II augmented this data with a further sample following more ambitious exposures focusing primarily on z ≃ 6 LBGs. Incorporating a sample of ESO VLT spectra, retrospectively selected using similar photometric criteria as those for the Keck sample from the FORS2 study of Vanzella et al. (2005Vanzella et al. ( , 2006Vanzella et al. ( , 2008Vanzella et al. ( , 2009, the current dataset amounts to a sample of 546 galaxies over the redshift range 3.5 < z < 6.3. Our earlier papers in this series concentrated primarily on the rate of occurrence of Lyman α (Lyα ) emission in our spectra (the "Lyα fraction"). The overall goal was to understand the significantly different evolutionary trends in the luminosity functions of LBGs and narrow-band selected Lyα emitters (LAEs, Ouchi et al. 2008) prior to the use of the Lyα fraction as a test of when reionization ended (Schenker et al. 2011). Paper I confirmed a result found by Shapley et al. (2003) at lower redshift, namely that Lyα emission is more frequent in lower luminosity LBGs rising to a high proportion ≃ 50% at M UV = −19. More importantly, the Lyα fraction was found to rise modestly with redshift over 3 < z < 6. By correlating the visibility of Lyα emission with UV continuum slopes derived from the HST photometry, it was argued that these trends in the visibility of line emission most probably arise from different amounts of dust obscuration. Reduced dust extinction in lower luminosity LBGs and those at higher redshift has also been deduced from studies of larger photometric samples (Bouwens et al. 2009, see also Reddy & Steidel 2009).
Paper I also discussed the possibility that the covering fraction of hydrogen may be lower in low-luminosity LBGs. Strong Lyα emission in luminous LBGs is often associated with low equivalent width interstellar absorption lines arising from a non-uniform covering fraction of neutral hydrogen (Quider et al. 2009;Shapley et al. 2003). This trend suggests that the high Lyα fraction in faint LBGs is partially due to a lower covering fraction, which would imply that Lyman continuum photons may more easily escape from intrinsically faint z ≃ 4 − 6 galaxies. Such a result would have great importance in understanding the role of z > 7 galaxies in maintaining cosmic reionization (Robertson et al. 2010).
The present paper represents our first analysis of the spectral properties of z ≃ 4 − 5 LBGs derived from composite spectra in the manner pioneered at z ≃ 3 by Shapley et al. (2003). The large database now amassed following the campaigns at Keck and the VLT makes a similar study now practical in the redshift range where there is evidence of increased short-term star formation and the mass assembly rate is particularly rapid. Via detailed studies of low-ionization absorption line and emis-sion line profiles, we aim to examine possible changes in the kinematics and covering fraction of neutral gas, which affects the strength of Lyα and the escape fraction of ionizing photons, as we approach the reionization era. As our redshift survey continues, in this paper we focus on a sample of galaxies with 3.5 < z < 4.5, selected as B-dropout LBGs. Combining our Keck sample with data from the VLT (see Paper I for details), herein we examine the spectral features and trends in composite spectra drawn from a sample of 131 galaxies.
A plan of the paper follows. We briefly review the spectroscopic observations and their data reduction in §2; much of the relevant discussion is contained in Paper I. In §3 we describe the selection of individual spectra that we consider appropriate for forming the composite mean spectrum at z ≃ 4 and the associated sources of uncertainty. §4 examines the mean spectrum in detail and introduces the various diagnostic features in the context of a physical model for LBGs of different masses and starformation rates. In §5 we compare spectroscopic trends grouped by observable properties such as mass and luminosity with those found at z ≃ 3 by Shapley et al. (2003). In §6 we discuss those trends which appear to be redshift-dependent, discussing implications for the role of early star-forming galaxies in cosmic reionization. Finally we summarize our results in §7.
Throughout this paper, we adopt a flat ΛCDM cosmology with Ω Λ = 0.7, Ω M = 0.3, and H 0 = 70 h 70 km s −1 Mpc −1 . All magnitudes in this paper are quoted in the AB system (Oke 1974).

OBSERVATIONS AND DATA REDUCTION
The rationale and procedures used to undertake our spectroscopic survey of LBGs over 3 < z < 7 were introduced in detail in Paper I and the interested reader is referred to that paper for further detail. Here we recount only the basic details. Our target LBGs were selected as B, V or i ′ -band "dropouts" based on deep photometry in the two GOODS fields. The photometric catalog used in our analysis will be described in detail in Stark et al. (2011, in preparation). The selection and photometric approach is largely similar to that described in Stark et al. (2009) and Paper I, but we highlight two key updates. First the selection is performed on the v2 GOODS ACS catalogs. Second, we utilize deep ground-based near-IR imaging in GOODS-N obtained from WIRCAM on CFHT (Wang et al. 2010) and deep HST Wide Field Camera 3 / IR imaging of GOODS-S from CANDELS (Grogin et al. 2011;Koekemoer et al. 2011).

Keck/DEIMOS
The majority of spectra we present are taken from an ongoing survey with the DEep Imaging Multi-Object Spectrograph (DEIMOS; Faber et al. 2003) on the Keck II telescope. In this paper we have used observations taken in 2008 April and 2009 March (masks GN081, GN082, GN083, GN094, and GN095 in Paper I), which targeted a total of 261 B-drop and 88 V -drop galaxies. These data were taken with the 600 lines mm −1 grating, covering the wavelength range 4850-10150Å with a resolution of ≃ 3.5Å. We do not use data at wavelengths λ > 9200Å which are affected by strong and variable absorption by atmospheric water vapor. The seeing was typically 0. ′′ 8 and ranged between 0. ′′ 5 and 1. ′′ 0 during the observations. All data were reduced and calibrated using a modified version of the IDL pipeline spec2d, developed specifically for DEIMOS by the DEEP2 survey team (Davis et al. 2003). The data were reduced as described in Paper I, with the addition of two important modifications. First, the continuum traces of all target galaxies and other objects occupying the same slit were carefully masked to exclude object flux from the sky background model. Second, the b-spline fit to the sky background was modified to include a 2nd-order polynomial fit to the spatial dimension. These modifications significantly improved the sky subtraction, particularly at long wavengths λ > 7000Å where bright sky lines can be problematic.
The reduced one-and two-dimensional spectra were visually inspected using the IDL program SpecPro (Masters & Capak 2011). Spectra which suffered from poor data quality were excluded from further analysis. The exclusions included unacceptable amounts of scattered light within the detector mosaic, defective CCD columns, contamination from bright nearby sources, or poor sky subtraction, for example arising from the location of the slit on the detector. A few low-redshift interlopers were also identified and excluded from further analysis. For the remaining spectra, redshifts were measured from either Lyα emission (where present) or interstellar absorption lines. Galaxies identified as hosting strong AGN based on the presence of Civ or other strong emission lines were excluded. We identified 5 stars, 2 dusty low-redshift galaxies (z ∼ 0.5), 4 AGN, and 134 star-forming galaxies with secure redshifts z > 3. After rejecting poor data, the final sample consists of 94 high-quality spectra with accurate redshifts. Examples of high redshift dropout spectra are shown in Figure 1. 2.2. Archival VLT/FORS2 Spectroscopy in GOODS-S To augment our sample of high redshift spectra, we have also made use of data from the FORS2 program of Vanzella et al. (2005Vanzella et al. ( , 2006Vanzella et al. ( , 2008Vanzella et al. ( , 2009 which targeted dropouts in the GOODS-S field. The characteristics of that survey in terms of resolution and spectral coverage are very similar to that undertaken at Keck with DEIMOS and details can be found in Paper I. Using the coordinates provided in the published FORS2 database we queried the version 2.0 ACS catalogs for GOODS-S and undertook our own photometric measures and dropout selection criteria in an identical fashion to that used for our Keck survey. The magnitude distribution of the FORS2 sample is generally weighted towards sources brighter than those in the overall Keck survey, but for the purposes of constructing the mean spectra discussed in this paper, the bulk of the individual spectra are of comparable brightness.

Redshift measurements
Care is needed in deriving accurate systemic redshifts from rest-frame UV spectra since the strongest features trace the kinematics of outflowing gas rather than that of the stars. Stellar absorption lines are usually too faint to be measured precisely given the signal to noise of the spectra. Typically the only features detected in individual spectra are Lyα and strong interstellar absorption lines such as Siiiλ1260, Oiλ1302+Siiiλ1304, Ciiλ1334, Siivλλ1393,1402, Siiiλ1526, and Civλλ 1548,1550. Absorption by outflowing gas results in a blueshift of interstellar absorption lines. Outflowing neutral hydrogen along the line of sight leads to a Lyα profile which displays very broad ( ∼ > 1000 km s −1 ) blueshifted absorption and net redshifted emission. The highionization Siiv and Civ lines also arise in P-Cygni stellar winds with broad blueshifted absorption. The magnitude of these offsets has been well-quantified for star-forming galaxies at z ≃ 2.3, and is typically −200 km s −1 for interstellar absorption lines and +500 km s −1 for Lyα emission .
To determine accurate redshifts for making composite spectra, we restrict our sample to those with redshifts measured from either Lyα emission (z Lyα ) or lowionization interstellar absorption lines (z IS ). Although, as discussed above, these are not at the systemic redshift, we follow well-established techniques to correct for the typical offsets ). We do not use redshifts based on Lyα absorption or high-ionization lines (Siiv and Civ) because of the complex and variable blueshifts of these features with respect to the systemic redshift. We define z Lyα as the centroid of the emission line and consider only spectra in which the line is detected at > 5σ significance. Interstellar absorption line redshifts require careful treatment to avoid spurious identification of sky line residuals. We consider only the lowionization Siiiλ1260, Oiλ1302+Siiiλ1304, and Ciiλ1334 features which are typically found in the highest signalto-noise regions of our spectra. Absorption features at longer wavelengths are less reliable due to the higher density of strong night sky lines, while shorter wavelength transitions are lost in the Lyα forest. To measure absorption line redshifts, we first estimate the redshift from Lyα (either in emission or absorption). We then fit Gaussian profiles to the spectrum near the expected position of the three features above. We require the best fit of all three lines to be consistent to within ±500 km s −1 with a combined significance of > 5σ. If these conditions are met, we define z IS as the weighted mean redshift of the three interstellar features.
With the criteria above, suitable redshifts are available for a total of 131 high-quality DEIMOS and FORS2 spectra. 91 redshifts are based on measures of z Lyα emission only, 31 have z IS only, and 9 have both z Lyα and z IS . The distribution of redshifts and absolute UV magnitudes for this sample is shown in Figure 2. Unless stated otherwise, further analysis in this paper is restricted to the 81 sources with redshift z < 4.5 and apparent magnitude z ′ AB < 26.0. We applied this additional magnitude criterion in order to ensure a well-defined continuum signal/noise in each individual spectrum. 53 of these 81 sources are drawn from the Keck survey and 28 from FORS2.

Sample bias
While the sample is 90% complete for dropout selected galaxies to an apparent magnitude z AB = 25, at fainter magnitudes Lyα emission or strong interstellar absorption features are required for reliable redshifts. This results in a bias towards stronger low-ionization absorp-FIG. 1: Examples of individual DEIMOS spectra for 3.5 < z < 4.5 LBGs. From top to bottom: a galaxy at z = 3.73 with Lyα emission and weak interstellar absorption features; a galaxy at z = 3.906 with Lyα and interstellar absorption; a faint galaxy at z = 4.07 with Lyα emission and no detectable absorption lines; and galaxy at z = 3.778 with Lyα emission whose spectrum is contaminated by scattered light and detector defects. The top three spectra are included in stacking analyses while the bottom spectrum is excluded based on poor data quality. The continuum signal-to-noise in each spectrum degrades noticeably at wavelengths λ > 7000Å due to increased OH sky emission. tion lines for the fraction (≃ 35%) without Lyα . We can quantify this bias through the detectability of the average low-ionization absorption line strength shown in the composite spectrum in Figure 3 (see §3). We find that we can measure redshifts for this average absorption line strength with 90% completeness at z AB = 24.7, declining through 50% at z AB = 25.2 to zero at z AB > 25.5. Of the galaxies with detected absorption line redshifts, those with z AB > 25.2 have line strengths only 10% stronger than for those with z AB < 25.2. Clearly this is a small effect.
The bias toward stronger Lyα emission for fainter galaxies is manifest in Figure 2 where a paucity of objects with W Lyα < 20Å can be seen at faint magnitudes. This bias was fully quantified in Paper I using Monte Carlo simulations. For galaxies at 3.5 < z < 4.5, a sample completeness of 95% is reached at W Lyα ∼ > 20Å for z ′ AB = 26, and W Lyα ∼ > 7Å for z ′ AB = 25. For the sample presented in this paper, the least biased subset is that with strong Lyα emission, followed by that with bright apparent magnitudes z ′ AB < 25. A final issue in considering composite spectra is that these are comprised of individual spectra across 3.5 < z < 4.5 with different rest-frame wavelength ranges. As our spectra generally cover the wavelength range 5000 < λ < 9200Å, galaxies at z = 3.5 contribute to the rest frame 1100 − 2050Å, while those at z = 4.5 contribute to 900 − 1650Å. At longer rest-frame wavelengths the composite spectrum is therefore largely contributed by galaxies at lower redshift. In addition there are wavelength-dependent sources of noise discussed in §3.2. In summary, the data used to construct the composite spectrum in Figure 3 correspond to a mean redshift z = 3.98 at λ rest = 1100Å, z = 3.95 at 1200Å, z = 3.90 at 1500Å, and z = 3.82 at 1650Å. For the wavelength range of interest in this work, the redshift bias ∆z < 0.15 is not particularly troublesome.

COMPOSITE SPECTRA
We now turn to the presentation of the composite spectra. We need to account for the difference between the redshifts determined using Lyα on the one hand and the low-ionization interstellar lines on the other hand and the systemic redshift prior to co-addition. We also seek FIG. 2: Redshifts and UV luminosities for the sample of 131 galaxies with suitably accurate redshifts measured either from Lyα or strong interstellar absorption lines. Points are color-coded according to the equivalent width of Lyα , WLyα . Galaxies brighter than an apparent magnitude z ′ AB < 26.0 (solid line) typically have well-detected continua suitable for forming a composite spectrum. Galaxies with redshift z < 4.5 (dashed line) are used for the composite spectra discussed in §3- §5, while §6 includes galaxies at all redshifts shown here.
to understand the signal/noise of the composite in terms of the statistical uncertainties and the variance among the individual spectra used to construct the composite. Composite spectra are constructed by shifting the individual spectra into the rest frame according to a deduced systemic redshift and then averaging the set. In general terms we will first identify a subsample based on their observable properties. Each spectrum in the sample is shifted to the adopted rest frame and interpolated to a common wavelength scale with a dispersion of 0.12Å. All spectra are normalized to have a median f ν = 1 in the range 1250-1500Å. Spectra taken with DEIMOS are smoothed to a resolution of 1.9Åto match the lower resolution of FORS2 data. The spectra are then averaged at each wavelength using a σ-clipped mean to reject outliers arising from sky subtraction residuals and cosmetic defects. An equal number of positive and negative outliers are rejected at each wavelength, totaling at most 30% of the data. The remaining data are averaged with an arithmetic mean.
Uncertainty in a composite spectrum will arise from both the finite signal to noise and the variance of the individual galaxies. For example, the variance in Lyα equivalent widths in our sample is much greater than the uncertainty measured in the individual spectra. It is especially important to quantify the sample variance for weak features that are generally not detected in individual spectra. We account for sample variance with a bootstrap technique. For each composite spectrum we create 100 alternate composites using the same number of spectra but drawn at random from the parent sample. Each alternate has an average 63% of the sample represented with 37% duplicates. Every measurement made on the composite spectrum is repeated for each of the 100 alternates. We then take the measurement error to be the standard deviation of the 100 alternate measurements, which reflects both the sample variance and finite signal to noise.
As discussed in §2.3 the most challenging issue is to determine the systemic redshift prior to shifting to the rest-frame. Here we follow the approach used by Shapley et al. (2003). As a first approximation we use the value of z Lyα (where available) to construct a composite spectrum. This enables us to locate the stellar photospheric line Ciiiλ1176 in the composite where we detect a velocity difference of −330 km s −1 with respect to Lyα . We can thus infer that Lyα emission in our sample is redshifted on average by ∆v Lyα = 330 km s −1 . In a similar fashion, stacking spectra using the redshift z IS based on low-ionization interstellar absorption results in a detection of Ciiiλ1176 with a velocity offset of +190 km s −1 . For comparison, at z ≃3 Steidel et al (2010) find ∆v Lyα = +445 km s −1 and ∆v IS = −164 km s −1 . To construct composite spectra, we use either z Lyα shifted by −330 km s −1 or z IS shifted by +190 km s −1 to approximate the systemic redshift of each galaxy. We use the Lyα -based redshift when available since it is typically determined with greater precision than z IS . Figure 3 shows the composite spectrum of 81 galaxies in our sample with 3.5 < z < 4.5 and apparent magnitude z ′ AB < 26.0 using this method.

Uncertainties in the Systemic Redshift
A natural concern is the extent to which these applied shifts might vary within the sample used to make the composite. This can be estimated from observations of higher signal to noise from spectra taken at lower redshift. Steidel et al. (2010) quantify the offset between z IS , z Lyα , and the systemic redshift z Hα in a sample of 89 galaxies at z ≃ 2.3. They find velocity offsets ∆v IS = −170 ± 115 km s −1 and ∆v Lyα = 485 ± 185 km s −1 relative to Hα. Assuming this is representative of our data, the uncertainty in our systemic redshift is therefore likely to be σ(v) ∼ 150 km s −1 . An upper limit on σ(v) can be estimated from the width of spectral lines in the composite spectrum. In particular the restwavelength of the stellar line Ciiiλ1176 in the composite provides a valuable measurement of the average offset from the systemic velocity and its width provides an upper limit on the effective spectral resolution. We measure a systemic velocity of 21 ± 101 km s −1 in the composite spectrum ( Figure 3) and a FWHM = 520 km s −1 (deconvolved from the instrumental resolution ≃ 450 km s −1 ). The uncertainty in the adopted redshifts about the true systemic stellar value is therefore ≤ 520 km s −1 FWHM or equivalently σ(z) ≤ 220 km s −1 , comparable in fact to what was achieved for individual spectra at lower redshift by Steidel et al. (2010). Figure 3 shows the composite spectrum of our sample as well as the 1σ error spectrum derived using the bootstrap technique discussed above. The error at each pixel is calculated as the standard deviation of all averaged data points (excluding outliers), divided by the square root of the number of data points. There are several wavelength-dependent factors contributing to the error spectrum in addition to the finite signal-to-noise of individual spectra. One factor is the intrinsic sample variance, seen clearly as a noise spike at the position FIG. 3: The composite spectrum of 81 galaxies in our sample with 3.5 < z < 4.5 and apparent magnitude z ′ AB < 26. The effective mean redshift for the sample averaged over wavelength is z = 3.9. The strongest spectral features are labelled. The gray filled region shows the ±1σ error at each point, determined from the scatter of individual spectra used to create the composite. The error spectrum peak at 1216Å is due to large scatter in the intrinsic distribution of Lyα equivalent widths. The error is lowest at ∼ 1300 − 1500Å where the continuum signal to noise ratio is ≃ 30. The error increases at shorter wavelengths where the instrument throughput is lower, and at longer wavelengths where sky emission is much stronger.

Error spectrum
of Lyα , and also evident for Ciiλ1334 and other absorption features. Another is the decreased instrument throughput at λ ∼ < 6000Å, leading to higher noise at shorter wavelengths. Similarly, stronger sky line emission causes increasing noise at longer wavelengths. Finally, the number of contributing spectra peaks at restframe λ ≃ 1150 − 1350Å , with increased noise at higher and lower wavelengths where fewer spectra are included. Our observed wavelength range 5000 < λ < 9200Å corresponds to a rest-frame 1000 < λ < 1800Å at the mean redshift z = 3.9 of our magnitude-limited sample. Ultimately we acheive a signal to noise ratio in the continuum of ≥ 10 between the Lyα line and rest-frame 1800Å and ≃ 5 in the Lyα forest, with a peak S/N ≃ 30 at 1350Å.

FEATURES IN THE COMPOSITE SPECTRUM
We now discuss the composite spectrum ( Figure 3) in more detail, focusing on the strong spectral features at 1215-1550Å where we have the best signal to noise. In this wavelength range we detect Lyα , Nvλλ1239,1243, Siiiλ1260, Siii*λ1265, Oiλ1302+Siiiλ1304 (blended), Siii*λ1309, Ciiλ1334, Siivλλ1394,1403, Siii λ1527, Siii*λ1533, and Civλλ1548,1550 at high significance. Our z = 3.9 composite spectrum is very similar to the composite of z = 3 LBGs presented in Shapley et al. (2003), which we show in Figure 4 for comparison. The absolute magnitude distribution of our sample is broadly similar to that at z = 3; both cover the range −22 < M UV < −20. Our discussion below follows closely that presented originally by Shapley et al. (2003) but we are also interested in whether there are differences seen over the redshift range 3 < z < 4.5. We will discuss these possible evolutionary trends in §5.
4.1. Lyman Break Galaxies: a physical picture It is helpful to begin by describing a possible physical picture of LBGs based on many analyses of the extensive observations at z ≃ 3 (see Shapley 2011 for a recent review). Typical L* LBGs at z = 3 have ultraviolet halflight radii ≃ 2.0 kpc, stellar masses ∼ 3 × 10 10 M ⊙ , and star formation rates ∼ 50 M ⊙ yr −1 (Bouwens et al. 2004;Ferguson et al. 2004;Shapley et al. 2001). The star formation surface density is sufficient to drive "superwinds" of outflowing gas, similar to those seen in local galax- FIG. 4: Composite spectrum of all 81 galaxies in our sample with z < 4.5 and apparent magnitude z ′ < 26, compared to the composite spectrum of 811 LBGs at z = 3 presented in Shapley et al. (2003). The higher redshift sample has a much stronger Lyα forest break, slightly redder UV spectral slope, slightly stronger Lyα emission, and weaker absorption lines. The inset shows a zoom-in of the region from 1200 − 1600Å, which contains most of the absorption lines of interest in this paper. Heckman 2002). Indeed, blueshifted interstellar absorption lines confirm there are outflows with typical velocities of ≃ 150 km s −1 , and in some cases as high as 800 km s −1 Pettini et al. 2002;Quider et al. 2009Quider et al. , 2010. However, the physical origin of these outflows remains unclear (e.g. Murray et al. 2010). It is thought these outflows produce an extended circumgalactic medium (CGM) of ejected material. Outflowing gas is found in both low and high ionization states (e.g. Siii and Siiv). Low-ionization transitions such as Siii and Cii are mostly associated with neutral hydrogen whereas high ionization lines occur in the fully ionized component. Based on trends in the strength of the various species, Shapley et al. (2003) suggest a geometry in which discrete clouds of neutral gas are embedded in a halo of ionized gas. Steidel et al. (2010) show that the neutral and ionized CGM components both extend to radii of at least 125 kpc.

Lyα
Lyα is the most prominent and diverse feature in our individual spectra. The line originates from hydrogen recombination in H II regions photoionized by massive stars and, as these stars dominate the adjacent contin-uum, its intrinsic equivalent width should be within the range W Lyα = 100 − 200Å for nearly all stellar populations (e.g. Forero-Romero et al. 2011). Our composite spectrum reveals a complex line profile with strong absorption extending blueward to v = −4000 km s −1 , and redshifted emission to v = +1000 km s −1 ( Figure 5) with a peak offset of v = +330 km s −1 relative to the systemic velocity. The net equivalent width W Lyα = 21 ± 3Å is only ∼ < 20% of the expected intrinsic value. The form of this Lyα profile has been readily understood in terms of the physical picture discussed in §4.1. Lyα emission produced at the systemic velocity can escape only along a line of sight free of neutral hydrogen or if it is shifted in velocity far from resonance. As photons escape they encounter blueshifted clouds of partially neutral gas, which absorb and re-emit isotropically. Photons backscattered from neutral clouds at small radii will appear redshifted and have a higher chance of escaping. Additionally, photons scattered from neutral gas at the edge of the CGM can escape into the ionized intergalactic medium and will be observed as blueshifted emission . Indeed, Figure 5 shows such a weak blueshifted emission peak in the composite spectrum. On the other hand, the reduced equivalent width compared to that expected intrinsically could be due to many effects including dust extinction, scattering at large radii where emission falls outside the spectroscopic slits and a non-zero escape fraction f esc of ionizing radiation. Since these effects cannot easily be disentangled, the absolute strength of the Lyα line must be interpreted with caution (Paper II, Schenker et al. 2011).
We can estimate the spatial extent of the neutral CGM (i.e. the radius at which escaping Lyα is last scattered) based on the strength of the diffuse blueshifted Lyα emission. The spatial profile of extended Lyα emission has been well-quantified at z ≃ 2.6 by Steidel et al. (2011) who find that diffuse haloes are a generic property of star-forming galaxies leading to a total Lyα flux ≃ 5× greater than that measured within the spectroscopic slits.
Average Lyα surface brightness profiles are well-fit at large radii by an exponential form, viz.
where b l ≃ 25 kpc. The blueshifted emission shown in Figure 5 has an equivalent width W = 1.5 ± 0.3Å measured within a slit aperture of ≃ 1. ′′ 0×1. ′′ 0, corresponding to 7.1 × 7.1 kpc 2 at z = 4. Assuming this approximates the peak value, we estimate Σ 0 = 1.5 ± 0.3Å per 50 kpc 2 and integrating Equation 1 yields a total equivalent width of W halo = 2πb 2 l Σ 0 . Using the measured value of Σ 0 , the scale length is given by Although the total Lyα flux is not measured, we can estimate its value from theoretical expectations as well as observations at lower redshift. The equivalent width in Figure 5 is W = 21.0Å (excluding the blueshifted emission component), thus we can write W halo = W tot − 21 A. As discussed above, the observed W tot is diminished by dust and the escape of ionizing radiation. Indeed, Steidel et al. (2011) measure W tot = 17 − 93Å in various subsamples of their data, all lower than the expected intrinsic value, W tot =100-200Å . Considering the composite DEIMOS spectra shown in Figure 5 and assuming an intrinsic W tot = 135Å (Forero-Romero et al. 2011) we derive upper limits of b l < 20 ± 2 kpc for galaxies with W Lyα > 0, and b l < 17±2 kpc for those with W Lyα > 20 A. This latter constraint is the most stringent, and also likely closest to the true value of b l . Steidel et al. (2011) find that their sample of Lyα emitters (LAEs, defined as W Lyα > 20Å) has the highest W tot (= 93Å) and largest scale length (b l = 28.4 kpc) of any sub-sample that they analyze. Our constraint of b l < 17 kpc for the LAEs therefore suggests that the characteristic size of Lyα haloes may be smaller at z = 4 than at z = 2 − 3. Due to the large inherent uncertainties, direct measurements of low surface brightness Lyα emission are needed to confirm this possibility.

Low-ionization metal transitions
According to the physical picture in §4.1, absorption in low-ionization transitions occurs in both the interstellar medium and outflowing clouds of cool gas. The former is at the systemic velocity while the latter is unblended transitions of Siii at 1260 and 1527Å whose equivalent width ratio is a valuable tracer of the optical depth . The profiles of both lines and their ratio are shown in Figure 6. Gas near the systemic velocity (|v| < 200 km s −1 ) is clearly optically thick, while for v < −200 km s −1 it is intermediate in optical depth. This difference is also seen in highresolution spectra of lensed galaxies (Pettini et al. 2002;Quider et al. 2010), and optically thin gas is seen at large galactocentric radii  suggesting that smaller, optically thin clouds are more easily accelerated to high velocity and large distances. We note that the column density at which both Siii transitions become optically thick (τ > 1) is N Siii = 1.9 × 10 12 cm −2 . Assuming log(Siii/H) = −4.86 as measured for the lensed galaxy cB58 (Pettini et al. 2002), we estimate that optically thin absorption occurs in clouds with hydrogen column densities N Hi ∼ < 10 17 cm −2 .

Fine structure transitions
A satisfying aspect of our composite spectrum is the successful identification of fine structure emission lines of Siii (see marked features in Figure 3). Siii ions in the CGM absorb photons in the resonance transitions and immediately re-emit a photon at approximately the same velocity. A photon absorbed at 1260Å will be re-emitted at either the same wavelength or to the fine structure transition Siii*λ1265. Likewise a photon absorbed at 1527Å will be re-emitted as either Siiiλ1527 or Siii*λ1533, and an absorbed Siiiλ1304 photon can be reemitted as Siii*λ1309. In all three cases the probability of emission in the resonant and fine structure transitions is approximately equal.
Since no absorption is seen in the fine structure transitions, we infer that atoms in the excited ground state will typically decay to the ground state before absorbing another photon. Since every absorbed photon is reemitted, the net equivalent width of the resonant and fine structure transitions is W Siii + W Siii* = 0. The precise equivalent widths depend on the initial absorption W Siii,abs and the optical depth. If the gas is optically thin, re-emitted photons will immediately escape giving W Siii,em = W Siii*,em = −0.5W Siii,abs . In the limit of optically thick gas, resonant photons will be continuously scattered until they emerge as Siii* (after 2 scatterings on average). In this case, W Siii,em = 0 and W Siii*,em = −W Siii,abs . Since the majority of absorption occurs in optically thick gas ( §4.3), we expect the equivalent width of Siii absorption lines to reflect the kinematics and covering fraction of neutral gas with minimal contamination from Siii re-emission.

The spatial extent of low-ionization absorption
From the picture above ( §4.3.1), we expect that the equivalent width of the fine structure emission lines should be equal and opposite to the resonant line equivalent width. This is not the case: W 1265 /W 1260 = −0.56± 0.36 and W 1533 /W 1527 = −0.52 ± 0.20. (W 1309 /W 1304 is more difficult to quantify since Siiiλ1304 is blended with Oiλ1302). There are two possible explanations. One is that scattered photons have larger path lengths and are subject to greater dust attenuation. Since the UV continuum slopes imply little differential extinction (mean E(B-V) = 0.10) and we expect only 2 scatterings for the average Siii* photon, this seems unlikely. More reasonably, the emitting region could be larger than that sampled by our slits. Although the 1. ′′ 0 × 1. ′′ 0 slit aperture samples most of the continuum light and its line of sight absorption, the CGM which absorbs and isotropically re-emits scattered photons extends to much larger radii (Steidel et al. , 2011. The strength of fine structure emission provides a direct measurement of the amount of absorption at small radii contained within the slit. The fine structure to resonant absorption line ratio suggests that a fraction 0.53 ± 0.17 (combining both λ1265 and λ1527 measures) of the total fine structure emission is contained within the extraction aperture of 1 arcsec 2 . The half-light radius of Siii* emission, and hence Siii absorption, thus corresponds to ∼ 0. ′′ 5 = 3.5 kpc at z = 4. This scale is only slightly larger than the median half-light radius r h = 0. ′′ 31 of galaxies in our sample, and would be reached in only 20 Myr at the typical outflow velocity (190 km s −1 ).
It is instructive to reconsider the z = 3 composite spectrum from Shapley et al. (2003) where the fine structure lines are clearly seen (but were not interpreted fully along the discussion above at the time). The Siii* lines are noticeably stronger in our z ≃ 4 composite, with an average W Siii* = 0.7 ± 0.2 at z ≃ 4 c.f. 0.30 ± 0.05Å at z = 3. Furthermore, the Siii absorption lines are weaker in the z = 4 composite, suggesting the absorption takes place at larger radii at z = 3. This difference cannot be explained through instrumental differences between the two surveys. Smoothing our z = 4 composite to match the 3.25Å resolution of Shapley et al. (2003) reduces the equivalent widths by < 5%. Likewise the wider slits used for the z = 3 data (1. ′′ 4, 11 kpc) c.f. the z = 4 data (1. ′′ 0; 7 kpc) is not the cause. If the CGM properties are similar, we expect the fine structure strength to constitute a larger fraction of the resonant absorption line strength in the z = 3 data, contrary to observations. Defining R F S = −(W 1265 + W 1533 )/(W 1260 + W 1527 ), we find R F S = 0.53 ± 0.17 for the z = 4 composite spectrum and R F S = 0.16 ± 0.04 for z = 3, indicating a smaller characteristic radius of fine structure emission at z = 4.
Inescapably, therefore, we conclude the circumgalactic gas around LBGs at z = 4 differs in two important ways from that at z = 3: there is greater low-ionization absorption at small radii at z = 4 and less total lowionization absorption. In addressing the origin of this effect, Figure 10 shows that the difference in equivalent width is at least partially linked to the kinematic offset from the Lyα emission. Thus it seems a higher covering fraction at small radii is required to produce the stronger fine structure emission. Higher resolution observations of individual bright, or perhaps gravitationally-lensed, galaxies at z ≃4 will ultimately be required to separate the relative contributions of covering fraction and kinematics in explaining this result.

High-ionization lines
The high-ionization lines Siiv and Civ arise both in interstellar gas and in stellar P-Cygni winds. The velocity centroid of Siiv, v = 140 km s −1 , is consistent with the low-ionization interstellar absorption lines suggesting that most absorption is interstellar in origin. This is supported by the absence of a significant redshifted emission component expected for a P-Cygni profile. In contrast, Civ is broader with a larger absorption velocity offset v = −370 km s −1 and redshifted P-Cygni emission indicating a large contribution from stellar winds. We also detect the Nvλλ1239,1243 P-Cygni feature, although the proximity to Lyα makes this feature difficult to study in detail.
We can determine the optical depth of highly ionized outflowing gas from the ratio of Siiv absorption lines. The ratio W Siivλ1394 /W Siivλ1403 = 2.0 for optically thin absorption, and 1.0 in the optically thick case. Shapley et al. (2003) find optically thin absorption in composite spectra of LBGs at z = 3, whereas we measure a ratio 1.4 ± 0.4 in the composite spectrum shown in Figure 3 indicating a significant contribution of optically thick absorption at z ∼ 4. The total equivalent width of the Siiv doublet is weaker by a factor 0.72 ± 0.12 in Figure 3 compared to the z = 3 composite of Shapley et al. (2003). The combination of higher optical depth and lower equivalent width suggests that the velocity range and/or covering fraction of the ionized gas traced by Siiv is lower at higher redshift.

Metallicity
The P-Cygni profile of Civ is sensitive to metallicity, and the combination of Civ and Heii equivalent widths constrains both the age and metallicity. Here we compare the equivalent widths measured in the composite spectrum ( Figure 3) with theoretical models in order to estimate the typical metallicity of galaxies in our sample. We note that the equivalent width of Civ contains significant interstellar absorption, so the value reported in Table 1 should be treated as an upper limit on the P-Cygni component. Assuming the interstellar absorption component of Civ is similar to that of Siiv (W ≃ 1.0Å), we take the P-Cygni absorption component to have equivalent width −1.6 ± 1.0Å with a conservative uncertainty. We compare this estimate and the measured equivalent width of Heii (reported in Table 1) with predictions from the stellar population synthesis code BPASS presented in Eldridge & Stanway (2009). We consider BPASS models which include binary evolution with continuous star formation rate, and determine the difference ∆W between observed and predicted equivalent width as a function of metallicity and relative carbon abundance. We restrict the stellar population age to be within the 1σ scatter of the median value for galaxies in our sample, determined to be 10 8.5±0.6 years from spectral energy density fits assuming constant star formation with a Kroupa (2002) initial mass function. Figure 7 shows the resulting relative error between BPASS models and measured equivalent width, which we define as Relative error = 0.5 ∆W Civ σ Civ 2 + ∆W Heii σ Heii 2 . (3) Relative error values ∼ < 1 are thus consistent with the data.
Models with Z ≤ 0.004 and somewhat depleted carbon abundance are in good agreement with the data. We note that observations of both local and high-redshift galaxies indicate typical relative carbon abundances (C/O) ≃ 2 − 5× lower than the solar value (equivalent to X C = 0.2 − 0.5 in Figure 7) for metallicities Z ∼ < 0.004 (Kobulnicky & Skillman 1998;Shapley et al. 2003;Erb et al. 2010;Eldridge & Stanway 2011). Solar metallicity models (Z = 0.020) do not fit the data. We conclude that the typical metallicity of galaxies in our sample is Z ∼ < 0.004 or ∼ < 0.2× solar metallicity. Bright galaxies at z < 3.8 in the FORS2 sample have gasphase metallicity 12 + log(O/H) = 7.7 − 8.5 or about 0.1 − 0.6× the solar value, in reasonable agreement (Maiolino et al. 2008). These values are consistent with the allowed metallicity range 0.1 < Z/Z ⊙ < 0.6 inferred for the Shapley et al. (2003) composite using the same BPASS models (Eldridge & Stanway 2011). Given that the equivalent widths of Civ and Heii in Shapley et al. (2003) agree well with those in Table 1, we expect the typical metallicity of galaxies to be similar in both samples.

SPECTROSCOPIC TRENDS
We now turn to an analysis of how the various spectroscopic features discussed in §4 are related to observable properties of LBGs as a prelude to considering how they might evolve with redshift. We will begin with the dependence of low-ionization absorption line strength with Lyα equivalent width. This is motivated by the common physical dependence of these features -both are governed by the kinematics and covering fraction of neutral circumgalactic gas -and also because this has been examined in detail for LBGs at z = 3 .
We construct composite spectra of three sub-samples of galaxies divided according to their Lyα equivalent width. Defining W LIS as the average equivalent width of λ1260, λ1303, λ1334, and λ1527Å features, Figure 8 shows weaker W LIS for galaxies with stronger Lyα . This was also noted for z = 4 LBGs by Vanzella et al. (2009). Both the z=3 and z=4 samples show this trend, although W LIS is weaker at z = 4 for galaxies with strong Lyα emission. This difference can be attributed  Figure 3 (see text for details). Relative error values ∼ < 1 indicate good agreement. Model results for a range of metallicity Z and carbon depletion factor XC as a function of age are described in Eldridge & Stanway (2011). The values XC = 0.5 (0.1) correspond to a carbon abundance reduced by a factor of 2 (10) relative to solar abundance ratios. Lines are plotted for the median age of our sample (10 8.5±0.6 years) determined from spectral energy density models with the same initial mass function and star formation history used in BPASS, and error bars indicate the allowed range for ages within 1σ of the median. Measurements of C/O abundance at both low and high redshift indicate carbon depletion factors XC ≃ 0.2 − 0.5 in galaxies with Z < 0.2 Z⊙ (e.g. Erb et al. 2010). For this range of XC , models with Z ∼ < 0.004 (equivalent to Z = 0.2 Z⊙) are in good agreement with the data. Solar metallicity models are unable to reproduce the observed equivalent widths.
to the luminosity-dependent trend of stronger W Lyα and weaker W LIS in fainter galaxies (Paper I; Shapley et al. 2003;Vanzella et al. 2009). If we consider a subset of galaxies in our sample with similar absolute magnitudes to those observed at z = 3 (−21.5 < M UV < −21.0), we recover the same normalization ( Figure 8). Our sample is 90% complete at the corresponding apparent magnitudes, so we expect a negligible bias.
Having established the correlation of W LIS with W Lyα , we now examine the dependence with other demographic properties. We divide the full sample into two bins of equal size according to each property of interest. The results are shown in Figure 9. We briefly review the trend of each property with W LIS and W Lyα and discuss the physical origin. Many of these trends were seen in Shapley et al. (2003) and Vanzella et al. (2009) and are clearly inter-related due to correlations between the demographic properties. Less luminous galaxies have stronger W Lyα and weaker W LIS illustrating that higher star formation rates drive larger amounts of neutral gas into the CGM with higher velocity and/or covering fraction. Defining the ultraviolet spectral slope as β = 5.30 · (i AB − z AB ) − 2.04 for B-drops (Bouwens et al. 2009), we also find that bluer galaxies with lower β have stronger W Lyα and weaker W LIS . This trend was also noted in Paper I, which showed that LBGs with strong Lyα emission have sys-tematically bluer β. Since neutral gas also presumably contains dust, the same gas which gives rise to W LIS also reddens the continuum. Stellar masses M * are measured for 60% of our sample for which there is unconfused Spitzer/IRAC photometry. As expected from the trends with luminosity, lower mass galaxies have stronger W Lyα and weaker W LIS . Finally, measuring half-light radii r h from the GOODS ACS data with Sextractor, we find smaller galaxies have stronger W Lyα and weaker W LIS . Could these trends could be due to selection effects or sample bias? The trend of stronger Lyα for less luminous LBGs is of particular concern, since fainter galaxies will require a larger W Lyα for detection. To address this, we consider only the 32 galaxies with apparent magnitudes z ′ AB < 25 for which the spectroscopic sample is 90% complete. This sub-sample is divided into two equal bins and the results confirm that trends seen in the larger sample also hold in brighter galaxies unaffected by sample bias.
The trends shown in Figure 9 are generally consistent with the overall trend of W LIS with W Lyα (Figure 8). The W Lyα variations mostly arise from the distribution of neutral gas (traced by W LIS ), with relatively little effect from other demographic properties examined (M UV , β, M * , r h ). However, demographic properties do have some effect. The strongest deviation seen in Figure 9 is that with β, which shows stronger W LIS than would be expected at a given W Lyα for galaxies with red UV slopes. Noted also by Shapley et al. (2003), this suggests that outflowing neutral gas contains dust which reddens the continuum. Also, more luminous galaxies (i.e. those with higher star formation rates) have stronger W LIS at a given W Lyα (Figure 8). This is likely due to increased absorption at large velocities in galaxies with higher SFR, as observed at z ≃ 1.4 (Weiner et al. 2009). For all other demographics, the composite spectra have values of W LIS within 1σ of that expected purely based on W Lyα .
In summary, the trends seen in Figure 9 arise almost entirely because of variations in the neutral gas covering fraction and/or kinematics, which are themselves correlated with the demographic properties.

Kinematics
We have now established that outflowing neutral gas is the dominant factor in determining both W Lyα and W LIS . However, variations from the trend of W LIS with W Lyα are apparent, particularly with β and M UV as discussed above. The most obvious mechanism for this behaviour is a systematic difference in the covering fraction f c and kinematics of neutral gas. Assuming is the outflow velocity), a higher f c and lower velocity range can conspire to give a constant W LIS . However, changing f c and v will also affect the transmission of Lyα photons resulting in a different W Lyα . It is therefore of interest to consider how to distinguish between the covering fraction and kinematics of neutral gas.
The kinematics of neutral gas can be roughly parameterized by the velocity dispersion and centroid of lowionization interstellar absorption lines. We have measure the velocity ∆v of low-ionization absorption lines with respect to Lyα as a proxy for outflow velocity. ∆v is strongly correlated with W Lyα and W LIS in LBGs at z = 3  in the sense that larger ve-  (2003) at mean z = 3, divided into quartiles of WLyα . Galaxies in the z = 3 sample have typical luminosities −21.5 < MUV < −21.0 corresponding to the blue triangles. Galaxies of the same luminosity lie on the same correlation between WLyα and WLIS at both z = 3.9 and z = 3. Fainter galaxies have weaker low-ionization absorption lines at fixed WLyα .
FIG. 9: Equivalent width of low-ionization absorption lines compared to that of Lyα , divided according to observable demographic properties as described in the text. Points binned by WLyα at z = 3 and z = 4 are the same as in Figure 8. The sample was divided into two bins for each demographic property (MUV , M * , β, and r h ), with WLyα and WLIS measured from composite spectra of the galaxies in each bin.
locities are associated with stronger interstellar absorption and weaker Lyα emission. Figure 10 shows ∆v measured from the same composite spectra used to determine demographic trends (e.g. Figure 9), as well as the results from Shapley et al. (2003). All z ≃ 4 composites are consistent (within 1σ) with the relation measured at z = 3 as well as the mean ∆v measured for B-dropout galaxies by Vanzella et al. (2009). We measure a velocity dispersion for each composite, and find that each is within 1.2σ of the effective spectral reso- lution σ = 290 km s −1 (measured for the stellar [Ciii] feature in Figure 3). We are therefore unable to detect trends in outflow kinematics with demographic properties in composite spectra. Higher signal-to-noise data, higher spectral resolution, or detailed studies of individual galaxies are required to address trends in kinematics and covering fraction of neutral gas at z = 4.

THE EVOLVING CGM
We now turn to a discussion of the redshift evolution of the neutral CGM surrounding typical LBGs. The most useful probes are Lyα and low-ionization absorption lines which trace the kinematics and covering fraction of neutral gas, and fine structure emission lines which provide a constraint on the spatial extent of the absorbing gas ( §4.3.2).
The common dependence of Lyα and low-ionization absorption lines on the neutral CGM results in a strong correlation between W LIS and W Lyα . Various physical properties of LBGs are correlated with both W LIS and W Lyα , but in such a way that the relation between W LIS and W Lyα remains nearly constant (Figure 9). Furthermore, the W LIS -W Lyα relation at fixed M UV does not change significantly with redshift between z = 3 and z = 4 (Figure 8).
To further examine evolutionary trends with redshift, we now divide our spectroscopic sample into two bins of redshift at fixed M UV , now including galaxies at all redshifts (no longer restricted to z < 4.5 as in previous sections). We consider galaxies with absolute magnitude −21.5 < M UV < −20.5, chosen to be representative of the sample in Shapley et al. (2003). There are 64 galaxies in our sample within this M UV range (see Figure 2). We construct composite spectra of galaxies with redshift above and below the median z = 4.1 and measure the equivalent width of Lyα and low-ionization lines (both resonant absorption and fine structure emission) in each composite. The results are given in Table 2 along with the demographic properties of galaxies in each subsample. The quantities M UV , M * , and r h are measured from photometry while β is determined from a direct fit to the ultraviolet continuum in the composite spectrum.
We define f λ ∝ λ β and fit the rest frame 1300 − 1700 A to determine β, with uncertainty quantified using the bootstrap method described in §3. Aside from redshift, the demographics of each sub-sample are quite similar. The higher redshift galaxies have slightly higher average M * , smaller r h , and smaller (bluer) β. W Lyα is consistent for both to within the sample variance, and is also consistent with the value W Lyα = 14.3Å measured for the composite spectrum of z = 3 LBGs in Shapley et al (2003). The most striking difference is in the strength of the low-ionization absorption lines, which are significantly weaker at higher redshifts ( Figure 11). The variation in W LIS is not explained by systematic differences in W Lyα or demographic properties, hence we seek an alternate explanation.
We first examine whether the evolution in W LIS could arise as a result of different equivalent width distributions for Lyα . Although the mean W Lyα across our two redshift subsamples is similar, the lower redshift subsample has a broader distribution and contains more galaxies with Lyα in absorption (W Lyα < 0). This is reflected in the larger sample variance in W Lyα at lower redshift (Table 2). We evaluate the effect of this potential bias on W LIS by constructing a composite spectrum from a subset of the z < 4.1 galaxies with intermediate W Lyα = 0 − 30Å, resulting in a consistent mean W Lyα with sample variance reduced by a factor of 2.5. This composite has W LIS = −1.5Å, 0.2Å higher than when the full range of W Lyα is used, but still considerably lower than the value W LIS = −1.0Å measured for the higher redshift galaxies. We therefore conclude that differences in the W Lyα distribution are insufficient to explain the observed variation in absorption line strength with redshift.
There are several possible physical explanations for the evolution of W LIS with redshift shown in Figure 11. One possibility is that the kinematics and/or covering fraction of neutral gas are systematically different. For example, an outflowing wind with fixed input energy and momentum will reach higher velocity at lower redshifts due to lower density of the IGM. However, this effect should be stronger between z = 3.0 − 3.8 than from z = 3.8 − 4.7, whereas the decrement in W LIS is much stronger from z = 3.8 − 4.7 ( Figure 11). Furthermore, we measure a higher offset between the velocity centroid of Lyα emission and low-ionization absorption for z > 4.1 galaxies (∆v = 660 ± 80) than for z < 4.1 (∆v = 550 ± 40). Both measurements are consistent with the trend shown in Figure 10. Based on the trend of W LIS with ∆v seen at z = 3 , we would then expect stronger W LIS at higher redshift. Kinematics are thus unable to explain the difference in W LIS , at least with the information currently available. Data with higher spectral resolution and signal-to-noise are needed to fully address differences in the covering fraction and kinematics of neutral gas.
Another obvious potential cause of weaker W LIS , or equivalently decreased W Lyα at a fixed W LIS , is the presence of neutral hydrogen with very low column density of heavy elements. In this case we would expect optically thin absorption line profiles. This could be directly tested as shown in §4.3.1, but the signal-to-noise of Siiiλ1527 in the high redshift composite is too low to constrain the optical depth. We do see some evidence that the column density of neutral gas is lower at higher redshifts based on Lyα line profiles in the composite spectra, shown in Figure 12. The Lyα absorption trough at ∼ 1200 − 1215Å is significantly weaker in the higher redshift composite. This absorption arises at least in part from damping wings of high column density gas with N Hi ∼ > 10 20 cm −2 associated with the low-ionization metal absorption lines (e.g. Pettini et al. 2000Pettini et al. , 2002. The higher redshift galaxies are therefore characterized by lower typical N Hi or/and a lower covering factor of high column density gas. Higher signal-to-noise data at z ≃ 5 is required to determine whether this affects the optical depth of low-ionization metal transitions.
Weaker low-ionization lines could also be caused by a systematically higher ionization state at higher redshifts resulting in lower column densities of low-ionization gas.
This would also explain the weaker damped Lyα absorption trough seen at higher redshifts (Figure 12). This scenario can be tested by measuring the equivalent widths of higher-ionization silicon transitions, in particular Siiiiλ1206 and Siivλλ1393,1402. We find that the Siiv lines are weaker at higher redshift by a factor of 0.68 ± 0.25 in W Siivλλ1393,1402 , consistent with the decrement in W Siii . We do not detect the Siiiiλ1206 line in the high redshift composite, with a 1σ upper limit of W Siiii < 0.8 times that of the lower redshift composite. Siiiiλ1206 is marked in Figure 12 and is clearly stronger at lower redshift. We therefore find no evidence of a significant change in the ionization state of outflowing gas at higher redshifts.
A final possibility for lower column density of heavy elements is that Lyα is scattered by "cold streams" of nearly metal-free neutral gas, predicted by simulations to accrete onto galaxies at small radii and with increasing rates at higher redshift (e.g. Dekel et al. 2009;Faucher-Giguère et al. 2011). If they exist, such streams could be identified as Hi absorption systems with no corresponding absorption from heavy elements. Steidel et al. (2010) find little evidence for the presence of such streams around LBGs at z = 2 − 3. A similar study at higher redshift would be very interesting but also extremely challenging with current observational facilities.
Could the difference in W LIS with redshift be caused in part by a different physical extent of circumgalactic gas? In §4.3.2 we showed that fine structure emission lines can be used to measure the amount of absorption taking place within the size of the spectroscopic aperture. This quantity does appear to change with redshift, in the sense that absorption takes place at smaller characteristic radius in LBGs at higher redshift. The composite spectra of LBGs at z = 4 has weaker absorption lines and stronger fine structure emission than the z = 3 composite of Shapley et al. (2003) (Figure 4), in qualitative agreement with trends seen at z = 3 ; see their Figure 9). We have only a weak constraint on the redshift evolution at z > 4 due to limited signal-to-noise in the high redshift composite. The ratio of total Siii fine structure emission W Siii * in the highand low-redshift composites is 1.1 ± 0.5 (taken as FIG. 11: Equivalent width of low-ionization absorption lines measured from composite spectra of LBGs at different redshifts. The DEIMOS and FORS2 data presented in this paper are separated into two sub-samples of equal size as described in the text. We also show the equivalent result at z = 3 from the composite spectrum of Shapley et al. (2003), with redshift distribution described in Steidel et al. (2003). The average low-ionization line profile of each composite is shown in the inset. Low-ionization absorption lines are significantly weaker for galaxies in the highest redshift composite. All three samples have consistent mean luminosity and WLyα . with values given in Table 2), while the corresponding ratio of absorption line strength W LIS is 0.6 ± 0.2. The inferred fraction of low-ionization absorption within the slit aperture is thus a factor 1.8±1.0 larger in the higher redshift sample, consistent with no evolution. Higher signalto-noise is required to accurately constrain the spatial extent of low-ionization absorption at z > 4.
To summarize, we find weaker low-ionization absorption lines in LBGs at higher redshift with fixed M UV . The difference in W LIS is not consistent with the trends observed for Lyα and various demographic properties examined in §5. We have discussed several possible causes for this discrepancy including variation in the kinematics, covering fraction, optical depth, ionization state, and spatial extent of absorbing gas. The data show that the ionization state has no significant effect, but we are unable to conclusively address other possible causes. Spectra of galaxies at z > 4 with higher signal-to-noise are required to more accurately constrain the optical depth and spatial extent. Additionally, spectra of individual galaxies taken with higher spectral resolution will be required to independently determine the covering fraction and kinematic structure of absorbing gas.

Galaxy evolution
We argue in §5 that Lyα equivalent width is determined primarily by the neutral CGM, which is correlated with various demographic galaxy properties (Figure 9). It is well established that these demographic properties vary with redshift and we expect W Lyα to vary accordingly. Large photometric surveys have shown that Lyman break galaxies at increasingly higher redshifts z > 3 have lower luminosities (Bouwens et al. 2007(Bouwens et al. , 2011a, bluer UV spectral slopes (Bouwens et al. 2009(Bouwens et al. , 2011b The Lyα absorption trough is weaker in the higher redshift composite indicating a lower incidence of neutral gas with high column density. Siiii absorption is also weaker in the higher redshift composite. and smaller sizes (Bouwens et al. 2004;Ferguson et al. 2004). This is in accordance with inside-out galaxy growth: galaxies increase in size and stellar mass as they evolve with time, while increasing metallicity and dust content reddens the ultraviolet continuum. Simultaneously, star formation drives large-scale outflows of gas which reach larger distances and are accelerated to larger velocities at later times (e.g. Murray et al. 2010). Galaxies which are more evolved (i.e. larger, more massive, redder) should therefore have a CGM characterized by larger spatial extent, larger velocity range, and higher covering fraction of neutral gas. Observationally this results in weaker W Lyα , stronger W LIS , and relatively weaker fine structure emission. These are precisely the trends observed at both z = 3 ) and z = 4 (this work).

Lyα in the epoch of reionization
We reiterate that Lyman break galaxies at increasingly higher redshifts z > 3 have lower luminosities, bluer UV spectral slopes, smaller stellar masses, and smaller sizes. Notably, all trends in the demographics of galaxies at higher redshift are correlated with stronger W Lyα and weaker W LIS (Figure 9). We therefore expect typical galaxies at higher redshifts to have, on average, stronger Lyα emission. Earlier results from this survey confirm that strong Lyα emission is more frequent in galaxies at higher redshift (Paper I; Paper II). We find no evidence of this trend reversing. In fact, galaxies with extremely small size, low mass, and blue β tend to be the strongest Lyα emitters (e.g. Erb et al. 2010). We do, however, expect the average Lyα emission strength to decrease significantly at increasing redshifts in the epoch of reionization due to neutral hydrogen in the IGM (e.g. Haiman & Spaans 1999).
We have devoted considerable discussion to the properties of Lyα in part because Lyα is of great interest as a tracer of cosmic reionization. Several authors have now presented evidence that reionization was incomplete at z ≃ 7 based on a rapidly decreasing fraction of galaxies with strong Lyα emission at z ∼ > 6.5 (Schenker et al. 2011;Ono et al. 2011;Pentericci et al. 2011). Although trends at z = 3 − 4 suggest that the galaxies observed by these authors should have a higher fraction of strong Lyα emission, we have presented evidence that the relation between W Lyα and W LIS is systematically different at z > 4 ( Figure 11). This implies a systematic difference in the spatial, kinematic, or optical depth structure of neutral circumgalactic gas compared to galaxies at lower redshift. The physical origin of this evolution and its effect on W Lyα at z > 4 will need to be understood in order to fully interpret the results of Lyα surveys at higher redshifts in the context of reionization.

SUMMARY
The rest-frame ultraviolet spectra of star-forming galaxies contains a wealth of information about the properties of the circumgalactic medium. In this paper we have presented an analysis of several features which trace the CGM with a focus on the properties of neutral gas. We find that the trends observed at lower redshift (z = 3; Shapley et al. 2003) also hold at z = 4 with approximately the same normalization ( Figure 8). However, we find evidence for rapid evolution at z > 4 with lower W LIS at fixed W Lyα and M UV , suggesting a systematic difference in the spatial distribution, kinematics, ionization state, or optical depth of circumgalactic gas at higher redshifts. We determine that the ionization state is not responsible for the observed evolution but are unable to distinguish between kinematics, covering fraction, optical depth, or the spatial extent of neutral gas as the likely cause. We are collecting additional spectra of LBGs at z > 4 with our ongoing survey, including high spectral resolution observations of bright lensed galaxies from which we can disentangle the kinematic profile and covering fraction of neutral gas. These data will allow us to address the precise magnitude and physical origin of evolution in circumgalactic gas properties.
As a final note, we emphasize that the neutral CGM is of great interest in the context of reionization of the universe. Neutral gas in the circumgalactic medium ab-sorbs ionizing radiation, thereby inhibiting the ability of galaxies to reionize the universe. The escape fraction of ionizing photons is one of the most important and uncertain factors in determining the contribution of star-forming galaxies to reionization (Robertson et al. 2010). We have shown that typical galaxies at higher redshift have weaker low-ionization absorption lines based on their demographic trends, and presented new evidence that absorption lines are systematically weaker at z > 4 even for fixed demographic properties. This is likely caused by a lower covering fraction and/or velocity range of neutral gas, and we will address the physical origin of this evolution with future data from our ongoing survey. Determining the redshift evolution of neutral gas covering fraction in LBGs will be of great interest for interpreting surveys of Lyα emission in the context of reionization and addressing the role of star-forming galaxies in reionizing the universe.