THE BIVARIATE SIZE–LUMINOSITY RELATIONS FOR LYMAN BREAK GALAXIES AT z ∼ 4–5

LBG color-selection criteria are designed for selecting star-forming galaxies at high redshift using two features in the spectral energy distributions (SEDs) of star-forming galaxies (e.g., Steidel & Hamilton 1992; Steidel et al. 1996; Giavalisco 2002). One feature is the spectral break at rest-frame 912 Å due to the absorption of UV photons by the stellar atmosphere of early-type stars and the neutral hydrogen in the interstellar medium (ISM) surrounding the star-forming regions. The second feature is the blanketing of Lyman-series lines by the intervening neutral hydrogen in the intergalactic medium (IGM) along the lines of sight to galaxies (Madau 1995; Meiksin 2006). Therefore, it is possible to use the colors between a set of broadband filters to target a specific redshift range.

We adopt the color-selection criteria in Giavalisco et al. (2004a), which are optimized to select LBGs at z ∼ 4 and z ∼ 5. The criteria for the B₄₃₅-dropouts are

$$\begin{eqnarray} (B_{435} - V_{606}) & \ge 1.2 + 1.4 \ (V_{606} - z_{850})\ \wedge \nonumber\\ (B_{435} - V_{606}) & \ge 1.2\ \wedge \ (V_{606} - z_{850}) \le 1.2 \end{eqnarray} \tag{ 1 }$$

and the criteria for the V₆₀₆-dropouts are

$$\begin{eqnarray} & \lbrace [(V_{606}-i_{775}) > 1.5 + 0.9 \ (i_{775}-z_{850})]\ \vee \nonumber\\ & [(V_{606}-i_{775}) > 2.0]\rbrace \ \wedge \ (V_{606}-i_{775}) \ge 1.2\ \wedge \nonumber\\ & (i_{775}-z_{850}) \le 1.3, \end{eqnarray} \tag{ 2 }$$

where ∧ is the logical "AND" operator and ∨ is the logical "OR" operator. The colors are computed using the isophotal magnitudes MAG_ISO reported by SExtractor in order to obtain consistent colors. All magnitudes quoted afterward are the MAG_AUTO from SExtractor unless otherwise stated.

We select galaxies that satisfy the above color criteria down to z₈₅₀ = 26.5 mag in GOODS and z₈₅₀ = 28.5 mag in HUDF, while also demanding that the signal-to-noise ratio (S/N) in the z₈₅₀-band be greater than 5. For V₆₀₆-dropouts we also require all sources to have B₄₃₅-band (S/N)_B ⩽ 5, as at z ∼ 5 the B₄₃₅-band is entirely blue-ward of the Lyman break and should not have significant detection. This S/N criterion is more liberal than those in the literature (e.g., Giavalisco et al. 2004a). The reason for adopting a liberal cut here is to increase the completeness of our sample, as our simulations show a non-negligible fraction of the true z ∼ 5 sources whose S/N in B₄₃₅-band are slightly higher than 2 (the value commonly used in the literature) because the SExtractor-estimated S/N becomes unreliable for the non-detected sources. Among the 52 sources with 2 < S/N < 5 in B₄₃₅-band, 6 are broken-up diffraction spikes, 46 have S/N < 3, 8 have 3 < S/N < 4, and only 2 have 4 < S/N < 5. Furthermore, all three of the sources with spectroscopic redshifts that have S/N > 2 in B₄₃₅-band are at z_spec > 4.4 and not low-z interlopers. Because the choice in the S/N limit is photometry-dependent and has some intrinsic uncertainty, we decide to use a more liberal threshold in the B₄₃₅-band S/N to increase our completeness without introducing a large number of interlopers. We also exclude sources with SExtractor-reported CLASS_STAR ⩾0.95 down to z₈₅₀ ⩽26 mag in GOODS and z₈₅₀ ⩽28 mag in the HUDF to eliminate field stars. The resultant sample still contains two active galactic nucleus (AGN) candidates (IAU ID J033216.86−275043.9 and J033243.84−274918.1; Wolf et al. 2004); we therefore exclude these two sources from further analyses. We next visually inspect each selected dropout source to remove other image artifacts that are not removed in the image-reduction process, and this removes 16 and 19 sources from the B₄₃₅-dropout and V₆₀₆-dropout sample, respectively, all of which are broken diffraction spikes from bright stars.

At the bright end of our samples (z₈₅₀ ⩽24.5 mag in GOODS and z₈₅₀ ⩽25.5 mag in the HUDF), the majority of sources have spectroscopic redshifts. For the few sources that do not, all of them have photometric redshifts from Dahlen et al. (2010; we describe the redshift information of our sample in Section 2.2). Therefore, we eliminate the bright sources that are confirmed low-z interlopers (defined as z < 3 for the B₄₃₅-dropout sample and z < 4 for the V₆₀₆-dropout sample) based on available redshift information. Doing so removes 1 B₄₃₅-dropout and 3 V₆₀₆-dropouts from GOODS. In fainter magnitude bins (z₈₅₀ >24.5 mag in GOODS and z₈₅₀ >25.5 mag in the HUDF) we do not have redshift information for most sources, so we will model the interloper fractions in Section 2.3 and include them as an additional component of our model.

In addition, if a given source exists in both the GOODS and HUDF source lists, we use the values measured in the HUDF because it is deeper. With the above criteria we have a total of 1063 B₄₃₅-dropouts and 465 V₆₀₆-dropouts. The number of dropouts in the two GOODS fields (excluding the deeper HUDF), 384 B₄₃₅-dropouts and 196 V₆₀₆-dropouts in GOODS-N and 454 B₄₃₅-dropouts and 157 V₆₀₆-dropouts in GOODS-S, are consistent within the cosmic variance estimated by the cosmic variance calculator of Trenti & Stiavelli (2008), which is about 15%. The numbers are summarized in Table 1.

There have been several campaigns to measure spectroscopic redshifts (z_spec) in the GOODS fields that cover the sources of interest here (Steidel et al. 1999; Cristiani et al. 2000; Cohen et al. 2000; Cohen 2001; Steidel et al. 2003; Wirth et al. 2004; Daddi et al. 2009a, 2009b; Barger et al. 2008; Vanzella et al. 2009; Popesso et al. 2009; D. Stern et al. 2013, in preparation). In addition, Dahlen et al. (2010) produced the deepest photometric-redshift (z_phot) catalog in GOODS-S (the photometric-redshift catalog of Cardamone et al. 2010 in GOODS-S, although derived with more bands, does not go as deep as Dahlen et al. 2010.) In Figure 1, we show the redshift distribution of our selected sources (we only show spectroscopic redshifts with good or acceptable qualities if such flags exists in the catalogs.) In general the color-selection criteria are able to select B₄₃₅-dropouts within the redshift range 3.4 ≲ z ≲ 4.4 and V₆₀₆-dropouts within 4.4 ≲ z ≲ 5.6. We do not use the redshift information to identify interlopers and exclude them from analysis in the faint end because the spectroscopic redshifts are incomplete and the photometric redshifts have larger errors. Instead we try to model the interloper fractions in Section 2.3 and include their contribution as an additional component in our model.

Zoom In Zoom Out Reset image size

High-redshift, star-forming galaxies have a higher fraction of Lyα-emitters (LAEs) than local star-forming galaxies (see, e.g., Reddy et al. 2008; Finkelstein et al. 2009; Cowie et al. 2010; Ciardullo et al. 2012). The fraction of star-forming galaxies at z ≳ 5 with large rest-frame Lyα equivalent widths could also be a function of redshift and luminosity (see, e.g., Stark et al. 2010; Schaerer et al. 2011). At z ∼ 5.6–5.8, Lyα is redshifted into the i₇₇₅-band, therefore a galaxy with strong Lyα emission (rest-frame equivalent width W_Lyα ≳ 50 Å) will have redder V₆₀₆−i₇₇₅ and bluer i₇₇₅-z₈₅₀ colors than the colors determined by the continuum alone, making them more likely to be selected as V₆₀₆-dropouts instead of i₇₇₅-dropouts (see, e.g., Malhotra et al. 2005; Dow-Hygelund et al. 2007; Stanway et al. 2007, 2008; Bouwens et al. 2007; Stark et al. 2010). Among our V₆₀₆-dropouts with spectroscopic redshifts, 6 sources out of 40 (15%) have z_spec ⩾ 5.5, representing a secondary peak of z_spec distribution. It is possible that these sources have strong Lyα-emission-perturbed broadband colors. However, at z ⩾ 5.5 there is yet to be a complete census of Lyα equivalent widths, so we can not model the effects of Lyα emission on our sample selection accurately. A simple estimate of the possible fraction of V₆₀₆-dropouts that are z ∼ 6 galaxies with strong Lyα emission—using the best-fit exponential distribution of Lyα equivalent widths of LAEs at z ∼ 3.1 (Ciardullo et al. 2012) and the number densities of i₇₇₅-dropouts from Bouwens et al. (2007)—suggests that these galaxies make up at most about 5% of our V₆₀₆-dropout sample in GOODS. Therefore, we do not consider z ∼ 6 star-forming galaxies with strong Lyα emission to contribute significantly to the V₆₀₆-dropout sample.

Sources other than high-z star-forming galaxies could also be selected by the LBG criteria. One type of interlopers are sources with intrinsically red SEDs such as dwarf stars (of spectral types M, L, and T) and z ⩽ 2 quiescent galaxies (Ryan et al. 2005; Pirzkal et al. 2005, 2009; Steidel et al. 1999; Yoshida et al. 2006). Another type of interlopers are sources with intrinsic SEDs outside of the target range, but with colors close enough to the selection boundaries that they are scattered into the color-selection window due to photometric error. To minimize the first type of interlopers, we identify low-z interlopers with available redshift information in the bright end in Section 2.1. At the faint end, identifying individual interlopers of either type becomes intractable due to the limited spectroscopic information. In this section, we estimate the interloper contributions from photometric scatter (hereafter "photometric interlopers").

We estimate the photometric interloper fractions at z₈₅₀ >24.5 mag in GOODS and z₈₅₀ >25.5 mag in the HUDF, in magnitude bins of width 0.5 mag, with a set of simulations (the "photometric interloper simulations"). Our simulations are very similar to the ones performed in Bouwens et al. (2006; where they estimated the contamination fraction in the HUDF) and Iwata et al. (2007). In the photometric interloper simulations we use all sources with S/N ⩾ 5 in z₈₅₀-band and S/N ⩾ 3 in B₄₃₅-, V₆₀₆-, and i₇₇₅-bands that are not selected as dropouts. Only sources with either spectroscopic redshift or photometric redshift much lower than z = 3 are used in order to minimize the contribution of sources that are just below the selection redshift range. The sources selected this way have reasonably constrained SED shapes that span a wide range. For each selected source we then rescale its z₈₅₀-band magnitude to a randomly selected value between 24.5 and 26.5 mag (in GOODS) and between 25.5 and 28.5 mag (in HUDF). After rescaling the magnitudes in other bands by the same factor, we randomly perturb its magnitude in each band according to the appropriate photometric error in the given band and magnitude bin (assuming Gaussian error in magnitude). We calculate the number of these sources that fall within the color-selection window in each magnitude bin, and also calculate their fractional contribution among the real dropouts, defined as

$$\begin{equation*} f_{\mathrm{interloper}} = \frac{\mathrm{{\rm number}\ of\ expected\ photometric\ interlopers}}{\mathrm{{\rm number}\ of\ selected\ dropouts}} \end{equation*}$$

The same procedure is repeated 5000 times to collect good statistics.

The results of the photometric interloper simulations are presented in Table 2 and Figure 2. The estimated overall photometric interloper fractions are below 8%, and for the V₆₀₆-dropouts in the HUDF the photometric interloper fraction is especially low, possibly owing to the extremely deep i₇₇₅-band data in HUDF that constrain the i₇₇₅-z₈₅₀ color well. In our photometric-interloper simulations, we do not reject sources using the B₄₃₅-band S/N threshold for V₆₀₆-dropouts, so the true photometric interloper fractions are likely even lower. In general, the photometric-interloper fractions stay relatively constant over all magnitude bins probed.

Zoom In Zoom Out Reset image size

As verification of the photometric interloper fractions from simulations, we compare the estimated photometric interloper fractions f_interloper with those estimated from available spectroscopic redshifts (z_spec) alone. Among the B₄₃₅-dropout sample in GOODS, 90 sources have secure spectroscopic redshifts at z₈₅₀ >24.5 mag, 1 of them have z_spec < 3 (1.1%). Similarly, among the V₆₀₆-dropout sample in GOODS, 37 sources have secure spectroscopic redshifts at z₈₅₀ >24.5 mag, and none of them have z_spec < 4. In the HUDF, only 4 B₄₃₅-dropouts and 4 V₆₀₆-dropouts with z₈₅₀ >25.5 mag have z_spec, and none of them are low-z interlopers. All of the interloper fractions estimated from available spectroscopic redshifts are consistent with what we find from the photometric interloper simulations. In all cases, the estimated interloper fractions from the simulations are higher than those estimated from the sources with spectroscopic redshifts. Because this set of simulations accounts for objects scattered into the color-selection window, but does not use other criteria such as stellarity or image flags to reject sources, the interloper fractions estimated are likely upper limits.

In this study we combine the UV luminosity function and size distribution to form the bivariate size–luminosity distribution in rest-frame 1500 Å. We adopt the Schechter function (Schechter 1976) as our model luminosity function, as is customary in most high-z UV luminosity-function studies and was shown to be a decent parameterization (e.g., Steidel et al. 1999; Bouwens et al. 2007; Reddy & Steidel 2009):

$$\begin{equation} \phi (L)dL = \phi ^* \left(\frac{L}{L^*}\right)^{\alpha }\exp \left(-\frac{L}{L^*}\right)\frac{dL}{L^*}, \end{equation} \tag{ 4 }$$

where ϕ(L)dL is the volume number density of galaxies within the luminosity range (L, L + dL), L* is the characteristic luminosity, and α is the power-law slope at the faint end, and ϕ* is the normalization.

For the size distribution of LBGs at a given luminosity, we adopt the log-normal distribution:

$$\begin{equation} p(R_e)dR_e=\frac{1}{\sigma _{\ln R_e}\sqrt{2\pi }} \exp \left[-\frac{\ln ^2(R_e/\bar{R_e})}{2\sigma _{\ln R_e}^2}\right]\frac{dR_e}{R_e}, \end{equation} \tag{ 5 }$$

where p(R_e)dR_e is the probability density that a galaxy at luminosity L has effective half-light radius between (R_e, R_e + dR_e), and $\bar{R_e}$, $\sigma _{\ln \mathrm{R_e}}$ are the peak and width of the distribution, respectively. The log-normal form of size distribution is motivated by the disk-formation theory of Fall & Efstathiou (1980), in which gas approximately conserves its specific angular momentum when it forms in the centers of DM halos, and disk size is proportional to the spin parameter λ and the halo virial radius R_vir: R_d∝λ R_vir. In this picture, the size distribution of the disks should reflect the underlying λ distribution of the halos, which is found to be well approximated by a log-normal distribution independent of redshift or halo mass (Barnes & Efstathiou 1987; Warren et al. 1992; Bullock et al. 2001).

To connect Equations (4) and (5), we follow the arguments of de Jong & Lacey (2000; see also Choloniewski 1985) that adopt a power-law relation to connect the luminosity L and the peak of the size distribution $\bar{R}$:

$$\begin{equation} \bar{R_e}(L) = R_0\left(\frac{L}{L_0}\right)^\beta. \end{equation} \tag{ 6 }$$

The peak of the size distribution $\bar{R}$ at a nominal luminosity L₀ is defined as $R_0 \equiv \bar{R_e}(L_0)$. We choose L₀ to correspond to M₀ = −21, or L_ν = 9.12 × 10⁻¹² erg s⁻¹ Hz⁻¹ at 1500 Å. Using Equation (6) to connect Equations (4) and (5), and expressing the distribution in terms of absolute magnitude M, we come to the full bivariate expression:

$$\begin{eqnarray} && \Psi (M,\log _{10} R_e)\,dM\,d\log _{10} R_e =0.4\ln (10) \nonumber\\ &&\quad \times 10^{-0.4(\alpha +1)(M-M^*)}\exp [-10^{-0.4(M-M^*)}] \nonumber\\ &&\quad \times \frac{\ln (10)}{\sigma _{\ln R_e}\sqrt{2\pi }} \exp \left\lbrace \frac{-[ \log _{10}(R_e/R_0) + 0.4\beta (M-M_0)]^2}{2\sigma _{\ln R_e}^2/\ln ^2(10)} \right\rbrace \nonumber\\ &&\quad \times dM\,d\log _{10} R_e. \end{eqnarray} \tag{ 7 }$$

In Equation (7), the bivariate distribution is characterized by a set of five free parameters $\bm {{P}} \equiv (\alpha , M^*, R_0, \sigma _{\ln R_e}, \beta)$ (again we determine ϕ* after we determine the shape of the distribution). We show in Figure 7 the distribution in the (M, log₁₀R_e) plane corresponding to the fiducial parameters $\bm {{P}}_0 = (-1.7, -21.0, 0\buildrel{\prime \prime}\over{.}21, 0.7, 0.3)$ down to apparent magnitude of 26.5 mag. The regions with brighter shades in Figure 7 denote higher number densities. We choose the values $\bm {{P}}_0$ because these are the initial guesses of the parameters in the fitting process. If the intrinsic distribution of the galaxy population is indeed described by Equation (7) with $\bm {{P}}_0$, one would expect the observed points to follow the top panel of Figure 7 if the sample is complete and the measurement error is zero. However, there are always sample incompleteness and there are always measurement errors. In the following sections we detail the procedure we use to incorporate sample incompletenesses and measurement biases and errors into our model. The transformed distribution is shown in the bottom panel of Figure 7, which is the expected distribution if the intrinsic distribution is characterized by $\bm {{P}}_0$. We compare the transformed distribution with the observed distribution to derive the best-fit parameters.

Zoom In Zoom Out Reset image size

4.2.1. Description of Simulation

4.2.2. Calculating Selection Completeness

The dropout simulations quantify the LBG selection function P(M, R_e, z) within the absolute-magnitude bin (M, M + ΔM), effective-radius bin (log₁₀R_e, log₁₀R_e + Δlog₁₀R_e), and redshift bin (z, z + Δz). By definition in each bin of (M, R_e, z),

$$\begin{equation} P(M, R_e, z) \equiv \frac{\mathrm{{\rm number}\ of\ detected\ and\ selected\ sources}}{\mathrm{total\ {\rm number}\ of\ input\ sources}}. \end{equation} \tag{ 8 }$$

In other words, P(M, R_e, z) is the probability that an LBG in the given (M, R_e, z) bin is detected in the image and selected in our sample. We adopt narrow bin sizes of (ΔM, Δlog₁₀R_e, Δz) = (0.5, 0.2, 0.1) so that the result is insensitive to the input distribution in either M, R_e, or z (as long as the input distribution is not pathological).

Figure 8 illustrates the dependence of the dropout-selection function P(M, R_e, z) on absolute magnitude M and effective half-light radius R_e derived from the simulations in the GOODS-depth images. The top row shows the P(M, R_e, z) for the B₄₃₅-dropout selection; the bottom row shows the same for the V₆₀₆-dropout selection. For both B₄₃₅-dropouts and V₆₀₆-dropouts, the selection function is skewed more toward the low-redshift end in fainter bins as well as in larger-effective-radius bins, showing that the redshift distribution of the dropout sample depends on both absolute magnitude and effective half-light radius. The plots for the HUDF-depth images are qualitatively similar, but they extend to fainter absolute magnitudes.

Zoom In Zoom Out Reset image size

4.2.3. Applying P(M, R_e, z)

As the first step toward building the expected observed size–luminosity distribution, we incorporate the dropout-selection incompleteness into Equation (7) using P(M, R_e, z). We calculate the expected distribution θ with the dropout-selection completenesses as

$$\begin{eqnarray} \theta (m, R_e) &= \int _{-\infty }^{M_{\mathrm{lim}}}\int _{0}^{\infty } \Psi (M, R_e)\,P(M, R_e, z)\,\frac{dV}{dz}\,\delta _M dMdz \nonumber\\ \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \delta _M &\equiv \delta \left[M-\mathrm{DM}(z)-K-m\right], \end{eqnarray} \tag{ 10 }$$

where dV/dz is the volume element contained in the surveyed area and redshift slice (z, z + dz), and δ_M is the Dirac delta function that facilitates the transformation from 1500 Å absolute magnitude M to apparent magnitude m. DM(z) is the distance modulus to redshift z and K is the k-correction term that converts from 1500 Å to the observed bandpass. The integration limit in absolute magnitude M_lim is chosen to include the faintest dropout candidate possible that could be selected in our sample. Although we integrate over all redshifts, P(M, R_e, z) is only non-zero within the redshift range targeted by the selection criteria. By multiplying the volume element dV/dz, the distribution θ predicts the number counts in the (m, log₁₀R_e) plane.

When implementing the correction, the integration becomes a discrete sum over small bins of (M, R_e, z) because we only have a discrete sampling of P(M, R_e, z). Therefore, we define the bivariate model Ψ on a discrete grid with (dM, dlog₁₀R_e) = (0.02, 0.02). Together with P(M, R, _ez) defined in bins of (ΔM, Δlog₁₀R_e, Δz) = (0.5, 0.2, 0.1), we calculate the distribution θ also on a discrete grid with (dm, dlog₁₀R_e) = (0.02, 0.02). The details of such calculations are elaborated on in Appendix A.1.

4.3.1. Description of GALFIT Simulation

4.3.2. GALFIT Transfer Function Kernels

Using the output simulation catalog, we calculate the differences between the input and output magnitudes and effective half-light radii for all sources. We denote the input magnitude and effective half-light radius as m and R_e, and the output (measured) values as m' and R'_e. The differences between the input and the measured magnitudes and sizes are defined as δm ≡ m' − m and δlog₁₀R_e ≡ log₁₀R'_e − log₁₀R_e, respectively. Then in the i'th magnitude bin (in which $m_{i^{\prime }} < m < m_{i^{\prime }}+\Delta m$) and the j'th size bin (in which $(\log _{10} R_e)_{j^{\prime }} < \log _{10} R_e < (\log _{10} R_e)_{j^{\prime }} + \Delta \log _{10} R_e$), we calculate the transfer function kernel $T_{i^{\prime }j^{\prime }}^\mathrm{GF}$ and completeness $C^{\mathrm{GF}}_{i^{\prime }j^{\prime }}$ as

$$\begin{eqnarray} T_{i^{\prime }j^{\prime }}^{\mathrm{GF}} &\equiv \mathcal {P}(\delta m, \delta \log _{10} R_e) \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} C^{\mathrm{GF}}_{i^{\prime }j^{\prime }} &= \int \int \mathcal {P}(\delta m, \delta \log _{10} R_e)\,d\delta m\,d\delta \log _{10} R_e, \end{eqnarray} \tag{ 12 }$$

where $\mathcal {P}(\delta m, \delta \log _{10} R_e)$ is the probability density function (PDF) in the plane of (δm, δlog₁₀R_e). The GALFIT-measurement completeness $C^\mathrm{GF}_{i^{\prime }j^{\prime }}$ accounts for the fact that GALFIT does not always give satisfactory fits. (Here $C^\mathrm{GF}_{i^{\prime }j^{\prime }}$ does not include the SExtractor-detection incompleteness.) We use the same criteria as described in Section 3 to calculate the fraction of sources poorly fit by GALFIT, namely, those with GALFIT-returned errors [R_e]_err/R_e ⩾ 0.6 and χ²_ν > 0.4 in the GOODS-depth images (and χ²_ν > 0.5 in z₈₅₀-band), and [R_e]_err/R_e ⩾ 0.6 and χ²_ν > 5.0 in the HUDF-depth images in both i₇₇₅- and z₈₅₀-band.

The transfer function kernels $T_{i^{\prime }j^{\prime }}^{\mathrm{GF}}$ in different bins are shown in Figure 9. The panels on the lower left have high surface brightnesses, while the panels on the upper right have low surface brightnesses. In general, surface brightness has an important effect on the behavior of the GALFIT measurement, namely kernels in low-surface brightness bins are more spread out (larger uncertainties) and center farther away from the origin (larger biases). The errors in the measured magnitudes and effective half-light radii are correlated as expected: an overestimation of flux (δm < 0) leads to an overestimation of effective half-light radius (δlog₁₀R_e > 0). GALFIT performs well overall; even in low-surface brightness bins, GALFIT still shows little bias compared with SExtractor.

Zoom In Zoom Out Reset image size

4.3.3. Comparison with SExtractor

One can also calculate the transfer function kernels for SExtractor-measured values. The kernels derived for SExtractor are shown in the right half of Figure 9. One can see the biases of SExtractor-measured values: SExtractor tends to overestimate the effective half-light radius in high-surface brightness bins (lower left) and underestimate the effective half-light radius in low-surface brightness bins (upper right). This is the main reason we decide to use GALFIT instead of SExtractor for size measurement.

We also show the plots of the measured magnitudes and effective half-light radii for very narrow ranges of input magnitudes around m = 24, 25, and 26 mag in Figure 10. One can also see here that SExtractor in general systematically returns magnitudes that are fainter than the input values. This indicates that SExtractor systematically misses a portion of the light from simulated sources. Although GALFIT measurements show similar scatter, there is no obvious bias in Figure 10. These comparisons demonstrate the advantages of profile fitting, at least in the case of smooth galaxies. In future efforts, we plan to quantify the biases for realistic—and possibly clumpy—galaxies drawn from hydrodynamical simulations.

Zoom In Zoom Out Reset image size

4.3.4. Applying T^GF

We apply the GALFIT transfer function kernels $T^{\mathrm{GF}}_{i^{\prime }j^{\prime }}$ to θ(m, R_e) in Equation (9). This is done by convolving θ(m, R_e) within the i'j'th bin of (m, R_e) with the transfer function kernel pertinent to that bin, and summing the contribution from all bins together:

$$\begin{eqnarray} &&\psi _{i^{\prime }j^{\prime }}(m^{\prime }, \log _{10} R_e) = \theta _{i^{\prime }j^{\prime }}(m, \log _{10} R_e) \ast T^\mathrm{GF}_{i^{\prime }j^{\prime }} (\delta m,\delta \log _{10} R_e) \nonumber\\ \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \psi (m^{\prime },\log _{10} R^{\prime }_e) &= \sum _{i^{\prime }j^{\prime }}\,\psi _{i^{\prime }j^{\prime }}(m^{\prime },\log _{10} R^{\prime }_e). \end{eqnarray} \tag{ 14 }$$

The distribution ψ(m', log₁₀R'_e) is the expected distribution in measured magnitude and effective half-light radius corresponding to the set of parameters $\bm {{P}}$. Therefore, ψ(m', log₁₀R') accounts for all the systematic effects related to dropout-selection and GALFIT measurement. We define the final distribution on the grid of bin size (dm', dlog₁₀R'_e) = (0.02, 0.02), and the detailed mathematical formulation of Equations (13) and (14) on the grid is included in Appendix A.2. After applying the dropout-selection function P(M, R_e, z) (Section 4.2.3) and GALFIT-transfer-function kernels T^GF, we arrive at the expected distribution of observed dropouts. The distribution corresponding to the fiducial set of parameters $(\alpha , M^*, R_0, \sigma _{\ln R_e}, \beta)= (-1.7, -21.0, 0\buildrel{\prime \prime}\over{.}21, 0.7, 0.3)$ is shown in the bottom panel of Figure 7.

We use the maximum-likelihood estimation (MLE) method to recover the best-fit parameters $\bm {{P}}\equiv (\alpha , M^*, R_0, \sigma _{\ln R_e}, \beta)$. Our formulation is equivalent to that proposed by Sandage et al. (1979; see also Trenti & Stiavelli 2008; Su et al. 2011). In general, given a set of observed data points (x₁, x₂, ..., x_n) and a model Γ(x) that describes the distribution of the observed points, the likelihood of the model Γ(x) is

$$\begin{equation} \mathcal {L} \equiv \Gamma (x_1)\cdot \Gamma (x_2)\cdots \Gamma (x_n) = \prod _{i=1}^{n}\Gamma (x_i). \end{equation} \tag{ 15 }$$

Here we try to determine the set of parameters $\bm {{P}}$ that produces the model ψ(m', log₁₀R'_e) that will maximize the likelihood function $\mathcal {L}$ for the observed magnitudes and sizes. Here we also need to include the component of the interlopers I(m', log₁₀R'_e) described in Section 2.3 (and calculated in Appendix A.3) that accounts for the number counts of possible low-z interlopers due to photometric scatter. Therefore, for each dropout sample we combine the log-likelihood from both fields (GOODS+HUDF) and find the set of parameters that maximizes the combined log-likelihood. In short, for a given set of parameters $\bm {{P}}$ we first calculate the overall normalization constant ϕ* through the equation

$$\begin{eqnarray} N_{\mathrm{G}} + N_{\mathrm{H}} &=& \phi ^*\sum _{s,l} [\psi ^{\mathrm{G}} (m^{\prime }_s, \log _{10} (R^{\prime }_e)_l)+I^{\mathrm{G}}(m^{\prime }_s,\log _{10} (R^{\prime }_e)_l)\nonumber\\ &&+\psi ^{\mathrm{H}}(m^{\prime }_s, \log _{10} (R^{\prime }_e)_l)+I^{\mathrm{H}}(m^{\prime }_s,\log _{10} (R^{\prime }_e)_l)] \nonumber\\ &&\times dm^{\prime }\,d\log _{10} R^{\prime }_e, \end{eqnarray} \tag{ 16 }$$

where N_G and N_H are the total number of dropouts in GOODS and HUDF, respectively. The volume surveyed in each field has already been included in ψ(m', log₁₀R'_e) (see Section 4.2.3 and Equation (9)). We then calculate the sum of the log-likelihood of both fields:

$$\begin{eqnarray} \log _{10} \mathcal {L}_{\mathrm{G}}(\bm {{P}}) &=& N_{\mathrm{G}} \,\log _{10} \phi ^*(\bm {{P}})\nonumber\\ && + \sum _{i=1}^{N_{\mathrm{G}}} \log _{10} [\psi ^{\mathrm{G}}(\bm {{P}};\,m^{\prime }_{i}, \log _{10} (R_e)_i)\nonumber\\ && + I^{\mathrm{G}}(m^{\prime }_i, \log _{10} (R^{\prime }_e)_i)] \end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray} \log _{10} \mathcal {L}_{\mathrm{H}}(\bm {{P}}) & = & N_{\mathrm{H}} \,\log _{10} \phi ^*(\bm {{P}})\nonumber\\ && + \sum _{i=1}^{N_{\mathrm{H}}} \log _{10} [\psi ^{\mathrm{H}}(\bm {{P}};\,m^{\prime }_{i}, \log _{10} (R_e)_{i})\nonumber\\ && + I^{\mathrm{H}}(m^{\prime }_i,\log _{10} (R^{\prime }_e)_i)] \end{eqnarray} \tag{ 17b }$$

$$\begin{equation} \log _{10} \mathcal {L}_{\mathrm{tot}}(\bm {{P}}) = \log _{10}\mathcal {L}_{\mathrm{G}}(\bm {{P}}) + \log _{10}\mathcal {L}_{\mathrm{H}}(\bm {{P}}). \end{equation} \tag{ 17c }$$

We use the scipy package optimize to find the global maximum of $\log _{10} \mathcal {L}_{\mathrm{tot}}$. Within the package optimize, the function fmin_l_bfgs_b uses the L-BFGS-B algorithm (Zhu et al. 1997) to search for global minimum. The best-fit parameters for each dropout sample are presented in Section 5.

We compare our derived best-fit values of the Schechter function with those of other studies at similar redshifts in Figure 15. Our derived α and ϕ* are in general consistent within 1σ of those from most studies; however, our M* is about 0.2–0.4 mag fainter than previous studies. The fainter M* is possibly driven by the low completeness near the faint end, where we potentially miss many galaxies due to low surface brightness. The inclusion of photometric interlopers is likely not the origin of the discrepancy: without photometric interlopers, the best-fit (α, M*, ϕ*) becomes (− 1.64, −20.48, 2.21 × 10⁻³ Mpc⁻¹) at z ∼ 4 and (− 1.67, −20.51, 1.70 × 10⁻³ Mpc⁻¹) at z ∼ 5, statistically indistinguishable from our nominal best-fit parameters. Extending our methodology to include more fields (e.g., HUDF parallel fields; Thompson et al. 2005; Oesch et al. 2007) will better constrain the luminosity function parameters by increasing the sample size.

Zoom In Zoom Out Reset image size

However, we have reason to believe that the cosmic variance introduces very small biases in our best-fit parameter values. Trenti & Stiavelli (2008; see also Robertson 2010) has a comprehensive discussion about the effect of cosmic variance on the determination of luminosity functions at high redshift, especially in narrow, pencil-beam fields such as HUDF. In particular, in their Section 5 they examined how the number density of high-redshift dropouts will bias the best-fit Schechter function parameters. They demonstrate that within a single field, M* tends to be fainter when the field is under-dense; however, when one combines a shallower, larger field with a deeper, smaller field, and derives the luminosity function parameters using the unbinned maximum-likelihood method as in Sandage et al. (1979), biases in both α and M* can be reduced. Because our fitting method is the one proposed by Sandage et al. (1979) generalized to include multiple fields, we do not expect our best-fit α and M* to be strongly biased due to cosmic variance.

The peak of the size distribution at the same rest-frame 1500 Å absolute magnitude, M₁₅₀₀ = −21 mag, is 1.34^+0.099_−0.108 kpc for the B₄₃₅-dropout sample (assuming z = 4) and is 1.19^+0.21_−0.16 kpc for the V₆₀₆-dropout sample (assuming z = 5). If instead we calculate the peak of the size distribution at the best-fit M* of each redshift, we will find that the size distribution peaks at 1.24 kpc at z = 4 (at M* = −20.60 mag), and at 1.07 kpc at z = 5 (M* = −20.53 mag), showing slightly stronger evolution than at fixed luminosity.

The conventional picture of Fall & Efstathiou (1980; see also Fall 1983; Dalcanton et al. 1997; Mo et al. 1998; Fall 2002) states that the size of the gaseous disk R_d, formed in the center of DM halo that acquired angular momentum through tidal interactions, is proportional to the virial radius R_vir of the host DM halo and the spin parameter λ of the halo: R_d∝λR_vir. Assuming that the peak of the size distribution R₀(z) at a given luminosity L₀ traces the peak of the log-normal λ distribution, and the λ distribution correlates weakly with redshift or halo mass (Barnes & Efstathiou 1987), one comes to the conclusion that R₀(z)∝R_vir(z). One can further relate the virial radius R_vir(z) of the halo to the circular velocity V_c and the mass M_H of the halo through (Ferguson et al. 2004)

$$\begin{equation} R_{\mathrm{vir}} = \frac{V_c}{10H(z)} = \left(\frac{GM_\mathrm{H}}{100}\right)^{1/3}\ H^{-2/3}(z), \end{equation} \tag{ 18 }$$

so that R_d∝H⁻¹(z) if R_vir at a fixed luminosity traces constant circular velocity, or R_d∝H^−2/3(z) if R_vir traces constant halo mass, and H(z) is the Hubble parameter. In the concordance cosmology model, the Hubble parameter as a function of z is $H(z) = H_0\sqrt{\Omega _m(1+z)^3+\Omega _\Lambda }$, so at z ⩾ 2, H⁻¹(z) ∼ (1 + z)^−3/2 and H^−2/3(z) ∼ (1 + z)⁻¹. Between z = 4 and 5, both relations predict only 20% and 31% changes in size, respectively. The peak size at M₁₅₀₀ = −21 mag evolves only by 13% from our results, favoring only slightly the (1 + z)⁻¹ relation. One needs data at lower redshifts to constrain the size-evolution scenario, where these two scaling relations diverge.

Studies of size evolution of disk-like (or late-type) galaxies with similar luminosities have generally shown larger fractional changes in size at z ⩽ 3. At z ≲ 1.1, using GEMS (Rix et al. 2004) data, Barden et al. (2005) found that the mean rest-frame V-band surface brightness of late-type galaxies with M_V ≲ −20 mag brightens by about 1 mag from z ∼ 0 to z ∼ 1. This corresponds to a 40% decrease in mean size of the galaxies with similar luminosities. Trujillo et al. (2006), combining FIRES (Franx et al. 2000), GEMS, and SDSS (York et al. 2000) data, found that at rest-frame V-band luminosity L_V > 3.4 × 10¹⁰h⁻²₇₀ L_☉ (or M_V ⩽ −21.53 mag) the mean effective half-light radius of late-type galaxies is only about 30% of that of the local SDSS sample. At high redshift, Ferguson et al. (2004) found the mean half-light radius of galaxies selected in GOODS from z ∼ 1.4 to z ∼ 5 follows the H⁻¹(z) ∼ (1 + z)^−3/2 scaling relation. On the other hand, both Bouwens et al. (2004) and Oesch et al. (2010) found the mean half-light radii of LBGs, selected in HUDF and HUDF ACS-parallel fields, to roughly trace (1 + z)⁻¹ from z ∼ 3 to z ∼ 8, although Oesch et al. (2010) noted that the peak of the size distribution evolves little between z ∼ 4–8. From the above observations of star-forming galaxies, the fractional change in size at similar luminosities is significant between z = 0–3, but smaller at z ≳ 3.

We plot the mean sizes of LBGs at z ≳ 2 (in rest-frame UV) from a compilation of previous papers, as well as our best-fit R₀ at M₁₅₀₀ = −21 mag, as a function of redshift in Figure 16. Also included in the figure are the curves H(z)⁻¹ and H(z)^−2/3, both normalized to the derived R₀ at z = 4 of this work. Due to different measurement methods, definitions of size, corrections for measurement biases and uncertainties, and corrections for incompletenesses, systematic differences in size at the same redshift exist. (The definition of effective half-light radius that we use is also meant to enclose half the total flux from a galaxy, but it does try to compensate for the missing flux at the edge of a galaxy, unlike the half-light radius measured by SExtractor as adopted in all other papers included here.) For example, Ferguson et al. (2004) derived the mean sizes of LBGs at z ∼ 1–5 with rest-frame UV luminosity between 0.7L*_{z = 3} and 5L*_{z = 3}, where L*_{z = 3} is the characteristic luminosity of LBGs at z = 3 from Steidel et al. (1999). They did not correct for sample incompleteness and measurement bias, and they reported the average size instead of the peak of the size distribution as in this paper (the mean of a log-normal distribution is larger than the peak due to the tail at the high end). Comparisons between high-z LBGs and local galaxies are also hampered by different rest-frame wavelengths probed. In order to compare with the rest-frame optical sizes of local galaxies, one needs resolved images of most LBGs at z ≳ 3 in the rest-frame optical. However, currently such data are not available, and the James Webb Space Telescope (JWST) will be necessary for this line of investigation. We also wish to extend our analyses to lower-redshift samples where high-resolution rest-frame UV data are available (down to z ∼ 2; see, e.g., Hathi et al. 2010) to obtain consistent size measurements across a longer redshift baseline.

Zoom In Zoom Out Reset image size

We find that the width of the size distribution, σ_ln R, is 0.83^+0.046_−0.044 at z ∼ 4 and 0.90^+0.15_−0.065 at z ∼ 5. This seems to indicate that the width of the size distribution at fixed luminosity is slightly larger at z ∼ 5, and both are larger than the width of the λ distribution (σ_ln λ ∼ 0.5) of DM halos found in the N-body simulations (Barnes & Efstathiou 1987; Warren et al. 1992; Bullock et al. 2001).

We plot the distribution of sizes of the sample—both measured by GALFIT and by SExtractor—in Figure 17. The observed widths of size distribution, without any statistical correction for biases and uncertainties, are σ_ln R ≳ 0.7 for the GALFIT-measured effective half-light radii. In comparison, the SExtractor-measured half-light radii shows narrower distributions: the widths of the best-fit log-normal distributions are σ_ln R ≳ 0.4. Although the widths of the SExtractor-measured size distribution are more in line with the theoretical expectations (σ_ln λ ∼ 0.5), we believe that it is biased because (1) SExtractor is severely limited by the PSF size, so it is strongly biased (i.e., the measured size is bigger than the true size) when the galaxies are barely resolved, and (2) for large galaxies SExtractor is also strongly biased due to the tendency to miss the tail of the light distribution. We therefore believe that the intrinsic width of the size distribution is significantly wider than what SExtractor-measured distribution indicates, and our statistical treatment indeed shows that the intrinsic size distribution is wider than the that measured by SExtractor.

Zoom In Zoom Out Reset image size

One must also consider the fact that GALFIT (or any size-measurement method) can artificially broaden the observed size distributions due to measurement uncertainties. However, our results show that although GALFIT might artificially broaden the intrinsic size distributions, the size distributions after accounting for sample-selection incompleteness are even wider than the GALFIT-measured distributions. Therefore, sample incompletenesses, especially at the small- and large-size ends, are more important in shaping the observed widths of the size distributions than whichever measurement method one uses.

The widths of the size distribution of local galaxies are in general smaller than what are found here for LBGs. When considering galaxies divided by their morphological types drawn from SDSS (z ∼ 0.1), Shen et al. (2003) found that σ_ln R ∼ 0.3 at the bright end, while the width increases to σ_ln R ∼ 0.5 at the faint end. This is true when considering either late-type galaxies alone or early-type galaxies alone. Furthermore, both de Jong & Lacey (2000) and Courteau et al. (2007) derived that σ_ln R ∼ 0.3 for their local spiral galaxies. All of the above results show that the size distribution for local disk galaxies is narrower than the λ distribution of DM halos, and mechanisms to systematically remove disk galaxies residing in small-λ and large-λ halos such as bar instability and supernova feedback are suggested. Given the wide size distribution derived for LBGs in this work, it is possible that LBGs sample the whole spectrum of λ, and some aspects of baryon physics further broaden the size distribution. In addition, there are also hints of wider size distributions than that of late-types or early-types alone when we consider both types together (e.g., Fall 1983; Romanowsky & Fall 2012). Because it is still unknown whether the majority of LBGs are rotation-supported disks, the relation between the width of the size distributions of LBGs and that of λ is still not clear.

Adopting the power-law form of the size–luminosity relation $\bar{R} \propto L^\beta$, we find that β = 0.22^+0.058_−0.056 for the B₄₃₅-dropouts and β = 0.25^+0.15_−0.14 for the V₆₀₆-dropouts. The trend between the peak of the size distribution $\bar{R}(L)$ and luminosity L, without accounting for incompletenesses at the low-surface brightness end, is steeper due to selection biases against low-surface brightness galaxies. The observed β, without accounting for incompleteness, is 0.32 for the B₄₃₅-dropouts and 0.30 for the V₆₀₆-dropouts. We also plot in Figure 18 the peak of the log-normal size distribution $\bar{R}(L)$ and the width σ_ln R in different apparent magnitude bins, when α, M*, and β are held fixed to the best-fit values. The different scaling relations between the mean size $\bar{R}$ and luminosity L are also plotted. After accounting for incompletenesses, the size–luminosity relation is still inconsistent with being flat.

Zoom In Zoom Out Reset image size

The observed size–luminosity (RL) relation of low-z disk galaxies in the optical bands has been shown to have similar power-law index β as our result for the high-z galaxies. We compile the size–luminosity relations for local galaxies in Table 4. For example, de Jong & Lacey (2000) observed that R_d∝L^0.25 for their local disk galaxy sample with types later than Sb. Shen et al. (2003) observed that R_e∝L^0.26 in the faint end of their late-type galaxy sample (defined as n < 2.5) drawn from SDSS, where R_e is their best-fit Sérsic effective half-light radius, although they also observed that at the bright end of their late-type galaxies the observed slope is significantly steeper (β ∼ 0.5). Courteau et al. (2007) observed that R_d∝L^0.32 for their local disk galaxy sample, including both field and cluster galaxies, and also observed a type-dependent trend in β such that later-type disk galaxies have shallower slope (β = 0.25 for their Sd-type sample). On the other hand, the observed RL relation for local early-type galaxies generally have steeper slopes: β ∼ 0.5–0.6 in general (Shen et al. 2003). We plot the RL relations from local galaxies along with the derived RL relations for out LBGs in Figure 19. The agreement between the slopes of the local RL relations and those for the LBGs in this study is interesting. Although most galaxies in our LBG sample have n < 2.5, which will be classified as late-type galaxies based on Sérsic index only, whether they are disk galaxies is still uncertain. For example, LBGs have ellipticity distributions inconsistent with being drawn from an optically thin disk (Ravindranath et al. 2006; Law et al. 2012). On the other hand, the normalization of the RL relation between different redshifts reveals the amount of size evolution between these two epochs. The similarities between the size–luminosity relation of the local late-type galaxies and high-z LBGs suggests that their formation histories share some similarities, or even that LBGs are the progenitors of local disk galaxies. Furthermore, if one scales the best-fit RL relation at z ∼ 4 according to H(z)^−2/3, the resultant RL relation matches the local relations very well (the dark-gray solid line in Figure 19).

Zoom In Zoom Out Reset image size

de Jong & Lacey (2000; see also Fall 1983; Dalcanton et al. 1997; Mo et al. 1998) offered a derivation of the expected value of β for disk galaxies. We saw in Section 6.2 that in the standard picture of disk formation, R_d∝λR_vir, and the virial radius R_vir is proportional to the circular velocity V_c at a given redshift (Equation (18)). The Tully–Fisher relation (Tully & Fisher 1977) of local disk galaxies in the i-band shows that V_c∝L^β, where β ∼ 0.29 (Courteau et al. 2007). This translates to R_d∝λV_c∝λL^β. The best-fit β of our sample, 0.22 for the B₄₃₅-dropouts and 0.25 for the V₆₀₆-dropouts, are slightly smaller than the prediction by the simple argument, but are still consistent with the TF slope of local disk galaxies.

Here we elaborate on the details of applying the dropout-selection function P(M, R_e, z) on the size–luminosity distribution (Equation (7)) on a grid of (M, log₁₀R_e). We also explain how we apply the dropout-selection function in the two fields separately, and how we combine the two fields to calculate the total likelihood.

For clarity first we define several variables as follows.

• 1.  
  M: the intrinsic 1500 Å absolute magnitude.
• 2.  
  m and R: the intrinsic apparent magnitude and effective half-light radius (we will omit the subscript e in R_e for clarity of notations).

Also the superscripts G stands for GOODS and H stands for HUDF in the equations that follow.

In the GOODS-depth images, the kernel can be expressed as the selection function P^G(M, log₁₀R, z), defined as the fraction of sources with input M, log₁₀R, and z that is detected and selected as a dropout by the color criteria. In this study we calculate the selection function P^G(M, log₁₀R, z) in bins of (M, log₁₀R, z), such that P^G_ijk is the completeness with absolute magnitude in the range (M_i, M_i + ΔM), size in the range (log₁₀R_j, log₁₀R_j + Δlog₁₀R), and redshift in the range (z_k, z_k + Δz). The bin size (ΔM, Δlog₁₀R, Δz) is selected to be (0.5, 0.2, 0.1), and we implicitly assume that the selection function stays the same within the bin. Usually the dropout selection function P^G_ijk is only greater than zero within Δz ∼ 0.5 around the central redshift (z = 4 for B₄₃₅-dropouts and z = 5 for V₆₀₆-dropouts).

The first step toward building the expected bivariate distribution of observed LBGs is to apply the dropout selection function P^G_ijk to the theoretical distribution (Equation (7)). We define Ψ, sans the normalization factor ϕ*, on a grid of (M_p, log₁₀R_q) where p and q are the indices of the "pixels." We choose the pixel size to be (dM, dlog₁₀R) = (0.02, 0.02). We apply the dropout-selection function in the following way: for each kernel within the range M_i ⩽ M < M_i + ΔM and log₁₀R_j ⩽ log₁₀R < log₁₀R_j + Δlog₁₀R, we calculate the contribution to the number density of LBGs predicted by Ψ in this bin, weighted by the volume in each redshift interval and the dropout selection function P^G_ijk. In GOODS,

$$\begin{eqnarray} \theta ^{\mathrm{G}}_{ij}(m_r, \log _{10} R_q) &=& \sum _{p,k} \Psi _{ij}(M_p, \log _{10} R_q)\,P^{\mathrm{G}}_{ijk}\, \frac{dV}{dz_k}\, \nonumber\\ &&\times \delta(M_p - \mathrm{DM}(z_k) - K - m_r) \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} &\equiv \sum _{ij} \Psi _{ij} \circledast P^{\mathrm{G}}_{ijk}, \end{eqnarray} \tag{ A2 }$$

where

$$\begin{eqnarray} && \Psi _{ij}(M_p,\log _{10} R_q) \nonumber\\ && \equiv \left\lbrace \begin{array}{@{} l @{\quad }l@{}}\Psi (M_p,\log _{10} R_q) &\mathrm{if}\ M_i\le M_p < M_i+\Delta M\ \mathrm{and}\\ & \log _{10} R_j \le \log _{10} R_q < \log _{10} R_j \\ & +\, \Delta \log _{10} R\\ 0 &\mathrm{otherwise} \end{array}\right. \qquad \end{eqnarray} \tag{ A3 }$$

and dV/dz_k is the volume element within the surveyed area and the redshift range (z_k, z_k + Δz). The Dirac delta function δ(M_p − DM(z_k) − K − m_r) facilitates the conversion from absolute magnitude M to apparent magnitude m so that the result is a distribution θ^G_ij defined on the grid of (m, log₁₀R) with the same pixel size (dm, dlog₁₀R) = (0.02, 0.02). In the argument of the delta function, DM(z_k) is the distance modulus to redshift z_k (calculated using the assumed cosmological parameters), and K is the k-correction term that converts from rest-frame 1500 Å magnitude to the apparent magnitude in the chosen band. Since we choose to use the bandpasses that are closest to rest-frame 1500 Å, the k-correction term is in general very small. Thus Equation (A1) is the contribution to the distribution of (m, log₁₀R) from the model Ψ within M_i ⩽ M < M_i + ΔM and log₁₀R_j ⩽ log₁₀R ⩽ log₁₀R_j + Δlog₁₀R. Lastly, Equation (A2) is a shorthand notation for the application of the dropout-selection function to the model Ψ.

The resultant distribution of number density, combining Ψ in all absolute-magnitude and effective-radius bins and after applying the dropout selection function, is

$$\begin{equation} \theta ^{\mathrm{G}}(m_r, \log _{10} R_q) = \sum _{i,j} \theta ^{\mathrm{G}}_{ij}(m_r, \log _{10} R_q). \end{equation} \tag{ A4 }$$

Similarly, we construct the distribution after applying the dropout selection in the HUDF using the same procedure:

$$\begin{eqnarray} \theta ^{\mathrm{H}}_{ij}(m_r, \log _{10} R_q) &=& \sum _{p,k} \Psi _{ij}(M_p, \log _{10} R_q)\, P^{\mathrm{H}}_{ijk}\,\frac{dV}{dz_k}\, \nonumber\\ &&\times \delta(M_p - \mathrm{DM}(z_k) - K - m_r) \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} &\equiv \Psi _{ij} \circledast P^{\mathrm{HUDF}}_{ijk} \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} \theta ^{\mathrm{H}}(m_r, \log _{10} R_q) &= \sum _{i,j} \theta ^{\mathrm{H}}_{ij}(m_r, \log _{10} R_q). \end{eqnarray} \tag{ A7 }$$

The second set of transfer-function (TF) kernels quantifies the GALFIT-measurement bias, scatter, and completeness of GALFIT giving reasonably good fits. It relates the intrinsic apparent magnitude and effective half-light radius to the observed apparent magnitude and effective half-light radius. Therefore we define

• 1.  
  m', R': the measured apparent magnitude and effective half-light radius;
• 2.  
  δm ≡ m_measured − m_intrinsic = m' − m;
• 3.  
  δlog₁₀R ≡ log₁₀R_measured − log₁₀R_intrinsic = log₁₀R' − log₁₀R.

Again the superscripts G stands for GOODS and H stands for HUDF in the equations that follow.

Similar to the case of dropout selection TF kernels, we use a bin size of 0.5 in m and 0.2 in log₁₀R and calculate GALFIT TF kernels in each intrinsic magnitude and size bin, assuming that the kernel property is uniform within the bin. The GALFIT TF kernels in each bin are PDF defined on a pixel grid of (δm, δlog₁₀R) within pixel width (dm, dlog₁₀R) = (0.02, 0.02). The value of the GALFIT TF kernel at (δm, δlog₁₀R) is the probability that the measured magnitude of the object is off by an amount of δm, and the size is off by an amount δlog₁₀R due to measurement. In the i'th apparent magnitude bin and the j'th size bin, in which $m_{i^{\prime }} \le m < m_{i^{\prime }} + \Delta m$ and $\log _{10} R_{j^{\prime }}\le \log _{10} R < \log _{10} R_{j^{\prime }}+ \Delta \log _{10} R$, the GALFIT TF kernel in GOODS is denoted as $T^{\mathrm{GF},\mathrm{G}}_{i^{\prime }j^{\prime }}(\delta m, \delta \log _{10} R)$.

Since each kernel $T^{\mathrm{GF},\mathrm{G}}_{i^{\prime }j^{\prime }}(\delta m, \delta \log _{10} R)$ only applies in its own bin, we calculate the contribution to the expected distribution of observed magnitude and size in each bin as

$$\begin{eqnarray} \psi ^{\mathrm{G}}_{i^{\prime }j^{\prime }}(m^{\prime }_s,\log _{10} R^{\prime }_l) &=& \sum _{q,r}\theta ^{\mathrm{G}}_{i^{\prime }j^{\prime }} (m_r,\log _{10} R_q)\,T^{\mathrm{GF},\mathrm{G}}_{i^{\prime }j^{\prime }} \nonumber\\ &&\times (m^{\prime }_s-m_r, \log _{10} R^{\prime }_l-\log _{10} R_q) \nonumber\\ \end{eqnarray} \tag{ A8a }$$

$$\begin{eqnarray} \equiv &\ \theta ^{\mathrm{G}}_{i^{\prime }j^{\prime }} \ast T^{\mathrm{GF},\mathrm{G}}_{i^{\prime }j^{\prime }}, \, \, \, \end{eqnarray} \tag{ A8b }$$

where

$$\begin{eqnarray} && \theta ^{\mathrm{G}}_{i^{\prime }j^{\prime }}(m_r,\log _{10} R_q) \nonumber\\ && \equiv \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\theta ^{\mathrm{G}}(m_r,\log _{10} R_q) & \mathrm{if}\ m_{i^{\prime }}\le m _r< m_{i^{\prime }}+\Delta m\ \mathrm{and}\\ & \log _{10} R_{j^{\prime }} \le \log _{10} R _q< \log _{10} R_{j^{\prime }}\\ &+\, \Delta \log _{10} R\\ 0 & \mathrm{otherwise.} \end{array}\right. \qquad \end{eqnarray} \tag{ A9 }$$

The distribution of (m', log₁₀R'), after applying the GALFIT transfer functions in all bins, is the sum of $\psi ^{\mathrm{G}}_{i^{\prime }j^{\prime }}$:

$$\begin{equation} \psi ^{\mathrm{G}}(m^{\prime }_s, \log _{10} R^{\prime }_l) = \sum _{i^{\prime },j^{\prime }} \psi ^{\mathrm{G}}_{i^{\prime }j^{\prime }}(m^{\prime }_s, \log _{10} R^{\prime }_l). \end{equation} \tag{ A10 }$$

Similarly in the HUDF, for the same dropout sample, one can calculate the expected distribution of (m', log₁₀R') as

$$\begin{eqnarray} \psi ^{\mathrm{H}}_{i^{\prime }j^{\prime }}(m^{\prime }_s,\log _{10} R^{\prime }_l) &=& \sum _{q,r}\theta ^{\mathrm{H}}_{i^{\prime }j^{\prime }}(m_r,\log _{10} R_q)\,T^{\mathrm{GF},\mathrm{H}}_{i^{\prime }j^{\prime }} \nonumber\\ &&\times (m^{\prime }_s-m_r, \log _{10} R^{\prime }_l-\log _{10} R_q) \nonumber\\ \end{eqnarray} \tag{ A11a }$$

$$\begin{eqnarray} &\equiv \ \theta ^{\mathrm{H}}_{i^{\prime }j^{\prime }} \ast T^{\mathrm{GF},\mathrm{H}}_{i^{\prime }j^{\prime }} \end{eqnarray} \tag{ A11b }$$

$$\begin{equation} \psi ^{\mathrm{H}}(m^{\prime }_s, \log _{10} R^{\prime }_l) = \sum _{i^{\prime },j^{\prime }} \psi ^{\mathrm{H}}_{i^{\prime }j^{\prime }}(m^{\prime }_s, \log _{10} R^{\prime }_l). \end{equation} \tag{ A12 }$$

Equations (A8b) and (A11b) say that the application of TF kernels is equivalent to convolution of θ and T^GF on a discrete grid in each bin. And Equations (A10) and (A12) say that the final distribution is the sum of the contributions from all (m', log₁₀R') bins.

We estimate the contribution from interlopers due to photometric scatter by combining the estimated interloper fraction as a function of apparent magnitude and size in Figure 2. We do not know exactly the size distribution of the interlopers, so we use the size distribution of all sources in the field as an approximation. In the GOODS fields, if the size distribution of all non-dropout sources is $\mathcal {S}^{\mathrm{G}}(R^{\prime }_l)$, and the interloper number density as a function of apparent magnitude is $\mathcal {I}^{\mathrm{GOODS}}(m^{\prime }_s)$, then the interloper contribution is

$$\begin{equation} I^{\mathrm{G}}(m^{\prime }_s, R^{\prime }_l) = \mathcal {I}^{\mathrm{G}}(m^{\prime }_s)\,\mathcal {S}^{\mathrm{G}}(R^{\prime }_l). \end{equation} \tag{ A13 }$$

Similarly, in HUDF

$$\begin{equation} I^{\mathrm{H}}(m^{\prime }_s, R^{\prime }_l) = \mathcal {I}^{\mathrm{H}}(m^{\prime }_s)\,\mathcal {S}^{\mathrm{H}}(R^{\prime }_l). \end{equation} \tag{ A14 }$$

After accounting for dropout-selection and GALFIT-measurement, the expected distributions on a discrete pixel grid in the GOODS and HUDF are constructed as

$$\begin{eqnarray} \psi ^{\mathrm{G}}(m^{\prime }_s,\,\log _{10} R^{\prime }_l) &= \sum _{i^{\prime },j^{\prime }} \left[ \sum _{i,j} \Psi _{ij} \circledast P^{\mathrm{G}}_{ijk} \right] \ast T^{\mathrm{GF},\mathrm{G}}_{i^{\prime }j^{\prime }} \nonumber\\ \psi ^{\mathrm{H}}(m^{\prime }_s,\,\log _{10} R^{\prime }_l) &= \sum _{i^{\prime },j^{\prime }} \left[ \sum _{i,j} \Psi _{ij} \circledast P^{\mathrm{H}}_{ijk} \right] \ast T^{\mathrm{GF},\mathrm{H}}_{i^{\prime }j^{\prime }}. \qquad \end{eqnarray} \tag{ A15 }$$

We add the interloper contributions (Appendix A.3) I^G and I^H to ψ^G and ψ^H to compare with the measured magnitudes and sizes of our LBG sample and derive the best-fit parameters.