Angular clustering and host halo properties of [OII] emitters at z >1 in the Subaru HSC survey

Teppei Okumura,1,2^*
Masao Hayashi,3
I-Non Chiu,1
Yen-Ting Lin,1

Ken Osato,4
Bau-Ching Hsieh,1
Sheng-Chieh Lin1,5

1 Institute of Astronomy and Astrophysics, Academia Sinica, No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan
2Kavli Institute for the Physics and Mathematics of the Universe (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
3 National Astronomical Observatory, Mitaka, Tokyo 181-8588, Japan
4Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 75014 Paris, France
5Department of Physics and Astronomy, University of Kentucky, 505 Rose Street, Lexington, KY 40506, USA.

cosmology: observations — galaxies: formation —
galaxies: halos — large-scale structure of
universe — methods: statistical

We study the angular correlation function of star-forming galaxies and
properties of their host dark matter halos at $z>1$ using the
Hyper-Suprime Cam (HSC) Subaru Strategic Program (SSP) survey.  We use [OII] emitters identified
using two narrow-band (NB) filters, NB816 and NB921, in
the Deep/UltraDeep layers, which respectively cover large angular
areas of $16.3 ~{\rm deg}^2$ and $16.9~{\rm deg}^2$.  Our sample
contains 8302 and 9578 [OII] emitters at $z=1.19$ (NB816) and
$z=1.47$ (NB921), respectively.  We detect a strong clustering
signal over a wide angular range, $0.001 < \theta < 1$ [deg], with the
bias $b=1.61^{+0.13}_{-0.11}$ ($z=1.19$) and $b=2.09^{+0.17}_{-0.15}$
($z=1.47$). We also find a clear deviation of the correlation from a
simple power-law form.  To interpret the measured clustering signal,
we adopt a halo occupation distribution (HOD) model that is
constructed to explain the spatial distribution of galaxies selected
by a star formation rate.  The observed correlation function and
number density are simultaneously explained by the best-fitting HOD
model.  From the constrained HOD model, the average mass of halos
hosting the [OII] emitters is derived to be $\log{M_{\rm
eff}/(h^{-1}M_\odot)}=12.70^{+0.09}_{-0.07}$ and
$12.61^{+0.09}_{-0.05}$ at $z=1.19$ and $1.47$, respectively, which
will become halos with the present-day mass, $M\sim 1.5 \times
10^{13}h^{-1}M_\odot$.  The satellite fraction of the [OII] emitter
sample is found to be $f_{\rm sat}\sim 0.15$.  All these values are
consistent with the previous studies of similar samples, but we obtain
tighter constraints
even in a larger parameter space due to the larger sample size from the HSC.
The results obtained for host halos of [OII] emitters in this paper
enable the construction of mock galaxy catalogs
and the systematic forecast study of cosmological constraints
from upcoming emission line galaxy surveys
such as the Subaru Prime Focus Spectrograph survey.

§ INTRODUCTION

Observation of the large-scale structure of the Universe via the
distribution of galaxies in galaxy redshift surveys provides a
powerful tool to reveal the nature of dark matter and dark energy
[Peebles, 1980, Weinberg et al., 2013].  As dark energy started to become
a dominant energy component toward today, the Universe experienced a
transition from decelerating to accelerating expansion at redshift
around unity.  It is thus crucial to investigate the large-scale
spatial distribution of galaxies over a broad redshift range from
$z=0$ to $z>1$.

The relation between galaxies and their host dark matter halos has
been extensively investigated using the prescription of the halo
occupation distribution (HOD) modeling
[Jing et al., 1998, Seljak, 2000, Peacock & Smith, 2000, Ma & Fry, 2000, Scoccimarro et al., 2001, Berlind & Weinberg, 2002, Cooray & Sheth, 2002].
Galaxies selected by the colors, specific star formation rates (sSFR), and stellar masses, have been mainly used as tracers of the
large-scale structure in wide-angle galaxy redshift surveys.  There
are a lot of preceding studies of the HOD modeling for such
populations, and it is known that a simple model proposed by
Zheng et al., 2005 can describe the spatial distribution of these
galaxies (e.g.,
[Zehavi et al., 2005, Blake et al., 2008, Zheng et al., 2009, Abbas et al., 2010, Zehavi et al., 2011, White et al., 2011, Wake et al., 2011, Coupon et al., 2012, de la Torre et al., 2013, Reid et al., 2014, Guo et al., 2014, Koda et al., 2016, Ishikawa et al., 2020]).
Most of the work is, however, limited at redshift less than unity
because these types of galaxies
become harder to target at higher redshifts [Takada et al., 2014, DESI Collaboration et al., 2016]:
for early/red-type galaxies one needs long exposure time to get the $4000$Å   break and
for color-selected galaxies the accuracy of the photometric redshift becomes worse at such redshifts.
On the other hand, statistical properties of dark matter
halos have been studied using the HOD of high-$z$ quasars and
galaxies, i.e., Lyman-break galaxies or Lyman-$\alpha$ emitters at
$z>2$, for relatively narrow angular regions
[Bullock et al., 2002, Moustakas & Somerville, 2002, Hamana et al., 2004, Conroy et al., 2006, Richardson et al., 2012, Kayo & Oguri, 2012, Durkalec et al., 2015, Harikane et al., 2016, Ishikawa et al., 2017, Harikane et al., 2018].

Due to the observational limit of color/sSFR/stellar mass-selected galaxies at higher redshifts, recently the focus of large galaxy surveys has
turned to emission line galaxies (ELGs) that can be targeted
at $z>1$. Such surveys include the Subaru
FastSound survey [Tonegawa et al., 2015, Okada et al., 2016], extended Baryon
Oscillation Spectroscopic Survey (eBOSS: [Dawson et al., 2016]), the
Hobby-Eberly Telescope Dark Energy Experiment (HETDEX:
[Adams et al., 2011]), Subaru Prime Focus Spectrograph (PFS:
[Takada et al., 2014]), Dark Energy Spectroscopic Instrument (DESI:
[DESI Collaboration et al., 2016]), the Spectro-Photometer for the
History of the Universe, Epoch of Reionization, and Ices Explorer
(SPHEREx: [Doré et al., 2014]), and Euclid [Laureijs et al., 2011].
However, statistical properties of the host halos for ELGs are not as simple as those
for galaxies selected by colors, sSFR and stellar masses because unlike them
there does not necessarily exist an ELG in the center of a halo at the massive end.  Geach et al., 2012 proposed a
model of HOD for ELGs taking into account the fact that ELGs selected
via star formation rates do not necessarily reside in the centers of
massive halos. Then they placed a constraint on the HOD parameters for
H$\alpha$ emitters from the Hi-Z Emission Line Survey (HiZELS:
[Geach et al., 2008]).  The HOD of the HiZELS sample has been reanalyzed
by Cochrane et al., 2017 and Cochrane et al., 2018.  Although a simple
HOD modeling has been performed for H$\alpha$ emitters of the
FastSound survey at $z\sim 1.4$ in Okumura et al., 2016, the sample
was so sparse that the clustering signal has been consistent with the
HOD of central galaxies only.  Kashino et al., 2017 analyzed the
H$\alpha$-selected galaxies from the FMOS-COSMOS survey at $z\sim 1.6$
and constrained the HOD model, but they still adopted the simple model
of Zheng et al., 2005.  The detailed HOD modeling has been performed
by Hong et al., 2019 for Ly$\alpha$ emitters at $z\sim 2.67$ selected
from the NOAO Deep Wide-Field Survey.

Since [OII] emitters ($\lambda =$ 3726, 3729 Å) are one of the main tracers of the large-scale
structure employed in upcoming cosmological surveys, it is of crucial
importance to investigate properties of halos which host them.  So
far, however, there are few observational studies of the HOD modeling
for [OII] emitters.  From the HiZELS survey Khostovan et al., 2018
measured the clustering of H$\beta$ + [OIII] as well as [OII]
emitters, but they did not find deviation of the correlation function
from a power-law form.  Rather than the HOD, properties of [OII]
emitters have been investigated using a semi-analytic model of galaxy
formation (e.g.,
[Contreras et al., 2013, Gonzalez-Perez et al., 2018, Favole et al., 2020, Gonzalez-Perez et al., 2020, Avila et al., 2020]).
Guo et al., 2019 analyzed the clustering of
spectroscopically-identified [OII] emitting galaxies at $0.7 < z <
1.2$ from the eBOSS survey.  They then investigated properties of the
host halos using the conditional stellar mass function.  Recently
there are also attempts to determine the HOD of ELGs by directly
identifying emission lines in cosmological hydrodynamical simulations
[Hadzhiyska et al., 2020, Osato & Okumura, 2021].  Properties of the luminosity
function of [OII] and other ELGs are also being actively investigated (e.g.,
[Comparat et al., 2016, Saito et al., 2020, Hayashi et al., 2020, Gao & Jing, 2020]).

Summary of the [OII] ELG data used in this paper.${}^*$

4c|NB816 ($1.178<z<1.208$)        4cNB921 ($1.453<z<1.489$)

Field      Area (${\rm deg}^2$)     $N_g $     $n_g (\ihmpct)$     $N_{\rm sub}$     Area (${\rm deg}^2$)     $N_g$     $n_g(\ihmpct)$    $N_{\rm sub}$

COSMOS     1.37     566    $4.48 \times 10^{-3}$    9     5.80     3134     $4.30 \times 10^{-3}$      41

DEEP2-3     $4.98$     2485     $5.40 \times 10^{-3}$     35     4.92     2107     $3.41 \times 10^{-3}$     35

ELAIS-N1     $4.81$     2927     $6.59 \times 10^{-3}$     34     4.81     3098     $5.13 \times 10^{-3}$     34

SXDS+XMM-LSS     $5.17$      2324     $4.86 \times 10^{-3}$     37     1.33     1239     $7.44 \times 10^{-3}$     9

Total     $16.33$     8302     $5.50 \times 10^{-3}$    115     16.86     9578     $4.53 \times 10^{-3}$     119

$*$ The value of the area shows the one after masking with bright
object masks, $N_g$ is the number of the galaxies used in this paper
after making the magnitude and flux cuts, $n_g$ is the number
densities of the galaxies, and $N_{\rm sub}$ is the number of
sub-regions using jackknife resampling (see section
<ref>). The last row shows the values for the sum of
the four fields.

In this paper we present a detailed HOD modeling for [OII] emitters
identified by Narrow-Band (NB) filters in the Subaru Hyper-Suprime Cam
(HSC) survey.  NB imaging surveys of ELGs allow us to cover a wide and
homogeneous field of view.  Thus it provides a suitable sample to
probe the large-scale spatial distribution of emission-line galaxies.
We, for the first time, constrain the general HOD model proposed by
Geach et al., 2012 for populations of [OII] emitters at $z=1.19$ and
$z=1.47$ in the HSC survey.  We then discuss properties of dark matter
halos which host the [OII] emitters based on the constrained HOD

The structure of this paper is as follows.  Section <ref>
describes the HSC survey and the [OII] emitter sample.  Section
<ref> presents the measurements of the angular
correlation function and its covariance matrix.  Power-law and
linearly-biased dark matter model fittings are shown in section
<ref>. We also discuss dependences of
the correlation function amplitude on different stellar mass and emission line luminosity thresholds.
The detailed HOD modeling is performed for the
measured correlation functions in section <ref>.  Conclusions
are given in section <ref>.  We compare the angular
correlation functions measured from four individual fields in appendix
We present the HOD parameter constraints without using the information of
the measured abundance of [OII] emitters in appendix <ref>.

Throughout this paper, we assume a flat $\Lambda$CDM cosmology with
the cosmological parameters based on the results of Planck CMB
measurements [Planck Collaboration et al., 2016]:
$\Omega_m=1-\Omega_\Lambda = 0.307$, $\Omega_b=0.0486$, $h=0.677$,
$n_s=0.967$ and $\sigma_8=0.816$.

§ DATA

§.§ Hyper-Suprime Cam Survey

The HSC Subaru Strategic Program (SSP) is an ongoing imaging survey since 2014 [Aihara et al., 2018]
with a $1.77~{\rm deg}^2$ field-of-view imaging camera installed on the Subaru Telescope
[Miyazaki et al., 2012, Miyazaki et al., 2018, Furusawa et al., 2018, Komiyama et al., 2018].

In the HSC survey the deep imaging is conducted in five broad-band
(BB) filters as well as four Narrow-Band (NB) filters in the Deep (D)
and UltraDeep (UD) layers over $28~{\rm deg}^2$ in total.  The imaging
reduction and catalog construction are carried out by Bosch et al., 2018
and Huang et al., 2018.  In 2019, the second public data release
(PDR2) of the HSC SSP data has been made
[Aihara et al., 2019].

Angular distribution of [OII] emitters at $z=1.19$ (top) and $z=1.47$ (bottom) with line flux greater than $3\times
10^{-17}{\rm erg ~s}^{-1}{\rm cm}^{-2}$ and magnitude
brighter than $23.5~ {\rm mag}$ in the NB filters, NB816  and NB921,
respectively, for the individual four fields.  The large empty
areas where no galaxy is found are mostly the regions masked for
bright objects.  See figure 4 of Hayashi et al., 2020 for details
of the mask information.

§.§ Narrow Band filters and [OII] emitter sample

We use the data from two NB filters, NB816 and NB921,
constructed by Hayashi et al., 2020 based on the D and UD layers of
Subaru-SSP from the PDR2.  These data are significantly updated from
the version of PDR1 [Aihara et al., 2018, Hayashi et al., 2018].  The coverage
of NB816 and NB921 increases to $26~{\rm deg}^2$ covered
by 13 pointings and $28~{\rm deg}^2$ covered by 14 pointings,
respectively, from the PDR1 [Aihara et al., 2019].  We apply the bright
star masks modified from Coupon et al., 2018 and Aihara et al., 2019
because Hayashi et al., 2020 noticed that the mask size for some
bright stars is not large enough to remove false detections around
prominent stellar halos.  As a result, the total area is reduced to
$16.3~{\rm deg}^2$ (NB816) and $16.9~{\rm deg}^2$ (NB921),
as shown in table <ref>.

The survey area in the D layer, where the NB data are available,
consists of four separate fields: E-COSMOS, DEEP2-3, ELAIS-N1, and
XMM-LSS.  Furthermore, each of E-COSMOS and XMM-LSS fields encompasses the
UD layer covered by a single pointing of HSC, named UD-COSMOS and SXDS,
respectively.  Since UD-COSMOS and E-COSMOS are jointly processed,
we call the combination the COSMOS field.

The complete description of the [OII] emitter sample used in this study
is provided in Hayashi et al., 2018 and Hayashi et al., 2020.
While we briefly summarize our sample in the following, we refer the reader
to these papers for more details.
The ELGs are selected using the NB data together with data from the five
BB filters.
[OII] emitters are
observed by NB816 and NB921 filters in two narrow redshift
ranges, $1.178<z<1.208$ and $1.453<z<1.489$, respectively.
To identify the [OII] emitters, we use spectroscopic redshifts if available and photometric redshifts otherwise.
If a galaxy has a spectroscopic redshift outside of the ranges above, we remove
it as a contaminant. When the photometric redshift is used, we have taken into account
its uncertainty for the identification of [OII] emitters.
For galaxies not identified with spectroscopic or photometric redshifts,
we use the color-color diagrams to distinguish [OII] emitters from
other possibilities.
Galaxies selected as NB [OII] emitters are expected to be star-forming galaxies [Ly et al., 2012, Hayashi et al., 2013].
Thus, our sample is limited by star forming rate (SFR).
While our [OII] emission lines are affected by the dust extinction,
the majority of typical star-forming galaxies is properly selected
because the cosmic SFR density in our catalogs is consistent with other preceding studies
[Hayashi et al., 2020].

We make magnitude and flux cuts
on our data, the model magnitude[We use a composite model magnitude
named .] $m\leq 23.5$ and the estimated line flux of $\geq
3\times 10^{-17}{\rm erg ~s}^{-1}{\rm cm}^{-2}$. We refer the reader to Hayashi et al., 2020 for the details of the line flux estimation.
The limiting magnitude of the NB filters is much deeper than $23.5$ [mag] so that fainter [OII] emitters have been included in our sample, and accordingly [OII] emitters which have the line flux with $(1-2)\times 10^{-17}{\rm erg ~s}^{-1}{\rm cm}^{-2}$ have been selected.
However, we have made such conservative magnitude and flux cuts to assure the completeness of the sample.
With these ranges,
the number density indeed increases monotonically with decreasing the
flux or increasing the magnitude (see
[Hayashi et al., 2018, Hayashi et al., 2020])
and the mean completeness of our [OII] sample becomes $0.974$ for NB816
and $0.956$ for NB921.
The numbers of [OII] emitters are
then reduced to 8302 and 9578 for NB816 and NB921 data,
respectively.  See table <ref> for the numbers of the
sample for each field.  The line flux cut corresponds to the observed
line luminosity thresholds of $L \geq 2.56\times 10^{41} \rm{[erg~s^{-1}]}$
($z=1.19$) and $L \geq 4.30 \times 10^{41} \rm{[erg ~s^{-1}]}$
($z=1.47$).  The median luminosity with the $2.5$th–$97.5$th percentiles
is $\overline{L} = 4.49^{+12.67}_{-1.82} \times 10^{41}
\rm{[erg ~s^{-1}]}$ at $z=1.19$ and $\overline{L} = 7.88^{+16.95}_{-3.35} \times 10^{41}
\rm{[erg ~s^{-1}]}$ at $z=1.47$.
Stellar masses for the [OII] emitters are estimated by
spectral energy distribution fit with five HSC BB data in Hayashi et al., 2020.
The median stellar mass of our [OII] emitter sample with the $2.5$th–$97.5$th percentiles
is $\overline{M}_* = 8.58^{+132.54}_{-7.13}\times 10^{9}M_\odot$
at $z=1.19$ and $\overline{M}_* = 1.72^{+12.40}_{-1.41}\times 10^{10}M_\odot$
at $z=1.47$.

Not all the galaxies in the catalog are real [OII] emitters, but some
are fake lines due to noise and some others are real emission lines
but not [OII] emitters.
The contamination rate is investigated by applying the photo-$z$ and
color selections to galaxies confirmed with spectroscopic redshifts.
The fraction of such non-[OII] emitters has
been significantly reduced in the current sample compared to that in
the PDR1 catalog because (1) the method of selecting emission-line galaxies has been improved from using the  magnitude to the fixed aperture magnitude [Hayashi et al., 2020] and (2)
the pipeline used for the data process, , has been upgraded in the PDR2
[Bosch et al., 2018, Aihara et al., 2019].  The analysis by Hayashi et al., 2020 estimated the
fraction of non-[OII] emitters in our sample is less than $10\%$ for
most cases and $\sim 20\%$ at most, by applying photometric redshift and color selections to galaxies confirmed with spectroscopic redshifts. This contaminant fraction can be
further reduced by a machine-learning based technique for photometric
redshift estimates (see e.g., [Hsieh & Yee, 2014]).  However, since the
current contaminant fraction is low enough for our analysis, we do not
include this process. The angular distribution of the [OII] emitters is shown in
figure <ref>.  Table <ref> summarizes
the properties of the data.

§.§ Random catalog

The HSC-SSP PDR2 provides us with the catalog of random points across
the survey area as one of the value-added products (see subsection 5.4
of [Aihara et al., 2019]). The number density of the random points is 100
per square arcmin. We apply the same basic flags as those for
the emission-line galaxies except for the flags related to the source
detection (i.e., merge_peak, sdssshape_flag, cmodel, and psfflux) to
select the random points (see subsection 2.1.1 of
[Hayashi et al., 2020]). We also apply the same bright star masks as
those used for the selection of the emission-line galaxies to the
random points (see [Hayashi et al., 2020] for the details of the masks
and also subsection 6.6.2 of [Aihara et al., 2019], as well as
[Coupon et al., 2018]). Consequently, the random points used in this study are
distributed over the regions identical to where the emission-line
galaxies are surveyed.

§ ANGULAR CORRELATION FUNCTION

In this paper we analyze the spatial distribution of [OII] emitters
using the angular correlation function.  The angular correlation
function is related to the spatial correlation function in three
dimension through the relation [Peebles, 1973, Peebles, 1980]:
w(θ) = ∫dx_1 p(x_1) ∫dx_2 p(x_2) ξ(r),
where $r$ is the 3-d separation between two galaxies, $r =
[x_1^2+x_2^2-2x_1x_2\cos{\theta}]^{1/2}$, $x_i=x(z_i)$ is the comoving
distance to a galaxy at redshift $z_i$, and $p(x)$ is the radial
selection function normalized as $\int^\infty_0dxp(x)=1$.

§.§ Correlation function estimation

We first measure the angular correlation function of [OII] emitters at
$z=1.19$ and $z=1.47$ from the four fields, COSMOS, DEEP2-3, ELAIS-N1,
and SXDS+XMM-LSS.  We then combine the four correlation functions at
each redshift for the statistical analysis.

We adopt the Landy-Szalay estimator [Landy & Szalay, 1993] to measure the
angular correlation function,
where $DD_k$, $RR_k$, and $DR_k$ are respectively the normalized
counts of data-data, random-random, and data-random pairs at given
angular separation $\theta$.
The random catalogs contain large numbers of points, which are more than 600 times as dense as our [OII] emitter catalogs to ensure that the shot noise due to the finite random points is negligible.
The subscript $k$ denotes the $k$th
field $(1\leq k \leq 4)$ and hat means an estimated quantity which can
be different from the true one [see equation (<ref>)
below].  To combine the correlation function measurements from the
four fields, we weigh each measurement by a function $W_k$,
1+ŵ(θ) = ∑_k=1^4W_k^2(θ) ( 1+ŵ_k(θ)) /∑_k=1^4W_k^2(θ).
Thus the combined correlation function is expressed as
ŵ(θ) = ∑_k=1^4W_k^2(θ) ŵ_k(θ)/∑_k=1^4W_k^2(θ).
We adopt the inverse-variance weighting,
$W_k(\theta)=\sigma_k^{-1}(\theta)$, where $\sigma_k^2(\theta)$ is the
diagonal component of the covariance matrix in the $k$-th field (see
next subsection).

Angular correlation functions of [OII] emitters measured from
the entire four fields at $z=1.19$ (top) and $z=1.47$ (bottom). Note
that the vertical axis mixes logarithmic and linear scalings.  The
red points are the measurement without the integral constraint
correction, $\hat{w}$.  The errorbars are estimated from jackknife
resampling.  The blue solid curve is the best-fitting model of the
power-law model, $(1-f_{\rm
fake})^2(w(\theta;A_w,\beta)-w_\Omega(A_w,\beta))$, with
$(A_w,\beta,f_{\rm fake})=(1.10,0.534,0.140)$ for $z=1.19$ and
$(1.17,0.538,0.139)$ for $z=1.47$.  The data enclosed by the two
blue vertical lines are used for this analysis.  The red solid curve
is the best-fitting model of the linearly-biased dark matter model,
$(1-f_{\rm fake})^2(b^2w_m(\theta)-w_\Omega(b))$, with $(b,f_{\rm
fake})=(1.60, 0.140)$ for $z=1.19$ and $(2.08,0.140)$ for
$z=1.47$ when the data enclosed by the two red vertical lines are used.  The red dashed curve is the same model but without
subtracting the integral constraint denoted by IC in the legend,
$(1-f_{\rm fake})^2b^2w_m(\theta)$.  The non-linear dark matter
correlation function, $w_m(\theta)$, is shown by the red dashed
curve as a reference.

The correlation function measured above is underestimated by a
constant due to the finite survey area.  The effect is known as the
integral constraint [Peebles & Groth, 1976] and calculated as

w_Ω= 1/Ω∫_ΩdΩ_1dΩ_2 w(θ),

where the integral is performed over the solid angle of the survey,
$\Omega$.  The integral constraint is calculated by a given model of
the correlation function $w(\theta;\Theta)$, where $\Theta$ is a set
of model parameters (see sections <ref> and <ref>).

Note also that the clustering amplitude is further reduced compared to
the true clustering of [OII] emitters due to the contamination of
non-[OII] emitters and noise lines. As seen in section
<ref>, we assume that the correlation of non-[OII] emitters
is negligible.  Then the estimated correlation function is related to
the true correlation function $w$ as

ŵ(θ) =  (1-f_fake)^2(w(θ) - w_Ω),

where $f_{\rm fake}$ is the fraction of non-[OII] emitters in our
sample, $0.1 \lesssim f_{\rm fake} \lesssim 0.2$ (see section <ref>).

The measured correlation function, combined over the four fields,
$\hat{w}$, is shown as the red points in the upper and lower panels of
Figure <ref> for $z=1.19$ and $z=1.47$, respectively.  The
correlation function becomes negative at scales, $\theta \sim 1 $ deg ($r\sim 50h^{-1}$Mpc),
due to the integral constraint $w_\Omega$. For the model parameter
fitting below, this effect is properly taken into account.  In this
figure we show various lines, which are the predictions based on the
power-law model and the linearly-biased dark matter model.  We discuss
these in detail in section <ref>.

The measurements of the angular correlation functions for individual
fields are presented for a consistency check in appendix

§.§ Covariance matrix

We use the Jackknife resampling technique to estimate the covariance
error matrix (see, e.g., [Lupton, 1993]).  The covariance of the
angular correlation function for $k$-th field, $C_{k,ij}\equiv C\left(
\hat{w}_k(\theta_i), \hat{w}_k(\theta_j)\right)$, is then estimated as

C_k,ij    =

where $\hat{w}^\ell_k(\theta_i)$ is the correlation function in the
$i$-th angular bin in $\ell$-th jackknifed realization, $N_{\rm sub}$
is the number of realizations in the $k$-th field and
$\bar{\hat{w}}_k(\theta_i)=N_{\rm sub}^{-1}\sum^{N_{\rm
sub}}_{\ell=1}\hat{w}^\ell_k(\theta_i)$.  The values of $N_{\rm
sub}$ for each field are shown in table <ref>.  In
total, we use 115 and 119 jackknifed realizations for the $z=1.19$ and
$1.47$ data, respectively.
We use the $k$-means algorithm
[ https://github.com/esheldon/kmeans_radec ] to divide our
survey regions to subregions with the same angular areas.
Each subsample includes a region continuous on the sky, $0.38~{\rm deg}$ on a side, which corresponds to the comoving size of $17h^{-1}$Mpc at $z=1.19$ and $20h^{-1}$Mpc at $z=1.47$.

Just like equation (<ref>), the covariance matrixes from the four
fields can be combined to $C_{ij}\equiv C\left( w(\theta_i),
w(\theta_j)\right)$ as

C_ij = ∑_k=1^4 W_k^2(θ_i)W_k^2(θ_j)C_k,ij/∑_k=1^4W_k^2(θ_i)∑_k=1^4W_k^2(θ_j) .

The square roots of the diagonal components of the covariance matrix,
$C_{ii}^{1/2}$, are shown as the error bars in figure <ref>.

§ RESULTS

§.§ Setup

In this section we consider two simple models for the angular
correlation function, the power-law and linearly biased dark matter
models in sections <ref> and <ref>, to
model the clustering amplitude of [OII] emitters.  The HOD modeling is
performed in section <ref>.

We adopt the small-angle approximation, known as the Limber's
approximation [Limber, 1954, Peebles, 1980].  With this, the
relation between the angular and spatial correlation functions given
by equation (<ref>) can be simplified to:

w(θ) = 2∫^∞_0 dx p^2(x) ∫^∞_0 dh ξ( r=√(x^2θ^2+h^2) ),

where $x=(x_1+x_2)/2$ and $h=x_2-x_1$.  The validity of the Limber's
approximation for the angular correlation with large angle separations
has been extensively discussed by, e.g., Simon, 2007 and
Crocce et al., 2011.  The approximation breaks down on large scales
when the width of the redshift distribution, $z_{\rm min}\leq z \leq
z_{\rm max}$, is too narrow.  The width of our NB filters is $\Delta x
\equiv x_{\rm max} - x_{\rm min} \simeq 50h^{-1}{\rm Mpc}$, where
$x_{\rm min} = x(z_{\rm min})$ and $x_{\rm max} = x(z_{\rm max})$, and
our analysis focuses on the angular scales of $\theta <1~{\rm deg}$.
We thus can safely rely on this approximation.

Although determination of the true radial selection function is not straightforward,
it is known to have a form similar with the transmission curve to some extent (see, e.g., figure 3 of [Hayashi et al., 2015]).
However, the
width $\Delta x \simeq 50h^{-1}{\rm Mpc}$ is narrow enough that we can
simply assume the radial selection function to be a constant,

p(x)= {
1/Δx    for   x_min ≤x ≤x_max,

0     otherwise.

Using the transmission curves available for the NB filters[https://subarutelescope.org/Observing/Instruments/HSC/sensitivity.html], we test the effect of choices of $p(z)$ on the calculation of the angular correlation functions.
We compare the angular correlation function computed using $p(z)$ from the transmission curve with that from the constant $p(z)$ in equation (<ref>). We then find that the correlation functions $w(\theta)$ have exactly the same shape with the amplitude different by $\leq 7\%$ for both NB816 and NB921. Thus, constraints on the parameters for the correlation function are not affected by the choice of $p(z)$ but it changes only the normalization, $f_{\rm fake}$, by at most $4\%$, much smaller than our $1-\sigma$ prior on $f_{\rm fake}$, $\sim 40\%$ (see the following paragraphs of this subsection). In the following analysis, we thus adopt the simple constant radial selection function [equation (<ref>)].

Given a model for the correlation function with a set of parameters
$\Theta$, $w(\theta; \Theta)$, the integral constraint
[equation (<ref>)] can be estimated as

w_Ω(Θ) = ∑_iw(θ_i;Θ)RR(θ_i)/∑_iRR(θ_i),

where $\theta_i$ is the $i$th separation bin and the sum is calculated
over the entire bins.  The normalized random-random pair count for the
combined field, $RR$, is computed in the same way as the correlation
function in equation (<ref>),

RR(θ) = ∑_k=1^4W_k^2(θ) RR_k(θ)/∑_k=1^4W_k^2(θ).

As discussed in section <ref>, not all the lines detected
as [OII] emitters are actually [OII] lines, but some are other
emission lines or just noises.  The fraction of fake [OII] lines is
denoted as $f_{\rm fake}$.  Then the observed correlation function has
contributions not only from the auto-correlation of [OII] emitters but
also from the auto-correlation of fake lines and their
cross-correlation, which are suppressed by factors of $f_{\rm fake}^2$
and $2f_{\rm fake}(1-f_{\rm fake})$, respectively (see e.g.,
[van den Bosch et al., 2013, Okumura et al., 2015, Okumura et al., 2017] for detailed
theoretical schemes for the tracer decomposition).  The contamination
fraction is estimated to be around $f_{\rm fake}\sim 0.14$ by applying
photometric redshift and color selections to galaxies confirmed with
spectroscopic redshifts [Hayashi et al., 2020].  Thus, even though the
non-[OII] emitters have a non-negligible auto-correlation, the
(intrinsically small) clustering amplitude is suppressed by $0.14^2
\sim 0.02$.  We therefore can safely assume that the observed
correlation is dominated by the auto correlation of [OII] emitters
(see, e.g., [Okumura et al., 2016, Kashino et al., 2017] for a similar
treatment)[If $f_{\rm fake}$ becomes large so that the
contribution from non-[OII] emitters is non-negligible, a
sophisticated treatment of the line contaminant for the cosmological
analysis is required (e.g., [Leung] et al., 2017, Addison et al., 2019])..  The
amplitude of the correlation function is thus simply scaled by a
factor of $(1-f_{\rm fake})^2$.  We treat $f_{\rm fake}$ as a free
parameter with a prior of $f_{\rm fake}=0.140\pm 0.060$ to cover
the uncertainties described in section <ref> within the $1-\sigma$ confidence level.
We chose this uncertainty to be very conservative and much larger than
it is to avoid inducing biased constraints on model parameters by imposing
strong (and possibly incorrect) priors.
The imposed prior is coincidentally equivalent to that adopted in Kashino et al., 2017.

Finally, once the model with parameters for the correlation function
is given as $w(\theta;\Theta)$, the theoretical prediction to be
compared to the observed correlation function, $\hat{w}$, is obtained
with an additional parameter $f_{\rm fake}$ through equation
(<ref>).  To constrain model parameters, the $\chi^2$
statistic is calculated for the correlation function, $w$, as

χ_w^2 (Θ)= ∑_i=1^N_bin∑_j=1^N_bin Δ_i C^-1_ij Δ_j,

where $N_{\rm bin}$ is the number of angular separation bins used in
the analysis, $\Delta_i = \hat{w}_{\rm obs}(\theta_i)-\hat{w}_{\rm
th}(\theta_i;\Theta)$ is the difference between the observed
correlation function [equation (<ref>)] and the theoretical
prediction [equation (<ref>)] with $\Theta$ being a set of
parameters to be constrained, and $C^{-1}_{ij}$ is the inverse
covariance matrix [equation (<ref>)] multiplied by a
correction for the finite realization effect, $(N_{\rm sub}-N_{\rm
bin}-2)(N_{\rm sub}-1)^{-1}$ [Hartlap et al., 2007].

In order to perform a maximum likelihood analysis, we use the Markov
chain Monte Carlo (MCMC) sampler
[Foreman-Mackey et al., 2013].

§.§ Power-law constraints

It has been known that shapes of the spatial and angular correlation
functions approximately follow a power law, originally found by
Totsuji & Kihara, 1969 (see also [Peebles, 1980]), as

w(θ; A_w, β) = A_w (
θ/1 arcmin

The $\chi^2$ statistic is calculated with $\Theta=(A_w,\beta,f_{\rm
fake})$.  Because a power-law model is known to fail to fit an
observed correlation function at small- and large-scale limits, we
compute the $\chi^2$ statistic for a relatively narrow angular
separation range, $0.004 < \theta< 0.4~[{\rm deg}]$, denoted by the
blue vertical lines in figure <ref>. Thus the number of bins
is $N_{\rm bin}=10$.

Constraints on the power-law model parameters $(A_w,\beta)$
and the fake line fraction parameter $f_{\rm fake}$ from the angular
clustering of [OII] emitters for $z=1.19$ (blue) and $z=1.47$
(orange). Contours show the $68\%$ and $95\%$ confidence levels. The
diagonal panel show the posterior probability distribution of each
parameter. A Gaussian prior is adopted for the fake line fraction,
$f_{\rm fake}=0.140\pm 0.060$, which is shown as the black curve compared with the one-dimensional
posterior. The posterior of $f_{\rm fake}$ for $z=1.47$ (orange) is almost entirely behind the one for $z=1.19$.

Parameter constraints on single power-law and linearly-biased dark matter models${}^*$

2c|NB816       2cNB921

Power-law parameter      Prior     Best-fit     Posterior PDF     Best-fit     Posterior PDF

$A_w$ (at $1$ arcmin)       None     $1.10$       $1.12^{ + 0.21} _{- 0.16}$       $1.17$     $1.20^{ + 0.20}_{ - 0.17}$

$\beta$                        None            $0.534$     $0.533 ^{+ 0.057}_{ - 0.055}$     $ 0.538$     $0.537^{+ 0.044}_{ - 0.043}$

$f_{\rm fake}$    $0.140\pm0.060$     $0.140$     $0.149^{ + 0.059}_{ - 0.059}$     0.139     $0.149^{ + 0.060}_{ - 0.059}$

$r_0(h^{-1}{\rm Mpc})$     $-$           $4.08$     $4.12^{ + 0.50}_{ - 0.41}$     $ 4.55$     $4.61^{+ 0.51}_{- 0.45}$

Linearly-biased DM Parameter      Prior     Best-fit     Posterior PDF     Best-fit     Posterior PDF

$b$                     None     $1.60$     $1.61 ^{+ 0.13}_{ - 0.11}$     2.08     $2.09^{ + 0.17}_{ - 0.15}$

$f_{\rm fake}$    $0.140\pm0.060$     $0.140$     $0.145^{ + 0.059}_{ - 0.059}$     0.140     $0.145^{ + 0.060}_{ - 0.059}$

$*$ The values in the column of “Prior” quote the mean and standard
deviation of the Gaussian priors. The column of “Best-fit” shows
the parameter set which gives the minimum value of $\chi^2$. In the
column of “Posterior PDF”, the central value is a median and the
error means $16-84$ percentiles after other parameters are
marginalized over.

Figure <ref> shows the joint constraints on the
power-law parameters and the fake line fraction, $(A_w, \beta, f_{\rm
fake})$ with blue and orange contours for $z=1.19$ and $z=1.47$,
respectively[To show parameter constraints using the
MCMC sampling here and in the rest of this paper, we use a python package
[Bocquet & Carter, 2016].].

The parameter set with the minimum value of $\chi^2$
for $z=1.19$ is $(A_w, \beta, f_{\rm fake})=(1.10, 0.534, 0.140)$
while that for $z=1.47$ is $(A_w, \beta, f_{\rm fake})= (1.17, 0.538,
0.139)$.  With these sets of the best-fitting parameters, we derive
the integral constraints as $w_\Omega = 0.110$ and $w_\Omega = 0.114$,
for lower and higher redshifts, respectively.  The best-fitting
power-law model prediction for the measured correlation function,
$(1-f_{\rm fake})^2(A_w\theta^{-\beta}-w_\Omega)$, is shown by the
blue line in figure <ref>.  The constrained parameters are
shown in table <ref>.

When the two parameters for a power-law model of the angular
correlation, $(A_w,\beta)$, are known, one can determine its
3-dimensional clustering, $\xi(r) = (r/r_0)^{-\gamma}$, as
[Peebles, 1980, Efstathiou et al., 1991, Simon, 2007]

A_w    =    r_0^γ( 10800/π )^γ-1
B( 1/2,γ-1/2 )∫^∞_0 dx p^2(x) x^1-γ

=   r_0^γ10800^γ-1(x_max^γ-1-x_min^γ-1)/π^γ-1(2-γ)(x_max-x_min)^2
B( 1/2,γ-1/2 ),

β   =    γ-1,

where $B$ is the beta function.  The final expression of $A_w$ is
derived using the constant radial selection function [equation
(<ref>)].  As a result, the clustering length is
constrained to be $r_0= 4.12^{+0.50}_{-0.41}$ at $z=1.19$ and $r_0=
4.61^{+0.51}_{-0.45}$ at $z=1.47$. These values correspond to the host
halo masses of $M\gtrsim 10^{12} M_\odot /h$ and are roughly
consistent with the preceding results of ELG clustering (see e.g.,
[Kashino et al., 2017]).  Here we do not make a detailed comparison with
previous works nor argue the clustering evolution further because the
power-law model is not so accurate to be used to interpret the precise
correlation function measurement anyways.  In section <ref>,
we discuss physical properties of halos hosting ELGs based on an HOD
model that is much more sophisticated than the power-law model.

§.§ Linearly biased dark matter model

Next, we consider another simple model, the non-linear correlation
function with the linear galaxy bias factor [Kaiser, 1984],

where $\xi_m$ is the underlying matter correlation function.  We use
the fitting formula of Takahashi et al., 2012 to compute the nonlinear
matter power spectrum, $P_m(k)$, and then Fourier transform it to obtain $\xi_m$.
Given equation (<ref>), a similar form is derived for
the angular correlation function as

where $w_m$ is the angular correlation function of matter, calculated
using equation (<ref>) with $\xi(r)$ replaced by the
3-dimensional mater correlation function, $\xi_m(r)$.  In this model
there are two free parameters to be determined, $\Theta=(b,f_{\rm
fake})$.  This model is valid for a wider range of angular scales
than the single power-law model.  The data at the scales $0.0025 <
\theta < 1~[{\rm deg}]$ are used to calculate the $\chi^2$ statistic
($N_{\rm bin}=13$).

The joint constraints on $(b,f_{\rm fake})$ are shown in figure
<ref>.  We find the best-fitting value of the bias
parameter, $b= 1.61 ^{+ 0.13}_{ - 0.11}$ and $b= 2.09^{ + 0.17}_{ -
0.15}$, respectively for the lower and higher redshift measurements.
These bias values can be compared to the effective bias $b_{\rm eff}$
to be determined by a HOD modeling [equation (<ref>)].  The
integral constraint is calculated by these best-fitting parameter sets
for the correlation in the combined field as $w_\Omega = 0.0270$ and
$w_\Omega = 0.0277$ for $z=1.19$ and $1.47$, respectively.

The best-fitting linearly-biased dark matter model, $(1-f_{\rm
fake})^2(b^2w_m(\theta)-w_\Omega)$, is shown as the red solid curve
in figure <ref>.  To see the contribution of the integral
constraint, we show the same model without $w_\Omega$ term, $(1-f_{\rm
fake})^2b^2w_m(\theta)$, as the red dashed curve. As expected, the
model correlation function is positive over all the scales studied

As a reference, the dark matter correlation function, $w_m$, is shown
by the black dashed curve in figure <ref>.

Constraints on the linear galaxy bias $b$ of the nonlinear
dark matter model and the fake line fraction parameter $f_{\rm
fake}$ for $z=1.19$ (blue) and $z=1.47$ (orange). Contours show
the $68\%$ and $95\%$ confidence levels. The diagonal panel show the
posterior probability distribution of each parameter. A Gaussian
prior is adopted for the fake line fraction, $f_{\rm fake}=0.140\pm
0.060$, which is shown as the black curve compared with the one-dimensional
posterior. The posterior of $f_{\rm fake}$ for $z=1.47$ (orange) is almost entirely behind the one for $z=1.19$.

The relatively large uncertainties on the bias parameter comes from
our conservative prior on $f_{\rm fake}$.
To see this, let us directly constrain the observed amplitude by
considering the relation:
ŵ(θ;b)    =    w(θ;b)-w_Ω(b)

=    w(θ;b) - ∑_iw(θ_i;b)RR(θ_i)/∑_iRR(θ_i),
where $\wt{b}$ is the bias parameter of the observed field, $\wt{w}(\theta,\wt{b}) = \wt{b}^2w_m(\theta)$.
Equation (<ref>) becomes equivalent to equations (<ref>) and (<ref>)
if we set $\wt{b} \equiv (1-f_{\rm fake})b$.
By adopting this 1-parameter linearly-biased dark matter model, we obtain
tighter constraints as $\wt{b} = 1.379^{+0.042}_{-0.043}$ and $\wt{b} = 1.786 ^{+0.052}_{-0.054}$ at $z=1.19$ and $z=1.47$, respectively, without any prior.

In the next subsection we further analyze the [OII] emitter samples with different luminosity and stellar mass thresholds using the linearly-biased DM model.

Angular correlation functions of [OII] emitters at $z=1.19$ (top) and $z=1.47$ (bottom) with different peak luminosity cuts. The unit of $L$ in the legend is $[{\rm erg~s}^{-1}]$. For the visualization purpose, the integral constraint is taken into account, and thus the vertical axis presents $\hat{w}+(1-f_{\rm fake})^2w_\Omega$. The amplitudes of the correlation functions shown by the blue, green, yellow, brown and purple points are multiplied by $2^n$ where $n=1,2,\cdots ,5$, respectively. The dashed curves are the best-fitting models of the linearly-biased dark matter model, $(1-f_{\rm fake})^2b^2 w_m(\theta)$.

Constraints on the measured clustering amplitude of [OII] emitters, $\wt{b}=(1-f_{\rm fake})b$, as a function of the minimum peak luminosity $L_{\rm min}$ (left) and stellar mass $M_{*,{\rm min}}$ (right). The blue and orange points are the results at $z=1.19$ and $z=1.47$, respectively.

§.§ Line luminosity and stellar mass dependences of [OII]-emitter clustering

It is interesting to see how the amplitude of a measured correlation function
depends on properties of  the [OII] emitters.
Here we consider two basic properties, the luminosity of the [OII] emission line and the stellar mass of the
host galaxies, and split our sample into subsamples based on these thresholds.

Figure <ref> shows the angular correlation functions
measured with different line luminosity cuts.
For clarity, we add the integral constraint correction so that all the correlation functions
become positive at all scales probed.
For each correlation function, we perform a model fitting using the linearly biased dark matter model
considered in section <ref>.
The best-fitting model for each measurement, $(1-f_{\rm fake})^2b^2w_m(\theta)$,
is presented by the dashed curve.
As discussed in section <ref>, this simple model fails to explain the measured
correlation function at small scales.
To see the dependence of the clustering amplitude on the luminosity,
we constrain the parameter $\wt{b}$ in equation (<ref>) as a function of the minimum line luminosity $L_{\rm min}$.
The result is shown in the left panel of figure <ref>.
For each redshift, we find that the bias $b$ increases with increasing line luminosity, by $\sim 10\%$.
If the fake line fraction, $f_{\rm fake}$, changed significantly with the luminosity, this luminosity dependence would be modulated.

We repeat the same analysis by splitting our [OII] emitter sample into subsamples with
stellar mass cuts $M_{*,{\rm min}}$. We measure the angular correlation function of [OII] emitters with the stellar mass of $M_* \geq M_{*,{\rm min}}$.
Since the fitting result for $w(\theta)$ is similar to that with the luminosity cuts in figure <ref>,
we do not plot the figure.
We constrain the measured clustering amplitude $\wt{b}$ as a function of $M_{*,{\rm min}}$, as shown in the right panel of figure <ref>.
A trend similar to the line luminosity limited samples is found:
[OII] emitters with higher stellar masses cluster more strongly than those with lower masses.
Almost no mass dependence is seen in the clustering at $z=1.47$.
As discussed in section 3.6 of Hayashi et al., 2018,
multi-band data at longer wavelength, such as near-infrared data, would be
required to accurately estimate the stellar masses at $z\sim 1.5$.
Thus, our result for the stellar mass dependence of the clustering at $z=1.47$ needs to be interpreted with caution.
Overall, however, the dependence of clustering of [OII] emitters on stellar masses is not as strong as
that on line luminosities,
which is qualitatively consistent with the result of Khostovan et al., 2018.

§ DARK MATTER-[OII] EMITTER CONNECTIONS

The halo occupation distribution (HOD) modeling is an empirical and
parametric way to describe the connection between galaxies and their
host halos through the modeling of the abundance and clustering of
[Jing et al., 1998, Peacock & Smith, 2000, Seljak, 2000, Ma & Fry, 2000, Scoccimarro et al., 2001, Cooray & Sheth, 2002, Berlind & Weinberg, 2002, Zheng et al., 2005, van den Bosch et al., 2013, Hikage et al., 2013, Okumura et al., 2015].
In this section we use HOD modeling to investigate physical properties
of halos which host [OII] emitters.

§.§ Formalism of halo model

The mean number of galaxies residing in a halo of mass $M$, denoted as
$\left\langle N(M) \right\rangle$, is decomposed into central and
satellite components,

\begin{eqnarray}
&& \left\langle N(M) \right\rangle =  \left\langle N_{\rm cen}(M)\right\rangle + \left\langle N_{\rm sat}(M)\right\rangle.
\end{eqnarray}

Once the HOD model is specified, the average number density of
galaxies is calculated as

n_g=∫⟨N(M) ⟩n(M)dM,
where $n(M)dM$ is the halo mass function.

In a halo model approach, the power spectrum of galaxies, $P(k)$ where
$k$ is the wavenumber, can be decomposed into the one- and two-halo

The one-halo term can be further decomposed into contributions from
the clustering of central-satellite and satellite-satellite pairs
hosted by the same halos, namely,

P_1h(k)    =    1/n_g^2∫n(M)dM [ 2⟨N_cenN_sat ⟩u(k,M) .

. + ⟨N_sat(N_sat-1) ⟩u ^2 (k,M)],

where $u(k,M)$ describes the Fourier transform of the density profile
of satellite galaxies within dark matter halos, which is assumed to
follow the matter density profile. We assume that $N_{\rm sat}$
follows the Poisson statistics such that $\langle N_{\rm sat}(N_{\rm
sat}-1) \rangle = \langle N_{\rm sat}\rangle^2$.  We also assume
that the occupation numbers of central and satellite galaxies are
independent, i.e., $\langle N_{\rm cen} N_{\rm sat}\rangle = \langle
N_{\rm cen}\rangle \langle N_{\rm sat}\rangle$.  The two-halo term can
be similarly decomposed into central-central, central-satellite and
satellite-satellite galaxy pairs hosted by distinct halos,

P_2h(k)    =    1/n_g^2∫n(M)dM ∫n(M')dM'

×[ ⟨N_cen ⟩+ ⟨N_sat⟩u(k,M)]

×[ ⟨N_cen ⟩+ ⟨N_sat⟩u(k,M')]

where $P_{hh}(k;M,M')$ is the cross-power spectrum of halos with
masses $M$ and $M'$.
The halo mass function is calculated using the
model derived by Tinker et al., 2010.  The density profile of halos is
assumed to have the NFW profile [Navarro et al., 1997] with the
concentration given by Duffy et al., 2008.
We compute $P_{hh}$ as a product of the nonlinear matter power spectrum [Takahashi et al., 2012] and the halo bias,
P_hh(k;M,M') = b_h(M)b_h(M')P_m(k),
where $b_h(M)$ is the bias for a halo of mass $M$.
We use the large-scale halo bias $b_h(M)$ proposed also by Tinker et al., 2010.

The three-dimensional correlation function of galaxies is calculated
by a Fourier transform of the power spectrum,

ξ(r) = 1/2π^2 ∫^∞_0 dkk^2P(k) sinkr/kr.
Then the angular correlation function is finally obtained through
equation (<ref>).

Constraints on parameters of Geach HOD model, $(\log{M_{\rm
c}}, \log{M_{\rm min}}, \alpha, \sigma_{\log{M}}, F_s,
F_c^A,F_c^B) $ and $f_{\rm fake}$ for $z=1.19$ (blue) and $z=1.47$
(orange). Two-dimensional contours show the $68\%$ and $95\%$
confidence levels after the other six parameters marginalized
over. The diagonal panels show the posterior probability distribution
of each parameter. Gaussian priors are assumed for $\alpha$ and
$f_{\rm fake}$, as depicted by the black sold curves in the 1-d
posterior panels.

In addition to the number density, given a set of HOD parameters, one
can determine various physical quantities such as the effective bias,

b_eff = 1/n_g∫b_h(M) ⟨N(M) ⟩n(M)dM,
the effective halo mass,

M_eff = 1/n_g∫M⟨N(M) ⟩n(M) dM,

and the satellite fraction,

f_sat = 1/n_g ∫⟨N_sat(M) ⟩n(M) dM.

Priors and constraints of the HOD parameters for Geach model${}^*$

2c|NB816       2cNB921

Parameter      Prior     Best-fit     Posterior PDF     Best-fit     Posterior PDF

$\log{M_{\rm c}} / (h^{-1}M_\odot)$          None                   $11.75$     $12.04^{+0.57}_{-0.32}$         $11.93$     $11.91^{+ 0.19} _{- 0.18}$

$\log{M_{\rm min}} / (h^{-1}M_\odot)$     None                    $12.46$     $12.61^{+0.32}_{-0.43}$        $12.47$     $12.57^{+0.29}_{-0.35}$

$\sigma_{\log{M}}$                                  $[0,1]$                 $0.06$     $0.40^{+0.28}_{-0.31}$            $0.13$     $0.17^{+ 0.15}_{- 0.10}$

$\alpha$                                                   $1.00\pm 0.20$      $1.06$     $1.03^{+0.12}_{-0.13}$             $1.23$     $1.12^{+ 0.13}_{- 0.12}$

$F_c^A$                                                   $[0, 0.5]$             $0.13$     $0.26^{+0.16}_{-0.12}$             0.14     $0.26^{+ 0.15}_{- 0.13}$

$F_c^B$                                                   $[0,1]$             $0.95$     $0.37^{+0.39}_{-0.28}$            0.90     $0.53^{+ 0.32}_{- 0.34}$

$F_s$                                                      $[0, 1]$                  $0.98$     $0.54^{+0.31}_{-0.34}$             0.73     $0.55^{+ 0.31}_{- 0.33}$

$f_{\rm fake}$                                          $0.140\pm 0.060$     $0.172$     $0.140^{+0.048}_{-0.053}$     0.128     $0.104^{+ 0.043}_{- 0.041}$

$\chi^2/\nu$    ($\nu=7$)         2.17        0.856

Inferred quantity     Measurement     Best-fit     Posterior PDF     Best-fit     Posterior PDF

$\log{n_g / (h^{-1}{\rm Mpc})^{-3}}$    $-2.259\pm 0.068$ (NB816)     $-2.313 $     $-2.316^{+0.071}_{-0.073}$   $-2.401$     $-2.381^{+0.074}_{-0.073}$

$-2.344\pm 0.070$ (NB921)

$f_{\rm sat}$     1c |$\cdots$                                                 $0.308$     $0.159^{+0.120}_{-0.047}$     $0.242$     $0.159^{ + 0.109}_{ - 0.049}$

$b_{\rm eff}$     1c |$\cdots$                                                 $1.797$     $1.700^{+0.084}_{-0.111}$     $2.029$     $1.981^{+0.072}_{-0.068}$

$\log{M_{\rm eff}} / (h^{-1}M_\odot) $     1c |$\cdots$           $12.831$     $12.703^{+0.093}_{-0.071}$     $12.699$     $12.609^{+0.085}_{-0.051}$

$*$ In the “Prior” column the ranges specified in brackets are for
uniform priors while the others we quote the mean and standard
deviation of the Gaussian priors. The column of “Best-fit” shows
the parameter set which gives the minimum value of $\chi^2$. In the
column of “Posterior PDF”, the central value is a median and the
error means $16-84$ percentiles after other parameters are
marginalized over.
The measured number density $\log{n_g}$ includes non-[OII] emitters,
and thus its best-fitting values differ from the measure ones by $\log{(1-f_{\rm fake})} \sim -0.06$.

HOD with the best-fitting parameters of Geach model at
$z=1.19$ (top) and $z=1.47$ (bottom).  The blue and yellow lines
show the average numbers of centrals from the Gaussian and error
function terms (see equation (<ref>)), while the
green line shows that of satellites. The red solid line is the sum
of the blue, yellow and green lines, the average number of all the
[OII] emitters. The dark and light red shaded regions indicate the
$68\%$ and $95\%$ confidence intervals.  The red dashed curve is the
same as the red solid curve but for the best-fitting HOD of
Zheng et al., 2005.  For comparison, the black curve shows the HOD
of the Geach model which gives the minimum $\chi^2$ value.

The angular correlation functions of [OII] emitters at $z=1.19$ (top) and
$z=1.47$ (bottom).  Note that the vertical axis mixes logarithmic
and linear scalings.  The red points are the measurement, $\hat{w}$,
the same as those in figure <ref>.  The blue and green
solid curves are the one- and two-halo terms with the best-fitting
Geach HOD model, respectively, and the red solid curve is their
sum. For comparison, we also show the best-fitting Zheng HOD model
as the red dashed curve.  The data enclosed by the two red vertical
lines are used for this HOD analysis.

§.§ Halo occupation distribution model for [OII] emitters

In this paper we adopt the HOD model developed by Geach et al., 2012
to describe the population of ELGs (see also
[Kim et al., 2011, Contreras et al., 2013]). In this model, the central HOD is
described by two components:

\begin{eqnarray}
\left\langle N_{\rm cen}(M)\right\rangle &=& F_c^B(1-F_c^A) \exp{\left[ -\frac{\log{(M/M_{\rm c})}^2}{2(\sigma_{\log{M}})^2} \right]} \nn \\
&& + F_c^A \left [1+{\rm erf} \left(\frac{\log{(M/M_{\rm c})}}{\sigma_{\log{M}}}\right) \right], \label{eq:hod_geach_cen}
\eey
where $F_c^{A,B}$ are normalization factors. The first component
describes the Gaussian distribution of centrals around halos of
average mass $M_{\rm c}$ with the dispersion $\sigma_{\log{M}}$, and
the second component describes the standard mass-limited step function
form proposed by Zheng} {et~al.}, 2005 (see section \ref{sec:hod_zheng}
below).  The satellite HOD is given by
\be
\left\langle N_{\rm sat}(M) \right\rangle= F_s \left [1+{\rm erf} \left(\frac{\log{(M/M_{\rm min})}}{\delta_{\log{M}}}\right) \right] \left( \frac{M}{M_{\rm min}}\right)^\alpha, \label{eq:hod_geach_sat}
\ee
where $F_s$ is the mean number of satellites per halo at the
transition mass $M_{\rm min}$ which corresponds to the characteristic
mass above which halos can contain satellites, $\delta_{\log{M}}$ is
the width of the transition from zero satellites per halo to the power
law, and $\alpha$ is the slope of the power law which gives the mean
number of satellites for $M>M_{\rm min}$.  We refer the reader to
Geach} {et~al.}, 2012 and Contreras} {et~al.}, 2013 for a more detailed
explanation of this model.

We fix the least-important parameter, $\delta_{\log{M}}$, to
$\delta_{\log{M}}=1$.  Thus the number of free parameters in the Geach HOD
model together with the fake line fraction is eight, $\Theta = (M_{\rm
c}, M_{\rm min}, \sigma_{\log{M}}, \alpha, F_c^A, F_c^B, F_s, f_{\rm
fake})$.  The parameter $\alpha$ is known to have a value around
unity.  We thus apply a Gaussian prior on $\alpha$ as $\alpha = 1.00
\pm 0.20$.  As with the analysis in section \ref{sec:result}, we
further impose a Gaussian prior on the fake line fraction parameter as
$f_{\rm fake}= 0.140 \pm 0.060$.  Uniform priors are applied for the
other 6 parameters.  Table \ref{tab:hod_geach} summarizes the priors
on all the eight parameters.  We use the Python package
\texttt{HALOMOD} [Murray} {et~al.}, 2013] to calculate the HOD model
prediction for the angular correlation function.

\subsection{HOD parameter constraints} \label{sec:constraint}
We present constraints on the HOD model using the measured angular
correlation function of the [OII] emitters. We compute the $\chi^2$
statistic for the HOD model constructed from the observed angular
correlation function, $\chi_w^2$, and number density, $\chi_n^2$, as
\be
\chi^2 (\Theta) = \chi_w^2 (\Theta) +   \chi_n^2 (\Theta),
\ee
where $\chi_w^2$ is given by equation (\ref{eq:chi2_w}) and $\chi_n^2$
is given by
\be
\chi_n^2 (\Theta) = \frac{\left[\log{n_g^{\rm obs}}-\log{n_g^{\rm th}(\Theta)}+\log{(1-f_{\rm fake})}\right]^2}{\sigma^2_{\log{n_g}}},
\ee
where $\Theta = (\log{M_{\rm c}}, \log{M_{\rm min}}, \alpha,
\sigma_{\log{M}}, F_s, F_c^A,F_c^B,f_{\rm fake})$.  Since the observed
number of [OII] emitters includes non-[OII] contaminants characterized
by the fraction $f_{\rm fake}$, we need to take into account a factor
of $(1-f_{\rm fake})$ difference between the observed number density
and its theoretical prediction computed by equation (\ref{eq:ng_hod}).
As shown in table \ref{tab:four_fields}, the observed number density
is $n_g^{\rm obs} = 5.50 \times 10^{-3} (h/{\rm Mpc})^3$ at $z=1.19$
and $n_g^{\rm obs} = 4.08 \times 10^{-3} (h/{\rm Mpc})^3$ at $z=1.47$.
We set the uncertainty of $\log{n_g^{\rm obs}}$ as $\sigma_{\log{n_g}}
= 0.03 |\log{n_g^{\rm obs}}|$ following the measurement uncertainty of
the luminosity function [Hayashi} {et~al.}, 2020].  Due to the flexibility
of the model, the HOD enables us to fit an observed correlation
function over a broader range than the power-law and nonlinear dark
matter models considered in the previous section.  We thus use the
data at $0.0015<\theta<1~[{\rm deg}]$ and the number of angular
separation bins is $N_{\rm bin} = 14$.  Together with the number
density, the degree of freedom is $\nu = N_{\rm bin} +1 - 8 = 7$.

The resulting constraints on the HOD parameters of the Geach model are
shown in figure \ref{fig:result_hod_contreras}.  This and the following
two figures are the main results of this paper.  The eight dimensional
posterior distributions are visualized in two-dimensional contours
with the other six parameters marginalized over.  The top row of each
column shows the one-dimensional posterior.  The constraints are
summarized in table \ref{tab:hod_geach}.  The normalizations
parameters, $F_c^{A,B}$ and $F_s$, are poorly constrained, which are
expected and consistent with Geach} {et~al.}, 2012 and Hong} {et~al.}, 2019.
Overall, constraints from the $z=1.47$ sample are tighter than those
from the $z=1.19$ one.  Particularly, the marginalized distribution
for $\sigma_{\log{M}}$ for the $z=1.19$ sample is much wider and even
becomes bimodal.  Interestingly, this bimodality was seen by
Hong} {et~al.}, 2019 who analyzed the clustering of Ly$\alpha$ emitters
using the Geach HOD model.  Although there are many free parameters
with few priors, we obtain meaningful constraints due to
our large sample. Particularly
strong constrains are obtained for the two mass parameters, $M_{\rm
c}$ and $M_{\rm min}$ as seen in table \ref{tab:hod_geach}, because
$M_c$ can be determined by the observed number density, $n_g$, and
$M_{\rm min}$ primarily degenerates with $M_{\rm c}$.  We do not see a
strong evolution of HOD from $z=1.47$ to $1.19$, which is consistent
with our ongoing work based on cosmological hydrodynamical simulations
[Osato} \& {Okumura}, 2021].

The posterior distribution of the HOD is shown in figure
\ref{fig:result_hodfig_contreras}.  The dark and light shaded regions
show the 68 and 95\% confidence intervals for the total HOD,
respectively, and red solid curve is the median. The black solid curve
is the HOD with a set of parameters which give the minimum $\chi^2$
value in the eight dimensional parameter space. The HOD with the
minimum $\chi^2$ looks somewhat different from the median of the HOD,
particularly the central HOD. This is a similar trend with the finding
of Hong} {et~al.}, 2019. They found that two independent algorithms give
very different best-fitting HOD parameters (See Model\#1 and Model\#2
of their figure 7). This small but non-negligible discrepancy largely
comes from the fact that the normalization parameters, especially
$F_c^B$, and the scatter of the mass, $\sigma_{\log{M}}$, are poorly
constrained.  Even though the HODs are different, the predicted
angular correlation functions with these HOD parameter sets become
very similar.

The best-fitting Geach HOD model prediction, $(1-f_{\rm
fake})^2(w(\theta;\Theta)-w_\Omega)$, is shown as the red solid
curve in figure \ref{fig:result_xi_contreras}. The contributions from
the 1-halo and 2-halo components are shown as the blue and green
curves, respectively.  Obviously, agreement of the measured
correlation function with the HOD model prediction is much more
remarkable than that with the power-law or linearly biased dark matter
model, compared to figure \ref{fig:w_all}.

% Section 5.4
\subsection{Derived physical parameters for host halos}

Using the constrained HOD parameters, one can infer the parameters
which characterize the properties of [OII] emitters and their host
dark matter halos, such as the galaxy number density, effective halo
bias and mass, and satellite galaxy fraction, through equations
(\ref{eq:ng_hod}), (\ref{eq:b_hod}), (\ref{eq:m_hod}), and
(\ref{eq:fsat_hod}), respectively.  The posterior distribution for
these parameters is shown in figure \ref{fig:result_params_b_m_n_f}
and summarized in table \ref{tab:hod_geach}.
% Number density
The number density is constrained following the imposed prior. Note
again that the constrained number density is different from the
observed number density by a factor of $(1-f_{\rm fake})$, $n_g^{\rm
th} = (1-f_{\rm fake})n_g^{\rm obs}$.
% bias
The effective bias parameter is determined at $z=1.19$ and $1.47$ as
$b_{\rm eff}=1.701^{+0.083}_{-0.110}$ and $b_{\rm
eff}=1.981^{+0.072}_{-0.068}$, respectively.  They are fully
consistent with the linear bias parameter from much simpler analysis
in section \ref{sec:nl_matter}, $b=1.61^{+0.13}_{-0.11}$ and

% Figure 10
\begin{figure}%[bt]
\begin{center}
\vspace{-.2cm}
\hspace{-0.2cm}
\FigureFile(20.0mm,45mm){result_params_NB0816_beff_cont_v7_29998_60_zhe_v6_12990_60.pdf}
\hspace{-0.05cm}
\FigureFile(20.0mm,45mm){result_params_NB0816_Meff_cont_v7_29998_60_zhe_v6_12990_60.pdf}
\hspace{-0.05cm}
\FigureFile(20.0mm,45mm){result_params_NB0816_ng_cont_v7_29998_60_zhe_v6_12990_60.pdf}
\hspace{-0.18cm}
\vspace{-.35cm}
\FigureFile(21.2mm,45mm){result_params_NB0816_fsat_cont_v7_29998_60_zhe_v6_12990_60.pdf}
\\
\hspace{-0.2cm}
\FigureFile(20.0mm,45mm){result_params_NB0921_beff_cont_v7_28019_60_zhe_v6_12997_60.pdf}
\hspace{-0.05cm}
\FigureFile(20.0mm,45mm){result_params_NB0921_Meff_cont_v7_28019_60_zhe_v6_12997_60.pdf}
\hspace{-0.05cm}
\FigureFile(20.mm,45mm){result_params_NB0921_ng_cont_v7_28019_60_zhe_v6_12997_60.pdf}
\hspace{-0.18cm}
\FigureFile(21.2mm,45mm){result_params_NB0921_fsat_cont_v7_28019_60_zhe_v6_12997_60.pdf}
\vspace{-0.3cm}
\end{center}
\caption{Posterior distribution for the parameters derived from the
best-fitting parameters for the Geach (blue) and Zheng (orange) HOD
models.  From the left to right, the posteriors for the effective
bias, effective halo mass, the galaxy number density and the
satellite fraction are shown. The upper and lower panels show
results for $z=1.19$ and $z=1.47$, respectively.
The number density presented here is different from the observed number density
shown in table \ref{tab:four_fields} by a factor of $(1-f_{\rm fake})$, $(1-f_{\rm fake})n_g^{\rm obs}$
(see table \ref{tab:hod_geach} for the constrained values for $f_{\rm fake}$).
Gaussian priors assumed for $\log{n_g}$ are depicted by the black solid curves, where the best-fitting value of $f_{\rm fake}$ is adopted to multiply by $(1-f_{\rm fake})$.
The posterior of $\log{n_g}$ for $z = 1.47$ (orange) is almost entirely behind the one for $z = 1.19$.
\label{fig:result_params_b_m_n_f}
\end{figure}

%Effective halo mass
The effective masses of halos hosting [OII] emitters are derived to be
$\log{M_{\rm eff}/(h^{-1}M_\odot)}=12.703^{+0.091}_{-0.069}$ and
$12.609^{+0.085}_{-0.051}$ at $z=1.19$ and $1.47$, respectively.
Figure \ref{fig:mass_accretion} shows the constraints on $M_{\rm eff}$
as a function of redshift obtained from our HSC [OII] emitters
together with the previous studies at $0<z<4$ which were already
presented in Kashino} {et~al.}, 2017.  The constraints include the
studies of photo-$z$ galaxies from the CFHT Legacy Survey at $z\sim 0.3$,
$0.5$ and $0.7$ [Coupon} {et~al.}, 2012], VIMOS-VLT Deep Survey
(VVDS) at $z\sim 0.55$ and $1$ [Abbas} {et~al.}, 2010], NEWFIRM Medium Band
Survey at $z\sim 1.1$ and $1.5$ [Wake} {et~al.}, 2011], H$\alpha$ emitters
from the FMOS-COSMOS Survey at $z\sim 1.6$ [Kashino} {et~al.}, 2017] and
from the HiZELS survey at $z\sim 2.2$ [Geach} {et~al.}, 2012], and VUDS at
$z\sim 2.5$ and $3.5$ [Durkalec} {et~al.}, 2015].  These samples have
number densities and stellar masses similar to our [OII] emitter
Here we also plot the average mass assembly history of halos with
different present-day masses, as derived by Behroozi} {et~al.}, 2013
based on an $N$-body simulation.  According to this prediction, our
measurements of $M_{\rm eff}$ at two redshifts are well explained by
the mass-assembly history with $M(z=0)=1.5\times 10^{13}M_\odot / h$
as depicted by the blue solid curve.  The effective halo mass of [OII]
emitters at $z=1.47$ is slightly higher than the blue curve, which
reflects the fact that [OII] emitters of our $z=1.47$ sample are
intrinsically brighter than those of the lower-z sample (see section
\ref{sec:data_nb}).
Given the accurate constraints on HOD, one can in principle infer the
mass of the host halo for an individual emission line galaxy
[Oguri} \& {Lin}, 2015]. It is, however, beyond the scope of this paper and
will be investigated in future work.

%Satellite galaxy fractions
For our [OII] emitter sample, the satellite galaxy fraction is
constrained as $f_{\rm sat}=0.158^{+0.114}_{-0.047}$ and
$0.159^{+0.109}_{-0.049}$ at $z=1.19$ and $1.47$, respectively.  The
preceding studies revealed that the satellite galaxy fraction of a
given population would depend on the redshift, number density and
stellar mass [Coupon} {et~al.}, 2012, Guo} {et~al.}, 2014, Guo} {et~al.}, 2019].  At least over the
redshift we studied here, the satellite fraction of [OII] emitters in
our sample does not significantly evolve with redshift.
The analysis of Favole} {et~al.}, 2016 constrained the satellite fraction
of ELGs with the host halo mass $M\sim 10^{12}h^{-1}M_\odot$ as
$f_{\rm sat}=0.225\pm 0.025$.
Avila} {et~al.}, 2020 constrained the value of $f_{\rm sat}$ for ELGs
from the eBOSS survey, and the best fitting value of their baseline
model is $f_{\rm sat} = 0.22$, consistent with our
measurements. However, their constraints are significantly
model-dependent and varying the baseline model changes the
best-fitting value between $0.18\leq f_{\rm sat}\leq 0.70$.
Moreover, the selection of ELGs in the eBOSS survey is quite
different from ours based on the NB. We thus
cannot make a quantitative comparison with their result.
Guo} {et~al.}, 2019 further studied the relation among the satellite
fraction, stellar mass and halo mass for ELGs from the eBOSS survey.
They showed that the host halo mass of ELGs with $f_{\rm sat} \sim
0.15$ is $\log{M/(h^{-1}M_\odot)}\sim 12.6$ is a strong function of
the stellar mass. It provides good agreement with our clustering
measurements of [OII] emitters from the HSC survey.

% Figure 11
\begin{figure}[bt]
\begin{center}
\vspace{-.3cm}
\FigureFile(82mm,82mm){mass_accretion.pdf}
\end{center}
\caption{Average mass of halos hosting ELGs as a function of redshift.
The blue stars are the result from our analysis of [OII] emitters at
$z=1.19$ and $z=1.47$.  The other points are the results of the
previous HOD analysis by Abbas} {et~al.}, 2010, Wake} {et~al.}, 2011,
Geach} {et~al.}, 2012, Coupon} {et~al.}, 2012, Durkalec} {et~al.}, 2015, and
Kashino} {et~al.}, 2017 (see the legend).  The solid curves show a
prediction of mass accretion histories for halos with
$M(z=0)=10^{15},10^{14},10^{13}$ and $10^{12}M_\odot$ from the top
to bottom [Behroozi} {et~al.}, 2013]. The blue curve represents the one
with $M(z=0)=1.5\times 10^{13}~[ h^{-1}M_\odot]$.  }
\label{fig:mass_accretion}
\end{figure}

% Figure 12
\begin{figure*}%[bt]
\begin{center}
\vspace{-.5cm}
\FigureFile(125mm,125mm){result_hod_zheng_v6_12990_60_v6_12997_60.pdf}
\end{center}
\caption{Constraints on HOD parameters of Zheng's model, $(\log{M_{\rm
c}}, \log{M_{\rm min}}, \alpha, \sigma_{\log{M}}) $ and $f_{\rm
fake}$ for $z=1.19$ (blue) and $z=1.47$ (orange). Contours show
the $68\%$ and $95\%$ confidence levels. The diagonal panels show the
posterior probability distribution of each parameter. Gaussian
priors are assumed for $\alpha$ and $f_{\rm fake}$, as depicted by
the black curves in the panels of the 1-d posterior distributions.
\label{fig:result_hod_zheng}
\end{figure*}

\subsection{Comparison to a simpler HOD}\label{sec:hod_zheng}
While the Geach HOD model is general and flexible, as a trade-off, the
HOD parameters are constrained more poorly and more degenerated with
each other.
Furthermore, it is not clear whether different halo occupation functions provide the same physical interpretation for a measured correlation function.
A similar attempt has been made by studying one of the simplest HOD model
[Zehavi} {et~al.}, 2005] in Hong} {et~al.}, 2019 (see their Appendix B).  Here
we investigate whether the predicted correlation function, halo
occupation function, and the host halo parameters of [OII] emitters
can be different between different HOD models.

To see this, we consider the commonly-adopted model of
Zheng} {et~al.}, 2005 (see also [Zheng} {et~al.}, 2007]).  The model contains five parameters in total and
the central and satellite HODs are, respectively, given by
\begin{eqnarray}
&& N_{\rm cen}(M) =
\frac{1}{2}\left[ 1+ {\rm erf}\left(\frac{\log{M / M_{\rm min}}}
{\sigma_{\log{M}}}\right) \right],  \\
&& N_{\rm sat}(M) =
\frac{1}{2}\left[ 1+ {\rm erf}\left(\frac{\log{M / M_{\rm min}}}
{\sigma_{\log{M}}}\right) \right]
\left(\frac{M-M_0}{M_1}\right)^{\alpha},
\end{eqnarray}
where $M_{\rm min}$ is the characteristic mass to host a central
galaxy, $M_1$ is a mass for a halo with a central galaxy to host one
satellite, $M_0$ is the mass scale to truncate satellites,
$\sigma_{\log{M}}$ is the characteristic transition width, and
$\alpha$ is the slope of the power law for the satellite HOD, the same
as $\alpha$ in the Geach model.  The parameter $M_0$ is a
poorly-constrained parameter.  In order to make this analysis simpler,
we impose a relation between $M_0$ and $M_1$ following
Conroy} {et~al.}, 2006 (see also [Kashino} {et~al.}, 2017]),
\begin{equation}
\log{M_0} / (h^{-1}M_\odot) = 0.76 \log{M_1}/(h^{-1}M_\odot) + 2.3.
\end{equation}
Thus the number of free parameters in the Zheng HOD model together
with the fake line fraction is five, $\Theta = (M_{\rm min}, M_1,
\sigma_{\log{M}}, \alpha, f_{\rm fake})$.  As is the case of the
analysis with the Geach HOD model, we apply the Gaussian priors on
$\alpha$ and $f_{\rm fake}$ as $\alpha= 1.00\pm 0.20$ and $f_{\rm
fake} = 0.140\pm 0.060$ , and flat priors on the other 3 parameters.
We use the data with the same angular separation range as the analysis
of the Geach HOD model (section \ref{sec:constraint}),
$0.0015 <\theta < 1$ [deg], the degree of freedom is
$\nu = N_{\rm bin} +1 - 5= 10$.

% Table 4
\begin{table*}[bt!]
\caption{Priors and constraints of the HOD parameters for Zheng model${}^*$}
\begin{center}
\begin{tabular}{ll | ll | ll}
\hline %\hline
& & \multicolumn{2}{c|}{{\it NB816}}  &  \multicolumn{2}{c}{{\it NB921}}   \\
Parameter  & \multicolumn{1}{l |}{Prior} & Best-fit & Posterior PDF & Best-fit & Posterior PDF \\
\hline %\hline
$\log{M_{\rm min}} / (h^{-1}M_\odot)$ & \multicolumn{1}{l |}{None}          & $12.16$ & $12.14^{+0.19}_{-0.18}$ &         $11.88$ & $11.95^{+0.15}_{-0.10}$ \\
$\log{M_1} / (h^{-1}M_\odot)$             & \multicolumn{1}{l |}{None}       & $13.070$ & $13.062^{+0.084}_{-0.090}$ & $13.000$ & $13.017^{+0.075}_{-0.078}$ \\
$\sigma_{\log{M}}$                              & \multicolumn{1}{l |}{$[0,1]$}            & $0.68$ & $0.65^{+0.21}_{-0.32}$ &             $0.18$ & $0.34^{+ 0.25}_{- 0.23}$  \\
$\alpha$                                               & \multicolumn{1}{l |}{$1.00\pm 0.20$}     & $0.99$ & $0.98^{+0.11}_{-0.12}$ &             $1.11$ & $1.07^{+ 0.10}_{- 0.11}$ \\
$f_{\rm fake}$                                      & \multicolumn{1}{l |}{$0.140\pm 0.060$} & $0.123$ & $0.133^{+0.048}_{-0.048}$  &    $0.115$ & $0.098^{+ 0.041}_{- 0.043}$  \\
$\chi^2/\nu$ \ \ ($\nu=10$)& & $1.86$& &$1.42$ & \\
\hline %\hline
Inferred quantity & \multicolumn{1}{l |}{Measurement}  & Best-fit & Posterior PDF & Best-fit & Posterior PDF \\
\hline
%\multicolumn{2}{l |}{
$\log{n_g / (h^{-1}{\rm Mpc})^{-3}}$ & $-2.259\pm 0.068$ ({\it NB816})   & $-2.327$   & $-2.325^{+0.075}_{-0.073}$ &    $ -2.369$ & $-2.379^{+ 0.075}_{- 0.074}$ \\
& $-2.344\pm 0.070$ ({\it NB921}) & &  &      &  \\
$f_{\rm sat}$ & \multicolumn{1}{c |}{$\cdots$}                                            & $0.100$ & $0.102^{+0.013}_{-0.012}$ &     $0.098$ & $0.094^{ + 0.012}_{ - 0.011}$ \\
$b_{\rm eff}$ & \multicolumn{1}{c |}{$\cdots$}                                            & $1.659$   & $1.668^{+0.102}_{-0.089}$ &           $2.034$ & $1.983^{+0.077}_{-0.095}$ \\
$\log{M_{\rm eff}} / (h^{-1}M_\odot) $ & \multicolumn{1}{c |}{$\cdots$}      & $12.645$ & $12.647^{+0.052}_{-0.050}$ & $12.591$ & $12.567^{+0.043}_{-0.046}$  \\
\hline
\end{tabular}
\end{center}
\label{tab:hod_zheng}
\begin{tabnote}
$*$ In the ``Prior'' column the ranges specified in brackets are for
uniform priors while the others we quote the mean and standard
deviation of the Gaussian priors. The column of ``Best-fit'' shows
the parameter set which gives the minimum value of $\chi^2$. In the
column of ``Posterior PDF'', the central value is a median and the
error means $16-84$ percentiles after other parameters are
marginalized over.
The measured number density $\log{n_g}$ includes non-[OII] emitters,
and thus its best-fitting values differ from the measure ones by $\log{(1-f_{\rm fake})} \sim -0.06$.
\end{tabnote}
\end{table*}

The resulting constraints on the HOD parameters of Zheng's model are
shown in figure \ref{fig:result_hod_zheng} and summarized in table
\ref{tab:hod_zheng}.  Overall, the constraints on the Zheng HOD model are
tighter than those on the Geach model due to the fewer number of
parameters.  The posterior distribution of the HOD $N(M)$ is shown as
the red dashed curve in figure \ref{fig:result_hodfig_contreras}.
Interestingly, the two models predict very similar HODs, particularly
the central ones, for the observed angular correlation function of
[OII] emitters.  The discrepancy seen at $M>10^{13}M_\odot/h$ is
reasonable because such massive halos are rare as hosts of [OII]
emitters and since the varied range of $F_s$ in the Geach model is
$0\leq F_s \leq 1$, the larger number of $\left\langle N_{\rm
sat}\right\rangle$ is allowed in the model (see the factor of $1/2$
in the Zheng model).
Furthermore, since the Zheng model has
$\left\langle N_{\rm cen}\right\rangle =1$ at the high mass end,
it would make sense for the Zheng model to have the lower number of
$\left\langle N_{\rm sat}\right\rangle$ when the total number density is fixed.
Nevertheless, the two $N(M)$'s are consistent with each other within $2-\sigma$.
The best-fitting correlation function is shown as the red dashed curve
in the upper and lower panels of figure \ref{fig:result_xi_contreras}
for $z=1.19$ and $z=1.47$, respectively.
As expected, the overall shape of
the correlation function is very similar with that of the Geach model
depicted by the red solid curve.
However, one can see that the Geach model shows a better fit to the
measured $w(\theta)$ at the smallest scales, which would imply that
the ELGs in our sample prefer the HOD model which allows
massive halos to have no central ELG.

The posterior distribution of the derived physical parameters is shown
by the function colored orange in figure
\ref{fig:result_params_b_m_n_f}.  Since the number density is
primarily constrained by our prior, the two models predict the same
values.  We find that the effective bias $b_{\rm eff}$ and halo mass
$\log{M_{\rm eff}}$ determined from the Zheng HOD model are also
consistent with those from the Geach model within the statistical
uncertainties.  On the other hand, the satellite fraction is derived
to be $f_{\rm sat} =0.102^{+0.013}_{-0.012}$ and
$0.093^{+0.012}_{-0.011}$ at $z=1.19$ and $1.47$, respectively.  These
values are slightly lower than those determined based on the Geach HOD
model, respectively $0.158^{+0.114}_{-0.047}$ and
$0.159^{+0.109}_{-0.049}$, though the differences are not significant.
This again reflects the fact that the Geach model adopts a more
flexible parameterization and the larger amplitude of $\left\langle
N_{\rm sat}\right\rangle$ is allowed.  Hence the 1-d posterior of
$f_{\rm sat}$ obtained from the Geach model has a long tail toward the
high $f_{\rm sat}$ end.

Overall, the differences of halo parameters for [OII] emitters
determined from the Geach and Zheng HOD models are not significant
compared to the current statistical uncertainties.  However, such
small differences need to be carefully treated in the future
clustering analysis ELGs from large galaxy surveys such as the Subaru
PFS, DESI, and Euclid surveys.

% Section 6
\section{Summary}\label{sec:conclusion}

In this paper we have studied the clustering of emission line galaxies
at $z>1$ and physical properties of dark matter halos which host them.
For this purpose we used [OII] emitters detected at $z=1.19$ and
$1.47$ using two narrow band filters, {\it NB816} and {\it NB921},
respectively, in the Subaru HSC survey [Hayashi} {et~al.}, 2020].  We then
measured the angular correlation functions of 8302 ($z=1.19$) and 9578
($z=1.47$) of [OII] emitters.  Using a simple model of the
nonlinear correlation function with the linear galaxy bias factor, we
measured the bias of [OII] emitters as $b=1.61^{+0.13}_{-0.11}$ and
$2.09^{+0.17}_{-0.15}$ at $z=1.19$ and $1.47$, respectively.
We also found that the bias monotonically increases with the line luminosity.

We, for the first time, performed an HOD analysis for the measured
correlation functions of [OII] emitters based on a model developed to
describe the population of galaxies selected by a star forming rate
[Geach} {et~al.}, 2012].  We varied eight parameters simultaneously with
only a few priors.  Nevertheless, the two HOD parameters related to
the host halo mass have been well constrained given
the large sample from the HSC survey.

Based on the constrained HOD parameters, we have derived parameters
which describe properties of the host halos of [OII] emitters, such as
the effective bias, effective halo mass and satellite fraction.  The
effective biases are determined as $b_{\rm
eff}=1.701^{+0.083}_{-0.110}$ and $1.981^{+0.072}_{-0.068}$ at
$z=1.19$ and $1.47$, respectively.  They are consistent with the
linear bias values determined based on a simple nonlinear dark matter
model.  The effective halo masses and satellite fractions are derived
as $\log{M_{\rm eff}/(h^{-1}M_\odot)} \simeq 12.7~ (12.6)$ and $f_{\rm
sat}\simeq 0.16~(0.16)$ at $z=1.19~(1.47)$.  They are consistent
with the result of Guo} {et~al.}, 2019 who studied the relation among
$f_{\rm sat}$, $M_{\rm eff}$ and the stellar mass using the ELG sample
from the eBOSS survey.
Furthermore, the determined effective halo masses are in good
agreement with previous studies of similar tracers with the
mass-assembly history with $M(z=0)=1.5\times 10^{13}h^{-1}M_\odot$.

As stated in section \ref{sec:introduction}, ELGs will be a main
tracer in large galaxy surveys at high redshifts.  Particularly, [OII]
emitters will be targeted by many ongoing/future surveys such as the
Subaru PFS survey [Takada} {et~al.}, 2014].  The result presented in this
paper is useful to construct a mock catalog for cosmological analysis.
The physical properties of halos hosting [OII] emitters revealed in
this paper can also be applicable to the Subaru PFS, DESI, and other
forthcoming surveys.

% Figure 13
\begin{figure*}%[bt]
\begin{center}
\vspace{-.8cm}
\FigureFile(172mm,172mm){w_nlmatter_individual.pdf}
\end{center}
\caption{Angular correlation functions of [OII] emitters of each of
the four fields measured at $z=1.19$ (top) and $z=1.47$ (bottom).
In each column, the measurement from each of the individual four
fields is presented.  The black points are the raw measurement
without the correction of the integral constraints, $\hat{w}_k$
($1\leq k \leq 4$).  The red points are after the integral
constraint (denoted by IC in the legend) is taken into account with
the linearly-biased dark matter model, $\hat{w}_k+(1-f_{\rm fake})^2w_{\Omega,k}$,
with the best-fitting parameters obtained
for the combined field, $(b,f_{\rm fake})=(1.60, 0.140)$ $(z=1.19)$
and $(b,f_{\rm fake})=(2.08, 0.140)$ $(z=1.47)$.  For the reference,
the angular correlation function measured for the combined field,
$\hat{w}+(1-f_{\rm fake})^2w_\Omega$, is shown as the blue dashed
curve.  }
\label{fig:w_individual}
\end{figure*}

In this analysis, both the shape of the angular correlation function
and number density are constrained by the HOD modeling.  However,
there is a claim that in 3-dimensional analysis the HOD would fail to
explain the redshift-space clustering and number density
simultaneously [Reid} {et~al.}, 2014, Saito} {et~al.}, 2016].  To investigate this, we
need a spectroscopic survey, and we can do it with the PFS, DESI or
{\it Euclid} surveys.

The NB data of the HSC-SSP PDR2 provide emission lines of not only
[OII] but also H$\alpha$ and [OIII], ranging from $z\sim 0.4$ to
$z\sim 1.6$ [Hayashi} {et~al.}, 2020].  Redshift evolution of physical
properties such as the halo mass for ELGs using these multiple
emission lines will be studied in our future work.

% Table 5

\begin{table*}[bt!]
\caption{Priors and constraints of the HOD parameters for Geach and Zheng models obtained using the clustering information only${}^*$
\begin{center}
\begin{tabular}{ll | l | l}
\multicolumn{4}{c}{Geach HOD model} \\
\hline %\hline
&          & \multicolumn{1}{c |}{{\it NB816}} & \multicolumn{1}{c}{{\it NB921}} \\
Parameter  & Prior & Posterior PDF & Posterior PDF \\
\hline %\hline
$\log{M_{\rm c}} / (h^{-1}M_\odot)$     &  None               & $11.50^{+0.72}_{-0.28}$     & $11.70^{+ 0.25} _{- 0.20}$ \\
$\log{M_{\rm min}} / (h^{-1}M_\odot)$ & None                & $11.40^{+1.25}_{-1.00}$    & $11.89^{+0.76}_{-1.36}$ \\
$\sigma_{\log{M}}$                              & $[0,1]$             & $0.27^{+0.33}_{-0.17}$      & $0.22^{+ 0.16}_{- 0.13}$  \\
$\alpha$                                               & $1.00\pm 0.20$  & $0.96^{+0.11}_{-0.09}$       & $1.10^{+ 0.13}_{- 0.13}$ \\
$F_c^A$                                               & $[0, 0.5]$         & $0.22^{+0.14}_{-0.07}$        & $0.17^{+ 0.18}_{- 0.06}$  \\
$F_c^B$                                               & $[0,1]$          & $0.56^{+0.32}_{-0.39}$        & $0.61^{+ 0.27}_{- 0.34}$  \\
$F_s$                                                  & $[0, 1]$            & $0.54^{+0.32}_{-0.30}$      & $0.52^{+ 0.32}_{- 0.34}$  \\
$f_{\rm fake}$                                      & $0.140\pm 0.060$  & $0.124^{+0.054}_{-0.051}$ & $0.091^{+ 0.054}_{- 0.046}$  \\
\hline %\hline
Inferred quantity & Measurement & Posterior PDF  &  Posterior PDF \\
\hline
$\log{n_g / (h^{-1}{\rm Mpc})^{-3}}$ & \multicolumn{1}{c |}{$\cdots$} &  $-1.47^{+0.89}_{-0.80}$ & $-1.96^{+1.18}_{-0.42}$  \\
$f_{\rm sat}$ & \multicolumn{1}{c |}{$\cdots$}                                        & $0.58^{+0.33}_{-0.43}$  & $0.38^{ + 0.54}_{ - 0.23}$ \\
$b_{\rm eff}$ & \multicolumn{1}{c |}{$\cdots$}                                       & $1.64^{+0.12}_{-0.10}$  & $1.94^{+0.13}_{-0.11}$ \\
$\log{M_{\rm eff}} / (h^{-1}M_\odot) $ & \multicolumn{1}{c |}{$\cdots$}  & $12.70^{+0.10}_{-0.09}$  & $12.63^{+0.12}_{-0.10}$  \\
\hline \\
\multicolumn{4}{c}{Zheng HOD model} \\
\hline %\hline
& & \multicolumn{1}{c|}{{\it NB816}}  &  \multicolumn{1}{c}{{\it NB921}}   \\
Parameter  & \multicolumn{1}{l |}{Prior} & Posterior PDF & Posterior PDF \\
\hline %\hline
$\log{M_{\rm min}} / (h^{-1}M_\odot)$ & \multicolumn{1}{l |}{None}         & $12.01^{+0.62}_{-0.61}$ &         $11.85^{+ 0.56}_{-0.26}$ \\
$\log{M_1} / (h^{-1}M_\odot)$             & \multicolumn{1}{l |}{None}       & $12.97^{+0.43}_{-0.64}$  &    $12.91^{+0.36}_{-0.31}$\\
$\sigma_{\log{M}}$                              & \multicolumn{1}{l |}{$[0,1]$}            & $0.58^{+0.30}_{-0.41}$ &            $0.32^{+0.41}_{-0.22}$  \\
$\alpha$                                               & \multicolumn{1}{l |}{$1.00\pm 0.20$}     & $0.94^{+0.11}_{-0.11}$ &         $1.02^{+0.11}_{-0.11}$ \\
$f_{\rm fake}$                                      & \multicolumn{1}{l |}{$0.140\pm 0.060$} & $0.123^{+0.060}_{-0.063}$  &  $0.077^{+ 0.055}_{- 0.047}$  \\
\hline %\hline
Inferred quantity & \multicolumn{1}{l |}{Measurement}  & Posterior PDF & Posterior PDF \\
\hline
$\log{n_g / (h^{-1}{\rm Mpc})^{-3}}$ & \multicolumn{1}{c |}{$\cdots$}    & $-2.24^{+0.56}_{- 0.44}$ &           $-2.27^{+0.28}_{-0.38}$ \\
$f_{\rm sat}$ & \multicolumn{1}{c |}{$\cdots$}                                       & $0.113^{+0.100}_{-0.044}$ &     $0.106^{+ 0.047}_{-0.040}$ \\
$b_{\rm eff}$ & \multicolumn{1}{c |}{$\cdots$}                                       & $1.64^{+ 0.15}_{-0.15}$ &           $1.91^{+0.13}_{-0.11}$ \\
$\log{M_{\rm eff}} / (h^{-1}M_\odot) $ & \multicolumn{1}{c |}{$\cdots$} & $12.633^{+0.102}_{-0.107}$ & $12.534^{+0.091}_{-0.077}$  \\
\hline
\end{tabular}
\end{center}
\label{tab:hod_nong}
\begin{tabnote}
$*$ In the ``Prior'' column the ranges specified in brackets are for
uniform priors while the others we quote the mean and standard
deviation of the Gaussian priors. In the
column of ``Posterior PDF'', the central value is a median and the
error means $16-84$ percentiles after other parameters are
marginalized over.
\end{tabnote}
\end{table*}

\section*{Acknowledgments }
TO thanks the members of the Subaru PFS Cosmology Working Group,
particularly Masahiro Takada, Eiichiro Komatsu, and Ryu Makiya for useful
correspondences during the regular telecon.
TO is grateful to Shogo Ishikawa for discussion.
We also thank the anonymous referee for the careful reading and suggestions.
TO acknowledges support
from the Ministry of Science and Technology of Taiwan under Grants
No. MOST 109-2112-M-001-027- and the Career Development Award,
Academia Sinica (AS-CDA-108-M02) for the period of 2019-2023.  KO is
supported by JSPS Overseas Research Fellowships.  This work is based
on data collected at Subaru Telescope, which is operated by the
National Astronomical Observatory of Japan.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical
communities of Japan and Taiwan, and Princeton University. The HSC
instrumentation and software were developed by the National
Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the
Physics and Mathematics of the Universe (Kavli IPMU), the University
of Tokyo, the High Energy Accelerator Research Organization (KEK), the
Academia Sinica Institute for Astronomy and Astrophysics in Taiwan
(ASIAA), and Princeton University. Funding was contributed by the
FIRST program from the Japanese Cabinet Office, the Ministry of
Education, Culture, Sports, Science and Technology (MEXT), the Japan
Society for the Promotion of Science (JSPS), Japan Science and
Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli
IPMU, KEK, ASIAA, and Princeton University.

\appendix
\section{Correlation functions in individual fields}\label{sec:w_individual}

In order to confirm the consistency of the angular correlation
functions among different fields, we present the measurement from each
of the individual four fields, $\hat{w}_k$, as the black points in
Figure \ref{fig:w_individual}.
The square roots of the diagonal components of the covariance matrix,
$C_{k,ii}^{1/2}$, are shown as the error bars.
The results for $z=1.19$ and $z=1.47$ are shown in the upper and lower
panels, respectively.  In order to compare the correlation functions
between different fields, we need to take into account the correction
of the integral constraint because the survey volume of each field is
To estimate the integral constraint for $k$-th field, $w_{\Omega,k}$,
the random-random count for the field, $RR_k$, is used as $RR$ in
equation (\ref{eq:ic_rr}), and
\be
w_{\Omega,k}(\Theta) = \frac{\sum_i{w(\theta_i;\Theta)RR_k(\theta_i)}} {\sum_i{RR_k(\theta_i)}}. \label{eq:ic_rr_k}
\ee
For $w(\theta_i;\Theta)$ in equation (\ref{eq:ic_rr_k}), we use the
best-fitting linearly-biased dark matter model for the combined field
obtained in section \ref{sec:nl_matter}, $(b,f_{\rm
fake})=(1.60,0.140)$ for $z=1.19$ and $(b,f_{\rm
fake})=(2.08,0.140)$ for $z=1.47$.
The red points show the result with the integral constraint
correction, $\hat{w}_k+(1-f_{\rm fake})^2 w_\Omega$.  For comparison,
the measurement from all four fields with the integral
constraint calculated using the same model is plotted as the blue
dashed curve. One can see that the measurement from each field is
largely consistent with each other.  Thus, throughout this paper we
perform the statistical analysis only for the data combined over all
the four fields.

\section{Constraints on HOD parameters using clustering only}\label{sec:hod_wo_nobs}
Given an HOD model for a certain galaxy population,
the average number density is calculated by equation (\ref{eq:ng_hod}).
One therefore needs to simultaneously analyze the observed correlation function and
number density to constrain the HOD model, as performed in section \ref{sec:constraint}.
In this appendix, we present the constraints on the HOD parameters based on Geach and Zheng models
without using the information of the measured number density of [OII] emitters for readers who are interested in the constraining power
coming from the clustering only (see, e.g., [Matsuoka} {et~al.}, 2011, Martinez-Manso} {et~al.}, 2015] for a similar attempt).

Table \ref{tab:hod_nong} shows these results without using the abundance information.
The halo mass parameters are strongly constrained by the abundance of [OII] emitters.
The constraints on these parameters are thus not as tight as those presented in tables \ref{tab:hod_geach} and \ref{tab:hod_zheng}, as expected:
the errors on $M_{\rm c}$ and $M_{\rm min }$ constraints without the abundance information in the Geach model are $\sim 1.5$ and $\sim 3$ times
larger than those with the simultaneous analysis of clustering and abundance, and the errors on $M_{\rm min}$ and $M_1$ constraints without the abundance information in
the Zheng model are $\sim 3$ and $\sim 4$ times larger.
There are not significant changes in the constraints of the other HOD parameters.
Interestingly, even though we do not use the number density to constrain the HOD parameters,
the posterior of the number density is consistent with the observed one.

\begin{thebibliography}{}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi

[Abbas} {et~al.}, 2010]
{Abbas}, U., {de la Torre}, S., {Le F{\`e}vre}, O., {et~al.} 2010, \mnras, 406,

[Adams} {et~al.}, 2011]
{Adams}, J.~J., {Blanc}, G.~A., {Hill}, G.~J., {et~al.} 2011, \apjs, 192, 5

[Addison} {et~al.}, 2019]
{Addison}, G.~E., {Bennett}, C.~L., {Jeong}, D., {Komatsu}, E., \& {Weiland},
J.~L. 2019, \apj, 879, 15

[Aihara} {et~al.}, 2018]
{Aihara}, H., {Armstrong}, R., {Bickerton}, S., {et~al.} 2018{\natexlab{a}},
\pasj, 70, S8

[Aihara} {et~al.}, 2018]
{Aihara}, H., {Arimoto}, N., {Armstrong}, R., {et~al.} 2018{\natexlab{b}},
\pasj, 70, S4

[Aihara} {et~al.}, 2019]
{Aihara}, H., {AlSayyad}, Y., {Ando}, M., {et~al.} 2019, \pasj, 71, 114

[Avila} {et~al.}, 2020]
{Avila}, S., {Gonzalez-Perez}, V., {Mohammad}, F.~G., {et~al.} 2020, \mnras,
499, 5486

[Behroozi} {et~al.}, 2013]
{Behroozi}, P.~S., {Wechsler}, R.~H., \& {Conroy}, C. 2013, \apj, 770, 57

[Berlind} \& {Weinberg}, 2002]
{Berlind}, A.~A., \& {Weinberg}, D.~H. 2002, \apj, 575, 587

[Blake} {et~al.}, 2008]
{Blake}, C., {Collister}, A., \& {Lahav}, O. 2008, \mnras, 385, 1257

[Bocquet \& Carter, 2016]
Bocquet, S., \& Carter, F.~W. 2016, The Journal of Open Source Software, 1,

[Bosch} {et~al.}, 2018]
{Bosch}, J., {Armstrong}, R., {Bickerton}, S., {et~al.} 2018, \pasj, 70, S5

[Bullock} {et~al.}, 2002]
{Bullock}, J.~S., {Wechsler}, R.~H., \& {Somerville}, R.~S. 2002, \mnras, 329,

[Cochrane} {et~al.}, 2018]
{Cochrane}, R.~K., {Best}, P.~N., {Sobral}, D., {et~al.} 2018, \mnras, 475,

[Cochrane} {et~al.}, 2017]
---. 2017, \mnras, 469, 2913

[Comparat} {et~al.}, 2016]
{Comparat}, J., {Zhu}, G., {Gonzalez-Perez}, V., {et~al.} 2016, \mnras, 461,

[Conroy} {et~al.}, 2006]
{Conroy}, C., {Wechsler}, R.~H., \& {Kravtsov}, A.~V. 2006, \apj, 647, 201

[Contreras} {et~al.}, 2013]
{Contreras}, S., {Baugh}, C.~M., {Norberg}, P., \& {Padilla}, N. 2013, \mnras,
432, 2717

[Cooray} \& {Sheth}, 2002]
{Cooray}, A., \& {Sheth}, R. 2002, \physrep, 372, 1

[Coupon} {et~al.}, 2018]
{Coupon}, J., {Czakon}, N., {Bosch}, J., {et~al.} 2018, \pasj, 70, S7

[Coupon} {et~al.}, 2012]
{Coupon}, J., {Kilbinger}, M., {McCracken}, H.~J., {et~al.} 2012, \aap, 542, A5

[Crocce} {et~al.}, 2011]
{Crocce}, M., {Cabr{\'e}}, A., \& {Gazta{\~n}aga}, E. 2011, \mnras, 414, 329

[Dawson} {et~al.}, 2016]
{Dawson}, K.~S., {Kneib}, J.-P., {Percival}, W.~J., {et~al.} 2016, \aj, 151, 44

[de la Torre} {et~al.}, 2013]
{de la Torre}, S., {Guzzo}, L., {Peacock}, J.~A., {et~al.} 2013, \aap, 557, A54

[DESI Collaboration} {et~al.}, 2016]
{DESI Collaboration}, {Aghamousa}, A., {Aguilar}, J., {et~al.} 2016, arXiv
e-prints, arXiv:1611.00036

[Dor{\'e}} {et~al.}, 2014]
{Dor{\'e}}, O., {Bock}, J., {Ashby}, M., {et~al.} 2014, arXiv e-prints,

[Duffy} {et~al.}, 2008]
{Duffy}, A.~R., {Schaye}, J., {Kay}, S.~T., \& {Dalla Vecchia}, C. 2008,
\mnras, 390, L64

[Durkalec} {et~al.}, 2015]
{Durkalec}, A., {Le F{\`e}vre}, O., {Pollo}, A., {et~al.} 2015, \aap, 583, A128

[Efstathiou} {et~al.}, 1991]
{Efstathiou}, G., {Bernstein}, G., {Katz}, N., {Tyson}, J.~A., \&
{Guhathakurta}, P. 1991, \apjl, 380, L47

[Favole} {et~al.}, 2016]
{Favole}, G., {Comparat}, J., {Prada}, F., {et~al.} 2016, \mnras, 461, 3421

[Favole} {et~al.}, 2020]
{Favole}, G., {Gonzalez-Perez}, V., {Stoppacher}, D., {et~al.} 2020, \mnras,
497, 5432

[Foreman-Mackey} {et~al.}, 2013]
{Foreman-Mackey}, D., {Hogg}, D.~W., {Lang}, D., \& {Goodman}, J. 2013, \pasp,
125, 306

[Furusawa} {et~al.}, 2018]
{Furusawa}, H., {Koike}, M., {Takata}, T., {et~al.} 2018, \pasj, 70, S3

[Gao} \& {Jing}, 2020]
{Gao}, H., \& {Jing}, Y.~P. 2020, arXiv e-prints, arXiv:2007.12878

[Geach} {et~al.}, 2008]
{Geach}, J.~E., {Smail}, I., {Best}, P.~N., {et~al.} 2008, \mnras, 388, 1473

[Geach} {et~al.}, 2012]
{Geach}, J.~E., {Sobral}, D., {Hickox}, R.~C., {et~al.} 2012, \mnras, 426, 679

[Gonzalez-Perez} {et~al.}, 2018]
{Gonzalez-Perez}, V., {Comparat}, J., {Norberg}, P., {et~al.} 2018, \mnras,
474, 4024

[Gonzalez-Perez} {et~al.}, 2020]
{Gonzalez-Perez}, V., {Cui}, W., {Contreras}, S., {et~al.} 2020, \mnras, 498,

[Guo} {et~al.}, 2014]
{Guo}, H., {Zheng}, Z., {Zehavi}, I., {et~al.} 2014, \mnras, 441, 2398

[Guo} {et~al.}, 2019]
{Guo}, H., {Yang}, X., {Raichoor}, A., {et~al.} 2019, \apj, 871, 147

[Hadzhiyska} {et~al.}, 2020]
{Hadzhiyska}, B., {Tacchella}, S., {Bose}, S., \& {Eisenstein}, D.~J. 2020,
arXiv e-prints, arXiv:2011.05331

[Hamana} {et~al.}, 2004]
{Hamana}, T., {Ouchi}, M., {Shimasaku}, K., {Kayo}, I., \& {Suto}, Y. 2004,
\mnras, 347, 813

[Harikane} {et~al.}, 2016]
{Harikane}, Y., {Ouchi}, M., {Ono}, Y., {et~al.} 2016, \apj, 821, 123

[Harikane} {et~al.}, 2018]
---. 2018, \pasj, 70, S11

[Hartlap} {et~al.}, 2007]
{Hartlap}, J., {Simon}, P., \& {Schneider}, P. 2007, \aap, 464, 399

[Hayashi} {et~al.}, 2013]
{Hayashi}, M., {Sobral}, D., {Best}, P.~N., {Smail}, I., \& {Kodama}, T. 2013,
\mnras, 430, 1042

[Hayashi} {et~al.}, 2015]
{Hayashi}, M., {Ly}, C., {Shimasaku}, K., {et~al.} 2015, \pasj, 67, 80

[Hayashi} {et~al.}, 2018]
{Hayashi}, M., {Tanaka}, M., {Shimakawa}, R., {et~al.} 2018, \pasj, 70, S17

[Hayashi} {et~al.}, 2020]
{Hayashi}, M., {Shimakawa}, R., {Tanaka}, M., {et~al.} 2020, \pasj, 72, 86

[Hikage} {et~al.}, 2013]
{Hikage}, C., {Mandelbaum}, R., {Takada}, M., \& {Spergel}, D.~N. 2013, \mnras,
435, 2345

[Hong} {et~al.}, 2019]
{Hong}, S., {Dey}, A., {Lee}, K.-S., {et~al.} 2019, \mnras, 483, 3950

[Hsieh} \& {Yee}, 2014]
{Hsieh}, B.~C., \& {Yee}, H.~K.~C. 2014, \apj, 792, 102

[Huang} {et~al.}, 2018]
{Huang}, S., {Leauthaud}, A., {Murata}, R., {et~al.} 2018, \pasj, 70, S6

[Ishikawa} {et~al.}, 2017]
{Ishikawa}, S., {Kashikawa}, N., {Toshikawa}, J., {et~al.} 2017, \apj, 841, 8

[Ishikawa} {et~al.}, 2020]
{Ishikawa}, S., {Kashikawa}, N., {Tanaka}, M., {et~al.} 2020, \apj, 904, 128

[Jing} {et~al.}, 1998]
{Jing}, Y.~P., {Mo}, H.~J., \& {Boerner}, G. 1998, \apj, 494, 1

[Kaiser}, 1984]
{Kaiser}, N. 1984, \apjl, 284, L9

[Kashino} {et~al.}, 2017]
{Kashino}, D., {More}, S., {Silverman}, J.~D., {et~al.} 2017, \apj, 843, 138

[Kayo} \& {Oguri}, 2012]
{Kayo}, I., \& {Oguri}, M. 2012, \mnras, 424, 1363

[Khostovan} {et~al.}, 2018]
{Khostovan}, A.~A., {Sobral}, D., {Mobasher}, B., {et~al.} 2018, \mnras, 478,

[Kim} {et~al.}, 2011]
{Kim}, H.-S., {Baugh}, C.~M., {Benson}, A.~J., {et~al.} 2011, \mnras, 414, 2367

[Koda} {et~al.}, 2016]
{Koda}, J., {Blake}, C., {Beutler}, F., {Kazin}, E., \& {Marin}, F. 2016,
\mnras, 459, 2118

[Komiyama} {et~al.}, 2018]
{Komiyama}, Y., {Obuchi}, Y., {Nakaya}, H., {et~al.} 2018, \pasj, 70, S2

[Landy} \& {Szalay}, 1993]
{Landy}, S.~D., \& {Szalay}, A.~S. 1993, \apj, 412, 64

[Laureijs} {et~al.}, 2011]
{Laureijs}, R., {Amiaux}, J., {Arduini}, S., {et~al.} 2011, arXiv e-prints,

[Leung} {et~al.}, 2017]
{Leung}, A.~S., {Acquaviva}, V., {Gawiser}, E., {et~al.} 2017, \apj, 843, 130

[Limber}, 1954]
{Limber}, D.~N. 1954, \apj, 119, 655

[Lupton}, 1993]
{Lupton}, R. 1993, {Statistics in theory and practice}

[Ly} {et~al.}, 2012]
{Ly}, C., {Malkan}, M.~A., {Kashikawa}, N., {et~al.} 2012, \apj, 757, 63

[Ma} \& {Fry}, 2000]
{Ma}, C.-P., \& {Fry}, J.~N. 2000, \apj, 543, 503

[Martinez-Manso} {et~al.}, 2015]
{Martinez-Manso}, J., {Gonzalez}, A.~H., {Ashby}, M. L.~N., {et~al.} 2015,
\mnras, 446, 169

[Matsuoka} {et~al.}, 2011]
{Matsuoka}, Y., {Masaki}, S., {Kawara}, K., \& {Sugiyama}, N. 2011, \mnras,
410, 548

[Miyazaki} {et~al.}, 2012]
{Miyazaki}, S., {Komiyama}, Y., {Nakaya}, H., {et~al.} 2012, in \procspie, Vol.
8446, Ground-based and Airborne Instrumentation for Astronomy IV, 84460Z

[Miyazaki} {et~al.}, 2018]
{Miyazaki}, S., {Komiyama}, Y., {Kawanomoto}, S., {et~al.} 2018, \pasj, 70, S1

[Moustakas} \& {Somerville}, 2002]
{Moustakas}, L.~A., \& {Somerville}, R.~S. 2002, \apj, 577, 1

[Murray} {et~al.}, 2013]
{Murray}, S.~G., {Power}, C., \& {Robotham}, A.~S.~G. 2013, Astronomy and
Computing, 3, 23

[Navarro} {et~al.}, 1997]
{Navarro}, J.~F., {Frenk}, C.~S., \& {White}, S.~D.~M. 1997, \apj, 490, 493

[Oguri} \& {Lin}, 2015]
{Oguri}, M., \& {Lin}, Y.-T. 2015, \apj, 801, 94

[Okada} {et~al.}, 2016]
{Okada}, H., {Totani}, T., {Tonegawa}, M., {et~al.} 2016, \pasj, 68, 47

[Okumura} {et~al.}, 2015]
{Okumura}, T., {Hand}, N., {Seljak}, U., {Vlah}, Z., \& {Desjacques}, V. 2015,
\prd, 92, 103516

[Okumura} {et~al.}, 2017]
{Okumura}, T., {Takada}, M., {More}, S., \& {Masaki}, S. 2017, \mnras, 469, 459

[Okumura} {et~al.}, 2016]
{Okumura}, T., {Hikage}, C., {Totani}, T., {et~al.} 2016, \pasj, 68, 38

[Osato} \& {Okumura}, 2021]
{Osato}, K., \& {Okumura}, T. 2021, in prep.

[Peacock} \& {Smith}, 2000]
{Peacock}, J.~A., \& {Smith}, R.~E. 2000, \mnras, 318, 1144

[Peebles}, 1973]
{Peebles}, P.~J.~E. 1973, \apj, 185, 413

[Peebles}, 1980]
---. 1980, {The large-scale structure of the universe} (Princeton, N.J.,
Princeton Univ. Press)

[Peebles} \& {Groth}, 1976]
{Peebles}, P.~J.~E., \& {Groth}, E.~J. 1976, \aap, 53, 131

[Planck Collaboration} {et~al.}, 2016]
{Planck Collaboration}, {Ade}, P.~A.~R., {Aghanim}, N., {et~al.} 2016, \aap,
594, A13

[Reid} {et~al.}, 2014]
{Reid}, B.~A., {Seo}, H.-J., {Leauthaud}, A., {Tinker}, J.~L., \& {White}, M.
2014, \mnras, 444, 476

[Richardson} {et~al.}, 2012]
{Richardson}, J., {Zheng}, Z., {Chatterjee}, S., {Nagai}, D., \& {Shen}, Y.
2012, \apj, 755, 30

[Saito} {et~al.}, 2020]
{Saito}, S., {de la Torre}, S., {Ilbert}, O., {et~al.} 2020, \mnras, 494, 199

[Saito} {et~al.}, 2016]
{Saito}, S., {Leauthaud}, A., {Hearin}, A.~P., {et~al.} 2016, \mnras, 460, 1457

[Scoccimarro} {et~al.}, 2001]
{Scoccimarro}, R., {Sheth}, R.~K., {Hui}, L., \& {Jain}, B. 2001, \apj, 546, 20

[Seljak}, 2000]
{Seljak}, U. 2000, \mnras, 318, 203

[Simon}, 2007]
{Simon}, P. 2007, \aap, 473, 711

[Takada} {et~al.}, 2014]
{Takada}, M., {Ellis}, R.~S., {Chiba}, M., {et~al.} 2014, \pasj, 66, 1

[Takahashi} {et~al.}, 2012]
{Takahashi}, R., {Sato}, M., {Nishimichi}, T., {Taruya}, A., \& {Oguri}, M.
2012, \apj, 761, 152

[Tinker} {et~al.}, 2010]
{Tinker}, J.~L., {Robertson}, B.~E., {Kravtsov}, A.~V., {et~al.} 2010, \apj,
724, 878

[Tonegawa} {et~al.}, 2015]
{Tonegawa}, M., {Totani}, T., {Okada}, H., {et~al.} 2015, \pasj, 67, 81

[Totsuji} \& {Kihara}, 1969]
{Totsuji}, H., \& {Kihara}, T. 1969, \pasj, 21, 221

[van den Bosch} {et~al.}, 2013]
{van den Bosch}, F.~C., {More}, S., {Cacciato}, M., {Mo}, H., \& {Yang}, X.
2013, \mnras, 430, 725

[Wake} {et~al.}, 2011]
{Wake}, D.~A., {Whitaker}, K.~E., {Labb{\'e}}, I., {et~al.} 2011, \apj, 728, 46

[Weinberg} {et~al.}, 2013]
{Weinberg}, D.~H., {Mortonson}, M.~J., {Eisenstein}, D.~J., {et~al.} 2013,
\physrep, 530, 87

[White} {et~al.}, 2011]
{White}, M., {Blanton}, M., {Bolton}, A., {et~al.} 2011, \apj, 728, 126

[Zehavi} {et~al.}, 2005]
{Zehavi}, I., {Eisenstein}, D.~J., {Nichol}, R.~C., {et~al.} 2005, \apj, 621,

[Zehavi} {et~al.}, 2011]
{Zehavi}, I., {Zheng}, Z., {Weinberg}, D.~H., {et~al.} 2011, \apj, 736, 59

[Zheng} {et~al.}, 2007]
{Zheng}, Z., {Coil}, A.~L., \& {Zehavi}, I. 2007, \apj, 667, 760

[Zheng} {et~al.}, 2009]
{Zheng}, Z., {Zehavi}, I., {Eisenstein}, D.~J., {Weinberg}, D.~H., \& {Jing},
Y.~P. 2009, \apj, 707, 554

[Zheng} {et~al.}, 2005]
{Zheng}, Z., {Berlind}, A.~A., {Weinberg}, D.~H., {et~al.} 2005, \apj, 633, 791

\end{thebibliography}

\end{document}