Abundance matching analysis of the emission line galaxy sample in the extended Baryon Oscillation Spectroscopic Survey

The extended Baryon Oscillation Spectroscopic Survey (eBOSS) is a newly completed six-year program (2014 – 2020; Dawson et al., 2016) within the Sloan Digital Sky Survey IV (SDSS-IV; Blanton et al., 2017). It uses the Sloan Foundation 2.5-m Telescope at the Apache Point Observatory (Gunn et al., 2006) with the same BOSS fiber spectrographs (Smee et al., 2013) to conduct observations. The goal of the eBOSS is to constrain cosmological models by measuring the large-scale structure of the universe, and it has successfully measured the cosmological distance scale at the percent level at four different redshifts (Alam et al., 2021). The eBOSS survey uses three different tracers to maximize redshift coverage: emission line galaxies (ELGs), luminous red galaxies (LRGs) and quasars. The eBOSS ELG program utilizes the DECam Legacy Survey (DECaLS; Dey et al., 2018) to perform target selection, because it has deeper $grz-$band photometry than the SDSS imaging (Raichoor et al., 2017). We use the final eBOSS ELG catalog from the SDSS-IV Data Release 16 (DR16) in this paper. It contains 173,736 reliable spectroscopic measurements of redshift between $0.6<z<1.1$, covering a total effective area of 727 deg² (Raichoor et al., 2021).

All abundance matching models should begin with observational measurements of the luminosity function or the mass function of the sample being modeled. Although we are modeling a sub-class of galaxies, this is still a critical piece of constructing our model. From the eBOSS ELG sample we measure the [$\rm OII$]3727 luminosity function, which quantifies the abundance of the ELG at a given [$\rm OII$] luminosity. The reason we choose [$\rm OII$] luminosity function instead of stellar mass function is that [$\rm OII$] flux can be obtained directly from the spectra, while stellar masses of ELGs are measured in indirect ways such as spectral energy distribution (SED) fitting.

The [$\rm OII$] flux of each ELG is extracted using the IDLSPEC2D eBOSS pipeline. The [$\rm OII$] luminosity of each ELG is then calculated as

where $f[\rm{OII}]$ is the [$\rm OII$] flux and $D_{L}(z)$ is the luminosity distance at redshift $z$. We show in Fig. 1 the observed [$\rm OII$] luminosity of the eBOSS ELG sample, measured in the luminosity range of $10^{40.5}\rm{erg~{}s^{-1}}$ $<L[\rm OII]<10^{42.5}\rm erg~{}s^{-1}$ with logarithmic bins of $\Delta\log L[\rm OII]=0.05$ dex. The errors are estimated using a bootstrap method. The Gaussian shape of the measured luminosity function is a product of the standard schechter form for such data and the eBOSS ELG selection function, which becomes inefficient at low luminosities.

Besides the one-point statistics of galaxies, i.e., the galaxy luminosity function, another important tool to study the statistics of galaxies is the two-point correlation function (2PCF) $\xi(r)$. It describes the excess probability to find a galaxy pair with separation $r$ over a random galaxy field (Peebles, 1980). Quantifying the 2PCF of galaxies is useful because it distinguishes models with different galaxy bias at scales of $r\gtrsim 1\rm Mpc$ and constrains the fraction of satellite galaxies at scales of $r\lesssim 1\rm Mpc$. Bias $b$ is defined as the ratio of clustering relative to clustering of dark matter, $b^{2}=\xi/\xi_{\rm DM}$. The bias of dark matter halo is not sensitive to the halo mass for small halos, however it is very sensitive to the halo mass for halos with higher mass. The large scale bias provides information on the overall halo mass scale for the galaxy sample, while the small scale clustering constrains the fraction of the sample that are satellite galaxies.

We use the Landy & Szalay estimator (Landy & Szalay, 1993) to measure the projected 2PCF $w_{p}(r_{p})$ for the eBOSS ELG sample. We correct the systematic effects on the clustering using the weighting scheme described in Raichoor et al. (2021). Each individual ELG is weighted by

where $w_{\textnormal{FKP}},w_{\textrm{sys}},w_{\textrm{cp}}$ and $w_{\textrm{noz}}$ are the FKP weights (Feldman et al., 1994), the imaging systematic weights, the close-pair weights, and the redshift failure weights, respectively. We refer readers to Raichoor et al. (2021) for details of ELG systematic weights. We choose 17 logarithmic bins from $0.34h^{-1}\textrm{Mpc}$ to $70h^{-1}\textrm{Mpc}$, covering both the one-halo term and two-halo term ranges. The errors are estimated with the jackknife resampling technique using 25 subsamples. Our measurements are presented in Fig. 2.

Estimating the covariance matrix directly from the data is noisy given the limited size of ELG sample. Thus to quantify the uncertainties in ELG clustering, we utilize the GLAM-QPM mocks from Lin et al. (2020). This set of mocks model the linear bias of eBOSS ELG sample at the scales we are interested. We first compute the covariance matrix $C_{ij}=\textrm{cov}(\xi_{i},\xi_{j})$ from 2000 GLAM-QPM mocks, where $\xi_{i}$ and $\xi_{j}$ are the clustering measurements at the $i^{\textrm{th}}$ and $j^{\textrm{th}}$ bin, respectively. Then we up-scale it by the jackknife errorbar of $\xi_{\rm ELG}$,

where $e_{i}$ and $e_{j}$ are the errors of the clustering measurements at the $i^{\textrm{th}}$ and $j^{\textrm{th}}$ bin, respectively.

In this paper we use the dark matter halo catalogs from the MultiDark Planck 2 simulation (MDPL2¹¹1https://www.cosmosim.org/cms/simulations/mdpl2/; Prada et al., 2012; Klypin et al., 2016). The MDPL2 simulation adopts a flat $\Lambda$CDM cosmology with cosmological parameters described in §1, and is run by the L-GADGET-2 code, which is a variant of cosmological code GADGET-2 (Springel, 2005). The box size of the simulation is $1h^{-1}\textrm{Gpc}$, and the mass resolution is $1.51\times 10^{9}h^{-1}M_{\odot}$. We choose the simulation output at redshift $z=0.819$, roughly equal to the ELG effective redshift $z_{\textrm{eff}}=0.845$. Halos are found using the RockStar algorithm (Behroozi et al., 2013). The peak maximum circular velocity $V_{\rm peak}$ and the peak halo mass $M_{\rm peak}$ are determined throughout the merger history. The attributes in the dark matter halo catalogs that are used in this paper and their definitions are presented in Table 2. The reason we choose MDPL2 simulation is two folds. First it has high mass resolution, allowing us to locate halos at low mass end (around $10^{11}h^{-1}M_{\odot}$). This is crucial in our study because ELGs are young galaxies that occupy lower-mass halos. Second, the volume of the simulation is large enough compared with eBOSS ELG footprint for us to get a good comparison to the data.

As described in §3.2, the 3 free parameters in our abundance matching model are $\sigma(\log L[\rm OII]|\log SFR)$, $L_{crit}$ and $\sigma_{\log L}$. Our implementation of the model is able to quickly produce mock galaxy catalogs and accurately compute small scale clustering for a given parameter set, making it ideal to be constrained using MCMC method. For each trial parameter set in the MCMC chain, we calculate a likelihood of the observed $L[\rm OII]$ luminosity function and $w_{p}(r_{p})$. We use the emcee package (Foreman-Mackey et al., 2013), which implements the affine-invariant ensemble sampler (Goodman & Weare, 2010), to perform the MCMC analysis.

For the $L[\rm OII]$ luminosity function, we use KL-divergence (also called relative entropy) as the metric. Suppose we have the normalized luminosity function from mock $n_{\rm model}(L)$ and the normalized luminosity function from ELG observation $n_{\rm ELG}(L)$, the KL-divergence of $n_{\rm model}(L)$ from $n_{\rm ELG}(L)$ is defined as:

Note that the luminosity function is normalized here, i.e.,

The KL-divergence measures the information gain when using $n_{\rm model}(L)$ to approximate $n_{\rm ELG}(L)$. With smaller KL-divergence, the mismatch between $n_{\rm ELG}(L)$ and $n_{\rm model}(L)$ is smaller. The KL-divergence is always non-negative, with $D_{\rm KL}(n_{\rm ELG}|n_{\rm model})=0$ if and only if $n_{\rm ELG}(L)=n_{\rm model}(L)$ everywhere.

We forbid parameter sets that yield too much discrepancy on $L[\rm OII]$ between model and ELG data by manually restricting $D_{\rm KL}(n_{\rm ELG}|n_{\rm mock})<0.01$. One can interpret this constraint as if we have a prior probability distribution $n_{\rm model}$ and we observe a posterior probability distribution $n_{\rm ELG}$, then the log-likelihood must be greater than -0.01. The reason we use this approach instead of using $\chi^{2}$ to constrain on $L[\rm OII]$ is that we have a relatively more accurate measurement on the $[\rm OII]$ luminosity than the two-point clustering. So if we combine the $\chi^{2}$ of $w_{p}$ with the $\chi^{2}$ of $L[\rm OII]$ as the full log-likelihood in the MCMC, the posterior distribution of parameters will be dominated by matching the luminosity function, which is not what we intended to do. Given the reason above, we choose a reasonable hard threshold cut for constraining the luminosity function and let the MCMC focus on matching the two-point correlation function.

In order to reproduce the small scale clustering of the eBOSS ELG sample, we constrain parameters using the projected correlation function $w_{p}(r_{p})$. Let $\xi_{\rm model}$ denote the projected correlation function from our model given a given set of parameters and $\xi_{\rm ELG}$ denote the projected correlation function from eBOSS ELG sample. We utilize the covariance matrix $C(\xi_{i},\xi_{j})$ from GLAM-QPM mocks to describe the distance between $\xi_{\rm model}$ and $\xi_{\rm ELG}$. Assuming the $\xi_{\rm model}$ has a multivariate normal distribution around its expected value, the log-likelihood of obtaining $\xi_{\rm model}$ is:

where $k$ is the dimension of $\xi_{\rm model}$, in our case $k=17$. This assumption can be validated through 2000 GLAM-QPM mocks, we measure cross-$\chi^{2}$ with 12 data points between $0.34h^{-1}\textrm{Mpc}$ and $70h^{-1}\textrm{Mpc}$ and the results indeed follow $\chi^{2}$ distribution with 12 degrees of freedom. Since $k$ and the covariance matrix $C$ are constants, the first term of $\log P(\xi_{\rm model})$ is a constant and doesn’t change if our model parameters change.

The main goal of this work is to construct a model for the one-point and two-point statistics of the eBOSS ELG sample with conditional abundance matching. We run 600,000 MCMC iterations to constrain the three free parameters described in §3. In Fig. 1 we present the comparison between the $[\rm OII]$ luminosity function of the eBOSS ELG sample and that from our abundance matching model. Within the space of models that pass the constraint of $D_{\rm KL}(n_{\rm ELG}|n_{\rm mock})<0.01$, we find a good match to the $w_{p}(r_{p})$ data. Note that although the discrepancy between the ELG luminosity function and our model is small, the errorbars of the ELG luminosity function, estimated using bootstrap method, are even smaller. This illustrates the point we made in §3.3 that if we add $\chi^{2}$ of $L[\rm OII]$ in MCMC to optimize the parameters, the model will mainly focus on matching $L[\rm OII]$ instead of $w_{p}$. So our approach focuses on the two-point clustering instead but also has a good constraint on the one-point statistics.

In Fig. 2 we show comparison of projected two-point correlation function between the eBOSS ELG sample and our model result, with the shaded regions as the 1$\sigma$ error. Overall the model and the ELG $w_{p}(r_{p})$ agree well at the 1$\sigma$ level. The best fit of $\chi^{2}$ is $10.0$ with 17 degrees of freedom for $w_{p}$. The corresponding p-value is not significant at $p<0.05$, indicating that there is no statistical significant between our fitted model and the observed data. The reduced $\chi^{2}$ is slightly below 1, which is likely due to the method that we up-scale the errorbar of $\xi_{\textnormal{ELG}}$ in Sec. 2.3. At scales of $r_{p}\gtrsim 10h^{-1}\rm Mpc$ our model produces relatively lower clustering, and at $r\lesssim 10h^{-1}\rm Mpc$ our $w_{p}$ is higher than the eBOSS ELG sample. This is likely due to the reason that our quenched fraction approach yields a slightly more satellites than the eBOSS ELG sample, so our model predicted $w_{p}$ is higher at the transition scale between 1-halo term and 2-halo term.

The constraints on the parameters are presented in Fig. 3. We find that the scatter in $\log L[\rm OII]$ at fixed SFR is $\sigma(L[\rm O_{II}]|SFR)=0.11_{-0.04}^{+0.03}$ dex. While most of the literature focus on estimating SFR using [OII] flux, there is limited discussion of the scatter in $[\rm OII]$ luminosity. However we could utilize the results from Moustakas et al. (2006), where they find that $\sigma(\textnormal{SFR}|L[\textnormal{OII}])\sim 0.3$ dex for various dust corrections. Since the $L[\rm OII]$ is proportional to SFR according to Eq. 5, and the empirical correction term is negatively correlated with SFR, the slope of $\log L[\rm OII]$ against SFR should be less than 1. Thus the $\sigma(L[\rm OII]|SFR)$ should be less than $\sigma(\rm SFR|L[OII])\sim 0.3$, which justifies our finding.

We find that the incompleteness coefficient in Eq. 6 is $\eta=0.12$ under the best MCMC fit. This means that the eBOSS ELG sample selects 12% of star-forming galaxies at mass range $10^{10}h^{-1}M_{\odot}\sim 10^{11.5}h^{-1}M_{\odot}$ from our mock catalog. This is in accordance with the HOD results of Guo et al. (2019), where they find that the eBOSS ELG target selection only select 1-10% of the star-forming galaxies at high mass end for various redshifts. Guo et al. (2019) also find that the completeness decreases at the low mass end. Accordingly, we use an error function to mimic the target selection at the low mass end and assume a universal incompleteness rate. Our best fit $\sigma_{\log L}=0.63$ dex, which model the sharpness of the $[\rm OII]$ selection due to the ELG target selection at the low $[\rm OII]$ end.

In our conditional abundance matching approach we connect the SFR of galaxies with the accretion rate of halos in order to model the galaxy assembly bias. Here we define galaxy assembly bias as the dependence of galaxy properties on a second halo property that is related with halo assembly history. To test the impact of assembly bias, we perform a parallel analysis that repeats the fiducial model, with the exception of removing any correlation between SFR and halo growth rate. In that case, SFR only depends on stellar mass. In Fig. 4 we show the comparison between our best fit model and the model without the assembly bias. The two models have the same number density and satellite fraction. From the left panel in Fig. 4 we can see that both two models can match the $[\rm OII]$ luminosity well, and from the right panel we find that the two-point correlation function without assembly bias is generally $\sim 20\%$ higher than the model with assembly bias. By measuring the $\chi^{2}$ of the model without assembly bias, we find that removing the assembly bias causes 2PCF to deviate from the best fit and fits worse than our fiducial model. This result is expected: if the assembly bias is not included in our model, the star-forming galaxies will more likely reside in older halos given a fixed halo mass bin from $10^{11}h^{-1}M_{\odot}$ to $10^{14}h^{-1}M_{\odot}$. Thus the two-point clustering of galaxies will be higher because older halos are more clustered.

We also test the impact of galaxy assembly bias on the void probability function (VPF). The VPF $P_{0}(r)$ is defined as the probability of a randomly placed sphere with radius $r$ containing no galaxy of a certain type. Compared with the two-point statistics, it is a useful tool to test whether halo occupation or galaxy bias changes from high to low density environment (Tinker et al., 2006). We make measurements for VPF in the range of $r\sim[1,21]h^{-1}\rm Mpc$ with linear binning of $\Delta r=2h^{-1}\rm Mpc$. For each radius we place $10^{6}$ spheres randomly in our mock catalog to obtain a relatively precise measurement of VPF. The errors on VPF are estimated using bootstrap method with 1000 resampling size. In Fig. 5 we show the difference of VPF with and without galaxy assembly bias. The model without assembly bias produces around 5% larger VPF at scale of $r\sim 20h^{-1}\rm Mpc$, as shown in the lower panel of Fig. 5. Note that the difference is more than twice of the errorbar, meaning that it is statistically significant, and can be used as a test with current and upcoming data.

We investigate the mean halo occupation function of our abundance matching model, as shown in Fig. 6. The centrals and the satellites are shown in black dashed line and black dotted line, respectively. The result indicates that ELG central galaxy can be modeled as a log-normal distribution at low mass end and a decreasing power-law at high mass end. According to our model, the occupation numbers of ELG at high halo mass end decreases, but not dramatically. We compare our result with Guo et al. (2019) at redshift $0.8<z<0.9$, as shown in the red line. Our amplitude of HOD function is higher than Guo et al. (2019) because we only show their results for redshift bin $0.8<z<0.9$. Guo et al. (2019) finds a significant decrease of the occupation function for centrals at high mass end, likely because they use a higher quenched fraction than we do. We also compare our result with Lin et al. (2020), and we find a similar satellite HOD. The difference in the central HOD is due to Lin et al. (2020) adopting a Gaussian function for central galaxies.

Recently, Avila et al. (2020) explored different forms of HOD model to study the eBOSS ELG sample. They have tested three different formulations for the central HOD: a step-wise erf function, a Gaussian function, and a Gaussian + decaying power-law at high mass end. The difference between those three formulas is mainly at the high mass end. With the same power-law for the satellite HOD, they find that the shape of the central HOD has a less significant impact on the two-point clustering than the satellite fraction. The data marginally rules out the Gaussian formula, and there is negligible difference between the step-wise erf function and the Gaussian + decaying power-law, which is in accordance with our result.

According to our best fit model, the satellite fraction is $19.30\pm 0.02$%, agreeing with an overall $\sim 20$ percent satellite fraction for all star-forming galaxies at $z=0.88$ (Tinker et al., 2013). Favole et al. (2016) also found that the satellite fraction of ELGs at redshift 0.8 is $22.5\pm 2.5\%$.