The halo bias for number counts on the light cone from relativistic N-body simulations

We present the halo number counts and its two-point statistics, the observable angular power spectrum, extracted for the first time from relativistic N-body simulations. The halo catalogues used in this work are built from the relativistic N-body code gevolution, and the observed redshift and angular positions of the sources are computed using a non-perturbative ray-tracing method, which includes all relativistic scalar contributions to the number counts. We investigate the validity and limitations of the linear bias prescription to describe our simulated power spectra. In particular, we assess the consistency of different bias measurements on large scales, and we estimate up to which scales a linear bias is accurate in modelling the data, within the statistical errors. We then test a second-order perturbative bias expansion for the angular statistics, on a range of redshifts and scales previously unexplored in this context, that is 0.4 ≤ z̅ ≤ 2 up to scales ℓ max ∼ 1000. We find that the angular power spectra at equal redshift can be modelled with high accuracy with a minimal extension of the number of bias parameters, that is using a two-parameter model comprising linear bias and tidal bias. We show that this model performs significantly better than a model without tidal bias but with quadratic bias as extra degree of freedom, and that the latter is inaccurate at z̅ ≥ 0.7. Finally, we extract from our simulations the cross-correlation of halo number counts and lensing convergence. We show that the estimate of the linear bias from this cross-correlation is consistent with the measurements based on the clustering statistics alone, and that it is crucial to take into account the effect of magnification in the halo number counts to avoid systematic shifts in the computed bias.


Introduction
At present and in the near future, many very large galaxy surveys are being and will be performed. The Dark Energy Survey (DES) recently released the Year 3 data products [1]; the Dark Energy Spectroscopic Instrument (DESI) [2] completed its Survey Validation observations, in preparation for the 5-year survey that just started. Future experiments such as the Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST) [3,4], the Euclid satellite [5,6], the Spectro-Photometer for the History of the Universe, Epoch of Reionization, and Ices Explorer (SPHEREx) [7], and the Square Kilometre Array (SKA) [8] are designed to map the galaxy distribution in the Universe on the largest scales, from early until present time.
One of the main challenges in cosmology is to develop a robust theoretical model for our observables that will allow us to extract information from ultra-large scales down to intermediate and small scales that we will observe with stunning precision, thanks to the vast volume of observational data available in the near future. In the concordance model for cosmology, the galaxy distribution that we observe originates from primordial perturbations generated in the early Universe, see e.g. [9]. Therefore, in order to correctly interpret observations from a galaxy survey, we need three key ingredients: i) an accurate description for the evolution of the matter density field and the geometry of the Universe, ii) an accurate description of the relation between the matter distribution and the distribution of galaxies (or other observable tracers), and iii) a way to account for all observational effects, including light-cone projection, selection effects and other systematics possibly present in our data. The first ingredient can be obtained either with analytical methods, within the framework of perturbation theory [10] and the Effective Field Theory of Large Scale Structure [11][12][13], or with numerical methods by running cosmological simulations, where the continuous dark matter distribution is represented via discrete N-body particles [14,15]. Concerning the second point, the relation between matter and galaxies is encoded in the so-called 'galaxy bias' [16]. The galaxy bias depends on the complex process of galaxy formation [17]. However, it is well established that galaxies reside inside dark matter halos -virialized objects that form through gravitational instability. Dark matter halos also constitute biased tracers of the large-scale dark matter distribution. Galaxy and halo bias can be studied and modelled using similar methods, with the simplification that halos are composed of dark matter and therefore baryonic effects do not affect halos to the same extent as galaxies. The third point encodes different effects. One important observational effect is that we perform our observations on our past light cone, therefore relativistic corrections such as redshift-space distortions [18] and lensing magnification [19] need to be modelled properly.
In this context, some of us have studied the angular statistics of the weak-lensing observables (convergence, ellipticity and rotation) in Ref. [57], and the number counts of N-body particles in Ref. [58], extracted directly from high-resolution N-body simulations performed with the relativistic code gevolution [49,50]. In this work we extend the analysis of Ref. [58] to study the relativistic number counts of dark matter halos. Dark matter halos are biased tracers of the LSS. Therefore, in this work, we tackle one of the main difficulties in interpreting galaxy surveys: the biasing problem.
In this paper, we present a detailed analysis of the halo bias, extracted directly from the observable angular power spectrum. Halos are collapsed density peaks of matter and therefore much easier to model than actual galaxies. We work with a complete population of halos and do not consider additional observational biases that would come with the selection function of a realistic galaxy survey, see e.g. [63][64][65]. For the first time, we investigate the nonlinear bias in the observable statistics in directly observable harmonic space, using relativistic simulations to ray-trace our past light cone. We test the validity and limitations of the linear bias model, and we compare our simulation data to a second-order perturbative bias expansion. Furthermore, we study the interplay of nonlinear biasing and redshift-space distortions (RSD) and the halo bias in the cross-correlation of number counts and lensing convergence. The angular harmonic space has the advantage that is constructed from the truly observable coordinates (redshift and angular positions), while in a 3D analysis in Fourier or configuration space one needs to assume a cosmological model to convert the observed redshifts and angular positions into distances. On the other hand, it is not convenient to use if the redshifts of the objects are known very precisely: in order to fully exploit the redshift information, one would need to split the sample in many thin redshift bins and run a full tomographic analysis (i.e. including also cross-correlations between unequal-redshift bins). This would make the analysis computationally prohibitive and would lead to large shot noise per redshift bin. For this reason, the specific application we have in mind for this analysis are photometric galaxy surveys, such as DES [1], LSST [3,4] and Euclid [5,6], that will extract cosmological information from a tomographic analysis of the LSS, cross-correlating galaxy clustering and cosmic shear (3 × 2pt analysis). Presently, the linear bias approximation is widely adopted both in the data analysis (e.g. in the DES fiducial analysis [1]) and in cosmological forecasts [66][67][68]. Our work is a first step towards a better understanding of biasing in the observable angular power spectrum, with the aim of improving the currently available models and helping the scientific community to extract information from a wider range of scales. For a photometric survey the radial resolution is 50 h −1 Mpc so that RSD which are purely radial can be treated with linear theory and velocity bias is not relevant. However, for density fluctuations which are studied also on small transversal scales a treatment beyond the linear bias is necessary as we shall demonstrate.
The rest of this paper is structured as follows. In Sec. 2 we describe the N-body code, the simulations, the halo catalogues, and the ray-tracing method used to obtain the observed positions of halos and particles on the light cone. In Sec. 2.1 we discuss the method we apply to construct the sky maps of number counts, and in Sec. 2.2 we describe our estimator to extract the angular power spectrum and its covariance. In Sec. 3 we test the linear bias model against our simulation data. In particular, in Sec. 3.1 we study the consistency of different methods to estimate the linear bias, comparing measurements on the light cone from the angular power spectrum and measurements on the snapshots, coming from the power spectrum in Fourier space. In Sec. 3.2 we study up to which scales we can robustly model the simulated spectra with the linear bias prescriptions. In Sec. 4 we explore minimal extensions of the linear bias model for the angular power spectrum, and we show how they improve the modelling of our simulation data. In Sec. 5 we discuss the cross-correlations of halo number counts and the lensing convergence. We show how to extract the linear bias from them, and we assess the importance of lensing magnification. In Sec. 6 we summarize our results and sketch future developments of this work. Additional details about our analysis are provided in the appendices.

The simulations
The simulation used in this work is the unity2 simulation described in Ref. [58] which has the following properties. It contains 5760 3 N-body particles in a volume of (4032 Mpc/h) 3 , which yields a mass resolution of 3 × 10 10 M /h. The grid resolution is 700 kpc/h which is equal to the mean particle separation, i.e. the simulation has the same number of N-body particles as grid points. The baseline cosmology is a ΛCDM model with two massive and one massless neutrino species. The massive states m 2 , m 3 are in normal hierarchy, with masses m 2 = 0.008689 eV, and m 3 = 0.05 eV. The fiducial cosmological parameters are based on the Planck 2013 analysis [69]. They are fixed to the values Ω cdm = 0.26858, Ω b = 0.049, h = 0.67, T CMB = 2.7255 K, A s = 2.215 × 10 −9 and n s = 0.9619. The initial conditions are generated from the linear transfer functions produced with the code class.
The simulation uses the public version 1.2 1 of the relativistic N-body code gevolution. This code solves Einstein's equations together with the geodesic equation for the particle ensemble, and evolves all six degrees of freedom of the metric in Poisson gauge under the weak-field assumptions. All scalar, vector and tensor perturbations are included at first order, while the first and second spatial derivatives are included to all orders for the scalar potentials. Since we are in ΛCDM cosmology, relativistic corrections are expected to be extremely small and our results should be in good agreement with a Newtonian calculation. Details on the code can be found in Refs. [49,50]. Our observer is chosen to occupy one of the corners of the periodic domain, and the data on the light cone are stored on the fly. Comoving positions and peculiar velocities of the particles are recorded when their trajectories intersect the background past light cone, while the metric perturbations are stored on spherical pixelised maps, that are concentric around the observer and spaced by a radial separation equal to the grid resolution. The pixelisation is handled by the HEALPix library [70]. These data are used by a post-processing routine that integrates the null geodesic equations backwards in time from the observer to each target object. The null geodesic equations are solved using a scheme which is non-perturbative for the scalar part of the gravitational potential φ. The gravitomagnetic vector potential and the transverse and traceless tensor modes are neglected in this post-processing step for simplicity, but they could easily be included. The method is identical to the one described in detail in Ref. [57], with the only difference that the target objects are not N-body particles but halos. In the context of ray tracing, the fact that we use a relativistic simulation has the benefit that the calculation is automatically self-consistent, and we do not need to consider subtleties of interpreting the data as in Newtonian simulations [71].
The halo population that we analyse in this work is produced from the particle data on the light cone using the halo finder rockstar [72]. This algorithm first identifies friends-of-friends groups in sixdimensional phase space and then computes halo properties based on spherical overdensity criteria in threedimensional comoving position space. Perturbations to the metric (i.e. relativistic effects) are ignored in this step, but these would only cause corrections of order φ ∼ 10 −5 to the matter density, similar to the kinetic energy ∼ v 2 which is also commonly ignored. However, other errors to the mass measurement are . The total number of halos (sources) for these catalogues are N sources ≈ 5.9×10 7 (unity2-lowz) and N sources ≈ 1.1×10 7 (unity2-highz) for the halo population with masses above M min = 10 12 M /h, whereas N sources ≈ 1.3 × 10 7 (unity2-lowz) and N sources ≈ 1.6 × 10 6 (unity2-highz) for the halo population with masses above M min = 5 × 10 12 M /h. The sky fractions of the catalogues are denoted by f sky . expected to be much larger. Two halo catalogues are produced for the same observer. For z 0.85 we have a catalogue with full-sky coverage , labelled unity2-lowz. At higher redshift, that is 0.85 z 3.5, a second catalogue, labelled unity2-highz, covers 1932 square degrees which correspond roughly to 5% of the sky.
From the full catalogues produced by rockstar we first remove all subhalos and then apply a selection based on halo mass. For the latter, we use the property M 200b as a proxy, which is the mass enclosed within a spherical overdensity of 200 times the background density. We note that the mass estimates from our simulation are biased low because of the limited resolution, and the numerical value of M 200b should therefore not be taken literally. Its value is not important, but rather the fact that we are selecting the most massive objects according to some threshold for some mass proxy.
In Fig. 1 we show the normalized redshift distribution for unity2-lowz and unity2-highz for two mass thresholds, M min = 10 12 M /h (red) and M min = 5 × 10 12 M /h (blue). For the more restrictive threshold M min = 5 × 10 12 M /h, a larger percentage of sources are found at lower redshift. Most of the results presented in this manuscript are based on the catalogue with that threshold.

From the halo catalogue to sky maps of number counts
From the halo catalogues described above, we extract sky maps of number counts in tomographic redshift bins. We follow a similar procedure as in Ref. [58], where it has been applied to extract the number counts of N-body particles, that is:

4.
When unity2-highz is analysed, we apply a Boolean mask on the map according to the fraction of sky f sky covered by the pencil beam.
In our analysis, the chosen half-bin width is σ z = 0.05 (1+z), which is similar to the typical redshift resolution of a photometric survey. Fig. 2 shows sky maps of number counts extracted from the halo catalogues and the corresponding particle catalogues, at two reference redshifts for unity2-lowz and unity2-highz.

Estimation of the angular power spectrum and its covariance
The angular power spectra are computed from the spherical harmonics coefficients of the map. For a full-sky map we use the HEALPix estimator anafast, while for maps that cover a smaller fraction of the sky we run the code PolSpice [73,74] which corrects for the effect of the mask.
Assuming that perturbations are Gaussian, the covariance of the angular power spectrum is where for i, j, p, q run over the list of cross-correlated maps, C est is the angular power spectrum estimated from the number counts maps, and ∆ is the multipole bandwidth, and f sky is the fraction of sky. The autocorrelations are strongly affected by shot-noise, i.e. C est,ij = C ij + δ ij /N i , for a redshift bin with N i objects per steradian. For the autocorrelation of a map, this translates into the error In order to remove the shot noise, we estimate the angular power spectrum at a fixed redshift bin using the jackknife resampling technique. We randomly split the objects in our redshift bin into two subsamples with roughly equal number of sources, we build maps of the number counts for the two sub-samples separately, and we estimate the autocorrelation as the cross-correlation of the two maps, see e.g. [75]. Using this method, we largely reduce the impact of the shot noise because, even if the sub-samples have half as many galaxies, their cross-correlation is not affected by shot noise coming from self-pairs.
We estimate the error on the angular power spectrum from Eq. (2.1), that is Comparing Eq. (2.2) and Eq. (2.3), we see that the jackknife estimator removes the shot noise from the mean, but it slightly increases the statistical fluctuations around the mean. In Appendix A, we show the effectiveness of the jackknife method to remove the noise from our measurements.

The linear bias model in harmonic space
In full generality, the relation between a biased field, such as the halo overdensity, and the underlying matter overdensity is an intricate function of the nonlinear density and tidal fields and their derivatives, see for example Ref. [16]. On large scales, where linear perturbation theory still remains valid, this relation simplifies to the linear bias model, i.e.
where δ h is the local halo overdensity, δ is the density contrast and b 1 is the linear bias, which simply depends on redshift. Note that δ h and δ are defined in a matter gauge and that a gauge correction appears when one uses the density contrasts defined in a different gauge [76]. The linear bias model has been widely used in the literature, both in data analysis and forecasts for future cosmological surveys [1,[66][67][68]. Nevertheless, it is well known that on small scales this model is not accurate and, for this reason, a stringent scale cut must be applied in order to avoid systematic biases in the analysis. While the validity of the linear bias model has been tested in Fourier space for the galaxy power spectrum and in configuration space for the 2-point correlation function, an analysis based on the angular power spectrum is lacking, and our work aims to fill this gap. With this purpose in mind, we extract the angular power spectra in tomographic redshift bins from unity2-lowz and unity2-highz. The angular power spectra are binned in band powers with ∆ = 5 for the full sky catalogue unity2-lowz and ∆ = 20 ≈ f −1 sky for unity2-highz. The angular power spectra are compared to their theoretical prediction, computed with the class code which implements the number counts in linear perturbation theory including relativistic effects. In our theoretical model, we include nonlinearities in the matter density using the HMCODE recipe implemented in class [77], with the model parameters fitted to the Cosmic Emulator dark matter only simulation [78].
In Sec. 3.1 we estimate the bias by comparing the halo and the particle number counts, while in Sec. 3.2 we extract the linear bias by fitting the halo angular power spectrum to our theoretical model, and we study at which scales the linear bias approximation breaks down.   green, orange and blue data points are the measurements from the autocorrelation of the particle maps, the crosscorrelation of halo and particle maps, and the autocorrelation of halo maps, respectively. The different amplitudes of three curves are due to the halo bias. Dashed lines represent the prediction from class using HMCODE for the angular power spectrum of matter (in black) and the biased angular power spectra, where the bias has been estimated through a fit as described in the text. Bottom panels: Bias functions constructed from the autocorrelation of halo maps (blue) and the cross-correlation of halos and particles (orange). A horizontal dashed black line shows the are the best-fit linear bias from auto and cross-correlation. The black vertical dotted line represents the scales where finite-grid resolution effects suppress the spectra at the level of 5%. As in Ref. [58], this is estimated to be 5(z) where Nyq is the multipole corresponding to the Nyquist frequency of the simulation. In the right panel, this line is not visible because atz = 1 it falls beyond the range of multipoles shown in the plot.

Linear bias: Fourier space versus harmonic space
As a first consistency check, we compare the linear halo bias estimated from two different statistics: the power spectrum in Fourier space, and the angular power spectrum in harmonic space. For both measurements, we apply the same mass selection, excluding from the halo catalogue all sources with masses below M min = 5 × 10 12 M /h.
We estimate the linear bias in Fourier space from the halo and the particle snapshots. The k-dependent bias in Fourier space is computed by taking the ratio of the cross-spectrum of halos and matter P hm with the matter power spectrum P m , at a fixed redshift, that is The linear bias is extracted by fitting this function to a constant in the range of scales between k min = 0.02 h/Mpc and k max = 0.08 h/Mpc. Note that k min is some two orders of magnitude larger than the horizon scale, so that the gauge dependence of the power spectrum, which is relevant mainly on scales of the order of the horizon and beyond, has a negligible effect. More details on the estimation of the linear bias from the snapshots are given in Appendix B.
In order to extract the linear bias from the statistics in harmonic space, we neglect here light-cone effects. We construct maps of halo number counts and particle number counts using the comoving positions of the sources and their background redshift. In a given redshift bin, we estimate the autocorrelation of the map of halo number counts, the autocorrelation of the map of particle number counts, and their cross-correlation. We construct the bias functions C theory fit density C theory fit GR C halos x particles C halos x halos P(k) halos x particles Figure 4: Comparison of different measurements of the redshift-dependent linear bias, for halos with masses above M min = 5 × 10 12 M /h. P (k) denotes measurements from power spectra on snapshots of a simulation in Fourier space (see Appendix B), whereas C denotes measurements from angular power spectra. In both cases, we show a measurement using the cross-correlation "halos × particles" (dashed blue and solid orange, respectively), and for C we additionally show a measurement using the auto-correlation "halos × halos" (green dashed). Instead of using the simulated number counts of particles as reference, we also constrain the bias using only the halo counts and a theoretical model that either accounts for local density fluctuations only, "theory fit density" (solid black), or includes the effects of redshift-space distortions and lensing on the light cone, "theory fit GR" (dashed grey), as explained in detail in Sec. 3.2.
where we remove the shot noise from the angular power spectra with the jackknife method described in Sec. 2.2 and Appendix A.
On large scales these -dependent functions are well described by a constant bias that we estimate with a fit in the range of multipoles ∈ [20, max ], where max = 80, 130 for unity2-lowz and unity2-highz, respectively. We apply a different cut at small scales for the two catalogues to ensure that enough data points are included in the fit for unity2-highz. By excluding low multipoles < 20, which are also most affected by sample variance, any residual gauge dependence of the calculation should be very small. We have verified that in the selected range of multipoles the linear bias model provides a good fit.
While the Fourier space analysis measures the linear bias at a fixed conformal time (or background redshift), the angular statistics measures the correlation of the overdensity fields projected onto a radial redshift bin:δ where dN/dz is the redshift distribution of the sources, while W encodes the redshift selection we apply to our catalogue. In our case, we apply a simple sharp cut in redshift (tophat), selecting all sources with background redshift in the range z ∈ [z − σ z ,z + σ z ]. The angular power spectra strongly depend on the bin size σ z , since radial modes are washed out for wide redshift bins. However, both the halo and the matter density fields are affected in the same way by the radial projection and, at least on linear scales, we have verified that the bias function depends only weakly on the bin size. Therefore, we expect the bias estimation from the harmonic and Fourier space spectra to be consistent with each others. In Fig. 3 we show the angular power spectra for the halo number counts, the particle number counts, and their cross-correlation, atz = 0.5 ( Fig. 3a) andz = 1 (Fig. 3b). The bottom panels show the functions b auto (blue) and b cross (orange), from which we extract the linear bias in harmonic space. For the scales below ∼ 150 the bias functions are roughly constant, and b auto overlaps with b cross as expected. Therefore, we estimate the best-fit linear bias b auto/cross 1 from auto and cross-correlation by fitting b auto/cross to a constant in the large-scale regime. On smaller scales, nonlinear biasing manifests itself with an excess of power in the simulation data, compared to the prediction from the linear bias model. Furthermore, we see that the bias functions b auto and b cross increase due to higher-order corrections in the bias expansion that affect auto and cross-correlations in different ways.
In Fig. 4 we compare the linear bias estimated from the snapshots using the Fourier space statistics, P (k), to the measurements in angular harmonic C , using either the autocorrelation "halos × halos" of halo maps or the cross-correlation "halos × particles." The measurements in harmonic space from the halo autocorrelation and the cross-correlation "halos × particles" give consistent results up toz ∼ 2. At higher redshift, the small discrepancy between the measurements increases to reach about ∼ 10% atz = 3.
This analysis refers to a halo population with masses larger than 5 × 10 12 M /h. For comparison, we repeat the analysis for a population of smaller objects (M min = 10 12 M /h). The redshift-dependent linear bias from the angular statistics for the two halo populations is shown on the left plot in Fig. 5. As expected, halos with the smaller mass threshold have a smaller linear bias. For the population of smaller objects the discrepancy between the best-fit bias from auto and cross-correlations, shown in the bottom panel of Fig. 5, is smaller compared to the more massive halos. On the right plot in Fig. 5 we show this discrepancy (in percent) as a function of the average best-fit bias from auto and cross-correlation. We see that for both populations the difference increases with the values of the bias in a similar way. As theoretically, measuring the linear bias via the halo-particle correlation or as the square root of the halo-halo divided by the particle-particle correlation should be equivalent, the difference we find for highly biased halos highlights that higher-order bias parameters are important even on large scales for these objects and thus contaminate the measurement of the linear bias. Our results suggest that the impact of these higher-order corrections grows roughly proportionally to the bias. It is also not surprising that linear bias fares worse at higher redshift. While the fluctuations themselves are smaller and therefore more linear at higher redshifts, massive halos are rarer and therefore more strongly biased. This does not only increase the linear bias but also the contribution from higher-order bias terms, for which the simple relation P hm /P m = P hh /P m no longer holds.

Testing the validity of the linear bias model for angular statistics
In the previous section, we estimate the linear bias by simply comparing the halo and the matter statistics. While this procedure is straightforward when the data come from a numerical simulation and particles information is available, the unbiased matter density cannot be observed. Therefore, in real observations the statistics of the biased number counts is fitted to a theoretical prediction, where the bias can be measured as one of the free parameters of the theoretical model.
In this section, we estimate the linear bias of our sample by fitting the simulation data to a oneparameter model. We want to assess whether we can accurately recover the linear bias estimated from the particle distribution using our theoretical prediction, whether the bias estimation is affected by light-cone effects such as redshift-space distortions (RSD), and up to which multipoles the linear bias prescription can be adopted without spoiling the quality of the fit.
The theoretical model used in this section is based on the linear number counts implemented in the class code [37], and the notation for the different contributions to the fluctuation of the counts follows Ref. [58]. We write the full halo number counts as where δ h is the local overdensity of halos, In the linear bias model it is simply proportional to the density contrast δ in comoving gauge. The term δ rsd includes the linear Kaiser RSD [18] and other subdominant Doppler corrections, where v = −∇V is the peculiar velocity in Poisson gauge (i.e. V is the velocity potential), r = τ o − τ is the conformal distance, H = a /a is the conformal Hubble rate and a prime denotes a derivative with respect to the conformal time τ . The evolution bias f evo describes the redshift evolution of the halo number density in a comoving volume. Since f evo affects only subdominant contributions to the number counts, we set it to zero in the analysis. The lensing magnification term [19,79,80] is denoted by δ κ , where ∆ Ω is the Laplace operator on the sphere. The last term in Eq. (3.5), δ gr , includes local and integrated combinations of the Bardeen potentials. These terms introduce a very small correction to the observable on the scales considered here and therefore we neglect them in the theoretical model. While we build our model from the linear computation of the number counts, nonlinearities in the matter density contrast, δ, are included using the HMCODE [77]. The linear bias b 1 is the only free parameter. We estimate it using the following procedures: 1. We analyse the maps of number counts built from the comoving position of the halos and we fit their spectra to the theoretical prediction with only the density term in the number counts, assuming that the linear bias is constant within the redshift bin. Therefore, our model is given by where C δ is the angular power spectrum sourced by only density perturbations for an unbiased tracer.
2. We use the fully relativistic maps, and we fit their spectra to the theoretical prediction which includes all relevant relativistic effects. The theoretical prediction is constructed computing the separate contributions of density C δ , the summed contributions of RSD and lensing magnification C (rsd, δκ) , and their cross-correlation with density for an unbiased tracer C δ×(rsd, δκ) . From those separate terms, our model is given by We fit the model to the data using the Levenberg-Marquardt algorithm, as implemented in the python routine curve fit, that minimizes the quantity where σ are the errors on our data points estimated as described in Sec. 2.2. We fit the simulated data for different values of max in the range max ∈ [50,1000]. For each of these fits, we quantify the quality of the fit from the best-fit Chi-Square χ 2 fit . We estimate the probability to find a χ 2 > χ 2 fit , under the hypothesis that our model is adequate: where PDF χ 2 (x, n dof ) is the χ 2 probability distribution function for n dof degrees of freedom, having n dof = max − min + 1 − n bias ) = ("number of points included in the fit" − "number of free parameters"). As there is only one free bias parameter to be fitted, we have here n bias = 1. The best-fit linear bias is chosen to be the value that, in the range of max considered, gives the largest P (χ 2 > χ 2 fit ). In Fig. 6 we show the angular power spectra extracted from our simulations for four representative redshifts. Simulation data are the blue (density maps) and orange (fully relativistic maps) points, while the best-fit theoretical models are represented by dashed and continuous black lines, respectively. The best-fit values for the linear bias that we obtain for all the redshifts considered here are plotted in Fig. 4 as "theory fit density" (black line), and "theory fit GR" (grey dashed line). We find that the best-fit values obtained by fitting the halo angular power spectrum to the theoretical model are in good agreement with the values measured directly from the cross-correlation with the particle maps. On small scales, the linear bias prescription breaks down; actually, the deviations between our simulation data and the best-fit theoretical model at a fixed scale ∼ 1000 increases with redshift (bottom panels of Fig. 6).
In order to quantify the quality of our fit and to test up to which nl our model describes reasonably well the simulation data, we choose two criteria to assess at which multipole the linear bias model becomes unreliable: (a) χ 2 criterion: the nonlinear cut nl is the minimum max at which P (χ 2 > χ 2 fit ) drops below the 5% threshold.
(b) 5% deviation criterion: nl is the minimum multipole such that the relative difference in percent, 100% × (C sim − C th )/C th , becomes larger than 5% for four and eight consecutive data points for unity2-lowz and unity2-highz, respectively.
The criterion a) takes into account the statistical uncertainty in the simulation data. On the other hand, the level of shot noise and its redshift dependence is a characteristic of the halo catalogue under consideration, and the cosmic variance error strongly depends on the survey geometry considered here. For these reasons, criterion b) gives a more general answer to the question of what are the smallest scales that can be modelled with the linear bias in the angular power spectrum.
In Fig. 7 we show the nonlinear angular scale estimated from our simulation as a function of redshift, for the two criteria described above. The nonlinear angular scale estimated from the density and fully relativistic maps gives very similar results, while the 5% deviation criterion gives a more conservative estimate of nl at large redshift. This is due to the large statistical uncertainty on the simulation measurements, the small fraction of sky available and the high level of shot noise at high redshift.
We compare the nonlinear angular scale from our simulation to the ansatz often used in the literature nl = r(z)k max , see for example Refs. [81,82]. Note that in our application k max is the maximum wavenumber at which the linear bias model is valid in Fourier space, assuming that the nonlinearities of matter are wellmodelled in the matter power spectrum by the HMCODE, while in Refs. [81,82] k max fixes the maximum wavenumber of validity for linear perturbation theory. In practice, for biased tracers, it is not realistic to shows the multipole at which resolution effects become larger than 5%. The corresponding angular scale is beyond the range shown for the other three redshifts. Bottom panels: relative difference between the best-fit theoretical model and the simulation data. A red vertical line shows the maximum multipole that is well modelled by a linear bias, using the 5% deviation criterion described in the text.
separate the nonlinearities of matter from the effect of biasing. Nevertheless, this approximation is often used in cosmological forecasts, see for example [66][67][68]. Therefore, in this section, we assess up to which scales this approach is legitimate. In our simulation results we cannot identify a clear redshift dependence from both our criteria. In unity2-lowz (z < 0.85) we see that the linear bias prescription that we adopt here is not valid anymore for > 200 − 400, while the nonlinear angular scale in unity2-highz (0.85 <z ≤ 3) ranges roughly between 150 < nl < 700 with large fluctuations. The general trend of our numerical results tells us that at high redshift the linear bias description cannot be extended up to smaller scales, unlike the ansatz nl = r(z)k max , and that including scales above nl ∼ 400 in a cosmological analysis without a consistent treatment of nonlinear biasing is too optimistic. We understand this result as follows. Matter non-linearities and biasing have an opposite redshift dependence: nonlinear corrections to the clustering of matter grow in the late-time Universe, while biasing has a larger impact at early times. In the clustering statistics of a biased tracer these two effects are mixed, and they lead to the noisy redshift dependence of non-linearities displayed in Fig. 7. At early times, the halo power spectrum is non-linear since halos are highly biased objects. At late times it is non-linear since matter fluctuations have entered the non-linear regime. This behaviour is certainly mass dependent and halos with much smaller masses will tend to be less biased, especially at higher redshift.

Beyond linear bias
In order to model the halo angular power spectrum accurately on scales > 400, we need to go beyond the linear bias prescription. A consistent formalism to describe nonlinear biasing has been developed in the past decade in the framework of perturbation theory, see for example Refs. [16,[83][84][85][86]. The perturbative bias expansion has been tested with numerical simulations, see e.g. [87], and recently applied to the analysis of cosmological surveys [88][89][90]. Furthermore, there are now public codes that compute nonlinear power spectra of biased tracers at one-loop in cosmological perturbation theory, such as class-pt [91], based on the linear Boltzmann code class, and PyBird [92]. While the perturbative approach to biasing provides a robust framework to model galaxy clustering observables in the mildly nonlinear regime, in full generality it also significantly increases the number of free parameters in the model. This is particularly problematic for a tomographic analysis where we need to marginalize over the extra bias parameters in each of the tomographic redshift bins. Therefore, it is crucial to find a balance between improving our modelling of biasing on small scales and limiting the number of nuisance parameters in order to preserve the information gain. For this reason, simulations have been used in the literature to investigate the relation between the linear bias and higher order bias parameters, see for example [93,94].
In this work, we test the accuracy of minimal extensions to the linear bias for the halo angular power spectrum. Following for example Ref. [88], we consider a second-order bias expansion in the basis of bias operators {δ, δ 2 , G 2 }, where δ is the matter density contrast, and G 2 is the tidal field operator. In comoving position space, i.e. neglecting light-cone effects, the perturbation of the halo number counts in this basis is Given this model, we compute the theoretical prediction for the angular power spectrum in two steps. First, we estimate the biased power spectrum in Fourier space P hh (z, k, µ) using the code class-pt (see Ref. [91] for the details on the implementation). Then, we compute the projected statistics, that is the angular power spectrum, using the flat-sky approximation [95,96]. For autocorrelations and tophat bins with half-bin width σ z , we need to compute 2 where j 0 is the spherical Bessel function, k = k 2 + 2 /r 2 (z) and µ = k /k.
Our theoretical model, based on Eq. 4.2, neglects wide-angle effects and the non-uniform redshift distribution of sources inside the bin. In Appendix C we test the impact of these two approximations on our study, and we find that up toz = 2 our approximation agrees within 5% with the full computation down to = 10. The largest discrepancies are found on large scales, and they are within the statistical uncertainties of our simulation measurements.
The simulation data are fitted to these models with the Markov Chain Monte Carlo (MCMC) method, making use of the Goodman & Weare's affine MCMC invariant ensemble sampler implemented in the python package emcee [97]. In Sec. 4.1 we report the results for the density maps, and in Sec. 4.2 we discuss the analysis of the maps from fully relativistic number counts.

Results from density maps
We first consider the maps of number counts in comoving position space. Neglecting RSD and relativistic effects, the halo power spectrum in Fourier space is only a function ofz and k. The full model for the Fourier space power spectrum, at fixed redshift, is then given by We point the reader to Ref. [91] for the exact form of the terms P 1-loop,SPT , I δ 2 , I G2 , F G2 , I δ 2 δ 2 , I G2G2 , I δ 2 G2 appearing in Eq. (4.3). We test four bias models against our simulation data: a 1-parameter model {b 1 }, that applies a linear bias on top of the linear power spectrum plus one-loop corrections from standard perturbation theory (SPT) P 1-loop,SPT ; a 2-parameter model {b 1 , b 2 } which also includes as free parameter the quadratic bias, but neglects the tidal bias b G2 ; a 2-parameter model {b 1 , b G2 } that instead includes as free parameter the tidal bias, but neglects the quadratic bias b 2 ; and a 3-parameter model {b 1 , b 2 , b G2 } that incorporates the full second-order bias expansion in Eq. 4.1. These models are obtained from Eq. (4.3), simply setting to zero the values of all parameters that we neglect in the corresponding bias expansion.
In our analysis, we fix the cosmological parameters to the input values in the simulation. This makes the fitting procedure very fast, as we can estimate the contributions to the angular power spectrum individually and only once for each redshift bin, solving integrals of the form with T X = {P lin + P 1-loop,SPT , I δ 2 , I G2 , F G2 , I G2G2 , I δ 2 G2 }. Once these terms are pre-computed, the model for the angular power spectrum and arbitrary values of the bias parameters can be quickly evaluated: where we omit the dependence on redshift of the bias parameters, and the (z, σ z )-dependence of the spectra for the sake of brevity. We fit this model to our simulation data for all redshift bins considered in Sec. 3, up toz = 2. For the extended bias models, we include in the fit a range of multipoles between min = 20 and max = 900. Furthermore, we use flat uninformative priors on the bias parameters, with one exception that will be discussed later.  Top panels: Blue data points represent the angular power spectra extracted from the simulated density maps. Red dashed lines are the theoretical from linear theory, while continuous red lines represents the best-fit for the {b1, bG2} model. Bottom panels: Percentage difference between the simulation data and the best-fit theoretical prediction, for linear theory and the four models described in Sec. 4.1.
In Fig. 8 we show the outcome of the MCMC fit at four representative redshifts, while in Table 1 (columns denoted by "density") we quantitatively compare the quality of the fit for the four bias models. As estimators for the quality of the fit, we consider the reduced χ 2 (χ 2 /n dof ) and the probability P (χ 2 > χ 2 fit ) to find a χ 2 > χ 2 fit , under the null hypothesis that we are using the correct model. We consider the best-fit value of the parameters to be the value at the 50th percentile.
At low redshift, that isz = 0.4, 0.5, the best-fit of the four models gives a very similar prediction for the angular power spectrum, and they all reproduce accurately the simulation data in the entire multipole range used for the fit. Surprisingly, at low redshift, even the minimal model {b 1 } provides a reasonably good fit to our data. Therefore, the one-loop corrections together with a linear bias provide a better model to the halo angular power spectrum than the prescription used in the previous section, that is a linear bias applied to the HMCODE matter power spectrum. The matter power spectrum used in the {b 1 } model does not include what is called the 'counterterm'in the literature [91], needed to model accurately the clustering of matter in perturbation theory. Therefore, the fact that the model with linear bias and one-loop corrections in standard perturbation theory works fairly well in the two lowest redshift bins could be simply a coincidence. At these redshifts the nonlinear bias has a similar amplitude but opposite sign as the contribution from the counterterm, and the two contributions roughly cancel. This cancellation does not occur anymore at higher redshift, starting fromz = 0.6, the quality of the fit for the 1-parameter model {b 1 } quickly degrades. Table 1: Accuracy of the best-fit (50th percentile value) for the four bias models described in the text. We report the reduced χ 2 , that is χ 2 /n dof , and the probability P (χ 2 > χ 2 fit ) to find a χ 2 > χ 2 fit , under the null hypothesis that we are using the correct model. Fits with P (χ 2 > χ 2 fit ) < 5% are considered not accurate, and highlighted in red. All the MCMC chains are run using uninformative flat priors, except for the model {b 1 , b 2 , b G2 } atz = 1, where we force b 2 being positive (negative) and report the accuracy of the fit for the two cases.
At higher redshift, this model is not able to reproduce the shape of the power spectrum on a wide range of scales. In order to accommodate the data on small scales it loses accuracy in the large-scale part of the spectrum, as shown in the bottom panels of Fig. 8 forz = 1, 1.4, 2. A similar trend emerges for the model {b 1 , b 2 } starting fromz ∼ 0.7. Overall, we find that both the 2-parameter model {b 1 , b G2 } and the 3-parameter model {b 1 , b 2 , b G2 } can describe the simulation data well in the whole range of redshifts and angular scales considered here, the former having the clear advantage of requiring only one nonlinear bias parameter.
As mentioned previously, we use uninformative flat priors to run the MCMC chains. The only exception is given by the model {b 1 , b 2 , b G2 } atz = 1, where we find that, for uninformative priors, the posterior distribution of b 2 has two local maxima, one at b 2 < 0 and one at b 2 > 0. The resulting 50th percentile value for b 2 does not provide a good fit to the data. Therefore, in the results we present here for this model and redshift bin we impose a flat prior b 2 > 0, which instead gives a good fit to the data. In Table 1 we report the information on the quality of the fit for both prior choices, b 2 > 0 and b 2 < 0.
In Fig. 9 we show the redshift dependence of the best-fit values for b 1 (top panel), b 2 (central panel) and b G2 (bottom panel). For b 1 , we compare the values obtained from models that use a perturbative bias expansion and the results of a linear fit, which only includes data on large scales. We see that the perturbative models providing a good fit on all scales, that is {b 1 , b G2 } and {b 1 , b 2 , b G2 }, find best-fit values of b 1 fully consistent with the measurements on large scales, while both {b 1 } and {b 1 , b 2 } give larger values of the linear bias. This means that these two models inaccurately predict a lower amplitude for the angular power spectra at low multipoles, which results in a larger value of the bias and an inaccurate fit. We also see from Fig. 9 that the models {b 1 , b G2 } and {b 1 , b 2 , b G2 } consistently find negative values for the parameter b G2 , with amplitude monotonically increasing with redshift. On the other hand, we also notice that in the 3-parameter model b 2 is mostly unconstrained atz > 1, and the best-fit values are consistent with zero

(z)
Linear bias -full sky Linear bias -flat sky  within the error bars. It seems that the quadratic contributions from the density are already well modelled by one-loop perturbation theory and introducing tidal bias is sufficient to model the data.

Results from relativistic maps
In this section, we apply the same methodology discussed in Sec. 4.1 to test the perturbative bias expansion against the fully relativistic number counts, and thus assess the robustness of our results when light-cone effects are taken into account. In our model, we neglect the effect of magnification and other relativistic effects. Using the class code we tested the impact of magnification in the autocorrelation of a linearly biased tracer, with bias given by the linear bias of our halo population. We find that forz ≤ 2 the impact of magnification is always below 1% and thus it would not affect the analysis in a significant way if included. RSD are included in the model using the simple linear Kaiser prescription [18]. The motivation for this choice is twofold. First, adding nonlinear RSD would extend the set of free parameters. Second, we know that RSD affects the angular power spectrum on large scales, where our simulation data have larger statistical errors due to cosmic variance, and that the effect is highly suppressed for wide redshift bins. While it was pointed out in Refs. [58,96,98] that nonlinear velocities affect the angular power spectrum on large scales for thin redshift bins, in our study we use relatively wide bins where the linear prescription for RSD is expected to work well. Therefore, to include the linear RSD model, we simply add the following two extra terms to Eq. (4.5), and , (4.8) Linear theory Linear theory density unity2-highz GR 10 1 10 2 10 3 100 Linear PT b1, bG2 (d) Figure 10: Outcome of our MCMC fit for the fully relativistic angular power spectra, at four representative redshifts. Top panels: Blue data points represent the angular power spectra extracted from the simulated maps. Red dashed (black dotted) lines are the theoretical prediction from linear theory including (neglecting) RSD, while continuous red lines represents the best fit for the {b1, bG2} model. Bottom panels: Percentage difference between the simulation data and the best-fit theoretical prediction, for linear theory and the four models described in Sec. 4.2.
where f (z) is the growth rate. We apply the same pipeline as described in Sec. 4.1 to the relativistic angular power spectra, using the updated model described above. In Fig. 10 we show the outcome of the fit for four representative redshifts, while the accuracy of the fit for the different models is reported in Table 1 (columns denoted by "full GR"). The analysis of the relativistic maps leads to the same results highlighted in the previous section: {b 1 } and {b 1 , b 2 } are inaccurate at redshiftz > 0.6, 0.7, as they are not able to fit the large-scale and small-scale parts of the spectra simultaneously; on the other hand, {b 1 , b G2 } and {b 1 , b 2 , b G2 } provide a good model up to scales ∼ 1000. We conclude that the model {b 1 , b G2 } is preferable for providing an accurate prediction with the minimum number of free parameters.
In Fig. 11 we compare the best-fit values of the bias parameters from the analysis of the density maps and the fully relativistic maps for the model {b 1 , b G2 }. We find that the bias measurements from the relativistic maps are fully consistent with the measurements from the density maps, and the linear bias agrees with the large scale bias estimated in Sec. 3.2. The agreement between the results with/without relativistic effects is not surprising: as we already mentioned, radial correlations are smeared over a wide bin. Therefore, the contribution of RSD is highly suppressed for the binning chosen in our work, and confined to scales with large cosmic variance. We expect this result not to hold for thin bins. Since this  Figure 11: Best-fit bias parameters for the model {b 1 , b G2 }. Blue data refer to analysis of the density maps, orange data refers to the analysis of the fully relativistic maps. For comparison, the black line shows the bias in linear theory, estimated only from the large scales in the angular power spectrum.
requires a separate analysis, we will address this aspect in future work.

Bias from cross-correlations
In the previous section, we analysed the correlations of number counts at equal redshiftz. Here we consider the cross-correlation between halo number counts and lensing convergence κ. This cross-correlation provides an alternative way to estimate the halo bias when we correlate maps at unequal redshift: the images of objects at z 2 are lensed by the halo and matter overdensities at z 1 < z 2 . Therefore, the convergence field at z 2 is correlated with the number counts at z 1 . By measuring the cross-correlation of the convergence field in our simulation with the halo number counts in comoving position space, C δ h κ (z 1 , z 2 ), and similarly the cross-correlation with the particle number counts, C δκ (z 1 , z 2 ), we can estimate the linear bias b 1 at z 1 as When considering the fully relativistic number counts, that is what we observe in galaxy surveys, there are extra terms that contribute to the halo and particle number counts. On top of RSD, also lensing contributes to the number counts through magnification δ κ , and this effect is highly correlated with the lensing convergence field because in linear theory δ κ = −2κ. Since the convergence is an integrated effect, the convergence maps at different redshifts are correlated. In this section, we study the impact of the magnification in the cross-correlation of number counts and convergence, and the effect on the linear bias estimation in entails.
We extract the convergence maps from our fully relativistic simulations for the unity2-highz catalogue and the binning in redshift considered in the previous sections. To this end we use the particle catalogue in order to have a good sampling of the field. Following the methodology described in detail in Ref. [57], we first use ray tracing to compute the area distance at the position of each particle. Then, we estimate the convergence as fluctuation of the area distance compared to its background value. Finally, we construct the map by taking in each pixel the average of the convergence over all the particles that fall into the pixel.
We compute the cross-correlation of the convergence maps with the halo number counts, C ∆ h κ (z 1 , z 2 ), and particle number counts, C ∆κ (z 1 , z 2 ), including all relativistic effects. We also estimate the correlation of the convergence field C κκ (z 1 , z 2 ).
In Fig. 12 (top panel on the left plot), we show the cross-correlation power spectra C ∆ h κ (z 1 , z 2 ) (blue), C ∆κ (z 1 , z 2 ) (orange), and the amplitude of the contribution from magnification to these spectra, that is Including (z1) (z2) Figure 12: Cross-correlation of number counts at z1 = 1 and z1 = 2 (left and right plots, respectively) with the convergence map at z2 = 3. Top panels: angular power spectra C ∆ h κ (z1, z2) (blue), C ∆κ (z1, z2) (orange), and the amplitude of the contribution from magnification to these spectra, that is 2C κκ (z1, z2) (green). Data points are the spectra extracted from the simulations, the dashed lines represent the linear theory prediction from the class code. Bottom panels: Bias estimation at z1 from this cross-correlation, in a case where we ignore the contribution of magnification (red lines) and in a case where we subtract that contribution (blue lines). The data points are computed from Eq. (5.2) and Eq. (5.3), respectively. The horizontal solid red/blue lines mark the respective linear bias estimated from the cross-correlation. A dashed black lines shows the linear bias estimated previously in Sec. 3.1 from the autocorrelation. For z1 = 2 (right plot) the contribution from magnification is comparable to the one from density on all scales in the spectrum C ∆κ (z1, z2), and thus it is not possible to estimate the linear bias in a meaningful way when neglecting this effect.
2C κκ (z 1 , z 2 ) (green), for z 1 = 1 and z 2 = 3. We estimate the bias at z 1 for a case where we neglect the presence of relativistic effects, that is and also for a case where we subtract the effect of lensing magnification, i.e.
The results are shown in Fig. 12 (bottom panel on the left plot). In the case where we neglect the presence of relativistic effects, the linear bias fitted to b , including scales below = 150, overestimates the value obtained in the previous sections from the autocorrelation. However, in the case where we subtract the contribution of magnification, we are able to recover the value of the linear bias accurately. For the configuration z 1 = 1 and z 2 = 3, the contribution of magnification to the particle angular power spectrum ranges between ∼ 20 − 70% on scales ∼ 30 − 1000, and the shift in the bias estimate is ∼ 20%. When we cross-correlate number counts at z 1 > 1, the contribution of magnification increases, while the density term becomes smaller. At z 1 = 2, magnification constitutes a significant fraction of both the spectra C ∆ h κ (z 1 , z 2 ) and C ∆κ (z 1 , z 2 ), as shown in Fig. 12 (right plot). Therefore, for these configurations, it is not possible to obtain a reasonable estimate of the linear bias neglecting the contribution of magnification. However, when magnification is taken into account in the bias measurement, we can recover the linear bias obtained from the autocorrelation at z 1 = 2.
In Fig. 13 we show the same cross-correlations, for number counts at z 1 = 3 and convergence at z 2 = 3 (left plot) and z 2 = 1 (right plot). For z 1 ≥ z 2 , the contribution from magnification is by far the dominant one in the cross-correlations C ∆ h κ (z 1 , z 2 ) and C ∆κ (z 1 , z 2 ), and there is no information on the halo bias. While we cannot measure the halo bias from these cross-correlations, they are a clean observable to detect cosmic magnification. A first detection of this effect from galaxy shear measurements cross-correlated with  Figure 13: Cross-correlation of number counts at z1 = 3 and the convergence map at z2 = 3 (left plot), and z2 = 1 (right plot). The data points show the measurements from simulation for C ∆ h κ (z1, z2) (blue), C ∆κ (z1, z2) (orange), and the amplitude of the contribution from magnification to these spectra, that is 2C κκ (z1, z2) (green). These three lines coincide within statistical fluctuations because magnification is the dominant contribution to both C ∆ h κ (z1, z2) and C ∆κ (z1, z2), for either of the configurations. The dashed black lines show the predictions from linear theory for C ∆κ (z1, z2), while the grey shaded region shows the 1σ statistical error for C ∆ h κ (z1, z2) around the theory prediction. This error is dominated by the shot noise of halos atz1 = 3.
background galaxy counts has been presented in Ref. [80]. In the near future, galaxy surveys like Euclid [5,6] and the Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST) [3,4] will be able to detect this effect with high significance and extract cosmological information from it [99].
In Fig. 12 and Fig. 13, we notice that the correlations dominated by magnification exhibit a power suppression on small scales when compared to the theoretical prediction. This is an artifact of finite resolution that occurs because magnification and lensing convergence are integrated effects and, therefore, the large error at low redshifts propagates into the estimated spectra at all redshifts. A detailed study of this issue can be found in Ref. [57], Appendix C.

Conclusions
In this paper we have computed for the first time the halo number counts on the light cone of a relativistic N-body simulation and the model independent two-point statistics, the angular power spectrum. Since halos are biased tracers of the LSS, this study is a first step towards building a robust model for the galaxy bias in the angular power spectrum of galaxies, the primary observable in the galaxy clustering analysis of photometric surveys.
We have measured the redshift dependence of halo bias for a halo population with masses above M min = 5 × 10 12 M /h, using different methods and estimators. We have computed the halo bias from the snapshots using cross-correlation of the density field for halos and particles, together with the matter power spectrum in Fourier space. We have compared these results to the bias values measured on the light cone using the angular power spectrum. The bias on the light cone has been extracted in tomographic redshift bins a) by comparing the correlations of halo density maps, particle density maps and their cross-correlation; b) by fitting the halo spectra to a theoretical model, both with and without relativistic effects. We have found that these different estimates give consistent results for the linear halo bias, i.e. the halo bias on large scales. On small scales, high order corrections become significant and the bias from the halo autocorrelation differs from the one obtained using the cross-correlation of halos and particles.
We have investigated the limitations of the linear bias model for our halo population. In particular, we have established at which multipoles a theoretical model with linear bias, and nonlinearities modelled in the matter power spectrum with the HMCODE, breaks down. This nonlinear multipole scale estimated in our simulation ranges between nl = 200 − 600, and it does not show a clear redshift dependence. The values of the nonlinear multipole scale that we obtain from our simulations are roughly compatible with the ansatz often used in the literature nl = r(z)k max , with k max = 0.1 − 0.2 h/Mpc. However, unlike this simple prediction, we find that at high redshift the linear bias description cannot be extended up to smaller scales. This reflects the fact that the bias becomes larger at high redshift, while nonlinearities of matter decay at early time. Nonlinear biasing mixes these two effects, thus it is not expected to decay at high redshift in the same way as the nonlinear clustering of matter does.
We have tested minimal extensions of the linear bias model, including second-order terms in the bias expansion. These extra terms are known in the literature as quadratic bias and tidal bias. We have found that a two-parameter model with linear bias and tidal bias (and quadratic bias set to zero) is able to fit the halo angular power spectra well at equal redshifts on a wide range of scales, up to ∼ 1000. This represents a great improvement compared to the linear bias model. We have also shown that this model is significantly preferred compared to adding a quadratic bias. A quadratic bias addition is not able to reproduce nonlinearities in our simulated angular power spectra atz ≥ 0.7.
Finally, we have shown that it is possible to estimate the linear bias from the cross-correlation of halo number counts and convergence, provided that the convergence field is at higher redshift than the number counts, and that the contribution of magnification is properly accounted for.
This work can be extended in several directions. We have tested the bias expansion against the halo angular power spectra at equal redshift. However, a significant amount of information can be extracted from the correlations at unequal redshift, see for example Ref. [38], and thus it is crucial to assess the robustness of a model against the full set of auto and cross-correlations. A detailed study of the cross-correlation beyond the linear bias model might require going beyond the flat-sky approximation that we have employed in our analysis as it provides a fast computation of our observable. We will address this issue in future work. We have found that RSD can be modelled reasonably well using Kaiser's linear prescription for the binning considered in our work. We expect this picture to change for thinner redshift bins, as several works have pointed out that nonlinear velocities contribute to RSD on large scales in the angular statistics for a typical bin size σ z ∼ 0.01 or smaller, see e.g. Refs. [58,96,98,100]. A detailed analysis of nonlinear RSD in the angular statistics for biased tracers is left for future research. In Fig. 14 we compare the jackknife estimator (JK) to the standard estimator of the spectrum (auto) that simply uses the autocorrelation of the number counts as is, for both halos and particles. We see that the standard estimator is affected by shot noise on all scales for halos, while the spectra of particles are visibly affected only on small scales. This is fully consistent with the level of Poisson noise, given by 1/N i for a redshift bin with N i galaxies per steradian, that is shown in the plot with a dashed horizontal line for halos, and a dotted horizontal line for particles. The noise level is ∼ 2 orders of magnitudes smaller for particles due to their higher average number density. The JK method provides an effective way to remove the shot noise contribution from the spectra, without assuming a prior knowledge of the noise. Fig. 14 shows that the JK method gives consistent results with the angular power spectrum extracted from the autocorrelation if the Poisson noise is manually subtracted from the spectrum.

B Linear bias on comoving slices
In order to measure the linear bias from the Fourier modes on equal-time hypersurfaces we run a separate simulation with the same parameters and mass resolution as unity2 but with a smaller volume of L 3 box = (1344 Mpc/h) 3 . The snapshots then have 1920 3 particles which is easier to handle and allows us to store them at a large number of redshift values z ∈ {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2.2, 2.6, 3}. The halo-matter cross-power spectrum and the matter auto-power spectrum in each snapshot is estimated using the publicly available code Pylians3. 4 This code implements a routine that deposits particles into a regular 3D grid with N 3 voxels using the cloud-in-cell mass-assignment scheme. Here, we have used a grid with N = 1000, which is sufficient for the intended measurements of the linear bias well below its Nyquist wavenumber of k Nyq = πN/L box ≈ 2.33 h/Mpc. We also note that L box is much smaller than the horizon scale and the gauge dependence of the measured power spectra is therefore expected to have a minuscule effect on our bias measurements. The auto-power spectra and cross-power spectra are then estimated from the density fields in Fourier space. The linear bias is finally estimated by dividing the halo-matter crosspower spectrum by the matter auto-power spectrum, see Eq. bias measurement is estimated from the standard deviation of the residuals of the data points of b(k,z) with respect to b fit (z), where N bins is the total number of bins which contain the individual estimates at the wavenumbers k i within the fitting range. The residuals for the bias measurement in each snapshot are shown in Fig. 15 and the best-fit bias value is shown in Fig. 4 over the entire redshift range.

C Validation of the flat-sky approximation
In Sec. 4 we test a model of the angular power spectrum that includes second-order terms in the bias expansion. In order to evaluate this model and run our MCMC chains efficiently, we neglect wide-angle effects and we assume a uniform redshift distribution of the sources in our tophat bin. In this appendix, we study the impact of these approximations on the angular power spectra of the number counts. For this purpose, we compare the angular power spectra for unbiased tracers in the flat-sky approximation to a) the angular power spectra estimated in full-sky with the code class, using a uniform redshift distribution of the sources within the tophat bins; b) the angular power spectra estimated in full-sky with the code class, for tophat bins and the redshift distribution of the halo population used in our analysis. In this test, we use the same bin width adopted in our analysis in Sec. 4, and nonlinear corrections are included in the power spectrum using the HMCODE recipe implemented in class.
In Fig. 16 we present this comparison at three representative redshifts:z = 0.5 (top row),z = 1 (middle row)z = 2 (bottom row). The left column shows the comparison of the flat-sky approximation with case a), while the right column shows the comparison with case b). We find that wide angle effects are below 4% at z = 0.5, and their impact decreases at higher redshift. On the other hand, the difference between the full computation and our approximation (which also includes the effect of a non-uniform dN/dz in the bin) increases with redshift. At z = 2 it is still below 5% on large scales, and it is still below the statistical errors of our simulated data. Therefore, we can safely use this approximation in our analysis for binsz ≤ 2.