New leaves of the tree: percolation analysis for cosmic web with discrete points
Jiajun Zhang Dalong Cheng Ming-Chung Chu



School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai, China
Department of Physics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong

School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai, China

Department of Physics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong


Percolation analysis has long been used to quantify the connectivity of the cosmic web. Most of the previous work is based on density fields on grids. By smoothing into fields, we lose information about galaxy properties like shape or luminosity. The lack of mathematical modelling also limits our understanding for the percolation analysis. In order to overcome these difficulties, we have studied percolation analysis based on discrete points. Using a Friends-of-Friends (FoF) algorithm, we generate the S-bb relation, between the fractional mass of the largest connected group (S) and the FoF linking length (bb). We propose a new model, the Probability Cloud Cluster Expansion Theory (PCCET) to relate the S-bb relation with correlation functions. We show that the S-bb relation reflects a combination of all orders of correlation functions. Using N-body simulation, we find that the S-bb relation is robust against redshift distortion and incompleteness in observation. From the Bolshoi simulation, with Halo Abundance Matching (HAM), we have generated a mock galaxy catalogue. Good matching of the projected two-point correlation function with observation is confirmed. However, comparing the mock catalogue with the latest galaxy catalogue from SDSS DR12, we have found significant differences in their S-bb relations. This indicates that the mock galaxy catalogue cannot accurately retain higher order correlation functions than the two-point correlation function, which reveals the limit of the HAM method. As a new measurement, S-bb relation is applicable to a wide range of data types, fast to compute, robust against redshift distortion and incompleteness, and it contains information of all orders of correlation functions.   



New leaves of the tree: percolation analysis for cosmic web with discrete points
    Ming-Chung CHU
    December 30, 2023
================================================================================



§ INTRODUCTION


Cellular or clumpy? That is the question. In the last century, when the first generation of redshift surveys were available, people were confused by the filamentary structure of galaxy distribution <cit.>. The existence of the cosmic web was qualitatively confirmed, but a quantitative description of it was not available until Zeldovich's percolation analysis <cit.>. We classify the percolation analyse for the cosmic web into two branches. The first one, dealing with discrete points, is called continuum percolation, whereas the second one, dealing with the density field on grids, is called site percolation. Continuum percolation analysis was also called cluster analysis in 1980s, because the same technique was used to define galaxy clusters. Cluster analysis can clearly tell the difference between a clumpy structure and web-like pattern <cit.>. However, Dekel & West <cit.> concluded that cluster analysis was insensitive to cosmological parameters, and it was also non-trivially affected by volume and point density. Thereafter, the site percolation analysis became popular, and properties of site percolation in large scale density fields have been studied in detail <cit.> . By smoothing points into field, the technique has also been used to look for properties of galaxy distribution <cit.>. We would like to point out that direct comparison between the percolation properties of density field and galaxy distribution is not fair, because galaxies don't trace the density field perfectly. If we can generate a mock catalogue of galaxies and compare it with a real redshift survey, the comparison will be fair and the result will give us some useful insights in properties of the cosmic web. 

On the other hand, the establishment of the ΛCDM model and rapid development of N-body simulation techniques provide us a general picture of the formation of the cosmic web <cit.>. We are now confident about the success of the ΛCDM model in large scales. However the three-point, four-point or even N-point correlation functions should also be tested <cit.>, though they are all very difficult to calculate and  model in theory. In principle, the two-point correlation function costs O(N^2) calculations, and the three-point correlation function costs O(N^3) calculations, so on and so forth. Although we have better ways to optimize the calculation, higher order correlation functions  still require unbearable calculation resources. Thus, several methods to characterize the cosmic web have been proposed, such as Minkowski Functional, Minimum Spanning Tree, and Genus Analysis <cit.>. These methods focus on the morphology and topology of the cosmic web. 

The comparison between simulation and observation is not straight forward. We are still lacking of knowledge about galaxy formation and evolution, and matching galaxies with simulated dark halos heavily depends on model <cit.>. Surprisingly, Halo Abundance Matching (HAM), a newly developed technique based on simple arguments, generates mock galaxy catalogues successfully <cit.>. HAM follows the  halo mass, halo circular velocity, halo potential well, galaxy circular velocity, galaxy luminosity, galaxy stellar mass, independent of galaxy formation model, and HAM is easy to understand.     

Shandarin et al. <cit.> developed and studied the percolation analysis, treating voids and superclusters on an equal footing and finding the nature and nurture of the formation of cosmic web. Using a similar method, the SDSS galaxy sample has been tested with the percolation analysis <cit.>. However we are still far from understanding the percolation transition in
the cosmic web. Density estimation is a prior step for site percolation analysis, but estimating  the density field from galaxy observation is not easy. We should take into account galaxies' different properties in the percolation analysis.  

Based on these considerations, we used the continuum percolation analysis to calculate the S-bb relation, where bb stands for the FoF linking length and S stands for the mass fraction of the largest group for the given linking length bb <cit.>. It is quite similar to the largest cluster statistics in previous percolation analyses <cit.>. However, we find that the continuum percolation analysis can distinguish among cosmological models with high resolution simulations, whose resolution was not available in 1985. We propose the Probability Could Cluster Expansion Theory (PCCET), to calculate the S-bb relation theoretically. We measure the S-bb relation in our mock galaxy catalogue and compare it with a real galaxy catalogue. We find that HAM cannot reproduce the S-bb relation as well as the projected two-point correlation function. 



§ METHODOLOGY



 §.§ S-bb Relation and Percolation Transition

For any discrete point distribution, we can always apply the FoF algorithm to find a group of points. The points can be massive points in N-body simulations, halos in simulated halo catalogues or galaxy catalogues in redshift space. We define S_1 and S_2 as:


    S_1(bb)=N_largest/N_total,


    S_2(bb)=N_second-largest/N_total,

 where
 N_total is the total number of points,
 N_largest (N_second-largest) is the number of points in the largest (second largest) group found by FoF, and normally, we refer to the S-bb relation the function of S_1(bb).
 It is obvious that S_1(bb) reveals quite similar information as the cluster analysis and site percolation analysis.
 Our simulation tests are all done using  <cit.>. Fig. <ref> is the S-bb relation for a random distribution as a comparison for later figures. It is not surprising that all three sets of data corresponding to different box sizes and number of particles give the same S_1(bb) and the transition threshold is around bb=0.88.



We define the transition threshold as the bb where S_2(bb) gets the maximum value. Another definition of the transition threshold was also proposed in <cit.>; in our case, it can be defined as the bb where the mass-weighted sum of the number of connected structures, μ^2=Σ_n N_nN_groups, where N_n is the number of particles in the n^th group and N_groups is the total number of connected structures, gets the maximum value. 

We measured the S_1, S_2 and μ^2 using the halo catalogue of the MDR1 simulation<cit.> with haloes with mass larger than 10^12M_⊙. The MDR1 simulation was performed in a box with 1h^-1Gpc sides, and so we separated the sample into 64 subsamples of boxes with 250h^-1Mpc sides to estimate the cosmic variance, shown as the errorbars in Fig. <ref>. The transition threshold bb_c is 0.74 and the power index α given by S_1=(bb-bb_c)^α was measured to be α=0.46, which is consistent with <cit.>. 
 
As shown in Fig. <ref>, both definitions of the transition threshold are good, and their FWHM are similar. We will discuss the meaning of the transition threshold in the definition of S_2 later for easier understanding. When the percolation transition happens, the largest structure grows tremendously as bb gets larger, and the second largest structure reaches its most massive point. After that, the mass of the second largest structure becomes smaller, not because the structure shrinks, but because the original second largest structure is connected to the largest one while another structure is recognized as the new second largest one. This connection and replacement reveals the nature of percolation: the largest structure becomes the one and the only one dominant structure, connecting most of the high density regions together and leaving only the far apart particles in voids alone.  
Thus, the definition of transition threshold based on either S_2 or μ^2 is appropriate.



While the ΛCDM (cold dark matter) model is not clearly distinguishable from the ΛWDM (warm dark matter) model, the ΛHDM (hot dark matter) model is clearly ruled out by observation of large scale structures. We would like to test our percolation analysis with these three models. The S-bb relations are expected to be indistinguishable for ΛCDM and ΛWDM models, but those for ΛHDM and ΛCDM models should be clearly different. This test is the bench mark test of percolation analysis.   We generate the initial conditions with  <cit.> and run the simulation from redshift 49 to redshift zero using the N-body simulation code <cit.>. The transfer function and additional thermal velocity for HDM and WDM are given by <cit.>. These three simulations share the same realization. Some basic parameters for the simulations are given in Table. <ref>. 


As shown in Fig. <ref>, the HDM model is quite different from CDM and WDM models in their matter 
power spectra and S-bb relations. The matter power spectrum is calculated by the publicly available
code <cit.>. It is well accepted that the HDM model has been ruled out by
measuring the power spectrum in both CMB <cit.> and large scale structures <cit.>.
We have shown here that the difference between HDM model and CDM model is significant enough
in the S-bb relation to be used as a supplementary. On the other hand, the unique
advantage of using the S-bb relation will be discussed in the following sections. 



 §.§ Probability Cloud Cluster Expansion Theory

The percolation analysis has been used for a long time to study the cosmic web. For site percolation there are analogies to similar questions in many other fields <cit.>. But for discrete points, there is still no available theory to link percolation with basic cosmological parameters. Here we propose the Probability Cloud Cluster Expansion Theory (PCCET) to connect the S-bb relation with correlation functions. The idea is as follows:



  * Consider a box with only one particle and take it as the only member in the largest group. So the S-bb relation will be a constant line S=1.

  * Add in another particle randomly, the probability of linking these two particles in a single group P(bb) being monotonically related to the linking length bb. The S-bb relation will be a step function jumping from S_1=0.5 to S_1=1.0 at a certain transition value of P: S=1/2+1/2P.

  * Add in the third particle randomly, the probability of linking any two of them being still P and the probability of linking all three of them shall be P^2. While from a certain particle as the starting point of the group, the probability of linking the other two particles shall be 2P but deducting P^2 to avoid over counting. The S-bb relation will then be a three-step function, with steps at 1/3, 2/3, and 1: S=1/3+1/3P+1/3P(2P-P^2).

  * Continuing on, we consider more and more particles, taking care of over counting and over deducting, and we finally find the relation S(P) when the number of particles goes to infinity. We hide all information about bb in the function P(bb).  

So in general, we can write S(P) as


    a_1=1/N, 
    
    a_n=a_n-1b_n-1 (n>1), 
    
    b_n=1-(1-P)^n, 
    
    S=∑_n=1^N a_n.


This series expansion is the reason why we call it cluster expansion theory. Cluster expansion is often used in condensed matter physics. It can successfully account for phase transitions in some typical problems. We borrow the idea of cluster expansion because of the similarity between percolation transition and phase transition. We have shown that S(P) is monotonically increasing with P and has a limit between 0 and 1. When N→∞, S(P) will become the Euler Function. The proof in detail is given in Sec. <ref>. 

In the framework of this theory, a particle cannot be treated as a discrete point. In fact, it is better to take the particle as a representation of the probability cloud. For a random distribution, the probability for a particle to be closer to another particle than distance r is inversely proportional to the volume of a sphere with radius r. However, we need a volume to set as the standard for this probability. Here we take all the particles in the theory as test points. We can first put in the test center. We imagine that every test center has its own cube with periodic boundary conditions. The length of each side of this test cube can be set as the average distance between particles in a real simulation. Then the distance to the center can be simply understood as the linking length bb. If bb is larger than 1/2, the influence of this test center will reach outside the test box. In such a case, we shall remove the influence region outside the test box and only consider the volume inside the box as the real probability distribution region. Shown in Eq. <ref> is the P-bb relation for a random distribution.






From the PCCET, we can find the relation between percolation and correlation function. Because of the cluster expansion, the 2-point correlation function is the correction to P as the probability of linking up 2 points, while the 3-point correlation function is the correction to P^2 as the probability of linking up 2 pairs of points, etc. In general, the n-point correlation function is the correction to P^(n-1) as the probability of linking up n-1 pairs of points. This shows that percolation actually contains the information of all orders of correlation functions. Our probability cloud cluster expansion theory also shows how every order of correlation function contributes to the S-bb relation, which represents the percolation phenomenon in the cosmic web.

Unfortunately, in the non-linear regime, high order correlation functions are also quite large and not necessarily smaller than the 2-point correlation function. So based on our theory, it is still not possible to calculate the whole S-bb relation analytically. It is important to use N-body simulations to set up constraints for the S-bb relation. After all, PCCET is still the only model explaining percolation. It can be used to find the S-bb relation in the linear-regime and it is also possible that in the future, we may find a way to use it in the non-linear regime to extract high order correlation functions.




 §.§ Halo Abundance Matching

From PCCET, we can see that the S-bb relation contains information of all orders of correlation functions. Thus we would like to see whether it provides new information in comparing simulation with observation. We first need to generate a mock galaxy catalogue and check the projected 2-point correlation function to see whether it matches the real galaxy catalogue. Then we will compare the S-bb relations of the mock and real galaxy catalogues. 


Based on the Bolshoi simulation <cit.> BDM halo merger tree, we follow the work in <cit.> using Halo Abundance Matching to generate a mock galaxy catalogue. Here is the process:


  * We take the Bolshoi simulation BDM halo catalogue at redshift 0 to extract the final positions and peculiar velocities of the galaxies. We assume that every BDM halo, including the host haloes and subhaloes, can host at most one galaxy at its center. 

  * From the merger tree of BDM haloes in the Bolshoi simulation, we find the maximum circular velocity V_max in the merger history with the halo circular velocity defined as the maximum value of the circular orbital velocity according to the density distribution of the halo.

  * We adopt the luminosity function from <cit.>, which is of the Schechter form. The luminosity function is:



    Φ(M)dM=0.4log_10Φ_*10^-0.4(M-M_*)(α+1)
    ×exp(-10^-0.4(M-M_*))dM,


 
where M is the absolute magnitude and Φ_*=0.0078, M_*-5log_10h=-20.83, α=-1.24, given by <cit.> using SDSS6 data.

  * We sort the haloes in the Bolshoi simulation according to their V_max. Record every halo's ranking. Thus, we can easily obtain the number density of haloes with V_max larger than a certain value by the ranking.

  * The integral of the luminosity function is an incomplete gamma function. We apply an one-one mapping from the number density of haloes with V_max larger than a certain value V to the integral of the luminosity function. Thus we obtain an one-one mapping between V and absolute magnitude M.

  * Combining the data together, we get the halo catalogue with the galaxy absolute magnitude. This is also the mock galaxy catalogue with scale 250h^-1Mpc. In this catalogue, we have positions, peculiar velocities and luminosities for the mock galaxies at redshift 0.



Halo Abundance Matching has been proven quite successful in recreating a galaxy catalogue which can reproduce the observed projected galaxy 2-point correlation function in the bins of absolute magnitude close to the steep transform range in the luminosity function, around Φ=Φ_*. Many works follow the idea of HAM <cit.>. There are slightly different methods to obtain a mock catalogue based on the same basic idea of HAM. Comparing to <cit.>, we did not take stochastic scattering into account. Since there are obvious bias to dimmer galaxies in stochastic scattering, we avoid this process. On the other hand, we take maximum circular velocity in the merger history of every halo as its monotonic mapping label. The physical insight of this picture is clear. Gas should fall into the center of a halo much faster than the merger rate, and once the gas has fallen into the halo center, it is so dense that the merger process cannot do much to the central dense region. So the galaxy luminosity depends more on the steepest potential well in the merger history, denoted by the maximum circular velocity V_max.


We follow the paper <cit.> to calculate the projected 2-point correlation function. Our results are shown in Fig. <ref>. We have confirmed the previous results given by <cit.> that the HAM reproduced sample exhibits the same projected 2-point correlation functions as the observation in the luminosity range [-20,-19], [-21,-20] and [-22,-21], with differences mostly within 1 or 2 σ. We take these samples for further percolation analysis comparison. We shall also notice that the bias of the projected 2-point correlation function depends on the luminosity of the galaxy. We shall have a different view of the bias in the S-bb relation.

 


§ RESULT



 §.§ S-bb Relation Evolution

We are interested in the evolution of the cosmic web, which we study through an ΛCDM cosmological simulation with 512^3 particles in a box with 1h^-1Gpc sides, initial conditions generated by BBKS power spectrum <cit.> and Zel'dovich approximation <cit.>. We apply Ω_m=0.28, Ω_Λ=0.72 and h=0.7 for this test. As shown in Fig. <ref>, the S-bb relation at high redshifts looks similar to that of random distribution (Fig. <ref>), while as the density contrast increases with time, the S-bb relation at low redshift becomes rather different, which has already been known in many previous works <cit.>. 






 §.§ Bias Determination



The error bars of the S-bb relation shown in fig. <ref> are determined as follows. Based on Bolshoi simulation<cit.>, we have only one realization of the HAM mock sample. So we use jackknife re-sampling to estimate the error bars. Here are the steps:


  * randomly choose 90% of the HAM mock sample for ten times, generating 10 different mock samples.

  * Measure the S-bb relations of these 10 mock samples.

  * Take the mean and standard variance of these 10 S-bb relations as the value and error bar.


From Fig. <ref> we can see that the lower luminosity samples' transition threshold bb is lower than those of the high luminosity samples. Based on this result, we may guess that the lower luminosity samples have higher bias in the 2-point correlation function than higher luminosity samples. However it is well known that more luminous galaxies have higher bias which seems to be contradictory to our results here. Note that the definition of bb is such that it is normalized to the average distance of points. So we shall look at the 2-point correlation function with the projected distances r_p also normalized to the average distance, which is shown in Fig. <ref>. The higher the luminosity of the samples are, 
the lower their normalized 2-point correlation function is. Therefore, a smaller transition threshold is well
expected.



Galaxies' real spatial distribution is different from their redshift space distribution due to the peculiar velocities of the galaxies. Using the HAM mock sample data with peculiar velocity, we can study the redshift distortion effect in the S-bb relation. We uniformly choose 14 different directions to look at the sample, all of which are about z=0.1 away from the sample center. The distance is translated by the cosmological parameters used in the Bolshoi simulation. As shown in Fig. <ref>, we find that redshift distortion has little effects on the S-bb relation. Because the S-bb relation shows mainly the topological information of the cosmic web, it is reasonable that it is not sensitive to local distortions. 





Our observation data comes from SDSS DR12 <cit.>. We take the North Galactic Cap (NGC) as our galaxy zoo region because of the continuous, large coverage of the NGC. Thus we have the spectrum of galaxies covering huge volume in the red-shift space in NGC. We obtain the data of galaxies in the range of 130<RA<230, 10<Dec<50, and 0.07<z<0.14, including their coordinates, spectral redshifts and Petrosian magnitudes in the r band. The volume of this region is quite similar with the volume of the Bolshoi simulation. Then we use the same cosmological parameters as the Bolshoi simulation to calculate the luminosity distances of the galaxies using their spectral redshifts. We also calculate the Petrosian absolute magnitudes of the galaxies. But we didn't correct their absolute magnitudes. So the magnitudes of the galaxies are in fact in the r^0.1 band as described in <cit.>, since the Luminosity Function we use to generate the HAM mock galaxy catalogue is also in the r^0.1 band. The K-correction for this sample is also not important because we have chosen the red shift range to be 0.07<z<0.14. According to the discussion in <cit.>, the K-correction for r^0.1 is no more than 0.1. So we just ignore K-correction for our observation sample. Now we have generated an observation galaxy catalogue with 3-D space positions and reasonable absolute magnitudes. The volume covered by the sample is also similar to the volume covered by the HAM mock sample.

Before SDSS DR12, there was no such large continuous angular coverage as well as deep redshift coverage. The public access to the SDSS DR12 allows us to play with a unprecedented, continuous, large scale, 3-D galaxy zoo. It is a good time to search for unusual statistical measurement of the large scale structures, not just the percolation analysis we proposed here but also many other ways of viewing the structures as a whole web.

However, we have to point out that the observational sample is in redshift space while the HAM mock sample is in real space. So the direct comparison of these two samples is not fair. We must take redshift distortion into consideration. Fortunately, our HAM mock sample also contains information of peculiar velocity, which means we can easily transfer the HAM mock sample into the redshift space.



We also need to consider the incompleteness of the observational sample. As shown in the DR12 paper <cit.>, that is mainly because of fiber collision. We can also simulate the incompleteness by excluding those haloes in the simulation that looks closer than 62" projected on the 2D plane due to fiber collision, which will exclude ∼ 5% halos mainly in high density region. From Fig. <ref> we can see that, this incompleteness has negligible effect for the S-bb relation.
Therefore, the S-bb relation is very robust against the observational selection bias introduced by fiber collision. We also find that, changing from the adopted cosmological parameters in both Bolshoi and MDR1 simulations to Planck cosmology, the S-bb relation will not be changed significantly. However, due to the limited boxsize of Bolshoi simulation, the difference between the S-bb relation of MDR1 and Bolshoi simulations is observable. That will cause the misidentification of the transition threshold by 0.04, but for bb>0.8, the S-bb relation difference is negligible. Therefore, using the HAM mock sample built from the Bolshoi simulation to compare with observation is fair, but we need to be aware of a small misidentification of the transition threshold, because of the limited size of the simulation.




 §.§ Comparison with SDSS DR12


Here we are comparing the observational galaxy sample and the HAM mock galaxy sample in two luminosity ranges -22<M<-21 and -21<M<-20, as shown in Fig. <ref>. We choose these two ranges because, firstly, we need generally enough number of galaxies to reduce the Poisson error and secondly, we require the galaxies to be detectable even at red shift z=0.14 to get rid of Malmquist bias.  Interestingly, although the projected 2-point correlation functions of the HAM mock sample and observed sample are quite similar, their S-bb relations are different by 2-3 σ, particularly the transition point in the S_2-bb relation. According to our PCCET, this is most probably because the HAM mock sample doesn't reproduce the high order correlation functions well enough. This is also supported by a recent work <cit.> comparing Illustris simulation and the data sample generated with the same two-point correlation function. The measured S-bb relations (called 'giant component-linking length' relation), for the two samples are significantly different. The transition threshold is also smaller for the real simulation sample. 



§ CONCLUSION AND DISCUSSION

We have proposed the S-bb relation as a branch of percolation analysis of the cosmic web. We have shown that the S-bb relation and its transition threshold reveal physical information of the cosmic web. It also has some advantages compared to the standard 2-point (n-point) correlation function:


  * The S-bb relation is based on the Friend-of-Friend (FoF) algorithm, which is easy to implement and costs very little computational resources.

  * Our theoretical model, PCCET, shows that the S-bb relation contains information of all orders of correlation functions.

  * The S-bb relation is robust against redshift distortion. 

  * The S-bb relation is robust against incompleteness, which is very useful in observation. We have shown that 90% completeness is quite acceptable and this has been reached in SDSS.

  * The S-bb relation gives information on the bias of the correlation function.  

However, it also has some limitations:


  * We need large sample of galaxies to look at the S-bb relation, with at least 100h^-1Mpc box side length.

  * We need large continuous region of observation to study the S-bb relation, especially for 3-D analysis.

  * It is not easy to extract every order of the correlation functions from the S-bb relation.

Cellular or clumpy? That's not the question any more. It is well known that the distribution of galaxies is cellular. How cellular is the cosmic web? That's the question. We suggest the S-bb relation as a new measurement, which is applicable to a wide range of data types from simulation particles to galaxy catalogue, fast to compute (as fast as 2-point correlation function measurement), robust (against redshift distortion and incompleteness), and it contains information of all orders of correlation functions. In the framework of PCCET, we have good understanding of the S-bb relation. We have compared the S-bb relation of the HAM mock galaxy catalogue with SDSS DR12 galaxy catalogue. A significant difference was found, though these two samples had similar projected 2-point correlation functions. The percolation analysis with discrete points can help us understand the cosmic web. At the minimal, it provides a quick check of whether our mock sample really represents the real distribution or not. 



§ APPENDIX



 §.§ Proof for PCCET

The series in Eq. <ref> is known q-Pochhammer function or Euler function in mathematics.
The q-Pochhammer symbol (a;q)_n is defined as 

    (a;q)_n     = ∏_k = 0^n-1(1-aq^k)   
        = (1-a)(1-aq)(1-aq^2)⋯(1-aq^n-1).

In particular, the Euler function ϕ(q) is defined as the infinite product

    ϕ(q) = (q;q)_∞,


    q = 1 - p.
                             
The S_n(p) can be expressed in terms of the q-Pochhammer symbols 

    S_n(p) = 1/n∑_k=0^n-1(q;q)_k,

where

    0 < q < 1.

Since 0<q<1, it is not difficult to realise the fact that

    1 > (q;q)_1 > (q;q)_2 > ⋯ (q;q)_n > 0.
                        
Therefore

    S_n(q) > 1/nn(q;q)_n-1 = (q;q)_n-1.

On the other hand,

    S_n(q) = 1/n∑_k = 0^n - 1(q;q)_k = 1/n∑_k = 0^m - 1(q;q)_k + 1/n∑_k = m^n - 1(q;q)_k.

Here we introduce an integer m such that

    m = INT(√(n)),

where INT(x) is defined as

    INT(x) = max{n ∈ Z | n ≤ x}.

Since as n increases, (q;q)_n decreases, and we have 

    S_n(q)     = 1/n∑_k = 0^m - 1(q;q)_k + 1/n∑_k = m^n - 1(q;q)_k  
        < m/n(q;q)_0 + n - m/n(q;q)_m 
        = m/n + n - m/n(q;q)_m.

In summary,

    (q;q)_n-1 < S(n) < m/n + n - m/n(q;q)_m.

m as a function of n has the following asymptotic properties

    lim_n →∞ m = ∞,   
       lim_n →∞m/n = 0.

For simplicity, we denote

    L_n(q) = (q;q)_n-1,       U_n(q) = m/n + n - m/n(q;q)_m.

Then 

    lim_n →∞ L_n(q) = (q;q)_∞ = ϕ(q), 
    lim_n →∞ U_n(q) = (q;q)_∞ = ϕ(q).

Since we always have

    L_n(q) < S_n(q) < U_n(q),

according to the Squeeze Theorem, we concluded that

    lim_n →∞ S_n(q) = ϕ(q).

In the form of a step function, we find the limit of S_n(p) to be:

    lim_n →∞ S_n(p) = 
              
                    0    if  p = 0, 
    ϕ(1-p)    if  0 < p < 1, 
    
                    1    if  p = 1.



§ ACKNOWLEDGEMENTS


J. Zhang thanks Prof. Y. Jing for giving quite useful advices, L. Yang and Z. Li for helping with measuring 2-point correlation function, Dylan for useful discussion, S. Liao for useful discussion and advice throughout this project, J. Chen for his great help in studying the PCCET mathematically in depth. 

“The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP).
The MultiDark database was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064.”
    
    The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de).

    The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center.



