An investigation of halo and galaxy assembly bias

Observing galaxy clustering is an effective way to gain knowledge of galaxy formation and evolution and to constrain cosmology. Cosmology determines halo population and clustering (bias), and galaxy-halo relation connects halo clustering to galaxy clustering. However, traditional halobased models of galaxy clustering only consider the halo mass dependence and ignore any effect caused by halo assembly history, which may cause systematics in constraining the galaxy-halo relation and cosmology. My work is focused on investigating halo and galaxy assembly bias, aiming at providing insights to improve models of galaxy clustering. First, the joint dependence of halo bias on halo mass and other halo assembly variables is investigated. It is found that halo bias increases outward from a global minimum on the plane of halo mass and one assembly variable. Based on the joint dependence, halo assembly bias is unlikely to be fully absorbed by any combination of halo variables. Furthermore, dependence of halo bias on one assembly variable can be independent from that on another one even if the two variables are correlated. Additionally, assembly effects on halo kinematics are shown to correlate with those of spatial clustering. Another component in halo models, the galaxy-halo relation, is studied with a hydrodynamic simulation. The correlations between a set of galaxy properties and halo properties are analyzed. Then a simple model is developed to relate galaxy assembly bias and halo assembly bias through the correlation coefficient between halo and galaxy properties, which provides an effective prescription to incorporate galaxy assembly bias into halo models of galaxy clustering. Finally, a specific aspect of the galaxy-halo relation, the distribution of galaxy luminosity and color at a given halo mass (i.e. the conditional color magnitude distribution, CCMD) in galaxy formation models is investigated. In a semi-analytical galaxy formation model, the CCMD of central galaxies is found to be comprised of a red and a blue population, in a pattern consistent with that inferred from modeling the observed luminosity and color-dependent galaxy clustering. An attempt has been made to map galaxy color and luminosity onto halo assembly variables to figure out the physical origin of the two populations.

AN INVESTIGATION OF HALO AND GALAXY ASSEMBLY BIAS by Xiaoju Xu A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Physics Department of Physics and Astronomy The University of Utah August 2019 Copyright © Xiaoju Xu 2019 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of Xiaoju Xu has been approved by the following supervisory committee members: Zheng Zheng , Chair(s) 13 May 2019 Date Approved Kyle Dawson , Member 13 May 2019 Date Approved Benjamin Bromley , Member 13 May 2019 Date Approved Pearl Sandick , Member 13 May 2019 Date Approved Jingyi Zhu , Member 14 May 2019 Date Approved by Peter E. Trapa , Chair/Dean of the Department/College/School of Physics and Astronomy and by David B. Kieda , Dean of The Graduate School. ABSTRACT Observing galaxy clustering is an effective way to gain knowledge of galaxy formation and evolution and to constrain cosmology. Cosmology determines halo population and clustering (bias), and galaxy-halo relation connects halo clustering to galaxy clustering. However, traditional halobased models of galaxy clustering only consider the halo mass dependence and ignore any effect caused by halo assembly history, which may cause systematics in constraining the galaxy-halo relation and cosmology. My work is focused on investigating halo and galaxy assembly bias, aiming at providing insights to improve models of galaxy clustering. First, the joint dependence of halo bias on halo mass and other halo assembly variables is investigated. It is found that halo bias increases outward from a global minimum on the plane of halo mass and one assembly variable. Based on the joint dependence, halo assembly bias is unlikely to be fully absorbed by any combination of halo variables. Furthermore, dependence of halo bias on one assembly variable can be independent from that on another one even if the two variables are correlated. Additionally, assembly effects on halo kinematics are shown to correlate with those of spatial clustering. Another component in halo models, the galaxy-halo relation, is studied with a hydrodynamic simulation. The correlations between a set of galaxy properties and halo properties are analyzed. Then a simple model is developed to relate galaxy assembly bias and halo assembly bias through the correlation coefficient between halo and galaxy properties, which provides an effective prescription to incorporate galaxy assembly bias into halo models of galaxy clustering. Finally, a specific aspect of the galaxy-halo relation, the distribution of galaxy luminosity and color at a given halo mass (i.e. the conditional color magnitude distribution, CCMD) in galaxy formation models is investigated. In a semi-analytical galaxy formation model, the CCMD of central galaxies is found to be comprised of a red and a blue population, in a pattern consistent with that inferred from modeling the observed luminosity and color-dependent galaxy clustering. An attempt has been made to map galaxy color and luminosity onto halo assembly variables to figure out the physical origin of the two populations. For my parents. CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi CHAPTERS 1. 2. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Cosmology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Dark matter halo formation and clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Press-Schechter model and excursion set theory . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Halo clustering and assembly bias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Defining halo in simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Galaxy-halo relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Statistical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Physical galaxy formation models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 2 4 5 5 6 7 DEPENDENCE OF HALO BIAS AND KINEMATICS ON ASSEMBLY VARIABLES 9 2.1 2.2 2.3 2.4 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation and halo bias calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Joint dependence of halo bias on mass and assembly variables . . . . . . . . . . . . . . . . . 2.4.1 Dependence of halo bias on mass and one assembly variable . . . . . . . . . . . . . . 2.4.2 Effective mass to absorb the assembly bias effect? . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Scale dependence of halo bias with assembly variable . . . . . . . . . . . . . . . . . . . 2.4.4 Dependence of halo bias on two assembly variables . . . . . . . . . . . . . . . . . . . . 2.5 Dependence of pairwise velocity and velocity dispersions on halo assembly . . . . . . . 2.5.1 Assembly effect on pairwise velocity of haloes . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Relation between spatial clustering and pairwise velocity under assembly effect 2.6 Summary and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. 9 10 12 13 13 16 21 26 28 28 34 37 GALAXY ASSEMBLY BIAS OF CENTRAL GALAXIES IN THE ILLUSTRIS SIMULATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1 3.2 3.3 3.4 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation and galaxy-halo catalogue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Relationship between stellar mass and halo properties . . . . . . . . . . . . . . . . . . . 3.4.2 Relationship between galaxy properties and halo properties . . . . . . . . . . . . . . . 3.4.2.1 At fixed Mh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2.2 At fixed Vpeak . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 40 41 42 43 46 46 49 3.4.2.3 Dependence on Mh and Vpeak . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4.3 Assembly bias of central galaxies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.5 Summary and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4. CONDITIONAL COLOR MAGNITUDE DISTRIBUTION FROM GALAXY FORMATION MODELS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 CCMD in hydrodynamic simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 EAGLE simulation and Gaussian mixture model . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 CCMD in EAGLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 CCMD from a semi-analytical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Galaxy formation model and CCMD components . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Relation between CCMD components and halo properties . . . . . . . . . . . . . . . . 5. 62 63 63 64 66 66 67 SUMMARY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 APPENDIX: MASS RESOLUTION EFFECT ON HALO BIAS AND PAIRWISE VELOCITY STATISTICS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 vi LIST OF FIGURES 2.1 Joint dependence of halo bias on halo mass and each assembly variable. The assembly variables in the four panels are peak maximum circular velocity Vpeak , halo formation scale factor aM/2 , halo concentration c, and halo spin λ, respectively. The four black curves in each panel mark the central 50 and 80 percent of the distribution of the corresponding assembly variable as a function of halo mass. . . . . . . . . . . . . . . . . . . . . 14 2.2 Effective halo mass to absorb assembly bias effect with halo spin. The top panel is the same as the bottom-right panel of Fig. 2.1, and dashed lines are added to illustrate how the effective mass is defined (see text). In the bottom panel, the black points show the halo bias as a function of halo mass, with the scatter at fixed halo mass coming from haloes with different spin parameters (assembly bias effect with spin). The red points show the dependence of halo bias on effective halo mass, and the scatter caused by the assembly effect with spin is substantially reduced. . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Same as Fig. 2.1, but with original halo mass replaced by effective halo mass. The effective halo mass absorbs the halo assembly bias with halo spin, manifested by the vertical contours in the bottom-right panel. However, assembly bias effect still exists with other assembly variables, as seen in the other three panels. . . . . . . . . . . . . . . . . . 19 2.4 Comparison between assembly bias effects under original halo mass (top) and effective halo mass (bottom). At fixed original or effective halo mass, haloes are grouped according to the percentiles of each assembly variable, and their bias values normalised by the mass-only dependent bias are shown as curves in each panel (dashed for the 10–25 percent percentile and solid for the 75–90 percent percentile). . . . . . . . . . . . . . 20 2.5 Scale dependence of assembly bias. The upper subplot shows the scale-dependent halo bias factor for different assembly variables. In each panel, solid curves are calculated through halo-mass cross-correlation functions and dotted curves through halo-halo auto-correlation functions. Each set of five curves are for haloes in different percentiles of each assembly variable: lower 25 percent (blue), lower-middle 25 percent (cyan), upper-middle 25 percent (green), upper 25 percent (red), and all haloes in the mass bin (black). Top and bottom panels are for haloes of log[Mh /(h−1 M )]=12 and 14, respectively. Bottom subplot is similar, but each halo bias curve is normalised by the corresponding large-scale value (average on scales larger than 10h−1 Mpc). . . . . . . . . . 22 2.6 Halo bias as a function of two assembly variables for haloes in the mass range of log[Mh /(h−1 M )] = 12.5 ± 0.5. In each panel, the contours mark the 68 percent and 95 percent distribution of the two assembly variables and the color scale shows the value of halo bias. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.7 Same as Fig. 2.6, but for haloes in the mass range of log[Mh /(h−1 M )] = 13.5 ± 0.5. . 25 2.8 Pairwise velocity and velocity dispersions of haloes. The left panel shows the scaledependent pairwise radial velocity as a function of halo mass. The middle and right panels are similar, but for pairwise radial and transverse velocity dispersions, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.9 Dependence of halo pairwise velocity and velocity dispersions on assembly variables for low-mass haloes (log[Mh /(h−1 M )] = 11.73). Right, middle, and light panels show pairwise radial velocity, pairwise radial velocity dispersion, and pairwise transverse velocity dispersion, respectively. In each panel, the curves are color-coded according to the value of the corresponding assembly variable. . . . . . . . . . . . . . . . . . . 31 2.10 Same as Fig. 2.9, but in each panel, the curves are color-coded according to the percentile of the corresponding assembly variable. . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.11 Dependence of halo pairwise velocity and velocity dispersions on halo mass and each of the assembly variables. left, middle, and right panels show pairwise radial velocity, pairwise radial velocity dispersion, and pairwise transverse velocity dispersion, respectively. The pairwise velocity and velocity dispersions are evaluated at r ∼ 10h−1 Mpc. The four black curves in each panel mark the central 50 and 80 percent of the distribution of the corresponding assembly variable as a function of halo mass. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.12 Relation between pairwise velocity and velocity dispersions and halo bias. Left, middle, and right panels show pairwise radial velocity, pairwise radial velocity dispersion, and pairwise transverse velocity dispersion versus halo bias, respectively. In each panel, the points are color-coded by halo mass and their size indicates the value of the corresponding assembly variable. The pairwise velocity and velocity dispersions are evaluated at r ∼ 10h−1 Mpc. The solid curve is computed from mass-only dependence of halo clustering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.1 Top-left: M∗ as function of Mh for central galaxies. The galaxies are color-coded according to log[Vpeak /(km s−1 )]. For galaxies in each bin of log Vpeak , the contours correspond to the 68.3 and 95.4 percent distribution, respectively. Top-right: M∗ as function of Vpeak for central galaxies, color-coded according to log[Mh /(h−1 M )]. Bottom-left: M∗ as function of Mh for central galaxies, with the mean relation color-coded according to the values of aM/2 . For clarity, the scatter in the mean relation is only shown for the bin with the highest aM/2 (latest forming haloes). Bottom-right: M∗ as function of Vpeak for central galaxies, color-coded according to aM/2 , with the shaded region illustrating the scatter for the bin with the highest aM/2 . Note the remarkable result that M∗ does not depends on aM/2 at fixed Vpeak (compared to the Mh case in the bottom-left panel). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2 Standard deviation in log M∗ as a function of Mh (solid) and Vpeak (dashed). The correspondence between Mh and Vpeak is from the mean relation Vpeak ∝ Mh1/3 in equation (3.1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 viii 3.3 Relation between each pair of galaxy and/or halo properties at log[Mh /(h−1 M )] ∼ 12. In each contour panel, the two contours show the central 68.3 and 95.4 percent of the distribution of the pair of properties. The panels with red contours (i.e. the top 6 rows and left 6 columns of contour panels) display the correlations between halo properties. Those with black contours (i.e. the bottom 4 rows and the left 7 columns of contour panels) show the correlations between galaxy and halo properties, and those with blue contours (i.e. the right 3 columns of contour panels) are for the correlations between pairs of galaxy properties. The number in each contour panel is the Pearson correlation coefficient for the pair of properties. The histogram at the top panel of each column is the probability distribution function of the variable of that column. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4 Same as Fig. 3.3, but at fixed log[Vpeak /(km s−1 )] ∼ 2.23. Note particularly the lack of correlation of M∗ with other assembly variables (including c, aM/2 , λ, MÛ h , MÛ h /Mh ), in contrast with the case in Fig. 3.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5 Pearson correlation coefficient ρ of each pair of galaxy and/or halo properties as a function of Mh (solid) and Vpeak (dashed). The galaxy and halo properties are marked to the far left of each row and at the bottom of each column. For clarity, we only label the values of Vpeak on the horizontal axis, and the values of Mh can 3 be inferred from Mh ∝ Vpeak from equation (3.1). As with Fig. 3.3, panels with red, black, and blue curves are for correlations between halo-halo, galaxy-halo, and galaxy-galaxy properties, respectively. In each panel, the dotted horizontal line indicates no correlation. Note particularly the lack of correlation of M∗ with other assembly variables (including c, aM/2 , λ, MÛ h , MÛ h /Mh ) for the Vpeak dependence case, in contrast with the Mh -dependent case. Also the correlations of SFR, sSFR, and color with halo assembly variables show more consistent behaviours in the Vpeak dependence case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.6 Illustration of the correlation between galaxy property and halo assembly property and the construction of the halo and galaxy samples for the study of galaxy assembly bias effect in Section 3.4.3. In the top panel, the ellipse denotes the joint distribution of halo property x and galaxy property y at fixed Mh or Vpeak , which is assumed to follow a 2D Gaussian distribution. A halo sample is constructed with haloes of the top f fraction of x (indicated by the red region), and a galaxy sample is constructed with central galaxies of the top f fraction of y (indicated by the region inside the dashed curve). Shown in the bottom panel are the probability distribution functions of halo property x for all the haloes at fixed Mh or Vpeak (black solid+dashed), the selected haloes (dashed), and the selected galaxies (red solid), respectively. . . . . . . . . 55 ix 3.7 Connection between halo and galaxy assembly effect. In the left panels, the quantities shown are the values of b(A, Mh )/b(Mh ) − 1 of different samples. For each sample, the quantity is the fractional difference of the bias factor of the sample selected based on property A with respect to that of all haloes at fixed Mh (i.e. δhb or δgb defined in Section 3.4.3), which is used to characterise the magnitude of the assembly bias effect. The thick dotted curves are for halo samples selected based on halo formation time, and the solid curves are for central galaxy samples selected based on color (top), SFR (middle), and sSFR (bottom). The thin dashed lines are the predictions from the simple model presented in Section 3.4.3 according to the correlation between galaxy and halo properties (δgb = ρδhb , with ρ the correlation coefficient). For clarity, jackknife error bars are only shown for the solid curves in each panel. The right panels are the same, but for the assembly bias effect as a function of Vpeak . See details in Section 3.4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1 Color magnitude distribution of galaxies at fixed halo mass bins from EAGLE. In each panel, black dots are galaxies, and red/blue contours are the two components of the Gaussian mixture model that fits data the best. Halo masses are at the bottom-right, in the form of log[Mh /(h−1 M )]. Except for the last three mass bins, fitted red and blue components cannot be separated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.2 AIC and BIC of Gaussian mixture models from the EAGLE simulation. In each panel, y-axis is the AIC/BIC estimator, and x-axis is the number of Gaussian components in the model. Solid line represents AIC, and dashed represents BIC. The smaller the estimator value, the better the fitting model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3 Same as Fig. 4.1, but with Guo2011 semi-analytical galaxy formation model. . . . . . . 68 4.4 Joint distribution of half mass scale and concentration of extreme red and blue components. Extreme color components are the red and blue galaxies found in Section 3.4.2, with a higher than 0.8 probability to belong to one of the color component. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.5 Distribution of galaxy major merger variables on main branch. Top left: loga distribution of each major merger of each red/blue galaxy along its main branch; the number on y-axis indicates each galaxy. Top middle: histogram of loga of first major merger of red/blue galaxies. Top right: histogram of loga of last major merger of red/blue galaxies. Bottom left: distribution of total number of major merger against loga of first major merger of red/blue galaxies. Bottom middle: distribution of total number of major merger against loga of last major merger. Bottom right: loga of last major merger against loga of first major merger. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 A.1 Same as Fig. 2.1, but the Bolshoi simulation is used for log[Mh /(h−1 M )] < 12.8 and MDR1 simulation for higher mass. The two vertical dashed lines in the each panel indicate the masses of haloes with 400 particles in the two simulations. . . . . . . . . . . . 75 A.2 Same as Fig. 2.10, but haloes from the Bolshoi simulation are used for the analysis. . . 76 A.3 Same as Fig. 2.11, but the Bolshoi simulation is used for log[Mh /(h−1 M )] < 12.8 and MDR1 simulation for higher mass. The two vertical dashed lines in the each panel indicate the masses of haloes with 400 particles in the two simulations. . . . . . . 78 x ACKNOWLEDGEMENTS I would like to thank my advisor Prof. Zheng Zheng. When I joined Zheng’s group, I had little understanding of both the field of study and the way scientific research works. Zheng patiently introduced the big picture of the research area, and showed me interesting works that we can do to contribute to the field. I learned a lot from his very insightful ideas. As an excellent advisor, Zheng paid more attention to train students to understand the underlying basic picture instead of specific knowledge. I benefited a lot from the order-of-maginitude problems in group meetings. Additionally, Zheng always gives great suggestions on scientific writing and presentation preparing, which make our work more understandable and impressive. I also thank Zheng for supporting me to attend conferences and encouraging me to interact with people from other places. I wish I could have joined the group earlier and had more time to learn from Zheng. I thank all our group members, Shiyu Nie, Kevin McCarthy, Li Yang, and previous members Haojie Xu and Greg Eng, for explaining their research and discussing my work. I also want to thank Prof. Kyle Dawson from whom I get suggestions for my work, and his group members, Vikrant Kamble and Hélion du Mas des Bourboux, from whom I learned a lot of knowledge of different fields. Last but not the least, I want to thank my parents and my husband for their love, care, and encouragement during my PhD study. Without their support, my PhD life could have been much harder. CHAPTER 1 INTRODUCTION Based on General Relativity, modern cosmology has been developed to enlarge our knowledge of the universe. The Friedmann equation describes the expansion of the universe with all the matter and energy components in it. The universe starts with a big bang, experiences a quick inflation at early age, and expands to what it looks like today. The ΛCDM (Λ-cold dark matter, cosmological constant dark energy and non-relativistic dark matter) model, which states that the universe is dominated by cold dark matter (CDM) and constant density dark energy (Λ), is the most successful model that can describe the expansion history after inflation. In this chapter, I will review the ΛCDM cosmology model briefly and discuss some tools that can be used to study cosmology. 1.1 Cosmology General Relativity says that the space-time curvature is connected to the energy distribution, and gravity is a result of space-time geometry changing from flat. Space tells matter how to move, and matter tells space how to curve. This is exactly the famous Einstein field equation: 1 8πG Rµν − Rgµν = 4 Tµν . 2 c (1.1) The left-hand side of the equation describes space-time curvature. Rµν is the Ricci tensor and R is scalar curvature defined by Ricci tensor. On the right-hand side, the energy-momentum tensor Tµν includes all the mass and energy components in space. Once all the components (e.g., baryonic matter, dark matter, dark energy) are considered, it predicts expansion of the universe and evolution of energy density (e.g., Friedmann equations). Under the principle that the universe is isotropic and homogeneous on large scales, the Friedmann–Robertson–Walker (FRW) metric can be used in the field equation to derive Friedmann equations: aÛ 2 kc2 8πG + 2 2 = ρ a 3 a R0 (1.2) In this equation, a is the scale factor of the universe, ρ is spatial density (includes matter, radiation, and Λ), k stands for spatial curvature, and R0 is the curvature radius at z=0. This equation describes 2 the expansion rate of the universe in terms of density and curvature. In the early universe, the dominant component is radiation, but density decays very quickly with expansion of the universe. Matter density decays more slowly than radiation, so it takes the stage after radiation. While matter density decreases with expansion, density of dark energy doesn’t decrease in the ΛCDM model. Therefore, dark energy gradually becomes the dominant component in the late universe. Structure forms and evolves under the competition of gravity and expansion. For example, dark matter will collapse into virialized objects called dark matter halos and baryonic matter condenses to form galaxies in the potential well of dark matter halos. Dark matter halos move according to the underlying gravitational field along the history of the universe, and produce a specific pattern of distribution. Galaxy clustering is similar to halo clustering on large scales. Observing this specific galaxy clustering pattern could be an effective way to gain knowledge of cosmology. 1.2 Dark matter halo formation and clustering To study galaxy clustering, it is natural to start with halo clustering, since galaxies form in dark matter halos. A successful halo formation model is supposed to predict the spatial density of halos, as well as spatial distribution of halos. Then the galaxy occupation model can be used to transfer halo clustering into galaxy clustering. The halo occupation model can be constrained with observational galaxy clustering, as long as they are working together with halo clustering predicted from certain cosmological parameters. Another possibility is that both halo occupation and cosmological parameters can be constrained at the same time with appropriate galaxy statistical observables. In this section, I will introduce analytical halo clustering models and galaxy-halo relation models. 1.2.1 Press-Schechter model and excursion set theory The initial condition of matter density is considered to be a Gaussian random field. In this field, overdense regions collapse under gravity and form dark matter halos. At first, volume of overdensity increases together with expansion of the universe, and then turns around to contract under its self gravity. Instead of contracts to a high density point, it will oscillate and finally stop at half the maximum radius to become a virialized (e.g., kinetic energy is half the potential energy) self-bound system. Theoretically, any region that has density contrast δ = ρ−ρ̄ ρ̄ > 0 (ρ is matter density and ρ̄ is average density) can collapse to form halos eventually. However, we can at most observe halos 3 formed at or before z=0, so there exists a threshold that the overdensity must exceed to collapse at some redshift. In the spherical collapse model, this density threshold is δsc = 1.686 by linearly extrapolating the initial overdensity to the redshift of collapse. The purpose of the analytical model of halo formation is connecting halo population at a later time to density distribution in the initial field, and deriving halo mass function which describes the number density of the halo at a fixed halo mass. The Press-Schechter model (Press & Schechter (1974)) demonstrates that the fraction of mass ended up in halos with mass greater than M equal to twice the fraction of the initial density field that exceeds the density threshold of spherical collapse on scale R. M and R are related by the mean density of the background universe. The factor of two is introduced to account for the situation that a region reaches the collapse threshold on a larger scale but under the threshold on the scale considered. The excursion set theory (EST, Bond et al. (1991)) smooths the overdensity field around a point with filters of different scales. The density series under these filters can be described by random walk. The distribution of the first upcrossing of the collapse threshold by random walk trajectory as a function of smoothing scales allows the inference of halo mass function. If we use a sharp k filter, there will be no correlation between random walk steps (Markov chain random walk). This can largely simplify the model. Same as the Press-Schechter model, the fraction of the mass that ended up in halos with mass greater than M is still linked to the fraction of the initial density field that exceeds the collapse threshold. But in EST, it only accounts for the first cross of the threshold if we consider the density as a function of different radius (or mass included in that radius, or equivalently the root mean square of mass enclosed in that radius). So it accounts for the missing factor of two in the Press-Schechter model. With the first crossing distribution found, the halo mass function at a given mass can also be established. A more complicated theory can also be developed when the random walk takes correlated steps (Maggiore & Riotto (2010); Paranjape et al. (2012)). Since the initial density field is in linear regime, density fluctuations on large scales cannot be big, so they require more time to reach the collapse threshold. On the other hand, fluctuations on small scales can be bigger and collapse earlier. So the low mass halos collapse first, and massive halos form later by merger of low mass halos. In EST framework, merger rates and accretion rates can be predicted. However, excursion set theory considers any arbitrary point in space, which may not be the center of a halo at later time. Excursion set peak (Hahn & Paranjape (2013); Paranjape et al. (2013)) can compensate for that by considering only peaks in initial density field, and is thus more consistent 4 with simulations. Both Press-Schechter and excursion set theory show discrepancies with N-body simulations on the halo mass function, by underpredicting less massive halos and overpredicting low mass halos. The reason for discrepancies is that both of them are in the framework of spherical collapse. The discrepancies can be alleviated by extending the spherical collapse model to the elliptical collapse model (Sheth & Tormen (2002)), where the threshold for collapse varies with scale (a moving barrier). The barrier for collapse is identical to that in spherical collapse for larger collapse scale (massive halo). It means that massive halos have more spherical shape than low mass halos. On the other hand, collapse barrier for smaller halos is higher than that in spherical collapse, so it suppresses the abundance of the low mass halo to be more consistent with N-body simulations. 1.2.2 Halo clustering and assembly bias Analytical halo formation models above are capable of predicting halo mass function, as well as the halo bias. The halo bias describes how dark matter halo distribution biases from the underlying matter distribution, ξhh (r) = b2 ξmm (r). ξhh (r) is the halo-halo auto-correlation function, which is the excess probability of finding a pair of halos at separation r comparing to that of a random distribution. Similarly, ξmm (r) is the mass autocorrelation function. Under the same matter correlation function, massive halos have higher correlation amplitude and halo bias, which means that their clustering is stronger than that of low mass halos. Mo & White (1996) derived the expression of halo bias as a function of peak height ν = δ/σ(M), with σ(M) the root mean square of the matter density field, and ν can be considered as a variable for halo mass. In this model, halo bias can be written as only a function of halo mass, without dependence on halo formation time and other variables related to assembly history, and this is consistent with the simulation they have. However, the dependence of halo bias factor on halo properties other than halo mass was discovered in new simulations (e.g., halo assembly bias, Gao et al. (2005); Gao & White (2007)). For example, at a fixed halo mass, low mass halos formed earlier have larger halo bias than those formed later. Low mass halos with higher concentration have a larger bias at a fixed halo mass, but for massive halos, this trend reverses. Halo assembly bias is a strong prediction of ΛCDM from N-body simulations, but it is not included in earlier halo formation models. Although there have been some efforts (Dalal et al. (2008); Castorina & Sheth (2013); Shi & Sheth (2017)) on explaining assembly bias theoretically from different approaches, it still 5 remains unclear at some level. Numerically, environment (matter overdensity or tidal anisotropy) is believed to be one of the most important causes of assembly bias, in the sense that halos in denser regions cluster more strongly (Han et al. (2019); Ramakrishnan et al. (2019)). Other variables also contribute, for example, halo shape, internal velocity anisotropy, and halo spin. Recently, Mansfield & Kravtsov (2019) demonstrated that splashback halos at the splashback radius of massive neighbors can also contribute to halo assembly bias. In Chapter 2, I will present a study of halo assembly bias, and also show the relation between halo kinematics and halo bias. 1.2.3 Defining halo in simulations One way to test theoretical predictions of halo statistics is to study them in N-body simulations. N-body simulation follows the evolution of a large number of dark matter particles under the gravitational influence. Dark matter particles are the test particles with a mass of 5-10 orders of magnitude of the solar mass. Particle mass reflects the mass resolution of the simulation; the lower particle mass, the higher resolution. N-body simulation starts with initial conditions that are set up by the power spectrum, and then it records information of particles at each time step. Finally, particles evolve into a non-linear structure. To identify individual halos with discrete particles, the Friend-of-Friend (FOF, Davis et al. (1985)) algorithm finds the neighbors of a particle within a linking length and links them together, then identifies them to be a halo if the average density of the region from these particles exceeds 178 or 200 times the average density in the simulation box. This is the most generally used method to identify halos in N-body simulations. However, this method has a disadvantage. If two halos are close, with a particle in between of them, it is possible that FOF algorithm will link these two halos to be one. To avoid this, in Rockstar (Behroozi et al. (2013)) halo finder, 6D phase space is used to link particles together instead of 3D space. FOF and Rockstar halo finder may not result in the same halo on the individual level, but the obtained halo mass functions are similar. 1.3 Galaxy-halo relation Although halo clustering is straightforward to study in N-body simulations, halos cannot be observed directly. One way to avoid this difficulty is to link galaxies to dark matter halos and predict galaxy population base on the halo population. However, assumptions need to be made to link them, for example, which galaxy properties depend on which halo properties. In this section, I introduce some methods to model galaxy-halo relations, from statistical models to physical models. In the 6 statistical models, there is no actual galaxy formation physics involved and galaxies are populated into dark matter halos with certain prescriptions. On the opposite side, physical models include some physical processes, through semi-analytical modeling or direct simulations. More details about galaxy-halo relations can be found in Wechsler & Tinker (2018). 1.3.1 Statistical models Abundance matching (Kravtsov et al. (2004); Reddick et al. (2013); Guo et al. (2016)) is one effective method to populate halos with galaxies. Galaxies are ranked based on galaxy properties like luminosity or stellar mass, and halos are ranked according to halo properties like halo mass or other mass indicators. Then the galaxy property and halo property is connected by matching the number densities of galaxies above a luminosity (or stellar mass) threshold and halos above a halo mass (or other indicators) threshold. For the central subhalo which is at the potential minimum of the host halo, the mass indicator can be halo mass itself. For satellite subhalos, their mass would be stripped once they are accreted onto the host halo, so the halo mass before accretion or the maximum circular velocity may be used in abundance matching as a mass indicator. Simple abundance matching models assume a monotonic relation between galaxy property and halo property, while others also model the scatter in this relation. An extended version is age matching, which matches galaxy color or star formation rate to halo age at fixed halo mass (Hearin & Watson (2013); Watson et al. (2015)). Halo occupation distribution (HOD; Berlind & Weinberg (2002); Zheng et al. (2005)) is another method to link galaxies and halos. It introduces the probability P(N| M) that N galaxies of a given type can be found in a halo of mass M. Each halo hosts one central galaxy at the halo center, and hosts one or more (or zero) satellite galaxies according to the host halo mass. A more massive host halo is more likely to host a more luminous central and larger number of satellites with a given luminosity. Galaxy spatial and velocity distribution can also be taken care of in HOD. The position of a central galaxy is at the potential minimum, and positions of satellites can be set according to NFW (Navarro-Frenk-White; Navarro et al. (1997)) profile or other more realistic satellite profiles measured from simulations. Constraining halo occupation can provide insights into galaxy formation physics, and constrain cosmological parameters. Based on the halo catalog from simulations, HOD populates galaxies in halos with some free parameters, and then constrains those parameters through comparing populated galaxy clustering and observed galaxy clustering. Galaxy clustering is a result of the combination of cosmology and HOD. However, previous research shows that different cosmology will predict distinguishable galaxy 7 statistics, and changes in HOD cannot compensate for that (Zheng et al. (2005)). Conditional luminosity function (CLF; Yang et al. (2003)) shares the same spirit with HOD by modeling galaxy luminosity function at fixed halo mass. At a fixed halo mass, HOD can be obtained by integrating CLF over luminosity. However, traditional HOD and CLF relate galaxies to halo only by halo mass without considering halo assembly history. There are models developed to include other halo properties, to take galaxy assembly bias and halo assembly bias into consideration (Hearin et al. (2016); Xu et al. (2018)). 1.3.2 Physical galaxy formation models Semi-analytical models (De Lucia & Blaizot (2007); Guo et al. (2011)) implement physical or empirical galaxy formation models on the subhalo merger tree of N-body simulations. Stars form by infalling and cooling of gas to a central region of their subhalos and evolve to form different components of the galaxy, for example, the stellar bulge and the disk which co-evolve with other gas components. Angular momentum evolution of stellar and gas components can also be modeled. Supernova and AGN feedbacks which eject energy and suppress cold gas from infalling may also be included. Semi-analytical models set up a number of free parameters and constrain them with observables like galaxy luminosity function, galaxy clustering, galaxy color, and others. Semianalytical models are effective ways to model galaxy formation inside dark matter halos, but some parts of them are based on empirical prescriptions, and the physical details remain unclear. Another kind of physical model is direct hydrodynamical simulations(Vogelsberger et al. (2014b); McAlpine et al. (2016)). Similar to N-body simulation, hydrodynamical simulation evolves particles or meshes according to gravity and hydrodynamical processes. Hydrodynamical simulations could have many kinds of particles, for example, dark matter, star, gas, black hole, to simulate several physical processes of different components in galaxies and halos. Hydrodynamical simulations are capable of reproducing observables like galaxy luminosity function, star formation rate, and Tully-Fisher relation. Different simulations implement different aspects of physical processes, so they show different levels of consistency with observations. Hydrodynamical simulations follow physical processes in more details than semi-analytical models. However, galaxy formation physics could involve processes on a large range of scales. So it is computationally expensive to have a large simulation box. The box size of the state-of-the-art hydrodynamical simulations is approximately 8 100Mpc on one side. This limits the ability of investigating massive objects, but these simulations can still be used for studies on appropriate mass scales. In Chapter 3, I will perform a study of galaxy-halo relations and galaxy assembly bias with a hydrodynamical simulation. In Chapter 4, I will introduce a study of galaxy conditional color-magnitude distribution with both a semi-analytical model and a hydrodynamical simulation. CHAPTER 2 DEPENDENCE OF HALO BIAS AND KINEMATICS ON ASSEMBLY VARIABLES This chapter originally appeared in a work that has been accepted by the MNRAS (Xu & Zheng (2018a)) published by Oxford University Press on behalf of the Royal Astronomical Society. Reprinted with permission, sections and figures have been renumbered and reformatted to match the dissertation. The coauthor of the paper is Zheng Zheng. 2.1 Abstract Using dark matter haloes identified in a large N-body simulation, we study halo assembly bias, with halo formation time, peak maximum circular velocity, concentration, and spin as the assembly variables. Instead of grouping haloes at fixed mass into different percentiles of each assembly variable, we present the joint dependence of halo bias on the values of halo mass and each assembly variable. In the plane of halo mass and one assembly variable, the joint dependence can be largely described as halo bias increasing outward from a global minimum. We find it unlikely to have a combination of halo variables to absorb all assembly bias effects. We then present the joint dependence of halo bias on two assembly variables at fixed halo mass. The gradient of halo bias does not necessarily follow the correlation direction of the two assembly variables and it varies with halo mass. Therefore in general for two correlated assembly variables, one cannot be used as a proxy for the other in predicting halo assembly bias trend. Finally, halo assembly is found to affect the kinematics of haloes. Low-mass haloes formed earlier can have much higher pairwise velocity dispersion than those of massive haloes. In general, halo assembly leads to a correlation between halo bias and halo pairwise velocity distribution, with more strongly clustered haloes having higher pairwise velocity and velocity dispersion. However, the correlation is not tight, and the kinematics of haloes at fixed halo bias still depends on halo mass and assembly variables. 10 2.2 Introduction Large-volume galaxy redshift surveys over a range of redshifts, including the Sloan Digital Sky Survey (SDSS; York et al. (2000)), the Two-degree Field Galaxy Redshift Survey (2dFGRS; Collessy (1999)), the SDSS-III (Eisenstein et al. (2011)) and SDSS-IV (Dawson et al. (2016)), and the WiggleZ Dark Energy Survey (Blake et al. (2011)), have transformed the study of large-scale structure, producing detailed distribution of galaxies in the universe as a function of their properties and resulting in galaxy clustering measurements with high precision. Galaxy clustering data from such surveys play an important role in understanding galaxy formation and evolution and in learning about cosmology, in particular in constraining dark energy and testing gravitational theories. The formation of galaxies involves poorly understood baryonic processes, causing a hurdle to interpret galaxy clustering data. In contrast, the formation and evolution of dark matter haloes are dominated by gravitational interactions and their properties are well understood with analytic models and N-body simulations (Press & Schechter (1974); Mo & White (1996); Tinker et al. (2008)). Over the past two decades, an informative way to interpret galaxy clustering has been developed and made wide applications, which is to link galaxies to the underlying dark matter halo population. The two commonly adopted descriptions of the galaxy-halo connection are the halo occupation distribution (HOD) and conditional luminosity function (CLF) frameworks (e.g. Berlind & Weinberg (2002); Yang et al. (2003); Zheng et al. (2005)), which have been successfully applied to galaxy clustering data to infer the relation between galaxy properties and halo mass (Zehavi et al. (2005, 2011)). The HOD/CLF framework chooses the halo mass as the halo variable, and the implicit assumption is that the statistical properties of galaxies inside haloes only depend on halo mass, not on halo environment or growth history. This assumption of an environment-independent HOD/CLF is partly motivated by the excursion set theory (Bond et al. (1991)). However, N-body simulations show that the clustering of haloes depends on not only halo mass, but also halo assembly history (Gao et al. (2005); Zhu et al. (2006); Jing et al. (2007)), halo structure (e.g. Wechsler et al. (2006); Gao & White (2007); Faltenbacher & White (2009); Paranjape & Padmanabhan (2017)), and halo environment (e.g. Harker et al. (2006); Salcedo et al. (2018)), which is termed as halo assembly bias. Halo properties that correlate with halo environment or assembly history are broadly referred to here as halo assembly variables, such as halo formation time, halo concentration, maximum circular velocity of halo, and halo spin. Studying the origin of halo assembly bias remains active ongoing effort (e.g. Borzyszkowski et al. (2017)). For massive haloes (with mass higher than the 11 nonlinear mass Mnl for collapse), the assembly bias appears to be a generic feature in the Extended Press-Schechter theory, related to the difference in the curvature of peaks of the same height in a Gaussian density field (e.g. Dalal et al. (2008)). For low mass haloes, the assembly bias is proposed to be originated from the strong influence of the environment, especially the tidal field, on the evolution of haloes (e.g. Hahn et al. (2007, 2009); Wang et al. (2011); Shi et al. (2015)). Whether halo assembly bias is inherited by galaxies is still under investigation, in both simulations (e.g. Chaves-Montero et al. (2016); Busch & White (2017); Garaldi et al. (2017)) and observations (e.g. Lin et al. (2016); Miyatake et al. (2016); Guo et al. (2017); Zu et al. (2017)). If galaxy properties are closely tied to halo growth history, the inherited assembly bias from the host haloes would require a modification of the current halo model of galaxy clustering to include such an effect. Otherwise, we would infer incorrect galaxy-halo connections and introduce possible systematics in cosmological constraints (e.g. Zentner et al. (2014); Hearin (2015); Zentner et al. (2019); but see McEwen & Weinberg (2016)). Along the path, a better characterisation of halo assembly bias is necessary, which motivates the work in this paper. In most previous studies, halo assembly bias is inferred by grouping each halo assembly variable at fixed halo mass into certain percentiles. While this can reveal the trend of assembly bias, it does not provide a clear view on the multivariate dependence of the halo bias, which needs the measurement of halo bias at certain values, not certain percentiles, of a given assembly variable. If halo assembly effect is to be incorporated into the halo modelling of galaxy clustering, a natural route is to describe halo bias in terms of multiple halo properties. We will present the multivariate dependence of halo bias and study whether there is a combination of halo variables to minimise assembly bias. In addition, we also investigate the scale dependence of halo assembly bias and the assembly effect on halo kinematics. The structure of the paper is as follows. In Section 2.3, we describe the simulation data used in this study and the measurement of halo bias. In Section 2.4, we present the joint dependence of halo bias on halo mass and each halo assembly variable, study an effective halo variable to minimise halo assembly bias, show the scale dependence of halo assembly bias, and describe the dependence of halo bias on two halo assembly variables. In Section 2.5, we analyse the assembly effect on halo kinematics. Finally, in Section 2.6, we summarise our results and discuss their implications. 12 2.3 Simulation and halo bias calculation For the work presented in this paper, we make use of the MDR1 MultiDark simulation (Prada et al. (2012); Riebe et al. (2013)) 1. The MDR1 N-body simulation adopts a spatially flat ΛCDM cosmology, with Ωm = 0.27, Ωb = 0.0469, h = 0.70, ns = 0.95, and σ8 = 0.82. The simulation has 20483 particles in a box of 1h−1 Gpc on one side, with particle mass 8.721 × 109 h−1 M . Dark matter haloes are identified with the Rockstar algorithm (Behroozi et al. (2013)), which finds haloes through adaptive hierarchical refinement of friends-of-friends groups using the six-dimensional phase-space coordinates and the temporal information. While the MDR1 simulation is adopted for presenting most of the results, we also make use of the Bolshoi simulation for some investigations that need high completeness in certain regions of the parameter space. The Bolshoi simulation (Klypin et al. (2011)) 2 with a 250h−1 Mpc box adopts the same cosmology as MDR1, but with 64 times higher mass resolution and 7 times higher force resolution. We use the z = 0 Rockstar halo catalogue. In addition to halo mass Mh , for which we adopt M200b (i.e. the mean density of a halo being 200 times the background density of the universe), we consider four other halo properties: Vpeak , peak maximum circular velocity of the halo over its accretion history; λ, halo spin parameter that characterises its angular momentum; aM/2 , cosmic scale factor when the halo obtains half of its current (z = 0) total mass, which characterises halo formation time; c, halo concentration parameter, which is the ratio of the halo virial radius Rvir to the scale radius rs . All the four quantities are related to the halo assembly history and therefore are dubbed as halo assembly variables. We measure halo bias for haloes above ∼ 4 × 1011 h−1 M (about 50 particles). The large simulation box also serves well our purpose of studying the joint dependence of halo bias on multiple halo properties, as it enables halo bias in fine bins of halo properties to be measured. To further reduce the uncertainty in the halo bias measurements, we derive halo bias from the two-point halo-mass cross-correlation function ξhm and the two-point mass auto-correlation function ξmm (e.g. Jing (1999); Jing et al. (2007)), i.e. b(r) = ξhm (r) . ξmm (r) 1 https://www.cosmosim.org/cms/simulations/mdr1/ 2 https://www.cosmosim.org/cms/simulations/bolshoi/ (2.1) 13 In practice, for a given sample of haloes and a pair separation bin r ± dr/2, we count the number of halo-mass particle pairs, Nhm (r ± dr/2), and that of mass particle-particle pairs, Nmm (r ± dr/2). With the periodic boundary condition, the corresponding pair counts from randomly distributed haloes and mass particles can be simply calculated as Nhm,ran (r ± dr/2) = nh nmV d 3r and Nmm,ran (r ± dr/2) = 2 V d 3 r/2, where n and n are the halo and mass particle number densities in the simulation box nm h m of volume V. The two-point correlation functions are then computed as ξhm (r) = Nhm (r ± dr/2) −1 Nhm,ran (r ± dr/2) (2.2) ξmm (r) = Nmm (r ± dr/2) − 1. Nmm,ran (r ± dr/2) (2.3) and 2.4 Joint dependence of halo bias on mass and assembly variables In this section, we first present the results of the joint dependence of halo bias on both mass and one of the assembly variables. Then we investigate whether an effective halo mass can be constructed to absorb the assembly effect in halo clustering. We also discuss how the halo assembly affects the scale dependence of halo bias. Finally, we present the dependence of halo bias on two halo assembly variables at fixed halo mass. 2.4.1 Dependence of halo bias on mass and one assembly variable To separate the effects of mass dependence and assembly dependence, we first choose the halo mass bin to have a narrow width, ∆ log Mh = 0.2. For each halo mass bin, we further divide haloes into bins of each assembly variable, and typically, there are about 20 assembly variable bins. We measure halo bias for haloes in each bin of mass and assembly variable. We take the average above 10h−1 Mpc as the large-scale halo bias factor, where the ratio between ξhm and ξmm has becomes flat. To maintain reasonable signal-to-noise ratios, we only consider bins with more than 100 haloes. The joint dependences of halo bias on mass and each assembly variable are shown in Fig. 2.1.3 3 While the plot does not show the level of uncertainties in the contours, it is clear that the main trend is not affected by noise. The contribution of noise distorts the otherwise smooth contours, showing up as jags in the contour curves or alterations of contour levels (such as the small fluctuation seen in the lower part of the bottom-right panel of Fig. 2.1). Rather than providing a corresponding plot of the noise levels, we suggest to use the fluctuations in the map as an estimate of the magnitude of noise. This rule of thumb also applies to curve plots, and in general, the deviation from a smooth 14 3.0 −0.1 3.0 2.5 −0.2 2.5 aM/ b 1.5 2.4 1.0 2.2 1.5 −0.4 1.0 −0.5 0.5 0.5 −0.6 12.0 Mh / h M ⊙ )] 12.5 log[ 13.0 13.5 ( −1 14.0 14.5 1.8 0.0 12.0 Mh / h M ⊙ )] 12.5 log[ 13.0 13.5 ( −1 14.0 14.5 3.5 1.6 3.5 3.0 1.4 3.0 −1.0 2.5 2.5 1.2 b 1.5 0.8 −1.5 λ 2.0 log c 2.0 1.0 0.0 b 2.0 log 2.0 −0.3 log peak 3.5 2.0 2.6 log[ 0.0 2 2.8 3.5 b 1 V /(km s− )] 3.0 1.5 −2.0 0.6 1.0 0.4 0.5 0.2 12.0 Mh /(h − M ⊙ )] 12.5 log[ 13.0 13.5 1 14.0 14.5 0.0 1.0 0.5 −2.5 12.0 Mh /(h − M ⊙ )] 12.5 log[ 13.0 13.5 1 14.0 14.5 0.0 Figure 2.1: Joint dependence of halo bias on halo mass and each assembly variable. The assembly variables in the four panels are peak maximum circular velocity Vpeak , halo formation scale factor aM/2 , halo concentration c, and halo spin λ, respectively. The four black curves in each panel mark the central 50 and 80 percent of the distribution of the corresponding assembly variable as a function of halo mass. 15 The black curves in each panel delineate the central 50 and 80 percent of the distribution of the corresponding assembly variable as a function of halo mass. The lowest mass bin in the investigation is ∼ 4 × 1011 h−1 M , corresponding to about 50 particles. Paranjape & Padmanabhan (2017) show that for haloes with less than 400 particles, the distribution of halo concentration is not converged (also see Trenti et al. (2010)). The spin distribution is also affected for haloes with low number of particles (e.g. Trenti et al. (2010); Benson (2017)). To see how the resolution effect changes the results, we perform a test (see Appendix A) of replacing haloes of log[Mh /(h−1 M )] < 12.8 with those in the Bolshioi simulation, which ensures that each of the lowest mass haloes in our analysis contains more than 3000 particles. We find that the overall trend shown in Fig. 2.1 still holds. Given the higher signal-to-noise ratio and larger range in assembly variables with the MDR1 simulation, we present the halo bias results using the MDR1 simulation except for the analysis in Section 2.4.2, where extending to much lower halo mass is needed. For halo pairwise velocity statistics, a caveat is made in Section 2.5. The top-left panel shows the joint dependence of halo bias on halo mass Mh and peak maximum circular velocity Vpeak . At low halo masses, halo bias increases as Vpeak increases, while at high halo masses, the trend is reversed. The transition mass is around log[Mh /(h−1 M )] = 13.0, which is about three times the nonlinear mass for collapse (log[Mnl /(h−1 M )] = 12.5 for the adopted cosmology). Haloes with low mass and high Vpeak have a bias factor comparable to that of high mass haloes, indicating that they could originate from highly stripped high mass haloes (e.g. Wang et al. (2011)). In the top-right panel, it can be seen that at fixed mass, haloes formed earlier (i.e. with lower aM/2 ) are more strongly clustered. The trend becomes weak for haloes above ∼ 1014 h−1 M , manifested by the almost vertical contours. Extremely old haloes (the bottom boundary in the panel) at all masses are all as strongly clustered as massive haloes. At fixed mass, old haloes tend to be more concentrated, i.e. with higher concentration parameters c. The dependence of halo bias on concentration (bottom-left panel) follows a similar trend as with formation time (see § 2.4.4 for the condition for this to be valid). The dependence becomes weak around 3Mnl , which is consistent with previous work of direct measurements by binning halo concentration (e.g. Wechsler et al. (2006); Jing et al. (2007)) and high precision inference with a curve gives us the level of noise. However, we do add error bars for a few curves in the curve plots of this paper, and in most cases, the noise level is low and the trend is not affected by the noise. 16 lognormal model and Separate Universe technique (Paranjape & Padmanabhan (2017)). We note that the latter shows concentration-dependent assembly bias as a function of the actual value of concentration, instead of percentiles of concentration. However, at higher halo masses, halo bias decreases with increasing halo concentration, clearly different from that in the aM/2 case. This difference implies that halo concentration c can depend on quantities other than formation time (aM/2 ), like the mass accretion rate (e.g. Wechsler et al. (2002)). In the bottom-right panel, the joint dependence of halo bias on halo mass and spin is shown. The distribution of halo spin (indicated by the black curves) is only weakly dependent on halo mass. Overall, the contours show an ordered pattern (except for the extremely low spin tail of low mass haloes), with haloes of higher spin more strongly clustered. Unlike the cases in the other three panels, the trend persists over all halo masses, which is consistent with the result using percentiles of spin distribution (Gao & White (2007)). Compared to studies of halo bias for haloes selected from percentiles of assembly variables, our results show clearly the joint dependence of halo bias as a function of halo mass and the value of each assembly variable. In the Mh –A (A being one of the assembly variables) plane, the common description of our results can be that halo bias increases outward from a point of global minimum. Such a description is not obvious if the investigation is limited to the main distribution range of the assembly variable. The previous results can all be understood by considering slices in the Mh –A plane and by noticing the shape and orientation of the contours. With the above description motivated by our calculation, it would be interesting to ask whether we can define a new halo variable to account for both the halo mass and assembly history dependences of halo bias. We do such an exercise below. 2.4.2 Effective mass to absorb the assembly bias effect? To see whether we can construct a halo variable from combinations of halo mass and assembly variables to minimise the assembly history effect on halo bias, we investigate an example with halo mass Mh and spin λ. As shown in Fig. 2.1, the distribution of halo spin is almost independent of halo mass and the joint dependence of halo bias on Mh and λ has a pattern easy to describe. Within the most part of the parameter space extended by log Mh and log λ, a given value of halo bias factor approximately corresponds to a straight line (see the dashed lines in Fig. 2.2 as 17 3.5 3.0 −1.0 2.5 −1.5 b log λ 2.0 1.5 13 14 −2.0 1 11 2 10 log −2.5 10.5 11.0 Mh /(h − M ⊙ )] 11.5 log[ 5 12.0 12.5 13.0 1 13.5 M 14.0 1.0 0.5 eff 0.0 14.5 original Mh effective Mh 4 b 3 2 1 0 11.0 11.5 12.0 Mh / h M ⊙ )] log[ 12.5 13.0 13.5 ( −1 14.0 14.5 15.0 Figure 2.2: Effective halo mass to absorb assembly bias effect with halo spin. The top panel is the same as the bottom-right panel of Fig. 2.1, and dashed lines are added to illustrate how the effective mass is defined (see text). In the bottom panel, the black points show the halo bias as a function of halo mass, with the scatter at fixed halo mass coming from haloes with different spin parameters (assembly bias effect with spin). The red points show the dependence of halo bias on effective halo mass, and the scatter caused by the assembly effect with spin is substantially reduced. 18 examples). We find that the set of straight lines can be well described as converged to a single point (log Mh,0 , log λ0 ), located towards the bottom-right corner of the left panel of Fig. 2.2. The slope (log Mh − log Mh,0 )/(log λ − log λ0 ) of each line can be defined as a new halo variable, which is a combination of halo mass Mh and spin λ. Halo bias is monotonically dependent on this variable. Equivalently, for each line, we can evaluate the halo mass at a fixed spin value and use this halo mass as the new variable. We choose the fixed spin value to be log λeff = −1.4, roughly the value averaged over all haloes, and define the corresponding halo mass as the effective halo mass Meff through log Meff = log Mh,0 + log λeff − log λ0 (log Mh − log Mh,0 ). log λ − log λ0 (2.4) To cover a large range in effective halo mass, we perform the analysis using a combination of haloes in both the MDR1 and Bolshoi simulations. As seen from Fig. 2.2, for an effec- tive mass of log[Meff /(h−1 M )] ∼ 11–12.5, we have contributions from haloes of mass below log[Mh /(h−1 M )] ∼ 11.6 with high spin. Furthermore, there is numerical effect for haloes with less than a few hundred particles (e.g. Paranjape & Padmanabhan (2017); Trenti et al. (2010); Benson (2017)) in the MDR1 simulation. To include low mass haloes and minimise the potential numerical effect, we calculate effective halo masses for both MDR1 and Bolshoi haloes, then present the results using MDR1 haloes for log[Meff /(h−1 M )] ≥ 12.7 and Bolshoi haloes for log[Meff /(h−1 M )] < 12.7. To verify that the effect of assembly bias in spin is absorbed into the effective halo mass, we measure halo bias in fine bins of log Meff and log λ. In the right panel of Fig. 2.2, the black circles show the dependence of halo bias on the original halo mass Mh , and the spread at fixed mass reflects the assembly bias effect in halo spin. The red dots represent the dependence of halo bias on the effective halo mass Meff . At fixed Meff , the scatter is much smaller than that seen in the b–Mh dependence, demonstrating that the assembly bias effect in spin has been well absorbed into the effective mass. Fig. 2.3 is similar to Fig. 2.1, but for the joint dependence of halo bias as a function of effective mass Meff and each of the assembly variables. The fact that the assembly bias effect in spin has been absorbed into Meff shows up as largely vertical contours in the log Meff –log λ plane (bottom-right panel). However, in the other three panels of Fig. 2.3, where the joint dependence of halo bias on Meff and the other assembly variables are shown, assembly bias still exists at fixed Meff . The assembly 19 3.5 2.8 2 1.0 1.8 b b 1.5 2.0 2.0 −0.3 1.5 −0.4 1.0 −0.5 0.5 1.6 0.5 −0.6 11.5 Mh / h M ⊙ )] 12.0 log[ 12.5 13.0 ( −1 13.5 14.0 14.5 1.8 0.0 11.5 Mh / h M ⊙ )] 12.0 log[ 12.5 13.0 ( −1 13.5 14.0 14.5 3.5 1.6 3.5 3.0 1.4 3.0 −1.0 2.5 2.5 1.2 b 1.5 0.8 −1.5 λ 2.0 log c 2.0 1.0 0.0 b peak 2.5 aM/ 2.0 2.2 log 3.0 −0.2 2.5 2.4 log[ −0.1 3.0 2.6 1.4 3.5 log 1 V /(km s− )] 3.0 1.5 −2.0 0.6 1.0 0.4 0.5 0.2 11.5 Mh / h M ⊙ )] 12.0 log[ 12.5 13.0 ( −1 13.5 14.0 14.5 0.0 1.0 0.5 −2.5 11.5 Mh / h M ⊙ )] 12.0 log[ 12.5 13.0 ( −1 13.5 14.0 14.5 0.0 Figure 2.3: Same as Fig. 2.1, but with original halo mass replaced by effective halo mass. The effective halo mass absorbs the halo assembly bias with halo spin, manifested by the vertical contours in the bottom-right panel. However, assembly bias effect still exists with other assembly variables, as seen in the other three panels. bias effect becomes weaker at the high Meff end. This is not surprising, given that the effect is weak at high Mh in Fig. 2.1 and the λ dependence appears to be the strongest one. In Fig. 2.4, we show the assembly bias effect as a function of the original halo mass (left panel) and the effective halo mass (right panel), by grouping haloes into 10–25 percent and 75–90 percent percentiles of different assembly variables. Error bars are calculated with the jackknife method. For clarity, we only plot error bars for b(M) (black solid) and aM/2 (dark blue). Under the effective halo mass, the assembly bias in halo spin almost disappears, and the curves of bias from the two percentile subsamples fall on top of each other (on the line of unity, corresponding to the average bias as a function of effective mass), while under the original halo mass, the assembly effect in halo spin is at the level of ∼15 percent. Above log[Meff /(h−1 M )] ∼ 13.5, the assembly effects in Vpeak and c are also reduced, but not in aM/2 . However, towards lower Meff , the assembly bias effect in other assembly variables becomes enhanced, compared to the case with the original halo mass. Therefore, the effective mass does not work all the way to reduce assembly bias effect in other variables and 20 1.4 b(A; M)/b(M) 1.3 1.2 -__ V aM/ 10%~25% peak 75%~90% 2 c λ 1.1 1.0 0.9 0.8 original Mh 0.7 0.6 12.0 12.5 Mh / h M ⊙ )] log[ 13.0 13.5 ( −1 14.0 14.5 1.4 b(A; M)/b(M) 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 11.0 effective 11.5 12.0 Mh Mh / h M ⊙ )] log[ 12.5 13.0 13.5 ( −1 14.0 14.5 15.0 Figure 2.4: Comparison between assembly bias effects under original halo mass (top) and effective halo mass (bottom). At fixed original or effective halo mass, haloes are grouped according to the percentiles of each assembly variable, and their bias values normalised by the mass-only dependent bias are shown as curves in each panel (dashed for the 10–25 percent percentile and solid for the 75–90 percent percentile). over the full mass range. With the joint dependence of halo bias on Meff and assembly variables, in principle, we can perform a similar variable transformation to largely reduce the assembly effect in another chosen variable, and the analysis can be repeated until the assembly effect in all the assembly variables of interest is recursively absorbed. It may be possible to find a principal direction or locus in the multidimensional space spanned by halo mass and assembly variables such that halo bias has the 21 strongest dependence on the corresponding variable combination. However, there is a limit on how far we would like to go following such a path. If the ultimate goal is to replace halo mass with the putative variable combination in halo modelling of galaxy clustering, it is desirable that the variable combination has a tight correlation with galaxies properties. It has been established that halo mass plays the dominant role in shaping galaxy formation and evolution and in determining the main properties of galaxies (e.g. White & Rees (1978); Birnboim & Dekel (2003); Kereš et al. (2005)). Although the combination variable Meff absorbs the assembly bias effect in halo spin, a fixed Meff spans a large range in halo mass. For example, haloes of effective mass Meff = 1012 h−1 M , which corresponds to the middle dashed line in the left panel of Fig. 2.2, can come from haloes of original mass Mh ∼ 1011 h−1 M with high spin or those of Mh ∼ 1013 h−1 M with low spin. The large range in original halo mass suggests that Meff would not be a good variable to choose for a tight correlation with galaxy properties. Together with the results in Fig. 2.3 and Fig. 2.4, it implies that it is unlikely to find the right combination of halo variables to completely absorb the assembly bias effect in every assembly variable and at the same time to keep a close connection to galaxy properties. 2.4.3 Scale dependence of halo bias with assembly variable While the assembly history of haloes affects large-scale clustering, it can also have an effect on small-scale clustering. In Fig. 2.5, the scale dependences of halo bias on assembly variables are shown at two halo masses, log[Mh /(h−1 M )]=12 and 14, respectively. The top eight panels show the original bias factors b(A; r) as a function of r, calculated by grouping haloes into various percentiles of each assembly variable A. To emphasise the shape difference on small scales, in the bottom eight panels, we also show halo bias factors normalised by the corresponding large-scale value (average values above 10h−1 Mpc), i.e. b(A; r)/hb(A; r > 10h−1 Mpc)i. For clarity, in each panel, error bars are only plotted for the curve corresponding to all haloes to indicate the uncertainties. In the non-normalised case, for Mh = 1012 h−1 M haloes, all halo bias factors (solid curves in the top panels) start to drop at ∼ 3h−1 Mpc towards small scales, which means that the two-point halo-mass cross-correlation function is shallower than the mass auto-correlation function on such scales. Note that the drop occurs on scales larger than the halo virial radius (Rvir ∼ 0.21h−1 Mpc). The panels show that halo assembly history influences how steep the bias factor drops towards small scale. Except for the halo spin case, all other cases have a similar trend that can be related through the correlation between halo assembly variables. Haloes that form earlier (smaller aM/2 ), thus on 22 b(A; r) 1.4 1.2 logM=12 logM=12 logM=12 1.0 logM=12 0.8 0.6 0.4 0.2 aM/2 dependence λ dependence Vpeak dependence c dependence 0.0 from ξhh /ξmm lower 25% upper-mid 25% 3.5 lower-mid 25% upper 25% All 3.0 2.5 2.0 1.5 logM=14 logM=14 logM=14 1.0 logM=14 Vpeak dependence aM/2 dependence λ dependence c dependence 0.5 0.1 1.0 10.0 0.1 1.0 10.0 0.1 1.0 10.0 0.1 1.0 10.0 r(h −1 Mpc) 1.2 1.0 logM=12 logM=12 logM=12 logM=12 b(A; r)/ b(A; r> 10h −1 Mpc) 0.8 0.6 0.4 0.2 0.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0.1 Vpeak dependence aM/2 dependence from ξhh /ξmm All lower 25% lower-mid 25% logM=14 Vpeak dependence 1.0 10.0 0.1 logM=14 aM/2 dependence 1.0 10.0 0.1 r(h Mpc) −1 λ dependence c dependence upper-mid 25% upper 25% logM=14 c dependence 1.0 10.0 logM=14 λ dependence 0.1 1.0 10.0 Figure 2.5: Scale dependence of assembly bias. The upper subplot shows the scale-dependent halo bias factor for different assembly variables. In each panel, solid curves are calculated through halomass cross-correlation functions and dotted curves through halo-halo auto-correlation functions. Each set of five curves are for haloes in different percentiles of each assembly variable: lower 25 percent (blue), lower-middle 25 percent (cyan), upper-middle 25 percent (green), upper 25 percent (red), and all haloes in the mass bin (black). Top and bottom panels are for haloes of log[Mh /(h−1 M )]=12 and 14, respectively. Bottom subplot is similar, but each halo bias curve is normalised by the corresponding large-scale value (average on scales larger than 10h−1 Mpc). 23 average with higher halo concentration c and higher Vpeak , have a steeper drop in halo bias towards smaller scales. The dotted curves in each panel are the measurements of halo bias using halo-halo auto-correlation functions, and the trend is similar to that derived from halo-mass cross-correlation functions. For the case with halo spin, the overall shape difference is smaller than seen in the other three cases. As halo spin shows an anti-correlation with halo concentration (e.g. Maccio et al. (2007); also Fig. 2.6 and Fig. 2.7), one would expect the trend to be that the scale-dependent bias for haloes with lower spin drops more steeply towards smaller scales. However, the result is opposite (see § 2.4.4 for an interpretation). The bottom panels of the non-normalised case in Fig. 2.5 show the results with Mh = 1014 h−1 M haloes. As the virial radius Rvir is ∼ 0.98h−1 Mpc, the rise on small scales reflects the mass distribution inside haloes. For all cases, the dotted curves for halo biases from halo-halo autocorrelation functions have a cutoff around r ∼ 2Rvir , a manifestation of the halo exclusion effect. Above ∼ 2Rvir , the shape of scale-dependent halo bias factors do not show a strong dependence on assembly variables (shown more clearly in the corresponding panels for the normalised case). Overall, the shape of the scale dependence of halo bias at small scales depends on the value of the assembly variable in consideration, as shown in the normalised case. The dependence is not strong, but it is discernible for low mass haloes. For high mass haloes, the dependence becomes weak (except for the spin case). In such a situation, the scale-dependent assembly bias at a given halo mass Mh can be well approximated as b(Mh , A; r) = hb(Mh , A)i f (Mh ; r), where hb(Mh , A)i is the large-scale bias with the value of assembly variable being A and f (Mh ; r) characterise the shape (being unity on large scales by construction). Sunayama et al. (2016) study the scale dependence of halo assembly bias for the case of Vmax , the maximum circular velocity at z = 0. They find that for low-mass haloes (Mh . Mnl ), the ratio of biases between high Vmax and low Vmax haloes exhibits a pronounced scale dependence at 0.5–5h−1 Mpc, and the scale dependence becomes weak towards higher halo mass. Our results, if put in a similar format, are in broad agreement with their findings. They also show that the scale dependence can be partially attributed to haloes previously residing in higher mass haloes but ejected to become host haloes at the epoch of interest, which is consistent with the proposed origin of assembly bias for low mass haloes (e.g. Wang et al. (2009); Wetzel et al. (2014)) If galaxy properties track the halo assembly history, the scale-dependent bias from the halo 24 2.50 2.50 1.4 2.25 1.4 2.00 1.2 0.8 1.00 0.6 0.25 −0.5 aM/ −0.4 −0.3 log −0.2 −0.1 0.0 V /(km s− )] 2.0 2.1 2.2 2.3 log[Vpeak 2.4 2.5 2.6 /(km s −1 )] 2.7 −0.5 2.50 2.25 2.00 2.4 b 1.25 2.3 1.00 0.75 2.2 1.75 1.50 −1.5 1.25 1.00 0.75 −2.0 0.50 2.1 2.00 −1.0 1.75 1.50 0.50 0.25 −0.5 aM/ −0.4 −0.3 log −0.2 −0.1 0.0 0.25 0.00 −2.5 0.2 0.4 0.6 2 −0.5 2.50 0.8 1.0 logc 1.2 1.4 1.6 −0.5 2.25 2.00 1.00 0.75 −2.0 0.50 1.75 λ 1.50 −1.5 1.25 b b 1.25 log logλ −1.5 2.00 −1.0 1.75 1.50 1.00 0.75 −2.0 0.50 0.25 −2.5 2.0 2.1 2.2 2.3 log[Vpeak 2.4 2.5 2.6 /(km s −1 )] 2.7 0.00 0.00 2.50 2.25 −1.0 0.00 b 1 0.25 2.25 2.5 peak 0.50 0.2 0.00 2.50 2.6 log[ 0.75 2 2.7 2.0 −0.6 1.00 0.4 logλ −0.6 1.25 0.6 0.50 0.2 1.50 0.8 0.75 0.4 1.75 1.0 logc 1.25 b log c 1.50 2.00 1.2 1.75 1.0 2.25 b 1.6 0.25 −2.5 −0.6 −0.5 −0.4 aM/ −0.3 log −0.2 −0.1 0.0 0.00 2 Figure 2.6: Halo bias as a function of two assembly variables for haloes in the mass range of log[Mh /(h−1 M )] = 12.5 ± 0.5. In each panel, the contours mark the 68 percent and 95 percent distribution of the two assembly variables and the color scale shows the value of halo bias. 25 2.50 1.2 1.75 1.25 1.00 0.6 2.00 1.75 1.0 1.50 logc 0.8 b log c 1.50 2.25 1.2 2.00 1.0 2.50 1.4 2.25 0.8 1.25 b 1.4 1.00 0.6 0.75 0.4 0.25 0.2 −0.5 −0.4 aM/ −0.3 log −0.2 −0.1 0.0 V /(km s− )] 2.50 2.6 log[Vpeak 2.7 2.8 2.9 /(km s −1 )] 3.0 2.50 2.25 1.25 b 2.7 1.00 2.6 0.75 1.75 1.50 −1.5 1.25 1.00 0.75 −2.0 0.50 2.5 2.00 −1.0 1.75 1.50 0.50 0.25 2.4 −0.6 −0.5 −0.4 aM/ −0.3 log −0.2 −0.1 0.0 0.25 0.00 −2.5 0.2 0.4 0.6 2 −0.5 2.50 0.8 1.0 logc 1.2 1.4 −0.5 2.25 2.00 1.00 0.75 −2.0 0.50 1.75 λ 1.50 −1.5 1.25 b b 1.25 log logλ −1.5 2.00 −1.0 1.75 1.50 1.00 0.75 −2.0 0.50 0.25 −2.5 2.4 2.5 2.6 log[Vpeak 2.7 2.8 2.9 /(km s −1 )] 3.0 0.00 0.00 2.50 2.25 −1.0 0.00 b 1 2.5 −0.5 2.00 2.8 peak 0.25 2.25 2.9 log[ 2.4 2 3.0 0.50 0.2 0.00 logλ −0.6 0.75 0.4 0.50 0.25 −2.5 −0.6 −0.5 −0.4 aM/ −0.3 log −0.2 −0.1 0.0 0.00 2 Figure 2.7: Same as Fig. 2.6, but for haloes in the mass range of log[Mh /(h−1 M )] = 13.5 ± 0.5. 26 assembly effect could be a way to reveal the assembly bias in galaxy clustering. However, the results in Fig. 2.5 show that the scale dependence from assembly effect is not strong, and the overall shape does not vary significantly with the magnitude of any assembly variables, at a level of ∼10 percent on scales around 3h−1 Mpc for low mass haloes (see the normalised curves in Fig. 2.5). Also on scales where the scale dependence has a relatively large amplitude, the one-halo contribution to the galaxy correlation function would tend to mask it. Therefore, in practice, it is probably difficult to discern the assembly effect from the scale dependence in galaxy clustering. 2.4.4 Dependence of halo bias on two assembly variables At fixed halo mass, Fig. 2.1 shows the dependence of halo bias on each assembly variable. While in general, assembly variables are correlated with each other, there are differences in the assembly bias trends with different assembly variables. To better understand the trends, we present the correlation between any two assembly variables and how halo bias depends on them at fixed halo mass. Each panel of Fig. 2.6 shows halo bias as a function of two assembly variables, measured with haloes of mass in the range log[Mh /(h−1 M )] = 12.5 ± 0.5 (i.e. about Mnl ). The contours enclose the 68 and 95 percent distribution of the two assembly variables, respectively, which also reveal the correlation between them. As can be seen in the top-left panel, halo concentration c and formation scale factor aM/2 are anti-correlated, i.e. older haloes are more concentrated. The gradient of the halo bias follows the correlation direction. In such a case, the trend of dependence of halo bias on one variable can be used to predict that on the other variable. For example, at the above mass, given that more concentrated haloes are more strongly clustered and generally older, we can infer that older haloes are more strongly clustered than younger haloes. This is indeed what is seen in Fig. 2.1. Besides c–aM/2 , the bias dependences on c–Vpeak (top-middle panel) and Vpeak –aM/2 (top-right panel) also follow a similar behaviour. However, as halo spin λ becomes involved, the above picture changes. For example, in the bottom-right panel, halo spin λ is positively correlated with aM/2 , meaning that haloes of lower spin are generally older. Unlike the previous cases, now the gradient of the halo bias is not along the correlation direction but about perpendicular to it, and the correlation can no longer be used to predict the bias dependence trend. That is, the fact that older haloes are more strongly clustered and have a lower spin does not imply that haloes of lower spin are more strongly clustered. In this 27 example, the contrary is true (haloes of lower spin are less clustered as seen in the panel and in Fig.2.1). The difference between the halo bias gradient and the correlation direction with two assembly variables can help resolve apparent puzzles in assembly bias trend with other variables, e.g. the dependence on halo age and the number of subhaloes (above some mass threshold). The occupation number of subhaloes is found to correlate with halo age in the sense that older haloes have fewer subhaloes (e.g. Gao et al. (2004)), resulting from dynamical evolution and destruction of subhaloes. As in the mass range considered here older haloes are more strongly clustered, the correlation would imply haloes of fewer subhaloes are more strongly clustered. However, the dependence on subhalo occupation number is found to be opposite to such a naive expectation (e.g. Croft et al. (2012)). Fig. 2.7 shows halo bias as a function of two assembly variables for massive haloes (log[Mh /(h−1 M )] = 13.5 ± 0.5, above Mnl ). Compared to the case with low mass haloes (Fig. 2.6), while the correlation between each two assembly variables does not change much, the direction of the gradient of halo bias can be substantially different. For example, the gradient direction in the c–aM/2 panel now becomes perpendicular to the correlation direction. As a consequence, the trend in assembly bias effect and its relation to the correlation also changes. For cases involving halo spin, the gradient direction has only a mild change, indicating that the origin of spin bias is different from others (e.g. Salcedo et al. (2018)). With the Separate Universe technique, Lazeyras et al. (2017) present the dependence of halo bias on two assembly variables in two halo mass bins, although the dependence is not compared with the correlation of each two variables. They show that the dependence changes with halo mass and the trend is weak if halo shape is used as one variable. Han et al. (2019) present analysis of the multidimensional dependence of bias on halo properties. Their results, if projected onto the space of two assembly variables at fixed halo mass, can be compared to ours. Overall, we see that the joint dependence of halo bias on two assembly variables does not necessarily follow the correlation between the variables. One should be cautious in inferring assembly bias trend in one variable based on its correlation to the other variable. A similar conclusion is reached by Mao et al. (2018) for the so-called “secondary bias” with cluster-size haloes. The pattern of the joint dependence varies with halo mass, which can be characterised by a rotation (e.g. in terms of the halo bias gradient) as halo mass increases. 28 2.5 Dependence of pairwise velocity and velocity dispersions on halo assembly The assembly bias from the spatial distribution or clustering of haloes is related to halo environment, which should also affect the velocity distribution of haloes. The velocity distribution of haloes is a major ingredient in modelling redshift-space clustering of galaxies. To yield insights on how halo assembly may affect redshift-space clustering, we investigate the halo assembly effect on halo velocity distribution. We present the results in terms of halo pairwise velocity, v12 = v2 − v1 , as a function of halo pair separation, r = |r12 | = |r2 − r1 |, where vi and ri (i=1, 2) are the velocities and positions of a pair of haloes. The pairwise radial and transverse velocities are calculated as v12,r = (v12 · n12 )n12 , and v12,t = v12 − v12,r , where n12 = r12 /|r12 | is the direction from one halo to the other of the pair. The 2 i − hv 2 2 2 corresponding velocity dispersions are σ12,r = hv12,r 12,r i and σ12,t = hv12,t i − hv12,t i , with the average over all halo pairs at a given separation r. Before moving on to present the results on pairwise velocity statistics, we point out the comparison between those from the Bolshoi and MDR1 simulations (see Appendix A). In general, the MDR1 results in similar dependence patterns of the pairwise velocity statistics on assembly variables, except for halo spin. The weak trend with spin seen in the MDR1 simulation does not show up in the Bolshoi simulation, which may be caused by noise in the Bolshoi simulation or unknown systematic effect in determining halo spin in low-resolution simulations. For consistency, we choose to present the MDR1 results with the above caveat for the spin results (more details in Appendix A). 2.5.1 Assembly effect on pairwise velocity of haloes In Fig. 2.8, we first show the halo mass-dependent pairwise velocity and velocity dispersions, as a function of halo pair separation. Consistent with previous studies (e.g. Zheng et al. (2002)), the pairwise infall velocity v12,r (left panel) increases as halo mass increases, and the increase is faster at smaller pair separation, reflecting the stronger gravitational influence from higher mass haloes. For high mass haloes (with mass above a few times Mnl ), the pairwise radial velocity continuously decreases toward large separation. For low mass haloes, a bump in the pairwise radial velocity shows up. The bump is around r ∼ 5h−1 Mpc for haloes of the lowest mass in the study and shifts to larger scales for haloes of higher mass. At large separation, the pairwise infall velocity tends to converge to an amplitude independent of halo mass, which reflects the fact that all haloes feel the 29 log M 14.53 14.33 14.13 13.93 13.73 | 13.53 13.33 13.13 12.93 12.73 12.53 12.33 12.13 11.93 11.73 800 700 600 500 400 300 200 100 0 log M 14.53 14.33 14.13 13.93 13.73 | 13.53 13.33 13.13 12.93 12.73 12.53 12.33 12.13 11.93 11.73 800 700 600 500 400 300 200 100 0 log M 14.53 14.33 14.13 13.93 13.73 | 13.53 13.33 13.13 12.93 12.73 12.53 12.33 12.13 11.93 11.73 σ 12 , t/ 2(kms−1 ) σ12, r (kms −1 ) v12, r (kms−1 ) 800 700 600 500 400 300 200 100 0 0.1 0.1 1.0 r h −1 10.0 1.0 r h −1 10.0 1.0 10.0 ( Mpc) ( Mpc) 0.1 r ( Mpc) −1 h Figure 2.8: Pairwise velocity and velocity dispersions of haloes. The left panel shows the scaledependent pairwise radial velocity as a function of halo mass. The middle and right panels are similar, but for pairwise radial and transverse velocity dispersions, respectively. 30 same large-scale gravitational potential field sourced by linear fluctuations. The pairwise radial velocity dispersion σ12,r (middle panel) decreases toward small halo pair separation. The dependence on halo mass is weak – in the range of r ∼ 1–5h−1 Mpc, lower mass haloes have slightly higher σ12,r . The dispersion σ12,t of the pairwise transverse velocity (right panel) has a similar trend. On large scales, the one-dimensional (1D) pairwise transverse velocity √ dispersion σ12,t / 2 (right panel) is about 10 percent lower than the radial one (middle panel). In Fig. 2.9, the dependence of pairwise radial velocity (top panel) and the radial and transverse velocity dispersions (middle and bottom panels) on assembly variables are shown for haloes of log[Mh /(h−1 M )] = 11.73. Below ∼ 1h−1 Mpc, the pairwise radial velocity does not show a strong dependence on any of the assembly variables. Above ∼ 1h−1 Mpc, the assembly effect becomes clear. Both radial and transverse pairwise velocity dispersions show a substantial dependence on assembly variables at all scales. To see the trend more clearly, Fig. 2.10 shows the dependences by grouping each assembly variable into four percentiles. For the radial pairwise velocity, the bump around r ∼ 5h−1 Mpc seen in Fig. 2.8 shows up for each group, and the scatter caused by assembly effect also reaches maximum at this scale. The assembly effect in pairwise velocity with halo spin is not as strong as with other assembly variables, opposite to the trend seen in halo bias (e.g. left panel of Fig. 2.4). For low-mass haloes presented here, similar to the spatial clustering, the trend of the dependence of pairwise velocity on assembly variables follows the correlations among assembly variables, except for the case with halo spin. At fixed mass, haloes that on average form earlier (with lower aM/2 , higher c, or higher Vpeak ) have higher pairwise velocity and velocity dispersions. This is in line with the evolution of those low-mass haloes being influenced by the surrounding environment, especially the tidal field (e.g. Hahn et al. (2007, 2009); Wang et al. (2011); Shi et al. (2015)) and with some of them being ejected haloes around massive haloes (e.g. Wang et al. (2009); Wetzel et al. (2014)). The deviation between the pairwise velocity trend with halo spin and the expectation from correlations with other assembly variables arises from the misalignment of the gradient of pairwise velocity and the direction of the correlation in the plane of spin and one other assembly variable, similar to what we discuss in § 2.4.4 for halo bias. For the joint dependence of pairwise velocity and velocity dispersions on halo mass and each √ assembly variable, we plot v12,r , σ12,r , and σ12,t / 2 at the scale r ∼ 5h−1 Mpc in Fig. 2.11. The trend is similar to the joint dependence for halo bias (Fig. 2.1), and at fixed mass, more strongly clustered 31 100 50 0 0.1 500 logM=11.73 400 300 200 100 r(h −1 Mpc) 1.0 10.0 300 200 150 v12, r 100 log c 0.38 0.53 0.67 0.82 | 0.97 1.12 1.27 1.43 400 300 200 100 50 0 0.1 1.0 ( Mpc) r h −1 150 v12, r 100 50 0 2(kms−1 ) , t/ 12 σ 500 logaM/2 -0.55 -0.51 -0.47 -0.43 -0.39 | -0.35 -0.31 -0.27 -0.23 -0.19 -0.15 -0.11 -0.07 -0.03 0.01 0.1 -2.05 -1.95 -1.85 -1.75 -1.65 -1.55 -1.45 -1.35 -1.25 | -1.15 0.1 logM=11.73 400 300 200 100 1.0 ( Mpc) r h −1 10.0 0 0 10.0 400 300 200 logM=11.73 r(h −1 Mpc) 1.0 100 10.0 0 peak 2.05 2.07 2.09 2.11 2.13 2.15 2.17 2.19 2.21 2.23 2.25 2.27 2.29 | 2.31 2.33 2.35 V 0.1 logM=11.73 1.0 − r h ( 1 Mpc) 10.0 logaM/2 -0.55 -0.51 -0.47 -0.43 -0.39 | -0.35 -0.31 -0.27 -0.23 -0.19 -0.15 -0.11 -0.07 -0.03 logM=11.73 r(h −1 Mpc) 0.1 1.0 10.0 500 400 log c 0.38 0.53 0.67 0.82 | 0.97 1.12 1.27 1.43 300 200 logM=11.73 0.1 1.0 r(h −1 Mpc) 0 500 0.1 300 200 logM=11.73 r(h −1 Mpc) 1.0 10.0 10.0 log λ 400 1.0 ( −1 Mpc) 0.1 r h logλ -2.05 -1.95 -1.85 -1.75 -1.65 -1.55 -1.45 -1.35 -1.25 | -1.15 log c 0.38 0.53 0.67 0.82 | 0.97 1.12 1.27 1.43 100 10.0 500 logλ σ12, r (kms −1 ) (kms−1 ) 200 0 10.0 300 250 1.0 r(h −1 Mpc) 500 logM=11.73 σ12, r (kms −1 ) (kms−1 ) 250 0 0.1 100 σ12,t / 2(kms−1 ) 150 -0.55 -0.51 -0.47 -0.43 -0.39 | -0.35 -0.31 -0.27 -0.23 -0.19 -0.15 -0.11 -0.07 -0.03 0.01 10.0 logM=11.73 2(kms−1 ) 200 logaM/2 1.0 ( Mpc) r h −1 0 200 , t/ v12,r (kms−1 ) 250 100 300 log 12 300 0.1 200 400 σ 0 300 2.05 2.07 2.09 2.11 2.13 2.15 2.17 2.19 2.21 2.23 2.25 2.27 2.29 | 2.31 2.33 2.35 2(kms−1 ) 50 400 500 logVpeak , t/ v12, r 100 500 logM=11.73 100 12 150 2.05 2.07 2.09 2.11 2.13 2.15 2.17 2.19 2.21 2.23 2.25 2.27 2.29 | 2.31 2.33 2.35 σ 200 logVpeak σ12, r (kms −1 ) (kms−1 ) 250 σ12,r (kms−1 ) 300 0 -2.25 -2.15 -2.05 -1.95 -1.85 -1.75 -1.65 -1.55 -1.45 -1.35 -1.25 -1.15 | -1.05 0.1 logM=11.73 1.0 ( −1 Mpc) 10.0 r h Figure 2.9: Dependence of halo pairwise velocity and velocity dispersions on assembly variables for low-mass haloes (log[Mh /(h−1 M )] = 11.73). Right, middle, and light panels show pairwise radial velocity, pairwise radial velocity dispersion, and pairwise transverse velocity dispersion, respectively. In each panel, the curves are color-coded according to the value of the corresponding assembly variable. 32 (kms−1 ) 250 200 1.0 logM=11.73 c dependence ( Mpc) r h −1 400 1.0 logM=11.73 λ dependence ( Mpc) r h −1 10.0 2(kms−1 ) , t/ 12 σ σ12,t / 2(kms−1 ) 2(kms−1 ) 1.0 logM=11.73 c dependence 1.0 0 0.1 500 logM=11.73 λ dependence 400 1.0 − r h lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 ( 1 Mpc) logM=11.73 λ dependence 300 200 100 0 0.1 100 r(h −1 Mpc) 200 50 r(h −1 Mpc) 200 lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 300 100 v12, r 400 lower-middle 25% upper-middle 25% upper 25% lower 25% 10.0 300 0 0.1 500 lower 25% lower-middle 25% upper-middle 25% upper 25% 0 0.1 500 logM=11.73 c dependence 100 10.0 150 0 0.1 1.0 , t/ v12, r 50 100 r(h −1 Mpc) 200 100 300 lower-middle 25% upper-middle 25% upper 25% lower 25% 10.0 300 150 0 0.1 200 0 0.1 400 logM=11.73 aM/2 dependence 300 100 500 lower 25% lower-middle 25% upper-middle 25% upper 25% lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 ( 1 Mpc) 12 200 400 1.0 − r h σ (kms−1 ) 250 500 logM=11.73 aM/2 dependence 200 logM=11.73 aM/2 dependence 1.0 10.0 r(h −1 Mpc) 0 0.1 2(kms−1 ) 300 1.0 100 r(h −1 Mpc) 300 100 0 0.1 0 0.1 400 150 50 100 500 lower-middle 25% upper-middle 25% upper 25% lower 25% 200 lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 1.0 lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 r(h −1 Mpc) , t/ 200 300 100 12 v12,r (kms−1 ) 250 logM=11.73 Vpeak dependence 1.0 10.0 r(h −1 Mpc) logM=11.73 peak dependence V σ v12, r 300 400 200 100 0 0.1 500 logM=11.73 Vpeak dependence 300 150 50 σ12, r (kms −1 ) 400 σ12,r (kms−1 ) 200 500 lower 25% lower-middle 25% upper-middle 25% upper 25% σ12, r (kms −1 ) (kms−1 ) 250 σ12, r (kms −1 ) 300 0 0.1 1.0 − r h lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 ( 1 Mpc) Figure 2.10: Same as Fig. 2.9, but in each panel, the curves are color-coded according to the percentile of the corresponding assembly variable. 33 σ ,r (km s− ) log[V /(km s − )] 12.0 Mh /(h − M ⊙ )] 12.5 13.0 log[ 500 13.5 1 14.0 14.5 0.0 13.0 log[ 13.5 1 14.0 14.5 1.8 12.0 Mh /(h − M ⊙ )] 12.5 13.0 log[ 13.5 1 14.0 14.5 1.8 450 150 1.0 0.8 200 0.6 150 12.0 Mh /(h − M ⊙ )] 12.5 13.0 log[ 13.5 1 14.0 14.5 2 aM/ 12 log 400 350 1.2 0 12.0 Mh /(h − M ⊙ )] 12.5 13.0 log[ 13.5 1 14.0 14.5 500 50 0 12.0 1 14.0 14.5 0 1 14.0 14.5 0 450 12 350 300 λ −1.5 250 200 −2.0 150 σ 150 12 200 −2.0 400 log ,r (km s −1 ) 250 v Mh /(h − M ⊙ )] log[ 13.5 13.5 500 300 λ λ log −1.5 log ,r (km s −1 ) 350 100 13.0 13.0 −1.0 350 200 12.5 Mh /(h − M ⊙ )] 12.5 log[ 400 250 12.0 100 0.2 400 50 150 0.4 450 300 −2.5 200 0.6 −1.0 150 250 0.8 500 450 −2.0 300 1.0 50 0.2 −1.0 −1.5 450 1.4 100 0.4 0 500 1.6 250 50 0.2 14.5 1.8 300 c 100 0.4 1 14.0 σ 0.6 13.5 350 1.2 log 200 13.0 400 v 0.8 Mh /(h − M ⊙ )] 12.5 log[ 1.4 12 250 50 12.0 ,r (km s −1 ) ,r (km s −1 ) c 300 0 450 350 1.0 100 −0.6 500 1.6 400 1.2 200 −0.5 50 −0.6 0 1.4 250 −0.4 100 500 1.6 300 150 σ 150 350 −0.3 c Mh /(h − M ⊙ )] 12.5 400 log 12.0 500 −0.2 200 −0.4 −0.5 50 0 2 (km s −1 ) ,r (km s −1 ) 2 aM/ 12 250 100 −0.6 14.5 450 300 −0.3 v 150 1 14.0 0.0 350 −0.2 200 −0.5 13.5 400 log 250 −0.4 13.0 −0.1 12 ,r (km s −1 ) 2 aM/ 300 −0.3 Mh /(h − M ⊙ )] 12.5 450 350 −0.2 50 12.0 log[ −0.1 400 log 2.0 500 450 −0.1 log 0 100 −2.5 50 12.0 Mh /(h − M ⊙ )] 12.5 log[ 13.0 13.5 1 14.0 14.5 0 12 2.0 σ ,t / 0.0 0 2 (km s −1 ) 14.5 12 1 14.0 σ ,t / log[ 13.5 100 50 2 (km s −1 ) 13.0 150 2.2 100 −2.5 12 Mh /(h − M ⊙ )] 12.5 σ ,t / 100 2.2 200 12 1 150 250 2.4 σ ,t / 12.0 peak peak 200 2.4 12 1 12 250 50 2.0 2 (km s −1 ) 1 1 1 100 2.2 300 2.6 250 150 350 300 2.6 200 400 2.8 350 300 2.4 450 400 2.8 350 500 3.0 450 400 2.6 peak 3.0 450 2.8 log[ 500 v ,r (km s− ) log[V /(km s − )] V /(km s− )] 500 3.0 50 12.0 Mh /(h − M ⊙ )] 12.5 log[ 13.0 13.5 1 14.0 14.5 0 Figure 2.11: Dependence of halo pairwise velocity and velocity dispersions on halo mass and each of the assembly variables. left, middle, and right panels show pairwise radial velocity, pairwise radial velocity dispersion, and pairwise transverse velocity dispersion, respectively. The pairwise velocity and velocity dispersions are evaluated at r ∼ 10h−1 Mpc. The four black curves in each panel mark the central 50 and 80 percent of the distribution of the corresponding assembly variable as a function of halo mass. 34 haloes display higher pairwise velocity and velocity dispersions. 2.5.2 Relation between spatial clustering and pairwise velocity under assembly effect The overall pattern of pairwise velocity and velocity dispersions in Fig. 2.11 shows a clear correlation with that of halo bias in Fig. 2.1. With the joint effect of halo mass and assembly, more strongly clustered haloes tend to move faster with higher velocity dispersions. Such a correlation motivates us to further study the relation between spatial clustering and halo velocity field. Under the assumption that halo pairs are conserved during evolution, a relation can be established between halo pairwise velocity and spatial clustering. On large scales, halo pairwise velocity is found to be proportional to halo bias (e.g. Sheth et al. (2001); Zhang & Jing (2004)), v12,r ∝ b. This would have interesting implications for the assembly effect. For example, in the low halo mass regime, if all older haloes originate from stripped massive haloes so that their stronger clustering is a manifestation of that of more massive haloes, their pairwise velocity should also follow that of the massive haloes. We would expect a tight correlation between v12,r and b for all subsets of haloes (divided according to mass and assembly history). In the top panels of Fig. 2.12, the values of v12,r around r ∼ 10h−1Mpc and halo bias b are plotted for halo subsamples in bins of halo mass and assembly variables. In each panel, the color indicates the halo mass and the size of the points denotes the magnitude of the assembly variable. We note that the magnitude of the assembly bias is not the focus here. At fixed mass, it is simply related to the value of the bias factor in a way seen in Fig. 2.1. The plot is meant to show the v12,r–b sequence induced by assembly effect within each halo mass bin (represented by one color) and as a collection from different mass bins. In each panel, the curve is the relation derived from the dependence on halo mass only, v12,r(M) versus b(M), which is close to the expected linear relation. On average, the v12,r–b relation from the points tracks the curve quite well. At fixed mass, more strongly clustered haloes originated from assembly effect tends to have pairwise velocity shifted toward the curve (e.g. the blue points with the lowest mass). However, the magnitude of the shift is not large enough to make the points fall on top of the curve. That is, within each halo mass bin, the v12,r–b relation from assembly effect shows deviations from the mass-only curve. In the v12,r–b relation as a whole from different mass bins, the assembly effect causes the scatter around the mass-only relation. The deviation and scatter from the expected v12,r–b v12, r (km s −1 ) v12, r (km s −1 ) 800 700 600 500 400 300 200 100 0 0.0 800 700 600 500 400 300 200 100 0 0.0 1.5 b 2.0 2.5 dependence __ cmass-only dependence 0.5 1.0 1.5 b 2.0 2.5 dependence __ mass-only dependence 1.0 1.5 b 2.0 log M 14.0 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 2.5 2.0 2.5 aM/ dependence __ mass-only dependence 1.0 1.5 b 2.0 2.5 dependence __ cmass-only dependence 0.5 1.0 1.5 b 2.0 2.5 dependence __ mass-only dependence 1.0 1.5 b 2.0 log M 14.0 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 3.0 800 700 600 500 400 300 200 100 0 0.0 2 (km s −1 ) σ12, t / σ12, r (km s −1 ) 3.0 800 700 600 500 400 300 200 100 0 0.0 log M 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 λ 0.5 3.0 800 700 600 500 400 300 200 100 0 0.0 log M 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 2 0.5 3.0 2 (km s −1 ) b 12 3.0 800 700 600 500 400 300 200 100 0 0.0 1.5 2.5 dependence __ mass-only dependence log M 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 Vpeak 0.5 1.0 1.5 b 2.0 2.5 M/ dependence __ amass-only dependence 0.5 1.0 1.5 b 2.0 2.5 dependence __ cmass-only dependence 0.5 1.0 1.5 b 2.0 2.5 dependence __ mass-only dependence 1.0 1.5 b 2.0 3.0 log M 14.0 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 3.0 log M 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 λ 0.5 3.0 log M 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 2 σ ,t / 3.0 800 700 600 500 400 300 200 100 0 0.0 log M 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 λ 0.5 3.0 1.0 2 (km s −1 ) __ mass-only dependence 1.0 800 700 600 500 400 300 200 100 0 0.0 log M 13.8 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 2 0.5 σ12, t / aM/ dependence 0.5 3.0 2 (km s −1 ) 2.5 800 700 600 500 400 300 200 100 0 0.0 log M 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 σ12, t / 2.0 ,r (km s −1 ) b 12 1.5 dependence __ mass-only dependence Vpeak σ 1.0 v 12 0.5 800 700 600 500 400 300 200 100 0 0.0 log M 13.6 13.4 13.2 13.0 12.8 12.6 12.4 12.2 12.0 11.8 11.6 σ12, r (km s −1 ) 800 700 600 500 400 300 200 100 0 0.0 dependence __ mass-only dependence Vpeak σ12, r (km s −1 ) 800 700 600 500 400 300 200 100 0 0.0 ,r (km s −1 ) v12, r (km s −1 ) 35 2.5 3.0 Figure 2.12: Relation between pairwise velocity and velocity dispersions and halo bias. Left, middle, and right panels show pairwise radial velocity, pairwise radial velocity dispersion, and pairwise transverse velocity dispersion versus halo bias, respectively. In each panel, the points are color-coded by halo mass and their size indicates the value of the corresponding assembly variable. The pairwise velocity and velocity dispersions are evaluated at r ∼ 10h−1 Mpc. The solid curve is computed from mass-only dependence of halo clustering. 36 relation implies that the environment plays a role more complicated than our naive expectation. In the middle and bottom panels of Fig. 2.12, the pairwise radial and transverse velocity dispersions show a considerably large scatter around the mass-only curve. The mass-only curve is almost flat, which means that for haloes selected by mass, the pairwise velocity dispersions are nearly independent of halo mass. Such a behaviour can be explained in linear theory by associating haloes with smoothed perturbation peaks (e.g. Bardeen et al. (1986); Sheth & Diaferio (2001)). Once haloes are split by their assembly history, the pairwise velocity dispersions no longer follow the mass-only curve, especially for low mass haloes. For example, the velocity dispersions for the lowest mass haloes (blue points) from assembly effect almost monotonically increase with increasing bias, substantially deviating from the mass-only curve. Low-mass haloes formed earlier (or with higher Vpeak or higher concentration) can have 1D pairwise velocity dispersions as high as 600–800 km s−1 , much higher than the mass-only value (∼300 km s−1 ). Interestingly, they are even higher than those of massive haloes. This is consistent with those haloes being in denser environment, where larger pairwise velocity dispersions are expected (Sheth & Diaferio (2001)). The other possible origin of the high velocity dispersion is that some low-mass haloes are stripped haloes ejected from massive haloes. The above results indicate that assembly effect in the kinematics of haloes, in particular in the pairwise velocity dispersions, does not follow a simple description and that the environment can play an important role in shaping it. As the assembly bias for low mass haloes is shown to be connected to their tidal environment (e.g. Hahn et al. (2009); Borzyszkowski et al. (2017); Paranjape et al. (2018)), studying the dependence of pairwise velocity statistics on the tidal environment (e.g. characterised by the tidal tensor) would be a useful step to help understand the assembly effect on the halo kinematics. The assembly effect on halo kinematics can affect the redshift-space clustering of haloes, raising the possibility of using redshift-space clustering to detect assembly effect with galaxy clustering. We analyse the redshift-space two-point correlation function and the multiple moments (monopole, quadrupole, and hexadecapole) of halo samples in mass and assembly variable bins. We find that the assembly effect on redshift-space clustering is mainly through the bias factor and that halo kinematics only introduce a minor effect. Similar to the scale-dependence case (§ 2.4.3), the halo assembly effect may be hard to be revealed in practice from redshift-space clustering, but it deserves a further, detailed study. 37 2.6 Summary and discussion Using dark matter haloes identified in the MultiDark MDR1 simulation, we investigate the joint dependence of halo bias on halo mass and halo assembly variables and the assembly effect on halo kinematics. For the halo assembly variables (peak maximum circular velocity Vpeak , halo formation scale factor aM/2 , halo concentration c, and halo spin λ) considered in this work, the joint dependence of halo bias on halo mass and each assembly variable can be largely described as that halo bias increases outward from a point of global minimum in the plane of mass and assembly variable. All previous results of halo assembly bias measured for certain percentiles of a given halo assembly variable can be inferred from the above dependence and the distribution of assembly variable. We explore the possibility of finding an effective halo variable to minimise the assembly bias by using a combination of halo mass and spin. While an effective halo mass constructed through the combination absorbs the dependence of assembly bias on halo spin, at fixed effective mass, halo bias still depends on other assembly variables. The investigation indicates that assembly bias is multivariate in nature and that it is unlikely for one halo variable to absorb every aspect of the assembly effect. From studying the joint dependence of halo bias on two assembly variables at fixed halo mass, we find that it is not necessarily true to predict the trend of assembly bias for one assembly variable solely based on that in the other assembly variable and their correlation. It only becomes possible if the gradient in halo bias follows the correlation direction. Whether the gradient and correlation directions align with each other relies on which two assembly variables to choose. It also depends on halo mass – with respect to the correlation direction, the gradient direction rotates as halo mass varies, showing the non-trivial nature of assembly bias. We also study the kinematics of haloes under the assembly effect by dividing haloes according to halo mass and assembly variables. In general, more strongly clustered haloes have higher pairwise radial velocity and higher pairwise velocity dispersions. For low-mass haloes showing higher bias caused by assembly effect, the pairwise radial velocity tends to approach that of massive haloes of similar clustering amplitude, while the pairwise velocity dispersions can be substantially higher than those of the massive haloes. The results supports the picture that the evolution of low mass haloes is influenced by the surrounding environment, especially the tidal field, and that some of low mass haloes could be ejected haloes around massive haloes. However, we do not find a simple description for the relation between halo kinematics and spatial clustering under assembly effect. 38 The assembly bias for low mass haloes (Mh / Mnl ) has a substantial scale dependence, showing as a drop on scales below ∼ 3h−1 Mpc. The scale dependence of assembly bias and the assembly effect on halo kinematics can potentially provide an approach to identify assembly effects in galaxy clustering data through the shape of the scale-dependent galaxy bias and redshift-space distortions. However, the effect can be subtle, which can be masked by the one-halo term of galaxy clustering and the scatter between the galaxy properties and halo assembly. Further study is needed to see how and whether this approach works with high precision galaxy clustering data by incorporating the description of galaxy-halo connection with assembly effect included. CHAPTER 3 GALAXY ASSEMBLY BIAS OF CENTRAL GALAXIES IN THE ILLUSTRIS SIMULATION This chapter originally appeared in a work that has been submitted to the MNRAS (Xu & Zheng (2018b)). Sections and figures have been renumbered and reformatted to match the dissertation. The coauthor of the paper is Zheng Zheng. 3.1 Abstract Galaxy assembly bias, the correlation between galaxy properties and halo properties at fixed halo mass, could be an important ingredient in halo-based modelling of galaxy clustering. We investigate the central galaxy assembly bias by studying the relation between various galaxy and halo properties in the Illustris hydrodynamic galaxy formation simulation. Galaxy stellar mass M∗ is found to have a tighter correlation with peak maximum halo circular velocity Vpeak than with halo mass Mh . Once the correlation with Vpeak is accounted for, M∗ has nearly no dependence on any other halo assembly variables. The correlations between galaxy properties related to star formation history and halo assembly properties also show a cleaner form as a function of Vpeak than as a function of Mh , with the main correlation being with halo formation time and to a less extent halo concentration. Based on the galaxy-halo relation, we present a simple model to relate the bias factors of a central galaxy sample and the corresponding halo sample, both selected based on assembly-related properties. It is found that they are connected by the correlation coefficient of the galaxy and halo properties used to define the two samples, which provides a reasonable description for the samples in the simulation and suggests a simple prescription to incorporate galaxy assembly bias into the halo model. By applying the model to the local galaxy clustering measurements in Lin et al. (2016), we infer that the correlation between star formation history or specific star formation rate and halo formation time is consistent with being weak. 40 3.2 Introduction It has been well established that galaxies form in dark matter haloes (White & Rees (1978)). As the first step to study galaxy formation and clustering, halo formation and clustering, which is dominated by gravity, have been extensively studied with analytic models (Press & Schechter (1974); Bardeen et al. (1986); Mo & White (1996); Sheth & Tormen (1999)) and cosmological N-body simulations (Springel et al. (2005); Prada et al. (2012)). It has been found that halo clustering depends not only on halo mass but also on halo assembly history or environment (Gao et al. (2005); Gao & White (2007); Paranjape et al. (2018); Xu & Zheng (2018a); Han et al. (2019)). This is called halo assembly bias, whose nature is still under investigation (Dalal et al. (2008); Castorina & Sheth (2013)). If galaxy properties are affected by halo formation and assembly history, halo assembly bias would translate to galaxy assembly bias. Operationally, galaxy assembly bias can be defined as that at fixed halo mass, the statistical galaxy content shows dependence on other halo variables or galaxy properties show correlations with halo assembly history. The widely adopted halo model (Cooray & Sheth (2002)) of interpreting galaxy clustering, such as the halo occupation distribution (Berlind & Weinberg (2002); Zheng et al. (2005)) or conditional luminosity function (Yang et al. (2003)), makes the implicit assumption of no galaxy assembly bias. Such methods have been successfully applied to galaxy clustering (Zehavi et al. (2005); Zheng et al. (2007); Xu et al. (2018)). However, if assembly bias is significant, neglecting it in the model would lead to incorrect inference of galaxy-halo connections and introduce possible systematics in cosmological constraints (e.g. Zentner et al. (2014, 2019); but see also McEwen & Weinberg (2016); McCarthy et al. (2018)). Conversely, observationally inferred galaxy assembly bias would help understand galaxy formation. The existence and strength of galaxy assembly bias are still a matter far from settled, either in theory or in observation. Galaxy assembly bias has been investigated in hydrodynamic or semi-analytic galaxy formation models (Berlind et al. (2003); Croton et al. (2007); Mehta (2014); Chaves-Montero et al. (2016); Zehavi et al. (2018); Contreras et al. (2019)), focusing on the effect on galaxy occupation function and galaxy clustering. The results seem to depend on the implementation details of star formation and feedback. Studying galaxy assembly bias from observation has the difficulty of determining halo mass, and the results are not conclusive (Yang et al. (2006); Berlind et al. (2006); Lin et al. (2016); Zu et al. (2017); Guo et al. (2017)). Given the potential importance of galaxy assembly bias in modelling galaxy clustering, in this paper we study the correlation between various 41 central galaxy properties and halo properties in the Illustris hydrodynamic simulation (Vogelsberger et al. (2014b)) at the halo level, aiming at providing useful insights in describing galaxy assembly bias. The structure of the paper is as follows. In Section 3.3, we introduce the simulation and the galaxy and halo catalogues. Then in Section 3.4, we investigate the relation between galaxy and halo properties, with primary galaxy-halo properties in Section 3.4.1 and general galaxy-halo properties in Section 3.4.2. In Section 3.4.3, we present a simple model to connect galaxy assembly bias with halo assembly bias. Finally, we summarise and discuss the results in Section 3.5. 3.3 Simulation and galaxy-halo catalogue In this work, we use galaxies and haloes from the state-of-the-art hydrodynamic galaxy formation simulation Illustris1 (Vogelsberger et al. (2014b); Nelson et al. (2015)) to study galaxy assembly bias, which is able to produce different type of galaxies seen in observation (Vogelsberger et al. (2014a); Genel et al. (2014)). In particular, we use the Illustris-2 simulation, which has a box size of 75h−1 Mpc on a side, and contains 9103 dark matter particles of mass 5 × 107 M and the same number of baryon particles of mass 1 × 107 M . The mass resolution is sufficient for our purpose of studying (central) galaxies in haloes that are more massive than a few times 1010 h−1 M . The simulation adopts a spatially-flat cosmology with the following parameters: Ωm = 0.27, Ωb = 0.0456, h = 0.70, ns = 0.963, and σ8 = 0.809. The haloes in the Illustris database are identified with the friends-of-friends (FoF) algorithm. As this algorithm is notorious for having the probability of bridging two separate haloes into one halo, we apply the phase-space halo finder Rockstar (Behroozi et al. (2013)) to identify haloes, with dark matter particles extracted from 30 snapshots of the simulation (from z=3.94 to z=0). We then build the merger tree using the consistent-tree algorithm (Behroozi et al. (2012)). We focus our study on the relation between central galaxies and dark matter haloes and the assembly effect at z = 0. Galaxies in the Illustris-2 simulation are assigned to the Rockstar haloes. For each galaxy, if its distance dgh to the centre of a halo is smaller than the virial radius Rvir of the halo, it is assigned to this halo. In rare cases, a galaxy may be assigned to more than one haloes, as haloes are not spherical. We then put the galaxy into the halo that corresponds to the lowest ratio of dgh /Rvir . With galaxies assigned to haloes, for each halo, we define the most massive galaxy inside 1 http://www.illustris-project.org 42 0.2Rvir from the halo centre to be the central galaxy. If there is no galaxy inside 0.2Rvir , we simply choose the most massive galaxy inside Rvir as the central galaxy. With the catalogue of central galaxies and the associated host haloes, we study the relation between galaxy properties and halo properties and the assembly effects. The halo properties we focus on are: (1) Mh , halo mass enclosed in a volume with mean density of 200 times the background density of the universe; (2) Vpeak , peak maximum circular velocity of the halo over its accretion history; (3) c, halo concentration parameter, defined as the ratio of halo virial radius to scale radius; (4) aM/2 , cosmic scale factor when the halo obtains half of its current (z = 0) total mass; (5) MÛ h , halo mass accretion rate near z = 0 (averaged between z=0 and z=0.197, about 2.4 Gyr, one dynamical time), in units of h−1 M yr−1 ; (6) MÛ h /Mh , specific halo accretion rate, in units of Gyr−1 . The central galaxy properties we consider include: (1) M∗ , stellar mass (sum of masses of star particles within twice the stellar half mass radius); (2) SFR, star formation rate within twice the stellar half mass radius; (3) sSFR, specific star formation rate, the ratio of SFR to M∗ ; (4) g − r, galaxy color defined by the g-band and r-band luminosity. 3.4 Results We aim at presenting the relation between central galaxies and haloes to learn about the correlation between halo formation and assembly and galaxy properties. In Section 3.4.1, we first study how galaxy stellar mass depends on the primary halo properties (Mh and Vpeak ). Then we investigate how various halo and galaxy properties are correlated in Section 3.4.2. Finally, in Section 3.4.3, we use a simplified model to describe the connection between galaxy and halo assembly bias factor. 43 3.4.1 Relationship between stellar mass and halo properties Galaxy stellar mass is a primary property inferred from observation. The relation between this primary galaxy property and certain primary halo property (e.g., Mh and Vpeak ) can be established based on subhalo abundance matching or modelling the stellar mass-dependent clustering, which encodes information about galaxy formation. Here we show the relation predicted by the Illustris simulation and study how tight stellar mass correlates with Mh and Vpeak . The top-left panel of Fig. 3.1 shows M∗ as function of Mh , color-coded with values of Vpeak . Galaxy stellar mass M∗ increases steeply with Mh at log[M∗ /(h−1 M )] < 12 and then slowly at log[M∗ /(h−1 M )] > 12, a trend similar to that inferred from observation (Behroozi et al. (2010); Leauthaud et al. (2012); Zu & Mandelbaum (2015)). The scatter in M∗ at fixed halo mass in the M∗ –Mh relation decreases with increasing Mh (solid curve in Fig. 3.2), varying from about 0.3 dex at the low-mass end to about 0.17 dex at the high-mass end. The scatter at the high-mass end is consistent with the value ∼ 0.16 dex inferred from galaxy clustering modelling (Tinker et al. (2017)). One source of the scatter can be the halo formation history (Tinker (2017)), which may affect the growth history of stellar mass (either from star formation or galaxy merging; e.g. Gu et al. (2016)). The color code in Vpeak in the top-left panel enables us to see how the scatter in the M∗ –Mh relation may be connected to halo assembly. On average, Vpeak and Mh are correlated, and the mean relation is found to be well described by Vpeak Mh = 170 12 10 h−1 M 1/3 km s−1 (3.1) in the Illustris simulation. However, there is scatter on top of the mean relation, and at fixed Mh , the distribution of Vpeak reflects that in the assembly history. As can be seen in the top-left panel, the assembly of haloes encoded in Vpeak does contribute to the scatter in the M∗ –Mh relation – at fixed Mh , galaxies residing in haloes of higher Vpeak tend to have higher stellar mass, especially at the low mass end. To see how well the scatter in M∗ can be attributed to the scatter in Vpeak , in the top-right panel of Fig. 3.1, we plot the M∗ –Vpeak relation. It follows a similar trend seen in the M∗ –Mh relation, steeper (shallower) dependence of M∗ on Vpeak at the low (high) Vpeak end, which is expected given the correlation between Vpeak and Mh . The M∗ –Vpeak relation appears to be tighter than the M∗ –Mh relation, in the sense that at fixed Vpeak , the scatter in M∗ is lower than that at the corresponding Mh (see Matthee et al. (2017) for a similar result in terms of z = 0 maximum halo circular velocity with 44 12 log[M * /(h −1M ⊙ ⊙] 11 logVpeak 1.81 1.87 1.93 1.99 2.05 2.11 2.17 2.23 2.29 2.35 2.44 2.65 10 9 8 7 10.0 10.5 11.0 11.5 12.0 12.5 h/( −1 ⊙ ⊙] log[M h M 13.0 13.5 14.0 11 log[M * /(h −1M ⊙ )] 12 logMh 10.6 10.8 11.0 11.2 11.4 11.6 11.8 12.0 12.2 12.4 12.7 13.7 10 9 8 7 1.8 2.0 2.2 2.4 2.6 log[Vpeak/(km s−1)] 2.8 3.0 Figure 3.1: Top-left: M∗ as function of Mh for central galaxies. The galaxies are color-coded according to log[Vpeak /(km s−1 )]. For galaxies in each bin of log Vpeak , the contours correspond to the 68.3 and 95.4 percent distribution, respectively. Top-right: M∗ as function of Vpeak for central galaxies, color-coded according to log[Mh /(h−1 M )]. Bottom-left: M∗ as function of Mh for central galaxies, with the mean relation color-coded according to the values of aM/2 . For clarity, the scatter in the mean relation is only shown for the bin with the highest aM/2 (latest forming haloes). Bottom-right: M∗ as function of Vpeak for central galaxies, color-coded according to aM/2 , with the shaded region illustrating the scatter for the bin with the highest aM/2 . Note the remarkable result that M∗ does not depends on aM/2 at fixed Vpeak (compared to the Mh case in the bottom-left panel). 45 log[Mh/(h −1M )] 10.6 11.0 11.4 11.8 12.2 12.6 13.0 13.4 13.8 14.2 Vpeak dependence Mh dependence 0.30 σlogM * 0.25 0.20 0.15 0.10 1.81 1.93 2.05 2.17 2.29 2.41 2.53 2.65 2.77 2.89 log[Vpeak/(kms−1)] Figure 3.2: Standard deviation in log M∗ as a function of Mh (solid) and Vpeak (dashed). The correspondence between Mh and Vpeak is from the mean relation Vpeak ∝ Mh1/3 in equation (3.1). the EAGLE simulation). The scatter varies from ∼0.28 dex at low Vpeak to ∼0.13 dex at high Vpeak (dashed curve in Fig. 3.2). In the M∗ –Vpeak plot (top-right panel), the contours are color-coded by Mh . Unlike the Mh case in the top-left panel, we find that M∗ does not show a clear dependence on Mh – at fixed Vpeak (i.e. by taking a vertical cut in the plot), M∗ in haloes of different Mh appears to follow similar mean and scatter. In the bottom-left panel, we show the effect of the other halo assembly property aM/2 on the M∗ –Mh relation. Each curve shows the mean M∗ –Mh relation for haloes in one aM/2 bin. The scatter around the mean relation is illustrated with the shaded region, only shown for the latest forming haloes to avoid crowdedness. There is a clear and substantial dependence of the M∗ –Mh relation on the assembly property aM/2 – at fixed Mh , haloes forming earlier (smaller aM/2 ) tend to host galaxies of higher M∗ . Such a trend is consistent with previous work based on the EAGLE simulation (Matthee et al. (2017)) and semi-analytic galaxy formation model (Zehavi et al. (2018)). Switching to the M∗ –Vpeak relation (bottom-right panel), we find that at fixed Vpeak there is no dependence of M∗ on the assembly property aM/2 , with the curves of different aM/2 bins all falling on top of each other. As further shown in Section 3.4.2, Vpeak is able to absorb the effect on M∗ from any other halo assembly variable. This is consistent with the results using z = 0 maximum halo circular velocity (Matthee et al. (2017)). We will discuss this remarkable result in Section 3.5. 46 The tighter correlation between M∗ –Vpeak , in comparison to M∗ –Mh , suggests that in galaxy clustering modelling, galaxy assembly bias effect can be partially accounted for by switching from a Mh –based model to a Vpeak -based model. That is, the galaxy-halo relation is parameterised as a function of Vpeak . This is in line with the finding by Chaves-Montero et al. (2016) using the EAGLE simulation, while their velocity quantity most strongly correlating with M∗ is slightly different, Vrelax , the highest value of the maximum circular velocity of a subhalo with a relaxation criterion imposed. Chaves-Montero et al. (2016) show that Vrelax is able to capture the majority (but not all) the galaxy assembly bias effect on galaxy clustering for stellar-mass-based samples. The tighter correlation between M∗ –Vpeak may also be the reason that when using subhalo abundance matching or its variants to model galaxy clustering, halo circular velocity-based models usually have good performances (Reddick et al. (2013); Guo et al. (2016)). In what follows, besides a further investigation of the dependence of stellar mass on halo assembly variables, we also extend the study to other galaxy properties. 3.4.2 Relationship between galaxy properties and halo properties With the 7 halo properties (Mh , Vpeak , c, aM/2 , λ, MÛ h , and MÛ h /Mh ) and the 4 central galaxy properties (M∗ , SFR, sSFR, and g − r color), we study the correlations between them. To aid the discussion, we also present the correlations among halo properties and those among galaxy properties. 3.4.2.1 At fixed Mh Fig. 3.3 shows the correlation between each pair of the halo and galaxy properties for central galaxies in haloes of a narrow mass bin, log[Mh /(h−1 M )] = 12.0 ± 0.1. In each contour panel, the contours indicate 68.3 and 95.4 percent of the distribution of the pair of properties. The number labelled in each panel is the Pearson correlation coefficient ρ of the two properties, indicating how strong the correlation is. It is calculated as ρ= hx yi − hxihyi , σx σy (3.2) where x and y denote the two properties, hi means average, and σx and σy are the standard deviations of x and y. The panel at the top of each column shows the marginalised distribution of the property labelled at the x-axis of the column. The panels with red contours (i.e. the top 6 rows and left 6 columns of contour panels) display 2.4 2.3 2.2 0.464 1.4 -0.029 0.656 -0.009 -0.610 -0.654 -0.064 -0.351 -0.380 0.382 0.088 -0.340 -0.332 0.358 0.226 -0.046 -0.413 -0.340 0.378 0.244 0.983 0.329 0.548 0.573 -0.402 -0.280 -0.344 -0.392 0.097 -0.060 0.109 0.164 -0.013 -0.034 -0.047 0.541 -0.070 -0.380 -0.193 0.420 0.142 0.153 0.165 0.079 0.881 0.107 0.442 0.208 -0.458 -0.184 -0.181 -0.200 -0.052 -0.734 9.5 10 .0 10 .5 −1 0 1 1.0 0.6 0.00 −0.25 −0.50 logλ logaM/2 logc logVpeak 47 −1 −2 Ṁ h 150 50 0.15 0.05 −0.05 10.5 10.0 9.5 logsSFR logSFR logM * Ṁ h/Mh −50 1 0 −1 −9 Ṁ h/Mh logM * -0.840 −1 1 −1 0 −9 0.2 5 0.5 0 0.7 5 5 .1 5 .0 .0 5 0 Ṁ h 0 −1 logλ − 0 logaM/2 −5 0 50 15 0 logc −2 1.0 logVpeak 1.4 −0 .50 −0 .25 0.0 0 logMh 0.6 .0 .1 2.2 2.3 2.4 12 11 .9 −11 0.75 0.50 0.25 12 g-r −10 logSFR logsSFR g-r Figure 3.3: Relation between each pair of galaxy and/or halo properties at log[Mh /(h−1 M )] ∼ 12. In each contour panel, the two contours show the central 68.3 and 95.4 percent of the distribution of the pair of properties. The panels with red contours (i.e. the top 6 rows and left 6 columns of contour panels) display the correlations between halo properties. Those with black contours (i.e. the bottom 4 rows and the left 7 columns of contour panels) show the correlations between galaxy and halo properties, and those with blue contours (i.e. the right 3 columns of contour panels) are for the correlations between pairs of galaxy properties. The number in each contour panel is the Pearson correlation coefficient for the pair of properties. The histogram at the top panel of each column is the probability distribution function of the variable of that column. 48 the correlation between halo properties. Within the small but finite halo mass bin, all but one halo property show almost no correlation with Mh (correlation coefficient close to zero). The exception is Vpeak , and the correlation is simply driven by the Vpeak ∝ Mh1/3 mean relation. At fixed Mh , any pair of halo properties shows a substantial correlation (with | ρ| above 0.2). The nearly perfect correlation (ρ = 0.983) between MÛ h and MÛ h /Mh is a consequence of fixed Mh . Overall, the correlation trend is that haloes of higher Vpeak are more concentrated, form earlier, spin more slowly, and have lower accretion rate, which have been seen in previous work (Jeeson-Daniel et al. (2011); Han et al. (2019); Xu & Zheng (2018a)). The panels with black contours (i.e. the bottom 4 rows and the left 7 columns of contour panels) show the correlation between halo and galaxy properties. The correlation between M∗ and Mh shows up because of the finite size of the halo mass bin. At fixed Mh , the central galaxy stellar mass M∗ correlates with all other halo properties – haloes of higher Vpeak , higher concentration, earlier formation, lower spin, and lower accretion rate tend to host more massive central galaxies. On the contrary, the SFR shows no strong correlation with any halo properties. The most significant one is with halo formation time (ρ ∼ 0.16), with on average higher SFR in haloes of later formation. Given the substantial correlation between M∗ and halo properties and the weak or lack of correlation between SFR and halo properties, the sSFR (≡ SFR/M∗ ) is expected to correlate well with halo properties, but in an trend opposite to and weaker than that with M∗ . This is indeed the case. The most significant correlation is with Vpeak or aM/2 (both with | ρ| ∼ 0.4). The correlation between sSFR and the average halo accretion rate over the past dynamic time is there but not strong (ρ ∼ 0.16). The correlation between galaxy color g − r and halo properties essentially follows the case of sSFR (with a sign change in ρ; redder galaxies having lower sSFR). The panels with blue contours (i.e. the right 3 columns of contour panels), the correlations between pairs of galaxy properties at fixed halo mass, are shown. The SFR positively correlates with M∗ (ρ ∼ 0.54), and the mean relation has a slope close to unity, SFR ∝ M∗ . It resembles the star-forming main sequence (Brinchmann et al. (2004); Speagle et al. (2014); Santini et al. (2017)). That is, even if we only consider central galaxies in haloes of fixed mass, the star-forming main sequence emerges. Note that this SFR–M∗ correlation is not driven by the correlation of SFR and M∗ with a common halo variable we consider here. In fact, from Fig. 3.3, it can be seen that their correlation with a common halo variable may lead to the opposite effect. For example, haloes of earlier formation tend to host central galaxies of higher M∗ and lower SFR, and naively, this would 49 imply an anti-correlation between SFR and M∗ , opposite to what is found here. While it is possible that the halo-level star-forming main sequence is related to a halo variable not considered here, it is more likely that the sequence is driven by baryonic physics, which may have complicated dependence on or decouple from halo formation history. Unlike the SFR, the sSFR shows little dependence on M∗ . However, the sSFR is tightly correlated with the SFR (ρ ∼ 0.88). Given that sSFR≡SFR/M∗ , the pattern in the mutual correlations among M∗ , SFR, and sSFR can be achieved if the SFR–M∗ correlation coefficient is close to the ratio of the scatters in log M∗ and log SFR2, which appears to be the case. Galaxy g − r color strongly correlates with sSFR and follows the same trends as sSFR in its correlations with M∗ and SFR. Overall, at fixed halo mass, galaxy properties other than SFR show significant correlations with one or more halo properties, manifesting galaxy assembly bias at the level of haloes. The correlations among galaxy properties, however, may largely result from baryonic physics, given that the trend cannot be simply explained by their correlation with halo properties. In Section 3.4.1, it is found that switching from Mh to Vpeak can remove the dependence of M∗ on other halo properties. We now extend the investigation to other galaxy properties. 3.4.2.2 At fixed Vpeak Fig. 3.4 is similar to Fig. 3.3, but the correlations are presented for haloes at a fixed Vpeak bin, log[Vpeak /(km s−1 )] = 2.23 ± 0.03. The correlations among halo properties (in panels with red contours) are similar to those in Fig. 3.3, and there are additional correlations between Mh and other halo properties. The galaxy-halo correlations are shown in panels with black contours. The finite bin size in Vpeak makes the M∗ –Vpeak correlation show up. Other than this (and the one with Mh ), M∗ does not correlate with any other halo properties at fixed Vpeak , reinforcing the result in Section 3.4.1. It indicates that the correlations of M∗ with halo properties seen at fixed Mh (Fig. 3.3) can be attributed to the M∗ –Vpeak correlation and the correlation of Vpeak with other halo properties. For the SFR, at fixed Vpeak , it correlates significantly with Mh and aM/2 , higher SFR in haloes of higher mass and later formation. The correlations between SFR and other halo properties are weak. The correlations between sSFR (or color) and halo properties closely follow the SFR case. 2 To see this, let x = log M , y = log SFR, and z = log sSFR = y − x. We can derive the relation among the correlation ∗ coefficients, ρ xz /ρyz = (ρ xy − σx /σy )/(1 − ρ xy σx /σy ). For ρ xz to be near zero, we have ρ xy ∼ σx /σy . logMh logVpeak logc logaM/2 logλ Ṁ h Ṁ h/Mh -0.333 -0.761 −1 0 1 5 5 -0.169 9.5 10 .0 10 .5 .1 .0 -0.189 0 -0.112 logM * −1 1 −1 0 −9 0.383 logSFR logsSFR −1 1 −1 0 −9 0.2 5 0.5 0 0.176 −1 0 1 −5 0 50 15 0 −0 .05 0.0 5 0.1 5 9.5 10 .0 10 .5 0.089 −1 0 1 0.205 0 -0.468 0.156 5 0.228 0.077 5 5 5 0 10 .5 .0 9. 5 .1 .0 .0 0 0 5 10 0 0 − 0 1 5 − 5 0.011 .0 0.511 0.050 − 0 -0.255 -0.014 −5 0 50 15 0 −0 .05 0.0 5 0.1 5 9.5 10 .0 10 .5 0.395 −5 0 50 15 0 −0 .05 0.0 5 0.1 5 −1 −2 0.175 −5 0 50 15 0 -0.179 −1 0.044 −2 0.018 −1 0.405 −2 -0.327 −1 −0 .50 −0 .25 0.0 0 1. 5 1. 0 0 5 0 1 5 5 0 0 0 − 1 − 2 .0 .2 .5 .5 .0 − 5 0 − 0 − 0 1 1 0.470 −2 1. −0 5 .50 −0 .25 0.0 0 1.0 -0.435 −1 −0 .50 −0 .25 0.0 0 1. −0 5 .50 −0 .25 0.0 0 .5 5 0 0 5 0 − 1 − 2 0 .2 .5 .0 − 0 − 0 1.5 1.0 5 0.5 2.2 -0.300 −2 0.056 −0 .50 −0 .25 0.0 0 -0.039 1.0 0.150 1. 5 0.394 0 .2 .2 0.072 1. 0 -0.028 0. 5 5 2 2 0.017 1.4 5 0.5 0 .0 .2 5 −0 .50 −0 .25 0.0 0 1. 5 1. 0 0. 5 2.2 11 .6 11 .8 12 .0 12 .2 2.2 0 -0.016 1.0 0.441 5 0.5 0.425 .2 0.353 .2 2 2 0.524 2 2 1 1 50 2.2 .0 .2 11 .6 11 .8 12 .0 12 .2 2.2 0 0.284 0 2 2 .8 .6 150 2.2 1 1 1 1 −1 0. 5 .0 1 1 0.611 0.6 -0.312 .2 .8 .6 1.5 1.0 5 0.5 2.2 11 .6 11 .8 12 .0 12 .2 2.2 0 logc -0.720 5 0.338 2.2 12 1 −0.05 1 logaM/2 2.2 5 11 .6 11 .8 12 .0 12 .2 2.2 0 1.0 2.2 −10 1 10.5 10.0 9.5 1 −50 12 1 0 −1 −9 11 .6 11 .8 12 .0 12 .2 2.2 0 logλ −2 .8 0.05 .6 Ṁ h 0.6 0.00 −0.25 −0.50 11 Ṁ h/Mh 1.4 2.2 −11 0.75 0.50 0.25 11 .6 11 .8 12 .0 12 .2 2.2 0 logM * 0.15 11 logsSFR logSFR 11 .6 11 .8 12 .0 12 .2 logVpeak 2.20 5 11 .6 11 .8 12 .0 12 .2 2.2 0 g-r 50 2.25 0.425 -0.017 -0.696 0.326 0.205 0.917 -0.089 0.735 0.908 -0.830 g-r Figure 3.4: Same as Fig. 3.3, but at fixed log[Vpeak /(km s−1 )] ∼ 2.23. Note particularly the lack of correlation of M∗ with other assembly variables (including c, aM/2 , λ, MÛ h , MÛ h /Mh ), in contrast with the case in Fig. 3.3. 51 In terms of the galaxy-galaxy correlations, the trends are similar to those seen at fixed halo mass, except that the sSFR and color now show clear correlations with M∗ . As a whole, using Vpeak as the halo variable largely removes the correlations between M∗ and other halo assembly properties, and the dependences of SFR, sSFR, and color on halo assembly variables follow each other. 3.4.2.3 Dependence on Mh and Vpeak The correlations shown in Fig. 3.3 and Fig. 3.4 are for haloes of log[Mh /(h−1 M )] ∼ 12.0 and log[Vpeak /(km s−1 )] ∼ 2.23. To obtain a full picture, in Fig. 3.5, we present the Mh and Vpeak -dependent Pearson correlation coefficients for the various pairs of galaxy and halo properties, by performing the calculation in different Mh and Vpeak bins, respectively. The panels correspond to those in Fig. 3.3 and Fig. 3.4, and the correlation shown in a panel of a given row and column is between the property as labelled at the far left of the row and that at the bottom of the column. In each panel, the solid (dashed) curve is the dependence of ρ on Mh (Vpeak ), with zero correlation marked by the black dotted curve. Note that only Vpeak is shown on the x-axis and the corresponding 3 Mh can be obtained according to Mh ∝ Vpeak from equation (3.1). The panels with red curves show the correlations among halo properties. If we limit to halo properties other than Mh and Vpeak , we find that the correlation of any pair of the assembly variables only weakly depends on Mh or Vpeak if any (manifested by the nearly flat curves) and that the correlation strength does not depend on whether we use Mh or Vpeak bins (manifested by the highly overlapped solid and dashed curves). For the galaxy-halo correlations (in panels with black curves), in terms of Mh dependence, the strongest correlation between galaxy and halo property is found between M∗ and Vpeak /c/aM/2 in low mass haloes, with | ρ| ∼ 0.5–0.6. It holds true in the full range of haloes considered here that using Vpeak largely removes the correlation between M∗ and any other halo assembly variable (dashed curves around zero). The only exception is that M∗ appears to be slightly anti-correlated with c in low-Vpeak haloes. With Vpeak , the dependences of correlation on Vpeak for SFR, sSFR, and color closely track each other, which is not the case for those on Mh . With Vpeak as the primary halo variable, star formation-related properties (SFR, sSFR, and color) mainly show dependence on halo formation time and then halo concentration. For galaxy properties (in panels with blue curves), the correlation between SFR and M∗ reaches logλ logaM/2 logc logVpeak 52 1 0 −1 1 Vpeak̇dependence Mḣdependence 0 −1 1 0 −1 1 0 g-r logsSFR logSFR logM * Ṁ h/Mh ρ Ṁ h −1 1 0 −1 1 0 −1 1 0 −1 1 0 −1 1 0 −1 1 0 −1 2.0 2.5 logMh 2.0 2.5 logVpeak 2.0 2.5 logc 2.0 2.5 2.0 2.5 log[Vpeak logaM/2 logλ 2.0 2.5 2.0 /(km s−1)] Ṁ h 2.5 Ṁ h/Mh 2.0 2.5 logM * 2.0 2.5 logSFR 2.0 2.5 logsSFR Figure 3.5: Pearson correlation coefficient ρ of each pair of galaxy and/or halo properties as a function of Mh (solid) and Vpeak (dashed). The galaxy and halo properties are marked to the far left of each row and at the bottom of each column. For clarity, we only label the values of Vpeak on the 3 horizontal axis, and the values of Mh can be inferred from Mh ∝ Vpeak from equation (3.1). As with Fig. 3.3, panels with red, black, and blue curves are for correlations between halo-halo, galaxy-halo, and galaxy-galaxy properties, respectively. In each panel, the dotted horizontal line indicates no correlation. Note particularly the lack of correlation of M∗ with other assembly variables (including c, aM/2 , λ, MÛ h , MÛ h /Mh ) for the Vpeak dependence case, in contrast with the Mh -dependent case. Also the correlations of SFR, sSFR, and color with halo assembly variables show more consistent behaviours in the Vpeak dependence case. 53 a maximum in haloes of ∼ 1012 h−1 M . It weakens in haloes of higher Mh or Vpeak , probably because galaxies move away from the star-forming main sequence and passive evolution starts to dominate. However, the tight sSFR–SFR correlation persists over the full Mh or Vpeak range. For color, the correlation with SFR and sSFR is weak in low-Mh or low-Vpeak haloes and becomes stronger in haloes of higher Mh or Vpeak . The results indicate that in the Illustris simulation, galaxy formation ties to Vpeak more closely than Mh . In comparison with the Mh -based results, we find that in terms of Vpeak , galaxy properties show a cleaner trend in the correlation with other halo assembly variable, such as the lack of correlation for M∗ and the similar correlation pattern for SFR/sSFR/color. It suggests that the Vpeak –based halo model would be a good choice for capturing galaxy assembly bias effect (e.g. with M∗ -based galaxy samples) and for studying galaxy assembly bias (e.g. with SFR/sSFR/color-based samples). With Vpeak as the primary halo variable in the model, halo formation time and concentration would be the main options for the secondary variable to describe the relation between haloes and star formation related quantities, with the former preferable. 3.4.3 Assembly bias of central galaxies With the set of galaxy and halo properties investigated in Section 3.4.2, we do not find a galaxy property that 100 percent correlates with a halo assembly property. For the Mh dependence, the strongest correlation has | ρ| ∼0.5–0.6, between M∗ and Vpeak /c/aM/2 . It means that halo assembly bias cannot be fully inherited by galaxies and be fully translated to galaxy assembly bias. Galaxy assembly bias should be different from halo assembly bias. For example, for haloes of the same mass, we can split them into two halo samples of low and high concentrations and then split central galaxies into two galaxy samples with low and high M∗ . There would be a difference in the clustering of the two halo samples, as well as in that of the two galaxy samples. Given that M∗ is not perfectly correlated with c, we expect that the difference in the galaxy samples is smaller than that in the halo samples. As connecting galaxy assembly bias to halo assembly bias at the halo level can be an important ingredient in incorporating assembly bias effect into clustering model, we develop a simplified model below to understand the connection between galaxy and halo assembly bias. Let us consider haloes at fixed Mh (or Vpeak ) and focus on one halo assembly variable x (e.g. aM/2 or concentration) and one galaxy property y (e.g. M∗ or SFR). Without losing generality, y is assumed to be positively correlated with x. The joint distribution of x and y is illustrated by an 54 ellipse in the top panel of Fig. 3.6. We can form a halo sample by selecting the fraction f of haloes with the highest x (to the right of the vertical dashed line, the region inside the dashed curve) and a galaxy sample by selecting the same fraction of central galaxies with the highest y (above the red-purple dividing line, the region in red). Then, what is the relation between the bias factors of the galaxy and halo sample? To proceed, we make the following assumptions – (1) The joint distribution of galaxy and halo properties follows a 2-dimensional (2D) Gaussian p(x, y), characterised by the centre (xc , yc ), standard deviations σx and σy , and the correlation ρ between x and y. We can take (xc , yc ) = (0, 0) by shifting x and y. A 2D Gaussian function is a reasonable approximation for the distributions seen in Fig. 3.3 and Fig. 3.4, which also follows the Taylor expansion of the logarithmic of the distribution function to the second order. (2) Galaxy property y has a strong dependence on the halo assembly property x, and only weakly on other halo assembly variables. (3) At the fixed Mh (or Vpeak ), halo bias factor is linear with respect to halo property x, b(x) = k x + bc , the first order approximation from Taylor expansion. The value bc is the bias factor at x = 0, which is also the average halo bias factor for haloes at mass Mh (or Vpeak ). The slope k is the first derivative of b with respect to x, k = ∂b/∂ x. For the top f fraction of haloes with the highest x and that of galaxies with the highest y, the distribution of halo property x is shown in the bottom panel of Fig. 3.6 as the dashed curve and red curve, respectively. They are the projections of the region inside the dashed curve and that in red. Clearly, the selected halo and galaxy samples differ in the mean halo property x, lower for the galaxy sample. The mean values can be calculated as ∫ +∞ ∫ +∞ dx −∞ dy xp(x, y) σx exp(−t 2 /2) tσx = ∫ +∞ hxix>tσx = ∫ +∞ ∫ +∞ dx −∞ dy p(x, y) dv exp(−v 2 /2) tσ t (3.3) x and ∫ +∞ tσy dy hxiy>tσy = ∫ +∞ tσy ∫ +∞ dx xp(x, y) −∞ ∫ +∞ dy −∞ ρσx exp(−t 2 /2) = ∫ +∞ , 2 /2) dx p(x, y) dv exp(−v t (3.4) where t is determined by having the correct fraction f , f = ∫ t +∞ 1 dv √ exp(−v 2 /2). 2π (3.5) The bias factors for the halo and galaxy samples are then bh = hbix>tσx = k hxix>tσx + bc (3.6) 55 Figure 3.6: Illustration of the correlation between galaxy property and halo assembly property and the construction of the halo and galaxy samples for the study of galaxy assembly bias effect in Section 3.4.3. In the top panel, the ellipse denotes the joint distribution of halo property x and galaxy property y at fixed Mh or Vpeak , which is assumed to follow a 2D Gaussian distribution. A halo sample is constructed with haloes of the top f fraction of x (indicated by the red region), and a galaxy sample is constructed with central galaxies of the top f fraction of y (indicated by the region inside the dashed curve). Shown in the bottom panel are the probability distribution functions of halo property x for all the haloes at fixed Mh or Vpeak (black solid+dashed), the selected haloes (dashed), and the selected galaxies (red solid), respectively. 56 and bg = hbiy>tσy = k hxiy>tσy + bc . (3.7) We can characterise the assembly bias effect by the fractional difference between the bias factor of the selected halo/galaxy sample and the average halo bias factor at Mh (or Vpeak ), δhb = (bh − bc )/bc and δgb = (bg − bc )/bc . Based on equations (3.3)–(3.7), we have δgb = ρδhb . (3.8) That is, the assembly bias effect of the galaxy sample is weaker than that of the halo sample, by a factor equal to the correlation coefficient of the galaxy and halo property. Only in the case that galaxy and halo properties are tightly correlated (with zero scatter; | ρ| = 1) does halo assembly bias effect completely translate to galaxy assembly bias effect. The connection in equation (3.8) is also valid for samples defined by the property range bounded by two percentiles (i.e. bin samples rather than threshold samples considered here). To test how well the simple model works, we choose a pair of halo and galaxy properties to construct the halo and galaxy samples. We then measure the two-point correlation functions (2PCFs) of the halo and galaxy samples in each Mh and Vpeak bin. To reduce the uncertainty, the large-scale bias factor of a given halo sample is derived from the ratio of the halo-matter two-point cross-correlation function and the matter auto-correlation function (Xu & Zheng (2018a)), bh = ξhm (r)/ξmm (r), averaged over scales of 5-18 h−1 Mpc. The bias factor for the galaxy sample is similarly derived. We consider samples based on halo formation time aM/2 and galaxy color/SFR/sSFR, which show Mh (Vpeak )-dependent correlation coefficient (Fig. 3.5). The results are shown in Fig. 3.7. In the top-left panel, the thick dotted red (blue) curves show the assembly bias quantity δhb for the 50 percent oldest (youngest) haloes as a function of Mh . The solid red (blue) curves are δgb for the 50 percent reddest (bluest) central galaxies. Both quantities decrease with increasing halo mass, i.e. the assembly bias effect becomes weaker for more massive haloes. The thin dashed curves are the same as the dotted curves but modulated by the correlation coefficient between color and aM/2 , i.e. ρδhb , the prediction for δgb from the simple model. The model works well in reproducing the halo mass-dependent galaxy assembly bias effect based on halo assembly bias and the galaxy-halo correlation. The top-middle and top-right panels are for galaxies selected according to SFR and sSFR. The bottom panels show the assembly bias effect as a function of Vpeak . In all the cases, δgb can be well described by ρδhb , which supports the effectiveness of the simple model. 57 0.1 0.0 −0.1 −0.2 −0.3 model model 11.0 11.5 12.0 12.5 log[Mh/(h −1M ⊙ )] 0.2 0.1 0.0 −0.1 −0.2 −0.3 −0.4 11.0 11.5 12.0 12.5 log[Mh/(h −1M ⊙ )] 0.2 0.1 0.0 −0.1 −0.2 −0.3 model model 10.5 11.0 11.5 12.0 12.5 log[Mh/(h −1M ⊙ )] 13.0 13.5 0.0 −0.1 −0.2 model model 1.8 2.0 2.2 2.4 log[Vpeak/(km s−1)] 2.6 50% %oungest haloes 50% oldest haloes 50% highe SFR 50% lowe SFR 0.4 0.3 0.2 0.1 0.0 −0.1 −0.2 model model −0.3 −0.4 13.5 50%⊙ oungest⊙haloes 50%⊙oldest⊙haloes 50%⊙higher⊙sSFR 50%⊙lower⊙sSFR 0.3 b(A; Mh)/b(Mh) − 1 13.0 0.1 0.5 model model 10.5 0.2 −0.4 13.5 50%⊙youngest⊙haloes 50%⊙oldest⊙haloes 50%⊙higher⊙SFR 50%⊙lo er⊙SFR 0.3 b(A; Mh)/b(Mh) − 1 13.0 0.3 −0.3 b(A; Vpeak)/b(Vpeak) − 1 10.5 50% younges haloes 50% oldes haloes 50% blues galaxies 50% reddes galaxies 0.4 1.8 2.0 2.2 2.4 log[Vpeak/(km s−1)] 0.5 2.6 50% %oungest haloes 50% oldest haloes 50% highe sSFR 50% lowe sSFR 0.4 b(A; Vpeak)/b(Vpeak) − 1 b(A; Mh)/b(Mh) − 1 0.2 0.5 b(A; Vpeak)/b(Vpeak) − 1 50% youngest haloes 50% oldest haloes 50% bluest galaxies 50% reddest galaxies 0.3 0.3 0.2 0.1 0.0 −0.1 −0.2 model model −0.3 −0.4 1.8 2.0 2.2 2.4 log[Vpeak/(km s−1)] 2.6 Figure 3.7: Connection between halo and galaxy assembly effect. In the left panels, the quantities shown are the values of b(A, Mh )/b(Mh ) − 1 of different samples. For each sample, the quantity is the fractional difference of the bias factor of the sample selected based on property A with respect to that of all haloes at fixed Mh (i.e. δhb or δgb defined in Section 3.4.3), which is used to characterise the magnitude of the assembly bias effect. The thick dotted curves are for halo samples selected based on halo formation time, and the solid curves are for central galaxy samples selected based on color (top), SFR (middle), and sSFR (bottom). The thin dashed lines are the predictions from the simple model presented in Section 3.4.3 according to the correlation between galaxy and halo properties (δgb = ρδhb , with ρ the correlation coefficient). For clarity, jackknife error bars are only shown for the solid curves in each panel. The right panels are the same, but for the assembly bias effect as a function of Vpeak . See details in Section 3.4.3. 58 The success of the model suggests that an easy recipe could be developed to incorporate galaxy assembly bias into the halo model. The contribution of central galaxies to the galaxy bias factor, in its full form in the simple model, is bg = b c + ρ ∂b xc + σx √ exp(−t 2 /2)/ f , ∂x 2π (3.9) where t = t( f ) is from equation (3.5). As an example, let us use the halo mass Mh as the primary variable in the halo model and consider a galaxy property (color) that correlates with halo formation time (aM/2 or logaM/2 ). In the halo model, besides the average halo bias bc , we also need to know how the halo bias changes with aM/2 (∂b/∂ x), the mean value of aM/2 (xc ), and the scatter in aM/2 (σx ), all as a function of Mh . As usual, we can construct fitting formulae for those four quantities based on N-body simulations. The quantities f and ρ belong to the description of the galaxy-halo relation, which can be parameterised. For f , it is simply the occupation fraction for haloes at Mh . For ρ, the results in Fig. 3.5 suggest that a quadratic form with Mh would suffice. To compute the galaxy bias factor for a galaxy sample, we only need to perform a 1D integral over halo mass. Certainly, it is straightforward to include assembly bias effect by populating dark matter haloes in N-body simulations. However, the above proposal has its virtue for analytic calculations in theoretical investigations. The simple model can also be applied to observation to infer the correlation between galaxy and halo properties. Lin et al. (2016) construct samples of central galaxies from the Sloan Digital Sky Survey data and study the assembly bias effect from the 2PCF measurements. Early and late galaxy samples are defined according to either star formation history (SFH) or sSFR. Weak lensing measurements are used to verify that the host haloes of the early and late galaxies are of similar halo mass (around 1012 h−1 M ). Lin et al. (2016) compare the difference in the early and late galaxy clustering to that in the early and late formed haloes, and do not find evidence for galaxy assembly bias. For SFH-based galaxy samples, they find the ratio of the early to late galaxy bias factor to be 1.00 ± 0.12 (see their fig. 5). If we take the mean of the bias factors of the two galaxy samples as the average halo bias and the uncertainty comes from two similar error bars added in quadrature, the measurement gives δgb = 0 with an uncertainty 0.085. For haloes around 1012 h−1 M , δhb ∼ 0.25, with early- and late-forming haloes (Fig. 3.5). The coefficient of the correlation between SFH and halo formation time is then inferred to be ρ = δgb /δhb = 0.00 ± 0.34. For the sSFR-based samples, the ratio of the early to late galaxy bias factor is 1.07 ± 0.14 (their fig.5). We infer δgb = 0.034 ± 0.099, and with δhb ∼ 0.25 the coefficient of the correlation between sSFR and halo formation time is 59 constrained to be ρ = 0.14 ± 0.40. For both cases, the correlation between galaxy and halo property is consistent with being small. It implies that galaxy SFH and sSFR at most only loosely track halo formation. Similar measurements with large samples can reduce the uncertainty in the inferred correlation, which would help test galaxy formation models (e.g. by comparing to those in Fig. 3.5). 3.5 Summary and discussion Properties in galaxies residing in haloes of the same mass may have a dependence on certain aspects of halo assembly or formation history, which is termed as galaxy assembly bias. Studying galaxy assembly bias and its relation to halo properties can help improve the halo model of galaxy clustering and yield insights in galaxy formation and evolution. Using the Illustris cosmological hydrodynamic galaxy formation simulation, we investigate the central galaxy assembly bias effect through studying the relation among a set of galaxy and halo properties. The main results can be summarised as follows. (1) Central galaxy stellar mass M∗ has a tighter relation with Vpeak than with Mh , manifested by the smaller scatter in M∗ at fixed Vpeak than that at fixed Mh . Once the assembly effect of M∗ on Vpeak is included, M∗ shows nearly no correlation with any other halo assembly properties. (2) The correlations between halo assembly properties and other galaxy properties also appear cleaner if studied at fixed Vpeak , which reveal that galaxy SFR, sSFR, and color mainly correlate with halo formation time (and to a less extent with halo concentration). (3) A simple model is presented to show the relation between galaxy and halo assembly bias, which is linked by the correlation coefficient of the galaxy and halo property in consideration. The Illustris simulation produces a relation between central galaxy stellar mass and halo mass (M∗ –Mh ) similar to that inferred from observation. We find that the scatter in the relation is closely related to halo assembly properties. For example, at fixed Mh , haloes of higher Vpeak or earlier formation tend to host galaxies of higher M∗ . If we choose Vpeak to be the primary halo variable, the scatter in M∗ at fixed Vpeak is reduced compared to that at fixed Mh . Remarkably, once switched to Vpeak , M∗ appears to have nearly no dependence on other halo assembly properties, at least for those considered in our study (including halo concentration, formation time, spin, accretion rate, and specific accretion rate). The property Vpeak , which is an indication of the maximum potential depth over the assembly history of haloes, is able to capture almost all the assembly effect in galaxy stellar 60 mass. The results are in broad agreement with the study using the EAGLE simulation in terms of the z = 0 maximum halo circular velocity (Matthee et al. (2017)). The reason for Vpeak to be the fundamental property in determining M∗ is likely related to the accretion and response of baryons in the gravitational potential. As for the scatter, the simulation noise in Illustris does not contribute much (Genel et al. (2019)). It could be related to chaotic or stochastic processes in star formation and feedback (Matthee et al. (2017); Genel et al. (2019)). Further study is needed to understand the cause of the correlation between M∗ and Vpeak and the origin of the scatter. We present the correlation between each pair of galaxy and halo properties in terms of the Pearson correlation coefficient. Besides M∗ , the other galaxy properties (SFR, sSFR, and color) show a more consistent and clear trend with Vpeak than with Mh . At fixed Vpeak , those other galaxy properties related to star formation are found to mainly correlate with halo formation time and concentration, with stronger correlation with the former. The relatively nice behaviour in the correlations with galaxy properties in the Vpeak -based investigation suggests that it would be advantageous to use Vpeak as the primary variable in the halo model of galaxy clustering, in particular in modelling stellar-mass-based samples. To further model SFR-, sSFR-, or color-selected samples of galaxies, our investigation suggests to introduce halo formation time as the secondary halo variable (and to a less extent, halo concentration). This is in line with the age-matching model of Hearin & Watson (2013), who assumes monotonic mapping between galaxy color and some variant of halo formation time at fixed galaxy luminosity (or stellar mass). That is, there exists a perfect correlation between the galaxy and halo property. Our investigation shows, however, that the correlation coefficient between galaxy SFR/sSFR/color and halo formation time should be included as one important ingredient in the model. The Illustris simulation is able to reproduce the observed star-forming main sequence reasonably well (Sparre et al. (2015)). Here we find that at fixed Mh or Vpeak the relation between SFR and M∗ of central galaxies follows the star-forming main sequence. Interestingly, compared with the SFR–M∗ relation, SFR shows a tighter correlation with sSFR (ρ ∼ 0.9) over the full halo Mh or Vpeak range. It is necessary to test its validity with observations and study its origin by tracking SFR and stellar mass growth of individual galaxies in simulations. The correlation between SFR and sSFR or between any pair of galaxy properties (M∗ , SFR, sSFR, and color) considered here cannot be explained solely by their common dependence on one halo assembly variable. Baryonic processes in galaxy formation (like star formation and feedback) likely play a major role in shaping 61 such correlations. For the effect of assembly bias on galaxy clustering, at fixed Mh or Vpeak , we come up with a simple model to relate the bias factors of a galaxy sample and the corresponding halo sample, which are connected by the correlation coefficient of the galaxy and halo properties used to define the two samples. It gives a reasonable description for the samples constructed with the simulation. It suggests a simple prescription to incorporate galaxy assembly bias into the halo model. By applying the simple model to the galaxy clustering measurements in Lin et al. (2016), we infer that the correlation between SFH/sSFR and halo formation time is consistent with being weak (ρ ∼ 0–0.14). The simple model can be further tested with other hydrodynamic simulations, like EAGLE (Schaye et al. (2015)) and IllustrisTNG (Nelson et al. (2018)), and semi-analytic models, which can also provide further insights on parameterising the correlations between galaxy and halo properties. While our study in this paper focuses on central galaxies, we plan to carry out similar investigations for satellite galaxies to complete the picture of galaxy assembly bias at the halo level to help improve the halo model. CHAPTER 4 CONDITIONAL COLOR MAGNITUDE DISTRIBUTION FROM GALAXY FORMATION MODELS In previous sections, I perform analysis of galaxy-halo relation in a hydrodynamic simulation. In this chapter, I will introduce a particular aspect of galaxy-halo relation study, which is focused on galaxy color. It is well known that a galaxy color magnitude diagram can be decomposed into two main populations, a blue cloud and a red sequence. This result can be considered as an integral effect over all halo mass. The reason that causes red and blue components is galaxy quenching, namely blue star-forming galaxies shut down their star formation and become red passive galaxies. In this chapter, I will focus on red and blue components of galaxies at fixed halo mass, in a hydrodynamic simulation and a semi-analytical galaxy formation model. This investigation can help extract more information encoded in the color magnitude diagram, as well as deepen our understanding of galaxy formation and evolution. 4.1 Introduction There are generally two types of processes that lead to galaxy quenching, environmental effects and feedback effects. Both of them can cause lack of cold gas to suppress star-forming activities. Environmental quenching can be complicated. For satellite galaxies, their cold gas can be removed by tidal stripping and ram pressure when they are accreted onto massive halos (Gunn & Gott III (1972); Read et al. (2006); Van Den Bosch et al. (2008)), and their star formation can be slowly shut down by suppression of gas inflow (strangulation, Peng et al. (2015)). For low mass central galaxies in dense environment, they are surrounding other massive host halos and the gas is more likely to associate with those massive ones. So gas flow into low mass centrals can be shut down and star formation activity slows down when it deplete its own gas reservoir (Von Der Linden et al. (2010); van de Voort et al. (2017); Pintos-Castro et al. (2019)). Massive central galaxies live in massive host halos that have high virial temperature. Gas inside the shock radius is in equilibrium with halo, 63 so it is hard for gas to cool and accrete onto central galaxies (Birnboim & Dekel (2003); Woo et al. (2012)). Besides all of these causes related to environment, stellar processes can also impact star formation rate: for example, the AGN and supernova feedback, which eject energy to heat up cold gas or remove cold gas by wind (Benson (2010); Maiolino et al. (2012)). Since galaxies form and evolve in dark matter halos, it is reasonable to consider that galaxy quenching and color bimodality are related to halos. Halo mass is believed to be one of the main causes of quenching (Zu & Mandelbaum (2016)). To further investigate and model its dependence on other halo properties, it is effective to decompose it along the halo mass direction. Recently, Xu et al. (2018) inferred the color magnitude distribution of galaxies as a function of halo mass (i.e. the conditional color-magnitude distribution or CCMD) from modeling the color and luminosity-dependent clustering of SDSS galaxies with HOD method. They found that the CCMD of central galaxies can be described by a superposition of two almost orthogonal components, a redder component with narrow color distribution but wide luminosity distribution and a bluer component with wide color distribution but narrow luminosity distribution. Those components may have a physical origin related to galaxy formation and evolution. Understanding this physical origin can provide insights on galaxy evolution, and relating it to halo properties can help to improve halo modeling to better explain observed galaxy clustering. In this chapter, we focus on exploring CCMD of central galaxies in a hydrodynamic simulation and a semi-analytical galaxy formation model. We briefly introduce the hydrodynamic simulations and our modeling method in Section 4.2.1, and show results in Section 4.2.2. In Section 4.3.1, we introduce the semi-analytical model and present corresponding results in Section 4.3.2. 4.2 4.2.1 CCMD in hydrodynamic simulations EAGLE simulation and Gaussian mixture model EAGLE (Schaye et al. (2015)) is one of the state-of-the-art hydrodynamic simulations with box size of (100Mpc)3 and 15043 dark matter particles and 15043 baryonic matter particles. The simulation adopts a ΛCDM cosmology with Ωm = 0.307, Ωb = 0.04825, h = 0.6777, σ8 = 0.8288, and ns =0.9611. The particle mass for dark matter and baryonic matter are 9.7 × 106 M and 1.8 × 106 M . Baryonic physical processes included in the simulation are radiative cooling, star formation and feedback, and supermassive black hole formation and feedback. The simulation can reproduce many observations, for example, the evolution of galaxy stellar mass function, the 64 color-magnitude diagram, galaxy size, and HI column density distribution. Halos are identified by FOF, and galaxies are identified by the star particles that are associated with halos. Xu et al. (2018) modeled the red galaxies and blue galaxies on the color-magnitude diagram as two-dimensional Gaussian at fixed halo mass, for centrals and satellites separately. Their modeling can reproduce galaxy color distribution, luminosity function, and clustering. To investigate the Gaussian components more directly from color magnitude distribution, we implementm a 2D Gaussian mixture model on color magnitude distribution of central galaxies in EAGLE. The Gaussian mixture model assumes data distribution is the mixture of several components of 2D Gaussian, and obtains parameters of these Gaussian components by maximizing the log-likelihood function. We use an advanced version of Gaussian mixture model in python, the Bayesian Gaussian mixture model, which maximizes the lower bound of model evidence. In our situation, the general Gaussian mixture model and Bayesian Gaussian mixture produce similar results based on the color-magnitude diagram, but the later can capture Gaussian components better in the whole halo mass range. 4.2.2 CCMD in EAGLE Fig. 4.1 shows the CCMD in the EAGLE simulation. Each panel shows the color-magnitude distribution of central galaxies in bins of log halo mass (at bottom of each panel). Black dots are galaxies and red/blue contours are redder/bluer components found by the Gaussian Mixture model. For the first three bins, the red sequence can be seen, but the algorithm cannot find it. The reason is that the galaxy luminosity function is not complete at this halo mass, and low luminosity galaxies cannot be resolved. It can be seen from the truncation at the low luminosity end of the red sequence. For log halo mass in the range from 11.1 to 12.3, the two components can only be found at lower g-r color, since the red sequence is weak. Only at higher halo mass can the red and blue components be separated. However, the last two bins are wider in mass range because of the small number of dark matter halos, so the result may not be robust for fixed halo mass. With all of the above, we conclude that the evidence of 2D Gaussian red and blue components in the CCMD in the EAGLE simulation is not robust. We also investigate the possibility that the color magnitude distribution has more than two Gaussian components. For this purpose, we use the Akaike information criterion (AIC) and Bayesian information criterion (BIC) to estimate the goodness of model fitting. They are defined by balance between the number of parameters of the model and the maximum likelihood. When increasing the 65 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 (10.5,10.7) (10.7,10.9) (10.9,11.1) (11.1,11.3) (11.3,11.5) (11.5,11.7) (11.7,11.9) (11.9,12.1) (12.1,12.3) (12.3,12.5) (12.5,12.9) (12.9,14.5) 0.0 0.9 0.8 g−r 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 −12.5 −15.0 −17.5 −20.0 −22.5 −25.0 −12.5 −15.0 −17.5 −20.0 −22.5 −25.0 −12.5 −15.0 −17.5 −20.0 −22.5 −25.0 −12.5 −15.0 −17.5 −20.0 −22.5 −25.0 Mr Figure 4.1: Color magnitude distribution of galaxies at fixed halo mass bins from EAGLE. In each panel, black dots are galaxies, and red/blue contours are the two components of the Gaussian mixture model that fits data the best. Halo masses are at the bottom-right, in the form of log[Mh /(h−1 M )]. Except for the last three mass bins, fitted red and blue components cannot be separated. total number of parameters of a model, the maximum likelihood of the model can also increase. So the best fitting model should be the one with relative large likelihood, but not too many parameters. This is done by minimizing AIC and BIC. The main difference between AIC and BIC is that BIC takes the total number of samples into consideration. Fig. 4.2 shows AIC (solid) and BIC (dashed) for modeling the CCMD as Gaussian components, as a function of the number of components. Halo mass for each bin is the same as in Fig. 4.1. It is clear that for most of the halo mass bins (not including the three lowest ones), AIC and BIC both reach minimum at two Gaussian components. Although for one mass bin it prefers one component and the other prefers four components, it is not reasonable to change the number of parameters along halo mass. As mentioned above, the 66 −700 −750 −800 −850 −900 −950 −1000 −3200 AIC, BIC −3250 −3300 −3350 −3400 −400 −420 −440 −460 −480 −500 −520 −540 0 2 4 6 8 10 0 2 4 6 8 10 0 2 4 6 number of components 8 10 0 2 4 6 8 10 Figure 4.2: AIC and BIC of Gaussian mixture models from the EAGLE simulation. In each panel, y-axis is the AIC/BIC estimator, and x-axis is the number of Gaussian components in the model. Solid line represents AIC, and dashed represents BIC. The smaller the estimator value, the better the fitting model. two Gaussian components cannot be separated in many mass bins in the CCMD. So in conclusion, conditional color magnitude distribution in EAGLE cannot be modeled very well with Gaussian components. 4.3 4.3.1 CCMD from a semi-analytical model Galaxy formation model and CCMD components Besides the hydrodynamic simulation, we also explore the CCMD with semi-analytical galaxy formation model. The model we use here is the Guo2011 (Guo et al. (2011)) galaxy formation model, which is built based on the subhalo merger tree of the Millennium simulation (Springel et al. 67 (2005)). The Millennium simulation has the comoving boxsize of (500Mpc/h)3 with 21603 dark matter particles of mass 8.6 × 108 h−1 M . The simulation adopts the cosmology based on WMAP first-year measurement (Spergel et al. (2003)), and the cosmological parameters are Ωm = 0.25, Ωb = 0.045, h = 0.73, σ8 = 0.9, and ns =1. Dark matter halos and subhalos are identified by Friend-of-Friend and SUBFIND (Springel et al. (2001)) algorithm. Physical processes of galaxy formation included are cooling, star forming, angular momentum evolution, black hole growth, and AGN and supernova feedback. The model can recover galaxy luminosity function, metal abundances, Tully-Fisher relation, and galaxy clustering in observation. Fig. 4.3 is the same as Fig. 4.1, but with the galaxy formation model changed to the Guo2011 semi-analytical model. Except for the first four halo mass bins, all the other bins show a red and a blue Gaussian component, and these two components are approximately orthogonal to each other, which is consistent with Xu et al. (2018). It is reasonable to consider the separation of red and blue components as a result of halo assembly history, and it can be linked to different halo assembly variables, for example concentration (ratio of maximum circular velocity and velocity at virial radius of halo, instead of the ratio of halo virial radius to scale radius), Vpeak , half mass scale aM/2 or earlier formation history parameter aM/3 (scale factor of the universe when halo reaches 1/3 of its final mass), scale at which the halo reaches its Vpeak , or parameters related to halo/galaxy major merger. Our final purpose is to find halo assembly variables that are responsible for the appearance of red and blue components. 4.3.2 Relation between CCMD components and halo properties From Section 3.4.2, galaxy color is most related to halo age (half mass scale) and concentration. To explore the relation, we take the galaxies that have high probability (greater than 0.8) to belong to red or blue component found by Gaussian mixture modeling in Section 4.3, and show the joint distribution of half mass scale and concentration of those galaxies in Fig. 4.4. With these kinds of extreme samples, it is easier to see whether the half mass scale and concentration are the cause of red and blue components. The red and blue components have a big difference in the distribution of half mass scale and concentration for the first four bins where the two components cannot be separated very well in Fig. 4.3. On the other hand, in the bins of higher mass where the two orthogonal components can be found, their distribution of half mass scale and concentration are distinguishable. This may imply that concentration and half mass scale are not enough to reveal the origin of red and 68 0.8 0.6 0.4 0.2 (10.9,11.1) (11.1,11.3) (11.3,11.5) (11.5,11.7) (11.7,11.9) (11.9,12.1) (12.1,12.3) (12.3,12.5) (12.5,12.7) (12.7,12.9) (12.9,13.1) (13.1,13.3) (13.3,13.5) (13.5,13.7) (13.7,13.9) 0.0 0.8 0.6 g−r 0.4 0.2 0.0 0.8 0.6 0.4 0.2 0.0 1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.8 0.6 0.4 0.2 (13.9,14.1) (14.1,14.3) (14.3,14.5) 0.0 −15 −20 −25 −15 −20 −25 −15 −20 Mr 0.0 −25 −15 −20 0.0 −25 −15 −20 −25 Figure 4.3: Same as Fig. 4.1, but with Guo2011 semi-analytical galaxy formation model. blue components in the CCMD. We also define red and blue galaxies by color cut, with g − r > 0.75 being red and g − r < 0.45 being blue. However, red and blue samples defined this way still cannot be separated by the distribution of half mass scale and concentration, and this remains the case if concentration is replaced with Vpeak . Galaxy color may also relate to merger history of halos, since galaxy merger is believed to trigger star formation activity. We define a small sample of extreme red and blue galaxies by color at halo mass of 1012 h−1 M , then follow the main branch of subhalo merger tree of each galaxy to find the major mergers. Major merger is the kind of merger in which the mass ratio of the infall satellite subhalo to central subhalo is greater than 1/3. Then we find the time of galaxy merger corresponding to each subhalo major merger. The sample includes 20 red galaxies and 16 blue galaxies; 5 of these blue galaxies don’t have major merger along their formation history. We show the distribution of 69 0.25 0.20 0.15 0.10 0.05 0.00 −0.05 (10.9,11.1) (11.1,11.3) (11.3,11.5) (11.5,11.7) (11.7,11.9) (11.9,12.1) (12.1,12.3) (12.3,12.5) (12.5,12.7) (12.7,12.9) (12.9,13.1) (13.1,13.3) (13.3,13.5) −0.10 0.25 0.20 0.15 0.10 0.05 logc 0.00 −0.05 −0.10 0.25 0.20 0.15 0.10 0.05 0.00 −0.05 −0.10 0.25 0.20 (13.5,13.7) (13.7,13.9) 1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.15 0.10 0.05 0.00 −0.05 −0.10 −1.0 (13.9,14.1) −0.8 −0.6 −0.4 (14.1,14.3) −0.2 −1.0 0.0 −0.8 −0.6 −0.4 (14.3,14.5) −0.2 −1.0 0.0 −0.8 −0.6 −0.4 0.0 −0.2 −1.0 0.0 −0.8 −0.6 −0.4 0.0 −0.2 −1.0 0.0 −0.8 −0.6 −0.4 −0.2 0.0 loga Figure 4.4: Joint distribution of half mass scale and concentration of extreme red and blue components. Extreme color components are the red and blue galaxies found in Section 3.4.2, with a higher than 0.8 probability to belong to one of the color component. major mergers of galaxies in this sample (with at least one major merger) in Fig. 4.5. We examine three quantities of major mergers: the time of the first major merger, the time of the last major merger, and the total number of major mergers along the history. In the top-left panel of Fig. 4.5, we show the time of each major merger of each red/blue galaxy (red/blue dots). Each grid line parallel to the x-axis indicates one galaxy, and dots on grid indicate major merger. Overall, blue galaxies have more mergers at later time than red galaxies. The top-middle/right panel shows the histogram of log scale factor of first/last major merger for galaxies in the top-left panel. From the 1D distribution, it is hard to distinguish red and blue galaxies. Bottom panels are 2D distributions of log scale factor of first major merger, last major merger, and total number of major mergers. The same 6 31 30 29 28 27 26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 7 5 6 4 5 3 4 3 2 2 1 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 loga of majormerger 0 0.0 total majormerger on mainbranch 3.5 3.0 2.5 2.0 1.5 1.0 1 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 loga of first majormerger 4.0 total majormerger on mainbranch 8 0.0 3.5 −0.1 3.0 2.5 2.0 1.5 loga of first majormerger 0.0 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0.0 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0.0 loga of last majormerger 4.0 −0.2 −0.3 −0.4 −0.5 −0.6 1.0 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.0 loga of last majormerger red and blue samples 70 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 loga of last majormerger 0.0 loga of first majormerger Figure 4.5: Distribution of galaxy major merger variables on main branch. Top left: loga distribution of each major merger of each red/blue galaxy along its main branch; the number on y-axis indicates each galaxy. Top middle: histogram of loga of first major merger of red/blue galaxies. Top right: histogram of loga of last major merger of red/blue galaxies. Bottom left: distribution of total number of major merger against loga of first major merger of red/blue galaxies. Bottom middle: distribution of total number of major merger against loga of last major merger. Bottom right: loga of last major merger against loga of first major merger. as 1D distribution, red and blue galaxies cannot be separated into two distinguishable components. So we conclude that halo major merger alone is not a main origin of red and blue components in the CCMD. Halo properties related to assembly history at earlier time are also possible to correlate with galaxy color at later time. Instead of half mass scale aM/2 , we use aM/3 . Instead of Vpeak , we use aVpeak , the scale factor when halo reaches its Vpeak . The distribution of log aM/3 and logaVpeak of red and blue samples are more distinguishable than previous halo properties, with red galaxies tending to form 1/3 of their final mass and reaching Vpeak earlier. However, aM/3 and aVpeak distribution 71 cannot fully separate red and blue galaxies; more information from halo assembly history is needed. In the above preliminary investigation, we show that the information of aM/2 , concentration, and Vpeak cannot reveal the difference between red and blue components in CCMD. Moreover, considering major merger variables alone is not sufficient to separate the two components. The combination of aM/3 and aVpeak is helpful for separating the two components, but more halo assembly variables are necessary. So we intend to explore red and blue components in 3D or higher dimensional halo property space, to see whether a combination of halo properties can be found to explain the origin of these two color components with the help of clustering analysis. CHAPTER 5 SUMMARY In this dissertation, we first present a study of halo assembly bias. From the joint dependence of halo bias on halo mass and one assembly variable, we find that the bias factor increases outward from a global minimum on the joint dependence plane. Our results of one variable dependence are consistent with previous studies measuring bias with percentiles. The possibility of constructing a new halo variable to reduce assembly effect is also explored, and we conclude that it is not possible to absorb halo assembly bias by one single halo variable. Besides one assembly variable dependence, we also study halo bias as functions of two halo variables at fixed halo mass. We show that bias dependence on one variable cannot be inferred correctly from bias dependence on another variable with the correlation between these two variables, if the correlation direction of the two variables is different from bias gradient. We also investigate assembly effects in halo kinematics by measuring halo pairwise velocity. We find that halos with higher bias also have higher pairwise velocity and dispersion. Galaxy assembly bias is related to halo assembly bias through galaxy-halo connection. We study galaxy halo relation and establish a simple model to describe galaxy assembly bias based on Illustris hydrodynamic simulation. We find that central galaxy stellar mass correlates more with halo peak maximum circular velocity rather than halo mass. Moreover, relations between other galaxy properties and halo properties are shown to be more explicit at fixed peak maximum circular velocity. Then we model galaxy assembly bias based on galaxy-halo relation. With a few assumptions, we predict that the fractional galaxy assembly bias is smaller than the fractional halo assembly bias by a factor of correlation coefficient between galaxy property and halo property. This prediction is consistent with galaxy assembly bias trends in simulation. Our model can also be used to model galaxy bias and to constrain galaxy-halo relation in observation. For understanding an important part of galaxy-halo relation, we focus on the conditional color magnitude distribution (CCMD) in both hydrodynamic simulation and the semi-analytical galaxy formation model. Orthogonal red and blue components in the CCMD of central galaxies have been 73 found in the latter but not the former. Our purpose is to find the physical origin of these red and blue components from properties of host dark matter halos. We find that the parameters like halo formation time, concentration, peak maximum circular velocity, merger history parameters, and their combinations cannot lead to a good separation of the features in the CCMD. The combination of earlier halo formation history and time of reaching peak maximum circular velocity is shown to be helpful for separating red and blue components. In the future, investigations of higher dimensional halo property space are necessary to explain the origin of red and blue components in the CCMD. APPENDIX MASS RESOLUTION EFFECT ON HALO BIAS AND PAIRWISE VELOCITY STATISTICS In the main analysis in this paper, we primarily use haloes from the MDR1 simulation, which has particle mass of 8.721 × 109 h−1 M . Paranjape & Padmanabhan (2017) show that for haloes with less than 400 particles, halo concentration may not be determined accurately, usually lower than the true value (their Fig. 2; also see Fig. 9 of Trenti et al. (2010)). The distribution of other assembly variables (such as spin) for haloes of low number of particles can also be affected (e.g. Fig. 9 in Trenti et al. (2010) and Fig. 3 in Benson (2017)). Such a numerical effect would potentially affect the pattern of the dependence of halo bias and velocity statistics on assembly variables for halos of mass below 3.5 × 1012 h−1 M in the MDR1 simulation. To test the effect, we make use of the Bolshoi simulation, with the same cosmology but a 64 times higher mass resolution. The mass of haloes with 400 particles in the Bolshoi simulation is 5.5 × 1010 h−1 M , well below the minimum mass in our analysis. First, we perform similar calculations as in Fig. 2.1, but using haloes in the Bolshoi simulation for log[Mh /(h−1 M )] < 12.8. The results are presented in Fig. A.1. By comparing with those in Fig. 2.1, we see that the trend in the joint dependence of halo bias on halo mass and one assembly variable from the Bolshoi simulation is similar to that from the MDR1 simulation in the mass range considered in the paper. The test assures that we can present the main results based on MDR1 simulation. The large volume of MDR1 compared to Bolshoi has the advantage of making the halo bias calculation less noisy and enabling the inclusion of a larger range in assembly variable at fixed mass, which is clearly seen when comparing Fig. 2.1 and Fig. A.1. For the above reason, the results on bias (except for those in Section 2.4.2; see below) from the MDR1 simulation are shown in the paper. We caution, however, that for investigations or applications that require accurate values of halo bias, the mass resolution effect needs to be accounted for. For the investigation with effective halo mass in Section 2.4.2, we need to combine the Bolshoi and MDR1 simulation. This is not mainly driven by the above numerical effect, but the need for 75 3.0 −0.1 3.0 2 2.2 1.5 2.0 1.6 MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 −0.7 −0.8 Bolshoi-400mp 3.5 3.0 MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 2.5 b 0.8 1.5 0.6 1.0 0.4 Bolshoi-400mp MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 0.5 0.0 −1.5 λ 2.0 log c 2.0 0.0 3.0 −1.0 2.5 1.0 0.5 3.5 1.2 log 1.0 −0.6 0.0 1.6 0.2 1.5 −0.5 0.5 1.8 1.4 2.0 −0.4 1.0 Bolshoi-400mp 2.5 −0.3 aM/ 2.0 b 2.4 1.8 −0.2 b peak 3.5 2.5 2.6 log[ 0.0 b 2.8 3.5 log 1 V /(km s− )] 3.0 1.5 −2.0 1.0 −2.5 Bolshoi-400mp MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 0.5 0.0 Figure A.1: Same as Fig. 2.1, but the Bolshoi simulation is used for log[Mh /(h−1 M )] < 12.8 and MDR1 simulation for higher mass. The two vertical dashed lines in the each panel indicate the masses of haloes with 400 particles in the two simulations. low mass haloes. As can be seen from the left panel of Fig. 2.2, for effective halo mass in the range log[Meff /(h−1 M )] ∼ 11–12.5, we need to reach haloes with high spin and low mass (below log[Mh /(h−1 M )] ∼11.6), which are covered incompletely in the MDR1 simulation. A combination of MDR1 and Bolshoi is therefore necessary for the effective mass analysis (see Section 2.4.2 for details). Then we analyse the pairwise velocity statistics of haloes of log[Mh /(h−1 M )] ∼11.73 in the Bolshoi simulation. Fig. A.2 shows the dependence of pairwise velocity and velocity dispersion on assembly variable, which is to be compared to Fig. 2.10 with the MDR1 simulation. For the cases with Vpeak , aM/2 , and c (the left three columns), the trends from the Bolshoi simulation well track those from the MDR1 simulation. There are slight amplitude shifts for the corresponding curves from the two simulations, which can be partly explained by the sample variance effect (i.e. difference in the long-wave length modes in the two simulations). The mass resolution may also contribute. Imagine that we have a simulation with the same box size and initial condition as Bolshoi but with MDR1 mass resolution. For Bolshoi haloes of a given concentration, those of the corresponding 76 300 (kms−1 ) 250 200 1.0 logM=11.73 c dependence ( Mpc) r h −1 400 1.0 logM=11.73 λ dependence ( Mpc) r h −1 10.0 2(kms−1 ) , t/ 12 σ σ12,t / 2(kms−1 ) 2(kms−1 ) 1.0 logM=11.73 c dependence 1.0 0 0.1 500 logM=11.73 λ dependence 400 1.0 − r h lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 ( 1 Mpc) logM=11.73 λ dependence 300 200 100 0 0.1 100 r(h −1 Mpc) 200 50 r(h −1 Mpc) 200 lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 300 100 v12, r 400 lower-middle 25% upper-middle 25% upper 25% lower 25% 10.0 300 0 0.1 500 lower 25% lower-middle 25% upper-middle 25% upper 25% 0 0.1 500 logM=11.73 c dependence 100 10.0 150 0 0.1 −1 , t/ v12, r 50 100 r(h Mpc) 1.0 200 100 0 0.1 200 lower-middle 25% upper-middle 25% upper 25% lower 25% 10.0 300 150 logM=11.73 aM/2 dependence 300 0 0.1 400 lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 ( 1 Mpc) 12 200 400 1.0 − r h σ (kms−1 ) 250 500 logM=11.73 aM/2 dependence 100 500 lower 25% lower-middle 25% upper-middle 25% upper 25% 0 0.1 ( Mpc) r h −1 200 logM=11.73 aM/2 dependence 1.0 10.0 r(h −1 Mpc) 100 2(kms−1 ) 300 1.0 300 100 0 0.1 0 0.1 400 150 50 100 500 lower-middle 25% upper-middle 25% upper 25% lower 25% 200 lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 1.0 lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 r(h −1 Mpc) , t/ 200 300 100 12 v12,r (kms−1 ) 250 logM=11.73 Vpeak dependence 1.0 10.0 r(h −1 Mpc) logM=11.73 peak dependence V σ 300 σ12, r (kms −1 ) v12, r 0 0.1 400 200 100 50 500 logM=11.73 Vpeak dependence 300 150 σ12,r (kms−1 ) 200 400 σ12, r (kms −1 ) (kms−1 ) 250 500 lower 25% lower-middle 25% upper-middle 25% upper 25% σ12, r (kms −1 ) 300 0 0.1 1.0 − r h lower 25% lower-middle 25% upper-middle 25% upper 25% 10.0 ( 1 Mpc) Figure A.2: Same as Fig. 2.10, but haloes from the Bolshoi simulation are used for the analysis. 77 MDR1-resolution haloes would have a scatter (e.g. Paranjape & Padmanabhan (2017); Trenti et al. (2010)). That is, for the same percentile, haloes from the two simulations would not have the exact correspondence, and the pairwise velocity amplitude is expected to be slightly different. For the above three assembly variables, the pattern of the dependence of pairwise velocity statistics on the value of assembly bias from the Bolshoi simulation is similar to that from the MDR1 simulation, as can be seen by comparing Fig. A.3 and Fig. 2.11. For the case with halo spin λ, the pairwise velocity statistics from the percentile analysis seem to be different between the MDR1 and Bolshoi simulations (right panels of Fig. 2.10 and Fig. A.2). With the MDR1 simulation, there is a clear trend of higher pairwise velocity and velocity dispersion for haloes of higher spin, although it is substantially weaker than those with the other three assembly variables. With the Bolshoi simulation, however, we do not see such a clear trend (also see the corresponding panels of Fig. A.3). If what Bolshoi suggests is close to the truth, given the expected scatter in determining the spin parameter for a low-resolution simulation (e.g. Trenti et al. (2010); Benson (2017)), we would not expect to see any trend with the MDR1 simulation. However, here we have the opposite results. We note that the error bars with the Bolshoi simulation analyses are substantially larger than those with the MDR1. It is possible that the sample variance effect in the Bolshoi simulation masks the spin dependence, i.e. the noise level is too high to reveal the intrinsically weak trend. The other possibility is that with low-resolution simulations, the spin measurement is more easily affected by environment and hence halo assembly, making the scatter mentioned above not random. Further investigations using simulations of the same initial condition and different resolutions, like those in Trenti et al. (2010), would help resolve the issue of the spin-dependence of the pairwise velocity statistics. To summarise, from comparing the halo statistics with MDR1 and Bolshoi simulations, for most results presented in this paper, the MDR1 simulation works well in revealing the patterns related to assembly variables, which has the benefit of better statistics from the much larger volume. For the effective halo mass analysis, we have to include the Bolshoi simulation to reach haloes of much lower mass for completeness consideration. For halo pairwise velocity statistics, the spin case appears to be different between the two simulations, weak trend with the MDR1 simulation and the lack of trend with the Bolshoi simulation. The cause of the apparent discrepancy is not entirely clear, which could be noise in the Bolshoi simulation or some unknown environment-dependent systematic effect in halo spin measurement. In the main text, we choose to present the spin results from the MDR1 78 Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 1.6 MDR1-400mp Mh / h M ⊙ )] log[ ( −1 2 (km s −1 ) 1 peak 1 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 50 0 200 1.8 1.6 150 100 Bolshoi-400mp MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 50 0 0.0 500 0.0 500 450 −0.1 450 −0.1 450 −0.8 Mh / h M ⊙ )] ( −1 ,r (km s −1 ) 150 50 −0.7 0 −0.8 MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 2 250 −0.5 200 150 −0.6 100 Bolshoi-400mp 300 −0.4 aM/ log −0.5 200 350 −0.3 50 −0.7 0 −0.8 σ ,t / MDR1-400mp 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ 100 Bolshoi-400mp MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 50 0 500 1.8 500 1.8 500 1.6 450 1.6 450 1.6 450 Bolshoi-400mp MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 500 Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 0 Mh / h M ⊙ )] log[ ( −1 λ 200 −2.0 400 350 300 −1.5 250 200 −2.0 150 100 −2.5 Bolshoi-400mp MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 50 0 0 450 −1.0 σ 150 50 500 250 100 50 MDR1-400mp 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 300 −1.5 v 12 200 MDR1-400mp 100 Bolshoi-400mp 350 log λ log 250 −2.5 Bolshoi-400mp 0.2 0 ,r (km s −1 ) ,r (km s −1 ) 300 150 50 400 350 −2.0 150 0.4 450 −1.0 400 −1.5 200 500 450 −1.0 250 0.6 σ ,t / 0.2 0 c 100 0.4 50 log 150 2 (km s −1 ) log[ ( −1 0.6 300 0.8 σ ,t / Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 200 1.0 λ MDR1-400mp 250 0.8 350 1.2 log Bolshoi-400mp ,r (km s −1 ) c 100 0.4 300 1.0 12 150 log 0.6 v 200 12 250 0.8 350 1.2 400 1.4 σ ,r (km s −1 ) c 300 1.0 0.2 1.4 350 1.2 400 12 400 2 (km s −1 ) 1.8 1.4 log 2 100 Bolshoi-400mp 250 −0.6 v −0.6 −0.4 12 150 300 aM/ −0.5 200 log 250 350 −0.3 400 −0.2 σ ,r (km s −1 ) −0.4 12 300 −0.7 −0.2 350 −0.3 400 2 (km s −1 ) 500 400 12 Bolshoi-400mp V /(km s− )] ,r (km s −1 ) 100 250 2.0 σ ,t / 1.8 150 0.0 aM/ 2 0 2.0 300 2.2 −0.1 −0.2 log 50 200 350 2.4 12 MDR1-400mp 250 2.2 400 2.6 12 Bolshoi-400mp 300 450 2.8 100 −2.5 Bolshoi-400mp MDR1-400mp Mh / h M ⊙ )] 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 log[ ( −1 12 1.6 100 350 log[ 1.8 150 2.6 σ 2.0 400 2.4 peak 200 v peak log[ 250 2.2 2.8 500 3.0 450 V /(km s− )] 300 2.4 log[ 350 ,r (km s −1 ) 2.6 12 1 V /(km s− )] 400 500 3.0 450 2.8 12 500 3.0 50 0 Figure A.3: Same as Fig. 2.11, but the Bolshoi simulation is used for log[Mh /(h−1 M )] < 12.8 and MDR1 simulation for higher mass. The two vertical dashed lines in the each panel indicate the masses of haloes with 400 particles in the two simulations. 79 simulation, but put the caveat that there is an apparent difference compared to those from the Bolshoi simulation. REFERENCES Bardeen, J. M., Szalay, A., Kaiser, N., & Bond, J. 1986, The Astrophysical Journal, 304, 15 Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, The Astrophysical Journal, 717, 379 Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, The Astrophysical Journal, 762, 109 Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., Busha, M. T., Klypin, A. A., & Primack, J. R. 2012, The Astrophysical Journal, 763, 18 Benson, A. 2010, Physics Reports, 495, 33 Benson, A. J. 2017, Monthly Notices of the Royal Astronomical Society, 471, 2871 Berlind, A. A., Kazin, E., Blanton, M. R., Pueblas, S., Scoccimarro, R., & Hogg, D. W. 2006, arXiv preprint astro-ph/0610524 Berlind, A. A., & Weinberg, D. H. 2002, The Astrophysical Journal, 575, 587 Berlind, A. A., Weinberg, D. H., Benson, A. J., Baugh, C. M., Cole, S., Davé, R., Frenk, C. S., Jenkins, A., Katz, N., & Lacey, C. G. 2003, The Astrophysical Journal, 593, 1 Birnboim, Y., & Dekel, A. 2003, Monthly Notices of the Royal Astronomical Society, 345, 349 Blake, C., Kazin, E. A., Beutler, F., Davis, T. M., Parkinson, D., Brough, S., Colless, M., Contreras, C., Couch, W., Croom, S., et al. 2011, Monthly Notices of the Royal Astronomical Society, 418, 1707 Bond, J., Cole, S., Efstathiou, G., & Kaiser, N. 1991, The Astrophysical Journal, 379, 440 Borzyszkowski, M., Porciani, C., Romano-Díaz, E., & Garaldi, E. 2017, Monthly Notices of the Royal Astronomical Society, 469, 594 Brinchmann, J., Charlot, S., White, S., Tremonti, C., Kauffmann, G., Heckman, T., & Brinkmann, J. 2004, Monthly Notices of the Royal Astronomical Society, 351, 1151 Busch, P., & White, S. 2017, Monthly Notices of the Royal Astronomical Society, 470, 4767 Castorina, E., & Sheth, R. K. 2013, Monthly Notices of the Royal Astronomical Society, 433, 1529 Chaves-Montero, J., Angulo, R. E., Schaye, J., Schaller, M., Crain, R. A., Furlong, M., & Theuns, T. 2016, Monthly Notices of the Royal Astronomical Society, 460, 3100 Collessy, M. 1999, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 357, 105 Contreras, S., Zehavi, I., Padilla, N., Baugh, C., Jiménez, E., & Lacerna, I. 2019, Monthly Notices of the Royal Astronomical Society, 484, 1133 81 Cooray, A., & Sheth, R. 2002, Physics Reports, 372, 1 Croft, R. A., Matteo, T. D., Khandai, N., Springel, V., Jana, A., & Gardner, J. P. 2012, Monthly Notices of the Royal Astronomical Society, 425, 2766 Croton, D. J., Gao, L., & White, S. D. 2007, Monthly Notices of the Royal Astronomical Society, 374, 1303 Dalal, N., White, M., Bond, J. R., & Shirokov, A. 2008, The Astrophysical Journal, 687, 12 Davis, M., Efstathiou, G., Frenk, C., & White, S. 1985, The Astrophysical Journal, 292, 371 Dawson, K. S., Kneib, J.-P., Percival, W. J., Alam, S., Albareti, F. D., Anderson, S. F., Armengaud, E., Aubourg, É., Bailey, S., Bautista, J. E., et al. 2016, The Astronomical Journal, 151, 44 De Lucia, G., & Blaizot, J. 2007, Monthly Notices of the Royal Astronomical Society, 375, 2 Eisenstein, D. J., Weinberg, D. H., Agol, E., Aihara, H., Prieto, C. A., Anderson, S. F., Arns, J. A., Aubourg, É., Bailey, S., Balbinot, E., et al. 2011, The Astronomical Journal, 142, 72 Faltenbacher, A., & White, S. D. 2009, The Astrophysical Journal, 708, 469 Gao, L., Springel, V., & White, S. D. 2005, Monthly Notices of the Royal Astronomical Society: Letters, 363, L66 Gao, L., & White, S. D. 2007, Monthly Notices of the Royal Astronomical Society: Letters, 377, L5 Gao, L., White, S. D., Jenkins, A., Stoehr, F., & Springel, V. 2004, Monthly Notices of the Royal Astronomical Society, 355, 819 Garaldi, E., Romano-Díaz, E., Borzyszkowski, M., & Porciani, C. 2017, Monthly Notices of the Royal Astronomical Society, 473, 2234 Genel, S., Bryan, G. L., Springel, V., Hernquist, L., Nelson, D., Pillepich, A., Weinberger, R., Pakmor, R., Marinacci, F., & Vogelsberger, M. 2019, The Astrophysical Journal, 871, 21 Genel, S., Vogelsberger, M., Springel, V., Sijacki, D., Nelson, D., Snyder, G., Rodriguez-Gomez, V., Torrey, P., & Hernquist, L. 2014, Monthly Notices of the Royal Astronomical Society, 445, 175 Gu, M., Conroy, C., & Behroozi, P. 2016, The Astrophysical Journal, 833, 2 Gunn, J. E., & Gott III, J. R. 1972, The Astrophysical Journal, 176, 1 Guo, H., Li, C., Zheng, Z., Mo, H., Jing, Y., Zu, Y., Lim, S., & Xu, H. 2017, The Astrophysical Journal, 846, 61 Guo, H., Zheng, Z., Behroozi, P. S., Zehavi, I., Chuang, C.-H., Comparat, J., Favole, G., Gottloeber, S., Klypin, A., Prada, F., et al. 2016, Monthly Notices of the Royal Astronomical Society, 459, 3040 Guo, Q., White, S., Boylan-Kolchin, M., De Lucia, G., Kauffmann, G., Lemson, G., Li, C., Springel, V., & Weinmann, S. 2011, Monthly Notices of the Royal Astronomical Society, 413, 101 Hahn, O., Carollo, C. M., Porciani, C., & Dekel, A. 2007, Monthly Notices of the Royal Astronomical Society, 381, 41 Hahn, O., & Paranjape, A. 2013, Monthly Notices of the Royal Astronomical Society, 438, 878 82 Hahn, O., Porciani, C., Dekel, A., & Carollo, C. M. 2009, Monthly Notices of the Royal Astronomical Society, 398, 1742 Han, J., Li, Y., Jing, Y., Nishimichi, T., Wang, W., & Jiang, C. 2019, Monthly Notices of the Royal Astronomical Society, 482, 1900 Harker, G., Cole, S., Helly, J., Frenk, C., & Jenkins, A. 2006, Monthly Notices of the Royal Astronomical Society, 367, 1039 Hearin, A., & Watson, D. 2013, Monthly Notices of the Royal Astronomical Society, 435, 1313 Hearin, A. P. 2015, Monthly Notices of the Royal Astronomical Society: Letters, 451, L45 Hearin, A. P., Zentner, A. R., van den Bosch, F. C., Campbell, D., & Tollerud, E. 2016, Monthly Notices of the Royal Astronomical Society, 460, 2552 Jeeson-Daniel, A., Vecchia, C. D., Haas, M. R., & Schaye, J. 2011, Monthly Notices of the Royal Astronomical Society: Letters, 415, L69 Jing, Y. 1999, The Astrophysical Journal Letters, 515, L45 Jing, Y., Suto, Y., & Mo, H. 2007, The Astrophysical Journal, 657, 664 Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, Monthly Notices of the Royal Astronomical Society, 363, 2 Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, The Astrophysical Journal, 740, 102 Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., Klypin, A. A., Gottlöber, S., Allgood, B., & Primack, J. R. 2004, The Astrophysical Journal, 609, 35 Lazeyras, T., Musso, M., & Schmidt, F. 2017, Journal of Cosmology and Astroparticle Physics, 2017, 059 Leauthaud, A., Tinker, J., Bundy, K., Behroozi, P. S., Massey, R., Rhodes, J., George, M. R., Kneib, J.-P., Benson, A., Wechsler, R. H., et al. 2012, The Astrophysical Journal, 744, 159 Lin, Y.-T., Mandelbaum, R., Huang, Y.-H., Huang, H.-J., Dalal, N., Diemer, B., Jian, H.-Y., & Kravtsov, A. 2016, The Astrophysical Journal, 819, 119 Maccio, A. V., Dutton, A. A., Van Den Bosch, F. C., Moore, B., Potter, D., & Stadel, J. 2007, Monthly Notices of the Royal Astronomical Society, 378, 55 Maggiore, M., & Riotto, A. 2010, The Astrophysical Journal, 711, 907 Maiolino, R., Gallerani, S., Neri, R., Cicone, C., Ferrara, A., Genzel, R., Lutz, D., Sturm, E., Tacconi, L., Walter, F., et al. 2012, Monthly Notices of the Royal Astronomical Society: Letters, 425, L66 Mansfield, P., & Kravtsov, A. 2019, arXiv preprint arXiv:1903.02007 Mao, Y.-Y., Zentner, A. R., & Wechsler, R. H. 2018, Monthly Notices of the Royal Astronomical Society, 474, 5143 Matthee, J., Schaye, J., Crain, R. A., Schaller, M., Bower, R., & Theuns, T. 2017, Monthly Notices of the Royal Astronomical Society, stw2884 83 McAlpine, S., Helly, J. C., Schaller, M., Trayford, J. W., Qu, Y., Furlong, M., Bower, R. G., Crain, R. A., Schaye, J., Theuns, T., et al. 2016, Astronomy and Computing, 15, 72 McCarthy, K. S., Zheng, Z., & Guo, H. 2018, arXiv preprint arXiv:1810.05183 McEwen, J. E., & Weinberg, D. H. 2016, arXiv preprint arXiv:1601.02693 Mehta, K. 2014, PhD thesis, The University of Arizona Miyatake, H., More, S., Takada, M., Spergel, D. N., Mandelbaum, R., Rykoff, E. S., & Rozo, E. 2016, Physical Review Letters, 116, 041301 Mo, H., & White, S. D. 1996, Monthly Notices of the Royal Astronomical Society, 282, 347 Navarro, J. F., Frenk, C. S., & White, S. D. 1997, The Astrophysical Journal, 490, 493 Nelson, D., Pillepich, A., Genel, S., Vogelsberger, M., Springel, V., Torrey, P., Rodriguez-Gomez, V., Sijacki, D., Snyder, G. F., Griffen, B., et al. 2015, Astronomy and Computing, 13, 12 Nelson, D., Springel, V., Pillepich, A., Rodriguez-Gomez, V., Torrey, P., Genel, S., Vogelsberger, M., Pakmor, R., Marinacci, F., Weinberger, R., et al. 2018, arXiv preprint arXiv:1812.05609 Paranjape, A., Hahn, O., & Sheth, R. K. 2018, Monthly Notices of the Royal Astronomical Society, 476, 3631 Paranjape, A., Lam, T. Y., & Sheth, R. K. 2012, Monthly Notices of the Royal Astronomical Society, 420, 1429 Paranjape, A., & Padmanabhan, N. 2017, Monthly Notices of the Royal Astronomical Society, 468, 2984 Paranjape, A., Sheth, R. K., & Desjacques, V. 2013, Monthly Notices of the Royal Astronomical Society, 431, 1503 Peng, Y., Maiolino, R., & Cochrane, R. 2015, Nature, 521, 192 Pintos-Castro, I., Yee, H., Muzzin, A., Old, L., & Wilson, G. 2019, arXiv preprint arXiv:1904.00023 Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, Monthly Notices of the Royal Astronomical Society, 423, 3018 Press, W. H., & Schechter, P. 1974, The Astrophysical Journal, 187, 425 Ramakrishnan, S., Paranjape, A., Hahn, O., & Sheth, R. K. 2019, arXiv preprint arXiv:1903.02007 Read, J. I., Wilkinson, M., Evans, N., Gilmore, G., & Kleyna, J. T. 2006, Monthly Notices of the Royal Astronomical Society, 366, 429 Reddick, R., Wechsler, R., Tinker, J., & Behroozi, P. 2013, The Astrophysical Journal, 771, 30 Riebe, K., Partl, A. M., Enke, H., Forero-Romero, J., Gottlöber, S., Klypin, A., Lemson, G., Prada, F., Primack, J. R., Steinmetz, M., et al. 2013, Astronomische Nachrichten, 334, 691 Salcedo, A. N., Maller, A. H., Berlind, A. A., Sinha, M., McBride, C. K., Behroozi, P. S., Wechsler, R. H., & Weinberg, D. H. 2018, Monthly Notices of the Royal Astronomical Society, 475, 4411 84 Santini, P., Fontana, A., Castellano, M., Di Criscienzo, M., Merlin, E., Amorin, R., Cullen, F., Daddi, E., Dickinson, M., Dunlop, J. S., et al. 2017, The Astrophysical Journal, 847, 76 Schaye, J., Crain, R. A., Bower, R. G., Furlong, M., Schaller, M., Theuns, T., Dalla Vecchia, C., Frenk, C. S., McCarthy, I., Helly, J. C., et al. 2015, Monthly Notices of the Royal Astronomical Society, 446, 521 Sheth, R. K., & Diaferio, A. 2001, Monthly Notices of the Royal Astronomical Society, 322, 901 Sheth, R. K., Diaferio, A., Hui, L., & Scoccimarro, R. 2001, Monthly Notices of the Royal Astronomical Society, 326, 463 Sheth, R. K., & Tormen, G. 1999, Monthly Notices of the Royal Astronomical Society, 308, 119 —. 2002, Monthly Notices of the Royal Astronomical Society, 329, 61 Shi, J., & Sheth, R. K. 2017, Monthly Notices of the Royal Astronomical Society, 473, 2486 Shi, J., Wang, H., & Mo, H. 2015, The Astrophysical Journal, 807, 37 Sparre, M., Hayward, C. C., Springel, V., Vogelsberger, M., Genel, S., Torrey, P., Nelson, D., Sijacki, D., & Hernquist, L. 2015, Monthly Notices of the Royal Astronomical Society, 447, 3548 Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, The Astrophysical Journal Supplement Series, 214, 15 Spergel, D. N., Verde, L., Peiris, H. V., Komatsu, E., Nolta, M., Bennett, C., Halpern, M., Hinshaw, G., Jarosik, N., Kogut, A., et al. 2003, The Astrophysical Journal Supplement Series, 148, 175 Springel, V., White, S. D., Jenkins, A., Frenk, C. S., Yoshida, N., Gao, L., Navarro, J., Thacker, R., Croton, D., Helly, J., et al. 2005, Nature, 435, 629 Springel, V., White, S. D., Tormen, G., & Kauffmann, G. 2001, Monthly Notices of the Royal Astronomical Society, 328, 726 Sunayama, T., Hearin, A. P., Padmanabhan, N., & Leauthaud, A. 2016, Monthly Notices of the Royal Astronomical Society, 458, 1510 Tinker, J., Kravtsov, A. V., Klypin, A., Abazajian, K., Warren, M., Yepes, G., Gottlöber, S., & Holz, D. E. 2008, The Astrophysical Journal, 688, 709 Tinker, J. L. 2017, Monthly Notices of the Royal Astronomical Society, 467, 3533 Tinker, J. L., Brownstein, J. R., Guo, H., Leauthaud, A., Maraston, C., Masters, K., Montero-Dorta, A. D., Thomas, D., Tojeiro, R., Weiner, B., et al. 2017, The Astrophysical Journal, 839, 121 Trenti, M., Smith, B. D., Hallman, E. J., Skillman, S. W., & Shull, J. M. 2010, The Astrophysical Journal, 711, 1198 van de Voort, F., Bahé, Y., Bower, R., Correa, C., Crain, R., Schaye, J., & Theuns, T. 2017, Monthly Notices of the Royal Astronomical Society, 466, 3460 Van Den Bosch, F. C., Aquino, D., Yang, X., Mo, H., Pasquali, A., McIntosh, D. H., Weinmann, S. M., & Kang, X. 2008, Monthly Notices of the Royal Astronomical Society, 387, 79 85 Vogelsberger, M., Genel, S., Springel, V., Torrey, P., Sijacki, D., Xu, D., Snyder, G., Bird, S., Nelson, D., & Hernquist, L. 2014a, Nature, 509, 177 Vogelsberger, M., Genel, S., Springel, V., Torrey, P., Sijacki, D., Xu, D., Snyder, G., Nelson, D., & Hernquist, L. 2014b, Monthly Notices of the Royal Astronomical Society, 444, 1518 Von Der Linden, A., Wild, V., Kauffmann, G., White, S. D., & Weinmann, S. 2010, Monthly Notices of the Royal Astronomical Society, 404, 1231 Wang, H., Mo, H., & Jing, Y. 2009, Monthly Notices of the Royal Astronomical Society, 396, 2249 Wang, H., Mo, H., Jing, Y., Yang, X., & Wang, Y. 2011, Monthly Notices of the Royal Astronomical Society, 413, 1973 Watson, D. F., Hearin, A. P., Berlind, A. A., Becker, M. R., Behroozi, P. S., Skibba, R. A., Reyes, R., Zentner, A. R., & van den Bosch, F. C. 2015, Monthly Notices of the Royal Astronomical Society, 446, 651 Wechsler, R., & Tinker, J. 2018, Annual Review of Astronomy and Astrophysics, 56, 435 Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, The Astrophysical Journal, 568, 52 Wechsler, R. H., Zentner, A. R., Bullock, J. S., Kravtsov, A. V., & Allgood, B. 2006, The Astrophysical Journal, 652, 71 Wetzel, A. R., Tinker, J. L., Conroy, C., & Bosch, F. C. v. d. 2014, Monthly Notices of the Royal Astronomical Society, 439, 2687 White, S. D., & Rees, M. J. 1978, Monthly Notices of the Royal Astronomical Society, 183, 341 Woo, J., Dekel, A., Faber, S., Noeske, K., Koo, D. C., Gerke, B. F., Cooper, M. C., Salim, S., Dutton, A. A., Newman, J., et al. 2012, Monthly Notices of the Royal Astronomical Society, 428, 3306 Xu, H., Zheng, Z., Guo, H., Zu, Y., Zehavi, I., & Weinberg, D. H. 2018, Monthly Notices of the Royal Astronomical Society, 481, 5470 Xu, X., & Zheng, Z. 2018a, Monthly Notices of the Royal Astronomical Society, 479, 1579 —. 2018b, arXiv preprint arXiv:1812.11210 Yang, X., Mo, H., & Bosch, F. C. v. d. 2003, Monthly Notices of the Royal Astronomical Society, 339, 1057 Yang, X., Mo, H., & van den Bosch, F. C. 2006, The Astrophysical Journal Letters, 638, L55 York, D. G., Adelman, J., Anderson Jr, J. E., Anderson, S. F., Annis, J., Bahcall, N. A., Bakken, J., Barkhouser, R., Bastian, S., Berman, E., et al. 2000, The Astronomical Journal, 120, 1579 Zehavi, I., Contreras, S., Padilla, N., Smith, N. J., Baugh, C. M., & Norberg, P. 2018, The Astrophysical Journal, 853, 84 Zehavi, I., Zheng, Z., Weinberg, D. H., Blanton, M. R., Bahcall, N. A., Berlind, A. A., Brinkmann, J., Frieman, J. A., Gunn, J. E., Lupton, R. H., et al. 2011, The Astrophysical Journal, 736, 59 86 Zehavi, I., Zheng, Z., Weinberg, D. H., Frieman, J. A., Berlind, A. A., Blanton, M. R., Scoccimarro, R., Sheth, R. K., Strauss, M. A., Kayo, I., et al. 2005, The Astrophysical Journal, 630, 1 Zentner, A. R., Hearin, A., van den Bosch, F. C., Lange, J. U., & Villarreal, A. 2019, Monthly Notices of the Royal Astronomical Society, 485, 1196 Zentner, A. R., Hearin, A. P., & van den Bosch, F. C. 2014, Monthly Notices of the Royal Astronomical Society, 443, 3044 Zhang, H.-Y., & Jing, Y.-P. 2004, Chinese Journal of Astronomy and Astrophysics, 4, 507 Zheng, Z., Berlind, A. A., Weinberg, D. H., Benson, A. J., Baugh, C. M., Cole, S., Davé, R., Frenk, C. S., Katz, N., & Lacey, C. G. 2005, The Astrophysical Journal, 633, 791 Zheng, Z., Coil, A. L., & Zehavi, I. 2007, The Astrophysical Journal, 667, 760 Zheng, Z., Tinker, J. L., Weinberg, D. H., & Berlind, A. A. 2002, The Astrophysical Journal, 575, 617 Zhu, G., Zheng, Z., Lin, W., Jing, Y., Kang, X., & Gao, L. 2006, The Astrophysical Journal Letters, 639, L5 Zu, Y., & Mandelbaum, R. 2015, Monthly Notices of the Royal Astronomical Society, 454, 1161 —. 2016, Monthly Notices of the Royal Astronomical Society, 457, 4360 Zu, Y., Mandelbaum, R., Simet, M., Rozo, E., & Rykoff, E. S. 2017, Monthly Notices of the Royal Astronomical Society, 551