Joint Ly emitters - quasars reionization constraints
We present a novel method to investigate cosmic reionization, using joint spectral information on high redshift Lyman Alpha Emitters (LAE) and quasars (QSOs). Although LAEs have been proposed as reionization probes, their use is hampered by the fact their Ly line is damped not only by intergalactic H but also internally by dust. Our method allows to overcome such degeneracy. First, we carefully calibrate a reionization simulation with QSO absorption line experiments. Then we identify LAEs ( and EW Å) in two simulation boxes at and and we build synthetic images/spectra of a prototypical LAE. The surface brightness maps show the presence of a scattering halo extending up to 150 kpc from the galaxy center. For each LAE we then select a small box of around it and derive the optical depth along three viewing axes. At redshift , we find that the Ly transmissivity , almost independent of the halo mass. This constancy arises from the conspiracy of two effects: (i) the intrinsic Ly line width and (ii) the infall peculiar velocity. At higher redshift, , where the transmissivity is instead largely set by the local H abundance and consequently increases with halo mass, , from 0.15 to 0.3. Although outflows are present, they are efficiently pressure-confined by infall in a small region around the LAE; hence they only marginally affect transmissivity. Finally, we cast LOS originating from background QSOs passing through foreground LAEs at different impact parameters, and compute the quasar transmissivity ( ). At small impact parameters, cMpc, a positive correlation between and is found at , which tends to become less pronounced (i.e. flatter) at larger distances. Quantitatively, a roughly 10 increase (from to ) of is observed in the range . This correlation becomes even stronger at . By cross-correlating and , we can obtain a H density estimate unaffected by dust. At , the cross-correlation is relatively weak, whereas at we find a clear positive correlation. We conclude by briefly discussing the perspectives for the application of the method to existing and forthcoming data.
keywords:intergalactic medium - cosmology: theory - diffuse radiation - reionization - numerical simulation
Lyman Alpha Emitters (LAEs) are galaxies showing prominent emission in the transition of the hydrogen atom resulting in a Ly emission line at Å. At redshift the line can be observed in optical/IR bands, thus allowing searches of distant galaxies traced by such radiation. What is the source of Ly photons? Short-lived, massive stars produce large amounts of UV photons with energy eV = 1 Ryd which ionize the surrounding gas. H-atoms recombine on a short time scale in the dense interstellar medium and 2/3 of the recombination cascade ends up in Ly photons, powering the observed luminosity. The rapid production of dust associated with star formation process, and most noticeably at high redshift by supernova explosions (Todini & Ferrara, 2001), might however lead to important attenuation of the Ly emission line. In spite of this drawback, Ly is still considered the best available tracer of high-redshift star-forming galaxies. This was realized already more than four decades ago (Partridge & Peebles, 1967), but a full-scale application of the method had to await for technological progresses and became routinely used only towards the end of the last century (Hu et al., 1998).
Since then, several hundreds of LAEs over redshift have been detected using narrowband imaging and spectroscopy (e.g., Hu et al. 1998; Rhoads et al. 2003; Malhotra & Rhoads 2004; Taniguchi et al. 2005; Kashikawa et al. 2006; Ouchi et al. 2008, 2009, 2010; Castellano et al. 2010; Pentericci et al. 2011). Redshift corresponds the end of the epoch of reionization, the last major cosmic phase transition in which the gas turned from a neutral to an ionized state. Albeit the general features of the reionization process are now understood and constrained by available data (Choudhury & Ferrara, 2006; Bolton & Haehnelt, 2007; Iliev et al., 2009; Mitra et al., 2011b, a), many physical details as e.g., the initial mass function (IMF), escape fraction of Ly, UV, and ionizing continuum photons, intrinsic source spectral energy distribution, and interstellar (ISM) and intergalactic (IGM) medium inhomogeneities, remain only vaguely known. There are strong hopes that LAEs can allow us to clarify these issues as several authors (e.g., Miralda-Escude & Rees 1998; Santos 2004; Haiman & Cen 2005; Dijkstra et al. 2007; Dayal et al. 2008, 2009) have pointed out.
However, before we can fully exploit the LAE potential in this sense, one has to deal with the fact that reionization itself affects the propagation of Ly photons. In fact, in addition to the aforementioned presence of dust, interstellar and intergalactic H atoms produce an effective opacity caused by their large scattering cross-section to Ly photons. As a result, inferring the intrinsic Ly luminosity of a source from the observed one it often quite difficult. Further complications arise: Dayal et al. (2011) have shown that a degeneracy between the IGM ionization state and the clumping of ISM dust is imprinted in the LAE visibility. Stated differently, a wide range of IGM H fractions, , can reproduce the observed Ly luminosity function with different Ly escape fractions, , from the galaxy, i.e. increased transmission through a more ionized IGM is compensated by higher dust absorption inside the galaxy. As a result, disentangling from the transmissivity through the IGM is very challenging.
To attack this problem we propose here a novel strategy. Adelberger et al. (2003, 2005) introduced a new method to measure IGM metallicity and neutral fraction around a target galaxy using background QSOs or galaxies. Such sources are separated from the target galaxy by an angular distance small enough that their absorption spectra trace the IGM around the galaxy. They measured the spatial distribution of metals and HI using the absorption spectra of background QSOs/galaxies at redshift . We exploit this successful method to investigate the transmissivity of LAEs through the IGM during reionization. Several tens of QSOs are detected at redshift higher than (e.g.Fan et al. 2001, 2003, 2006a; Songaila 2004; Willott et al. 2007, 2010) and new QSOs are still being discovered (Mortlock et al., 2011) in the United Kingdom Infrared Telescope (UKIRT) Infrared Deep Sky Survey (UKIDSS).111http://www.ukidss.org With hundreds of LAEs in the samples at , it is possible to derive joint constrains on the IGM during cosmic reionization.
The main goal of this paper is to break dust-HI the degeneracy affecting the visibility of LAEs by isolating through the additional information provided by near-passing QSO lines of sight the IGM contribution to the damping of the line. To achieve this goal we will rely on a set of high-resolution, radiative transfer cosmological simulations.
In Section 2, we describe how we fine tune the H density field before running Ly radiative transfer. First we obtain the H density field in a volume through UV radiative transfer simulations at redshift . We normalize the UV photoionization flux so that the synthetic absorption spectra match the observed Ly transmissivity deduced from QSO targets. Next, in Section 3, we present our LAE model study, a numerical aspect of LAE for an individual halo. We present the transmissivities for both LAE and QSO and find a relation between them at different redshifts (Section 4). We conclude with a critical discussion of the results in Section 5. In appendix, we compare two Ly radiative transfer methods, i.e. a full radiative transfer and a simplified model often used in the literature.
2 Neutral hydrogen density field
The first step of the method consists in building an accurate description of the IGM H density field at the epoch at which high-redshift LAE are observed. This is a crucial step as, due to its resonant radiation transfer of the Ly line is very sensitive to the precise value of . The H distribution is obtained through cosmological simulations described in the following; however, it is of outmost importance that the simulated results are calibrated against available absorption line observations before they can be used for our purposes. In the remainder of this Section, we describe how we take care of this aspect.
We start from two (modified, for reasons explained below) snapshots (at ) of the S1 cosmological reionization simulation presented in Baek et al. (2010) from which we take all relevant gas physical properties. Such simulation results have been obtained by post-processing the output of a GADGET-2 (Springel, 2005) cosmological hydrodynamical SPH simulation with the radiative transfer code LICORICE on an adaptive grid. Such grid has been constructed in such a way that each cell contains at most SPH particles. This value gives a minimum cell size of 66.88 comoving kpc at z6, which is smaller than the size of a typical, virialized halo at z6 ( ckpc). The S1 run simulates a volume with baryon and dark matter particles, which gives a mass resolution of for dark matter and for baryons. We assumes WMAP3 cosmological parameters: , , (Spergel et al., 2007). Reionization starts at and ends at as shown in Fig.2 in Baek et al. (2010). Dark matter halos are identified with a friends-of-friends (FOF) algorithm. For each halo, we obtain its dark matter, , and gas, , mass and star formation rate (SFR) , with the method described in Baek et al. 2010. All details concerning the source modeling (SFR, IMF, SED, photoionization rate) are described in Baek et al. 2010. We keep track of stars formed inside each SPH particle according to a Schmidt-law, so to at least qualitatively account for photons produced by sub-resolution objects. If the star fraction of a particle is 0.1% (corresponding to a mass of ), we considered it as a UV source particle. The number of source particles outside resolved halos is about 70% in number, corresponding to 20% of the total luminosity of ionizing sources. At , the simulated mean star formation rate density is , in good agreement with current data (Hopkins & Beacom, 2006; Bouwens et al., 2011). Fig. 1 shows the SFR as a function of the halo mass of the galaxies identified as LAEs at redshift . A tight correlation between the two quantities is clearly observed, with a dependence of the SFR on halo mass approximately as , with in the range .
It is assumed that stars form according to a Salpeter Initial Mass Function (IMF) in the range . The galaxy Spectral Energy Distribution (SED) has been calculated by integrating the stellar properties given in Meynet & Maeder (2005) and Hansen & Kawaler (1994) over the IMF. We use 1000 frequency bins between 13.6 eV and 100 eV for stellar type UV sources, and a photon packet propagates under periodic boundary condition. We take into account supernova feedback by simply injecting in the surrounding gas particles an amount of erg per unit solar mass of stars formed of which 20% (80%) is in thermal (kinetic) form. All SPH neighbor particles are kicked with a velocity depending on the distance from the center. As we will explain later on, matching the observed QSO transmissivity data, requires an ionizing photon escape fraction .
Even a small IGM H neutral fraction, is sufficient to produce a Ly scattering optical depth
larger than unity for mean IGM conditions at the redshift of main interest in this paper, . Indeed the simulated volume-weighted H neutral fraction is found to be in the range . Hence, due to the sensitivity of on it is necessary to determine the latter quantity with a precision . The LICORICE code that we use for photo-ionizing radiative transfer uses a ray-tracing Monte Carlo method on adaptive grid. The snapshot extracted from S1 shows Monte Carlo noise fluctuations on which exceed the required precision. To overcome this problem and also to normalize the UV flux with QSOs observations by varying , we post-process again the S1 snapshot. With an increased (by a factor ) number of photon packets, corresponding to photon packets during a 10 Myr evolution time; on average photon packets pass through each radiative transfer cell. This procedure enables us to achieve a precision on and hence to meet the required standard as shown in the next paragraph.
The two panels in Fig. 2 show the probability distribution function (PDF) of and photoionization rate, , in the high precision run. The PDF has a peak at and closely approximates a log-normal distribution. The volume-weighted mean is . At this redshift, the reionization is already completed except for very high density clumps (visible as a high tail of the PDF) where the recombination rate is boosted. From the comparison of the PDFs in the right panel of Fig. 2 we can also conclude that the IGM is in photoionization equilibrium. The blue curve, , represents the photoionization rate distribution for the particles and it has been obtained from the full radiative transfer simulation. The red curve , is instead the photoionization rate reconstructed from the simulated values of , gas temperature, , and density, , assuming ionization-recombination equilibrium, so that
where is the temperature-dependent recombination rate. The volume-weighted photoionization rate is . Both values and are consistent with the work of Mitra et al. 2011b which synthesizes and analyzes the observational constraints on reionization.
Fig. 3 and Fig. 4 show the distributions at , the latter obtained by rerunning the snapshot of S1 at with escape fraction for comparison. At , the simulation show that reionization is essentially complete and the neutral density field is relatively homogeneous. At , instead, large portions (25%) of the simulated volume are still essentially neutral and the ionization field is patchy with the HII regions confined around sites of intense galaxy/star formation (i.e., inside-out topology).
We caution that the mean we obtain at might be somewhat high. This is likely caused by the fact that in the original S1 simulation, appreciable star formation appears only around (set by the simulation resolution). The resulting integrated e.s. optical depth is 0.062, to be compared with the WMAP7 fiducial range observed value (Jarosik et al., 2011). On the other hand it is useful to remember that the Gunn-Peterson effect in high- quasars only gives lower limits for . The above tension does not represent a major issue for our study as in practice it only modifies slightly the relation between redshift and corresponding mean .
The volume-weighted mean we find is in excellent agreement with the data (e.g. Fan et al. 2006b) ; however, there is extra information contained in the PDF distribution (Fig. 2, left panel) than can be used to calibrate the model. To this aim, we use the observed Ly absorption line spectra of 17 QSOs from Songaila (2004) and Fan et al. (2006b). Let us define the transmission at a given redshift, as the average ratio of observed flux to the intrinsic (unabsorbed) one:
where the average is over a certain wavelength interval along the line of sight (LOS).
In order to compute from our hydro simulations, we interpolate the H density field on a fixed grid, using the SPH smoothing length of each particle and compute the optical depth . We produce 3000 synthetic spectra by piercing through the simulation box with randomly oriented LOS, using the same spectral resolution () and wavelength interval (15Å) to compute Eq. (3) as assumed in the data (Songaila 2004). The best match to the data is obtained by iteration, recursively varying the assumed value.
In Fig. 5, we compare the LOS-averaged simulated with that obtained from 17 observed quasar spectra. Since the emitted photons are all on the blue side of the Ly resonance, their transmissivity is very sensitive to the H distribution. The form of histogram varies sensitively with escape fraction; we found that the best fit with using our source model (initial mass function, spectral energy distribution, stellar lifetime etc). We verified that for synthetic spectra converges to the final distribution to %. This comparison assures that our modeling of the IGM H density field is accurate enough to robustly predict the Ly resonant line transfer from sources in the simulation volume.
3 Simulating LAEs
Ly photons experience both frequency and spatial diffusion during the propagation and modeling the Ly spectrum emerging from galaxies requires a high level of accuracy due to the resonant nature of the Ly line. With increasing computational power, several groups have developed Ly transfer codes (e.g.Zheng & Miralda-Escudé 2002; Ahn et al. 2002; Dijkstra et al. 2006; Verhamme et al. 2006; Tasitsiomi 2006; Semelin et al. 2007; Pierleoni et al. 2009; Laursen et al. 2009; Faucher-Giguère et al. 2010). Most codes are based on Monte Carlo (MC) methods, capable of treating arbitrary geometries and various scattering processes. In this Section we describe our LAE model with radiative transfer and discuss numerical aspects of the LAE modeling. Ly radiative transfer in this work is done by LICORICE, which has been widely tested and used in our previous works (Semelin et al. 2007; Baek et al. 2009, 2010; Vonlanthen et al. 2011). However, differently from the ionizing UV radiative transfer, the Lyman line sector of LICORICE uses a fixed grid.
3.1 LAE physical model
We assume that the main source of Ly emission in LAEs is star formation, and compute the intrinsic luminosity of Ly from their SFR (see Sec. 2). About of ionizing photons are converted to Ly photons assuming case-B recombination (Osterbrock, 1989), so the intrinsic Ly luminosity can be expressed as,
where is the ionizing photon production rate. The intrinsic luminosity can be attenuated by dust absorption of Ly photons in the LAE interstellar medium (ISM). As a result, only a fraction can escape into the IGM. Secondly, the frequency distribution of Ly photons is broadened by Doppler effects induced by galaxy rotation. As a consequence, the spectrum emerging from the LAE has the form
where and is the circular rotation velocity of the LAE. If star formation is a quiescent process, it is for realistic halo and disc properties (Mo et al. 1998; Santos 2004), where is the halo rotational velocity at radius within which the collapsed region has a mean over density of 200 times the background. We use in this work. Finally, due to damping by H in the IGM, only a fraction, of photons escaping out of the LAE actually reaches us. Hence, the observed bolometric Ly luminosity is,
In the following we will fix , i.e. we neglect internal dust absorption. This is because our aim is to isolate the effects of IGM on the Ly transmissivity. A detailed study of the effects of grains can be found in Dayal et al. (2011). As already mentioned, they found a dust absorption/IGM transmissivity degeneracy that we aim at breaking with this work through an independent determination of . Before we proceed, we pause to discuss some interesting side issues.
3.2 Spectral imaging
In order to produce images and spectra of LAEs, we store the data as a three-dimensional array similarly to Zheng & Miralda-Escudé (2002). Two dimensions correspond to the sky plane (thus giving imaging information), while the third one represents the frequency axis (spectral information) like in an IFU data cube. From the beginning and at each scattering, we compute the probability, , for the photon to escape along the direction perpendicular to the sky plane; the photon is then added to the element of the array corresponding to the projected position and frequency. This procedure is necessary because only a very small fraction of the Monte Carlo photons escape exactly in the direction of the image (i.e. observer). The direction into which the photon is scattered is deduced from a dipolar phase function , where is the angle between the incident and outgoing directions. The probability is added to the array with weight , where is the solid angle subtended by the image pixel size, and is the cosine of the angle between the incident photon and the direction perpendicular to the image plane. The optical depth is computed from the gas density along the LOS to the image, with the frequency that the photon would have if it had been scattered in that direction. The surface brightness of each pixel of the constructed image is
7 where is the bolometric luminosity, the luminosity distance to the galaxy, the number of photon packets, the unit vector of the incident photon and is the one to the observer/image. The sum is performed over all photons and all of their scatterings. LICORICE uses the comoving frequency and varying expansion factor as a function of the propagation time of the photon (Baek et al., 2009). This method is more accurate than the usual approach which uses the same expansion factor in the whole cosmological simulation box and avoids first order errors in , which are not negligible when the propagation time reaches several tens of Myr. To determine frequency diffusion after each scattering, we take into account peculiar velocities; the Hubble flow is also included as explained above.
Galaxies with and observed equivalent width Å are identified as LAEs. For each identified LAE in the simulation box, we select a smaller volume of centered on it. Such volume size is chosen in such a way that all Ly photons on the blue side of the line are redshifted and hence not scattered anymore; we have checked that on sub-volumes is not seen to vary any longer. We then interpolate all physical properties of the gas particles in a grid of weighting on the kernel, , i.e. a spherically symmetric spline function which depends on the smoothing length of SPH simulation (see Monaghan & Lattanzio 1985). A resolved halo contains several tens of SPH particles within 10-100 kpc, so we can obtain higher resolution by interpolating on a finer grid as shown in Fig. 6. The SPH density and velocity distribution can be sensibly improved by interpolating on progressively finer grid using SPH smoothing length. We check the total transmitted flux varying grid size as shown in Fig. 6. Above grid linear size , begins to converge. The pixel size is 16 kpc (physical) or at ; the probabilities are collected in a spectral array of 240 bins spaced by 0.1 Å in wavelength, thus yielding a total spectral image depth of 24 Å (rest frame).
3.3 A prototypical LAE
To illustrate the physical features of a prototypical LAEat we select a halo with , intrinsic luminosity , and rotation velocity . These properties are consistent with those typically inferred for LAEs. The H column density integrated over the line-of-sight through the center of the simulation box () is . Assuming a static gas with homogeneous temperature K, this value gives a Ly optical depth at line center .
Fig. 7 shows an overview of the results. Panel (a) is the Ly surface brightness map obtained by collapsing data along the spectral array dimension. In practice, the current detection threshold of LAE narrowband imaging is (e.g. Shimasaku et al. 2006; Ouchi et al. 2008; Zheng et al. 2010). Only the central pixel () exceeds such detection limit. This is because we are assuming a point source located at the galaxy center and therefore the photons reaching us are predominantly those escaping directly with little contribution from scattered radiation. At larger distances instead, all the flux is produced by Ly scattering by H in the halo; its surface brightness dims with radius due to the progressive decrease of H density. The halo SB is potentially at reach of future instruments and can be used to gain unique information about the IGM density structure around the emitter.
The spectra of the source at three different spatial locations are also shown in Fig. 7. The black dotted curve in each sub-panel shows the intrinsic Ly profile emerging from the galaxy before being filtered by the IGM; the solid curves refer to the transmitted spectrum. The spectrum of the central pixel is very asymmetric as photons in the blue side of Ly are completely suppressed as they are redshifted through the line core, while photons on the red side are fully transmitted beyond Å.
Let us now analyze the spectra in two external regions, kpc, and kpc. The SB observed at those locations is due to scattered photons escaping from the central pixel. Moving outward from the source, both and the magnitude of the peculiar velocity decrease with consequent reduction of the frequency diffusion; in addition, the damping of the line core becomes less prominent. Both spectra show two sharp peaks in analogy with the case of an expanding sphere (Semelin et al., 2007); the left peak (blue side of Ly line) is more suppressed than the right one because of Hubble expansion, which starts to dominate over the infall velocity for distances kpc. The separation between the two peaks (Å) is mainly set by the gas velocity field rather than by the H density since is relatively small. In fact, with K and , the analytic solution of Neufeld (1990), holding for a monochromatic source and an extremely optically thick system222See Dijkstra et al. (2006) for the static, homogeneous sphere case., gives a separation of Å at redshift , i.e. smaller than what is obtained here.
Panel (c) of Fig. 7 shows the spherically averaged Ly SB profile as a function of galactocentric radius. The SB profile steeply () drops until kpc, from where it flattens into a plateau extending to kpc, followed by a final descent with slope similar to the inner one. The SB dimming is clearly related to the decrease of the H column density with radius. The presence of spherically symmetric Ly halos around high redshift sources was pointed out by Loeb & Rybicki (1999), who predicted the existence of a roughly uniform surface brightness region of size kpc for the cosmology adopted here. Our simulation, allowing a more precise description of the H density distribution, is in rough agreement with the above conclusion if the uniform surface brightness is identified with the plateau extension; the SB profile is however not constant within the Ly halo as it steepens towards the center.
The above agreement might be somewhat coincidental for the following reasons. In our simulation, the radius at which the Ly optical depth approaches unity, hence allowing photons to free-stream away from the source, is kpc, a factor 2 smaller than predicted by the analytical model of Loeb & Rybicki (1999), which was obtained by assuming a completely neutral, average density, homogeneous gas with no peculiar velocity added to the Hubble flow. In contrast, our results show that the typical in the halo is and the mean gas over density is . If these values are used, becomes much smaller than the analytical prediction. However, the further inclusion of peculiar velocities, restores the rough agreement with the Loeb & Rybicki (1999) model.
The previous argument demonstrates that infalling gas is the dominant process causing a reduction of transmissivity: when infall is artificially turned off, in fact, transmissivity for the selected LAE increases from 0.25 to 0.35. In principle, in addition to infall, additional peculiar velocities might be induced by outflows, which can modify the visibility of LAEs. For example, Dijkstra & Wyithe (2010) and Dijkstra et al. (2011), by using the spherically symmetric thin shell model introduced by Verhamme et al. (2006, 2008), found that the transmissivity can be even through a fully neutral IGM. Although outflows are included in our simulations, their effect on visibility appears to be rather minor. The only visible effect of the galactic outflow emerging from our prototypical LAE is to decelerate the infalling gas, as it can be seen from the change in the velocity profile slope within the central 50 kpc (see Fig. 8). Stated differently, it appears that the outflow expansion is severely quenched by the ram pressure of the infalling gas, thus affecting only mildly the final transmissivity. A simple estimate reinforces this conclusion. The radius at which the infalling gas kinetic input rate becomes equal to the mechanical luminosity, , of the outflow can be computed from the following equality
In the previous equation is the infalling gas over density with respect to the background density ; we have assumed 1 SN, of total energy erg, every 100 solar masses of stars formed and a kinetic energy conversion efficiency . It is clear that the infall can quench the outflow on very small scales, and therefore its influence on the escape of Ly photons is marginal. Very likely though the actual situation could be far more complex than it is possible to mode here: hydrodynamical instabilities might perturb the wind/infall interface; the wind might be asymmetric and or bipolar. These complications require a dedicated investigation that is beyond the scope of this study and we leave them for future work.
The results shown in this Section are based on full radiative transfer of Ly photons. However, this procedure results in relatively large shot noise due to Monte Carlo sampling, making the comparison with the QSO transmissivity (see below) challenging. For this reason in the following we will use the model, most often used in literature, to compute rather than the full radiative transfer. We have performed a careful comparison (discussed in the Appendix) between the two methods, from which we conclude that the differences are not significant. Therefore we use the exponential model as the fiducial case in the rest of the paper.
4 Lyman Alpha Transmissivity
The main aim of this paper is to put joint constraint on the LAE Ly transmissivity using information coming both from the LAE itself (i.e. in emission) and from the absorption spectra of background quasars. To perform this task, we cast thousands of LOS around each LAE halo to produce synthetic absorption QSO spectra with varying impact parameter. By analyzing these two sets of data in combination, we show in the following that it is possible to gain key information on the environment in which LAEs are embedded and to break the dust absorption/HI scattering degeneracy that plagues the use of LAEs as reionization probes. We discuss the properties of the Ly transmissivity in LAEs and QSO separately first; next we combine the two result sets.
4.1 Lyman Alpha Emitters
Using the friends-of-friends (FOF) algorithm, we identified 3070 halos at redshift , corresponding to a halo number density . For each halo, as explained above, we select a sub-volume of centered on it and interpolate all physical properties in a grid of cells. We use the model (see Appendix A) and we average over 3 LOS along the axis from the halo center to the edge of simulation box. We identified objects meeting the LAE selection criterion ( and EWÅ).
|Model||RT method||PV||Ly profile|
Fig. 9 shows as a function of the identified LAE host halo mass . The various radiative transfer parameters for each curve/model are described in Tab. 1. depends on several physical conditions, such as the LAE luminosity, clustering, Ly profile width, gas infall, surrounding . The first three properties are positively correlated with while the remaining ones tend to suppress ; the net effect therefore involves a complex interplay among them. We run several simulations to understand how these elements affect . L3, implementing the model, is our fiducial run; L2 uses instead full radiative transfer. As discussed in Appendix A, transmissivities of L2 are slightly larger than L3 but show similar pattern.
Except for the first three low-mass bins, , almost independent on halo mass. This is an indication of a balance between transmissivity-enhancing and suppressing factors. For example, massive galaxies have larger SFR (see Fig. 1) and consequently higher photoionization rates. As a result, their H regions are larger and more transparent to Ly photons (larger ). However, massive galaxies have also more prominent infall rates which tend to suppress transmissivity. This physical interpretation is supported by the results from run L1, in which we turned off peculiar velocity effects. Consistently with the above, increases with halo mass. However, we find that the positive correlation between and is not as much driven by the increased photoionization rate, but rather by the broader intrinsic line widths. In fact, at where the neutral HI density is low and very homogeneous, the ionizing flux from the LAE itself does not alter the ionization state of the IGM, and therefore the dependence on the halo mass is cancelled. Instead, if we run the L4 model without peculiar velocities and constant Ly line width, we recover again a flat slope. This is because broader Ly line widths allow the escape of more red-wing photons. We conclude that remains approximately constant with halo mass because peculiar velocities and Ly line width effects roughly balance each other.
For comparison, we run L5 at with same additional parameters as L3. At this redshift, reionization is still relatively patchy and the effects of the HII regions carved by the LAEs themselves becomes dominant. Hence, more massive galaxies are embedded in larger HII bubbles allowing a higher fraction of Ly photons to reach us, i.e. increases with halo mass.
The next piece of information that we need to compute is the transmissivity, , of QSOs whose LOS passes through the environment of a foreground LAE at a given impact parameter, (see Fig. 10). is computed following eq. 3, but we raised the spectral resolution to to capture properties of the IGM around a LAE as much as possible. At one spectral bin size corresponds to Å in the rest frame; for comparison, at , 1 comoving Mpc corresponds to Å (rest frame) i.e. 17 spectral bins. We select the transmitted flux of 17 spectral bins in the range where is the redshift of the LAE. This range corresponds to positions in cMpc around the LAE center along the LOS. For each LAE, we cast 12,000 LOS along varying , and compute averaged values of in the following four impact parameter ranges: .
Differently from , at shows a strong dependence on (Fig.11). This effect is caused by the gas infall around these halos, whose velocity increases rapidly with halo mass, and particularly beyond . In addition, for massive halos () is less dependent on impact parameter than for smaller halos.
At , however, transmissivity drops sharply to below , since reionization is not yet completed () and many smaller halos are still embedded in the dense neutral IGM.
As we have emphasized several times, gas peculiar velocities are an important physical factor in the determination of . We define , the modulus of the infall velocity, as
where is the number of LOS, is the number of pixels around a LAE, is the peculiar velocity of the pixel and is velocity of the central pixel which contains the halo to the direction of the LOS. We take along the three axis and which corresponds to Mpc (comoving) around the halo, i.e. a few times the virial radius of typical LAEs, where gas infall should become evident.
Fig. 12 shows infall velocities of LAEs as a function of halo mass at redshift 5.7 and 6.6. increases with halo mass at both redshifts. Indeed, as we will see shortly, is a key parameter to understand as well as .
Photons from the QSO are in the blue side of Ly line and they are completely scattered out of the LOS by even tiny H amounts. Usually, the around a LAE is large enough to strongly suppress the quasar flux; however, peculiar velocities (infall and outflows) can shift the core of the Voigt profile and prevent the suppression. Different from , a large infall velocity increases since the core of the Voigt profile is shifted from the central Ly frequency. As we show in Fig. 10, photons in the red (blue) side of interact with the gas having a velocity component parallel (antiparallel) to their propagation direction. In both cases, the shift of the Voigt profile core produces a increase.
Panel (a) of Fig. 13 shows all data points for at different impact parameters , while panel (b) shows the median of the distribution, both as a function of . Values are not presented since they overshoot the observational precision (). In this case a clear correlation between and exists. The slope is sharpest at the smallest impact parameter and at large infall velocities. This is expected from the simple argument that the environmental properties of the gas deeper into the LAE gravitational potential, where the component of the peculiar velocity along the LOS and are higher (see Fig.6), deviate more from the mean IGM ones. At larger impact parameters the slope of the median becomes flatter and closer to the mean IGM value that we obtained in the Sec. 2. The medians become independent of beyond 200 , exceeding the mean IGM transmissivity by a factor of about 3.
For QSOs, the peculiar velocity is the most dominant factor to determine their transmissivity. is computed in the wavelength range Å (rest frame, covers the 1 cMpc bracketing the center) corresponding to a Doppler shift velocity of 135 . The magnitude of for massive halo is about 150-200 , well above the velocity required to shift the core of the Voigt profile out of the above wavelength interval. Thus even if the is highest within the smaller impact parameter ( cMpc), the infall velocity is large enough to produce an increase of . Since the infall velocity is correlated with halo mass, we observe a similar tendency - the independence of on impact parameter - in panels (b) of Fig. 11 and Fig. 13.
Panel (c) of Fig. 13 shows similar data for for the closest impact parameter, which shows a steeper positive correlation between and . At this redshift, 25% of the volume is still neutral and some of the halos are embedded into completely neutral patches. As mentioned, at higher redshift the global ionization state of the gas dominates the transmissivity over peculiar velocity effects and the size of H bubble plays a more important role on than . The data of at is more scattered than the data at due to the inhomogeneous distribution and the smaller number of resolved halos. The median of below drops very sharply and even the median is under for . We interpret this trend as due to the fact that halos with small host less massive/luminous galaxies, possibly embedded in smaller HII regions leading to a complete suppression.
4.3 Lyman Alpha Emitters/Quasars correlation
In the previous Sections, we have separately analyzed the properties of and depending on the relevant physical parameters. We now ask the question: is there a relationship between the two transmissivities ? To clarify this point, we plot the - correlation in Fig. 14. A general property of the two transmissivities is that is sensitive to various physical parameters (e.g.infall velocity, clustering, neutral fraction) while is not. This is because the transmitted photons from QSOs are all in the blue side of the Ly line, hence they experience strong scattering in the Ly profile core as they are redshifted by cosmic expansion. On the other hand, almost all transmitted photons from LAE are already in the red side of the line having small cross section to H atoms when they escape from the galaxy. Only photons near to the core of the Ly line are affected by infall or H fraction.
At , we do not find a strong correlation between the two transmissivities. Most of the data points loiter in 0.2-0.3, whereas shows a large scatter in . is weakly anti-correlated with and this erases the possible correlation with , which instead is strongly correlated with the infall velocity. However, at there is a clear positive correlation between the two transmissivities, although the dispersion around it is relatively large. as shown in panel (b) of Fig. 14. Almost data below drop under which is the detection criteria, but they increase very sharply from to between and . LAEs with low are embedded in the almost neutral IGM patches, and the photons from background QSO passing in their vicinity have a very high probability to be scattered out of the LOS.
5 Summary and discussion
The evolution of the luminosity function (LF) of high- galaxies and, in particular, of LAEs is one of the most promising tools to study the interplay between early galaxy formation and cosmic reionization. However, before solid conclusions can be drawn from such a method, a number of delicate aspects must be worked out quantitatively. One of the most challenging ones is the presence of a degeneracy between Ly line damping by dust in the source and/or by intergalactic H . As Dayal et al. (2011) have pointed out, a wide range of H fractions is compatible with the observed Ly luminosity function since the effects of a largely neutral IGM are compensated by a higher escape fraction of Ly photons from galaxies containing only small amounts of dust. With the aim of breaking this dust/H degeneracy we have proposed a possible solution which relies on an independent determination of exploiting the presence of a background luminous source (typically a QSO but the use of Gamma Ray Bursts can be also conceived) whose LOS pierces through and samples the target LAE surrounding matter. The theoretical basis of the method, which we have applied at two different redshifts and 6.6 where large LAE and QSO samples are available, can be summarized as follows.
First we normalize the UV ionizing flux and the H field by using results from QSO absorption line experiments. We produce 3000 synthetic QSO spectra at redshift with randomly selected position and direction in a simulation volume. Next we obtain the transmissivity for each LOS, using the same spectral resolution and wavelength range as for a sample of 17 observed QSOs (Songaila, 2004; Fan et al., 2006a). We find that the transmissivity distribution changes very sensitively with the mean UV photoionization rate; the best fit to the data is obtained for a volume-averaged .
Before running Ly radiative transfer for all identified LAEs in the simulation, we have presented the individual case of a prototypical LAE to illustrate some relevant physical features behind the Ly line transfer. The observed surface brightness () is the sum of the central point source and a scattering halo extending up to 150 kpc from the galaxy center. Such halo is just below current detection threshold and at reach of future experiments. Our results also indicate that Ly photons in the red wing of the line are largely transmitted at this redshift. We also find that infalling gas is the dominant process causing a reduction of transmissivity, whereas outflows seem to be efficiently quenched by infall ram pressure, thus confining their effects in a relatively small region around the galaxy ( kpc).
For each resolved LAE in the simulation, we select a small box of around it and derive the optical depth along the three LOS parallel to the three viewing axes. At redshift , we find that , almost independent of the halo mass. This constancy arises from the conspiracy of two different physical effects: (i) the intrinsic Ly line width and (ii) the infall peculiar velocity. The Ly width depends on the galaxy rotation velocity, which increases with halo mass; the broader is the line, the larger is the fraction of red wing photons that can be transmitted, thus increasing . The infall velocity also increases with halo mass, but it acts to suppress , as it blue-shifts the line in the rest frame of the infalling H atoms. Evidently, the two effects balance each other almost perfectly, yielding a roughly constant . At higher redshift, , where the transmissivity is instead largely set by the local H abundance and by the ability of the galaxy to carve a sufficiently large H region around itself. As the SFR and the photoionization rate increases toward large masses, consequently increases with halo mass from 0.15 to 0.3.
As a next step, we cast thousands of LOS originating for background QSOs passing through foreground LAEs at different impact parameters. Differently from the case of LAE, photons from the quasars are in the blue side of the Ly line and are scattered out of the LOS by H in the vicinity of the LAE. Again, we have emphasized the importance of gas infall motions which increase the quasar transmissivity ( ). At smaller impact parameters, cMpc, a positive correlation between and is found at , which tends to become less pronounced (i.e. flatter) at larger distances. Quantitatively, a roughly 10 increase (from to ) of is observed in the range . The correlation is even stronger at .
At , we do find a relatively weak correlation between and . This is because most of the data points loiter in 0.2-0.3, whereas shows a large scatter in . However, at there is a clear positive correlation between the two transmissivities, although the dispersion around it is relatively large. The median is below for , and it increases very sharply from to between and , showing a very strong sensitivity of to . The correlation signal is diluted by the fact that the properties of the environment of many LAEs (i.e the most clustered ones) is not completely determined by the radiation and peculiar motions caused by the central object alone. Thus a QSO line of sight through such a clustered environment results in contaminating high data. Properly selecting relatively isolated LAEs would boost the correlation signal.
The proposed method therefore appears to be promising to determine the evolution of the physical state of the IGM with redshift and to allow a detailed study of reionization. Obviously, its feasibility relies on the availability of a sufficient number of high redshift LAEs and QSOs. In this sense, we live in an age full of promises. Currently, about 60 QSOs at 5.7 6.5 have been detected in total. Of these, Fan et al. (2006a) discovered 19 QSOs with and Songaila (2004) discovered 5 QSOs at in the Sloan Digital Sky Survey333http://www.sdss.org/ (SDSS) using moderate resolution spectra, 5000. Additional 19 QSOs are discovered by Willott et al. (2007, 2009, 2010) in the Canada-France High- Quasar Survey (CFHQS) with moderate resolution . In the VLA FIRST survey444http://sundog.stsci.edu, Becker et al. (2006, 2010) have obtained 17 QSOs with high () and moderate (5000) spectral resolution. Recently Mortlock et al. (2011) discovered the highest QSO at . New QSOs are still being discovered and identifying QSOs at even higher redshifts is the next challenge. Ongoing near-IR sky surveys like the UKIDSS (Lawrence et al. 2007) will be completed in few years. The Large Area Survey (4000 ) aims at finding quasars. Finally, the ESO VISTA Telescope surveys (Sutherland, 2009) will be even more efficient building on a larger camera probably capable of finding sources.
On the LAE side, a large sample already exists at . Up to now 54 (45) LAEs are confirmed by spectroscopic observation at redshift (6.5) in the Subaru Deep Field (Kashikawa et al., 2011). The number of confirmed LAEs are 70% (81%) of their photometric candidates at (6.5). The Subaru/XMM-Newton Deep Survey (SXDS) has also found a large number of LAEs at . In a sky Ouchi et al. (2008) found 401 photometric candidates, 17 of which spectroscopically confirmed. The survey locations of SDF (, ) and SXDS (, ) are within the survey area of VLA FIRST and SDSS which provides a first direct opportunity to test our method.
The survey field for LAE is about 1 in which hundreds of candidates and tens of confirmed LAE have been found. On the other hand QSOs are found over wide area over thousands of and the probability that we can find them next to the target LAE is low. A much more efficient strategy to find QSO-LAE pairs would be to perform a deep spectroscopic sky survey around a known high- QSO. The quasar J1335+3533 can be an optimal candidate for that, as it is suitably located close to the SDF. We will consider these follow-up ideas in future work.
Appendix A Full radiative transfer vs. exponential models
Zheng et al. (2010) have studied LAEs using full Ly radiative transfer in cosmological reionization simulations. They compare their results with previous works (e.g. McQuinn et al. 2007; Iliev et al. 2008; Dayal et al. 2011), where the transfer of Ly is treated in an approximate manner: in such works, in fact, the intrinsic line profile is damped by , where is the scattering optical depth at frequency along the LOS. Zheng et al. (2010) emphasize that the model does not provide an accurate description of the observed Ly spectrum, because it does not account for frequency (no frequency change occurs for any Ly photon) and space (no surface brightness information) diffusion. To investigate further this point, in addition to the full radiative transfer model discussed so far, we have reconstructed the LAE spectra also with the exponential model. We can then evaluate the differences between the two.
The black curve in Fig.15 shows the spectrum from the central pixel of the prototypical LAE at presented in Sec. 3.3. The apparent shift of the peak with respect to the intrinsic profile is about 5 Å. This is consistent with the result of Zheng et al. (2010) (see their Fig.6) for the same halo mass, . The red curve in the same Figure has been obtained from the model.
We find that at all frequencies such model underestimates the flux, but different from Zheng et al. (2010), we find the difference between the two methods on to be %. Most of the scatterings occurring in the central pixel shift the photon frequency in the range within , corresponding to the line core and thus they are efficiently damped. Average over all the LAEs identified in the simulation, the transmissivity computed with full radiative transfer is larger by about ; some halos (0.5% of total) are embedded in a dense HI region, where . Full radiative transfer of these halos gives , while the model suppresses the line almost completely.
We have analyzed the possible reasons for the discrepancy between Zheng et al. (2010) and our findings. We came to the conclusion that it originates from the combination of two factors: the different line width and IGM properties. Zheng et al. (2010) used line widths
This is about 3 times smaller than what adopted here. With this line width, almost all the photons are initially within the core of the line profile and the flux of the model is greatly suppressed; as a result only a very small amount of photons can escape without any scattering and the main contribution to the observed flux is coming from scattered photons. In our case, on the other hand, most photons with Å are transmitted without any scattering and only a very small contribution is coming from scattered photons. Even when we enlarged the integrating flux area to kpc the difference with respect to the model remains %. The very small contribution from scattered photons is related to a different HI density field. As discussed above we have carefully calibrated the IGM transmissivity through an iterative optimization with QSO absorption spectra. This procedure gives on average a 4 times lower than adopted by Zheng et al. (2010), which results in a reduced spatial/frequency diffusion.
Finally, the neutral density at the very center of halo can be underestimated due to the lack of mass resolution. The cell size of Ly radiative transfer is 1.6 times smaller than Zheng et al. (2010) but the mass resolution is 1700 times larger than theirs. Hence, we have less particles per halo, potentially leading to underestimating the recombination rate and . is determined by the amount of transmitted photons rather than scattered photons. However, we verify that even with 10 times higher H density in the central cell, the photons with Å are transmitted without scattering since they are far from the core of the line profile. The fraction of photons transmitted without scattering events is rather sensitive to the peculiar velocity of the infalling gas. As we shown in Fig.8, however, the infall profile is well captured by remapping SPH particles on a finer grid.
We acknowledge useful discussions with B. Ciardi, A. Mesinger, P. Laursen and other DAVID members.
- Adelberger et al. (2005) Adelberger, K. L., Shapley, A. E., Steidel, C. C., Pettini, M., Erb, D. K., & Reddy, N. A. 2005, ApJ, 629, 636
- Adelberger et al. (2003) Adelberger, K. L., Steidel, C. C., Shapley, A. E., & Pettini, M. 2003, ApJ, 584, 45
- Ahn et al. (2002) Ahn, S.-H., Lee, H.-W., & Lee, H. M. 2002, ApJ, 567, 922
- Baek et al. (2009) Baek, S., di Matteo, P., Semelin, B., Combes, F., & Revaz, Y. 2009, A&A, 495, 389
- Baek et al. (2010) Baek, S., Semelin, B., Di Matteo, P., Revaz, Y., & Combes, F. 2010, A&A, 523, A4+
- Becker et al. (2010) Becker, G. D., Rauch, M., & Sargent, W. L. W. 2010, ApJ, 698, 1010
- Becker et al. (2006) Becker, G. D., Sargent, W. L. W., Rauch, M., & Simcoe, R. A. 2006, ApJ, 640, 69
- Bolton & Haehnelt (2007) Bolton, J. S., & Haehnelt, M. G. 2007, MNRAS, 382, 325
- Bouwens et al. (2011) Bouwens, R. J., et al. 2011, Nature, 469, 504
- Castellano et al. (2010) Castellano, M., et al. 2010, A&A, 524, A28+
- Choudhury & Ferrara (2006) Choudhury, T. R., & Ferrara, A. 2006, MNRAS, 371, L55
- Dayal et al. (2008) Dayal, P., Ferrara, A., & Gallerani, S. 2008, MNRAS, 389, 1683
- Dayal et al. (2009) Dayal, P., Ferrara, A., Saro, A., Salvaterra, R., Borgani, S., & Tornatore, L. 2009, MNRAS, 400, 2000
- Dayal et al. (2011) Dayal, P., Maselli, A., & Ferrara, A. 2011, MNRAS, 410, 830
- Dijkstra et al. (2006) Dijkstra, M., Haiman, Z., & Spaans, M. 2006, ApJ, 649, 14
- Dijkstra et al. (2011) Dijkstra, M., Mesinger, A., & Wyithe, J. S. B. 2011, MNRAS, 414, 2139
- Dijkstra & Wyithe (2010) Dijkstra, M., & Wyithe, J. S. B. 2010, MNRAS, 408, 352
- Dijkstra et al. (2007) Dijkstra, M., Wyithe, J. S. B., & Haiman, Z. 2007, MNRAS, 379, 253
- Fan et al. (2006a) Fan, X., Carilli, C. L., & Keating, B. 2006a, ARA&A, 44, 415
- Fan et al. (2001) Fan, X., et al. 2001, AJ, 122, 2833
- Fan et al. (2006b) —. 2006b, AJ, 132, 117
- Fan et al. (2003) —. 2003, AJ, 125, 1649
- Faucher-Giguère et al. (2010) Faucher-Giguère, C.-A., Kereš, D., Dijkstra, M., Hernquist, L., & Zaldarriaga, M. 2010, ApJ, 725, 633
- Haiman & Cen (2005) Haiman, Z., & Cen, R. 2005, ApJ, 623, 627
- Hansen & Kawaler (1994) Hansen, C. J., & Kawaler, S. D. 1994, Stellar Interiors. Physical Principles, Structure, and Evolution. (Astronomy and Astrophysics Library)
- Hopkins & Beacom (2006) Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142
- Hu et al. (1998) Hu, E. M., Cowie, L. L., & McMahon, R. G. 1998, ApJL, 502, L99+
- Iliev et al. (2009) Iliev, I. T., Pen, U.-L., McDonald, P., Shapiro, P. R., Mellema, G., & Alvarez, M. A. 2009, Ap&SS, 320, 39
- Iliev et al. (2008) Iliev, I. T., Shapiro, P. R., McDonald, P., Mellema, G., & Pen, U.-L. 2008, MNRAS, 391, 63
- Jarosik et al. (2011) Jarosik, N., et al. 2011, ApJS, 192, 14
- Kashikawa et al. (2006) Kashikawa, N., et al. 2006, ApJ, 648, 7
- Kashikawa et al. (2011) —. 2011, ApJ, 734, 119
- Laursen et al. (2009) Laursen, P., Razoumov, A. O., & Sommer-Larsen, J. 2009, ApJ, 696, 853
- Lawrence et al. (2007) Lawrence, A., et al. 2007, MNRAS, 379, 1599
- Loeb & Rybicki (1999) Loeb, A., & Rybicki, G. B. 1999, ApJ, 524, 527
- Malhotra & Rhoads (2004) Malhotra, S., & Rhoads, J. E. 2004, ApJL, 617, L5
- McQuinn et al. (2007) McQuinn, M., Hernquist, L., Zaldarriaga, M., & Dutta, S. 2007, MNRAS, 381, 75
- Meynet & Maeder (2005) Meynet, G., & Maeder, A. 2005, A&A, 429, 581
- Miralda-Escude & Rees (1998) Miralda-Escude, J., & Rees, M. J. 1998, ApJ, 497, 21
- Mitra et al. (2011a) Mitra, S., Choudhury, T. R., & Ferrara, A. 2011a, ArXiv e-prints
- Mitra et al. (2011b) —. 2011b, MNRAS, 413, 1569
- Mo et al. (1998) Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319
- Monaghan & Lattanzio (1985) Monaghan, J. J., & Lattanzio, J. C. 1985, A&A, 149, 135
- Mortlock et al. (2011) Mortlock, D. J., et al. 2011, Nature, 474, 616
- Neufeld (1990) Neufeld, D. A. 1990, ApJ, 350, 216
- Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae and active galactic nuclei (University Science Books, Sausalito, CA)
- Ouchi et al. (2009) Ouchi, M., et al. 2009, ApJ, 696, 1164
- Ouchi et al. (2008) —. 2008, ApJS, 176, 301
- Ouchi et al. (2010) —. 2010, ApJ, 723, 869
- Partridge & Peebles (1967) Partridge, R. B., & Peebles, P. J. E. 1967, ApJ, 147, 868
- Pentericci et al. (2011) Pentericci, L., et al. 2011, ArXiv e-prints
- Pierleoni et al. (2009) Pierleoni, M., Maselli, A., & Ciardi, B. 2009, MNRAS, 393, 872
- Rhoads et al. (2003) Rhoads, J. E., et al. 2003, AJ, 125, 1006
- Santos (2004) Santos, M. R. 2004, MNRAS, 349, 1137
- Semelin et al. (2007) Semelin, B., Combes, F., & Baek, S. 2007, A&A, 474, 365
- Shimasaku et al. (2006) Shimasaku, K., et al. 2006, Pub. Astron. Soc. Japan, 58, 313
- Songaila (2004) Songaila, A. 2004, AJ, 127, 2598
- Spergel et al. (2007) Spergel, D. N., et al. 2007, ApJS, 170, 377
- Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
- Sutherland (2009) Sutherland, W. 2009, in Science with the VLT in the ELT Era, ed. A. Moorwood, 171–+
- Taniguchi et al. (2005) Taniguchi, Y., et al. 2005, Pub. Astron. Soc. Japan, 57, 165
- Tasitsiomi (2006) Tasitsiomi, A. 2006, ApJ, 645, 792
- Todini & Ferrara (2001) Todini, P., & Ferrara, A. 2001, MNRAS, 325, 726
- Verhamme et al. (2008) Verhamme, A., Schaerer, D., Atek, H., & Tapken, C. 2008, A&A, 491, 89
- Verhamme et al. (2006) Verhamme, A., Schaerer, D., & Maselli, A. 2006, A&A, 460, 397
- Vonlanthen et al. (2011) Vonlanthen, P., Semelin, B., Baek, S., & Revaz, Y. 2011, ArXiv e-prints
- Willott et al. (2007) Willott, C. J., et al. 2007, AJ, 134, 2435
- Willott et al. (2009) —. 2009, AJ, 137, 3541
- Willott et al. (2010) —. 2010, AJ, 139, 906
- Zheng et al. (2010) Zheng, Z., Cen, R., Trac, H., & Miralda-Escudé, J. 2010, ApJ, 716, 574
- Zheng & Miralda-Escudé (2002) Zheng, Z., & Miralda-Escudé, J. 2002, ApJ, 578, 33