# nIFTy Cosmology: the clustering consistency of galaxy formation models

Arnau Pujol, Ramin A. Skibba, Enrique Gaztañaga, Andrew Benson, Jeremy Blaizot, Richard Bower, Jorge Carretero, Francisco J. Castander, Andrea Cattaneo, Sofia A. Cora, et al.

► To cite this version:

Arnau Pujol, Ramin A. Skibba, Enrique Gaztañaga, Andrew Benson, Jeremy Blaizot, et al. nIFTy Cosmology: the clustering consistency of galaxy formation models. Monthly Notices of the Royal Astronomical Society, 2017, 469 (1), pp.749-762. 10.1093/mnras/stx913. hal-01645208

HAL Id: hal-01645208

https://hal.science/hal-01645208v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# nIFTy cosmology: the clustering consistency of galaxy formation models

Arnau Pujol,¹,²* Ramin A. Skibba,³,⁴ Enrique Gaztañaga,¹ Andrew Benson,⁵

Jeremy Blaizot,⁶,⁷,⁸ Richard Bower,⁹ Jorge Carretero,¹⁰ Francisco J. Castander,¹

Andrea Cattaneo,¹¹ Sofia A. Cora,¹²,¹³ Darren J. Croton,¹⁴ Weiguang Cui,¹⁵

Daniel Cunnama,¹⁶,¹⁷ Gabriella De Lucia,¹⁸ Julien E. Devriendt,¹⁹ Pascal J. Elahi,²⁰

Andreea Font,²¹ Fabio Fontanot,¹⁸ Juan Garcia-Bellido,¹⁵,²² Ignacio D. Gargiulo,¹²,¹³

Violeta Gonzalez-Perez,⁹,²³ John Helly,⁹ Bruno M. B. Henriques,²⁴,²⁵

Michaela Hirschmann,²⁶ Alexander Knebe,²²,²⁷ Jaehyun Lee,²⁸ Gary A. Mamon,²⁶

Pierluigi Monaco,¹⁸,²⁹ Julian Onions,³⁰ Nelson D. Padilla,¹² Frazer R. Pearce,³⁰

Chris Power,³¹ Rachel S. Somerville,³² Chaichalit Srisawat,³³ Peter A. Thomas,³³

Edouard Tollet,¹¹ Cristian A. Vega-Martínez¹² and Sukyoung K. Yi²⁸

*Affiliations are listed at the end of the paper*

Accepted 2017 April 11. Received 2017 April 7; in original form 2017 February 7

## ABSTRACT

We present a clustering comparison of 12 galaxy formation models [including semi-analytic models (SAMs) and halo occupation distribution (HOD) models] all run on halo catalogues and merger trees extracted from a single Λ cold dark matter N-body simulation. We compare the results of the measurements of the mean halo occupation numbers, the radial distribution of galaxies in haloes and the two-point correlation functions (2PCF). We also study the implications of the different treatments of orphan (galaxies not assigned to any dark matter subhalo) and non-orphan galaxies in these measurements. Our main result is that the galaxy formation models generally agree in their clustering predictions but they disagree significantly between HOD and SAMs for the orphan satellites. Although there is a very good agreement between the models on the 2PCF of central galaxies, the scatter between the models when orphan satellites are included can be larger than a factor of 2 for scales smaller than $1h^{-1}$ Mpc. We also show that galaxy formation models that do not include orphan satellite galaxies have a significantly lower 2PCF on small scales, consistent with previous studies. Finally, we show that the 2PCF of orphan satellites is remarkably different between SAMs and HOD models. Orphan satellites in SAMs present a higher clustering than in HOD models because they tend to occupy more massive haloes. We conclude that orphan satellites have an important role on galaxy clustering and they are the main cause of the differences in the clustering between HOD models and SAMs.

**Key words:** methods: numerical – galaxies: haloes – cosmology: theory.

## 1 INTRODUCTION

In Λ cold dark matter (ΛCDM) cosmology, gravitational evolution causes dark matter to cluster around peaks of the initial density field and to collapse into virialized objects (i.e. dark matter haloes). Structures form hierarchically, such that smaller haloes merge to form larger and more massive haloes. All galaxies are thought to

form as a result of gas cooling at the centre of the potential well of dark matter haloes. When a halo and its 'central' galaxy are accreted by a larger halo, it becomes a subhalo and its galaxy becomes a 'satellite' galaxy. However, due to tidal stripping and the gravitational interaction of subhaloes with their environment (other subhaloes, the gravitational potential of the halo centre, etc.), sometimes can be disrupted and the galaxy, if it survives, becomes an 'orphan' galaxy. In addition to mergers, haloes also grow by smooth accretion and galaxies grow by in situ star formation when fuel (i.e. cold gas) is available (Cooray & Sheth 2002; Sheth &

<footnote>* E-mail: arnaupv@gmail.com</footnote>
Tormen 2002; van den Bosch 2002; De Lucia et al. 2004; Gill
et al. 2004; van den Bosch, Tormen & Giocoli 2005; Diemand,
Kuhlen & Madau 2007; Giocoli, Tormen & van den Bosch 2008).

In this paradigm of hierarchical structure formation, there is a
correlation between halo formation, their abundances and the sur-
rounding large-scale structure where more massive haloes tend to
reside (Mo & White 1996; Sheth & Tormen 2002). Most galaxy for-
mation models implicitly assume that the properties of a galaxy are
determined primarily by the mass and formation history of the dark
matter halo within which it formed (White & Rees 1978; Cole 1991;
Lacey & Silk 1991; White & Frenk 1991; Baugh et al. 1999; Benson
et al. 2001). Thus, the correlation between halo properties and en-
vironment (matter density, substructure, etc.) induces a correlation
between galaxy properties and environment.

There are multiple statistical quantities used to study large-scale
structure, and here we focus on the two-point correlation func-
tion (2PCF), the radial distribution and the mean occupation num-
ber of galaxies. Clustering studies have shown that a variety of
galaxy properties (such as luminosity, colour, stellar mass, star for-
mation rate and morphology) are dependent on the environment
and halo properties across a wide range of scales. Galaxy forma-
tion models in simulations are crucial to study the connection be-
tween galaxies and haloes, and hence it is important to understand
the consistency or differences between different galaxy formation
models.

Galaxy formation is a complex, non-linear process, driven
by the interplay of many different physical mechanisms (e.g.
Benson 2010). The goal of galaxy formation models is to estimate
the statistical properties of the galaxy population given some set of
assumptions and thereby to better understand the physical processes
involved. One fruitful approach has been to utilize semi-analytic
models (SAMs) of galaxy formation (e.g. Cole et al. 2000; Hatton
et al. 2003; Baugh 2006, 2013; Cattaneo et al. 2006; Cora 2006;
Croton et al. 2006; De Lucia & Blaizot 2007; Monaco, Fontanot
& Taffoni 2007; Lo Faro et al. 2009; Benson 2012; Henriques
et al. 2013; Lee & Yi 2013; Gonzalez-Perez et al. 2014; Gargiulo
et al. 2015), in which a statistical estimate of the distribution of
dark matter haloes and their merger history – either coming from
cosmological simulations or extended Press–Schechter/Lagrangian
methods – is combined with simplified yet physically motivated
prescriptions of processes such as star formation, gas cooling, feed-
back from supernovae and active galactic nuclei (AGN), etc. that
allows one estimate the distribution of galaxy properties. New mod-
els are starting now to also use observations of galaxy clustering to
constrain their parameters (van Daalen et al. 2016).

An alternative approach to SAMs is (analytic) dark matter halo
occupation models, which determine the halo occupation of galax-
ies based on the properties of their parent halo. Usually observations
of clustering are used to constrain this occupation. This approach is
used to study the link between galaxy formation and halo assembly
(see Cooray & Sheth 2002; Mo, van den Bosch & White 2010 for
a review). Halo models of galaxy abundances and clustering gen-
erally consist of halo occupation distribution (HOD, or conditional
luminosity functions; e.g. Seljak 2000; Scoccimarro et al. 2001;
Berlind & Weinberg 2002; Cooray & Sheth 2002; Yang, Mo &
van den Bosch 2003; Kravtsov et al. 2004; Cooray 2006; Guo
et al. 2016) and (sub)halo abundance matching (HAM or SHAMS;
Conroy, Wechsler & Kravtsov 2006; Vale & Ostriker 2006; Hearin
et al. 2013; Reddick et al. 2013; Guo et al. 2016). Such models
are useful for exploring the relations between galaxy formation and
dark matter halo assembly in the context of the large-scale structure
of the Universe.

Subhaloes closer to the halo centre tend to accrete earlier on
(Gao et al. 2004), and therefore tidal stripping (a process that is
in part numerical but also physical) has more time to act on these
subhaloes. Hence, they are more frequently disrupted and this is the
reason why subhaloes are antibiased with respect to the dark matter
(see e.g. Ghigna et al. 2000; Diemand, Moore & Stadel 2004; Pujol
et al. 2014) in these regions. In simulations, lack of mass resolu-
tion causes the disappearance of a subhalo, causing the galaxy to
become an orphan. Sometimes the halo finder will merge a subhalo
with its parent halo, but the subhalo can reappear when its member
particles bounce out of the halo. There are different ways to follow
the positions of these orphan galaxies, and these can lead to different
clustering of galaxies, especially on small scales (Gao et al. 2004;
Wang et al. 2006; Guo et al. 2011; Budzynski et al. 2012; Lee
et al. 2014). In Gao et al. (2004), they used high-resolution resimu-
lations of galaxy clusters and analysed the radial density profiles of
both subhaloes and galaxies from a SAM. This study showed that by
including orphan galaxies the radial density distribution was very
close to that of the dark matter, as inferred in the observational data.
They also argued that increasing the resolution would not improve
the situation. This result has been confirmed in Guo et al. (2011),
who showed that orphans are still dominating the central regions of
galaxy clusters when increasing the resolution of the simulations.
They showed that by tracking the position of the most bound par-
ticle at the time of disruption convergence between simulations of
different resolution was achieved. In Wang et al. (2006), an HOD
approach was used, but using the number and positions of galaxies
from a SAM, and they showed that orphan galaxies are needed to
reproduce the clustering signal at small scales, also confirmed by
recent studies (Budzynski et al. 2012). Finally, Kang et al. (2012),
Guo et al. (2013) and Henriques et al. (2013) showed that cosmol-
ogy, within the current precision, has no impact in the clustering
when compared to the differences from galaxy formation physics,
even on large scales.

The trajectory and lifetime of orphan galaxies can be determined
from different approaches in SAMs. On one side, some SAMs im-
mediately merge galaxies with the central galaxy when the subhalo
is lost, and then they have no orphan satellites by construction. Other
SAMs define an analytical orbit for the orphan galaxies according to
the position and velocity of the galaxies when they became orphan.
The radius of the orbit is then continuously decreased until it merges
with the central galaxy. Finally, other SAMs define the position and
velocity of orphan galaxies directly from the dark matter particle
that was the most gravitationally bound from the disrupted subhalo.
And other SAMs (Guo et al. 2011) use a combination of both the
analytical orbits with the dark matter particle trajectories. All these
different treatments of orphan galaxies can have consequences on
the abundance and distributions of such galaxies, especially at small
scales.

SAM and HOD have important differences on the treatment of
orphan galaxies. First of all, while SAMs make use of the merger
trees to derive the initial trajectories of orphan galaxies, the HOD
models define the galaxy distribution from the present distribution of
haloes, without using information from their evolution. Moreover,
classical HODs do not account for the presence of substructures –
they are built on top of dark matter haloes and the population of
satellites is just distributed according to a Navarro–Frenk–White
(NFW) model.

Given the variety of galaxy formation models that are used
in simulations nowadays, it is important to study the differences
that arise from the different treatment of galaxy formation physics
in each model. In fact, many efforts have been done comparing
<table><caption>**Table 1.** The list of galaxy formation models used in this paper. The first column shows the acronyms used for each of the models. The second column specifies whether the model is a SAM or HOD. The third column specifies if the model has implemented a treatment of orphan satellites. The fourth column specifies if the model has calculated the positions or orbits of the orphan satellites for this work. Finally, a list of the references is shown in the last column.</caption><thead><tr><th>Model</th><th>Type</th><th>Orphans</th><th>Orphan positions</th><th>Reference</th></tr></thead><tbody><tr><td>DLB07</td><td>SAM</td><td>Yes</td><td>No</td><td>De Lucia & Blaizot (2007)</td></tr><tr><td>GALACTICUS</td><td>SAM</td><td>Yes</td><td>No</td><td>Benson (2012)</td></tr><tr><td>GALFORM</td><td>SAM</td><td>Yes</td><td>No</td><td>Gonzalez-Perez et al. (2014)</td></tr><tr><td>GALICS 2.0</td><td>SAM</td><td>No</td><td>–</td><td>Hatton et al. (2003), Cattaneo et al. (2006, 2017)</td></tr><tr><td>LGALAXIES</td><td>SAM</td><td>Yes</td><td>Yes</td><td>Henriques et al. (2013)</td></tr><tr><td>MICE</td><td>HOD</td><td>Yes</td><td>Yes</td><td>Carretero et al. (2015)</td></tr><tr><td>MORGANA</td><td>SAM</td><td>Yes</td><td>No</td><td>Monaco et al. (2007), Lo Faro et al. (2009)</td></tr><tr><td>SAG</td><td>SAM</td><td>Yes</td><td>Yes</td><td>Cora (2006), Gargiulo et al. (2015)</td></tr><tr><td>SAGE</td><td>SAM</td><td>No</td><td>–</td><td>Croton et al. (2016)</td></tr><tr><td>SKIBBASUBS</td><td>HOD</td><td>Yes</td><td>Yes</td><td>Skibba et al. (2006), Skibba & Sheth (2009)</td></tr><tr><td>SKIBBAHOD</td><td>HOD</td><td>Yes</td><td>Yes</td><td>Skibba et al. (2006), Skibba & Sheth (2009)</td></tr><tr><td>YSAM</td><td>SAM</td><td>Yes</td><td>Yes</td><td>Lee & Yi (2013)</td></tr></tbody></table>

different galaxy formation models and their physical prescriptions (Somerville & Primack 1999; Fontanot et al. 2009, 2011, 2012; Kimm et al. 2009; De Lucia et al. 2011; Contreras et al. 2013; Kang 2014; Somerville & Davé 2015; Guo et al. 2016). This study focuses on the differences in the galaxy clustering for a large variety of models run on the same simulation and with the same merger tree.

In an attempt to put together a large representation of the models from the literature in an extensive comparison study, Knebe et al. (2015, hereafter K15) presented 14 models (12 SAMs and 2 HOD models) using the same simulation input (halo catalogues and merger trees) and analysed the consistency between the models looking at the stellar mass function, the star formation, stellar-to-halo mass relations, stellar mass fractions or abundance of galaxies per halo. The present paper is a complement of K15, where we study the consistency between several galaxy formation models on the clustering and the distribution of galaxies in haloes. We analyse the consequences of the differences between the models on the distribution of the galaxies inside and outside haloes. We do this by comparing the 2PCF of galaxies, the halo occupation number and the radial distribution of galaxies in haloes. We also analyse orphan satellites separately in order to focus on the consequences from the different treatments of orphan satellites between HOD models and SAMs.

This paper is organized as follows. In the next section, Section 2, we describe the dark matter halo simulation and the orphan treatments of the galaxy formation models. Then we describe our methodology in Section 3. We present our results, including comparisons of halo occupation numbers, radial distribution of galaxies in haloes and galaxy clustering in Section 4. Finally, we end by summarizing and discussing our results in Section 5.

## 2 SIMULATION DATA

For this study we use a dark matter halo catalogue generated from a GADGET-3 N-body simulation (Springel 2005) of a $62.5 h^{-1}$ Mpc side box. We use $270^3$ particles with a particle mass resolution of $9.31 \times 10^8 h^{-1}$ M$_\odot$, producing an output of 62 snapshots. From each snapshot we generate a halo catalogue using the SUBFIND (Springel et al. 2001) code that generates haloes and subhaloes

from dark matter overdensities. We used the code MERGERTREE to generate the merger trees of the haloes.<fnref n="1" />

From the simulation we obtained several mass definitions for the haloes that can be used for the galaxy formation models. The mass definitions used are detailed in Appendix A. Some properties of the galaxy formation models can be sensitive to the mass definition and to the galaxy formation models. We discuss the mass definition criteria in Section 3.

We use several galaxy formation models together with this dark matter only simulation for the comparison analysis. Some of them are SAMs of galaxy formation, while others are based on the HOD model. We refer to K15 for a detailed description of these models and some comparisons between them. In this section we enumerate the models, their acronyms and references in Table 1, and briefly describe the treatment and merging of the orphan satellites of each, since this is one of the most relevant aspects for this comparison analysis.

All the models were originally calibrated to reproduce a given set of observations. However, each model uses different observational data and simulated cosmologies to calibrate its parameters, as stated in the corresponding papers describing the models. It is worth stressing that our strategy forces all models to the same underlying merger tree, therefore we do not expect the original calibrations to be optimal. We have seen in K15 that this leads to model-to-model variations larger than if they were all calibrated for this particular simulation. Nonetheless, we are interested in the general agreement between the different galaxy formation models.

### 2.1 Treatment of orphan satellites

The treatment of the orphan satellites (satellites with no associated dark matter subhalo) has a direct impact on galaxy clustering. In this section we give a brief overview on how models deal with the orphan population (if any) and we refer the interested reader to K15 for more details on the modelling of other physical processes.

#### 2.1.1 DLB07, GALACTICUS and GALFORM (SAM)

In these models, when a subhalo disappears (it is stripped below the resolution of the parent simulation), a merger time is assigned

<footnote n="1">¹ MERGERTREE forms part of the AHF package (Knollmann & Knebe 2009).</footnote>
to its galaxies according to some variations of the Chandrasekhar formula and galaxies are merged when this time is over. These galaxies are assumed to continue orbiting within their parent halo until dynamical friction causes it to merge with the central galaxy. Positions and velocities of orphan satellites are assumed to be traced by those of the most bound particles of substructures at the last time they were identified. This information was not provided for the simulation used in this study. Therefore, the positions of orphan satellites in these models cannot be used for the clustering analysis presented here.

## 2.1.2 GALICS 2.0 (SAM)

In this model, the effects of the merging time-scale are degenerate with those of supernova and AGN feedback and the shock-heating scale. Then, the same effects in the stellar mass function due to the contribution of orphan satellites can also be obtained without them by lowering the efficiency of supernova feedback or the shock-heating mass. This model contains a free parameter that sets the dynamical friction efficiency, producing orphan satellites only when this parameter is larger than 0 (otherwise galaxies always merge when haloes and subhaloes merge). As a good fit to observations (not including clustering observations) was obtained without the need of orphan satellites, we decided to use the simplest solution and set the dynamical friction efficiency parameter to 0. Because of this, GALICS 2.0 has no orphans.

## 2.1.3 LGALAXIES (SAM)

In the original form of LGALAXIES the positions of orphans are followed by tracking the most bound particle of their host dark matter halo just before it was tidally disrupted. The orphan satellite is then placed not at the current position of the particle with which it is identified, but at a position whose (vector) offset from the central galaxy is reduced from that of the particle by a factor of $(1 - \delta t/t_{\text{friction}})$, where $\delta t$ is the time since the dynamical friction clock was started. This time dependence is based on a simple model for a satellite with 'isothermal' density structure spiralling to the centre of an isothermal host on a circular orbit (Guo et al. 2013; Henriques et al. 2013).

Since the dynamical information of most bound particles is not available for the current simulation, for this work LGALAXIES simply decays the positions of orphans from their value at the time they become orphans. Instead of $(1 - \delta t/t_{\text{friction}})$, a factor of $2 \times \sqrt{1 - \delta t/t_{\text{friction}}}$ is used in order to obtain satellite profiles and small-scale clustering that roughly resemble those from the default model.

## 2.1.4 MICE (HOD)

The galaxy population in haloes is determined from the halo mass, independently of their substructure. In the original implementation of the model, the luminosity function and the colour–magnitude diagrams are determined from observations (Blanton et al. 2005). Then, the galaxies are split into centrals and satellites. Using a modified NFW profile for the satellite distribution inside haloes, the occupation of galaxies as a function of halo mass is calibrated in order to reproduce the 2PCF of galaxies from observations. The modification of the NFW profile corresponds to a slightly steeper distribution that improves the clustering consistency with observations (Zehavi et al. 2011).

For this particular project, once the number of satellite galaxies in a halo is set, each satellite galaxy is assigned to a different subhalo. When there are more satellites than subhaloes, the excess of satellites are considered orphans and populate the halo according to a modified NFW profile. This is not the approach used in the original implementation of the model, where all the satellite galaxies are distributed with the modified NFW profile. In this case, we will be able to study the clustering of galaxies that follow subhaloes (the non-orphan satellites) and those that are consistent with the original implementation of the HOD model (the orphan satellites).

Although originally this model was implemented and calibrated in the MICE simulation (Crocce et al. 2015; Fosalba et al. 2015a,b) to make it consistent with clustering observations, in this simulation we used the same parameters obtained from the calibration in the MICE simulation, and hence the model is not necessarily reproducing the clustering observations in this study. In order to reproduce clustering observations, we would need to recalibrate the parameters of the model in this simulation and cosmology (the original simulation used the parameters $\Omega_m = 0.25$, $\Omega_\Lambda = 0.75$ and $\sigma_8 = 0.8$, while this simulation has $\Omega_m = 0.272$, $\Omega_\Lambda = 0.728$ and $\sigma_8 = 0.807$, to mention some parameters).

## 2.1.5 MORGANA (SAM)

This model has been originally designed to work with merger trees generated by the Lagrangian code PINOCCHIO (Monaco, Theuns & Taffoni 2002), and hence some adjustments have been needed in order to use it interface with SUBFIND-based merger trees.

As MORGANA does not explicitly follow the evolution of substructures, only central galaxies are linked to a dark matter structure. Whenever a dark matter halo becomes a substructure, its galaxies become satellites and each of them receive a merging time (computed from the Taffoni et al. 2003 prescriptions), which is defined independently from substructure evolution. As these merger times are estimated statistically, the merger of a satellite galaxy with the central object is decoupled from its parent substructure survival, i.e. satellite galaxies may merge before their host substructure is lost (giving rise to a population of substructures whose galaxies have already disappeared; the other models assign a residual merger time when the subhaloes is lost, assuring that no merger takes place before the subhalo disappears) or, vice versa, after it. In both cases they are considered as orphans, and placed at the centre of the host main halo (meaning that we do not track the trajectories of the orphans). As this assignment of the position of the orphans is not physical when describing small-scale clustering, we exclude them for the analysis of this model. The other non-central galaxies are considered non-orphan satellites for this analysis (differently than in K15, where all satellite galaxies of this model are considered orphans), and their position is then defined from their subhaloes.

Moreover, merging times computed from the Taffoni et al. (2003) prescriptions are typically shorter than those estimated from N-body simulations (De Lucia et al. 2010). Both effects have important implications on the satellite number density, showing a lower halo occupation number as we discuss later.

## 2.1.6 SAG (SAM)

When the subhaloes are no longer identified due to the mass loss form the merging with a larger structure, their galaxies become orphans. The trajectory of orphan satellites is calculated to be a circular orbit with a velocity determined by the virial velocity of
the host subhalo and initially located at a halocentric distance given by the virial radius of the subhalo. The decaying radial distance is estimated from the dynamical friction, with position and velocity components randomly generated. The orphan satellites finally merge with the central galaxy of the substructure in which they reside according to the dynamical friction time-scale (Binney & Tremaine 1987). Because of this time-scale, orphan satellites can be found inside the biggest substructure of the halo (where the central galaxy resides) or inside the substructure of another (satellite) galaxy.

### 2.1.7 SAGE (SAM)

When a halo/central galaxy system is captured by something larger to become a subhalo/satellite galaxy, the *expected* average merger time of the system is calculated using the Binney & Tremaine (1987) dynamical friction formula. The subhalo/satellite is then tracked with time until its subhalo-to-baryonic mass ratio falls below a critical threshold, taken as 1.0. At this point the current survival time as a subhalo/satellite is compared to the expected merger time calculated at infall. If the subhalo/satellite has survived longer than average we say it is more resistant to disruption and the satellite is merged with the central in the usual way. Otherwise the satellite is disrupted and its stars are added to a new ‘intracluster’ mass reservoir. As a consequence, SAGE does not produce an orphan galaxy population, since the decision about (and implementation of) the ultimate fate of a satellite is always made before (or when) its subhalo is lost in the merger tree.

This model is an update of Croton et al. (2006), and the suppression of orphan galaxies and the satellite treatment are some of the changes. Satellite galaxies in Croton et al. (2006) were found to be too red, mainly because of the instantaneous hot gas stripping that was causing a premature suppression of star formation. In SAGE satellite galaxies are treated more like central galaxies, in the sense that hot gas stripping now happens in proportion to the subhalo mass stripping. The lack of orphans makes the model to be resolution dependent, since the population of satellite galaxies depends on the resolution and detection of subhaloes in the simulations. However, most modern cosmological simulations have sufficiently high resolution to recover the galaxy population down to the limit of that typically probed by current surveys (Croton et al. 2016).

### 2.1.8 SKIBBASUBS (HOD)

As in MICE, the occupation of galaxies in haloes is determined as a function of the halo mass, independently of the substructure, and calibrated to recover luminosity function, colour–magnitude diagrams and the clustering 2PCF (Blanton et al. 2003; Zehavi et al. 2005). Satellite galaxies are distributed in subhaloes, and the exceeding galaxies are considered orphan satellites and populate the haloes according to the NFW profile. For an additional comparison in this study, we also constructed a catalogue where all the satellite galaxies follow the NFW profile (independently of the substructure of the haloes), that we call SKIBBAHOD. The difference between SKIBBASUBS and SKIBBAHOD is that the satellite galaxies have a different density profile. In SKIBBASUBS the satellite galaxies follow the subhaloes (except the orphan satellites) while in SKIBBAHOD all satellite galaxies follow an NFW profile. This is useful to understand the importance of the different treatments of the satellite distribution and how they affect clustering. In this study SKIBBAHOD will be treated as a reference for an NFW-based model.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Comparison of the orphan fraction (with respect to all the galaxies) of the different galaxy formation models as a function of the stellar mass. Each line corresponds to a different model. SAGE and GALICS 2.0 did not consider orphans in this work.</figcaption></figure>

As in MICE, the parameters have been calibrated in another simulation and cosmology, and hence the model does not necessarily reproduce clustering observations even they were fit to do it in the original simulation.

### 2.1.9 YSAM (SAM)

All subhaloes are tracked even after the halo finder loses them in the central dense region of a main halo, and populated with a resident galaxy. In these cases, galaxies only merge when they are closer to the centre of the halo than $0.1R_{vir}$. If a substructure disappears before reaching the central region of its host halo, the galaxy is considered orphan and YSAM calculates its mass (Battin 1987) and orbit (Binney & Tremaine 1987) analytically until approaching the very central regions. This has a large impact on the lifetime of subhaloes and galaxy merging time-scale (Yi et al. 2013).

## 2.2 Orphan fraction

In Fig. 1 we show the orphan fraction as a function of stellar mass for the galaxy formation models that have orphan satellites. This figure is similar to fig. 12 in K15 where the same fraction has been plotted, but as a function of halo mass. The orphan fraction is defined as

$$f_{\text{orph}}(M_{*}) = \frac{N_{\text{orph}}(M_{*})}{N(M_{*})}, \qquad (1)$$

where $N_{\text{orph}}(M_{*})$ is the number of orphan satellites in the catalogue with stellar mass $M_{*}$ and $N(M_{*})$ is the total number of galaxies with the same stellar mass. We see that the fraction decreases with mass. This trend is expected since small subhaloes are more easily affected by tidal stripping, so orphan satellites tend to originate from small subhaloes. As the mass of the orphan satellites is strongly related to the mass of the subhaloes at the time of accretion, their masses then tend to be small. Because of this, orphan satellites are more important at small masses than at large masses, and this implies that the role of orphan satellites on galaxy clustering will mostly be important for low-mass thresholds (Gao et al. 2004; Wang et al. 2006; Guo et al. 2011; Budzynski et al. 2012; Lee et al. 2014). The large orphan fractions of DLB07, GALACTICUS and SAG at $M_{*} > 3 \times 10^{11} h^{-1} M_{\odot}$ are not significant due to the low number of galaxies with these masses, that makes very few orphan satellites (less than five) represent a large fraction.
We can see that the fraction of orphan satellites depend strongly on the different galaxy formation model, with the scatter between the models being large. The models that show the lowest orphan fractions (for high masses) are YSAM and MORGANA, as expected from their treatment of orphan satellites. On one hand, YSAM tracks the galaxies in the substructures even after its mass has decreased below the resolution level. These galaxies are still considered non-orphan galaxies, and because of this it is more difficult for a galaxy to become orphan in this model. On the other hand, in MORGANA the galaxy merger times are shorter than those estimated from N-body simulations, which means that orphan satellites merge more quickly with central galaxies. Therefore there are fewer orphans for this model. Interestingly, GALACTICUS and GALFORM show the highest orphan fractions. Several studies (Contreras et al. 2013; Campbell et al. 2015; Simha & Cole 2016) show that the analytical equation used in GALFORM allows galaxies to orbit around the central galaxy for longer than other approximations, also causing a more centrally concentrated distribution of satellite galaxies. Note that the models that track the positions of the orphan satellites, and hence the models that we will use to measure the clustering of orphan satellites, are the ones that present the lowest orphan fraction at small masses. This means that the results of the distribution of orphan satellites that we show in this study might have a stronger impact on galaxy clustering for the rest of the models.

## 3 METHODOLOGY

To study the clustering between the models, we will use the 2PCF, which describes the excess of probability dP over a random distribution of finding pairs of galaxies at a given separation r:

$$dP = n(1 + \xi(r))dV, \qquad (2)$$

where $\xi(r)$ represents the 2PCF at a separation r and n is the number density of galaxies. There are several estimators of the 2PCF (Kerscher, Szapudi & Szalay 2000; Coil 2013), but for our study we use the estimator described by the following formula:

$$\xi(r) = \frac{DD(r)}{RR(r)} - 1, \qquad (3)$$

where $\xi(r)$ is the 2PCF as a function of scale, $DD(r)$ is the number of data pairs separated a distance r between them and $RR(r)$ is the number of random pairs at the same distance. $DD(r)$ and $RR(r)$ are normalized by $n_D(n_D - 1)$ and $n_R(n_R - 1)$, respectively, where $n_D$ and $n_R$ are the numbers of data and random points used. This estimator is equivalent to Landy & Szalay (1993) when the random sample is large enough, as it is here, where we use $n_R = 10^6$. Because of the size of the simulations, we calculate $\xi(r)$ up to $R = 6 h^{-1}$ Mpc, since the measurement becomes noisy for larger scales. Because of the resolution, the minimum scale for the study of $\xi(r)$ is $R = 300 h^{-1}$ kpc (Guo et al. 2011).

To calculate the errors of the 2PCF we use the jackknife method (Norberg et al. 2009). We divide the simulation box into 64 cubic subvolumes, and we measure the 2PCF 64 times excluding each time one of the subsamples. We obtain the error from these measurements using the unbiased standard deviation according to the following formula:

$$\Delta\xi(r) = \sqrt{\left(\frac{N_{JK} - 1}{N_{JK}}\right) \sum_{i=1}^{N_{JK}} (\xi_i(r) - \bar{\xi}(r))^2}, \qquad (4)$$

where $N_{JK}$ is the number of jackknife subsamples used and $\xi_i(r)$ corresponds to the measurement of $\xi(r)$ excluding the ith subsample.

The error bars, computed with jackknife, give an idea of the scatter that we would expect from different realizations of the same volume and number density. However, in this study we compare different models run on the same haloes, and then these errors do not reflect the uncertainties of the scatter between the models. Differences between models are systematic and could be significant even when they are below the error bars. Then, differences below the error bar must be taken with care in this study since the same comparison applied to a larger volume could reduce the error bars but not necessarily the difference between the models.

## 4 RESULTS

In this section we present the model comparison of the mean halo occupation number, the 2PCFs and the radial distribution of the galaxies in the haloes, with an emphasis on the orphan galaxy distribution. We use the galaxy catalogues at redshift $z = 0$ and apply different stellar mass cuts in order to see the mass dependence of the convergence and differences of the models.

### 4.1 Mean halo occupation number

In this section we study the mean number of galaxies populating haloes within a given range in mass. In the halo model paradigm it is usually assumed that haloes cluster according only to their masses. Hence, the distribution of galaxies provides a window into the clustering.

In both SAM and HOD models the galaxy populations are commonly characterized by central and satellite galaxies. According to these models, every halo can be occupied by at most one central galaxy, and only the haloes that contain a central galaxy can have a non-zero number of satellites:

$$\langle N|M, M_* \rangle = \langle N_{\text{cen}}|M, M_* \rangle (1 + \langle N_{\text{sat}}|M, M_* \rangle), \qquad (5)$$

with

$$\langle N_{\text{cen}}|M, M_* \rangle < 1, \qquad (6)$$

where $\langle N|M, M_* \rangle$ is the mean number of galaxies N of stellar mass $M_*$ that populate haloes of mass M, and $N_{\text{cen}}$ and $N_{\text{sat}}$ are the number of central and satellite galaxies, respectively. According to all this, for $\langle N|M \rangle \lesssim 1$ the contribution of the HOD mainly comes from the central galaxies and for $\langle N|M \rangle > 1$ the contribution mainly comes from the satellite galaxies.

In Fig. 2 we show the comparison of the mean halo occupation numbers of the models. This measurement corresponds to the mean number of galaxies per halo as a function of halo mass. Left-hand panel shows galaxies with $M_* > 10^9 h^{-1} M_{\odot}$, and right-hand panel shows galaxies with $M_* > 10^{10} h^{-1} M_{\odot}$. This figure gives similar information as figs 13 and 14 from K15, where the number of galaxies is normalized by the halo mass to explore the specific frequency of galaxies as a function of halo mass.

We can see a scatter where each model starts populating galaxies for low stellar masses, which is a consequence of the different implementations of cooling, reionization and stellar feedback. In particular, the minimum mass where all the haloes are populated (i.e. where $\langle N(M) \rangle = 1$) changes a factor of 3 between the models for galaxies of a stellar mass threshold of $M_* > 10^9 h^{-1} M_{\odot}$. This scatter decreases to a factor of 2 when the $M_* > 10^{10} h^{-1} M_{\odot}$ cut is applied, except for GALACTICUS.

If we focus on the right-hand panel, we see a strong difference in GALACTICUS for $M_* > 10^{10} h^{-1} M_{\odot}$, where these massive galaxies also populate very small haloes. This is due to the excess of galaxies
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Mean halo occupation number as a function of halo mass for the different models using all the galaxy types. Left-hand panel shows galaxies with $M_{*} > 10^{9} h^{-1} M_{\odot}$, while the right-hand panel shows galaxies with $M_{*} > 10^{10} h^{-1} M_{\odot}$.</figcaption></figure>

at these masses for this model, which can be seen as a bump around $M_{*} \sim 2-3 \times 10^{10} h^{-1} M_{\odot}$ in the stellar mass function from fig. 2 in K15. This comes from the fact that the galaxy formation model has been calibrated using another simulation to match observations. Changing the simulation without recalibrating the stellar mass function has a significant impact on GALACTICUS (see fig. 6 of K15). We also like to mention that the fact that SKIBBAHOD and SKIBBASUBS are not identical due to some stochastic components of the models. We also note that SAGE and GALICS 2.0 show the lowest occupation number at high masses in the left-hand panel. This is expected since these models do not have orphan satellites by construction.

In order to study the contributions of the different galaxy types, we show in Fig. 3 the mean occupation number of galaxies split into galaxy types. Top panels show central galaxies, middle panels show satellite galaxies (orphan and non-orphan) and bottom panels show orphan satellites. The same stellar mass cuts as in Fig. 2 has been applied here for left- and right-hand panels.

Given the fact that all the haloes are populated by a central galaxy, the mean occupation number for a given stellar mass cut is directly related to the fraction of central galaxies that are more massive than the given stellar mass cut. Models whose central galaxies are less massive will present lower occupation numbers when a stellar mass cut is applied. If we focus on the top left-hand panel where the central galaxies are shown for a stellar mass cut of $10^{9} h^{-1} M_{\odot}$, we see that SAG, GALACTICUS, GALICS 2.0 and DLB07 are the models that show higher occupation numbers, while LGALAXIES is the model that has the lowest values. This is consistent with table 3 of K15, where we see that SAG, GALACTICUS, GALICS 2.0 and DLB07 are precisely the models that present more central galaxies above $10^{9} h^{-1} M_{\odot}$, and LGALAXIES presents the lowest number.

If we focus on the middle panels, where satellite galaxies are plotted, we see that the models show a large scatter at small masses, but the number of galaxies per halo increases with mass with a similar slope for larger masses. This agreement in the slope of the relation means that the galaxy formation models distribute satellite galaxies in haloes in a similar way. We note the different behaviour shown by MORGANA, which shows a significantly lower occupation number for small haloes. This is due to the decoupled modelling of satellite galaxies with respect to substructures that can leave some naked substructures without any satellite galaxy in it. The reason of this treatment of satellite galaxies is that the model was originally designed to be run in PINOCCHIO simulations, where

only haloes (but not subhaloes) are obtained from the output of the simulations.

The differences between HOD models and SAMs are stronger on the bottom panels, where we show the mean halo occupation number for orphan satellites. We can see that the slope in the mean halo occupation numbers of HOD models is much shallower than SAMs for $M_{*} > 10^{12} h^{-1} M_{\odot}$, showing a difference between HOD models and SAMs of one order of magnitude higher at $M_{*} \approx 10^{14} h^{-1} M_{\odot}$ than at $M_{*} \approx 10^{12} h^{-1} M_{\odot}$. This is because massive haloes have many substructures and hence the HOD models occupy them with as many non-orphan satellites as possible. If the total occupation number is not high enough, then the number of orphan satellites in these haloes is low. On the contrary, orphan satellites in SAMs originate from the disruption of subhaloes, and this usually happens in high-density environments with strong gravitational interactions. Consequently, in massive haloes, many subhaloes can interact with the environment and suffer tidal stripping. Hence, we expect that the orphan occupation increases quickly with halo mass, as we can see from these panels.

Finally, note also that SKIBBAHOD has a higher amplitude than SKIBBASUBS, and the slope is similar to some SAMs at high enough halo masses. This is due to the fact that SKIBBAHOD populates the haloes only with orphan satellites by construction, for which the total number is significantly higher than SKIBBASUBS. Also see that most of the SAMs show a good agreement in this mean halo occupation number of orphan satellites for $M_{*} > 10^{9} h^{-1} M_{\odot}$, especially for DLB07, GALACTICUS, GALFORM and SAG, and there is a good agreement between LGALAXIES and YSAM too. However, the differences become more significant for $M_{*} > 10^{10} h^{-1} M_{\odot}$. MORGANA shows the lowest occupation number due to the shorter merger times implemented in this model.

## 4.2 Radial distributions

In this section we compare the radial distributions measured in all the models. We do the measurement from the following equation:

$$n(R/R_{200}) = \frac{N_g}{(4\pi/3)[(R + \Delta R)^3 - R^3]}, \quad (7)$$

where $n(R/R_{200})$ is the number density of galaxies in the radial annulus $R$ to $R + \Delta R$, with $R$ referring to the radial distance to the halo centre, and $N_g$ is the number of galaxies between $R$ and $R + \Delta R$.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Mean halo occupation number as a function of halo mass for the different models and galaxy types. Left-hand panels show galaxies with $M_* > 10^9 h^{-1} M_{\odot}$, while the right-hand panels show $M_* > 10^{10} h^{-1} M_{\odot}$. Top panels show central galaxies, middle panels satellite galaxies and bottom panels orphan satellites.</figcaption></figure>

$R_{200}$ is the radius that encloses 200 times the critical density. So, this basically describes the density of galaxies as a function of the radial distance of the halo centre.

Fig. 4 shows the radial distribution of the different models, applying the same stellar mass thresholds of $M_* > 10^9 h^{-1} M_{\odot}$ (left) and $M_* > 10^{10} h^{-1} M_{\odot}$ (right) as in the previous figure. Top panels show all the satellites, while bottom panels show only orphan satellites.

Each panel shows the residual as $n(R)/\bar{n}(R) - 1$, where $\bar{n}(R)$ is the median of the distribution at each R bin.

In the top panels we only show those models that have implemented a treatment for orphan satellites (or have not done it by construction). The differences between most of the models are lower than $2\sigma$ for $M_* > 10^9 h^{-1} M_{\odot}$ and lower than $1\sigma$ for $M_* > 10^{10} h^{-1} M_{\odot}$. We see that the scatter is smaller than a factor
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of the radial distributions of the different models. Left-hand panels show galaxies with $M_{*} > 10^{9} h^{-1} M_{\odot}$, while in the right-hand panels with $M_{*} > 10^{10} h^{-1} M_{\odot}$. In the top panels we show all the satellite galaxies, while bottom panels show only orphan satellites. Each panel includes the residual with respect to the median of the distribution in each $R/R_{200}$ bin.</figcaption></figure>

of 2 for all scales smaller than $0.2R/R_{200}$ and for all the models except SKIBBAHOD, GALICS 2.0 and SAGE. SKIBBAHOD shows the highest radial distribution, a factor of 2 higher than the median. This is because all the satellites are orphans, and hence all the galaxies follow a NFW profile instead of following substructures. The fact that most of the models agree with the radial distribution of SKIBBAHOD except for the smallest scales is consistent with previous studies (Gao et al. 2004; Kang 2014; Pujol et al. 2014; van Daalen et al. 2016). On the other hand, SAGE and GALICS 2.0 show much flatter radial distributions than the rest of the models. This is because these models have no orphans, and hence only non-orphan galaxies contribute to these distributions. These two models then show the contribution of non-orphan satellite galaxies to the radial distributions, showing agreement with the rest of the models only at the largest scales.

We see that the HOD models (MICE, SKIBBAHOD and SKIBBASUBS) present a steeper slope of the radial distributions than SAMs in all the cases. Given the similarities between the top and bottom panels, we see that these differences basically come from the different treatments of orphan satellites, which dominate the smallest scales (satellites only contribute to the large scales of the panels). This indicates that the orphan satellites have a very important role in the distribution of galaxies in haloes, consistent with the conclusions at Gao et al. (2004). While HOD models distribute orphan satellites without any information from substructure or evolution (in SKIBBAHOD

this is the case for all satellites), orphan satellites from SAMs are a consequence of subhalo disruption, and hence the positions of orphan satellites are correlated with substructure. Moreover, orphan satellites in SAMs are limited in the densest regions, close to the halo centre, where orphans merge quickly with the central galaxy. This exclusion effect in the inner parts of the halo is one reason of why SAMs show a lower orphan density at the smallest scales.

## 4.3 Two-point galaxy correlation functions

We now compare the 2PCF between the different models. We have applied the same stellar mass thresholds used previously, and we also study the different galaxy types separately. Again, we show the residuals with respect to the median of the measurements in order to see the scatter between the models. Since the stellar mass functions are different between the models, the number density of galaxies for the same stellar mass cut can be different. We have also studied the differences using number density cuts instead of stellar mass cuts and we obtain the same results, meaning that the difference between the models is not due to the differences between their number densities.

In Fig. 5 we show the 2PCFs of galaxies for all the models that have computed the orphan positions (and those that did not do it by construction), using all the galaxies. For both stellar mass cuts we find a good agreement, with most of the models consistent within the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of the 2PCF measurements of the different galaxy formation models using a stellar mass threshold. Left-hand panel shows galaxies with $M_{*} > 10^{9} h^{-1} M_{\odot}$, while right-hand panel corresponds to $M_{*} > 10^{10} h^{-1} M_{\odot}$. Both panels use all the galaxies of the models.</figcaption></figure>

error bars (although we must be careful when interpreting the error bars, as discussed in Section 3). This is an encouraging result, since it highlights a consistency between the models even when most of them did not use any observations of clustering to constrain their parameters. However, we note that GALICS 2.0 and SAGE show a lower clustering at small scales, showing a factor of 2 lower at the smallest scales. Again, this is a consequence of the fact that these models do not have orphan satellites. These models then show the impact that excluding orphan satellites can have on the galaxy clustering predictions.

We split the galaxy samples according to their type in Fig. 6. We show in the top and middle panels the 2PCF of galaxies for central and satellite galaxies, respectively. In the cases of central galaxies the models have very good agreement, while some differences appear for satellite galaxies. The scatter between models in the middle panels is larger for smaller stellar masses, since galaxies are more dominated by orphan galaxies. In particular, all the satellites in SKIBBAHOD are orphans, causing a large difference with the rest of the models at the smallest scales due to the orphan radial distributions discussed previously. We also note a high clustering signal for SKIBBASUBS, and a lower signal for GALICS 2.0 (which has no orphans) and SAG. For most of the models, the scatter between the models is in general lower than a 25 per cent for the satellite galaxies and lower than 20 per cent for the central galaxies.

In the bottom panels of Fig. 6 we focus on the models that have orphan satellites in order to study their distribution. For both stellar mass thresholds we see a strong and significant difference between the HOD models and SAMs. The HOD models show a lower clustering amplitude, and they all agree between them, while SAMs agree between them but with a higher amplitude and different shape than the HODs. The clustering of 2PCF on small scales depends on two main factors, the halo occupation number and the density profile of galaxies in haloes, especially in massive haloes. Although the HOD models show a slightly steeper radial distribution of orphan satellites with respect to SAMs, they also show a much flatter halo occupation number distribution, as indicated in the lower panel of Fig. 3. This implies that orphan satellites in HODs populate less massive haloes than in SAMs. This is the main reason for the lower 2PCF of orphan satellites in HOD models compared to SAMs, since small-scale clustering is strongly affected by the occupation numbers in massive haloes. This difference between the clustering of orphan satellites of SAMs and HOD models is large, reaching an order of magnitude at scales of $\approx 1 h^{-1} Mpc$. The impact of orphan

satellites in the agreement between models will depend in general on the orphan fraction of the galaxy samples, and we have seen that it also depends on the galaxy formation model. In this analysis the models that computed the orbits and positions of the galaxies are those with the lowest orphan fractions, meaning that the impact of orphan satellites on the clustering of the rest of the models might be stronger.

## 5 DISCUSSION AND CONCLUSIONS

In this paper, we present a comparison of the clustering and halo occupation statistics of 12 different galaxy formation models. We use a dark matter only N-body simulation and run SAMs and HOD-based models with the same dark matter and merger trees input, and compare the results of mean halo occupation numbers, radial distributions of galaxies in haloes and 2PCFs. The goal of this paper is to study the clustering and distribution of galaxies in haloes, and to understand the roles of different galaxy types, in particular of orphan satellites (satellites that are not assigned to any dark matter subhalo). This work is part of a series of papers comparing galaxy formation models that started with K15.

The most important results of the study can be summarized as follows.

1. The slope in the mean occupation number of orphan satellites as a function of halo mass is much shallower in HOD models than in SAMs, due to the different treatments of orphan satellites between both approaches. Orphan satellites in SAMs originate from the disruption of subhaloes, and this happens more often in massive haloes. However, in this study most of the HOD models populate satellites in subhaloes and only when there are more satellites than subhaloes these extra galaxies are considered as orphans. As massive haloes have many subhaloes, the number of orphan satellites in these HOD models is not as high as in SAMs.

2. HOD models have a steeper radial distribution of orphan satellites in haloes than SAMs. This is because HOD models distribute orphan satellites following a NFW profile, independently of the substructure and evolution of the haloes. This allows HOD models to populate with more orphan satellites in the inner and denser regions of the haloes than SAMs. SAMs are constrained to where subhaloes have been disrupted, and this causes a lower density in the innermost regions, where subhaloes quickly merge into the central structure, as well as galaxies merging with other
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of the 2PCF measurements of the different galaxy formation models using a stellar mass threshold. Left-hand panels show galaxies with $M_{*} > 10^{9} h^{-1} M_{\odot}$, while right-hand panels show galaxies with $M_{*} > 10^{10} h^{-1} M_{\odot}$. In the top panels we show the central galaxies, middle panels show all the satellite galaxies and bottom panels show orphan satellites.</figcaption></figure>

galaxies for some models. When comparing all the satellites, the different models have a scatter of two times the measurement uncertainty (due to the limited volume used) in their radial distribution of galaxies for $M_{*} > 10^{9} h^{-1} M_{\odot}$. The scatter is smaller for higher thresholds, since orphan satellites become less relevant.

(3) Using all the galaxies above a certain mass threshold for the measurements of 2PCFs (see Fig. 5), we see a scatter of a factor

of 2 between the models. However, part of this scatter is due to the lower clustering found for the models that do not have orphan galaxies by construction, an indication of the importance of orphan galaxies on galaxy clustering. We find a good agreement between the models for central galaxies and at large scales for all the galaxy selections. Using a larger volume in simulations would allow us to measure linear bias, and this would be a valuable extension to this work.
(4) HOD models and SAMs have significant differences in their clustering of orphan satellites (see bottom panels of Fig. 6). Both SAMs and HODs show good agreement for models of the same kind, but SAMs have a higher 2PCF than HOD models. This is due to the differences on the halo occupation numbers of orphan satellites between both schemes. Although HOD models show a steeper radial distribution for orphan galaxies than SAMs, they statistically occupy less massive haloes. The clustering at small scales is strongly affected by the halo occupation of massive haloes, and because of this the orphan satellites in SAMs show a higher 2PCF at these scales than HOD models.

It is important to notice that the models used have not been recal-
ibrated for this particular simulation. The agreement between the
models could be improved by calibrating the models in the simu-
lation where the comparison has been done or even using the same
observational constraints (Knebe et al., in preparation). This study
is limited by the resolution of the simulation. A higher resolution
simulation would allow us to study smaller scales, and would also
allow to detect subhaloes in inner regions of the haloes. This could
have an impact on both the satellite distributions at small scales. In
addition, a comparison with hydrodynamic simulations would be
useful to study the baryonic effects on both galaxy and dark matter
clustering. It has been shown that baryons affect the dark matter dis-
tribution at small scales (Tissera et al. 2010; Cui et al. 2012, 2016;
Sawala et al. 2013; Cui, Borgani & Murante 2014). Finally, an-
other interesting extension would be the comparison of SAMs with
new implementations that take into account observations of galaxy
clustering to constrain their parameters (van Daalen et al. 2016).

ACKNOWLEDGEMENTS

The authors would like to express special thanks to the Instituto de Física Teórica (IFT-UAM/CSIC in Madrid) for its hospitality and support, via the Centro de Excelencia Severo Ochoa Program under Grant no. SEV-2012-0249, during the 3 weeks workshop 'nIFTy cosmology' where this work developed. We further acknowledge the financial support of the 2014 University of Western Australia Research Collaboration Award for 'Fast Approximate Synthetic Universes for the SKA', the ARC Centre of Excellence for All Sky Astrophysics (CAASTRO) grant number CE110001020 and the two ARC Discovery Projects DP130100117 and DP140100198. We also recognize support from the Universidad Autonoma de Madrid (UAM) for the workshop infrastructure.

Funding for this project was partially provided by the Span-
ish Ministerio de Ciencia e Innovación (MICINN), Consolider-
Ingenio CSD2007-00060, European Commission Marie Curie Ini-
tial Training Network CosmoComp (PITNGA-2009-238356). We
acknowledge support from the European Commission's Frame-
work Programme 7, through the Marie Curie International Re-
search Staff Exchange Scheme LACEgal (PIRSES-GA-2010-
269264). AP acknowledges support from beca FI and 2009-SGR-
1398 from Generalitat de Catalunya, project AYA2012-39620 and
AYA2015-71825 from MICINN, and from a European Research
Council Starting Grant (LENA-678282). RAS acknowledges sup-
port from the NSF grant AST-1055081. FJC acknowledges support
from the Spanish Ministerio de Economía y Competitividad project
AYA2012-39620. SAC acknowledges grants from CONICET (PIP-
220), Argentina. DJC acknowledges receipt of a QEII Fellowship
from the Australian Government. PJE is supported by the SSimPL
programme and the Sydney Institute for Astronomy (SIfA) via
ARC grant, DP130100117. WC and CP acknowledge support of

ARC DP130100117. FF acknowledges financial support from the grants PRIN INAF 2010 ‘From the dawn of galaxy formation’ and PRIN MIUR 2012 ‘The intergalactic medium as a probe of the growth of cosmic structures’. JG-B is supported by Spain through the MINECO grant FPA2015-68048, as well as the Consolid- Ingenio 2010 Programme of MICINN under grant PAU CSD2009- 00060 and the Severo Ochoa Programme SEV-2012-0249. VG-P acknowledges support from a European Research Council Start- ing Grant (DEGAS-259586). This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the Na- tional E-Infrastructure. The work of BMBH was supported by a Zwicky Prize fellowship and by Advanced Grant 246797 ‘GAL- FORMOD’ from the European Research Council. MH acknowl- edges financial support from the European Research Council via an Advanced Grant under grant agreement no. 321323 NEOGAL. AK is supported by the Ministerio de Economía y Competitividad and the Fondo Europeo de Desarrollo Regional (MINECO/FEDER, UE) in Spain through grants AYA2012-31101 and AYA2015-63810- P and the Consolider-Ingenio 2010 Programme of the Spanish Ministerio de Ciencia e Innovación (MICINN) under grant Mul- tiDark CSD2009-00064. He also acknowledges support from the Australian Research Council (ARC) grant DP140100198. PM has been supported by a FRA2012 grant of the University of Trieste, PRIN2010-2011 (J91J12000450001) from MIUR and Consorzio per la Fisica di Trieste. NDP was supported by BASAL PFB-06 CATA, and Fondecyt 1150300. Part of the calculations presented here were run using the Geryon cluster at the Center for Astro- Engineering at U. Catolica, which received funding from QUIMAL 130008 and Fondequip AIC-57. CP acknowledges support of the Australian Research Council (ARC) through Future Fellowship FT130100041 and Discovery Project DP140100198. RSS thanks the Downsborough family for their generous support. PAT acknowl- edges support from the Science and Technology Facilities Council (grant number ST/L000652/1). SKY acknowledges support from the National Research Foundation of Korea (Doyak 2014003730). Numerical simulations were performed using the KISTI supercom- puter under the programme of KSC-2013-C3-015.

The authors of the paper contributed in the following ways: AP, RAS and EG lead the study and wrote the paper. AK, JG-B and FRP organized the second week of the nIFTy workshop from where this study began. The authors listed in Table 1 performed their galaxy formation models using their code, in particular AB, FJC, AC, SC, DC, GDL, FF, VG-P, BH, JL, PM, RAS, RS, CVM and SY actively ran their models with the assistance of JH, MH and CS. WC, DC, PJE, CP and JO assisted with the analysis and data format issues. All authors had the opportunity to proof read and comment on the paper.

REFERENCES

Battin R. H., 1987, An Introduction to the Mathematics and Methods of As-
trodynamics. AIAA Education Series. American Institute of Aeronautics
and Astronautics Inc., New York

Baugh C. M., 2006, Rep. Progress Phys., 69, 3101

Baugh C. M., 2013, Publ. Astron. Soc. Aust., 30, 30

Baugh C. M., Benson A. J., Cole S., Frenk C. S., Lacey C. G., 1999, MNRAS,
305, L21
Benson A. J., 2010, Phys. Rep., 495, 33

Benson A. J., 2012, New Astron., 17, 175

Benson A. J., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2001, MNRAS, 327, 1041

Berlind A. A., Weinberg D. H., 2002, ApJ, 575, 587

Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton Univ. Press, Princeton, NJ

Blanton M. R. et al., 2003, ApJ, 592, 819

Blanton M. R., Lupton R. H., Schlegel D. J., Strauss M. A., Brinkmann J., Fukugita M., Loveday J., 2005, ApJ, 631, 208

Budzynski J. M., Koposov S. E., McCarthy I. G., McGee S. L., Belokurov V., 2012, MNRAS, 423, 104

Campbell D. J. R. et al., 2015, MNRAS, 452, 852

Carretero J., Castander F. J., Gaztañaga E., Crocce M., Fosalba P., 2015, MNRAS, 447, 646

Cattaneo A., Dekel A., Devriendt J., Guiderdoni B., Blaizot J., 2006, MNRAS, 370, 1651

Cattaneo A., Blaizot J., Devriendt J., Mamon G. A., Tollet E., Guiderdoni B., Kucukbas M., Thob A., 2017, MNRAS, in press

Coil A. L., 2013, in Oswalt T. D., Keel W. C., eds, Planets, Stars, and Stellar Systems, Vol. 6. Springer-Verlag, Dordrecht, p. 387

Cole S., 1991, ApJ, 367, 45

Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168

Conroy C., Wechsler R. H., Kravtsov A. V., 2006, ApJ, 647, 201

Contreras S., Baugh C. M., Norberg P., Padilla N., 2013, MNRAS, 432, 2717

Cooray A., 2006, MNRAS, 365, 842

Cooray A., Sheth R., 2002, Phys. Rep., 372, 1

Cora S. A., 2006, MNRAS, 368, 1540

Crocce M., Castander F. J., Gaztañaga E., Fosalba P., Carretero J., 2015, MNRAS, 453, 1513

Croton D. J. et al., 2006, MNRAS, 365, 11

Croton D. J. et al., 2016, ApJS, 222, 22

Cui W., Borgani S., Dolag K., Murante G., Tornatore L., 2012, MNRAS, 423, 2279

Cui W., Borgani S., Murante G., 2014, MNRAS, 441, 1769

Cui W. et al., 2016, MNRAS, 458, 4052

De Lucia G., Blaizot J., 2007, MNRAS, 375, 2

De Lucia G., Kauffmann G., Springel V., White S. D. M., Lanzoni B., Stoehr F., Tormen G., Yoshida N., 2004, MNRAS, 348, 333

De Lucia G., Boylan-Kolchin M., Benson A. J., Fontanot F., Monaco P., 2010, MNRAS, 406, 1533

De Lucia G., Fontanot F., Wilman D., Monaco P., 2011, MNRAS, 414, 1439

Diemand J., Moore B., Stadel J., 2004, MNRAS, 352, 535

Diemand J., Kuhlen M., Madau P., 2007, ApJ, 667, 859

Fontanot F., De Lucia G., Monaco P., Somerville R. S., Santini P., 2009, MNRAS, 397, 1776

Fontanot F., De Lucia G., Wilman D., Monaco P., 2011, MNRAS, 416, 409

Fontanot F., Cristiani S., Santini P., Fontana A., Grazian A., Somerville R. S., 2012, MNRAS, 421, 241

Fosalpa B., Gaztañaga E., Castander F. J., Crocce M., 2015a, MNRAS, 447, 1319

Fosalpa B., Crocce M., Gaztañaga E., Castander F. J., 2015b, MNRAS, 448, 2987

Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004, MNRAS, 355, 819

Gargiulo I. D. et al., 2015, MNRAS, 446, 3820

Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 2000, ApJ, 544, 616

Gill S. P. D., Knebe A., Gibson B. K., Dopita M. A., 2004, MNRAS, 351, 410

Giocoli C., Tormen G., van den Bosch F. C., 2008, MNRAS, 386, 2135

Gonzalez-Perez V., Lacey C. G., Baugh C. M., Lagos C. D. P., Helly J., Campbell D. J. R., Mitchell P. D., 2014, MNRAS, 439, 264

Guo Q. et al., 2011, MNRAS, 413, 101

Guo Q., White S., Angulo R. E., Henriques B., Lemson G., Boylan-Kolchin M., Thomas P., Short C., 2013, MNRAS, 428, 1351

Guo H. et al., 2016, MNRAS, 459, 3040

Hatton S., Devriendt J. E. G., Ninin S., Bouchet F. R., Guiderdoni B., Vibert D., 2003, MNRAS, 343, 75

Hearin A. P., Zentner A. R., Berlind A. A., Newman J. A., 2013, MNRAS, 433, 659

Henriques B. M. B., White S. D. M., Thomas P. A., Angulo R. E., Guo Q., Lemson G., Springel V., 2013, MNRAS, 431, 3373

Kang X., 2014, MNRAS, 437, 3385

Kang X., Li M., Lin W. P., Elahi P. J., 2012, MNRAS, 422, 804

Kerscher M., Szapudi I., Szalay A. S., 2000, ApJ, 535, L13

Kimm T. et al., 2009, MNRAS, 394, 1131

Knebe A. et al., 2015, MNRAS, 451, 4029 (K15)

Knollmann S. R., Knebe A., 2009, ApJS, 182, 608

Kravtsov A. V., Berlind A. A., Wechsler R. H., Klypin A. A., Gottlöber S., Allgood B., Primack J. R., 2004, ApJ, 609, 35

Lacey C., Silk J., 1991, ApJ, 381, 14

Landy S. D., Szalay A. S., 1993, ApJ, 412, 64

Lee J., Yi S. K., 2013, ApJ, 766, 38

Lee J. et al., 2014, MNRAS, 445, 4197

Lo Faro B., Monaco P., Vanzella E., Fontanot F., Silva L., Cristiani S., 2009, MNRAS, 399, 827

Mo H. J., White S. D. M., 1996, MNRAS, 282, 347

Mo H., van den Bosch F. C., White S., 2010, Galaxy Formation and Evolution. Cambridge Univ. Press, Cambridge

Monaco P., Theuns T., Taffoni G., 2002, MNRAS, 331, 587

Monaco P., Fontanot F., Taffoni G., 2007, MNRAS, 375, 1189

Norberg P., Baugh C. M., Gaztañaga E., Croton D. J., 2009, MNRAS, 396, 19

Pujol A. et al., 2014, MNRAS, 438, 3205

Reddick R. M., Wechsler R. H., Tinker J. L., Behroozi P. S., 2013, ApJ, 771, 30

Sawala T., Frenk C. S., Crain R. A., Jenkins A., Schaye J., Theuns T., Zavala J., 2013, MNRAS, 431, 1366

Scoccimarro R., Sheth R. K., Hui L., Jain B., 2001, ApJ, 546, 20

Seljak U., 2000, MNRAS, 318, 203

Sheth R. K., Tormen G., 2002, MNRAS, 329, 61

Simha V., Cole S., 2016, MNRAS, preprint (arXiv:1609.09520)

Skibba R. A., Sheth R. K., 2009, MNRAS, 392, 1080

Skibba R., Sheth R. K., Connolly A. J., Scranton R., 2006, MNRAS, 369, 68

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Somerville R. S., Primack J. R., 1999, MNRAS, 310, 1087

Springel V., 2005, MNRAS, 364, 1105

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726

Taffoni G., Mayer L., Colpi M., Governato F., 2003, MNRAS, 341, 434

Tissera P. B., White S. D. M., Pedrosa S., Scannapieco C., 2010, MNRAS, 406, 922

Vale A., Ostriker J. P., 2006, MNRAS, 371, 1173

van Daalen M. P., Henriques B. M. B., Angulo R. E., White S. D. M., 2016, MNRAS, 458, 934

van den Bosch F. C., 2002, MNRAS, 331, 98

van den Bosch F. C., Tormen G., Giocoli C., 2005, MNRAS, 359, 1029

Wang L., Li C., Kauffmann G., De Lucia G., 2006, MNRAS, 371, 537

White S. D. M., Frenk C. S., 1991, ApJ, 379, 52

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

Yang X., Mo H. J., van den Bosch F. C., 2003, MNRAS, 339, 1057

Yi S. K., Lee J., Jung I., Ji I., Sheen Y.-K., 2013, A&A, 554, A122

Zehavi I. et al., 2005, ApJ, 630, 1

Zehavi I. et al., 2011, ApJ, 736, 59

## APPENDIX A: IMPACT OF HALO MASS DEFINITION

Halo mass is one of the properties that galaxy formation models use to determine the population and properties of galaxies. Because of this the models might produce different results if they use different definitions of halo mass. In order to make fair comparisons
<figure><img src="image_9.png" /><figcaption>**Figure A1.** 2PCF for the GALACTICUS model for galaxies with $M_{*} > 10^{10} h^{-1} M_{\odot}$. Each line represents a different halo mass definition used for the model.</figcaption></figure>

in our analysis we need to study the dependence on the halo mass definition.

In Fig. A1 we show the 2PCF for different halo mass definitions using the GALACTICUS model for galaxies with $M_{*} > 10^{10} h^{-1} M_{\odot}$. We only show this model because it used all the different mass definitions and also because we expect the other models to show similar behaviour. Although we only present one model, the different mass definitions reveal only very small changes to the clustering compared with the differences between the models. Thus, the clustering of the models is not affected by the halo mass definition significantly. This result is also independent of the stellar mass selection used. Hence, the results of our study do not depend on the masses used and we focus on few definitions. In this paper we use $M_{200m}$ (defined as the mass enclosed in a radius within the density 200 times the mean density) in all the models that used this mass to obtain a catalogue, and $M_{200c}$ (defined as the mass enclosed in a radius within the density is 200 times the critical density) or $M_{FOF}$ (defined from the total number of particles belonging to the FOF group) for the other models, which do not use $M_{200m}$. The results and conclusions of the paper do not depend on the mass definitions used.

<footnote n="1">¹Institut de Ciències de l'Espai, IEEC-CSIC, Campus UAB, E-08193 Bellaterra, Barcelona, Spain</footnote>

<footnote n="2">²CosmoStat Laboratory, DRF/IRFU/SEDI-Service d'Astrophysique, CEA Saclay, F-91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="3">³Department of Physics, Center for Astrophysics and Space Sciences, University of California, 9500 Gilman Drive, San Diego, CA 92093, USA</footnote>

<footnote n="4">⁴University of California, Santa Cruz, Science Communication Program, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="5">⁵Carnegie Observatories, 813 Santa Barbara Street, Pasadena, CA 91101, USA</footnote>

<footnote n="6">⁶Université de Lyon, F-69003 Lyon, France</footnote>

<footnote n="7">⁷Université Lyon 1, Observatoire de Lyon, 9 avenue Charles Andrè, F-69230 Saint-Genis Laval, France</footnote>

<footnote n="8">⁸CNRS, UMR 5574, Centre de Recherche Astrophysique de Lyon, École Normale Supérieure de Lyon, F-69007 Lyon, France</footnote>

<footnote n="9">⁹Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="10">¹⁰Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra, Barcelona, Spain</footnote>

<footnote n="11">¹¹GEPI, Observatoire de Paris, CNRS, 61 Avenue de l'Observatoire, F-75014 Paris, France</footnote>

<footnote n="12">¹²Instituto de Astrofísica de La Plata (CCT La Plata, CONICET, UNLP), Paseo del Bosque s/n, B1900FWA La Plata, Argentina</footnote>

<footnote n="13">¹³Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque s/n, B1900FWA La Plata, Argentina</footnote>

<footnote n="14">¹⁴Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="15">¹⁵Departamento de Física Teórica, Módulo 15, Facultad de Ciencias, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="16">¹⁶South African Astronomical Observatory, PO Box 9, Observatory, Cape Town 7935, South Africa</footnote>

<footnote n="17">¹⁷Department of Physics and Astronomy, University of the Western Cape, Cape Town 7535, South Africa</footnote>

<footnote n="18">¹⁸INAF – Astronomical Observatory of Trieste, via Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="19">¹⁹Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="20">²⁰Sydney Institute for Astronomy, A28 School of Physics, The University of Sydney, NSW 2006, Australia</footnote>

<footnote n="21">²¹Astrophysics Research Institute, Liverpool John Moores University, IC2 Liverpool Science Park, 146 Brownlow Hill, Liverpool L3 5RF, UK</footnote>

<footnote n="22">²²Instituto de Física Teórica, Universidad Autónoma de Madrid (IFT-UAM/CSIC), E-28049 Madrid, Spain</footnote>

<footnote n="23">²³Institute of Cosmology and Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth PO1 3FX, UK</footnote>

<footnote n="24">²⁴Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching b. München, Germany</footnote>

<footnote n="25">²⁵Institute for Astronomy, ETH Zurich, CH-8093 Zurich, Switzerland</footnote>

<footnote n="26">²⁶Institut d'Astrophysique de Paris (UMR 7095: CNRS and UPMC), 98 bis Bd Arago, F-75014 Paris, France</footnote>

<footnote n="27">²⁷Astro-UAM, UAM, Unidad Asociada Consejo Superior de Investigaciones Científicas</footnote>

<footnote n="28">²⁸Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 120-749, Korea</footnote>

<footnote n="29">²⁹Dipartimento di Fisica, Università di Trieste, via Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="30">³⁰School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK</footnote>

<footnote n="31">³¹International Centre for Radio Astronomy Research, University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia</footnote>

<footnote n="32">³²Department of Physics and Astronomy, Rutgers University, 136 Frelinghuysen Road, Piscataway, NJ 08854, USA</footnote>

<footnote n="33">³³Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote>This paper has been typeset from a TEX/LATEX file prepared by the author.</footnote>