Cosmological Constraints from the Redshift Dependence of the Alcock–Paczynski Effect: Possibility of Estimating the Nonlinear Systematics Using Fast Simulations

The tomographic Alcock–Paczynski (AP) method is so far the best method in separating the AP signal from the redshift space distortion (RSD) effects and deriving powerful constraints on cosmological parameters using the ≲ 40 h − 1 Mpc clustering region. To guarantee that the method can be easily applied to the future large-scale structure surveys, we study the possibility of estimating the systematics of the method using the fast simulation method. The major contribution of the systematics comes from the nonzero redshift evolution of the RSD effects, which is quantified by ξ ˆ Δ s ( μ , z ) in our analysis, and estimated using the BigMultidark exact N-body simulation and approximate COmoving Lagrangian Acceleration (COLA) simulation samples. We find about 5%/10% evolution when comparing the ξ ˆ Δ s ( μ , z ) measured as z = 0.5/z = 1 to the measurements at z = 0. We checked the inaccuracy in the 2pCFs computed using COLA, and find it 5–10 times smaller than the intrinsic systematics of the tomographic AP method, indicating that using COLA to estimate the systematics is good enough. Finally, we test the effect of halo bias, and find ≲1.5% change in ξ ˆ Δ s when varying the halo mass within the range of 2 × 1012–1014 M⊙. We will perform more studies to achieve an accurate and efficient estimation of the systematics in the redshift range of z = 0–1.5.


Introduction
The large-scale structure (LSS) of the universe encodes an enormous amount of information about its expansion and structure growth histories. In the past two decades, large-scale surveys of galaxies have greatly enriched our understanding about the universe (Colless et al. 2003;Beutler et al. 2012;Blake et al. 2011aBlake et al. , 2011bYork et al. 2000;Eisenstein et al. 2005;Percival et al. 2007;Anderson et al. 2012;Alam et al. 2017). The next generation of LSS surveys, such as Dark Energy Spectroscopic Instrument (DESI), 5 Euclid satellite (EUCLID), 6 Large Synoptic Survey Telescope (LSST), 7 and Wide Field Infrared Survey Telescope (WFIRST), 8 will enable us to measure the z1.5 universe to an unprecedented precision, shedding light on the dark energy problem (Weinberg 1989;Riess et al. 1998;Perlmutter et al. 1999;Li et al. 2011;Weinberg et al. 2013).
The Alcock-Paczynski (AP) test (Alcock & Paczynski 1979) is a pure geometric probe of the cosmic expansion history based on the comparison of observed radial and tangential sizes of objects that are known to be isotropic. Under a certain cosmological model, the radial and tangential sizes of some distant objects or structures take the forms of D = D r z c H z ( )  and q D = + D r zD z 1 A ( ) ( ) , where Δz, Δθ are their observed redshift span and angular size, while D A , H are the angular diameter distance and the Hubble parameter computed using theories. In the case that incorrect models were assumed for computing D A and H, the values of Δr P and Δr ⊥ are wrongly estimated, resulting in geometric distortions along the line-ofsight (LOS) and tangential directions. This distortion can be quantified via statistics of the large-scale galaxy distribution, and has been widely used in galaxy surveys to place constraints on the cosmological parameters (Ryden 1995;Ballinger et al. 1996;Matsubara & Suto 1996;Outram et al. 2004;Blake et al. 2011a;Lavaux & Wandelt 2012;Alam et al. 2017;Mao et al. 2017;Ramanah et al. 2019).
The tomographic AP (Li et al. 2014(Li et al. , 2015 method is a novel application of the original AP test that uses the redshift evolution of the clustering anisotropies, which are sensitive to the AP effect while being relatively insensitive to the anisotropies produced by redshift space distortions (RSDs). This makes it possible to differentiate the AP distortion from the large contamination from the RSD effect. Li et al. (2016) first applied the method to the Sloan Digital Sky Survey (SDSS) Baryon Oscillation Spectroscopic Survey (BOSS) DR12 galaxies, and achieved∼30%-40% improvements in the constraints on the ratio of dark matter Ω m and dark energy equation of state (EOS) w when combining the method with the data sets of the Planck measurements of the cosmic microwave background (CMB; Planck Collaboration et al. 2016), type Ia supernovae (SNIa; Betoule et al. 2014), baryon acoustic oscillations (BAOs; Anderson et al. 2014), and the local H 0 measurement (Riess et al. 2011;Efstathiou 2014). In follow-up studies, Li et al. (2018) and Zhang et al. (2019) studied the constraints on models with the possible time evolution of the dark energy EOS, and showed that the method can decrease the errors of the parameters by 40%-50%. Li et al. (2019) forecast the performance of the tomographic AP method on future DESI data and found that when combining with Planck+BAO data sets the dynamical dark energy constraints are improved by a factor of 10.
The tomographic AP method is so far the best method in separating the AP signal from the RSD effects and extracting information from the - h 40 Mpc 1 clustering regions. Thus, it is essentially important to conduct more studies to guarantee that the method can be applied to future LSS surveys. The biggest caveat of the method is the systematics from the time evolution of the RSDs. On the scales explored by the tomographic AP method, the accurate quantification of RSDs can only be achieved via precise numerical simulations. Li et al. (2016) utilized the Horizon Run 4 (HR4) N-body simulation  to estimate the systematics when applying the method to the SDSS galaxies. They showed that, within the redshift range of 0.2-0.5/0.2-0.7, the systematics creates2%/6% time evolution in the anisotropic correlation function x m Dŝ ( ). Li et al. (2018) conducted more tests on the systematics and found that the shifts in the derived cosmological parameters caused by the systematics are well below the 1σstatistical error. Park et al. (2019) studied the systematics using simulations generated in five different cosmologies, and reported a nonnegligible cosmology dependence. More studies about the systematics have been performed in Li et al. (2015Li et al. ( , 2019, and their conclusions are consistent with Li et al. (2016Li et al. ( , 2018 and Park et al. (2019).
While many studies on the systematics of the method have been performed, it is necessary to conduct more checks regarding the systematics including: 1. First, the redshift coverage of current studies of the method are limited to z0.7. To meet the requirements of future surveys, this should be enlarged to z∼1.5. 2. Second, so far the systematics are mainly studied by using the Horizon Run N-body simulations (Kim et al. 2011, which are performed in the Ω m =0.26 ΛCDM cosmology. In different cosmologies, the influence of the RSD effect is also different. The cosmological dependence of the systematics is first studied in Park et al. (2019), and remains to be investigated in greater detail. 3. Finally, it would be necessary to have a fast method for the systematics estimation. The next-generation experiment will survey an unprecedented large volume of the universe. A fast and accurate method would be very helpful if we were to estimate the systematics of the tomographic AP method for these surveys.
In this paper, we explore the first and third issues listed above. In Section 2, we introduce the simulation materials used in the paper, including both the N-body and fast simulation samples. The methodology of the tomographic AP method is briefly reviewed in Section 3. In Section 4 we present our results, including the high-z measurements of the systematics, a comparison between N-body and fast simulation results, a test on the clustering scales, and a check on the halo bias effect. We conclude in Section 5.

Simulation
We use two sets of simulation samples, the BigMultiDark (BigMD) N-body simulation sample (Klypin et al. 2016) and also a set of fast simulation samples generated using the COmoving Lagrangian Acceleration (COLA) algorithm.
In both sets of samples, to mimic the RSDs caused by galaxy peculiar velocities, we perturb the positions of halos along the radial direction, using the following formula: where v LOS is the LOS component of the peculiar velocity of galaxies and z is the cosmological redshift. For testing the redshift evolution of the RSD effect, in both samples we use the outputs of halos at 12 redshifts, i.e., z={0, 0.102, 0.2013, 0.2947, 0.4037, 0.5053, 0.6069, 0.7053, 0.7976, 0.8868, 1, 1.445}, respectively.

The BigMD Simulation
The Multiverse simulations are a set of cosmological N-body simulations designed to study the effects of cosmological parameters on the clustering and evolution of cosmic structures. Among them, the BigMultidark simulation is produced using 3840 3 particles in a volume of h 2.5 Gpc  (Klypin et al. 2016). The size and number of particles of this simulation are sufficient for the purpose of the study in this work (although it is smaller than the Horizon Run N-body simulations, Kim et al. 2009, used in Li et al. 2014, 2016. Having both a large volume and a good resolution, this simulation is able to accurately reproduce the observational statistics of the current redshift surveys (Rodríguez-Torres et al. 2016). We use the publicly available halo catalog of BigMD simulation, created using the ROCKSTAR halo finder (Behroozi et al. 2013). ROCKSTAR is a halo finder based on adaptive hierarchical refinement of a friends-of-friends groups in six phase-space dimensions and one time dimension, allowing for robust tracking of substructure. To make the samples at different redshifts comparable to each other, we maintain a constant number density = n 0.001 ) in all snapshots. Both halos and subhalos are included in the analysis. In the Appendix, a rough comparison between the 2pCFs of the BigMD ROCKSTAR halos and the BOSS CMASS galaxies shows that the simulation sample can recover the observational results at a5% accuracy level.

The COLA Samples
Although powerful in constructing high-quality, realistic dark matter halo catalogs, N-body simulations are computationally expensive (Angulo et al. 2012;Fosalba et al. 2015;Heitmann et al. 2015;Potter et al. 2017), making them difficult to use for creating a large number of mock catalogs for current and future large surveys such as SDSS, LSST, EUCLID, and DESI. In order to circumvent this problem, some fast algorithms have been proposed to reproduce the large-scale statistics of N-body simulations. An incomplete list of these algorithms includes PTHalos (Manera et al. 2013(Manera et al. , 2015, In this work, we concentrate on the possibility of using the COLA algorithm (Tassev et al. 2013) as a replacement for the N-body method to quickly generate a large number of mocks for the estimation of systematics in the tomographic AP method. The second-order Lagrangian perturbation Theory (2LPT) adopted by the COLA code is fast in computation and still accurate enough in describing the large-scale dynamics. Due to the fact that 2LPT can be easily incorporated in any N-body code, COLA combines it within N-body simulations by adopting the 2LPT for time integration for large-scale dynamical evolution, and using a full-blown N-body code with the Particle-Mesh (PM) algorithm to deal with small-scale dynamics. Compared with the fastest simulation algorithms, COLA is better in simulating the structures on nonlinear scales (Chuang et al. 2015b), which is rather suitable for the science case in this work.
We generate 150 COLA simulations in the BigMD cosmology as reference samples to the BigMD simulation. For each COLA simulation, 600 3 particles in a box size of 512h Mpc 1 are used and the time steps are set at 20 ( Table 1). The resolution of the COLA simulations is 30% higher than the BigMD simulation, and the total volume of the 150 simulations, being » h 2720 Mpc 1 3 ( ), is also slightly larger than the BigMD box size. Besides, to be comparable with the BigMD simulation samples, we still use ROCKSTAR halo finder to build up the halo catalogs from the COLA particles.

Methodology
In what follows, we first briefly introduce the AP effect and its redshift evolution, and then discuss how to use it for estimating the cosmological parameter, as well as the necessity of conducting studies about its systematics.

The AP Effect
The AP effect (Alcock & Paczynski 1979) is known as geometric distortions when incorrect cosmological models are assumed for transforming redshift to comoving distance (see Li et al. 2014 for a detailed description). We probe the size of an object in the universe through measuring the redshift span Δz and angular size Δθ, which are represented as where the cosmological dependence enters via the Hubble parameter H and the angular diameter distance D A . In the special case of a flat universe with a constant dark energy EOS, they take the forms of where we have neglected the contribution from the radiation. Here = + a z 1 1 is the cosmic scale factor, H 0 is the present value of Hubble parameter, and r(z) is the comoving distance.
When we adopted the wrong cosmological parameters, the Dr  and Δr ⊥ in Equation (2) would take the wrong values, resulting in the distorted shape (AP effect) and the incorrectly estimated volume (volume effect). In the directions parallel and perpendicular to the LOS, the distortions are where "true" and "wrong" denote the values of quantities in the true and incorrectly assumed cosmologies, respectively. As a result, apparently the shape of the object is changed by a ratio of The above relationships mean that the AP and volume effect, once detected at any clustering scales, would lead to constraints on the cosmological parameters that control the cosmic expansion.

2pCF and Anisotropy
We use the 2pCF statistics to measure the anisotropies in the LSS clustering. The Landy-Szalay estimator (Landy & Szalay 1993) is adopted to calculate the 2pCF, where DD is the number of galaxy-galaxy pairs, DR is the number of galaxy-random pairs, and RR is the number of random-random pairs. All those numbers are measured with a dependence of (s, μ), where s is the distance between the galaxy pair and μ=cos(θ), with θ being the angle between the line joining the pair of galaxies and the LOS direction. We use the publicly available code CUTE (Alonso 2012) for computation of the 2pCF. In order to probe anisotropy, as was done in Li et al. (2015, 2016), we integrate the 2pCF over s and only focus on due to reasons explained in Li et al. (2015Li et al. ( , 2016. To remove the uncertainty from the clustering strength and galaxy bias, we further normalize the amplitude of x m Ds ( ) to only study its shape, i.e., A cut μ<μ max is imposed to reduce the finger-of-god (FOG) effect (Jackson 1972; and also the fiber collision effect when studying the observational data).

The Redshift Evolution
In addition to the AP effect, another source of apparent anisotropy in the galaxy distribution we observed is the RSD effect due to the peculiar velocity of galaxies, resulting in significant anisotropy even if the adopted cosmology is correct. Li et al. (2014) found that the anisotropy generated by the RSD effect is, although large, maintaining a nearly constant magnitude within a large range of redshift, while the anisotropies generated by AP varies with the redshift much more significantly. So they proposed to measure the AP effect using the redshift dependence of the distortion.
Due to the growth of structure, the galaxy peculiar velocities evolve with redshift, and thus the RSDs are not exactly constant with time. The small redshift evolution of RSD would cause redshift dependence in the LSS anisotropy, which is the main source of the systematics of this method.
In this work, we measure dx m D z z , , s i ĵ ( )from the BigMD and COLA simulations to quantify the systematics from the RSD effect. The redshift evolution of the clustering anisotropy can be described as where z i and z j are two different redshifts. The above quantity then represents the systematics when the measurement is done in the correct cosmology (the simulation cosmology). We further use to quantify the magnitude of the systematics contributed at a given redshift.

Results
In this section, we present the redshift evolution of x m Dŝ ( ) measured from the BigMD and COLA simulations at a redshift region of 0-1.445. A comparison was made between the results measured from the two sets of simulations. Figure 1 shows the 2D contour of the 2pCF ξ(s, μ), from the BigMD simulation and COLA simulation and at redshifts 0 and Figure 1. The 2D counter of 2pCF ξ(s, μ) from the BigMD and COLA samples, at redshifts 0 and 1. The contour lines are not horizontal due to the anisotropy produced by the RSD effect. The tilts on the left side and the right side are caused by the finger-of-god (FOG) and Kaiser effects, respectively. The slight difference, which is below a few percent between the BigMD and COLA results, indicates that the COLA has the ability to reproduce the clustering pattern at

Results of ξ(s, μ)
. Noted that the similarity in the region of 1−μ0.2 implies that COLA can correctly simulate the FOG effect. 1, respectively. Due to the RSD effect, the contour lines are not horizontal. The tilts at the left and right sides are caused by the FOG (Jackson 1972) and Kaiser effects (Kaiser 1987), respectively.
The similarity between the results from the BigMD and COLA simulations implies that COLA as an approximate algorithm can successfully reproduce the clustering pattern at The difference between the two sets of results is below a few percent. In particular, the similarity in the region of 1−μ0.2 suggests that COLA has the ability to correctly simulate the FOG effect.

x m ( ) in the Correct Cosmology
In Figure 2, we compare the x m D s ŝ ( ) measured from the COLA and BigMD samples, constructed using the underlying together with an amplitude normalization. We split the angular range of m -  0.06 1 1 into as many as 40 bins. The upper left panel shows x Dŝ measured at z=0, 0.102, 0.5053, 1.0, 1.445. Due to the RSD effect as well as its redshift evolution, all curves have≈30% tilt and are more tilted at higher redshift. In the upper right panel we find the results from COLA and BigMD merely have a1% difference at z<1, and1.5% at z>1. In the middle panels we plot the evolution against z=0 to better quantify the redshift evolution of the clustering anisotropy. In the region of small/ large 1−μ, the evolution becomes as large as 5% at z=0.5, and increases to 10% at z=1, 1.445. COLA can well reproduce the BigMD results at levels of0.5%, except in the case of z=1, where we find a 2%-3% discrepancy. In the lower panels, we plot x m D d dz ŝ ( ) at the three redshifts of z=0.051, 0.4545, 0.9434 to measure the redshift evolution from a different perspective. The estimation from COLA has5% error at z<0.5, and has a relatively large error of 15% at z∼1. true cosmology (i.e., the simulation cosmology). We compute them by integrating ξ(s, μ) within the range )caused by the redshift evolution of the RSD effect. In the region of small/large 1−μ, i.e., the region parallel/perpendicular to the LOS, the evolution becomes as large as 5/2.5% at z=0.5, and increases to 10/5% at z=1, 1.445. COLA can well reproduce the BigMD results at levels of0.5%, except in the case of z=1, where we find a 2%-3% discrepancy.
3. In the lower panels we plot x m D d dz ŝ ( ) at the three redshifts of z=0.051, 0.4545, 0.9434 to measure the redshift evolution from a different angle. The estimation from COLA has5% error at z<0.5, and has a relatively large error of 15% at z∼1.

x m ( ) in the Wrong Cosmologies
In Figure 3, we plot the measured x D Dŝ and x D d dz ŝ when assuming the "true" cosmology as well as a wrong cosmology Ω m =0.25, w=−1.0. In the wrong cosmology, the amount of difference between COLA and BigMD results is similar to that which was measured in the correct cosmology, while this amount of difference is smaller than the shift of x Dŝ and x D d dz ŝ induced by the wrong cosmology. Clearly, COLA is accurate enough for the purpose of systematics estimation in distinguishing this cosmology with the underlying true cosmology.

Different Clustering Scales
In previous studies of Li et al. (2015Li et al. ( , 2016Li et al. ( , 2018Li et al. ( , 2019 the authors only focus on the clustering scales of  ) would not enhance the power of the statistics, but rather decrease the signal-to-noise ratio and hence weaken the power of constraints. Second, on large clustering scales, one can directly compare the theoretical expectation and the measured results of the 2pCFs to achieve a more complete exploration of the information (not just the AP signal) encoded in the data. The tomographic AP method is designed to explore the nonlinear clustering regime, which is difficult for the traditional approaches.

Bias Effect
The effect of halo bias on x Dŝ has been briefly investigated in Li et al. (2016). Here we revisit this issue by using finer mass bins covering a larger range. Figure 5 shows the x m Dŝ ( ) at z=0, measured in seven mass bins ranging from 2×10 12 M e to 10 14 M e . We find that the difference among them is not apparent, mostly below 0.5%.  , respectively. We find that using large clustering scales results in a larger slope and evolution of slope in ξ Δs (μ). The statistical fluctuation (scattering) increases with increasing clustering scales.

Conclusion
The tomographic AP test is a novel statistical method entering nonlinear clustering scales of h 6 40 Mpc 1 -. Tight cosmological constraints have been achieved by applying the method to the SDSS data (Li et al. 2016(Li et al. , 2018Zhang et al. 2019). The method measures the redshift evolution of anisotropy, usually quantified by x m Dŝ ( ), to mitigate the effects from the RSDs while still being sensitive to the AP signal. In this work, we studied the systematics of this methods in detail, as an early step in preparation of its application to future LSS surveys.
Since the next-generation galaxy surveys will reach much higher redshift than the current ones, we study the redshift evolution of x m Dŝ ( ) to the z=1 region. We find persistent redshift evolution in the whole region, which manifests itself as a persistent, nonzero ) to quantify the difference between the high-redshift and z=0 results, we find an evolution of≈5% at z=0.5. The evolution reaches≈10% at z=1.
We also study the possibility of estimating the systematics via the cheap mocks generated by COLA. In most cases we studied, COLA can well reproduces the N-body results. In terms of estimating x m COLA can reach an accuracy of 0.5% at z=0, 0.5, while a 2%-3% discrepancy is found at z=1, 1.445. Finally, we studied the dependence of the performance on the clustering scales and galaxy bias.
Our investigation suggests that the application of the method should be easy for z0.5, where the systematics are relatively small, and also can be accurately estimated using fast simulations.
In the next few years, surveys such as DESI and Taipan 9 will probe the nearby universe with unprecedented precision. Their data sets would be ideal for the application of our method.
Applying this method to the region of z≈1 would be challenging. The systematics reach ≈10%, while the inaccuracy of COLA also reaches 2%-3%. To be well prepared for the cosmological analysis on the DESI, EULCID, WFIRST, and PFS surveys, this problem has to be resolved. To achieve a fast and accurate systematics estimation, one can improve the performance of the simulation by increasing the time step, setting high initial redshift, trying alternative algorithms, or making a correction to the COLA results using N-body results as baseline reference to decrease its inaccuracy.
This work shows that COLA is a possible tool for estimating systematics in multi-cosmologies. A fast and accurate enough systematics estimation method is important if one wants to conduct systematics correction in multiple cosmologies. In all previous works (Li et al. 2015(Li et al. , 2016(Li et al. , 2018, the study of systematics was done using a simulation running under only one set of cosmological parameters. The cosmological dependence of the systematics can be efficiently studied using the COLA algorithm. A possible way to further reduce systematics is to design better statistics that is more insensitive to the redshift distortions. For example, we can use the marked correlation function (Beisbart & Kerscher 2000;Gottloeber et al. 2002;Sheth & Tormen 2004;Sheth et al. 2005;White 2016;Satpathy et al. 2019), which assigns different weights to regions with different densities, to suppress the strong RSD coming from the most clustered regions. These issues have not yet been discussed in details in this paper and are worth further investigating in future research.
In general, the tomographic method requires precise redshift measurement from spectroscopic surveys. The method uses the clustering region at »h 5 Mpc 1 , requiring a redshift measurement with an accuracy of δz≈0.001(1 + z), which is difficult for most photometric surveys (e.g., DES and LSST).
The tomographic AP method is one of the best methods in deriving cosmological constraints using the - h 40 Mpc 1 clustering region. We will continue to work on this method so that we can safely use it to derive tight, robust cosmological constraints from the next-generation LSS surveys. Figure 6 illustrates a comparison between the 2pCFs measured from (1) two sets of BigMD ROCKSTAR halo samples, at the z=0.48 and 0.62 snapshots; (2) two sets of COLA ROCKSTAR halo samples, at the same redshifts; (3) two sets of subsamples selected from the BOSS DR12 CMASS galaxies, by imposing the redshift ranges of 0.430<z<0.511 and 0.572<z 6 <0.693; these two subsamples have effective redshifts of 0.48 and 0.62, respectively.
In the left panel we plot x m D s ŝ ( ) . The simulation samples achieve5% level of accuracy in recovering the observational measurements. The only exception is the leftmost point at 1−μ=0. Here (1) the FOG effect is very strong; (2) the observational result may be problematic due to the fiber collision effect. Anyway, this region is always abandoned when conducting cosmological analysis (Li et al. 2016 imposed a cut 1 − μ > 0.03).
Notice that the strong peak near 1−μ=0, produced by the FOG, is not detected in Figures 2-5, as we imposed a cut 1−μ<0.03 in these figures.
The redshift evolution of x m D s ŝ ( ) is plotted in the right panel. The observational results have much larger statistical error. A linear regression shows that results from the BigMD, COLA, and BOSS samples are in good consistency with each other.
In all plots, the BigMD and COLA results (the dashed and dotted curves) are so close to each other that we can hardly distinguish them by eye.