BASILISK II. Improved Constraints on the Galaxy-Halo Connection from Satellite Kinematics in SDSS
Kaustav Mitra Frank C. van den Bosch Johannes U. Lange




BASILISK II. Satellite Kinematics in SDSS]
      BASILISK II. Improved Constraints on the Galaxy-Halo Connection from Satellite Kinematics in SDSS

Mitra et al.]
   Kaustav Mitra^1E-mail: kaustav.mitra@yale.edu0000-0001-8073-4554,
   Frank C. van den Bosch^10000-0003-3236-2068,
   Johannes U. Lange^2,3,40000-0002-2450-1366




   ^1Department of Astronomy, Yale University, PO. Box 208101, New Haven, CT 06520-8101

   ^2Department of Physics, American University, 4400 Massachusetts Avenue NW, Washington, DC 20016, USA

   ^3Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA
 
   ^4Leinweber Center for Theoretical Physics, University of Michigan, Ann Arbor, MI 48109, USA







firstpage–lastpage
2024

[
    [
    
=====







§ ABSTRACT 
is a novel Bayesian hierarchical method for inferring the galaxy-halo connection, including its scatter, using the kinematics of satellite galaxies extracted from a redshift survey. In this paper, we introduce crucial improvements, such as updated central and satellite selection, advanced modelling of impurities and interlopers, extending the kinematic modelling to fourth order by including the kurtosis of the line-of-sight velocity distribution, and utilizing satellite abundance as additional constraint. This drastically enhances 's performance, resulting in an unbiased recovery of the full conditional luminosity function (central and satellite) and with unprecedented precision. After validating ’s performance using realistic mock data, we apply it to the SDSS-DR7 data. The resulting inferences on the galaxy-halo connection are consistent with, but significantly tighter than, previous constraints from galaxy group catalogues, galaxy clustering and galaxy-galaxy lensing. Using full projected phase-space information, breaks the mass-anisotropy degeneracy, thus providing precise global constraint on the average orbital velocity anisotropy of satellite galaxies across a wide range of halo masses. Satellite orbits are found to be mildly radially anisotropic, in good agreement with the mean anisotropy for subhaloes in dark matter-only simulations. Thus, we establish as a powerful tool that is not only more constraining than other methods on similar volumes of data, but crucially, is also insensitive to halo assembly bias which plagues the commonly used techniques like galaxy clustering and galaxy-galaxy lensing. 
 




methods: analytical —
methods: statistical —
galaxies: haloes — 
galaxies: kinematics and dynamics —
cosmology: dark matter







§ INTRODUCTION



According to the current cosmological paradigm, the vast majority of all galaxies form and reside in extended dark matter haloes <cit.>. Halo occupation modelling tries to use observational constraints on the population of galaxies in order to infer the statistical link between the galaxy properties (mainly their luminosity or stellar mass) and the properties of the dark matter haloes (mainly some measure of halo mass) in which they reside <cit.>.  The resulting `galaxy-halo connection' provides valuable insight regarding the formation and evolution of galaxies, and benchmarks to calibrate, compare and validate semi-analytic models <cit.> and simulations <cit.>. In addition, since it describes the link between the light we see and the mass that governs the dynamical evolution of the Universe, it is a powerful tool that allows astronomers to constrain cosmological parameters using the observed distribution of galaxies <cit.>.

Arguably, the most straightforward method to infer the galaxy halo connection, and one that has become extremely popular, is subhalo abundance matching <cit.>. It matches the ordered list of galaxies (typically ranked by stellar mass or luminosity) to that of subhaloes (typically ranked by their peak or infall mass)[Rather than abundance matching individual galaxies to subhaloes, one can also match the abundance of galaxy groups (identified using some group finder) to dark matter host haloes <cit.>.]. In this mostly non-parametric method, one usually allows for some amount of scatter (a free parameter) in the rank-order matching, to have realistic spread in the stellar mass - halo mass relation <cit.>. A key advantage of SHAM over other methods, discussed below, is that it only requires stellar mass (or luminosity) measurements of the galaxies. However, an important downside is that it relies crucially on the assumption that both the galaxy sample and the (sub)halo sample (typically taken from a N-body simulation) are complete. Hence, SHAM cannot be applied to subsamples of galaxies (i.e., samples of emission line galaxies, or galaxies selected by colour). In addition, even if the galaxy sample is complete, the (sub)halo catalogues used, which are typically extracted from numerical simulations, suffer from incompleteness due to artificial disruption <cit.> and failures of subhalo finders <cit.>. This can significantly impact the galaxy-halo connection inferred via SHAM <cit.>. 

These problems can be overcome using data that more directly constrains halo mass. The two most commonly used methods are galaxy clustering <cit.> and galaxy-galaxy lensing <cit.>. The former relies on the fact that more massive haloes are more strongly clustered <cit.>; hence, the clustering strength of a given population of galaxies informs the characteristic mass of the haloes in which they reside. Unfortunately, its reliability is hampered by the finding that halo clustering strength depends not only on mass but also on secondary halo properties <cit.>, something that is collectively referred to as halo assembly bias. Galaxy-galaxy lensing, which is a manifestation of weak gravitational lensing, uses the tangential shear distortions of distant background galaxies around foreground ones in order to constrain the halo masses of the latter <cit.>. Although, in principle, a fairly direct probe of halo mass, this method requires tedious shape measurements of faint background sources, which can be prone to effects like blending and intrinsic alignment. Typically the background sources lack spectroscopic redshifts, which can also cause systematic errors in the modelling of their measured shear distortions. In addition, on large scales the 2-halo term of the lensing shear is also impacted by the same assembly bias issues that plague clustering.

Another method that can be used to constrain the galaxy-halo connection, but which has hitherto been severely under-utilized, is satellite kinematics. It uses measurements of the line-of-sight velocities of satellite galaxies with respect to their corresponding central galaxy in order to constrain the gravitational potential, and hence the mass, of the host halo[Famously, the same principle was used by <cit.> in order to infer the presence of dark matter in the Coma cluster]. With the exception of large galaxy groups and clusters, individual central galaxies typically only have a few spectroscopically detected satellites. Consequently, it is common to combine the satellite velocity measurements from a large stack of central galaxies in order to estimate an average satellite velocity dispersion, which in turn is used to infer an average host halo mass using either a virial mass estimator or a simple Jeans model <cit.>.

It has often been argued that satellite kinematics is not a reliable mass estimator for any combination of the following reasons: (a) satellite galaxies are not necessarily virialized tracers <cit.>, (b) their orbits may well be anisotropic <cit.>, resulting in a well-known mass-anisotropy degeneracy <cit.>, and (c) the stacking that is used implies `mass-mixing' (i.e., combining the kinematics of haloes of different masses), which muddles the interpretation of the data. In addition, the selection of centrals and satellites from a redshift survey is unavoidably impacted by impurities, incompleteness and interlopers, further complicating the analysis. Despite these concerns, a number of studies have progressively improved satellite kinematics and have shown that it can yield reliable, as well as precise, constraints on the galaxy-halo connections[and on the masses of individual clusters <cit.>]. In particular, <cit.> demonstrated that by selecting centrals and satellites using iterative, adaptive selection criteria the impact of impurities and interlopers can be minimized.  <cit.> has shown that by combining different weighting schemes one can accurately account for mass mixing, and even constrain the scatter in the stellar mass-halo mass relation <cit.>. This was significantly improved upon by <cit.> who demonstrated that kinematics of satellite galaxies from a large redshift survey such as the Sloan Digital Sky Survey <cit.> can yield constraints on the galaxy-halo connection that are complementary to, and competitive with, constraints from galaxy clustering and/or galaxy-galaxy lensing. 
<cit.> were the first to analyse satellite kinematics while accounting for orbital anisotropy. Using a method first developed by <cit.> they were able to simultaneously constrain halo mass, halo concentration and orbital anisotropy, albeit without accounting for mass mixing.

<cit.> developed , a Bayesian hierarchical inference formalism that further improves on the ability of satellite kinematics to constrain the galaxy halo connection. Unlike previous methods, does not resort to stacking the kinematics of satellite galaxies in bins of central luminosity, and does not make use of summary statistics, such as satellite velocity dispersion. Rather, it leaves the data in the raw form and computes the corresponding likelihood. Consequently, it can simultaneously solve for halo mass and orbital anisotropy of the satellite galaxies, while properly accounting for scatter in the galaxy-halo connection. In addition, can be applied to flux-limited, rather than volume-limited samples, greatly enhancing the amount and dynamic range of the data. 

Paper I also tested and validated against mock data sets of varying complexity, and demonstrated that it yields unbiased constraints on the parameters specifying the galaxy-halo connections. However, in order to speed up the analyses, all those tests where performed using mock data samples that were only about 1/8 the size of the full SDSS sample analysed here.  When we ran on full-sized mocks instead, the model parameter uncertainties shrank considerably, as expected, revealing several significant discrepancies that turned out to be systematic. This necessitated a number of modifications to , which we present in the first half of this paper. Most notably, we introduce significant improvements to the treatment of interlopers (i.e., galaxies that are selected as satellites but that do not reside in the same dark matter halo as the central), allowing for both a population of splash-back galaxies <cit.> and a large scale infall population akin to the well-known <cit.> effect. In addition, we slightly modify the cylindrical selection criteria in order to improve the purity of our sample (i.e., reduce misclassification of satellites as centrals), we assure that the selection of secondaries around each individual primary is volume-limited, and we forward-model the contribution of impurities that arise from haloes in which the brightest galaxy is a satellite rather than the central. We also let go of the oversimplified assumption that the satellite velocity profile along any given line-of-sight is Gaussian, as was done in . Rather, we now use the fourth-order Jeans equations to model  the kurtosis of the line-of-sight velocity distribution (LOSVD). This enables more accurate modelling of the full phase-space distribution of satellite galaxies, and allows to break the mass-anisotropy degeneracy. Finally, we also replace fitting binned statistics of centrals with zero (detected) satellites, as done in , with a more general, Bayesian hierarchical modelling of the number of satellites around each central. Although  this data on satellite abundances does not yield direct kinematic constraints on halo mass, it greatly helps to constrain the overall galaxy-halo connection. 

The goal of this paper is threefold: (i) showcase the advancements in satellite kinematics methodology that we have introduced in , and highlight its improved performance when tested against realistic SDSS-like mock data; (ii) apply to SDSS DR7 data to simultaneously constrain the conditional luminosity functions of central and satellite galaxies, the satellite velocity anisotropy and satellite radial distribution, all with unprecedented precision, and compare those with previous constraints on halo occupation statistics; and (iii) establish as a powerful method to infer the galaxy-halo connection which is free of halo assembly bias effects, and that is even more constraining than commonly used techniques like galaxy clustering and galaxy-galaxy lensing, when applied on data of similar volumes.

Throughout this paper we adopt the flat Planck18 ΛCDM cosmology <cit.>, which has matter density parameter Ω_ = 0.315, power spectrum normalization σ_8=0.811, spectral index n_ = 0.9649, Hubble parameter h = (H_0/100) = 0.6736 and baryon density Ω_ h^2 = 0.02237.[These are the TT,TE,EE+lowE+lensing best-fit values assuming a base-ΛCDM cosmology.]




§ SAMPLE SELECTION





 §.§ Selecting central-satellite pairs



The first step in analysing satellite kinematics is to select a sample of centrals and their associated satellites from a redshift survey. Unfortunately, this selection is never perfect; one undoubtedly ends up selecting some bright satellites as centrals (we refer to these as `impurities') and not every galaxy selected as a satellite actually resides in the same host dark matter halo as the corresponding central (those that don't are referred to as `interlopers'). In what follows, we therefore use `primaries' and `secondaries' to refer to galaxies that are selected as centrals and satellites, respectively.

A galaxy at redshift z is considered a potential primary if it is the brightest galaxy in a conical volume of opening angle Θ_ ap^ pri≡/d_(z)  centered on the galaxy in question, and extending along the line-of-sight from z-(Δ z)^ pri to z+(Δ z)^ pri. Here d_(z) is the angular diameter distance at redshift z, and (Δ z)^ pri = (/c)   (1+z). The parameters  and  specify the primary selection cone. Following <cit.>, we select the primaries in a rank-ordered fashion, starting with the most luminous galaxy in the survey. Any galaxy located inside the selection cone of a brighter galaxy is removed from the list of potential primaries. All galaxies fainter than the primary and located inside a similar cone, but defined by  and , centred on the primary are identified as its secondaries.  Note that, although it is common to refer to these selection volumes as `cylinders', a convention we also adopted in , in actuality the selection volumes are frustums of cones. In order to rectify this confusing nomenclature, in this paper we refer to them as `selection cones' (see Fig. <ref>).

The four parameters , , , and  control the completeness and purity of the sample of primaries and secondaries. Increasing  and/or , boosts the purity among primaries (i.e., it reduces the number of satellites erroneously identified as centrals), but reduces the overall completeness. Similarly, decreasing  and/or  reduces the number of interlopers, but at the costs of a reduced number of satellites, which are the dynamical tracers of interest. Since brighter primaries typically reside in larger, more massive haloes, it is advantageous to scale the sizes of the selection cones with the luminosity of the primary <cit.>. In particular, we adopt = 0.6  σ_200, = 0.15  σ_200, and =  = 1000 σ_200. Here σ_200 is a rough measure of the satellite velocity dispersion in units of 200, which, following <cit.> and <cit.>, we take to scale with the luminosity of the primary as


    logσ_200 = 0.04 + 0.48 log L_10 + 0.05 (log L_10)^2 ,


where L_10 = L / (10^10), and σ_200 is allowed to take a maximum value of 4. The values of  and  correspond to roughly 1.65 and 0.4 times the halo virial radius, respectively, while the value for  is large enough to include the vast majority of all satellites around primaries of the corresponding luminosity. Note that the numerical values in these selection criteria are tuned in order to optimize the selection of primaries and secondaries against impurities and interlopers. In particular, they are slightly different from the values we adopted in . As detailed in Appendix <ref>, this is done in order to reduce the fraction of impurities that are neither true centrals, nor the brightest satellites in their corresponding host haloes. These impurities are particularly difficult to account for in our forward-modelling approach and can cause a small but systematic overestimate of the scatter in the relation between halo mass and central luminosity.



The SDSS redshift catalogue, to which we apply in this study, is a flux-limited survey. As emphasized in , an important advantage of over earlier studies of satellite kinematics, is that it is not limited to volume-limited subsamples, thereby greatly boosting the number of primaries and secondaries to be used in the analysis. However, in order to facilitate proper modelling of the number of secondaries (true satellites and interlopers) we need to assure that the selection of secondaries around each individual primary is volume limited. This is something that was not implemented in , but which turned out to be important in order to avoid a systematic bias in the inferred faint-end slope of the satellite luminosity function. This effect was not significant in the smaller mock data samples used to test in , but could no longer be overlooked using data sets comparable in size to the SDSS data used here.

In this paper, we limit our analysis to primaries in the luminosity range 9.504 ≤log(/[h^-2]) ≤ 11.104, corresponding to -19 ≤ M_r^0.1 - 5log h ≤ -23. Here M_r^0.1 is the absolute magnitude in the SDSS r-band K+E corrected to z=0.1. In addition, we only use data in the redshift range 0.02≤ z ≤ 0.20.  Note that the selection cone used to identify secondaries around a primary of luminosity L_ pri at redshift z_ pri extends from  to , given by


    = z_ pri -  Δ V_ max^ sec(L_ pri) / c (1 + z_ pri)      and
        = z_ pri +  Δ V_ max^ sec(L_ pri) / c (1 + z_ pri)  .


Hence, as depicted in Fig. <ref> we are only complete in the selection of secondaries with luminosities ≥ L_ min(). Here L_ min(z) is the minimum luminosity of galaxies at redshift z that make the apparent magnitude limit of our survey data (m_r = 17.6; see <ref>). In order to assure a complete, volume-limited selection of secondaries around each primary, secondaries fainter than L_ min() are discarded. In addition, in order to assure that the entire secondary selection cone around a given primary fits within the flux-limits of the SDSS data, we require that  > L_ min(). Finally, the redshifts of the primaries are restricted to 0.034 ≤ z_ pri≤ 0.184, such that the  and  of the most luminous primaries fit within the 0.02 ≤ z ≤ 0.20 limits of the entire sample. 

To elucidate this further, Fig. <ref> illustrates the bounds on the selection of primaries and secondaries.  The solid, vertical lines at z=0.02 and z=0.20 mark the minimum and maximum redshifts of the entire sample, while the dashed, vertical lines at z=0.034 and z=0.184 mark the redshift limits allowed for primaries. Dashed, horizontal lines at log L = 11.104 and log L = 9.504 mark the luminosity cuts for primaries. The solid circles, labelled A to H, represent hypothetical primaries of three different luminosities L_1, L_2 and L_3 (indicated by three different colours). The shaded rectangle associated with each primary indicates the allowed luminosity-redshift ranges of its corresponding secondaries. Since the redshift extent of the secondary selection cone scales with the luminosity of the primary, the shaded regions of fainter primaries have a smaller Δ z-extent, as is evident from the figure. Note that these shaded regions extend down to where the apparent magnitude of the secondaries at the back-end of the selection cone is equal to the magnitude limit of the survey, which can be  significantly lower than log L = 9.504, specifically for the primaries that are relatively nearby (e.g., primaries A and F).

Primary A is at the minimum allowed redshift for primaries, z_ min^ pri = 0.034, such that the `front' end of its secondary selection cone is equal to the minimum redshift of our survey data (i.e., =0.02). Similarly, primary C is located at the maximum allowed redshift for primaries, z_ max^ pri = 0.184, and has = 0.20.  Primary H is also special in that it has the highest redshift possible given its luminosity. Had it been any farther away, the far end of its secondary selection cone would stick outside of the SDSS flux limit, resulting in incompleteness. The three dashed and dotted curved lines, labelled as L_ min((L_1/2/3)), show the lower luminosity limits for secondaries as a function of , corresponding to each of the three different primary luminosities represented in the figure. For example, primaries of luminosity L_2 (like primary D and E) can not have secondaries fainter than the middle dashed curve in green. This ensures that their secondaries are individually volume limited around each of those primaries.






 §.§ Survey incompleteness



As any spectroscopic redshift survey, the SDSS data, from which our sample of primaries and secondaries derives, suffers from spectroscopic incompleteness due to fibre collisions and other failure modes <cit.>. Each galaxy in the survey is assigned a spectroscopic completeness, 𝒞_ spec, which indicates the fraction of spectroscopic targets in the angular region of the galaxy in question with a useful spectrum. In order to avoid primaries in regions with poor spectroscopic completeness, we remove all primaries with 𝒞_ spec < 0.8.

If a primary is close to the edge of the survey, such that its secondary selection cone sticks partially outside of the survey footprint, or if the secondary selection cone overlaps with a masked region, the number of secondaries may be underestimated. In order to account for this, we proceed as follows. For each primary we uniformly distribute ∼ 5× 10^4 particles in its secondary selection cone. We then compute the fraction, w_ app, of those particles that are located inside the angular footprint of the SDSS, accounting for both survey edges and masked areas. In what follows we use w_ app,i to denote the aperture completeness of galaxy i. In order to avoid primaries with a poor aperture completeness, we remove all primaries with w_ app < 0.8.

As demonstrated in <cit.>, it is important to correct satellite kinematics data for  fibre-collision induced incompleteness. In the SDSS, spectroscopic fibres cannot be placed simultaneously on a single plate for objects separated by less than ϑ_ fc≡ 55” <cit.>. Although some galaxies are observed with multiple plates, yielding spectroscopic redshifts even for close pairs, roughly 65% of galaxies with a neighbour within 55” lack redshifts due to this fibre collision effect. In order to correct the data for the presence of fibre collisions, we follow <cit.> and start by assigning each fibre-collided galaxy the redshift of its nearest neighbour <cit.>. Note that we only use these during the identification of primaries. Once the selection is complete, all fiber-collided primaries and secondaries are removed from the sample[As shown in <cit.>, including fibre-collided galaxies during the selection of primaries significantly reduces sample impurity.]. In addition, each galaxy is assigned a spectroscopic weight, w_ spec, that is computed as follows. For each galaxy we first count the number of neighboring galaxies, n, brighter than m_r = 17.6 within a projected separation less than 55”. Next, for all galaxies in the survey with n neighbours, we compute the fraction, f_ spec, of those neighbours that have been successfully assigned a redshift. Finally, all galaxies with n neighbours are then assigned a spectroscopic weight equal to w_ spec = 1 / f_ spec. 

In order to correct for aperture incompleteness and fibre collisions, down-weights the expectation value for the number of secondaries around primary i (see equation [<ref>] below), using the following correction factor:


    f_ corr,i = /∑_j=1^ w_ spec,ij  w_ app,i .


Here w_ spec,ij is the spectroscopic weight for secondary j associated with primary i. Since correcting for fibre collisions is extremely difficult on scales below the fibre-collision scale, we remove all secondaries with R_ < R_ cut() ≡ d_()  ϑ_ fc. Hence, the secondary selection volumes used in the end are conical frustums with a central hole with an opening angle of 55” (see Fig. <ref>). As shown in and <cit.>, this combined approach of down-weighting the model predictions for the number of secondaries and ignoring secondaries below the fibre-collision scale accurately accounts for incompleteness arising from fibre-collisions in the SDSS.



§ OBSERVABLES



Here we describe the various observables used by in order to constrain the galaxy-halo connection. These include (i) accessible 2D phase-space parameters of primary-secondary pairs (line-of-sight velocity and projected separation), which contains the information regarding the kinematics of satellite galaxies, (ii) statistics regarding the number of secondaries per primary (including primaries with zero secondaries), which helps to constrain the halo occupation statistics, and (iii) the galaxy luminosity function. The following subsections discuss each of these observables in detail.



 §.§ Satellite kinematics



For each primary-secondary pair in the sample we compute their projected
separation 


    = d_()  ϑ ,


and their line-of-sight velocity difference


    = c  ( - )/1 +  .


Here  and  are the observed redshifts of the primary and secondary, respectively, c is the speed of light, and ϑ is the angular separation between the primary and secondary on the sky.

As detailed in , the main data vector used in is given by 


    _ SK = ⋃_i=1^N_+ ( {,  | j=1,...,} | , , ) .


where the union is over all N_+ primaries with at least one secondary. Here  is the number of secondaries associated with primary i, and it is made  explicit that , , and  are only treated as conditionals for the data {,  | j=1,...,}. In other words, we consider ,  and  as `given' and shall not use the distributions of these quantities as constraints on our satellite kinematics likelihood. Rather, uses the number densities of all galaxies as additional constraints (see <ref>). The main reason for doing so is to make the method less sensitive to the detailed selection of primaries, which is difficult to model in detail. In particular, as discussed in , this approach makes fairly insensitive to details regarding the σ_200(L) relation (equation [<ref>]) used to define the selection cones.



 §.§ Number of secondaries



The data vector _ SK described above only contains primaries with at least one secondary. The complementary data vector _0 = ({,}  |  i=1,2,...,N_0) lists all N_0 primaries with zero spectroscopically detected secondaries. Even though _0 contains no kinematic data, it still provides additional constraints on the galaxy-halo connection, in particular regarding the occupation statistics of satellite galaxies. In we utilized this information by computing the fraction, P_0 = N_0 /(N_0 + N_+), of primaries, in a given bin of log and , that have zero secondaries. Here N_0 is the number of `lonely primaries' with zero detected secondaries, and N_+ is the number of primaries that have at least one secondary. As discussed in , this P_0 statistic provides valuable constraints on the galaxy-halo connection. However, upon closer examination we found that the binning used in this method causes small, but systematic errors in the inference. Using smaller bins was not able to solve this problem, which is why we ultimately opted for the following alternative, unbinned approach. 

In line with 's philosophy to leave the data as much as possible in its raw form, rather than computing P_0 on a (log, )-grid, we use the following raw data vector as constraint on the model:


    = ⋃_i=1^N_ NS(  | , ) .


Here the union is over (a random subset of) all = N_0 + N_+ primaries, independent of how many secondaries they have (i.e., including the primaries with zero secondaries). Since N_0 ≫ N_+, computing the likelihood  for all  primaries is much more time-consuming than computing the likelihood for the satellite kinematics data vector (eq. [<ref>]). Therefore, we only use a downsampled, random subset of N_ NS = (N_+) primaries, where each primary has a probability equal to N_+/ to be included. In the case of the SDSS data set described in <ref> this probability is 0.094. This downsampling assures that the computation of the likelihood for  has a CPU requirement that is comparable to that for _ SK. We emphasize that our constraints are primarily driven by the satellite kinematics data. Hence, this down-sampling of the satellite abundance data has no significant impact on our constraining power of the central galaxy-halo connection or the orbital anisotropy of the satellite galaxies. It only slightly broadens the posterior constraints for some of the parameters characterizing  
galaxy-halo connection of satellite galaxies.



 §.§ Galaxy Number Densities



The final observable that we use to constrain the galaxy-halo connection is the galaxy luminosity function, which provides important additional constraints on the CLF <cit.>, and therefore helps to tighten the posterior in our inference problem. We use the number density of galaxies in 10 bins of 0.15 dex in luminosity, ranging from 10^9.5 to 10^11. These are computed using the corresponding, volume-limited subsamples, carefully accounting for the SDSS DR-7 footprint. In what follows, we refer to the data vector representing these 10 number densities as _ LF. The covariance matrix of this data is computed using a jackknife estimator. In particular, we apply a recursive routine[ https://github.com/rongpu/pixel_partition] developed by <cit.>, that takes into account the survey mask and window, and iteratively constructs 𝒩 maximally compact, equal-area partitions of the survey footprint <cit.>. We adopt 𝒩=100 which is large enough to capture the covariance in the survey while also being small enough to assure that each subregion still hosts an adequate number of galaxies.[We apply a Hartlap correction factor <cit.> to the inverse of the covariance matrix to account for the relatively small number of jackknife samples, but note that this has a negligible impact.]



§ METHODOLOGY



We analyse the data described above using the Bayesian, hierarchical satellite kinematics code , which is described in detail in . Here we briefly summarize its salient features and introduce a few modifications that improve 's performance.

uses an affine invariant ensemble sampler <cit.> to constrain the posterior distribution, 


    P( | ) ∝( | )   P() .


Here = _ SK +  + _ LF is the total data vector,  is the vector that describes our model parameters, P() is the prior probability distribution on the model parameters, and (|) is the likelihood of the data given the model. The latter consists of three parts: the likelihood _ SK for the satellite kinematics data _ SK, the likelihood _ NS for the numbers of secondaries as described by the data vector , and the likelihood _ LF for the luminosity function data _ LF. In what follows we briefly describe the computation of each of these three different likelihood terms in turn. However, we first describe the model that we use to characterize the galaxy-halo connection.



 §.§ Galaxy-halo connection model





  §.§.§ Conditional luminosity function



The galaxy occupation statistics of dark matter haloes are modelled using the conditional luminosity function (CLF), Φ(L|M,z)   L, which specifies the average number of galaxies with luminosities in the range [L- L/2,   L +  L/2] that reside in a halo of mass M at redshift z <cit.>. In particular, we write that


    Φ(L|M,z) = Φ_(L|M) + Φ_(L|M) .


Here and throughout the rest of the paper, subscripts `c' and `s' refer to central and satellite, respectively, and we assume that the CLF is redshift independent, at least over the redshift range considered here (0.02 ≤ z ≤ 0.20).




The CLF of centrals is parametrized using a log-normal distribution (see blue, dashed curve in Fig. <ref>),


    Φ_ (L | M)  L = log e/√(2πσ_^2)exp[ -(log L - logL̅_/√(2)σ_)^2]  L/L.


The mass dependence of the median luminosity, L̅_, is parametrized by a broken power-law: 


    L̅_ (M) = L_0 (M / M_1)^γ_1/(1 + M / M_1)^γ_1 - γ_2.


which is characterized by three free parameters; a normalization, L_0, a characteristic halo mass, M_1, and two power-law slopes, γ_1 and γ_2. 

Motivated by the fact that several studies suggest that the scatter, σ_, increases with decreasing halo mass <cit.>, we allow for a mass-dependent scatter using


    σ_(M) = σ_13 + σ_ (log-13)


Hence, the scatter is characterized by two free parameters, a normalization, σ_13, that specifies the intrinsic scatter in log in haloes of mass = 10^13, and a power-law slope σ_. Note that this is slightly different from the parametrization adopted in .

For the satellite CLF we adopt a modified Schechter function (see red curve in  Fig. <ref>): 


    Φ_ (L | M) = ϕ_^*/L_^*( L/L_^*)^α_exp[ - ( L/L_^*)^2 ].


Thus, the luminosity function of satellites in haloes of a given mass follows a power-law with slope α_ and with an exponential cut-off above a critical luminosity, L_^*(M). Throughout we adopt


    L_^*(M) = 0.562  L̅_ (M) .


which is motivated by the results from galaxy group catalogues <cit.>. As in <cit.>, we assume a universal value for the faint-end slope of the satellite CLF, α_ s, independent of halo mass. Finally, the normalization ϕ_^*(M) is parametrized by 


    log[ ϕ_^*(M) ] = b_0 + b_1 log M_12 + b_2 (log M_12)^2.


where M_12 = M/(10^12). Note that this characterization of the CLF is very similar to that adopted in a number of previous studies <cit.>. All CLF parameters, along with parameters that characterize the satellite velocity anisotropy, and nuisance parameters used for interloper modelling, are listed in Table <ref>. It also includes the best-fitting values and 1σ confidence intervals for all the parameters, obtained by fitting the SDSS-DR7 data.



  §.§.§ Spatial distribution of satellites



Throughout we assume that the radial distribution of satellite galaxies is given by a spherically symmetric, generalized NFW (gNFW) profile


    n_ sat(r|M,z) ∝( r/  r_)^-γ( 1 + r/  r_)^γ - 3 .


Here  and γ are free parameters and r_ is the scale radius of the dark matter halo, which is related to the halo virial radius via the concentration parameter =  / r_. This gNFW profile has sufficient flexibility to adequately describe a wide range of radial profiles, from satellites being unbiased tracers of their dark matter halo (γ =  = 1), to cored profiles that resemble the radial profile of surviving subhaloes in numerical simulations (γ=0, ∼ 2). This also brackets the range of observational constraints on the radial distribution of satellite galaxies in groups and clusters <cit.>.



 §.§ Satellite Kinematics



The data vector for the satellite kinematics is given by eq. (<ref>) and contains the projected phase-space coordinates  and  of all secondaries (satellite galaxies plus interlopers) associated with the N_+ primaries (centrals plus impurities). We make the reasonable assumption that the data for different primaries is independent. Additionally, for a primary with more than one secondary, we assume that the phase-space distribution of the secondaries are not correlated with each other. The latter may not be entirely justified, given that satellites are often accreted in groups, which can bias halo mass estimates <cit.>. We emphasize, though, that the majority (71% in the case of the SDSS data discussed in <ref>) of primaries that contribute to the satellite kinematics data only have a single secondary. In addition, tests based on realistic simulation-based mocks (see <ref>) indicate that any potential correlations between satellites (subhaloes) that occupy the same host halo can safely be ignored (i.e., do not cause a significant systematic error in our inference). Hence, we have that


    _ SK   ≡(_ SK|) 
        = ∏_i=1^N_+ ∏_j=1^ P(,  | , , , ) .


Here, P(,  | , , ) is the probability that a secondary galaxy in a halo at redshift , with a primary of luminosity, , and with a total of  detected secondaries has projected phase-space parameters (, ). For true satellites, the probability is computed assuming that satellite galaxies are a virialized, steady-state tracer of the gravitational potential well in which they orbit (see <ref>). Throughout, we assume dark matter haloes to be spherical and to have NFW <cit.> density profiles characterized by the concentration-mass relation of <cit.> with zero scatter. Hence, host haloes are completely specified by their virial mass, , alone[Throughout this paper, we define virial quantities according to the virial overdensities given by the fitting formula of <cit.>.], which implies that we can factor the likelihood as


    _ SK    = ∏_i=1^N_+ ∫  P( | , , )  ×
                              ∏_j=1^ P(, |, , ) .


This equation describes a marginalization over halo mass, which serves as a latent variable for each individual primary, accentuating the hierarchical nature of our inference procedure. Note that the `prior' for halo mass is informed by , , and  according to the model .  Using Bayes theorem, we have 


    P(M|L,z,) = P(|M,L,z)   P(M,L,z)/∫ M   P(|M,L,z)   P(M,L,z) .


In what follows we discuss each of the conditional probability functions required to compute _ SK in turn.



  §.§.§ The probability P(|M,,)



The number of secondaries, N_, associated with a particular primary consists of both satellites (galaxies that belong to the same dark matter host halo as the primary), and interlopers (those that do not). Throughout we assume that the number of interlopers and the number of satellite galaxies are independent, and that both obey Poisson statistics. As shown in , this implies that


    P(N_|M,L,z) = λ^N__ tot ^-λ_ tot/N_! ,


where


    λ_ tot = f_ corr× [λ_ sat + λ_ int] ,


is the expectation value for the number of secondaries, corrected for fibre collision and aperture incompleteness using the correction factor of equation (<ref>), and with λ_ int(L,z) and λ_ sat(M,L,z) as the expectation values for the numbers of interlopers and satellites, respectively.

The expectation value for the number of satellites brighter than the magnitude limit L_ min(), in a halo of mass M at redshift , that fall within the aperture used to select secondaries around a primary of luminosity , is given by


    λ_ sat(M,,) = f_ ap(M,,)  ∫_L_ min^∞Φ_(L|M)   L


Note that L_ min is a function of  which in turn is a function of  and  (see <ref>). Here Φ_(L|M) is the satellite component of the CLF given by equation (<ref>) and f_ ap is the aperture fraction, defined as the probability for true satellites to fall within the secondary selection cylinder specified by  and . Given that  is much larger than the extent of the halo in redshift space, we have that


    f_ ap(M   , , ) = 
        4 π∫_0^n̅_ sat(r|M,) 
     [ζ(r,) - ζ(r,)]   r^2   r .


Here = (M,) is the virial radius of the halo in question, and = () and = R_ cut() are the outer and inner radii of the conical volume used to select the secondaries. The function n̅_ sat(r|M,z) is the average, radial profile of satellites around haloes of mass M at redshift z, normalized such that


    4 π∫_0^n̅_ sat(r|M,z)   r^2   r = 1 ,


and


    ζ(r, R) = 
    1      if  r ≤ R 
    
    1 - √(1 - R^2 / r^2)     otherwise.


More specific expressions for λ_ sat(M,L,z) and f_ ap(M,L,z) are provided in .

For the interlopers, one naively expects the abundance to be proportional to the number density of galaxies within the relevant range of luminosities and the volume of the secondary selection cone. However, being biased tracers of the mass distribution, galaxies are highly clustered which typically will boost the number density of galaxies in the vicinity of a bright primary. Moreover, this clustering strength is known to depend on halo mass, galaxy luminosity and redshift <cit.>, and to be affected by peculiar velocities, in particular due to large-scale infall <cit.>. We bypass the intricate complexities involved with modeling this clustering on small scales by modeling the expectation value for the number of interlopers as the product of an effective `bias', b_ eff, and the expectation value for the number of galaxies with L_ min() < L < in a randomly located conical selection volume, V_ cone(, ):


    λ_ int(,) = b_ eff× V_ cone×n̅_ gal .


where each term on the right-hand side is a function of {,}. Here


    n̅_ gal(, ) = ∫_L_ min^ L ∫_0^∞Φ(L|M)   n(M,)   M


is the average number density of galaxies at redshift  with luminosity in the range [L_ min,],
with  n(M,z) the halo mass function at redshift z, computed using the fitting function of <cit.>, and 


    V_ cone(,) = π [^2 - ^2]  2 /H()  (1+)^3


with H(z) the Hubble parameter[Note, there was a typo in Eq. (22) in , where the power-law index of (1+) was 2, rather than the correct 3]. The effective bias is modelled as


    b_ eff(,) = η_0  (/10^10.5)^η_1 ( 1+)^η_2


where η_0, η_1, and η_2 are three free nuisance parameters that fully specify our interloper bias model, and that are constrained simultaneously with all other physical parameters. This model has proved to be sufficiently flexible to accurately model the full complexity of interloper abundance in realistic simulation-based mock data (see Section <ref>).



  §.§.§ The probability P(M,,)



The function P(M,,) describes the probability distribution function of primaries as a function of host halo mass, luminosity and redshift, and can be written as


    P(M,,) = P(|M,)    n(M,)  (M,,)  .


where n(M,z) is the halo mass function <cit.> and (M,L,z) is a `completeness', to be defined below. As in , if we assume that all primaries are true centrals, then we have that P(|M,) = Φ_(|M). However, in reality some primaries are misidentified satellites, and such impurities need to be accounted for. In we argued that the impact of these impurities is sufficiently small that it can be ignored. Although this was indeed the case for the small mock data sets used there, the impact of impurities can no longer be ignored when using data sets similar in size to the SDSS data analysed here. In fact, detailed tests showed that they can systematically bias the inferred scatter in the relation between halo mass and central luminosity, and we therefore devised the following scheme in order to account for impurities.

The vast majority of all impurities in realistic SDSS-like mocks (such as the Tier-3 mock described in <ref>) are those satellite galaxies which happen to be the brightest galaxy in their halo (even brighter than their central). In what follows, we refer to these as Type-I impurities. Since primaries are by definition the brightest galaxies in their selection cones, such brightest-halo-galaxy (hereafter BHG) satellites typically end up being selected as primaries, rather than their corresponding central. In rare cases, a primary is neither a true central, nor a BHG satellite. We refer to these as Type-II impurities, which arise, for example, if the true central is the BHG but is absent from the SDSS survey data, either because of fiber collisions or because it falls  outside the window of the SDSS footprint. As detailed below, Type I impurities can be accounted for in our new theoretical modelling. However, Type-II impurities are virtually impossible to model accurately. Using detailed mock data sets, we therefore tuned our selection criteria in order to minimize the contribution of Type-II impurities. In particular, we found that we were able to significantly reduce their frequency by slightly enlarging the volume of the primary selection cone as described in <ref>; In particular, the new criteria reduce the fraction of Type-II impurities from ≳ 1% when using the old selection criteria used in , to ∼ 0.5% with our new selection criteria. More importantly, in mock data, the new selection criteria predominantly eliminate the presence of Type-II impurities that are extreme outliers of the average relation between halo mass and primary luminosity, and which are the main culprits for causing mild systematic errors in the inferred galaxy-halo connection (specifically in the scatter, σ_). Detailed tests with mock data, presented in <ref> below, show that our new primary selection criteria sufficiently suppress the impact of Type-II impurities that it allows for unbiased estimates of the galaxy-halo connection (at least for a survey the size of SDSS).

Therefore, in what follows, we ignore Type-II impurities and assume that primaries are either true centrals or BHG satellites (i.e., Type-I impurities). Hence, we have that

 
    P(|M,z)     =  P(=|M,z) P(L_ bs<|M,z)
     
       + P(<|M,z) P(L_ bs=|M,z).


Here P(L_ bs<L|M,z) is the probability that the brightest satellite in a halo of mass M at redshift z has a luminosity less than L, which is given by


    P(L_ bs<L|M,z) = exp[ - Λ(L|M,z) ] .


Here Λ(L|M,z) is the expectation value for the number of satellites brighter than L in a halo of that mass and redshift, which in turn is given by


    Λ(L|M,z) = ∫_L^∞ L'  Φ_(L'|M)


Differentiating  P(L_ bs<L|M,z) with respect to luminosity yields:


    P(L_ bs=L|M,z) = Φ_ s(L|M,z)   e^-Λ(L|M,z)


The two other terms that appear in equation (<ref>) are  P(=|M,z), which is simply equal to the central CLF, Φ_(|M), and its cumulative distribution, which is given by


    P(<L|M,z) = 
    ∫_0^L'  Φ_('|M)



The expression for P(|M,) given by equation (<ref>), when substituted in equation (<ref>), accurately forward models the impact of the vast majority of impurities. 

The final ingredient we need is an expression for the completeness (M,L,z), which is defined as the fraction of haloes of mass M at redshift z= with a central or brightest satellite of luminosity  that falls within the survey volume of the SDSS, and that is selected as a primary by our selection criteria.  In general we have that (M,L,z) = (M|L,z) _0(L,z). As is evident from equation (<ref>), the modelling in is independent of _0, which drops out (see also ). In other words, we only need to account for the halo mass dependence of the completeness. As shown in Appendix <ref>, this mass-dependence is already accounted for by our forward-modelling of the Type-I impurities. Hence, we set (M,L,z) = 1 throughout.



  §.§.§ The probability P(,|M,,)



In order to model the line-of-sight kinematics of the secondaries we proceed as follows. Since secondaries consist of both true satellites and interlopers, which have distinct phase-space distribution, we write


    P(, |M,L,z) =     f_ int  P_ int(, |L,z)   + 
        [1 - f_ int]   P_ sat(, |M,L,z)


with the interloper fraction defined as


    f_ int = f_ int(M,L,z) = λ_ int(L,z)/λ_ tot(M,L,z) ,


where λ_ int and λ_ tot have been individually defined in <ref>. We first describe how we compute P_ sat(, |M,L,z) (in <ref>) before detailing our treatment of interlopers (<ref>).



  §.§.§ The phase-space distribution of satellites:

 

In computing the joint 2D probability P_ sat(, |M,L,z), we assume that the baryonic matter of the central galaxy has a negligible impact on the kinematics of its satellite galaxies[We address the accuracy of this assumption, which is common to virtually every study of satellite kinematics, in a forthcoming paper (Baggen et al., in prep).], and we model the satellites as tracers in a pure dark matter halo which is fully characterized by its halo mass and concentration. Throughout, we use the median concentration-halo mass relation of <cit.>, and we emphasize that our modelling is fairly insensitive to the exact choice of the concentration-mass relation within reasonable bounds of its theoretical uncertainty. We also assume the central galaxy to be located at rest at the centre of the halo. As shown in , relaxing this assumption by allowing for non-zero velocity bias for centrals has negligible impact on 's inference. 

Under these assumptions we have that


    P_ sat(, |M,L,z) = P(|M,L,z)   P(|,M,z) .


with


    P(|M,L,z) = 2 π  Σ̅(|M,z)/f_ ap(M,L,z) .


Here f_ ap is defined in equation (<ref>), and


    Σ̅(|M,z) = 2 ∫_^r_ sp(M,z)n̅_ sat(r|M,z) r  r/√(r^2 - R^2) ,


is the projected, normalized number density distribution of satellite galaxies.

In , we made the simplified assumption that the line-of-sight velocity distribution, P(|,M,z), is a Gaussian, which is completely characterized by the line-of-sight velocity dispersion σ_ los(|M,L,z). However, there is no a priori reason why the LOSVD should be Gaussian. In fact, the detailed shape of the LOSVD contains valuable information regarding the velocity anisotropy <cit.>, which we aim to constrain using . In this work we therefore improve upon by extending our modelling of the kinematics to fourth-order and by describing P(|,M,z) as a generalised Gaussian with a projected velocity dispersion, σ_ los(|M,L,z), and a line-of-sight kurtosis, κ_ los(|M,L,z). The projected, line-of-sight velocity dispersion is related to the intrinsic, radial velocity dispersion, σ_r^2(r|M,z), according to 


    σ^2_ los( |M,z) = 2/Σ̅()   ∫_^r_ sp(M,z)[ 1 - β(r|M) ^2/r^2] ×
       n̅_ sat(r|M,z)  σ_r^2(r|M,z)  r   r/√(r^2 - ^2) ,


where σ_r^2(r|M,z) follows from the second order Jeans equation for a spherically symmetric NFW halo <cit.>. Here, the local anisotropy parameter


    β(r|M) ≡ 1 - σ^2_(r|M)2  σ^2_(r|M)


relates the tangential (σ_) and radial (σ_) velocity dispersions. For our fiducial model we assume that β is independent of both radius and halo mass, and we constrain this `average' velocity anisotropy using the satellite kinematics data. In  <ref> we discuss the implications of adopting more flexible models in which the anisotropy parameter is allowed to depend on halo mass. Note that the upper-integration limit of equations (<ref>) and (<ref>) is set to r_ sp(M,z) = f_ sp r_ vir(M,z), instead of r_ vir(M,z), to account for a population of splash-back galaxies (see <ref>). 

The projected, fourth moment of the LOSVD at projected separation  is given by


    v_ los^4(R_) =    2Σ(R_)∫_R_^r_ sp[1-2β R_^2/r^2 + 12β (1+β) R_^4 / r^4 ] ×
       v_r^4(r|M,z)   n_ sat(r|M,z) r   r√(r^2-R_^2) ,


where β is β(r|M) in general, and v_r^4(r|M,z) follows from the fourth-order spherical Jeans equation <cit.>, which for radius-independent anisotropy is given by:


    v_r^4(r   |M,z) = 3   Gr^2β n̅_ sat(r|M,z)×
       ∫_r^r_ sp dr'   (r')^2β n̅_ sat(r'|M,z)  σ_r^2(r'|M,z)  M(r')r'^2


Here M(r) is the enclosed mass of the spherical NFW halo inside radius r. Given the fourth-order line-of-sight velocity moment, we can compute the projected kurtosis as


    κ_ los(R_|M,z) =  v_ los^4  (R_|M,z)   /  σ_ los^4  (R_|M,z).




Finally, in order to account for non-zero redshift errors in the data, the line-of-sight velocity dispersion is modified according to σ_ los→√(σ^2_ los + 2 σ^2_ err), with σ_ err = 15 the typical SDSS redshift error <cit.>. Having computed both the velocity dispersion and kurtosis, we model the detailed shape of the LOSVD, P(|,M,z), using a symmetric (all odd moments are equal to zero), generalized form of the normal distribution, known as the <cit.> distribution[The Langdon distribution is often used to characterize the non-Maxwellian velocity distribution of electrons heated due to the inverse-Bremsstrahlung process <cit.>.]:


    P_() = 12 Γ(1/m) ma_m exp(-|/a_m|^m) .


Here the parameters a_m and m are related to the variance, σ^2, and the kurtosis, κ, according to 


    σ^2 =  a_m^2 Γ(3/m)Γ(1/m)     and    κ = Γ(5/m) Γ(1/m) Γ^2(3/m)


The reason for using this particular distribution function is purely one of convenience; P_() has a nice analytic closed form, is simple to compute, has all the features required of a probability distribution (normalized and positive-definite), and includes the Gaussian as a special case (m=2).

In , we use equation (<ref>) to compute a_m and m from σ^2_ los( |M,z) and κ_ los( |M,z), after which we compute 


    P(|,M,z) = P_()  Γ(1/m)Γ(1/m) - Γ(1/m, (/a_m)^m)


which is properly normalized such that its integral from - to  is unity.



  §.§.§ The phase-space distribution of interlopers
 


In we assumed that interlopers have a constant projected number density and a uniform distribution in line-of-sight velocities, so that P_ int(,  |L,z) =    /  [  (^2 - ^2)]. Here = (L) and = R_ cut(z) are the outer and inner radii of the conical volume used to select secondaries around primaries of luminosity L at redshift z and  is the corresponding line-of-sight depth (see <ref>).

However, as discussed in detail in , a subset of the interlopers are either infalling or splash-back galaxies and have kinematics that are very similar to the true satellites. Assuming that the velocity distribution of interlopers is uniform ignores this `kinematically coupled interloper population', which causes to overestimate the number of satellite galaxies. Although the resulting offsets were modest for the smaller mock samples studied in , they cause a significant, systematic bias (predominantly in the satellite CLF parameters b_0, b_1 and b_2) when using larger samples. This motivated us to develop a more sophisticated treatment for the phase-space distribution of interlopers. 

Based on a detailed assessment of interlopers in our mock data sets (see and <ref> for details), we  now model the interlopers as consisting of three fairly distinct populations: (i) a population of `splash-back galaxies' associated with the host halo of the primary, and extending out to a distance r_ sp from the primary, (ii) a roughly uniform background population of `true' interlopers that are uncorrelated with the primary, and (iii) an `infalling' population of interlopers, located outside of the splash-back radius. This infall motion, on large linear scales, is responsible for redshift space distortions in clustering data known as the Kaiser effect <cit.>.

We assume that the phase-space distribution of splash-back galaxies can be modelled similar to that of the satellites; i.e., they follow the same n_ sat(r|M,z), extrapolated to beyond the halo's virial radius, and their kinematics obey the same Jeans equations. The only difference is that they are located between the host halo's virial radius, r_ vir, and a splash-back radius r_ sp≡ f_ sp r_ vir. In order to account for a population of splash-back galaxies we simply change the upper-integration limit of equations (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>) from r_ vir(M,z) to r_ sp(M,z). Throughout we adopt f_ sp = 2, which is motivated by estimates of the splash-back radius in simulations <cit.>. In addition, detailed tests with mock data sets (see <ref>) show that yields unbiased estimates of the velocity anisotropy parameter β for f_ sp≳ 1.5. We find that not accounting for splash-back galaxies (i.e., setting f_ sp=1) results in a weak bias of β, without significantly affecting any of the other parameters. On the other hand, setting a much larger value for splash-back radius, like f_ sp=3, yields posteriors that are indistinguishable from those for f_ sp=2. Hence, our choice of f_ sp=2 is reasonable and our results are robust against modest changes in the adopted value of f_ sp.

We assume that both the uncorrelated `background' interlopers (bg) as well as the `infalling' interlopers (inf) have a uniform angular distribution on the sky, such that their phase-space distribution can be written as


    P_ int(, ) = 2 /(^2 - ^2) [P_ bg() + P_ inf()] ,


where the term in square brackets is normalized such that its integral from - to + is unity. Since the secondary selection volume is conical in shape, the backside has a larger volume than the front (see Fig. <ref>). Note that all the secondaries of any given primary have luminosities above a fixed threshold. Therefore, due to the conical selection volume, we expect the velocity distribution of the uncorrelated background interlopers, P_ bg(), to increase with . In particular, P_ bg() is proportional to the comoving volume of the corresponding velocity slice of the secondary selection cone, and we therefore adopt


    P_ bg() = 3/H(z') d^2(z')/d^3() - d^3()  (1+z_ pri)


Here  z' = z_ pri + (1+z_ pri) /c and d(z) is the comoving distance out to redshift z.



We have experimented with modelling the line-of-sight velocity distribution of the infalling population of interlopers, P_ inf(), using the linear <cit.> model, but this did not yield sufficiently accurate results. We therefore opted for a semi-empirical approach. Tests with mock data (see <ref>) show that P_ inf() is accurately fit by a Gaussian, P_ inf(Δ V) =  A_ inf ^-1/2(Δ V^2 / σ_ inf^2), with A_ inf and σ_ inf free parameters that vary from primary to primary. Rather than trying to devise an analytical model for these parameters, we use the following data-driven approach. Around each primary we select a set of `tertiary' galaxies in a conical volume similar to that used for the secondaries, but at larger projected distances from the primary. More specifically, the tertiary selection cone is specified by an inner projected radius R_ min^ int = 0.6  σ_200 and an outer projected radius R_ max^ int = 0.9  σ_200, and by the same redshift depth as the secondary selection cone (shown by the outermost hollow annular cone in Fig. <ref>). Tests with mock data (see Fig. <ref>) indicate that (i) the line-of-sight velocity distribution of these tertiaries is virtually indistinguishable from that of the infalling interlopers among the secondaries, (ii) the results are insensitive to the exact radii of the tertiary selection cone[A weak dependence of σ_ inf on projected radius is apparent in Fig. <ref>, and has also been noted by <cit.>. We have experimented with implementing such a R_-dependence and extrapolating this
to the radial interval of the secondary selection cone, but found that this had a negligible impact on the inference.], and (iii) less than one percent of the tertiaries are actual satellites of the primary. The latter indicates that for most primaries R_ min^ int lies well outside the virial radius of the host halo of the primary, as required. We assume that A_ inf and σ_ inf are each quadratic functions of log(L_ pri) and  (including the cross-term), and determine the corresponding 6 × 2=12 coefficients by simultaneously fitting the velocity distribution of all tertiaries around all primaries. These coefficients are then used in to model the line-of-sight velocity distribution of the infalling population of interlopers.

Note that this rather elaborate model for the phase-space distribution of interlopers has zero degrees of freedom. The only degrees of freedom for the interloper-modelling is with regard to their number density, which is modelled via the effective bias described by equation (<ref>). 



 §.§ Modelling the number of secondaries



As mentioned in <ref>, the data vector _ NS that describes the number of secondaries for a random subset of N_ NS primaries, including those with zero secondaries, contains valuable information regarding the occupation statistics of satellite galaxies, and hence, the CLF. 

Similar in spirit to how we compute the likelihood for the satellite kinematics by marginalizing over halo mass (cf. equation [<ref>]), the likelihood for _ NS given a model θ is given by


    _ NS = ∏_i=1^N_ NS∫  P(|, )   P( | , , ) .


Here P(|M,L,z) is given by equation (<ref>), while


    P(M|L,z) = P(L, M, z)/∫ M  P(L,M,z) ,


with P(L,M,z) given by equation (<ref>). 



 §.§ Modelling the galaxy luminosity function


 
The final observational constraint that we use to constrain the halo occupation statistics is the comoving number density of SDSS galaxies, n_ gal(L_1,L_2) in ten 0.15 dex bins in luminosity, [L_1,L_2]  covering the range 9.5 ≤log L/(h^-2) ≤ 11.0 (see <ref>). We include this data in our inference problem by defining the corresponding log-likelihood


    ln_ LF(_ obs|) = -1/2  [() - _ obs]^t  Ψ  [() - _ obs] .


Here _ obs is the data vector and () is the corresponding model prediction, computed from the CLF and the halo mass function using


    n_ gal(L_1,L_2) = ∫_L_1^L_2 L ∫_0^∞Φ(L|M)   n(M,z_ SDSS)   M ,


where z_ SDSS = 0.1 is a characteristic redshift for the SDSS data used,[We have verified that our results do not depend significantly on this choice.] and Ψ is the precision matrix, which is the inverse of the covariance matrix, with <cit.> correction (see <ref>).



 §.§ Numerical implementation



The fiducial model used by is characterized by a total of 16 free parameters: 6 parameters describing Φ_(L|M) (namely log M_1, log L_0, γ_1, γ_2, σ_13 and σ_), 4 parameters describing Φ_(L|M) (namely α_ s, b_0, b_1, and b_2), 1 parameters (β) to quantify the average velocity anisotropy of satellite galaxies, 2 parameters (γ and ) that describe the radial number density profiles of satellite galaxies (see equation [<ref>]), and 3 nuisance parameters  (η_0, η_1 and η_2) that specify the abundance of interlopers (see equation [<ref>]). We assume broad uniform priors on all parameters, except for β. The value of the anisotropy parameter β formally ranges from -∞, for maximal azimuthal anisotropy, to +1, for maximal radial anisotropy, which is difficult to probe with our MCMC sampler. Hence, in order to assure roughly equal amounts of parameter space for radially and azimuthally anisotropic models, we sample ℬ≡ -log(1-β), rather than β. In particular, we adopt uniform priors over the range -1 ≤ℬ≤ 1, which corresponds to -9 ≤β≤ 0.9.




Probing the posterior P(|) over the fiducial 16-dimensional parameter space requires close to a million likelihood evaluations, each of which involves thousands of numerical integrations (see for details). In order to make this problem feasible we perform the Bayesian inference under the assumption of a fixed normalized, radial number density distribution of satellite galaxies, i.e., fixed values for γ and . This has the advantage that f_ ap(M,L,z) and P_ int(, |L,z) are all independent of the model, , while P_ sat(, |M,L,z) only depends on a single anisotropy parameter (see <ref>). We compute P_ sat(, |M,L,z) for each central-satellite pair for 10 values of ℬ between -1 and 1 (or β between -9 and 0.9), and we interpolate it for intermediate values. Combined with the fact that we perform all integrations over halo mass using Gaussian quadrature with fixed abscissas (see ) implies that we only need to compute all these quantities once for each primary and/or secondary, which we then use to find the posterior P(|) in the 14-dimensional parameter space at fixed (γ,). As a consequence, for a single evaluation of the full likelihood 


    ln_ tot = ln_ SK + ln_ NS + ln_ LF ,


for the full SDSS data consisting of 18,373 primaries with at least one secondary and a total of 30,431 secondaries (see <ref>), it only takes of the order of 200 milliseconds using a single run-of-the-mill CPU. This is sufficiently fast, that it allows one to run different Monte-Carlo Markov Chains for different assumptions regarding  n̅_ sat(r|M,z), and to find the best-fit radial profile, marginalized over all other model parameters. First, we combine the posteriors from separate MCMC runs on 15× 15 grid in (γ, log)-space, each time marginalizing over the other 14 parameters, to constrain n̅_ sat(r|M,z) (see Appendix <ref>). Having determined the values of γ and  that maximize _ tot, we then run a MCMC sampler to infer the full posterior P(|) keeping γ and  fixed at these best-fit values. The MCMC sampler used to probe our 14 dimensional parameter space is the affine invariant stretch-move algorithm of <cit.>. Throughout we use 1,000 walkers and the proposal density advocated by <cit.>. This results in typical acceptance fractions between 0.3 and 0.4, and the MCMC chain is typically converged after about 500 steps (i.e., 5 × 10^5 likelihood evaluations). We have experimented at length with other initial guesses, and find the results to be extremely robust, and to always fully converge well under 1 million likelihood evaluations. Finally, throughout we adopt flat priors on all parameters, with very broad prior bounds that do not affect our inference.








§ VALIDATION WITH MOCK DATA





 §.§ Tier-3 mock data



We validate the performance of using the Tier-3 mock data introduced in . This mock sample is constructed using the z=0 halo catalogue of the high-resolution SMDPL simulation <cit.>, which uses 3840^3 particles to trace structure formation in a cubic volume of (400 )^3, adopting cosmological parameters consistent with <cit.>.

Each host halo in the catalogue with a mass M_ vir≥ 3 × 10^10 is populated with mock galaxies with luminosities L ≥ 10^8.5 according to a particular fiducial CLF model. Each central galaxy is given the position and velocity of its halo core, defined as the region that encloses the innermost 10% of the halo virial mass. Satellite galaxies are assigned the phase-space coordinates of the subhaloes with the highest peak halo masses. If the number of satellites, drawn from the input CLF, exceeds the number of resolved subhaloes in a specific host halo, we randomly assign the excess satellites the halo-centric positions and velocities of subhaloes hosted by other haloes of similar mass. Note that no assumption of quasi-equilibrium dynamics has been made in the mock making procedure. Therefore, our Tier-3 mock satellites obey the Jeans equations only as much as the live subhaloes do in the SMDPL simulation.

Once all haloes have been populated with mock galaxies, we construct a mock SDSS survey as follows.  First, we place a virtual observer at a random position within the simulation volume. We use this virtual observer to convert the (x, y, z) coordinates of each galaxy into a cosmological redshift, z_ cosm, and sky coordinates (using a random orientation). If necessary, the simulation box is repeated with random sets of right angled rotations until the entire cosmological volume out to z_ cosm = 0.20 is covered. Next, we overlay the SDSS DR7 footprint on the simulated sky, and only keep galaxies with m_r ≤ 17.6 that lie within the SDSS DR7 survey window. Redshift-space distortions are simulated by adding (1 + z_ cosm) v_ los / c to the redshift of each galaxy, with v_ los the galaxy's peculiar velocity along the line-of-sight. Spectroscopic redshift errors in the SDSS are simulated by adding a random Δ z from a Gaussian with scatter σ_ err = 15 <cit.>. Finally, we simulate the effect of fibre collisions induced spectroscopic incompleteness following the method of <cit.>. Once the mock spectroscopic survey is completed, we select primaries and secondaries using the selection cones described in <ref>, and assign spectroscopic weights to all secondaries using the method described in <ref>. Similar to what we do for the real data, we remove primaries with an aperture completeness w_ app< 0.8 and exclude secondaries that are located within 55” from their primary. Finally, we use the mock data to compute the comoving abundances of galaxies in the ten luminosity bins described in <ref> using the same method as used for the real SDSS data (i.e., by taking into account the SDSS DR7 footprint with its mask and window functions as well incompleteness caused by fiber collision).









 §.§ Inference from the Tier-3 mock



The Tier-3 mock data described above is used to test the performance of . As described in <ref>, we first determine the best-fit values of γ and , which characterize the radial number density profile of satellite galaxies, properly marginalized over all other model parameters. In the Tier-3 mock, satellites are placed in subhaloes which are known to have a radial profile that differs starkly from that of the dark matter. The true radial number density distribution closely follows the generalised NFW shape (equation <ref>) with (γ,) ≈ (0.0, 2.6). Hence, the radial profile of satellites, in our Tier-3 mock, is cored and has a scale-radius that is significantly larger than that of their dark matter host haloes. This feature of the DM-only simulation is consistent with many previous studies <cit.>. The best-fit values obtained by , when applied on the Tier-3 mock data is (γ,) = (0.0, 2.5), in almost perfect agreement with the profile inferred directly from the N-body simulation. Thus, in agreement with what we reported in , can accurately recover the radial profile of satellite galaxies. Next, we keep (γ,) fixed at the best-fit values, (0.0, 2.5), and run to infer the posterior distribution of the remaining 14 parameters. The brown lines in Fig. <ref> show the CLF constraints thus obtained for our fiducial model[No results are shown for the three nuisance parameters η_0, η_1 and η_2 that characterize the abundance of interlopers, or for the velocity anisotropy parameter β, which is discussed in detail in <ref>.]. Note that the posteriors of all CLF parameters are in excellent agreement with the input values, shown as blue dashed lines. For comparison, the posteriors indicated in green correspond to a model that allows for mass-dependence of the orbital anisotropy, and will be discussed in <ref>.

Panel (a) of Fig. <ref> demonstrates that accurately fits the galaxy luminosity function. The blue circles show the mean number density of galaxies in logarithmic bins of luminosity of width 0.15 dex (roughly the width of the blue vertical bars), while the brown band indicates the 95% confidence interval as inferred from 's posterior.



Fig. <ref> panel (b) shows the relationship between halo mass and median central galaxy luminosity. Note that the true input model (blue squares) is perfectly recovered by (brown band, indicating the 95 % confidence interval from the posterior). The pink dots indicate the true halo masses and luminosities of all primaries (both true centrals and impurities) in the Tier-3 satellite kinematics sample. The magenta triangles show their median values of log in bins of log M. Note that, due to selection bias, these do not agree with the true input model (blue squares). A more luminous primary has a larger secondary selection volume associated with it,  as well as a larger luminosity range between  and L_ min(). Therefore, at any given redshift, brighter primaries are more likely to have at least one secondary than less luminous primaries in haloes of the same mass, hence the brighter primaries get preferentially selected in the _ SK data. This causes the median luminosity of primaries in the satellite kinematics sample to be biased high with respect to that of true centrals. This effect is especially pronounced at the low halo mass end, where the expectation value for the number of secondaries is lowest.  Despite this strong and unavoidable selection bias, perfectly recovers the true input relation between central luminosity and halo mass, and not the biased relation! This indicates that , in its forward modelling, accurately accounts for selection bias and other systematics such as the presence of impurities. 

Panels (c) and (d) of Fig. <ref> show, respectively, the logarithmic scatter in central galaxy luminosity at fixed halo mass, σ_, and the normalization of the satellite CLF, ϕ_^*. Once again, blue dots show the true input values, while the brown shaded regions mark the 95% confidence interval of the posterior distribution inferred with . Finally, the brown-shaded histograms in panels (e) and (f) show, respectively, the posterior distributions for the faint-end slope of the satellite CLF, α_, and that of the orbital anisotropy parameter, β. The blue, vertical line in panel (e) indicates the true input value of α_, while the blue-shaded region in panel (f) show the range of mean velocity anisotropy of subhaloes in different halo masses in the SMDPL simulation used to construct the Tier-3 mock (discussed further in <ref>). As is evident, in each case the posterior constraints are in excellent agreement with the input values, which indicates that can put tight and accurate constraints on the more intricate aspects of the galaxy-halo connection, beyond the mere relation between halo mass and central luminosity. 

Fig. <ref> shows the interloper fraction, f_ int, as a function of primary galaxy luminosity in six redshift bins (indicated by different colours). Solid circles indicate the true interloper fractions in the Tier-3 mock sample with the error bars computed assuming Poisson statistics. The coloured bands show the corresponding posterior predictions as inferred by , which are in good agreement with the true interloper fractions. Hence, correctly distinguishes satellites from interlopers, at least in a statistical sense, and accurately recovers their relative prevalence as a function of luminosity and redshift of the primary.

Finally, Figs. <ref> and <ref> compare the posterior predictions for the HOD and CLF, respectively, to their true input used to construct the Tier-3 mock data (solid dots). In particular, Fig. <ref>, shows the average number of central (circles) and satellite (triangles) galaxies per halo, above four different luminosity thresholds (L_ th), as a function of host halo mass. Interestingly, the true satellite mean occupation in our mock data deviates significantly from a simple power-law (which is a common assumption in the literature), and accurately recovers that complexity in its shape across all luminosities and halo masses. Fig. <ref> plots the central (purple) and satellite (orange) CLFs for 8 different halo masses, as indicated in the top-right corner of each panel. Note that in all cases the halo occupation statistics are recovered with exquisite precision and accuracy.




 §.§ Velocity Anisotropy in the Tier-3 Mock



Unlike previous studies of satellite kinematics, which used satellite velocity dispersions in bins of primary luminosity, a unique aspect of is that it models the full probability distribution P(, |M,L,z). By modelling the full 2D phase-space information, rather than just the second moment of velocity, has the potential to break the mass-anisotropy degeneracy that hampers dynamical models which only rely on measurements of the satellite velocity dispersions. Hence, is able to simultaneously constrain the halo masses of the primaries, as well as the velocity anisotropy of its secondaries. 

The brown shaded contours in Fig. <ref> show the 68 and 95 percent confidence intervals on the orbital anisotropy parameter β as inferred by from the Tier-3 mock data. Note that the constraints (β = 0.11 ± 0.05) are remarkably tight and indicate a mild, radial anisotropy. Recall that the Tier-3 mock was constructed by placing satellite galaxies inside subhaloes in the SMDPL simulation. Hence, these satellite galaxies have the same orbital anisotropy as those subhaloes. The solid, black line indicates the average anisotropy parameter of subhaloes in the SMDPL simulation as a function of host halo mass. It, too, indicates a mild radial anisotropy in reasonable agreement with the constraints obtained with . However, this comparison is not entirely fair. After all, our secondary selection criteria typically only selects secondaries with projected separations R_ ap^ sec 0.4 (see <ref>). Therefore, it is more meaningful to compare our posterior constraints with the orbital anisotropy of subhaloes located in the inner regions of their host haloes. The blue, solid line in Fig. <ref> shows the orbital anisotropy of subhaloes with a 3D halo-centric radii less than 0.4. These are in better agreement with 's posterior constraints, especially given that the satellite kinematics data is dominated by haloes above 10^12. Hence, we conclude that , by using the full line-of-sight velocity distributions of the secondaries, is indeed able to break the mass-anisotropy degeneracy and obtain simultaneous, reliable constraints on both halo mass and orbital anisotropy. 

Thus far we have only considered models in which the orbital anisotropy is a   `universal' constant, independent of halo mass or halo-centric radius. However, a detailed analysis of the orbital anisotropy of dark matter subhaloes in the SMDPL simulation (and hence our Tier-3 mock satellites) reveals a rather complicated dependence on both halo mass and halo-centric radius (see Fig. C1 in ). Indeed, the blue and black curves in Fig. <ref> indicate some dependence on halo mass, albeit weak. 
We therefore also analysed the Tier-3 mock data using a more flexible model with a mass-dependent orbital anisotropy parameter, given by β(M) = 1 - 10^-ℬ(M), with
 

    ℬ(M) =
    ℬ_12,    if M < 10^12
    ℬ_14,    if M > 10^14
    ℬ_12 + 1/2 (ℬ_14-ℬ_12) log M_12,    otherwise


Here both ℬ_12 and ℬ_14 are free parameters for which we adopt uniform priors ranging from -1 to +1. Hence, by replacing our fiducial model, in which β is independent of halo mass, with this mass-dependent model we add one extra free parameter to the mix. 






Remarkably, we find that this extra degree of freedom has no discernible impact on the constraints of any of the other parameter. This is evident from Fig. <ref>, where the green contours show the posterior constraints from the mass-dependent β model, which are indistinguishable from those of our fiducial model (brown contours). The green shaded contours in Fig. <ref> show the 68 and 95 percent confidence intervals on  β(M) of the corresponding model. It reveals a weak hint for the orbital anisotropy to become more radially anisotropic in lower mass haloes, in agreement with the weak trend for the SMDPL subhaloes with r < 0.4. However, the uncertainties at the low mass end are rather large, and  's inference is consistent with a constant, mass independent β at the 2σ level. 










§ APPLICATION TO SDSS





 §.§ The Data



Having demonstrated the can accurately infer the galaxy-halo connection from kinematic data of primaries and secondaries that can be extracted from a galaxy redshift survey, we now apply to data from the SDSS. In particular, we use the New York University Value-Added Galaxy Catalog <cit.>, which derives from the Seventh Data Release of the SDSS <cit.>. More specifically, we use the VAGC  sample[<http://sdss.physics.nyu.edu/lss/dr72/bright/0/>], which includes ∼ 570,000 galaxies with a limiting Petrosian magnitude of m_r < 17.6. We use this data to identify primaries and secondaries using the selection criteria outlined in <ref>.  As already mentioned there, we limit our analysis to galaxies in the redshift range 0.02≤ z ≤ 0.20. In order to assure that the secondary selection cones fit entirely within this redshift range, the redshifts of primaries are restricted to 0.034 ≤ z ≤ 0.184 (see Fig. <ref>). Primaries are also limited to have luminosities in the range 9.504 ≤log[/(h^-2)] ≤ 11.104. We end up with a total of ∼ 165,000 primaries[For the computation of  _ NS, the number of primaries is downsampled by an order of magnitude to 16,491, as discussed in <ref>.], of which N_+=18,373 primaries have at least one secondary. The total number of secondaries, and thus the total number of primary-secondary pairs for which kinematic data is available, is 30,431. 


 
The upper left-hand panel of Fig. <ref> shows the luminosity distributions of all primaries (blue), primaries with at least one secondary (orange), and secondaries (green). Note that primaries with at least one secondary are brighter, on average, than those with zero secondaries. This simply follows from the fact that, on average, brighter centrals reside in more massive haloes, which host more satellites. Also, given a fixed halo mass the brighter primary is assigned a larger secondary selection volume, making it more likely to have a secondary. This latter effect, though, is subdominant.  Note also that there are no satellites with L < 10^9. As discussed in <ref>, this is a consequence of the apparent magnitude limit of the SDSS survey combined with the fact that we only allow for primaries brighter than 10^9.504. The upper right-hand panel shows the probability, P_0, that a primary of luminosity  contains zero secondaries. It is simply defined as the fraction of primaries, in a given luminosity and redshift bin, with zero secondaries, i.e., P_0 = N_0 / (N_0 + N_+), where N_0 and N_+ are the number of primaries in our sample with zero and at least one detected secondary, respectively. These probabilities have been computed using a 8× 6 uniformly-spaced grid in (log, ) covering the ranges [9.504,11.104] and [0.034,0.184], respectively. Different colours correspond to different redshift bins, as indicated. We emphasize that this binned data are not used in our inference; it is merely used here to illustrate how P_0 scales with luminosity and redshift. Errorbars are computed assuming Poisson statistics, and are smaller than the data points in most cases. Note that P_0 increases with decreasing luminosity and increasing redshift, as expected from the Malmquist bias resulting from the apparent magnitude limit of the spectroscopic SDSS data. The lower left-hand panel shows the multiplicity function, i.e., the number of primaries each with N_ sec secondaries. Note that most primaries have zero secondaries, and that there are very few primaries with more than 20 secondaries. Finally, the lower right-hand panel shows the distribution of the aperture completeness w_ app for all primaries. Note that the vast majority of primaries have w_ app=1, i.e., their entire secondary selection cone completely falls within the SDSS survey footprint. As mentioned in <ref>, primaries with w_ app < 0.8 have been removed from the sample. 


Figure <ref> plots the line-of-sight velocity difference  as a function of the luminosity of the primary galaxies (left-hand panel) and as a function of the projected distance between primary-secondary pairs (right-hand panel). Data points are colour-coded according to the redshift of the primary, as indicated. This constitutes the satellite kinematics data in SDSS, that we attempt to forward-model with in order to constrain the galaxy halo connection. The deficit of data points at small R_ p reflects that we have removed secondaries with a projected separation less than θ_ fc = 55” because of fiber collision issues (see <ref>). Similarly, the absence of data points for low  and large |Δ V| reflects the luminosity dependence of the secondary selection criteria (see <ref>). Evidently, that the velocity dispersion of secondaries is a strong function of primary luminosity, consistent with the expectation that more luminous centrals reside in more massive dark matter haloes. The low-density, high velocity wings of P() at any given  reflects the contribution of foreground and background interlopers, i.e., galaxies selected as secondaries that do not reside in the dark matter halo of the primary. 












 §.§ Results



We start our analysis of the SDSS data presented above by constraining the satellite radial number density profile (equation [<ref>]). Using a 15× 15 grid in (γ, log)-space, we obtain (γ,) = (0.94, 1.7) (see Appendix <ref> for details), in good agreement with the results of <cit.>
and <cit.>. Note that this central slope, γ, is significantly steeper than what we inferred for the Tier-3 mock data (γ=0.0, see <ref>), in which the satellites were positioned on subhaloes in numerical simulations. This indicates that the radial distribution of real satellite galaxies is more centrally concentrated than that of subhaloes in DM-only numerical simulations. This discrepancy is most likely is due to a combination of artificial disruption in simulations <cit.> and failures of the subhalo finders being used <cit.>.

Next, keeping  and γ fixed at 1.7 and 0.94, respectively, we run to constrain the posterior distribution of the remaining 14 parameters that characterize the CLFs of central and satellite galaxies, the interlopers, and the average orbital anisotropy of satellite galaxies. Once again, we adopt very broad non-informative priors for all parameters. Table <ref> lists the best-fit parameters plus their 68% confidence intervals thus obtained. We emphasize that, as shown in and Appendix <ref>, all these results are extremely robust to modest changes in  and γ.

Before showing the key results on the CLF, we first demonstrate that the best-fit model of is an excellent fit to the data. In order to illustrate this, we bin the data in 2D-bins of luminosity and redshift of the primaries. We emphasize that no such binning was used in the analysis; it is merely used here for the purpose of visualization of the data and its corresponding prediction from . The various panels in Fig. <ref> show the LOSVDs of primary-secondary pairs for bins in log (different rows) and  (different columns). We only show panels for which the luminosity lower bound of the {log,} bin falls above the flux limit at the redshift upper bound of that bin.  Blue dots and shaded histograms show the stacked data, while the red shaded bands show to 95% confidence intervals obtained using the inferred posterior distributions. In order to quantify the level of agreement between the data and the model, we proceed as follows. Let N_(,) be the actual number of secondaries in the SDSS data for each of the various {log,} bins. Using the best-fit model's predicted P(Δ V) for each {log,} bin we draw N_(,) values of Δ V, and compute the likelihood of this fake data representing the best-fit model. We repeat this 10^4 times. The red-shaded histogram in the inset-panel in the lower-right corner of Fig. <ref> shows the resulting distribution of likelihoods, which we then compare to the analogous likelihood for the actual SDSS data (blue, vertical dashed line). The fact that the latter is perfectly consistent with the distribution of likelihoods (red histogram) indicates that  the LOSVDs obtained from are in excellent agreement with the data, fitting not only the roughly Gaussian LOSVDs centered on Δ V=0, but also the extended wings due to the interlopers.

Fig. <ref> uses the same (log,) binning and panels[Note that bins with five or fewer primaries in the down-sampled 
_ NS data have been omitted.] as Fig. <ref>, but this time we plot the distributions of the numbers of secondaries per primary (red dots with Poisson errorbars).  More specifically, the x-axis indicates the number of secondaries, N_ sec, and the y-axis indicates log(1+N_ pri), where N_ pri is the number of primaries that each have N_ sec secondaries. In most cases the distributions clearly peak at N_ sec=0, as most primaries in SDSS DR7 do not have a spectroscopically detected secondary. The only exceptions are a few high-log L_ pri bins. Once again, in each panel the red shaded bands indicate 95% confidence intervals obtained using the posterior distributions inferred with .
The inset-panel in the lower-right compares the likelihood of the SDSS data given the best-fit model, to the distribution of expected likelihoods computed by drawing 10^4 random realisations of the best-fit multiplicity function. This time the likelihood of the SDSS data falls at the edge of the expected range, indicating that the fit to the data is not optimal. Indeed, upon closer inspection one can notice that the best-fit model overpredicts the multiplicity of primaries with N_ sec∼ 3-6, especially for some intermediate  and  bins. As we demonstrate in a forthcoming paper, this small discrepancy arises from certain limitations in the satellite CLF model, and can be resolved by adopting a slightly more flexible halo-occupation modelling without significantly affecting any of the main relations presented here.


Fig. <ref> shows several key halo mass dependencies that characterize the galaxy-halo connection inferred by from the SDSS data. In each panel, the brown shaded bands show the inferred 95 % confidence intervals, while the coloured symbols show best-fit constraints from previous SDSS-based studies[Where needed we have converted these results to our definition of halo mass.]. In particular, we compare our inference to the results from an analysis of galaxy group catalogues by <cit.>, to results based on a simultaneous analysis of galaxy clustering and galaxy-galaxy lensing by <cit.>, and to results of the most recent analysis of satellite kinematics by <cit.>. Note that the latter did not use the Bayesian hierarchical methodology of , but rather was based on the standard summary statistics of host-weighted and satellite-weighted  velocity dispersions as a function of binned primary luminosity.



Panel (a) plots the median central luminosity, L_, as a function of halo mass. As is evident, the constraints obtained with are in excellent agreement with previous results, though we emphasize that our constraints are significantly tighter.

Panel (b) plots the posterior constraints on σ_(M), characterizing the scatter in log L_ at fixed halo mass. Most studies in the past have assumed σ_ to be independent of halo mass, and inferred values that lie roughly in the range 0.15-0.2 dex <cit.>. , on the other hand, allows for a mass dependence as characterized by equation (<ref>). Yet, despite this extra degree of freedom, our inference is statistically consistent with a constant σ_ = 0.17 dex over the entire range of halo masses probed.  Note that this is different from <cit.> who, using a similar mass-dependent characterization of the scatter in the log L_ - log M relation, inferred that σ_ increases with decreasing halo mass, as depicted by the blue circles in Panel (b). As discussed in more detail in <ref>, the reason for this discrepancy can be attributed to the fact that all previous analyses, including <cit.>, invariably assumed the brightest halo galaxy to be the central.

Panel (c) shows the posterior constraints on the the faint-end slope, α_, of the satellite CLF. Throughout we have assumed a global, mass-independent α_ similar to <cit.>, <cit.> and most other previous work. Our inference that α_ = -0.87 ± 0.06 is in excellent agreement with <cit.>, and is largely consistent with the constraints obtained by <cit.> and <cit.> given their uncertainties. Note, though, that <cit.> inferred that α_ becomes significantly steeper at the massive end, reaching values as low as -1.5 for groups with an inferred halo mass M  3 × 10^14. In future work we plan to allow for a mass-dependent α_, to see if satellite kinematics reveal a similar mass dependence as that inferred from the galaxy group catalogue of <cit.>, and to study how this extra degree of freedom impacts the other parameters that characterize the galaxy-halo connection. 

Finally, panel (d) of Fig. <ref> shows the constraints on the normalization, ϕ_^*, of the CLF of satellite galaxies, as a function of host halo mass. The constraints obtained by , as depicted by the brown bands, are in fair agreement with previous constraints, especially if the uncertainties on the latter are taken into account (note that the coloured symbols only indicate the best-fit values). At the low mass end, the results of <cit.>, which are also based on satellite kinematics, seem to suggest significantly larger values of ϕ_^* (i.e., more satellites per halo). This is likely due to the fact that  <cit.> have assumed that interlopers have a uniform distribution of line-of-sight velocities, <cit.>. As discussed in <ref>, this oversimplified assumption implies that some of the splash-back galaxies and infalling interlopers, which have a LOSVD that resembles that of true satellites, will be incorrectly `counted' as satellite galaxies. Because of this, and since we have demonstrated that can accurately recover the interloper fraction as well as the CLF normalization, ϕ_^*(M), we reckon that our results for the abundance of satellite galaxies are likely to be more accurate. In all fairness, we emphasize that most previous studies of satellite kinematics were not aiming to accurately recover satellite abundances; rather they mainly focused on constraining halo masses as a function of primary galaxy luminosity. As discussed in detail in <cit.>, this aspect of the analysis is not significantly impaired by the oversimplified assumption that the LOSVD of interlopers is uniform. 




 §.§ Scatter in central galaxy luminosity



Different empirical estimates of the galaxy-halo connection, at present, broadly agree on the relation between central galaxy luminosity (or stellar mass) and halo mass. However, the constraints on the scatter in this relation, especially as a function of halo mass, has yet to attain similar convergence, and therefore is considered a key parameter at the forefront of empirical modelling <cit.> and is highly informative for testing physical models <cit.>.





's inference of the logarithmic scatter in central luminosity, σ_(M_ vir), shows no significant halo mass dependence despite having the freedom in the model. Fig. <ref> compares our constraints on σ_ (grey band) with estimates from previous studies (all error-bars and uncertainties in this plot are 68% confidence intervals). Crucially, our constraints disagree with <cit.>, who also used satellite kinematics extracted from the SDSS DR-7 data. They split their sample into red and blue centrals, and the scatter in the two sub-populations are shown by the shaded regions of corresponding colour. At the high halo mass end, their red-fraction of centrals approaches unity, and thus the red-shaded region should be a good approximation of the overall scatter. As is evident, it reveals weak but significant mass-dependence with dσ_/ dlog M ≈ -0.04, which, at first sight appears inconsistent with our results. However, <cit.>, as all other previous studies, have simply assumed that their primary, defined as the brightest galaxy in the selection cone, is always the central galaxy. Hence, their σ_ has to be interpreted as the scatter in the brightest halo galaxy, log L_ BHG,  rather than that in log L_.

, on the other hand, accounts for type-I impurities, which are BHG-satellites that are misclassified as primaries. In particular, we have demonstrated that by directly forward modelling these BHG-satellites, 's inferred σ_ is an unbiased estimate of the intrinsic luminosity scatter of true centrals. The probability that the BHG is a satellite, rather than a central, increases strongly with halo mass <cit.>. Therefore, the inferred scatter at the high mass end, from studies that did not account for type-I impurities, may have been biased. We can directly test this with , which makes it straightforward to compute the expected scatter in log L_ BHG and compare it to that in log L_. The black dashed curve in Fig. <ref> shows the predicted scatter in BHG luminosity, σ_ BHG, as a function of host halo mass for our best-fit CLF model. Note that σ_ BHG drops significantly below σ_ at the high-mass end. This is because it is mostly the fainter centrals that are `replaced' by a brighter satellite, which causes the distribution of BHG luminosities to be narrower than that of the true centrals. For M  10^13.5, the mass-dependence of the inferred σ_ BHG is in excellent agreement with <cit.> (recall that the vast majority of all centrals in massive haloes are red, and the comparison should thus be with the red-shaded region). Above 10^13, the black dashed line also shows improved agreement with previous results from an analysis of galaxy clustering and galaxy-galaxy lensing by <cit.>, who assumed that σ_ is mass-independent, and that from an analysis of a SDSS galaxy group catalogue by  <cit.>.

As is evident from Fig. <ref>, the various results disagree strongly at the low-mass end (M  10^13). This, however,  has to be interpreted with caution, as none of the constraints are particularly reliable there. For example, the results of <cit.> can not be trusted at the low-mass end, since their halo masses are estimated from the total group luminosity. Since this is dominated by the central luminosity in low mass haloes, their inferred scatter at the low mass end is guaranteed to be an underestimate <cit.>. In the case of <cit.>, most of the constraining power comes from haloes with M ≥ 10^13. As they assume a simple linear dependence of σ_ on log M, the constraints at the low-mass end mainly reflect extrapolation of the assumed linear relation. Our results are also affected this way, but less so, since we have used a flux-limited sample, rather than a more restricted volume-limited sample as in <cit.>. This allows for a better sampling of fainter centrals, that reside in lower-mass haloes. 



 §.§ Orbital Anisotropy of Satellite Galaxies in SDSS



The brown contours in Fig. <ref> show the 68 and 95 percentile constraints on the orbital anisotropy parameter, β, for satellite galaxies in the SDSS data as inferred from our fiducial model. We infer a significant radial anisotropy with β = 0.29^+0.05_-0.04. These global constraints on the average orbital anisotropy of satellite galaxies, across a large range of halo masses, are perfectly consistent with, but significantly tighter than, the results of <cit.> who also analysed the kinematics of satellite galaxies in SDSS data to infer  β=0.26 ± 0.09. We emphasize that, unlike , the analysis of Wojtak & Mamon did not account for mass mixing, and was based on a much smaller sample of primary-secondary pairs than used here.

Interestingly, our constraints on the orbital anisotropy are also consistent with the typical orbital anisotropy of subhaloes in numerical simulations of structure formation in a ΛCDM cosmology. In fact, the green contours show the constraints we obtain using a model in which the orbital anisotropy is allowed to depend on halo mass as given by equation [<ref>]. We find a weak indication that the orbits of satellite galaxies become more radially anisotropic towards lower halo mass. Most importantly, these results for the SDSS data (Fig. <ref>) are consistent with those for the Tier-3 mock data (Fig. <ref>), in which the satellite orbits reflect those of subhaloes in the ΛCDM-based SMDPL. The fact that the orbital anisotropy of satellite galaxies in the SDSS appears to be consistent with that of subhaloes in N-body simulations can be heralded as yet another success for the ΛCDM model.




Although the weak mass-dependence of the orbital anisotropy inferred here is intriguing, especially in light of the agreement with the Tier-3 results, we emphasize that these results have to be interpreted with caution. The reason is that we have excluded data on projected separations < 55” because of fibre collision issues. As a consequence, the range in radii probed, in terms of the halo virial radius, in low mass haloes is different than that probed in more massive haloes. Hence, any potential radial dependence of the orbital anisotropy of satellite galaxies can, in principle, masquerade as a mass dependence in our analysis. We intend to address this `degeneracy' in a forthcoming paper (Mitra et al., in prep) in which we consider models in which the orbital anisotropy is allowed to depend on halo-centric radius, as well as halo mass. In particular, we will consider Osipkov-Merritt model <cit.>, as well as more realistic simulation-inspired models such as those used by <cit.>.







§ SUMMARY AND CONCLUSION



In we presented , a novel, Bayesian hierarchical method for analysing the kinematics of satellite galaxies. Based on the spherically symmetric Jeans equations it models the kinematics of large ensembles of satellite galaxies associated with central galaxies that span a wide range in halo mass and luminosity. The halo masses of the individual centrals act as latent variables in a hierarchical Bayesian framework that uses the data to  constrain the detailed galaxy–halo connection as characterized by the CLF.  Unlike traditional methods for analysing satellite kinematics, does not make use of any summary statistic, such as velocity dispersions of satellite galaxies in central galaxy luminosity bins. Rather, it leaves the data in its raw form, which has the advantage that all data are used optimally while avoiding systematics that arise from binning. In addition, whereas traditional methods typically require volume-limited samples, can be applied to flux limited samples, thereby greatly enhancing the quantity and dynamic range of the data. And finally, is the only available method that simultaneously solves for halo mass and orbital anisotropy of the satellite galaxies, while properly accounting for `mass-mixing'.

In this paper we have presented a number of important improvements to , required for an unbiased recovery of all parameters when using large samples of data comparable to what can be achieved with existing SDSS catalogues. In particular, 




  * We introduced an improved selection of primaries and secondaries that assures that the secondaries associated with each individual primary are volume-limited, even-though the overall sample is still flux-limited. This facilitates a more accurate modelling of the abundance and velocity distribution of the secondaries.


  * We forward model the contribution of impurities among the primaries, where impurities are predominantly those satellites that are brighter than their corresponding centrals. 


  * We slightly modified the selection criteria of primaries to minimize the effect of other kinds of impurities that are extremely difficult to forward-model.


  * We extended the satellite kinematics model to higher-order, by using the fourth-order Jeans equation to compute the kurtosis of the LOSVD. Incorporating this, in the modelling of the full 2D phase-space distribution of satellites, allows to break the mass-anisotropy degeneracy, and to put tight constraints on the global, average velocity anisotropy of satellite galaxies.


  * We drastically improved the modelling of interlopers by (i) accounting for the fact that the selection volume of secondaries is conical rather than cylindrical, (ii) accounting for splash-back galaxies, and (iii) using a data-driven method to model the line-of-sight velocity distribution of the interlopers.


  * Instead of discarding primaries with zero secondaries, utilizes their information to further constrain the galaxy-halo connection. Congruous with the satellite kinematics methodology, we introduced a similar Bayesian hierarchical framework to model the abundance of secondary galaxies around each primary, which improves significantly on the stacking-based approach used in . This allows to put unprecedented constraints on the satellite CLF.



Using realistic mock data of similar quality and volume as the SDSS DR7, we have demonstrated that, with this improved methodology, can break the mass-anisotropy degeneracy, and simultaneously constrain the host masses and average orbital velocity anisotropy of satellite galaxies. In particular, achieves an unbiased recovery of all 10 CLF parameters that characterize the galaxy-halo connection covering almost four orders of magnitude in halo mass (from ∼ 10^11 to 10^15), and with unprecedented accuracy. In addition, it simultaneously recovers the orbital anisotropy parameter, β, the luminosity and redshift dependence of the interloper fraction, and the radial number density profile of satellite galaxies. It is worth emphasizing that the recovery is unbiased despite the fact that the selection of primaries and secondaries is (unavoidably) plagued by biases, incompleteness, and impurities.

We applied to a sample of 18,373 primaries and 30,431 secondaries extracted from the SDSS DR-7 data, yielding some of the tightest constraints on the galaxy-halo connection to date (Table <ref>). The model accurately reproduces both the abundance and line-of-sight velocity distributions of secondaries (Figs. <ref> and Fig. <ref>), and is in good agreement with previous constraints on the galaxy-halo connection derived from galaxy group catalogues, galaxy clustering, galaxy-galaxy lensing and previous analyses of satellite kinematics. 

Assuming that the orbital anisotropy of satellite galaxies is independent of halo mass and halo-centric radius, our analysis of SDSS data reveals a significant radial anisotropy of β = 0.29^+0.05_-0.04, in excellent agreement with, but significantly tighter than, previous results <cit.>. We also find a weak indication that β is slightly larger in lower mass haloes, in good agreement with the orbital anisotropy of subhaloes in dark-matter only simulations of structure formation in a ΛCDM cosmology <cit.>. Since satellite are believed to reside in subhaloes, this may be considered another success of the standard model for structure formation.

We find that the radial number density profile of satellite galaxies, n_ sat(r|M), is tightly constrained and well characterized by a generalized-NFW profile (equation [<ref>]) with a central cusp-slope γ=0.94 (compared to γ=1 for a pure NFW profile), and a characteristic scale radius that is roughly two times larger than what is expected for the dark matter. Within the uncertainties, this is consistent with several previous studies <cit.>. Consistent with paper I, we find our results to be extremely robust to modest changes in n_ sat(r|M); the only parameter that displays some dependence is the anisotropy parameter, β (see Appendix <ref>). This is to be expected given that both β and n_ sat(r|M) appear in the Jeans equation used to model the kinematics of the satellite galaxies.

Interestingly, we find no evidence for a significant halo mass dependence of the scatter in central luminosity, and at any given halo mass we find the luminosity scatter to be around σ_ = 0.17   dex. This is inconsistent with the latest analysis of satellite kinematics by <cit.>, also based on SDSS DR-7 data, who inferred that the scatter decreases with increasing halo mass, albeit only weakly with σ_/log M ∼ -0.04. As discussed in <ref>, this discrepancy arises primarily from the fact that <cit.> and all previous studies simply assumed the brightest galaxy in the halo to be the central. We, however, take into account the existence of brightest halo galaxy satellites, and forward-model the probability of misidentifying them as primaries. By doing so, we demonstrate that 's inference of σ_ is an unbiased recovery of the intrinsic luminosity scatter of true centrals. From our best-fit model we can predict what the scatter in brightest halo galaxy luminosity should be, as a function of halo mass, and that is in good agreement with the scatter inferred by <cit.> and other previous analyses.

For completeness, we point out that several studies that used stellar mass, rather than r-band luminosity, to characterize the galaxy-halo connection also inferred that the scatter in stellar mass of central galaxies decreases with increasing halo mass <cit.>. However, most of these inferences were only significant below the ∼ 2 σ level. Hence, observationally  it remains unclear whether or not the scatter in the galaxy-halo connection has a significant mass dependence. Taking our results at face value, it seems that scatter is fairly mass-independent, at least for log M  10^12, and that previous indications for a significant mass dependence at the massive end are likely a result of confounding true centrals with brightest halo galaxies. 

On the theory side, the situation is even more higgledy-piggledy, with a clear lack of consensus <cit.>. In general, semi-analytical models <cit.> predict a weak mass dependence with a small, negative value for σ_/log M, but the magnitude of the overall scatter is typically much larger than what is inferred observationally <cit.>. Hydrodynamical simulations of galaxy formation, typically predict a significantly lower scatter, at least for haloes with M  10^12, in much better agreement with observations. Typically, though, they predict that the scatter rapidly increases for M  10^12 <cit.>. Finally, empirical models such as the  <cit.> and  <cit.> seem to predict σ_(M) relations that fall roughly in between the predictions from semi-analytical models and hydrodynamical simulations.

To conclude, we have demonstrated that satellite kinematics extracted from galaxy redshift surveys contain a wealth of information regarding the statistical relation between galaxies and their associated dark matter haloes. The Bayesian hierarchical framework , developed here and in , is able to analyse such data in an unbiased way, yielding accurate constraints on the galaxy-halo connection over a wide range of halo mass, and with unprecedented precision. Hence, satellite kinematics is complementary to other techniques that are used to constrain the galaxy-halo connection, in particular galaxy clustering and galaxy-galaxy lensing. Importantly, by only probing the smallest, most non-linear scales (i.e., the 1-halo term) it is insensitive to halo assembly bias, which hampers an unambiguous interpretation of the 2-halo term in clustering and galaxy-galaxy lensing. Hence, it is to be expected that by combining all these methods, degeneracies can be broken which opens up new avenues to test our cosmological paradigm and our models for galaxy formation. To this end, we plan to use, and where necessary further develop, in future work. In particular, among others, we intend to explore additional degrees of freedom in the characterization of the galaxy-halo connection (for example, mass dependence in the faint-end slope of the satellite CLF and in the ratio of the characteristic luminosities of the centrals and satellites in haloes of any given mass), the impact of baryonic effects on the halo potential (which may introduce systematic errors in the inference from satellite kinematics), and the impact of scatter in the halo concentration-mass relation (and the expected correlation with the abundance of subhaloes/satellite galaxies). In addition, we are excited about the prospects of using to probe the galaxy-halo connection as a function of secondary galaxy properties, such as galaxy colour and/or size, and to put constraints on cosmological parameters by combining satellite kinematics with other observables. 







§ ACKNOWLEDGMENTS


We are grateful to the anonymous referee for an insightful referee report that has resulted in significant improvements of the manuscript. FvdB has been supported by the National Aeronautics and Space Administration through Grant No. 19-ATP19-0059 issued as part of the Astrophysics Theory Program and by the National Science Foundation (NSF) through grant AST-2307280, and received additional support from the Klaus Tschira foundation. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611, and at the KITP Santa Barbara, which is supported in part by the National Science Foundation under Grant No. NSF PHY-174895.

This work utilized, primarily for plotting purposes, the following python packages:  <cit.>,  <cit.>,  <cit.>, and  <cit.>.







§ DATA AVAILABILITY


The SDSS DR7 New York University Value Added Galaxy Catalog  sample, used in this work, is publicly available at <http://sdss.physics.nyu.edu/lss/dr72/bright/0/>. All derived data, such as kinematic data of primaries and secondaries, and the data corresponding to each plot, will be made available on reasonable request to the corresponding author.





mnras










§ OPTIMIZING SELECTION CRITERIA






As described in <ref>, the selection of primaries and secondaries makes use of conical selection volumes. In this appendix we describe modifications to the parameters characterizing these selection volumes that result in an improved purity of the sample. Following <cit.>, <cit.> and <cit.> the selection cones are characterized by  = a_ σ_200, = a_ σ_200,  = b_ σ_200, and = b_ σ_200 (see Fig. <ref>). Here σ_200 is an estimate for the satellite velocity dispersion in units of 200, which scales with the luminosity of the primary as logσ_200 = c_0 + c_1 log L_10 + c_2 (log L_10)^2, where L_10 =  / (10^10). In we adopted exactly the same parameters as <cit.>: a_=0.5, a_ = 0.15, b_ = 1000, b_ = 4000/σ_200 and (c_0, c_1, c_2) = (-0.04, 0.38, 0.29). As discussed in , this results in impurity fractions of ∼ 5 percent. These impurities cause a slightly bias estimate of the scatter in the galaxy-halo connection. In the case of the small mock data samples used in , the effect was not significant. However, when using an 8-times larger, full-size SDSS sample the systematic bias in scatter becomes >3 σ significant.

In <ref> we classified impurities as either Type-I (BHG satellites) or Type-II (neither a central nor a BHG satellite). The green dots in the left-hand panel of Fig. <ref> show the luminosities as function of host halo mass of the primaries selected from the Tier-3 mock data using the selection criteria used in . Blue and orange contours mark the 4 and 7σ_ ranges around the median relation between halo mass and central luminosity used to construct the mock data. Red and black circled dots mark impurities of Type-I and II, respectively. As is evident, Type-I impurities have luminosities that are comparable to those of true centrals at the same host halo mass. That is because Type-I impurities are BHG satellites which are brighter than their corresponding central galaxy (hence they must have a luminosity in the typical range of Φ_(L|M) for a true central to be fainter). Being the brightest one in the corresponding halo, a Type-I impurity is impossible to avoid in the selection procedure. However, as we have demonstrated in <ref>, we actually forward model the contribution of Type-I impurities.

Type-II impurities, though, are a much bigger concern. As is evident from Fig. <ref>, these can have luminosities that are much lower than that of a typical central at the corresponding halo mass (by as much as 7σ). Since lower luminosities are indicative of a lower halo mass, a too-large contribution of Type-II impurities can give rise to significant, systematic errors in the inference. Since the kinematic information from secondaries associate with Type-II impurities still reflect a high velocity dispersion consistent with the actual halo mass, the main effect of Type-II impurities is to cause a systematic overestimate in the scatter of central luminosities at a fixed halo mass (i.e., an overestimate of σ_). Since we are not aware of a reliable method to forward model the impact of Type-II impurities, it is prudent that we minimize their incidence by tuning our selection criteria accordingly.

After extensive testing with different mock data sets similar to the Tier-3 mock discussed in the main text, we finally settled on the following set of parameters: a_ = 0.6, a_ =0.15, b_ = b_ = 1000 and {c_0, c_1, c_2} = {0.04, 0.48, 0.05}. With these new selection criteria we are able to reduce the fraction of Type-II impurities from ∼ 1.1 % to ∼ 0.5 %. The impact of this reduction can be seen by comparing the two panels of Fig. <ref>. Note that in addition to dropping the fraction of Type-II impurities to sub-percent levels, the new selection criteria preferentially removes the most dramatic outliers and also drastically reduces the contribution of Type-II impurities at the high mass end, where the old selection criteria caused the fraction of Type-II impurities to be very high.

With these new and improved selection criteria we find that Type-II impurities no longer cause a significant overestimate of σ_. Although the new selection criteria reduces the number of primaries in the satellite kinematics sample by almost 40 percent, we find that this does not significantly compromise the precision with which can infer the galaxy-halo connection. The reason is that the main reduction of primaries occurs at the low-luminosity end, where most of the secondaries are interlopers that do little to constrain the halo occupation model. 



§ THE COMPLETENESS OF CENTRALS






The selection of centrals as primaries (<ref>) is not complete; i.e., not every central is selected as a primary. This incompleteness owes to two different reasons: (i) incompleteness in the SDSS redshift survey, for example due to fibre-collisions, or (ii) the central is located inside the selection cone of a brighter galaxy.  Let the completeness (M,L,z) be the fraction of centrals — of luminosity L, at redshift z, residing in haloes of mass M, in the survey volume of the SDSS — that are selected as primaries. We can write that (M,L,z) = (M|L,z) _0(L,z). As discussed in the main text (see <ref>), the modelling in is independent of _0, which drops out. In other words, we only need to account for any potential halo mass dependence of the completeness given by (M|L,z).

In order to gauge this mass dependence, we construct 100  mock SDSS redshift surveys similar to the Tier-3 mock survey discussed in the main text to which we apply our primary selection criteria. For each central galaxy in the mock SDSS volumes we assess whether it is selected as a primary. The combined results from all 100 mocks are shown as symbols with Poisson errorbars in Fig. <ref>. Different colours correspond to different luminosities of the centrals, as indicated. Here we have combined data on all centrals over the entire redshift range, but we emphasize that the redshift dependence is weak. A few trends are evident. First of all, the completeness is lower for fainter centrals. This simply reflects that fainter centrals are more likely to have a brighter galaxy in a neighbouring halo that happens to fall within the primary selection criterion. Modelling this would require accurate knowledge of the clustering of haloes (2-halo term) and is sensitive to assembly bias issues. Fortunately, we do not need to model this luminosity dependence. All we care about is the halo mass dependence as characterized by (M|L,z). 

As is evident from the data, the completeness for centrals of a given luminosity is roughly independent of halo mass at the low-mass end, but then drops drastically at the high-mass end. This transition occurs at higher mass for brighter centrals. The yellow-shaded region indicates, for each luminosity bin, the 5 to 95 percentile range of halo masses. Note that the exponential decline of (M|L,z) at large masses only affects the top ∼ 5 percent of centrals of a given luminosity. This is why, in , we decided to ignore this mass dependence all together. However, as we demonstrate below, it is actually fairly straightforward to model (M|L,z). As it turns out, this mass dependence owes almost entirely to the fact that a central of a given luminosity in a more massive halo is more likely to have a satellite that is brighter than itself. Recall that Φ_(L|M) is modelled as a log-normal distribution. Hence, the central galaxies in haloes of a given mass have a tail of excessively faint centrals. And since we assume that the luminosities of satellite galaxies are independent of that of their central, those faint centrals are more likely to have a brighter satellite, and thus to fail selection as a primary. As shown in Appendix A of <cit.>, the probability that the central galaxy is the brightest galaxy in its halo is given by


    P_ BC(M|L,z) = exp[-Λ_ BTC] .


Here Λ_ BTC is the expectation value for the number of satellites that are brighter than the central, which is given by equation (<ref>).

The solid lines in Fig. <ref> show the predictions for P_ BC(M|L,z), computed using equation (<ref>) for the same CLF model as used to construct the mocks. The absolute normalization for each luminosity bin, which represents _0(L,z), is tuned to match the mock data at the low-log M end. As is evident, equation (<ref>) accurately describes the halo-mass dependence of the completeness of primaries. Hence, the mass-dependence of the completeness of centrals can be modelled as (M|L,z) = P_ BC(M|L,z). 

Note that 1-(M|L,z) is the probability that a central is not the brightest galaxy in its halo, and thus the probability that the halo gives rise to a Type-I impurity. As discussed in the main text, we forward model the contribution of these Type-I impurities, which effectively means that we already account for the mass-dependence of the completeness of centrals depicted in Fig. <ref>. Indeed, P_ BC given by equation (<ref>) is identical to P(L_ bs<L|M,z) (equation [<ref>]) used in <ref> to forward model the Type-I impurities. 



§ THE RADIAL PROFILE OF SATELLITES







Throughout we assume that n_ sat(r|M) is characterized by a generalized-NFW form (equation <ref>) which has two free parameters: the inner logarithmic density slope, γ, and the concentration ratio ℛ = c_ vir / c_ sat which characterizes the scale radius of the number density profile. As discussed in <ref>, pre-computes and stores essential arrays which are then used in each step of the MCMC chain varying the CLF, anisotropy, and nuisance parameters. This pre-computation drastically speeds up , but requires that n_ sat(r|M), and thus {γ,}, are held fixed. Therefore, instead of keeping the radial profile free, we run separate MCMC chains for each assumed radial profile on a 15 × 15 grid of {γ, log}. We then combine the posteriors and likelihoods from each of these runs to compute the marginalized likelihood (γ,  | ).

The grey contours in the top-left panel of Fig. <ref> show the 68, 95 and 99 percent confidence intervals for γ and log thus obtained using the SDSS data described in <ref>. The best-fit values, indicated by the black cross, correspond to {γ, } = {0.94,  1.7}, which are the values we adopt for our detailed SDSS-analysis described in <ref>. However, the confidence intervals for γ and log reveal a significant degeneracy along a narrow ridge-line in γ- parameter space (see also ). To demonstrate the impact this degeneracy has on our inference, the coloured histograms in the other panels of Fig. <ref> show the posteriors on our CLF parameters and the anisotropy parameter, β, for 9 different combinations of γ and  (indicated by the circles of corresponding colour in the top-left panel) that roughly trace out the boundary of the 3σ confidence interval. The vertical black dashed line in each of these panels show the best-fit parameters inferred by with the best-fit {γ,  } combination, same as the values quoted in table <ref>.  As is evident, the inferred CLF parameters are extremely robust to changes in γ and  along the direction of this degeneracy. The only parameter that shows a weak dependence is the orbital anisotropy parameter β (bottom-right panel), which is to be expected from the fact that both n_ sat(r|M) and β appear in the expression for the line-of-sight velocity dispersion given by equation (<ref>). 

Note also that the constraints on γ and  are inconsistent with satellite galaxies following the same radial profile as the dark matter (i.e., γ =  = 1, indicated by the red square in the top-left panel) at >5σ significance. If we had assumed that satellite galaxies are an unbiased tracer of their host halo mass distribution, which is not uncommon in the literature when modelling the galaxy-halo connection, we would have obtained the posteriors indicated by the red histograms. Interestingly, most CLF parameters would still be consistent with the values inferred using our fiducial, best-fit model with {γ, } = {0.94,  1.7}. The main exceptions, though, is the orbital anisotropy parameter, which would be biased high (i.e., we would markedly overestimate the radial velocity anisotropy). A few other parameters like γ_2, α_, and  σ_, are also somewhat biased in red histograms. Thus, we would incorrectly infer a steeper L̅_ (M) relation, a shallower faint-end slope of satellite CLF, and a slight decreasing trend in the central luminosity scatter with host halo mass, if we wrongly assumed the satellites to follow the dark matter radial distribution. In conclusion, yields tight constraints on the radial number density profile of satellite galaxies, and whatever degeneracy remains between the central density slope and concentration has no significant impact on any inferred parameter. 




