William I Cowley

Karina I Caputi

Smaran Deshmukh

Matthew L N Ashby

Giovanni G Fazio

Olivier Le Fèvre

Johan P U Fynbo

Olivier Ilbert

Henry J Mccracken

Bo Milvang-Jensen

Rachel S Somerville

The Galaxy-Halo Connection for 1.5z5 as Revealed by the Spitzer Matching Survey of the

Keywords: methods, statistical -galaxies, evolution -galaxies, formation -galaxies, high-redshiftlarge-scale structure of Universe

1. INTRODUCTION A central ansatz of the Λ cold dark matter (ΛCDM) cosmological paradigm is that galaxies form from baryonic condensations within the potential well of a dark matter halo (e.g., White & Rees 1978). However, understanding the relationship between a dark matter halo and the galaxies it hosts is far from a trivial exercise. In particular, the issue this paper explores is how the stellar mass of a galaxy is related to the mass of its host dark matter halo, and how this relationship evolves with cosmic time.

A conceptually straightforward way in which to investigate the relation between galaxies and their host halos is a direct simulation of galaxy formation within a cosmological context. The current state of the art of such efforts comprises an N -body computation of the evolution of dark matter combined with either a hydrodynamical (e.g., Vogelsberger et Lacey et al. 2016) treatment of the baryonic processes involved. Both schemes have enjoyed considerable and groundbreaking successes in reproducing multiple observational datasets, thereby furthering our understanding of galaxy formation and evolution. However, they are both hampered by a lack of our knowledge regarding some of the complex physical processes involved, such as star and black hole formation, and associated feedback processes. This is largely related to the extreme computational challenge posed by attempting to resolve these processes in simulations of cosmologically significant volumes and time periods. As a result, some uncertainties regarding how accurately the physical processes involved in producing galaxies within dark matter halos are accounted for in these models remain.

Direct observational probes of the dark matter halos of galaxies include galaxy-galaxy lensing (e.g., Brainerd et al. 1996;Hoekstra et al. 2004; Leauthaud et al. 2010;Hudson et al. 2015) and the kinematics of satellite galaxies (e.g., Zaritsky et al. 1993 More et al. 2011). Galaxy-galaxy lensing uses distortions of the shapes and orientations of background galaxies caused by intervening mass along the line of sight and thus can be used to infer the foreground mass distribution. Satellite kinematics uses satellite galaxies as test particles in order to trace out the dark matter velocity field, and thus potential well, of the dark matter halo.

However, both of these techniques are mostly limited to low redshifts (z < 1), due to the difficulty of resolving individual galaxies at earlier epochs, and they rely on ensemble stacking of galaxies in order to extract a signal.

A simpler, though less direct, approach is to compare the observed abundance and clustering properties of galaxy samples with predictions from a phenomenological halo model (e.g., Neyman & Scott 1952;Berlind & Weinberg 2002;Cooray & Sheth 2002). This is a purely statistical description of how galaxies occupy halos and thus forgoes an understanding of the physical processes involved. However, this technique has been shown to provide a good description of the observed clustering of galaxies and can in principle be applied over very broad redshift ranges.

A drawback of this approach is that it relies on accurate a priori knowledge of the matter power spectrum, halo mass function, halo density profile, and the bias with which halos trace the underlying matter distribution. The relations used in practice are calibrated against numerical simulations, though there remains some doubt over the applicable regime for these calibrations. It also assumes that the bias with which halos trace the underlying matter distribution depends only on halo mass, ignoring potential effects such as 'halo assembly bias' (e.g., Zentner et al. 2014).

In recent years the halo occupation modeling technique has been utilized in many studies investigating the galaxy-halo connection (e.g., Zheng 2004;Zheng et al. 2007; Zehavi et al. 2011;Béthermin et al. 2014;Skibba et al. 2015), and some works have combined this technique with constraints from galaxy-galaxy lensing (e.g., Leauthaud et al. 2010;Coupon et al. 2015). However, interpreting the results from such studies in terms of physical galaxy properties is often complicated by the selection of the galaxies used. Many samples are selected by luminosity, resulting in uncertainty regarding the conversion to stellar mass.

Another difficulty is probing a sufficient range of halo masses such that the galaxy-halo relation can be meaningfully characterized. Increasing the depth of a survey, allowing lower mass galaxies to be investigated, often limits the survey area such that the number of massive galaxies is not sufficient for a robust clustering analysis. This trade-off has proven difficult to overcome for z 1.

Galaxy clustering at higher redshifts (up to z ∼ 7) has been investigated through use of the Lyman break 'dropout' selection method (e.g., Giavalisco et al. 1998; Harikane et al. 2016;Hatfield et al. 2017;Ishikawa et al. 2017) and Lyman α emission (e.g., Ouchi et al. 2017). However, as in these works galaxies are selected by broadband color and nebular emission respectively, the connection to properties such as stellar mass is even more uncertain.

Recently, McCracken et al. (2015) provided a clustering analysis of galaxies selected by stellar mass up to z ∼ 2 based on UltraVISTA DR1 near-infrared and ancillary COSMOS broadband data (McCracken et al. 2012). These authors were able to directly characterize the stellar-to-halo mass relation (hereafter SHMR) up to z ∼ 2.

Here, we use the unique Spitzer Matching survey of the UltraVISTA ultra-deep Stripes (SMUVS) galaxy catalog (Ashby et al. 2018;Deshmukh et al. 2018) to extend our understanding of the relationship between a galaxy's stellar mass and its host dark matter halo to higher redshifts than were previously possible. We do so through com-paring the observed clustering and abundances of stellar mass-selected samples of galaxies to predictions from a phenomenological halo occupation model. This paper is structured as follows: in Section 2 we describe the galaxy catalog used throughout this work and the construction of our galaxy samples. In Section 3 we describe the phenomenological halo model used to interpret our observations, and introduce a novel method to account for the use of photometric redshifts (this is discussed in more detail in Appendix A). We present our results in Section 4 (our main results are tabulated in Table 3, Appendix C). A brief discussion of some of the modeling assumptions made in our analyses is given in Section 5. We conclude in Section 6.

Throughout we assume a flat ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7, Ω b = 0.045, h = 0.7, σ 8 = 0.8, and n s = 0.95. All magnitudes are quoted in the Absolute Bolometric (AB) system (Oke 1974). (Taniguchi et al. 2007). The UltraVISTA data considered here correspond to the third data release 8 , which reaches an average depth of K s = 24.9±0.1 and H = 25.1±0.1 (2 arcsec diameter, 5σ). This paper forms part of a series of scientific studies that make use of SMUVS data, such as the search for strong Hα emitters at z = 4 -5 (Caputi et al. 2017) and the study of galaxy structural properties for z 5 (Hill et al. 2017).

A thorough description of the SMUVS multiwavelength source catalog construction and spectral energy distribution (SED) fitting is given in Deshmukh et al. (2018), however, we summarize the main details here.

Sources are extracted from the UltraVISTA HK s average stack mosaics using SExtractor (Bertin & Arnouts 1996). The positions of these sources were then used as priors to perform iterative point-spread function (PSF) fitting photometric measurements on the SMUVS 3.6 and 4.5 µm mosaics, using the daophot package (Stetson 1987).

For all of these sources, 2 arcsec diameter circular photometry on 26 broad, intermediate, and narrow bands U to K s is measured (Deshmukh et al. 2018). After cleaning for galactic stars using a B-J- [3.6] color selection (e.g., Caputi et al. 2011), and masking regions of contaminated light around the brightest sources, the final catalog contains ∼ 2.9 × 10 5 UltraVISTA sources with a detection in at least one IRAC band over an area of ∼ 0.66 square degrees.

The SED fitting is performed with all 28 bands (26 U through K s as well as Spitzer 3.6 and 4.5 µm) using the χ 2 minimization code LePhare 9 (Arnouts et al. Photometric redshifts and stellar mass estimates are obtained for > 99 percent of our sources. Using ancillary spectroscopic data in COSMOS to assess the quality of the obtained photometric redshifts, we found that the standard deviation, σ z , of |z phot -z spec |/(1+z spec ), based on ∼ 1.4 × 10 4 galaxies with reliable spectroscopic redshifts in the COSMOS field (see Table 1 Throughout we use the best-fit redshifts and stellar masses computed by LePhare.

Sample selection

We construct volume-limited stellar mass-selected samples of galaxies for this analysis. The stellar massredshift plane for the SMUVS catalog is shown in Figure 1. The 80 percent stellar mass completeness limit (based on 4.5 µm photometry) is calculated following the method described in Chang et al. (2013, see also Tomczak et al. 2014) and is shown as the gray line. All of our samples are above this completeness limit. The number of objects in each sample, as well as the median stellar mass, is summarized in Table 1.

We restrict our analysis to z > 1.

METHODS 3.1. The angular two-point correlation function

The angular autocorrelation function, w(θ), describes the excess probability, compared to a random (Poisson) distribution, of finding a pair of galaxies at some angular separation θ > 0. It is defined such that

δ 2 P 12 = η2 [1 + w(θ)] δΩ 1 δΩ 2 , ( 1
)
where η is the mean surface density of the population per unit solid angle, θ is the angular separation, and δΩ i is a solid angle element. Here, the angular correlation function is computed according to the standard Landy & Szalay (1993) estimator

w(θ) = DD(θ) -2DR(θ) + RR(θ) RR(θ) , (2)
where DD, DR and RR represent the number of datadata, data-random, and random-random pairs in a bin of angular separation respectively, and the random catalog is constructed to have the same angular selection as the data. We use a random catalog with ∼ 5 × 10 5 objects.

The errors on the two-point angular correlation function are estimated from the data using a jackknife approach (Norberg et al. 2009). We divide the SMUVS footprint into 60 approximately equal area jackknife regions, removing one region at a time we compute the covariance matrix as

C i,j = N -1 N N l=1 (w l i -wi ) × (w l j -wj ), ( 3
)
where N is the total number of jackknife regions, w is the mean correlation function ( l w l /N ) and w l is the estimate of w with the l th region removed. We compute w(θ) for each stellar mass threshold and redshift bin described in Table 1, for 12 evenly spaced logarithmic bins of angular separation in the range -3 < log 10 (θ/deg) < -0.6.

The halo occupation model

We use an established phenomenological halo occupation model (Zheng et al. 2007) to connect the observed galaxy angular correlation functions and abundances to host dark matter halos. In the halo model, the expected mean number of galaxies in a dark matter halo, N (M h ), the halo occupation distribution (HOD), is a sum of contributions from central galaxies, N c , and satellites, N s , such that

N (M h ) = N c (M h ) × [1 + N s (M h )]. (4)
The contribution from central galaxies (i.e., those at the center of the halo potential well) is modeled as a step function with a smooth transition:

N c (M h ) = 1 2 1 + erf log M h -log M h,min σ log M h , (5)
where M h,min is the mass at which 50 percent of halos host a single galaxy and σ log M h is the width of the central Note.

- (a) in log 10 (M /M ). For each stellar mass threshold and redshift bin we report the number of galaxies and the median stellar mass.

galaxy mean occupation. The contribution from satellite galaxies is modeled as a power law with a cutoff at low halo masses:

N s (M h ) = M h -M h,0 M h,1 αsat , (6)
where M h,0 is the cutoff mass scale, M h,1 characterizes the amplitude, and α sat describes the asymptotic slope at high halo mass. Thus, the HOD is described by five parameters: M h,min , M h,1 , M h,0 , α sat and σ log M h . A spatial correlation function, ξ(r), where r is the comoving spatial separation, is then constructed from the HOD 10 . This step assumes a number of relations, generally calibrated against numerical simulations, which are described in the following two paragraphs. The matter power spectrum is computed according to Eisenstein & Hu (1999) with the nonlinear correction of Smith et al. (2003) applied. Additionally, we implement the two halo exclusion model of Tinker et al. (2005), which improves on that presented in Zheng (2004). A thorough description of the construction of a similar halo model, upon which the one used in this work is based, is given in Appendix A of Coupon et al. (2012).

Given an HOD we can compute the following derived quantities which we will discuss later: the satellite fraction, f sat ,

f sat (z) = 1 -f cen (z) = 1 -N c (M h , z) n(M h , z) dM h / n gal (z), (7)
where n(M h , z) is the halo mass function and n gal is the 10 This is later projected for comparison with the measured angular correlation functions as described in Section 3.3.

galaxy number density,

n gal (z) = N (M h ) n(M h , z) dM h , (8)
and the effective large-scale galaxy bias, b gal ,

b gal (z) = b h (M h , z) N (M h ) n(M h , z) dM h / n gal (z), (9)

Projection and the effect of photometric redshift errors

In this work we measure angular correlation functions of galaxies in bins of redshift. This is a necessary consequence of our photometric redshifts, which are not accurate enough to measure the galaxy distribution in three dimensions. Thus a spatial correlation function, ξ(r), computed from an HOD, needs to be projected along the line of sight into two dimensions for comparison with our data. For this we use the Limber (1953) equation,

w(θ) = n gal (z) dV dz W (z) 2 dz dχ dz ξ(r, z) du n gal (z) dV dz W (z) dz 2 , (10
) where n gal is the number density of galaxies predicted by the HOD, dV /dz is the comoving volume element, dz/dχ = H 0 E(z)/c where E(z) = [Ω m (1 + z) 3 + Ω Λ ] 1/2 , and χ corresponds to the comoving radial distance to redshift z. The comoving line-of-sight separation, u, is defined by r = [u 2 + χ 2 2 ] 1/2 where 2 /2 = [ 1cos(θ)]. The W (z) term in Equation 10 relates to the redshift window that is being probed by the survey. If the redshifts were known precisely, then (ignoring further complications such as redshift space distortions) this would be a top-hat function equal to unity between the limits of the redshift range probed and zero elsewhere. However, with photometric redshifts, this is not the case, and the top-hat window should be convolved with an error kernel that is generally unknown a priori.

In order to mitigate this, here we assume a Gaussian error kernel and approximate a (1 + z) evolution in the error kernel dispersion, ∆ z , such that

W (z) = 1 2 erf z -z lo ∆ z (1 + z lo ) -erf z -z hi ∆ z (1 + z hi )
.

(11) Here z lo and z hi represent the lower and upper redshift limits of the photometric redshift bin respectively. Whilst the integral in Equation 10 is, in principle, over all redshift, it only has significant contributions from redshifts where W (z) is appreciably nonzero. We leave ∆ z as a free parameter in our fitting procedure, which is described in Section 3.4. This decision, and our modeling of the photometric redshift dispersion, is discussed in more detail in Appendix A.

Once projected according to Equation 10, we account for the integral constraint (e.g., Groth & Peebles 1977). This is the correction required as the measured angular correlation function will integrate to zero over the whole field, by construction of the estimator for w(θ) used. This results in the measured angular correlation function being underestimated by an average amount

σ 2 IC = 1 Ω 2 w true (θ) dΩ 1 dΩ 2 , (12)
where w true is the true angular correlation function and the angular integrations are performed over a field of area Ω.

Here we evaluate Equation 12according to the numerical method proposed by Roche & Eales (1999),

σ 2 IC = w true (θ) RR(θ) RR(θ) . ( 13
)
We take w true to be the angular correlation function predicted by our HOD model, and subtract σ 2 IC from it before comparing it to our observed correlation functions. We do this for each evaluation of w(θ) in our fitting procedure, which is described below. These corrections are typically < 10 -2 .

Fitting

We derive the best-fitting halo model corresponding to our clustering and abundance measurements using the python affine-invariant implementation for Markov chain Monte Carlo (MCMC), emcee 11 (Foreman-Mackey et al. 2013). Although our model in principle has six adjustable parameters, after some initial tests we determined that, despite a good signal being measured in our correlation functions, our data were insufficient to fully constrain all of them. This was based on consideration of the Bayesian information criterion (BIC ;Schwarz 1978; see e.g., Liddle 2007 for a discussion of this criterion). It is an approximation of the Bayesian evidence and is similar in construction to the Akaike (1974) information criterion. We therefore decided to fix a number of our model parameters to avoid overfitting our data. For M h,0 , we follow Conroy et al. (2006) There are thus three free parameters in our model: M h,min , M h,1 and ∆ z . We discuss our decision to leave ∆ z as a free parameter in more detail in Appendix A. We simultaneously fit both the observed clustering and number of galaxies by summing both contributions to the total χ 2 such that

χ 2 = i,j [w obs (θ i )-w mod (θ i )](C -1 ) i,j [w obs (θ j )-w mod (θ j )] + [N obs gal -n mod gal V mod ] 2 σ 2 √ N + σ 2 CV + σ 2 fit , (15)
where n mod gal is the number density of galaxies predicted by the HOD and the volume, V mod , is computed according to V mod = (Ω/4π) W (z) (dV /dz) dz where Ω is the area of our survey and W (z) is defined as in Equation 11. This therefore incorporates the effect of photometric redshift error on the observed number of galaxies through the ∆ z parameter. There are contributions to the error on the observed number of galaxies from Poisson noise, σ √ N , cosmic variance, σ CV , which is computed using the method presented in Moster et al. (2011), and a term that accounts for the error in the SED fitting procedure, σ fit . For this latter term we construct 100 mock catalogs by scattering each SMUVS galaxy within the probability distribution for its stellar mass and redshift given by LePhare (the construction of these mock catalogs are described in more detail by Deshmukh et al. 2018). We then apply the same photometric redshift and stellar mass cuts to these mock catalogs and take σ fit to be the standard deviation in the number of galaxies in each in each sample across the 100 mock catalogs. The value of these different sources of error for each sample is given in Table 3, the total error on the observed number of galaxies is typically ∼ 10 percent. This is comparable to the uncertainty quoted by Wake et al. (2011), who derived errors of ∼ 15 percent on their galaxy number densities using a different method for a slightly smaller area than surveyed here (∼ 0.4 deg 2 ). Additionally, we have scaled the inverse covariance matrix, C -1 , (see Equation 3) in order to account for bias from a finite number of jackknife samples, according to Hartlap et al. (2007).

We assume uninformative (i.e., flat) priors for our three free parameters (our fitting procedure is thus analogous to a maximum likelihood estimation). After a conservative 'burn-in' phase we run 40 walkers for 2500 loops, which results in the posterior distribution being sampled with 10 5 points. The chains appear 'well-mixed' by this point, and inspection of their autocorrelation and the Gelman-Rubin statistic (Gelman & Rubin 1992) indicates that they have converged. Our best-fit parameters are taken to be the sample in our chains that returns the minimum χ 2 , and the uncertainties represent the bounds of the 1σ contours in parameter space described by χ 2 < χ 2 min + ∆χ 2 , where ∆χ 2 = 3.53 for three free parameters. The value and uncertainties of the derived quantities (e.g., satellite fraction) are determined in the same way. Some examples of the likelihood distributions produced by our modeling and fitting procedure are shown in Figure 17, Appendix D.

RESULTS

In this section we display our main results. In Section 4.1 we present some of our measured angular correlation functions and the resulting best-fit halo models. We present the characteristic halo masses of, and satellite galaxy fractions derived from, our best-fit halo models in Sections 4.2 and 4.3 respectively. In Section 4.4 we show the clustering amplitude (comoving correlation length and large-scale bias) evolution inferred from our data, and in Section 4.4.1 we use our bias measurements to propose a technique for computing coevolutionary histories for dark matter and stellar mass assembly. Finally, in Section 4.5 we investigate the SHMR of our data. We show that with our derived SHMRs we can reproduce our stellar mass functions in Appendix B, highlighting the consistency of our analysis. Our main results are summarized in Table 3 in Appendix C.

Halo model analysis of angular two-point correlation functions

Here we present our measurements of the clustering of galaxies and results from the corresponding best-fit halo model. Some examples of these are shown in Figure 2.

In general, the fits to the observed clustering appear to be good. However, a purely visual inspection of the fit may be misleading. This is because the covariance between different angular bins and the abundance constraint used in the fitting are not visualized in Figure 2. Our goodness of fit is also reflected in our low values for the minimum reduced χ 2 (typically ∼ 1, see Table 3).

It should be noted that the amplitudes of the measured angular correlation functions do not necessarily display a monotonically increasing trend with increasing stellar mass, evident in the 1.5 < z < 2.0 panel of Figure 2. However, as higher stellar mass galaxies are less abundant one would expect them to reside in higher mass halos12 that are more biased and would thus have a greater correlation amplitude.

In our framework, this is accommodated by the ∆ z parameter in our fitting procedure. For example, at 1.5 < z < 2.0 the log 10 (M /M ) > 9.8 sample has a lower clustering amplitude on large scales ( 5 × 10 -3 deg) than the log 10 (M /M ) > 9.2 sample (see the top left panel of Figure 2). Therefore, the higher mass sample must have a higher best-fit value of ∆ z . Our modeling produces a value of ∆ z ∼ 0.09 relative to ∆ z ∼ 0.04 for the highand low-mass samples mentioned above respectively. We caution against overinterpreting these values in terms of |z photz spec |/(1 + z spec ) diagnostics, as our ∆ z parameter describes a global photometric redshift dispersion. It does not segregate out redshift outliers and ignores selection effects introduced by comparing to spectroscopic redshifts. A higher value of ∆ z projects the intrinsic spatial correlation function over a larger redshift range (according to Equation 10) and thus lowers the amplitude of the angular correlation function. This is also reflected in the amplitudes of the matter correlation functions, shown as the dashed lines in Figure 2.

Generally, though it is not always the case, we find that the samples with higher stellar mass in a given redshift bin tend to have higher values of ∆ z . This is not necessarily unexpected. Massive galaxies (M 10 10 M ) exhibit, on average, greater dust attenuation than less massive ones. This, in turn, leads to an increased degeneracy in the parameter space of the SED fitting between stellar age and dust reddening 13 . In general, however, we find that our best-fit values of ∆ z are low relative to the size of our photometric redshift bins, suggesting that our redshifts are accurate enough for the analysis performed in this study. Typically, we find ∆ z (1+z 50 )/(z hi -z lo ) 0.5 (where z 50 is the median redshift of the sample and z hi and z lo represent the limits of the bin), and all of our values for ∆ z (1 + z 50 )/(z hiz lo ) are consistent with being lower than unity (within 1σ, see Figure 14 and the discussion in Appendix A for more details).

Characteristic halo masses

Here we consider the characteristic halo masses, M h,min and M h,1 , which represent the mass at which 50 percent of halos host a single central galaxy in the sample, and at which each halo hosts an additional satellite galaxy respectively 14 . These quantities are shown in Figure 3 as a function of the median stellar mass of each sample, and redshift.

We can see that both halo masses are well constrained by our clustering and abundance measurements and that they form tight, approximately linear, relationships with stellar mass (in log space) with no significant evolution with redshift (though our 4.0 < z < 5.0 point is offset to lower halo masses, as is expected from the hierarchical nature of structure formation). We find the 'mass gap' between M h,min and M h,1 to be ∼ 10-20 which is broadly consistent with earlier studies (e.g., Zehavi et

Satellite fractions

Satellite galaxies are those that do not sit at the center of the potential well of a dark matter halo (as 'central' galaxies do), but rather formed in other dark matter halos that later merged with (or were accreted onto) their current host. They are on bound orbits, which decay due to the dynamical friction of dark matter acting on the sub-halo, and will eventually merge onto the central galaxy. Thus an increase in the satellite fraction can indicate halo merging activity, whilst a decrease can indicate satellites merging with the central galaxy.

The satellite fractions derived from our clustering and abundance measurements (Equation 7) are shown in Figure 4. We find that our inferred satellite fractions generally increase with cosmic time, as can be expected from 13 Interestingly, our values for ∆z are lower for 2.0 < z < 3.0 than 1.5 < z < 2.0, as our redshift bin is broader and the constraint of the age of the Universe (and thus of the stellar populations) at these higher redshifts reduces this degeneracy (e.g. Thomas et al. 2017). 14 In fact, the halo mass at which the HOD is equal to 2 is not exactly M h,1 , due to the cutoff halo mass, M h,0 , and the scaling of the satellite galaxy distribution by that of the central galaxy distribution. However, for our purposes here this distinction is unnecessary. hierarchical structure formation. This is supported by studies at lower redshift that indicate higher satellite fractions than we observe here, for example, Zehavi et al.

(2011) derive a value of f sat ≈ 0.3 at z = 0. Additionally, we observe a mild decreasing trend of f sat with increasing stellar mass for 2.0 < z < 3.0, but a much flatter relationship for 1.5 < z < 2.0. This is indicative of halo merging activity (which converts central galaxies into satellites) and/or in-situ stellar mass assembly in satellite galaxies (which moves satellite galaxies into higher stellar mass samples). However, distinguishing the relative contribution of these two potential processes would require further analysis of these data with more physically motivated models than the HOD framework adopted here.

Our data compare favorably with the satellite fractions inferred by McCracken et al. (2015). These authors performed a similar halo model analysis as performed in this work, though their initial photometry was based on the earlier UltraVISTA DR1 (McCracken et al. 2012). Additionally, they left all of the HOD parameters as free in their fitting procedure (whereas some are fixed here as was discussed in Section 3.4) and did not consider photometric redshift dispersion in their model, equivalent to

log 10 (M h /M ) M h,min M h,1 1.5 < z < 2.0 2.0 < z < 3.0 3.0 < z < 4.0 4.0 < z < 5.0 M h,min M h,1

Comoving correlation length and galaxy bias

Earlier studies of the clustering of galaxies at low redshift found that their correlation function could be adequately described by a power law, i.e., ξ(r) = (r/r 0 ) -γ (e.g., Davis & Peebles 1983;Hawkins et al. 2003). Here r 0 is the 'correlation length' that describes the separation at which the correlation function is equal to unity, and γ, which is typically ∼ 1.8, describes the slope.

However, this description of the correlation function of galaxies is somewhat problematic as more recent observations, at higher redshifts, have shown that galaxy clustering exhibits deviations from a simple power law that can be readily understood in terms of halo occupation models such as those used in this work (e.g., Zheng 2004). Moreover, the results from fitting a simple power law can be sensitive to the scales used in the fitting procedure and whether γ was fixed or not. However, it can still be informative, both for comparing with earlier work and as it roughly describes the amplitude of the observed correlation function.

Here, rather than attempt to fit a power law to our observed angular correlation functions, we compute this quantity directly from the spatial correlation functions predicted by our best-fit halo models. For this we define the correlation length such that ξ(r 0 ) ≡ 1. We show our inferred r 0 values in Figure 5. Here, we can see that r 0 generally increases monotonically with redshift and stellar mass. For example, at 1.5 < z < 2.0 galaxies with M ∼ 10 10.3 M display a correlation length of ∼ 5.6 h -1 Mpc, whereas at 4.0 < z < 5.0 the correlation length for galaxies of a similar stellar mass is ∼ 7.9 h -1 Mpc.

Comparing to some earlier measurements in Figure 5 9 we can see that our data occupy a hitherto relatively unexplored region of this parameter space. Marulli et al. (2013) fit a power law to the projected spatial correlation function measured from galaxies identified in the VIMOS public extragalactic redshift survey (VIPERS, Guzzo & The Vipers Team 2013) and Durkalec et al. (2015) performed a similar analysis with data from the VIMOS ultra-deep survey (VUDS). Both of these VIMOS surveys allow for the determination of spectroscopic redshifts. The Marulli et al. data appear broadly consistent with the trends identified in ours, though the redshift range of the two studies does not overlap. They find the same trend of increasing r 0 with stellar mass and their correlation lengths are at smaller scales than ours (for a fixed stellar mass) which would be expected as they are measured at lower redshift. At lower stellar masses (M ∼ 10 9.6 M ) we find a larger r 0 than Durkalec et al., though this may be due their r 0 being determined through a power-law fit, rather than a halo model as we do here, as the large-scale bias Durkalec et al. inferred from fitting a HOD model to their clustering measurements agrees well with ours, as is shown in Figure 6. Wake et al. (2011) find similar correlation lengths to ours for high-mass galaxies, though they derive their correlation lengths by fitting a power law to their observed angular clustering, which is different to the procedure we follow.

A perhaps more physically meaningful quantity, related to the comoving correlation length, is the largescale galaxy bias. This describes the difference in amplitude between the matter correlation function and that of the observed galaxy sample. We compute an effective large-scale bias from our best-fit halo models, according to Equation 9. These are shown in Figure 6 as a function of stellar mass and redshift. A striking monotonically increasing relationship between redshift and large-scale bias (at a fixed stellar mass) is evident here. The relationship between stellar mass and bias appears to be flatter at lower redshifts (z 3.0), and at lower stellar masses (though we note that the stellar mass range probed also changes with redshift).

A similar relationship between large-scale bias and stellar mass as is also seen at z ∼ 1.5 in the Martinez 3) appear lower than those we infer for 2.0 < z < 3.0, however, this may be due to the broader redshift range we consider, which would include more biased galaxies in the 2.5 < z < 3.0 range (for a fixed stellar mass threshold) than the McCracken et al. data. We also note that the halo model-derived bias values found by Wake et al. (2011) for z ∼ 1.6 -2.2 are higher than ours for 1.5 < z < 2.0. They find b gal ∼ 3 -4 for galaxies with log 10 (M /M ) > 10.5.

Again, this comparison with literature values indicates that we have filled in a large hitherto unexplored area of this plane, and our results show reasonable agreement where they coincide with earlier measurements.

Potential evolutionary paths

In this section, we propose a method to track the coevolution of dark matter halo and stellar mass assembly, based on our large-scale bias and stellar mass measurements presented above.

First, we generate descendant halo masses, according to our z = 0 halo mass function. To these we apply halo mass assembly histories, using the form proposed by McBride et al. (2009). These authors parametrized halo mass accretion histories from the Millennium Nbody simulations (Springel et al. 2005) according to

M h (z) = M h,z=0 (1 + z) β e -γz , (16)
where M h,z=0 is the mass of the z = 0 descendant. We note that other parametrizations for halo mass accretion histories in ΛCDM have been proposed in the literature (e.g., van den Bosch 2002; Wechsler et al. 2002). We generate a realistic ensemble of these histories according to the distributions of β and γ described in Appendix A of McBride et al. We then select halo histories with a bias equal to the effective galaxy bias we measure for a given galaxy population at a given redshift and track the evolution of this bias over cosmic time. We use the halo Using our galaxy population with log 10 (M /M ) > 10.8 at 1.5 < z < 2.0 as our starting point, we show the median evolution of the bias for such histories as the solid dark line in Figure 7. The dashed lines indicate the 16-84 percentile scatter which we can compute from the (∼ 10 3 ) halo histories that go through this locus on the bias-redshift plane. At high redshift, the scatter in the halo histories encompasses a broad range of progenitor halo masses.

For stellar mass evolution, we apply the cumulative number density approach proposed by Behroozi et al. (2013a). These authors found, using the abundance model of Behroozi et al. (2013b), that tracking the progenitors of galaxies selected by their stellar mass at a given redshift could be done simply by increasing the cumulative number density of the population by 0.16 × ∆z dex at previous redshifts 15 . Behroozi et al. attributed this simple power-law behavior to the roughly constant halo merger rate per unit halo per unit ∆z in ΛCDM cosmologies (e.g., [START_REF] Fakhouri | [END_REF]).

We do this now for our SMUVS data, starting again with log 10 (M /M ) > 10.8 galaxies at 1.5 < z < 2.0. At earlier redshifts, we relate the progenitor number density (based on the Behroozi et al. argument outlined above) to a stellar mass threshold and a median stellar mass. The evolution of the median stellar mass is shown as the colored line, plotted on top of the median halo mass accretion history described above, in Figure 7. We track the stellar mass evolution to z ∼ 4, at which point the stellar mass threshold of the progenitors falls below our 80 percent stellar mass completeness limit. Of interest is where our median halo history intersects with our measurements of galaxy bias at higher redshifts. Here, our median stellar masses computed via the progenitor number density method outlined above appear to agree broadly with those measured directly at that redshift for galaxy samples with a similar bias to our computed median halo mass accretion history. This is a potentially powerful result, as it indicates that clustering measurements of stellar mass-selected galaxy samples, such as those performed here, can be used to project consistent evolutionary paths for the dark matter and stellar mass assembly over significant proportions of the history of the Universe. It is also more straightforward to compute than other methods used to propose evolutionary pathways for galaxies based on clustering measurements (e.g., Conroy et al. 2008). We present this result again in the left panel of Figure 8, where the assembly histories are instead plotted as a function of lookback time and the bias measurements from Figure 7 have been converted into halo masses using the relation of Tinker et al. (2010). In the right panel of Figure 8 we show another example, applying this technique to galaxies with log 10 (M /M ) > 10.2 at 1.5 < z < 2.0.

Also of note here is that the assembly histories for stellar mass and dark matter are intrinsically different shapes. This is a challenge for most physical galaxy formation models to reproduce, as is discussed in detail by Mitchell et

al. (2014).

There are a number of caveats to the result in this section, and it should be noted that the agreement between the proposed assembly histories and our data is not perfect and that the scatter on both the halo and stellar mass histories is significant. Choosing halo histories based on a single bias value ignores the fact that a range of halo masses will be occupied, which is indeed implied by our halo occupation models. This would probably only increase the (already significant) scatter and have a relatively minor impact on the median history, as there is a steep fall-off of the halo mass function toward high halo masses and in our HODs toward low halo masses. The halo histories we generate are based on simulations performed with slightly different cosmological parameters than those assumed here16 , however, any change due to this is likely to be sub-dominant to the intrinsic scatter in the histories, which we have explicitly shown.

Nevertheless, we consider the agreement good enough for this technique, combining ΛCDM halo mass accretion histories and cumulative number density arguments for stellar mass evolution, based on stellar mass-selected clustering measurements, to be explored further in order to aid our understanding of the coevolution of the stellar and dark matter content of halos.

The stellar to halo mass relation

The SHMR is of great importance to galaxy formation models as it can be interpreted as the star formation history integrated over the lifetime of the halo, and thus as the efficiency with which baryons are converted into stars in halos of a given mass (assuming a constant fraction of a halo's mass is composed of baryons).

We are able to constrain this relation for our two lowest redshift bins (1.5 < z < 2.0 and 2.0 < z < 3.0) as shown in Figure 9, where we plot the ratio of the median stellar mass of a galaxy sample to the characteristic halo mass M h,min . The error bars shown here have been propagated from the uncertainty on M h,min derived from our MCMC fitting procedure. At higher redshifts, our clustering measurements lack the range of stellar/halo masses required to adequately characterize this relationship.

At 4.0 < z < 5.0 our data suggest an increase in the normalization of the SHMR. While this is a tentative result, we interpret this as subsequent stellar mass assembly being regulated by feedback processes, whilst the growth of a dark matter halo is not, reducing the normalization of this ratio over cosmic time. It also points toward the star formation efficiency being greater at higher redshifts. Note. -M peak is not a fitted parameter, but is derived from the others.

We describe the SMHR using the following functional form, proposed by Behroozi et al. (2013b; see their Equation (4)):

log 10 (M (M h )) = log 10 (εM t ) + f log 10 M h M t -f (0), (17)
where the function f (x) is described by

f (x) = -log 10 (10 αx + 1) + δ (log 10 (1 + exp(x))) γ 1 + exp(10 -x ) . ( 18
)
That is, a power law with slope -α for M h M t and an exponential transitioning into a sub-power law with slope γ for M h > M t . Other parameterizations of this relationship have been proposed, mainly four-parameter double power-laws with different high-and low-mass slopes (e.g., Yang et al. 2003;Moster et al. 2010), however Behroozi et al. found that the five-parameter form performed better based on their abundance-matching technique. We fit this relation to our data, with the best-fit parameters being given in Table 2; the 1σ uncertainty of the fit is shown by the gray regions in Figure 9. In Appendix B we show that we can accurately reconstruct the observed SMUVS stellar mass function using these bestfit relations, highlighting the consistency of our analysis. At z 3.0 our results would be complemented by combining them with a consistent analysis of data from larger and/or deeper surveys in order to fully resolve the SHMR. The peaked nature of the SHMR is commonly interpreted as the integrated effect of stellar and AGN feedback processes inhibiting star formation in the low-and high-mass regimes respectively within the hierarchical structure formation of ΛCDM (e.g., Benson et al. 2003;Bower et al. 2006). Given its peaked nature, the halo mass at which the SHMR is maximized, M peak , is of further interest. This indicates the halo mass at which the conversion of baryons into stars over the history of the halo has been most efficient. We derive this from the maximum of our best-fit SHMRs and give our values for M peak in Figure 10 (see also the rightmost column of Table 2). We can see that our estimations link up smoothly with previous estimates from similar studies at lower redshifts. Though we note that the Zehavi et al. (2011) and Coupon et al. (2012) studies are based on luminosityselected samples, rather than stellar mass-selected.

To describe the mild redshift dependence of M peak we fit a simple power law to the observational data shown in Figure 10 (e.g. Moster et al. 2010), and find that log 10 (M peak /M ) = 11.6 × (1 + z) 0.06 (shown as the solid line), is a reasonable description of the data, though this is of course somewhat dependent on the observational data chosen.

Our results for M peak are consistent with those from the sophisticated abundance matching technique of Behroozi et al. (2013b), also shown in Figure 10. This model constrains the SHMR based on halo merger trees from the Bolshoi simulations (Klypin et al. 2011) and observational estimates of the stellar mass function, the cosmic star formation rate density and the specific star formation rate -stellar mass plane for 0 < z < 8. Here we compute their value for M peak based on their best-fit SHMR (see their Section 5) and propagate through their uncertainties on the SHMR parameters. Their model agrees well with other observational estimates of M peak at low redshift (z 1), and is consistent with our measurements at higher redshifts, though understandably their model is less well constrained here due to the availability of observational data at these redshifts that could be used in their fitting procedure.

DISCUSSION OF MODELING ASSUMPTIONS

In this section, we briefly discuss some issues pertaining to the validity of some of the assumptions we have made in our analyses.

Throughout, we have assumed a five-parameter functional form for the HOD that is motivated by results from semi-analytical and hydrodynamical simulations (e.g., Zheng et al. 2005Zheng et al. , 2007)). As is discussed in Guo et al. (2016), this follows from assuming a lognormal distribution for central galaxy stellar mass at a fixed halo mass, and a power-law relation between the mean stellar mass of central galaxies and their host halo masses. In regimes where the SHMR deviates significantly from a power law, this halo occupation may not hold. Leauthaud et al. (2011) investigated this and found that changing the form of their assumed SHMR affected their best-fit HOD parameters only within their 1σ uncertainty. Thus while introducing different functional forms for the HOD may change the interpretation of some of the parameters the modeling results should not be significantly affected. We note that other forms for the galaxy HOD have been proposed in the literature (e.g., Geach et al. 2012; Gonzalez-Perez et al. 2018), though these are often designed to model populations that have been selected by properties that trace star formation, rather than stellar mass.

We include an additional free parameter, ∆ z , in order to encapsulate the dispersion due to the use of photometric redshifts in our analysis. Here we have assumed that this results in a Gaussian dispersion that evolves linearly with (1+z), leading us to assume the window function in Equation 11 when projecting our spatial correlation functions according to Equation 10. These assumptions are based on comparisons of our 'best-fit' LePhare photometric redshifts, with literature spectroscopic redshifts in COSMOS. The distribution of the resulting (z photz spec )/(1 + z spec ) roughly resembles a Gaussian. However, there is a small tail of outliers that we do not account for explicitly in our analysis. Benjamin et al. (2010) present a method that can in principle be used to estimate the impact that outlying photometric redshifts can have on clustering analyses. However, this method does not provide a unique solution for redshift bin contamination due to photometric redshift outliers because of degeneracies in the parameter space. We therefore refrain from implementing the method of Benjamin et al. and note that we do not place any prior constraint on the value of ∆ z determined by our fitting procedure. In addition, these outliers are a very small proportion (∼ 5 percent) of the sample for which we have reliable spectroscopic redshifts, so we expect the impact of these on our science results to be minimal and certainly subdominant to the more general redshift dispersion which we have accounted for. In Appendix A we show that we can reasonably reproduce the observed photometric redshift distribution, i.e., one constructed from the best-fit photometric redshifts computed by LePhare, from our best-fit values of ∆ z , indicating the consistency of this method, and that values of ∆ z (1 + z 50 )/(z hiz lo ) are typically 0.5.

In implementing our HOD model we make a number of further assumptions. One of which is that the distribution of satellite galaxies traces that of the dark matter profile within a halo, which is assumed to be of the Navarro et al. (1997) form. This assumption may affect our conclusions regarding the satellite fraction of galaxies. Again, it is motivated to some degree by hydrodynamical simulations (e.g., Nagai & Kravtsov 2005) though some studies have shown that this can be sensitive to the implementation of baryonic physics (e.g., Simha et al. 2012).

The halo bias parameterization we are using (Tinker et al. 2010) is calibrated against N -body simulations for z 2.5, so there is some uncertainty about the application to the higher redshifts studied here. Additionally, more detailed modeling of scale-dependent bias (e.g., Angulo et al. 2008) may be required at higher redshifts. For example, nonlinear clustering (e.g., Jose et al. 2016) can provide a better fit to galaxy clustering at high redshifts (z = 3 -5, Jose et al. 2017).

Finally, the halo occupation approach used here implicitly assumes that galaxy occupation statistics depend solely on the mass of the host dark matter halos. Whilst it is generally accepted that this is the halo property that most significantly influences the properties of the galaxies they host, over the last decade or so a body of evidence has emerged suggesting that it may also depend on additional properties such as the halo formation time (e.g., Gao et al. 2005;Croton et al. 2007;Zentner et al. 2014). This phenomenon is often referred to as 'halo assembly bias', and adjustments to the standard HOD formalism in order to account for this have been proposed (e.g., Hearin et al. 2016). However, it is still not clear if assembly bias is a significant effect, either in hydrodynamical simulations (e.g., Chaves-Montero et al. 2016), or observations (e.g., Lin et al. 2016), so we do not consider it here.

CONCLUSIONS

In this work we have applied a phenomenological halo model to the observed clustering and abundance of galaxies in the Spitzer Matching survey of the UltraVISTA ultra-deep Stripes (SMUVS) to understand the connection between galaxies and their host dark matter halos. SMUVS provides a unique combination of large area (∼ 0.66 square degrees of the COSMOS field) and unparalleled depth at 3.6 and 4.5 µm. This, combined with the latest ultra-deep UltraVISTA data and other ancillary data in the COSMOS field, allows for precise photometric estimates of redshift and stellar mass and robust statistics from such a large catalog (∼ 2.9 × 10 5 objects). As a result, our analysis is performed over an unprecedented redshift range (1.5 < z < 5.0) for volume-limited galaxy samples selected by stellar mass.

To interpret the observed clustering we utilize a wellestablished five-parameter halo model, though we fix the values of three of these parameters. Novel to our approach is the inclusion of the photometric redshift dispersion, ∆ z , as a free parameter in the fitting procedure. This parameter is designed to encapsulate the effect of dispersion in the photometric redshifts used in our anal-ysis and allows us to fully exploit the clustering information in our photometric galaxy catalog. This parameter is discussed in more detail in Appendix A.

From our best-fit halo models we can derive a number of interesting quantities. The characteristic halo masses, M h,min and M h,1 , are well constrained by our data and exhibit fairly tight linear relationships (in logarithmic space) with stellar mass. This is similar to what has been found in earlier studies with similar analyses (e.g., McCracken et al. 2015).

The satellite fraction, i.e., the fraction of galaxies in a sample that do not sit at the center of the potential well of their halo, generally increases with cosmic time. This is in broad agreement with other studies in the literature and in line with expectations of hierarchical structure formation.

The comoving correlation length, r 0 , [defined such that ξ(r 0 ) ≡ 1], and large-scale effective galaxy bias, b gal , show strong monotonic trends with increasing redshift for a fixed stellar mass over the entire redshift range probed.

Projecting ΛCDM halo mass accretion histories computed according to Appendix A of McBride et al. (2009) through our measurements of the large-scale bias, and combining this with the evolution of stellar mass content based on the cumulative number density (Behroozi et al. 2013a), we propose an evolutionary path for the dark matter and stellar mass content of galaxies for 1.5 < z < 4.0. This shows a favorable agreement with our independent determinations of the relationship between galaxy bias and stellar mass at these redshifts and appears to be an interesting method to investigate the coevolution of these quantities in the future.

Finally, we investigate the SHMR. This indicates the efficiency with which baryons have been converted into stars over the lifetime of a dark matter halo. We fit the SHMR parametrization of Behroozi et al. (2013b) for our two lowest redshift bins (1.5 < z < 2.0 and 2.0 < z < 3.0). In Appendix B we show that these SHMRs can reproduce the stellar mass function computed directly from the SMUVS catalog, highlighting the consistency of our analysis. The halo mass at which this relationship peaks indicates the halos for which the integrated star formation efficiency is greatest. We conclude that this quantity evolves mildly with redshift, i.e., log 10 (M peak ) ∝ (1 + z) 0.06 for z 3. Our inferred values of M peak are consistent with independent constraints on the evolution of this quantity inferred from abundancematching techniques (Behroozi et al. 2013b).

We end by acknowledging that the results presented in this work could be complemented in the future by deeper observations and/or larger area surveys. This would allow a number of trends/relationships identified here (and in earlier studies) to be fully investigated at higher redshift, which will provide important constraints for physical models of galaxy formation and evolution. Conversely, further modeling of our data with more physically motivated models than the HOD framework used in this study would allow us to disentangle the physical processes responsible for the trends observed. In this appendix, we discuss our modeling of the redshift dispersion caused by using photometric redshifts, which we describe through the parameter ∆ z and Equation 11. This assumes that the distribution of (z photz true ) is Gaussian with a width that scales with (1 + z true ). First, we show that under these assumptions, our modeling can accurately recover the redshift dispersion. To do this we use a simulated 2 deg 2 galaxy catalog from the semi-analytical model of galaxy formation, GALFORM (Lacey et al. 2016). We simulate photometric redshifts by applying a Gaussian scatter that scales with (1 + z true ) and then select galaxies based on their photometric redshifts in the range 2 < z phot < 3 and log 10 (M /M ) > 10.4. We do this for two values of standard deviation for the Gaussian scatter, ∆ z , of 0.02 and 0.08. For simplicity, we do not apply any scaling or error to the simulated stellar masses. We then measure the angular clustering and fit a HOD model to this following the same procedure we use for our SMUVS data as outlined in Section 3, for this we assume a 10 percent error on N gal as is typical of our SMUVS data. The resulting marginalized likelihood distributions from our MCMC fitting are shown in Figure 11. We can see from the bottom right panel that the best-fit MCMC value for ∆ z is remarkably close to the input value, and the other panels show that the impact on the other two halo model parameters, M h,min and M h,1 , is minimal.

The effect of increasing the photometric redshift dispersion on the measured clustering is shown in the left panel of Figure 12. We can see that increasing the value of ∆ z reduces the amplitude, as the clustering is being projected over a broader redshift range. However, the model is able to account for this such that the best-fit HOD that is recovered is the same, as can be seen in the inset panel. In the inset panel we also show the HOD predicted by the GALFORM model. We note that this is not precisely reproduced by the statistical HOD model, but consider further investigation of this beyond the scope of this work. The photometric redshift broadening is shown further in the right panel of Figure 12 for the case of ∆ z = 0.08. We see here that the true redshift distribution of galaxies selected by 2 < z phot < 3 is actually somewhat broader than these photometric redshift limits suggest and that this can be reasonably recovered by the HOD model, reminding the reader that the HOD model cannot account for evolution of the population within the redshift bin chosen. We conclude from these tests that our methodology can accurately account for the redshift dispersion arising from our necessary use of photometric redshifts.

We now discuss our decision to leave ∆ z as a free parameter. In principle, ∆ z could be measured from the data by comparison with spectroscopic redshifts. However, the spectroscopic coverage of our catalog is small (∼ 5 percent) and has a complicated selection function. Nevertheless, our data suggest σ z = 0.035 for galaxies in the redshift range covered by this study (1.5 < z < 5.0). We remind the reader that σ z is the standard deviation of |z photz spec |/(1 + z spec ), excluding outliers [objects with |z phot -z spec |/(1+z spec ) > 0.15]. We therefore fit a set of HOD models fixing ∆ z = 0.035 and compare the resulting best-fit models to our fits with ∆ z left as a free parameter. To do this, we use the Bayesian information criterion (BIC; Schwarz 1978; see e.g., Liddle 2007 for a discussion). This criterion is evaluated as follows:

BIC = -2 ln L max + k ln N, ( A1
)
where k is the number of model parameters and N is the number of data points (12 angular clustering data points and the observed number of galaxies). It is an approximation of the Bayesian evidence that should be a reasonable one given our uninformative priors. Models with the lowest value of BIC are favored, and though models with more free parameters can generally achieve a greater maximum likelihood, they are penalized by the k ln N term in Equation A1.

We consider the difference in BIC values, ∆BIC, between our ∆ z fixed and ∆ z free models in the left panel of Figure 13. As a 'rule of thumb', values of ∆BIC > 6 indicate that a model is 'strongly' favored (e.g., Appendix E of [START_REF] Fabozzi | Appendix E: Model Selection Criterion: AIC and BIC[END_REF]). Thus, we can see that in many cases the models in which ∆ z is a free parameter are strongly favored. Additionally, we assess the impact of ∆ z on the two remaining halo model parameters. The fractional change in these is typically 0.2 dex, as is shown in the right panel of Figure 13. Interestingly, a higher value of ∆ z appears to produce slightly higher values of M h,min (which will lead to higher values of derived quantities such as b gal and r 0 ) and slightly lower values of M h,1 (which will lead to higher values of f sat ). Given that ∆ z -free models are strongly favored by the BIC and the impact on the other model parameters is fairly small ( 0.2 dex) we are justified in our decision to leave ∆ z as a free parameter. Finally, we discuss whether our modeling returns reasonable values for ∆ z . In Figure 14 we show our values of ∆ z (1 + z 50 )/(z hiz lo ), which gives a measure of the redshift dispersion relative to the width of the redshift bin we are considering (here z 50 is the median stellar mass of a sample and z hi and z lo represent the limits of the photometric redshift bin). The values of this quantity are typically 0.5, which suggests that the bin-to-bin photometric redshift contamination in our work is minor once the general redshift dispersion has been accounted for. We note that our 1.5 < z < 2.0, log 10 (M /M ) > 10.4 sample appears to return an anomalously high value of ∆ z , however, none of our scientific conclusions would be changed by removing this sample from our study. Our sample at 4.0 < z < 5.0 also returns a high value for ∆ z , which we accept as a consequence of the extreme high-redshift nature of this sample in terms of stellar mass-selected clustering measurements.

A further check on the values of ∆ z returned by our model is to use them to reconstruct the photometric redshift distribution of our catalog for a given stellar mass threshold. We show this in Figure 15 over the redshift range 1.5 < z < 3.0 for two stellar mass thresholds. These appear to be a reasonable reconstruction of the redshift -Left panel: difference in the Bayesian information criterion (BIC) between models in which ∆z is fixed to a value of 0.035, as suggested by comparing our spectroscopic and photometric redshifts, to models where it is a free parameter. A circle indicates that the ∆z free model is favored, i.e., has the minimum BIC, whereas a square indicates that the ∆z fixed model is favored. The different colors represent different photometric redshift bins as indicated in the legend. For clarity, samples are displayed from left to right in order of increasing stellar mass threshold within each redshift bin. For reference, a dashed line is drawn at a value of ∆BIC = 6. A value of ∆BIC > 6 indicates that this model is 'strongly' favored. The right panel shows the fractional difference on the derived best-fit characteristic halo masses, M h,min (circles) and M h,1 (triangles), resulting from fixing ∆z = 0.035, compared to leaving it as a free parameter in our model fitting. For clarity, the results for each sample are displayed as described for the left panel.

distribution computed directly from the SMUVS catalog, which is based on the best-fit photometric redshift computed by LePhare, and so add confidence to the values of ∆ z determined by our method.

B. RECONSTRUCTING THE STELLAR MASS FUNCTION

In this appendix, we use our best-fit SHMRs from Section 4.5 (see Table 2) to reconstruct the stellar mass function for 1.5 < z < 2.0 and 2.0 < z < 3.0. For this, we begin with a halo mass function generated according to the prescription of Tinker et al. (2008) and incorporate a constant Gaussian scatter around our best-fit SHMR with a dispersion of 0.16 dex (More et al. 2009;Moster et al. 2010). We show the results in Figure 16, where we compare our SHMR-derived stellar mass function to those computed directly from the SMUVS catalog using the method described in Deshmukh et al. (2018). The two methods are in good agreement, highlighting the consistency of our analyses. One -Best-fit absolute redshift dispersion as a fraction of the photometric redshift bin width (here z 50 is the median stellar mass of a sample and z hi and z lo represent the limits of the photometric redshift bin). For clarity, the results for each of our samples are displayed as described for Figure 13. The black dashed lines indicate a value of ∆z = 0.035, suggested by comparing our spectroscopic and photometric redshifts. For reference, a horizontal gray dashed line is drawn at unity. can see clearly how the peaked nature of the SHMR converts the halo mass function to the Schechter function shape of the observed galaxy stellar mass function. We have ignored satellite galaxies in converting a halo mass function to a stellar mass function as shown in Figure 16. However, as our satellite fractions are generally small ( 0.1) this does not affect the conclusion of this section.

C. TABULATED RESULTS

In this appendix, we tabulate our main results from Section 4 in Table 3.

D. EXAMPLE LIKELIHOOD DISTRIBUTIONS

In this appendix, we show some examples of the likelihood distributions produced by our MCMC fitting procedure in Figure 17. These correspond to the samples shown in the upper right-hand panel of Figure 2. Note.-All masses in log 10 (M/M ). Error bars are 1σ. † Each of our fits has nine degrees of freedom. There are 12 clustering data points, plus the observed number of galaxies. The model has three free parameters. † † For this sample, we are unable to fully constrain all of the model parameters, therefore we do not quote a minimum χ 2 red .

9.5 10.0 10.5 11.0 11.5 12.0 12.5 log 10 (M /M ) -

Fig. 1 .

 1 Fig. 1.-The stellar mass -redshift plane for SMUVS galaxies. The dark blue lines indicate our stellar mass-selected volumelimited samples. The light gray line indicates the 80 percent stellar mass completeness limit. The color scale indicates the density of objects at that point on the plane (darker color indicates a higher density). 1999; Ilbert et al. 2006). We assume Bruzual & Charlot (2003) templates corresponding to a simple stellar population formed with a Chabrier (2003) stellar initial mass function and either solar or sub-solar (1 Z or 0.2 Z ) metallicity, and allow for the addition of nebular emission lines. Additionally, we assume exponentially declining star formation histories.Photometric redshifts and stellar mass estimates are obtained for > 99 percent of our sources. Using ancillary spectroscopic data in COSMOS to assess the quality of the obtained photometric redshifts, we found that the standard deviation, σ z , of |z phot -z spec |/(1+z spec ), based on ∼ 1.4 × 10 4 galaxies with reliable spectroscopic redshifts in the COSMOS field (see Table1in Ilbert et al., 2013 and references therein) is 0.026 (Deshmukh et al. 2018). This statistic is computed excluding outliers, defined as objects for which |z phot -z spec |/(1+z spec ) > 0.15, which comprise ∼ 5.5 percent of the spectroscopic catalog. These results compare favorably with other photometric surveys in the literature (e.g. Ilbert et al. 2013; Laigle et al. 2016) and highlight the high accuracy of our derived photometric redshifts.Throughout we use the best-fit redshifts and stellar masses computed by LePhare.

  in Ilbert et al., 2013 and references therein) is 0.026 (Deshmukh et al. 2018). This statistic is computed excluding outliers, defined as objects for which |z phot -z spec |/(1+z spec ) > 0.15, which comprise ∼ 5.5 percent of the spectroscopic catalog. These results compare favorably with other photometric surveys in the literature (e.g. Ilbert et al. 2013; Laigle et al. 2016) and highlight the high accuracy of our derived photometric redshifts.

  5 as the COSMOS field has a well-documented overabundance of rich structures (McCracken et al. 2007; Meneux et al. 2009; Mc-Cracken et al. 2015) at z ∼ 1 -1.5 that can nullify the halo modeling analysis performed in this work, and we wish to focus on the high-redshift nature of our catalog. Halo modeling analyses have been performed for z 1 in many previous studies (e.g. Leauthaud et al. 2011; Coupon et al. 2012; McCracken et al. 2015) which cover a larger volume at these redshifts than SMUVS.

  A Navarro-Frenk-White (NFW; Navarro et al. 1997) halo density profile, with the concentration relation of Bullock et al. (2001) is assumed. The halo mass function used is the parametrization of Tinker et al. (2008), with the high-redshift correction of Behroozi et al. (2013b), as well as the large-scale dark matter halo bias parametrization of Tinker et al. (2010). These Tinker et al. relations adopt the definition of a halo as a spherical overdensity of 200 relative to the mean cosmic density at the epoch of interest (e.g. Lacey & Cole 1994; Tinker et al. 2008).

  where b h (M h , z) is the large-scale halo bias parametrization of Tinker et al. (2010).

Fig. 2 .

 2 Fig. 2.-Stellar mass-selected angular correlation function measurements in SMUVS. Each panel represents a different redshift bin. Colors and symbols indicate different stellar mass limits as shown in the legends. Solid lines show the angular correlation function of our best-fit halo models. The dashed lines indicate the angular correlation function of dark matter; note that these have been projected incorporating the best-fit value for ∆z. The inset panels show the best-fit halo occupation distribution; dotted lines here show the contribution from satellite galaxies.

Fig. 3 .

 3 Fig. 3.-Characteristic halo masses M h,min (circles) and M h,1 (triangles) as a function of galaxy sample median stellar mass. Colors indicate different redshift bins, as shown in the upper legend.

Fig. 4 .

 4 Fig. 4.-Satellite fraction as a function of galaxy sample median stellar mass. Colors indicate different redshift bins, as shown in the legend. Observational data from McCracken et al. (2015, open squares) are also shown.

3 Fig. 5 .

 35 Fig. 5.-Comoving correlation length as a function of galaxy sample median stellar mass. Colors indicate different redshift bins, as shown in the legend. Observational data from Wake et al. (2011, open blue triangles), Marulli et al. (2013, open stars), and Durkalec et al. (2015, open red triangle) are also shown. We have approximately scaled the Wake et al. data from their quoted stellar mass threshold to a median stellar mass by applying a similar shift between these two quantities as found in the SMUVS data. This change is indicated by the horizontal error bars on the Wake et al. data.

  -Manso et al. (2015) data, with which our determinations agree excellently. Our bias values are also in reasonable agreement with McCracken et al. (2015), though they find a flatter trend at high stellar masses. This may be in part due to our inclusion of ∆ z as a free parameter, as calculations in which we fix ∆ z = 0.035 (as suggested by the σ z value for our catalog in the range 1.5 < z < 5.0) produce a flatter relationship in better agreement with the McCracken et al. data. Martinez-Manso et al. account for the redshift dispersion as they use the full redshift probability distributions for each galaxy in projecting the spatial correlation function produced by their HOD model. The bias values found by McCracken et al. from their 2.0 < z < 2.5 measurements (b gal

  Fig. 6.-Inferred large-scale galaxy bias as a function of stellar mass and redshift. Left panel: large-scale bias as a function of galaxy sample median stellar mass. Colors indicate different redshift bins, as shown in the legend. Other observational data are from Durkalec et (2015, open triangle), Martinez-Manso et al. (2015, open diamonds), and McCracken et al. (2015, open squares) are also shown. Right panel: large-scale bias as a function of redshift. The color scale indicates the galaxy sample median stellar mass. Solid gray lines indicate the bias evolution of halos of a constant mass (indicated in the panel), according to the prescription of Tinker et al. (2010).

  Fig. 7.-Same as for the right panel of Figure 6 but also indicating a potential evolutionary path for the host dark matter halos and stellar mass content of the progenitors of galaxies with log 10 (M /M ) > 10.8 at 1.5 < z < 2.0. The solid black line indicates the bias of the median halo mass accretion history, which is computed according to McBride et al. (2009), and the colored line on top of this indicates the median stellar mass of the progenitors, computed following the cumulative abundance argument of Behroozi et al. (2013a). The dashed black lines indicate the 16 -84 percentile scatter of the halo mass accretion histories.

- 2 Fig. 8 .

 28 Fig. 8.-Examples of our proposed evolutionary paths. Dark matter halo (black line) and stellar mass (red line) assembly histories computed as described in Section 4.4.1, as a function of lookback time. Halo masses are converted from bias measurements presented in Figure 7 using the relation of Tinker et al. (2010). The gray shaded region describes the 16 -84 percentile scatter of the halo mass accretion histories. The median of dark matter accretion histories has been scaled by 10 -2 for presentation purposes. Data from Figure 7 that intersect this history are also shown as points with error bars. The horizontal error bars indicate the range of the redshift bin that point is drawn from, in terms of lookback time. The dashed lines indicate the 1σ scatter on the stellar mass history. The starting point for the histories calculation is galaxies with log 10 (M /M ) > 10.8 in the left panel and > 10.2 in the right panel.

  Fig. 9.-Ratio between the median stellar mass and M h,min for each galaxy sample for the redshifts indicated in the panel. The solid lines in the 1.5 < z < 2.0 and 2.0 < z < 3.0 panels indicate the best-fit SHMRs of Behroozi et al. (2013b), with the 1σ (16 -84 percentile) uncertainties indicated by the gray shaded region.

13 Fig

 13 Fig. 10.-Location of the maximum in the SHMR (M peak ) as a function of redshift. Observational data from Zehavi et al. (2011, open star), Leauthaud et al. (2012, red downward triangles), Coupon et al. (2012, green triangles), Martinez-Manso et al. (2015, yellow diamond), and McCracken et al. (2015, light blue squares) are also shown. The solid line indicates a best-fit powerlaw relation in (1 + z) to the observational data shown. Simulation data from the abundance-matching model of Behroozi et al. (2013b, dashed black line) are also shown, with the gray shaded region indicating the 1σ uncertainty of this quantity.

  10 (M h,min /h -1 M )

Fig. 11 .

 11 Fig. 11.-Example one-dimensional (diagonal panels) and two-dimensional (off-diagonal panels) likelihood distributions produced from our MCMC fitting procedure for log 10 (M /M ) > 10.4 galaxies in the range 2.0 < z phot < 3.0, based on the clustering of a 2 deg 2 simulated galaxy catalog produced by the Lacey et al. (2016) GALFORM model. The different colors represent different values of input photometric redshift dispersion, ∆z, as indicated in the legend. The solid lines in the bottom right panel indicate the best-fit value of ∆z, which is in remarkable agreement with the input value. The contours in the off-diagonal panels indicate the 1, 2 and 3 σ regions.

1 Fig. 13

 113 Fig. 12.-Left panel: the angular correlation function measured from a 2 deg 2 simulated galaxy catalog produced by the Lacey et al. (2016) GALFORM model (points with error bars) for simulated galaxies with 2.0 < z phot < 3.0 and log 10 (M /M ) > 10.4. The solid lines indicate the best-fit halo model, according to our modeling procedure outlined in Section 3. The resulting HODs are shown in the inset panel; the HOD predicted by the GALFORM model is also shown here for reference (gray line). The different colors represent different values of input photometric redshift dispersion, as indicated in the legend. Right panel: shows the photometric and true redshift distributions (blue and red histograms respectively) for galaxies selected with 2.0 < z phot < 3.0 and log 10 (M /M ) > 10.4 with ∆z = 0.08. The best-fit HOD reconstruction of the true redshift distribution is shown as the dashed red line.

  Fig.14.-Best-fit absolute redshift dispersion as a fraction of the photometric redshift bin width (here z 50 is the median stellar mass of a sample and z hi and z lo represent the limits of the photometric redshift bin). For clarity, the results for each of our samples are displayed as described for Figure13. The black dashed lines indicate a value of ∆z = 0.035, suggested by comparing our spectroscopic and photometric redshifts. For reference, a horizontal gray dashed line is drawn at unity.

Fig. 15 .

 15 Fig. 15.-Example photometric redshift distribution of SMUVS galaxies (blue lines) and the best-fit HOD reconstruction (red lines) in the range 1.5 < z < 3.0 (indicated by the vertical dashed lines). Stellar mass thresholds used are log 10 (M /M ) > 9.4 (left panel) and log 10 (M /M ) > 10.2 (right panel).

1 dn

 1 Fig. 16.-The SMUVS stellar mass functions for 1.5 < < 2.0 (left panel) and 2.0 < z < 3.0 (right panel). The solid lines indicate this statistic computed using the best-fit SHMRs from Section 4.5 (see Table 2), with the propagated 16 -84 percentile errors indicated by the gray shaded region. For reference, the dashed line indicates the halo mass function of Tinker et al. (2008). The points with error bars indicate the stellar mass function computed directly from the SMUVS catalog as described by Deshmukh et al. (2018).

  10 (M h,min /h -1 M ) log 10 M > 9.4 log 10 M > 10.0 log 10 M > 10.6

Fig. 17 .

 17 Fig. 17.-Example one-dimensional (diagonal panels) and two-dimensional (off-diagonal panels) likelihood distributions produced from our MCMC fitting procedure for stellar mass-selected samples in the 2.0 < z < 3.0 redshift range. The contours in the off-diagonal panels represent the 1, 2, and 3 σ regions. The solid and dashed lines in the diagonal panels represent the best-fit values and 1σ errors, respectively. The different colors correspond to a different stellar mass threshold as shown in the legend. The samples correspond to the clustering and best-fit halo models shown in the upper right-hand panel of Figure 2.

  al. 2014; Schaye et al. 2015) or Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark 5 Institut d'Astrophysique de Paris, CNS & UPMC, UMR 7095, 98 bis Boulevard Arago, F-75014, Paris, France 6 Department of Physics and Astronomy, Rutgers University, The State University of New Jersey, 136 Frelinghuysen Road, Piscataway, NJ 08854 USA semi-analytical (e.g., Somerville et al. 2012; Henriques et al. 2015;

	1 Kapteyn Astronomical Institute, University of Groningen,
	P.O. Box 800, 9700 AV, Groningen, The Netherlands; cow-
	ley@astro.rug.nl 2 Harvard-Smithsonian Center for Astrophysics, 60 Garden
	Street, Cambridge, MA 02138, USA 3 Aix Marseille Université, CNRS, LAM (Laboratoire
	d'Astrophysique de Marseille) UMR 7326, F-13388, Marseille,
	France 4

  ; Brainerd & Specian 2003; van den Bosch et al. 2004; Norberg et al. 2008;

2 .

 2 THE SMUVS SURVEY 2.1. Survey overview The SMUVS program (PI: K. Caputi; Ashby et al. 2018) has collected ultra-deep Spitzer 3.6 and 4.5 µm data over the region of the COSMOS 7 field (Scoville et al. 2007) overlapping with three of the UltraVISTA ultradeep stripes (McCracken et al. 2012) with deep optical coverage from the Subaru Telescope

TABLE 1

 1

			Characteristics of each galaxy sample
		1.5 < z < 2.0	2.0 < z < 3.0	3.0 < z < 4.0	4.0 < z < 5.0
	Threshold (a)	N gal	M	(a) ,50	N gal	M	(a) ,50	N gal	M	(a) ,50	N gal
	9.00	19 637	9.49	-		-	-		-	-	-
	9.20	15 309	9.64	-		-	-		-	-	-
	9.40	11 232	9.86	18 362	9.73	-		-	-	-
	9.60	8 186	10.09	12 201	9.92	6 187	9.89	-	-
	9.75	-		-	-		-	-		-	1 768	10.02
	9.80	6 129	10.27	7 851	10.14	3 894	10.04	-	-
	10.00	4 634	10.42	5 194	10.35	-		-	-	-
	10.20	3 458	10.54	3 516	10.52	-		-	-	-
	10.40	2 398	10.66	2 357	10.66	-		-	-	-
	10.60	1 469	10.79	1 446	10.81	-		-	-	-
	10.80	702	10.95	740	10.93	-		-	-	-

  We also choose to fix σ log M h = 0.2 and α sat = 1.0 as these are standard values adopted throughout the literature and are supported by both theoretical and observational studies (see e.g.Wake et al. 2011; Martinez-Manso et al. 2015; Harikane et al. 2016, and references therein).

	, who propose,
	based on their simulations, the relationship
	log 10 M h,0 = 0.76 log 10 M h,1 + 2.3.	(14)

11

http://dan.iel.fm/emcee/

TABLE 2

 2

			SHMR best-fit parameter values
	redshift bin	log 10 (ε) log 10 (Mt/M )	α	δ	γ	log 10 (M peak /M )
	1.5 < z < 2.0 -1.81 +0.10 -0.20 2.0 < z < 3.0 -1.78 +0.06 -0.11	11.92 +0.20 -0.24 12.17 +0.19 -0.22	-1.64 +0.39 -0.31 3.30 +1.56 -1.50 0.59 +0.27 -0.28 -1.31 +0.48 -0.33 2.81 +1.51 -1.43 0.59 +0.28 -0.35	12.33 +0.07 -0.06 12.50 +0.10 -0.08

  This work is based in part on observations carried out with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. It also uses data products from observations conducted with ESO Telescopes at the Paranal Observatory under ESO program ID 179.A-2005, data products produced by TER-APIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium, observations carried out by NASA/ESA Hubble Space Telescope obtained and archived at the Space Telescope Science Institute; and the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. This re-search has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. The authors would like to thank Cedric Lacey and the GALFORM collaboration for use of the GALFORM model in Appendix A. This work has made use of the open-source Python packages: Numpy, Scipy, Matplotlib and IPython. The color-vision-impairedfriendly color schemes used throughout can be found at https://personal.sron.nl/ ~pault/.

	W.I.C., S.D., and K.I.C. acknowledge funding from the
	European Research Council through the award of the
	Consolidator Grant ID 681627-BUILDUP.
	APPENDIX
	A. MODELING THE PHOTOMETRIC REDSHIFT DISPERSION

TABLE 3 Tabulated

 3

	main results

This is in fact constrained to be the case by the assumption that our central galaxy HOD reaches an amplitude of unity, the fixed value of σ log M h and that the observed number of galaxies is used as a constraint in our fitting procedure.

In practice, we used the Behroozi et al. code publicly available at https://code.google.com/archive/p/nd-redshift/ to compute the progenitor number density and its 1σ scatter.

The Millennium simulations assumed Ωm = 0.25, Ω b = 0.045, Ω Λ = 0.75, h = 0.73, σ 8 = 0.9, and ns = 1.

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous referee for a comprehensive and constructive report that allowed us to improve the overall quality of the manuscript.


David Moreira

Purificacion Lopez
email: puri.lopez@u-psud.fr

Purificación López-García

Time series are critical to understand microbial plankton diversity and ecology

published or not. The documents may come

'As there are good and bad wine and fruit years, so there are rich and barren plankton years' [START_REF] Haeckel | Plankton-Studien[END_REF]. With these eloquent words, Ernst Haeckel exemplified more than one century ago the critical role that time series must play in the study of plankton, which cannot be fully understood just on the basis of punctual observations, as they may lead to wrong conclusions if done in moments when, for whichever reason, plankton does not behave normally. Indeed, 'to obtain a complete and more certain survey of the temporary variations of plankton composition requires an unbroken series of observations, carried on at one and the same place at least for the space of a full year-still better for several successive years-to obtain from the yearly and monthly oscillations a general average'. This idea is behind the remarkable proliferation of marine stations across the world during the XIX th and early XX th centuries. For decades, naturalists were observing marine plankton, describing new species, their interactions, and the functioning of the complex plankton ecosystem. This continuous surveillance also allowed discovering periodical rhythms in the abundance of particular species, leading to the elaboration of the concept of "periodic plankton" [START_REF] Hensen | Einige Ergebnisse der Plankton-Expedition der Humboldt-Stiftung[END_REF]. Because of the obvious easiness of observation and identification, early work largely concentrated on planktonic animals. By contrast, the study of microbial planktonic populations was more challenging but, after some decades of scrutiny, periodicity begun to be described, especially for some conspicuous species (see Fig. 1 for some examples), with strong emphasis on phytoplanktonic species because of their importance as primary producers in the global C cycle (e.g., diatoms [START_REF] Pratt | The phytoplankton of Narragansett Bay[END_REF]. Knowledge remained for a long time much more partial for the difficult-toobserve, smallest planktonic size fractions (pico-and nano-plankton, 0.2-3 μm and 3-20 μm, respectively), mostly composed of prokaryotes but also of a variety of tiny eukaryotic cells. Among those eukaryotes, some important algal species, such as the globally distributed Micromonas pusilla [START_REF] Butcher | Contribution to our knowledge of the smallest marine algae[END_REF], were described, though most diversity was considered to correspond to heterotrophic, likely bacterivorous, species. These organisms were recognized to play a significant role in the function of the plankton trophic chain by fuelling the "microbial loop" [START_REF] Azam | The ecological role of water-column microbes in the sea[END_REF]. However, the precise taxonomic identity of these tiny eukaryotic predators remained for the most part unknown and were collectively classified as "eukaryotic nanoflagellates" and treated as a sort of black box in ecological studies [START_REF] Fenchel | Ecology of heterotrophic microflagellates. II. Bioenergetics and growth[END_REF]. Only in the early 2000's, thanks to 18S rRNA gene PCR amplification and sequencing from environmental DNA, it started to become possible to fully appreciate the actual wealth of diversity that this black box hid [START_REF] Díez | Study of genetic diversity of eukaryotic picoplankton in different oceanic regions by small-subunit rRNA gene cloning and sequencing[END_REF][START_REF] López-García | Unexpected diversity of small eukaryotes in deep-sea Antarctic plankton[END_REF][START_REF] Moon-Van Der Staay | Oceanic 18S rDNA sequences from picoplankton reveal unsuspected eukaryotic diversity[END_REF]. Metabarcoding, initially based on time-consuming Sanger sequencing and now on high-throughput techniques, revealed an extensive diversity of small-sized eukaryotes in marine plankton, largely dominated by heterotrophic -including parasitic-species [START_REF] De Vargas | Ocean plankton. Eukaryotic plankton diversity in the sunlit ocean[END_REF][START_REF] Moreira | Molecular ecology of microbial eukaryotes unveils a hidden world[END_REF]. Not surprisingly, the relatively simple morphologies of these small eukaryotes concealed a vast cryptic phylogenetic diversity that has been intensively studied in recent years though several large-scale plankton sampling and metabarcoding efforts (e.g., [START_REF] De Vargas | Ocean plankton. Eukaryotic plankton diversity in the sunlit ocean[END_REF][START_REF] Pernice | Large variability of bathypelagic microbial eukaryotic communities across the world's oceans[END_REF]. However, these efforts, despite their wide geographical scope, lack the temporal dimension indispensable to really understand plankton community composition, dynamics, biogeography, and ecology, something that only the examination of time series can achieve.

Analysis of time series of small-sized plankton samples using molecular techniques, such as terminal restriction fragment length polymorphism and metabarcoding, had been applied for several years to study prokaryotic communities and, more rarely, eukaryotic ones. These studies have revealed seasonal rhythmicity of the marine and freshwater eukaryotic communities, notably by measuring pairwise Bray-Curtis similarity values between communities sampled at different time intervals (e.g., [START_REF] Kim | Monthly to interannual variability of microbial eukaryote assemblages at four depths in the eastern North Pacific[END_REF][START_REF] Simon | Marked seasonality and high spatial variability of protist communities in shallow freshwater systems[END_REF]. However, this averaging approach misses the fine-grain information contained in the temporal dynamics of each individual species or OTU. Analyzing that information is challenging because the structure of the planktonic eukaryotic communities is very complex, as they contain hundreds or thousands of different OTUs and most of these OTUs appear at first glimpse to exhibit a rather sporadic temporal behavior. [START_REF] Giner | Quantifying long-term recurrence in planktonic microbial eukaryotes[END_REF] have developed a new Recurrence Index (RI) to specifically address this question on a rich dataset of plankton samples collected monthly over 10 years at the Blanes Bay Microbial Observatory (Mediterranean Sea) and submitted to 18S rRNA gene metabarcoding analysis.

For both eukaryotes and prokaryotes, a recurrent observation in previous marine plankton metabarcoding analyses is that diversity patterns are most often dominated by rare OTUs occurring at very low frequencies and in single or very few samples [START_REF] Pedrós-Alió | Marine microbial diversity: can it be determined?[END_REF][START_REF] Sogin | Microbial diversity in the deep sea and the underexplored "rare biosphere[END_REF]. By applying their new RI index, Giner et al. corroborated this expectation as they found that only 13% and 19% of the pico-and nano-eukaryotic OTUs, respectively, showed seasonality.

However, these OTUs represented a large proportion of the sequence reads (39% and 37%, respectively), indicating that they are quantitatively important members of the small-size eukaryotic plankton community. The high abundance of these OTUs probably drives the seasonality pattern at the global eukaryotic community level inferred from the Bray-Curtis similarity analyses.

Among the OTUs with the highest RI values, there are several Mamiellophyceae green algae. Planktonic microalgae have traditionally been considered seasonal organisms that tend to bloom during spring and summer as a response to favorable conditions (e.g., increased daylight availability and warmer waters). Some heterotrophic eukaryotes also appear among the high-RI OTUs, in particular MALV-III marine alveolates. For these organisms, seasonality appears to be a conserved class-level trait since the most abundant OTUs of these groups have high RI. By contrast, non-seasonality can be the general trait of most OTUs in other taxonomic groups, both photosynthetic (e.g., Chlorodendrophyceae) and heterotrophic (e.g., cercozoa) ones. Similarly, the proportion of seasonal and non-seasonal OTUs varies in other groups. Thus, seasonality seems to be a trait that differs among classes and even among OTUs within classes. This result stresses that eukaryotic taxonomic groups must not be regarded as homogeneous assemblages but, in most cases, as bins of species with idiosyncratic temporal dynamics, most likely owing to their specific niche adaptations.

A general assumption to explain plankton seasonality is that algal growth, driven by seasonal changes in physicochemical parameters (light, temperature, nutrient availability) drives in turn the abundance of heterotrophic species (e.g., predators and parasites) that depend on them for growth. However, seasonality has also been observed for deep-sea protists, probably much less influenced by those changes in the phytoplanktonic species [START_REF] Dolan | An exploratory study of heterotrophic protists of the mesopelagic Mediterranean Sea[END_REF]. Likewise, for most eukaryotic heterotrophs, it remains to be elucidated whether they rely directly on the co-occurring algae or whether they rely on other types of prey, notably bacteria. This open question raises an important point concerning time series: sampling frequency. Aquatic microbial communities are highly dynamic and many microorganisms can exhibit fast population size fluctuations [START_REF] Mangot | Short-term dynamics of diversity patterns: evidence of continual reassembly within lacustrine small eukaryotes[END_REF]. Thus, given the complexity of plankton communities both in terms of structure and temporal oscillation of the composing species, it seems clear that even a monthly sampling strategy most likely misses part of such rapid plankton dynamics and still does not allow to fully infer the fine temporal interplay between autotrophic and heterotrophic plankton populations. The next challenge in this field will be to apply similar techniques as

  How diverse are marine planktonic protist communities? How much seasonality do they exhibit? For very long time, these two old and challenging questions in the field of plankton ecology could be addressed only for large-size protist species, based on cell counting under the microscope. The recent application of molecular techniques, notably massive marker-gene amplicon sequencing or metabarcoding approaches, has allowed investigating with unprecedented level of resolution the smallsized (<20 µm) planktonic eukaryotes too. An amazing diversity of these tiny organisms has been unveiled but details about their temporal dynamics have remained much more elusive. In this issue of Molecular Ecology, Giner et al. (2018) introduce a new Recurrence Index (RI) to specifically look for seasonality in time-series metabarcoding data. They inspected the temporal dynamics of all operational taxonomic units (OTUs) in a rich sequence dataset of pico-and nano-planktonic eukaryotes in samples collected monthly during 10 years. Although most OTUs did not show seasonality, some abundant ones did, which explains why some averaging methods can find seasonality at the less-detailed level of whole planktonic communities. Not surprisingly, the very complex small-sized eukaryotic plankton communities are composed of organisms with miscellaneous temporal dynamics.

Figure legend FIGURE 1

 legend1 Figure legend FIGURE 1 Large protists with seasonal dynamics: the tintinnid ciliate Stenosomella nivalis (a) and the dinoflagellate Ceratium furca (b) have abundance peaks during winter, whereas their counterparts Protorhabdonella curta (c) and Ceratium pentagonum (d) increase their abundance during summer (photographs by John Dolan).

Grant/Award Number: ANR-15-CE32-0003