title: Breaking the Dark Degeneracy with the Drifting Coefficient of the Field Cluster Mass Function
authors: Suho Ryu Jounghun Lee Marco Baldi
text: 
  =0.0in
  =cmbx10 scaled 1000 =cmr10 scaled 1000
  =cmr8 scaled 1000 =cmr7 scaled 1000
  et al.
  i.e.
  e.g.
  etc.
  JCAP
  
      ∑ m_ννΛ CDM  eV h^-1 M_⊙ h^-1 Mpc1Astronomy Program, Department of Physics and Astronomy, FPRD, 
          Seoul National University, Seoul 08826, Korea shryu@astro.snu.ac.kr, jounghun@astro.snu.ac.kr2Dipartimento di Fisica e Astronomia, Alma Mater Studiorum Università di Bologna, viale Berti Pichat, 
          6/2, I-40127 Bologna, Italy3INAF - Osservatorio Astronomico di Bologna, via Ranzani 1, I-40127 Bologna, Italy4INFN - Sezione di Bologna, viale Berti Pichat 6/2, I-40127 Bologna, Italy
          
          We present a numerical analysis supporting the evidence that the redshift evolution of the drifting coefficient of the field cluster mass function 
          is capable of breaking several cosmic degeneracies. 
          This evidence is based on the data from the  CoDECS and  DUSTGRAIN- pathfinder simulations performed separately for various 
          non-standard cosmologies including coupled dark energy, f(R) gravity and combinations of f(R) gravity with massive neutrinos 
          as well as for the standard ΛCDM cosmology. We first numerically determine the field cluster mass functions at various redshifts 
          in the range of 0≤ z≤ 1 for each cosmology. Then, we compare the analytic formula developed in previous works with the numerically 
          obtained field cluster mass functions by adjusting its drifting coefficient, β, at each redshift. It is found that the analytic formula 
          with the best-fit coefficient provides a good match to the numerical results at all redshifts for all of the cosmologies. The empirically 
          determined redshift evolution of the drifting coefficient, β(z), turns out to significantly differ among different cosmologies. 
          It is also shown that even without using any prior information on the background cosmology the drifting coefficient, β(z), can 
          discriminate with high statistical significance the degenerate non-standard cosmologies not only from the ΛCDM but 
          also from one another. It is concluded that the evolution of the departure from the Einstein-de Sitter state and spherically symmetric collapse 
          processes quantified by β(z) is a powerful probe of gravity and dark sector physics. 
          
          
          
          § INTRODUCTION
          
          
          A cosmic degeneracy refers to the circumstance that a standard diagnostic fails to distinguish between different cosmologies with high statistical significance. 
          For example, the cluster mass function, which is regarded as one of the most powerful probes of cosmology based on the large 
          scale structure, is unable to discriminate the effect of a low amplitude of the linear density power spectrum from that of massive neutrinos (ν) 
          (dubbed the σ_8-∑ m_ν degeneracy) in a νΛCDM (massive neutrinos ν + cosmological constant Λ + Cold Dark Matter) 
          cosmology. 
          Another example that the cluster mass function fails to discriminate from the ΛCDM cosmology is a coupled dark energy (cDE) model in which a scalar field 
          DE coupled to DM particles follows a supergravity potential <cit.>. 
          Since the cosmic degeneracy is caused by the limited sensitivity of a given standard diagnostic on which the degenerate models have almost the same effects, 
          what is required to break it is to overcome the limitation by utilizing prior information from other independent diagnostics. In the aforementioned example, 
          the σ_8-∑ m_ν degeneracy can be broken by prior information on the large-scale amplitude of the linear density power spectrum from the CMB 
          observations.
          
          There are, however, a few cosmic degeneracies which have been found more difficult to break even by combining the priors from several independent diagnostics. 
          A notorious example is the cosmic degeneracy between the ΛCDM + GR and the νCDM + MG cosmologies, where GR and MG stand for the 
          general relativity and modified gravity, respectively <cit.>. All different versions of the MG theory adopt a common tenet that the apparent acceleration of the present Universe 
          is caused not by the dominance of the anti-gravitational Λ at the present epoch but by the deviation of the gravitational law from the prediction of GR on cosmological scales. The consequence of this tenet is the existence of a long-range fifth force, which in turn has an effect of enhancing the density power spectrum 
          on the scales comparable to those affected by the suppression due to free streaming massive neutrinos <cit.>. 
          
          In the theory of f(R) gravity, the gravitational dynamics is defined by a modified Einstein-Hilbert action functional to which 
          an arbitrary function of the Ricci scalar, f(R), is introduced as a substitution for the Ricci scalar R itself of the original action in GR <cit.>. 
          Choosing as a viable MG the f(R) gravity in which the Ricci scalar term, R, in the Einstein-Hilbert action functional is replaced by an arbitrary function, 
          f(R), <cit.> numerically investigated a possible cosmic degeneracy between the ΛCDM + GR and the νCDM + f(R), 
          and demonstrated that the two cosmologies cannot be discriminated from each other by several standard diagnostics such as the nonlinear density power spectra, 
          halo bias and cluster mass functions <cit.>. The nonlinear growth rate functions, cluster velocity dispersions, and tomographic higher-order weak lensing statistics were proposed 
          in subsequent works as candidate diagnostics that could be capable of breaking this cosmic degeneracy <cit.>, also employing Machine Learning techniques <cit.>.
          
          Very recently, <cit.> have developed a new independent diagnostics based on the evolution of the drifting coefficient of the field cluster abundance and shown that 
          this new diagnostics is capable of distinguishing between those dynamical DE cosmologies which are degenerate with the ΛCDM case in their linear density 
          power spectra and cluster mass functions. In the follow-up work of <cit.>, it was also found that the aforementioned σ_8-∑ m_ν degeneracy can 
          be in principle broken by this new diagnostics. Our goal here is to explore whether or not this new diagnostics can break other cosmic degeneracies including that 
          between the ΛCDM + GR and the νCDM + MG models. 
          
          The key contents of the upcoming Sections are as follows. Section <ref> will present a succinct review of the analytic model for the field cluster abundance on which 
          the new diagnostics is based. Section <ref> will present a numerical evidence that the new diagnostics can distinguish between the standard 
          ΛCDM and a cDE cosmology that produce very similar cluster mass functions. Section <ref> will present a proof for the validity of the analytic formulae for 
          the field mass function and drifting coefficient in a MG gravity cosmology. Section <ref> will present a numerical evidence that the new diagnostics can break 
          the degeneracy between the standard ΛCDM+GR and the νCDM+MG cosmologies. Section <ref> will be devoted to summarizing the final results and 
          discussing the follow-up works.
          
          
          
          § A SUCCINCT REVIEW OF THE ANALYTIC MODEL
          
          
          According to the generalized excursion set theory that incorporates the non-spherical collapse conditions and non-Markovian random-walk 
          process <cit.>, the differential mass function of the cluster halos can be written as:
          
          
              d N(M, z)/d  ln M = ρ̅_m/M|d  ln σ^-1/d  ln M|f[σ(M, z)] ,
          
          where ρ̅_m is the mean mass density of the universe, σ (M) is the standard deviation of the linear density contrast field 
          smoothed on the mass scale M, and f(σ) is the multiplicity function characterized by two coefficients, β and D_B: 
          
          
              f(σ ; D_B,β)    ≈    f^(0)(σ ; D_B, β) + 
              f^(1)_β=0(σ ; D_B) + 
              f^(1)_β(σ ; D_B, β) + 
              f^(1)_β^2(σ ; D_B, β)  , 
              
              f^(0)(σ ; D_B, β)    =   δ_sc/σ√(1+D_B)√(2/π) 
              e^-(δ_sc+βσ^2)^2/2σ^2(1+D_B)  ,
              
              f^(1)_β=0(σ ; D_B)    =   -κ̃δ_sc/σ√(2a/π)[e^-aδ_sc^2/2σ^2
              -1/2Γ(0,a δ_sc^2/2 σ^2)]  , 
              
              f^(1)_β(σ ; D_B, β)    =   
              -β a δ_sc[f^(1)_β=0(σ; D_B)+κ̃ erfc(δ_sc/σ√(a/2))]  , 
              
              f^(1)_β^2(σ ; D_B, β)    =   β^2a^2δ^2_scκ̃{erfc(δ_sc/σ√(a/2))+
                    σ/aδ_sc√(a/2π)[e^-aδ_sc^2/2σ^2(1/2-a δ_sc^2/σ^2)+3/4aδ_sc^2/σ^2Γ(0,a δ_sc^2/2 σ^2)]}  ,
          
          with a≡ 1/(1+D_B), κ̃ = κ a, κ = 0.475, upper incomplete gamma function Γ(0, x), complementary 
          error function erfc(x), and critical density contrast for the spherical collapse δ_sc. 
          The diffusion coefficient, D_B, is a measure of the stochasticity of δ_c caused by the influence of the perturbing neighbors as well as by the 
          uncertainty in the identification of a cluster halo, while the drifting coefficient, β, quantifies the deviation of the critical density contrast for the realistic non-spherical 
          collapse, δ_c, from δ_sc at a given σ (M).
          
          Suggesting that for the isolated field clusters, the degree of the stochasticity of δ_c should be negligible, (i.e., D_B=0), <cit.> modified the 
          Corasaniti-Achitouv formalism to construct a single parameter model for the field cluster mass function as 
           
          
              dN_ I(M, z)/dln M = ρ̅/M|d  ln σ^-1/d  ln M|
              f[σ(M,z) ; D_B=0,β]  .
          
          In spite of having only single parameter, this analytic model has been shown to be very successful in describing the field cluster 
          mass functions in a wide redshift range not only for the standard ΛCDM cosmology but also for several non-standard 
          cosmologies including dynamical dark energy or massive neutrinos <cit.>.
          
          Although the exact value of δ_sc has been known to weakly depend on the background cosmology as well as on the redshift 
          <cit.>, <cit.> regarded δ_sc as a constant, setting it at the Einstein-de Sitter value of 1.686<cit.>, 
          as done in the original formulation of the generalized excursion set mass function theory <cit.>. 
          In reality, the gravitational collapse proceeds in a non-spherical way, for which the actual critical density contrast, δ_c, 
          departs from the idealistic spherical threshold, δ_sc. The cosmology dependence of δ_c is expected to overwhelm 
          that of δ_sc, given that the degree of the non-sphericity of the collapse process is closely linked with the anisotropy of the cosmic web, which in turn 
          possesses strong dependence on the background cosmology <cit.>. 
          Unlike δ_sc, however, the value of δ_c and its link to the initial conditions cannot be analytically derived from first principles 
          due to the complexity associated with the non-spherical collapse process <cit.>. 
          
          <cit.> imparted any redshift and cosmology dependences of the collapse threshold at a given mass scale to the empirical parameter, β, 
          and showed that the redshift evolution of the empirically determined β(z) has a universal form of an inverse sine hyperbolic function of z, 
          regardless of the background cosmology: 
          
          
              β(z) = β_A sinh^-1[1/q_z(z-z_c)]  ,
          
          where β_A, q_z and z_c, represent three adjustable parameters, whose best-fit values depend on the background cosmology. 
          Especially, the critical redshift, z_c, defined as β(z_c)=0 (i.e., δ_c(z_c)=δ_sc=1.686), has been found to depend most 
          sensitively on the background cosmology <cit.>, as it reflects not only how severely the real gravitational collapse deviates from the 
          spherical symmetry but also how rapidly the universe evolves away from the Einstein-de Sitter state.
          
          
          
          § EVOLUTION OF THE DRIFTING COEFFICIENT IN NON-STANDARD COSMOLOGIES
          
          
           §.§ Effect of Coupled Dark Energy on β(z)
          
          
          A cDE cosmology describes an alternative universe where the role of DE is played by a dynamical scalar field, ϕ, coupled to DM particles through energy-momentum exchange. The DE-DM coupling that causes the time-variation of DM particle mass <cit.>
          generates a long-range fifth force via which the growth of structures can be enhanced <cit.>.
          Categorized by the shape of DE self-interaction potential, V(ϕ), as well as by the strength of the DE-DM coupling, s(ϕ)≡ -dln m_ DM/dϕ, 
          a cDE cosmology has recently attained delving attentions since it has been found to provide a possible solution to the Hubble tension <cit.>.
          
          To investigate the effect of cDE on the redshift evolution of β(z), we utilize the data from the Large Coupled Dark Energy Cosmological 
          Simulations ( L-CoDECS) run by <cit.> with a modified version of the  GADGET3 code, a non-public developers version of the widely-used public code 
          GADGET-2 <cit.>. The  L-CoDECS is a series of N-body cosmological runs that simulate a standard ΛCDM and five different cDE cosmologies 
          on a periodic box of linear size 1  h^-1Gpc containing 1024^3 collisionless DM particles of 
          individual mass m_ DM=5.84× 10^10 as well as an equal number of collisionless baryon particles of m_ baryon=1.17× 10^10. 
          The initial conditions of the standard ΛCDM cosmology were chosen to meet the constraints from the Seven-Year Wilkinson Microwave Anisotropy 
          observations <cit.>.
          The five different cDE cosmologies are divided into three categories: the constant DM-DE coupling and exponential potentials (EXP001, EXP002, EXP003), 
          the exponential DM-DE coupling and exponential potential (EXP008e3) and the constant coupling and supergravity potential (SUGRA). 
          All cDE cosmologies simulated by the L-CoDECS were ensured to have a flat geometry, sharing the same values of the five key cosmological parameters, 
          h = 0.703, Ω_ CDM = 0.226, Ω_ DE = 0.729, Ω_b = 0.0451, A_s = 2.42 × 10^-9 and n_s = 0.966. 
          They differ from one another only in the potential shape and DM-DE coupling as well as in the linear density power spectrum amplitude, information on which 
          are provided in the first four columns of Table <ref>. 
          For more detailed description of the cDE cosmologies and the L-CoDECS [All data are available at the CoDECS website, http://www.marcobaldi.it/web/CoDECS.html], 
          we refer the readers to <cit.>.
          
          The  L-CoDECS simulations have been released with catalogs of gravitationally bound halos identified for each cosmology through a two-step process starting with a 
          Friends-of-Friends (FoF) algorithm with linking length parameter of l_c=0.2 followed by a gravitational unbinding procedure of each individual FoF halo using the 
           SUBFIND algorithm <cit.> that allows to associate spherical overdensity masses and radii to each gravitationally bound main substructure. Selecting the 
          cluster halos with virial masses M≥ 3× 10^13 from the halo catalog and applying again the FoF algorithm with linking length parameter of l_sc=2× l_c 
          to such halo sample, we identify the clusters of cluster halos as marginally bound superclusters. Here, the choice of l_sc=2× l_c for the FoF finding of the superclusters 
          is made to guarantee D_B=0 for the field cluster halos, as explained in <cit.>.
          
          Sorting out the field cluster halos as the superclusters which have only one member cluster halo, we create 
          a sample of the field cluster halos at each redshift in the range of 0≤ z≤ 1. The differential mass function of the field clusters halos is determined by 
          computing the number density of the field cluster halos, dN_ I/dln M, whose masses fall in the differential bin of the logarithmic mass, [ln M, ln M + dln M], 
          at each redshift. As done in <cit.>, we employ the Jackknife method to compute the errors in the determination of dN_ I/dln M. Splitting the sample of the field clusters 
          into 8 Jackknife subsamples of equal size, we obtain dN_ I/dln M separately from each subsample and then calculate the one standard deviation scatter of dN_ I/dln M 
          among the 8 subsamples as errors at each logarithmic mass bin. 
          
          Using the linear density power spectrum of each cDE cosmology provided within the  CoDECS public data release, we evaluate the linear density rms fluctuation, 
          σ(M), and the analytic mass function of the field cluster, Equation (<ref>), as well. 
          The best-fit value of the drifting coefficient, β in Equation (<ref>) is determined at each redshift by minimizing the following χ^2 under the assumption 
          that there is no correlation in dN_ I/dln M at different mass bins:
          
          
              χ^2(β)=∑_i=1^N_p[n^ n(ln M_i)-n^ a(ln M_i;β)]^2/σ^2_n_i  ,
          
          where n^ n(ln M_i) and n^ a(ln M_i;β) are the numerically obtained and analytically evaluated values of 
          dN_ I/dln M at the ith logarithmic mass bin, respectively,  N_ p is the number of the logarithmic mass bins and σ_n_i is the Jackknife error 
          on n^ n(ln M_i). The error on β is then determined as the square root of the inverse of the second derivative of χ^2 with respect to β as 
          σ_β=(d^2χ^2/dβ^2)^-1/2<cit.>.
          
          As done in <cit.>, once the values of β(z) are determined at various redshifts, we fit them to Equation (<ref>) via a nonlinear least square 
          regression procedure with the SciPy python code <cit.> to find the best-fit values of the three parameters, β_A, q_z and z_c and their associated 
          errors σ_β_A, σ_q_z and σ_z_c, respectively (see Table <ref>).  
          Then, we calculate the statistical significance of the differences in the three parameters among the cosmologies as 
          Δβ_A/σ_Δβ_A, Δ q_z/σ_Δ q_z and Δ z_c/σ_Δ z_c where 
          Δβ_A, Δ q_z and Δ z_c are the differences in the three parameters between two cosmologies, while 
          σ_Δβ_A, σ_Δ q_z and σ_Δ z_c correspond to the propagated errors in the determination of the differences. 
          
          Figure <ref> (Figure <ref>) plots the numerically determined mass functions of the field cluster halos (filled black circles) as well as the analytic model 
          (red solid line), Equation (<ref>), with the best-fit value of β for the six cosmologies at z=0 (z=1), respectively. In each panel, the analytic model for the 
          ΛCDM case (black dashed line) is also plotted to clearly show the differences. Although the analytic model, Equation (<ref>), succeeds in matching the numerical results at both of the redshifts for all of the cDE cosmologies, the field cluster mass functions are found to be incapable of telling apart with high statistical significance 
          the three cosmologies ΛCDM, EXP001 and SUGRA at both of the redshifts, z=0 and 1.
          
          Figure <ref> plots the redshift evolution of the empirically determined drifting coefficient, β(z) (filled black circles) as well as the fitting formula 
          (red solid lines) for the six cosmologies. In each panel, the fitting formula for the ΛCDM case (black dashed line) are also plotted to show the differences. As can be seen, 
          the fitting formula expressed in terms of the inverse sine hyperbolic function, Equation (<ref>), with the best-fit values of q_z, β_A and z_c indeed describes 
          quite well the behaviors β(z) for all of the six cosmologies. Note that the SUGRA can be distinguished by β(z) from the ΛCDM with high statistical 
          significance despite that the two cosmologies are mutually degenerate in the cluster mass functions. The statistical significance of the difference in the critical redshift parameter, 
          z_c, between the ΛCDM and the SUGRA cosmologies is found to be as high as 7.48 σ.
          Although β(z) distinguishes with high statistical significance the other cDE cosmologies except for the EXP001 from the ΛCDM, 
          it fails to break the degeneracy between the ΛCDM and the EXP001 cases, due to the extremely weak DM-DE coupling of the latter cosmology. 
          
          We have so far used prior information on the background cosmology for the determination of dN_ I/dln M and β(z). In other words, to examine if β(z) can break a 
          cosmic degeneracy between two different cosmologies, we assume that information on the shape of the linear density power spectrum are available. 
          In practice, however, this prior information is not available for the determination of β(z). Especially, if a background cosmology is indistinguishable from the ΛCDM 
          case by the standard diagnostics, then it may not be justified to make such a preemptive assumption about the shape of the linear density power spectrum. 
          The EXP001 corresponds to this case where no prior information on the background cosmology should be assumed to be available in practice, since the standard diagnostics 
          including the linear density power spectrum, mass function, etc., are unable to distinguish it from the ΛCDM case. 
          
          To deal with this degeneracy, we use the linear density power spectrum of the ΛCDM case, P(k;Λ CDM), for the computation of σ(M) in 
          dN_ I/dln M for the EXP001 case and compare the reevaluated analytic model with the numerical results to find the best-fit β(z). That is, 
          we redetermine β(z) for the EXP001 case without using prior information on P(k; EXP001). 
          Figure <ref> plots the analytical mass function of the field clusters (red solid lines) obtained by using P(k;Λ CDM) and compares it with the numerical results (black filled circles) for the EXP001 case at z=0 (top panel) and z=1 (bottom panel). As can be seen, in spite of no prior information on the background cosmology, the analytical mass function of the field clusters still describes quite well the 
          numerical results at both of the redshifts for the EXP001 case. 
          
          Figure <ref> plots the same as the top-right panel of Figure <ref> but without using prior information on P(k; EXP001). 
          As can be seen, the EXP001 turns out to yield larger differences in β(z) from the ΛCDM. The best-fit values of β_A, q_z and z_c for the EXP001 case listed in Table <ref> correspond to the ones obtained without using prior information 
          on P(k; EXP001). The statistical significance of the difference in z_c between the ΛCDM and the EXP001 is found to be as high as 2.53. 
          Figure <ref> summarizes the statistical significance of the difference in z_c among the three cosmologies, ΛCDM, EXP001, and SUGRA, 
          which are mutually degenerate in the field cluster mass functions. Although the degeneracy between the ΛCDM and the EXP001 can be broken by 
          β(z) only with 2.53 significance, we speculate that a larger data set would improve the significance. 
          
          
          
           §.§ Effect of f(R) Gravity on β(z)
          
          
          In the theory of f(R) gravity, the strength of a long range fifth force is quantified by the absolute value of the derivative of f(R) with respect to the 
          Ricci scalar R at the present epoch, | f_R0|≡| df/dR|_0. A larger value of | f_R0| corresponds to a stronger fifth force, 
          which would more severely enhance the small-scale density power <cit.>. 
          If neutrinos have a non-zero mass in a f(R) gravity cosmology, however, the suppressing effect of the free streaming neutrinos on the small-scale density power spectrum 
          could compensate the enhancing effect of the fifth force, resulting in a suppression of the deviations from the standard ΛCDM + GR cosmology that each of these two scenarios would individually imprint on structure formation. In other words, the linear density power spectra may not be capable of distinguishing a certain combination of 
          f_R0 with ∑ m_ν from the standard ΛCDM + GR cosmology, since they could have zero net effect on the amplitude of small-scale density perturbations <cit.>. 
          
          To investigate if β(z) can also break the cosmic degeneracy between ΛCDM+GR and νCDM+f(R), we use a subset of the data from 
          the  DUSTGRAIN- pathfinderN-body simulation suite that were conducted by <cit.> on a box of volume 750^3 h^-3Mpc^3 
          for various νCDM+f(R) cosmologies as well as the ΛCDM+GR cosmology. The  DUSTGRAIN- pathfinder simulations were performed with the  MG-GADGET code <cit.>– another modified version of  GADGET-3 implementing an adaptive mesh solver for the spatial fluctuations of the f_R scalar degree of freedom – to trace the evolution of 
          768^3 DM particles of mass 8.1× 10^10 h^-1M_⊙. To simulate the νCDM+f(R) cosmologies, the 
           DUSTGRAIN- pathfinder adopted the widely-used realisation of f(R) proposed by <cit.> and a particle-based implementation of massive neutrinos developed by <cit.>. 
          Collapsed structures were identified through a FoF finder with a linking length parameter of l_c=0.16 followed by the unbinding procedure implemented in the  SUBFIND code to identify the halo center and its spherical overdensity mass and radius for all gravitationally bound objects in each cosmology, similarly to what described above for the  CoDECS simulations. For a detailed description of the technical details of the  DUSTGRAIN- pathfinder simulations, 
          see <cit.>. 
          
          Among the various cosmologies simulated by the  DUSTGRAIN- pathfinder, we consider three different CDM+f(R) (namely, fR4, fR5 and fR6 corresponding to 
          | f_R0|=10^-4, 10^-5 and 10^-6, respectively) and three different νCDM+f(R) (namely fR4+0.3, fR5+0.15 and fR6+0.06 
          corresponding to =0.3, 0.15 and 0.06, respectively) as well as the standard ΛCDM + GR (from here on, GR) with initial conditions 
          set at the Planck values <cit.>. These 7 different cosmologies were ensured to be flat, described by the common key cosmological parameter values of 
          h= 0.67, Ω_m = 0.31, Ω_ DE = 0.67, Ω_b = 0.0481, A_s = 2.2 × 10^-9 and n_s = 0.97. 
          The first four columns of Table <ref> list the values of | f_R0|, ∑ m_ν, σ_8 for 
          each of the seven cosmologies considered in the present work.
          
          We first examine whether or not the analytic model for the field cluster mass function, Equation (<ref>), is valid for the three CDM+f(R) cosmologies. 
          Analyzing the FoF halo catalogs in the redshift range 0≤ z≤ 1 and following the same procedure described in Section <ref>, we numerically determine 
          dN_ I/dln M for the GR, fR4, fR5 and fR6 cases. 
          To evaluate the analytic model, Equation (<ref>), and compare it with the numerically determined dN_ I/dln M to derive β(z) for 
          each of the three f(R) gravity cosmologies, we use the  MGCAMB code <cit.>.
           
          Figure <ref> (Figure <ref>) depicts the same as Figure <ref> (Figure <ref>) but for the fR4, fR5 and fR6 cases, 
          revealing that the analytic model matches quite well the numerical results even for the f(R) gravity models. As expected, the fR4 (fR6) yields the most (least) abundant field 
          clusters in the entire mass range. No statistically significant difference is found in dN_ I/dln M between the GR and the fR6 cases, indicating their mutual degeneracy 
          in the field cluster mass functions. Figure <ref> plots the same as Figure <ref> but for the fR4, fR5 and fR6 cases. As can be seen, despite that the 
          field cluster mass functions fail to distinguish between the GR and the fR6 cases, the field cluster drifting coefficient, β(z), can break the degeneracy, showing a substantial 
          difference between the two cosmologies. 
          
          It is worth noting the distinct redshift dependence of the difference in β(z)  between the GR and each f(R) gravity cosmology.  
          The fR4 case yields an almost redshift-independent shift of  β(z) from that of the GR case, while the other two cases show different redshift-dependent shifts 
          between each other. That is, for the fR5 (fR6) case, the largest deviation of β(z) from that of the GR case occurs at the low (high) redshift ends. 
          A qualitative explanation for this trend is provided in the following.
          For the fR4 case, the fifth-force is basically always unscreened at low redshifts, which implies that the haloes of all masses are equally affected by the 
          fifth force, and that there is no sharp transition between screened and unscreened regions in the cosmic web. 
          
          Whereas, for the fR5 and fR6 cases, as the screening properties imply that the massive halos are screened, while less massive halos are not. 
          This introduces a mass dependence in the deviation of the halo mass function from that of the GR case. In particular, for the fR6 case, the high-mass tail of the halo mass function 
          should be almost unaffected and thus there should be an enhancement in the number of smaller-mass halos. This may have a different impact on the bias of halos at different 
          masses, and consequently an impact on the definition of the field clusters (i.e. isolated massive objects) which could induce a different evolution of β(z).
          A more quantitative analysis is required to fully understand the distinct differences in β(z) between the GR and each f(R) cosmology, 
          which is, however, beyond the scope of this paper.
          
          
          
           §.§ Combined Effect of f(R)+ν on β(z)
          
          
          Now that the validity of the analytic model of the field cluster mass function for the f(R) gravity cosmology is confirmed, we repeat the whole 
          process but for the fR4+0.3, fR5+0.15 and fR6+0.06 cosmologies, which were shown to be degenerate with the GR in the standard statistics 
          including the nonlinear density power spectrum, cluster mass functions and halo bias <cit.>. 
          Figure <ref> (Figure <ref>) depicts the same as Figure <ref> (Figure <ref>) but for the 
          fR4+0.3, fR5+0.15 and fR6+0.06 cosmologies. As can be seen, the analytic model is still quite valid in matching the numerically obtained field cluster 
          mass functions even for the f(R)+ν cosmologies. At z=0, the three f(R)+ν cosmologies show no difference from the GR case in the field cluster mass 
          functions. At z=1, the differences in dN_ I/dln M between the f(R)+ν and the GR cases are slightly larger but still not statistically significant. 
          Figure <ref> plots the same as Figure <ref> but for the fR4+0.3, fR5+0.15 and fR6+0.06 cosmologies.
          As can be seen, the fR4+0.3 cosmology yields a substantial difference in β(z) from the GR case, in spite of their mutual degeneracy 
          in the standard statistics. Yet, both of the fR5+0.15 and the fR6+0.06 cosmologies show almost no difference in β(z) from the GR case. 
          
          As done in Section <ref>, we redetermine dN_ I/dln M for both of the fR5+0.15 and fR6+0.06 cases without using prior information 
          on the shapes of their power spectra, which are plotted in Figure <ref>. As can be seen, the analytic model, Equation (<ref>), still agrees quite 
          well with the numerically obtained field cluster mass functions for both of the cases at both of the redshifts, despite that P(k;GR) is substituted for P(k; fR5+0.15) 
          and P(k; fR6+0.06). The drifting coefficient, β(z), redetermined without using prior information is plotted in Figure <ref>, which reveals 
          that the three cosmologies yield much larger differences in β(z). 
          
          For each cosmology, we determine the best-fit values of β_A, q_z and z_c by fitting Equation (<ref>) to β(z) obtained without priors. Then, we calculate 
          the statistical significance of the differences in the three parameters among the three cosmologies, which are shown in Figure <ref>. As can be seen, without using prior information 
          on the linear density power spectra of the f(R)+ν cosmologies, the statistical significance of the differences in z_c between the GR and the fR5+0.15 and between the 
          fR6+0.06 and the fR5+0.15 are as high as 3.48 and 3.22, respectively.
          
          Meanwhile, for the fR6+0.06 case, it turns out to be not z_c but β_A that is able to distinguish it from the GR case with Δβ_A/σ_Δβ_A=2.03. The lower statistical significance of the differences in β(z) between the GR and the fR6+0.06 is likely to be at least partially due to the large errors caused by the relatively small 
          box size of the DUSTGRAIN-pathfinder simulations. Given the distinct behaviors of β(z) between the the GR and the fR6+0.06 shown in Figure <ref>, 
          we suspect that if a halo sample from a larger simulations were used, the statistical significance would increase.
          The best-fit values of β_A, q_z and z_c for each of the seven cosmologies simulated by the DUSTGRAIN-pathfinder are shown in Table <ref>. 
          For the fR6+0.06 and fR5+0.15 cosmologies that are degenerate with the GR case in the standard statistics, what is listed in Table <ref> is the best-fit values 
          obtained without using priors in the shapes of the linear density power spectra.
          
          
          
          § SUMMARY AND DISCUSSIONS
          
          
          The new diagnostic developed by <cit.> traces the redshift evolution of the drifting coefficient, β(z), which is a single parameter of the analytic model 
          for the field cluster mass function derived by <cit.> in the frame of the generalized excursion set theory <cit.>.
          Motivated by the recent finding that this diagnostic can in principle break the σ_8-∑ m_ν degeneracy <cit.>, 
          we have studied whether or not the same diagnostic can break the degeneracy between the non-standard and the standard ΛCDM+GR cosmologies by 
          utilizing the data from the CoDECS and DUSTGRAIN-pathfinder simulations. 
          For this study, we have considered eleven different non-standard cosmologies which include 5 different cDE 
          (EXP001, EXP002, EXP003, EXP008e3 and SUGRA), 3 different f(R) gravity (fR4, fR5, fR6), and 3 different f(R) gravity+ν cosmologies 
          (fR4+0.3eV, fR5+0.15eV, fR6+0.06eV). 
          
          Among the cDE and f(R) gravity cosmologies, the EXP001 and fR6 have been known to be very similar to the ΛCDM+GR in the linear density 
          power spectra and cluster mass functions at z=0, due to their extremely weak DM-DE coupling and fifth force, respectively. 
          The three f(R) gravity+ν cosmologies have been known to be degenerate not only with the ΛCDM+GR but also among one another in the 
          standard diagnostics that include the cluster mass functions, halo bias, and nonlinear density power spectrum <cit.>.
          Analyzing the catalogs of the FoF bound objects for each cosmology, we have identified the field cluster halos and determined their mass functions at each redshift in the range 
          of 0≤ z≤ 1. The best-fit value of the drifting coefficient β has been found by adjusting the analytic model of <cit.> to the numerically determined field cluster 
          mass functions. The fitting formula for β(z) proposed by <cit.> have been used to assess the statistical significance of the 
          differences in β(z) among the cosmologies. This analysis has lead us to find the following.
          
            *  
          The analytic model of <cit.> for the field cluster mass functions with the best-fit values of β agrees excellently well with the numerical results at all redshifts for 
          all of the non-standard cosmologies. 
          
            * 
          The empirical formula of <cit.> for β(z) works fairly well for all of the non-standard cosmologies. 
          
            * 
          Despite that they produce very similar (field) cluster mass functions, the ΛCDM and the SUGRA cosmologies substantially differ in β(z) from each other. 
          
            * 
          The degeneracy between the ΛCDM and the EXP001 in the (field) cluster mass functions can be broken by β(z) with 2.53σ significance 
          without using any prior information on the linear density power spectrum, P(k; EXP001). 
          
            * 
          The degeneracy between the ΛCDM+GR and the fR4+0.3eV in the linear density power spectra and (field) cluster mass functions can be broken by 
          β(z) with high statistical significance. 
          
            * 
          The degeneracy among the ΛCDM+GR, fR5+0.15eV and fR6+0.05eV cosmologies in the standard diagnostics can be broken by β(z) with minimum 
          2.01σ significance, without using any prior information on the linear density power spectra. 
          
          
          To understand the advantage of using β(z) as a cosmology discriminator, it may be worth comparing β(z) with the standard diagnostics such as the linear density 
          power spectrum, nonlinear density bi spectrum and cluster mass function.
          As for the linear density power spectrum, it deals with isotropically averaged densities and thus fail to capture independent information contained in the anisotropic 
          nonlinear cosmic web about the background cosmology <cit.>. 
          As for the nonlinear density bi spectrum that treats the nonlinear anisotropic density field, it is not readily observable, suffering from 
          highly nonlinear halo bias and redshift space distortion effects. 
          Regarding the cluster mass function, although it is free from the halo bias and redshift space distortion effect, it varies most sensitively with the value of σ_8.
          If two different cosmologies share an identical value of σ_8 (e.g., ΛCDM and SUGRA), the cluster mass function is apt to fail in telling them apart. 
          
          Meanwhile, the field cluster drifting coefficient, β(z), deals with the non-spherical collapse occurring in the anisotropic cosmic web that contains additional information on the 
          initial conditions. It is free from the halo bias and redshift space distortion effect, directly quantifying how the background cosmology deviates from the Einstein-de Sitter state which 
          sensitively depends on the evolution of the energy contents of the universe.
          Notwithstanding, we have yet to find a direct link of β(z) to the initial conditions, which weakens its power as a probe of gravity and dark sector physics.
          The very fact that the inverse sine hyperbolic function provides a fairly good approximation to the empirically determined β(z) for all of the cosmologies hints that 
          it should be beyond a mere fitting formula. Our future work is in the direction of theoretically deriving β(z) from first principles, providing a physical explanation for 
          why β(z) behaves as an inverse sine hyperbolic function of z and establishing its direct link to the initial conditions. 
          
          Another advantage is high observational applicability of β(z). As shown in Figure <ref>,  the fR5+0.15eV and fR6+0.06eV cases substantially differ 
          from the GR case in the low-z values of β(z) (z<0.5). In other words, it does not require a large sample of the high-z clusters with z>0.5 to distinguish between the 
          f(R)+ν and the GR cases with β(z) in practice.  
          Yet, to distinguish between the fR5+0.15eV and the fR6+0.06eV cases as well as between the f(R) gravity and the cDE cosmologies with β(z), it indeed requires  
          a large cluster sample from a wide range of redshifts.  The upcoming large-scale deep surveys such as the LSST (Large Synoptic Survey Telescope) or EUCLID 
          that is expected to cover the redshift range up to z∼ 2<cit.> will improve prospects for β(z) as a cosmological discriminator. 
          
          
          
          We thank an anonymous referee for useful comments and suggestions. 
          JL acknowledges the support by Basic Science Research Program through the National Research Foundation (NRF) of Korea 
          funded by the Ministry of Education (No.2019R1A2C1083855). 
          MB acknowledges support by the project "Combining Cosmic Microwave Background and Large Scale Structure data: an Integrated Approach for Addressing Fundamental Questions in Cosmology", funded by the MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2017 - grant 2017YJYZAH. The N-body simulations described in this work have been performed on the Hydra supercomputer at RZG and on the Marconi supercomputer at Cineca thanks to the PRACE allocation 2016153604 (P.I. M. Baldi).
          
          
          000[Amendola(2000)]ame00 
          Amendola, L. 2000, , 62, 043511 
          [Amendola(2004)]ame04 
          Amendola, L. 2004, , 69, 103524 
          [Amendola et al.(2013)]euclid 
          Amendola, L., Appleby, S., Bacon, D., et al. 2013, Living Reviews in Relativity, 16, 6
          [Baldi et al.(2010)]bal-etal10 
          Baldi, M., Pettorino, V., Robbers, G., et al. 2010, , 403, 1684 
          [Baldi(2012a)]codecs1 
          Baldi, M. 2012a, , 420, 430
          [Baldi(2012b)]codecs2 
          Baldi, M. 2012b, , 422, 1028
          [Baldi et al.(2014)]bal-etal14 
          Baldi, M., Villaescusa-Navarro, F., Viel, M., et al. 2014, , 440, 75
          [Bond & Myers(1996)]BM96 
          Bond, J. R., & Myers, S. T. 1996, , 103, 1
          [Buchdahl(1970)]buc70
          Buchdahl, H. A. 1970, , 150, 1
          [Clifton et al.(2012)]mgreview 
          Clifton, T., Ferreira, P. G., Padilla, A., et al. 2012, , 513, 1
          [Corasaniti & Achitouv(2011a)]CA11a 
          Corasaniti, P. S. & Achitouv, I. 2011a, , 106, 241302
          [Corasaniti & Achitouv(2011b)]CA11b 
          Corasaniti, P. S. & Achitouv, I. 2011b, , 84, 023009
          [Di Valentino et al.(2020)]cde_hubble20 
          Di Valentino, E., Melchiorri, A., Mena, O., et al. 2020, Physics of the Dark Universe, 30, 100666
          [Eke et al.(1996)]eke-etal96 
          Eke, V. R., Cole, S., & Frenk, C. S. 1996, , 282, 263
          [Giocoli et al.(2019)]gio-etal19 
          Giocoli, C., M., Baldi, M., & Moscardini, L. 2018, , 481, 2813
          [Garcia-Farieta et al.(2019)]gf-etal19
          García-Farieta, J. E., Marulli, F., Veropalumbo, A., et al. 2019, , 488, 1987
          [Gunn, & Gott(1972)]GG72 
          Gunn, J. E., & Gott, J. R. 1972, , 176, 1 
          [Hagstotz et al.(2019)]hag-etal19 
          Hagstotz, S., Costanzi, M., Baldi, M., et al. 2019, , 486, 3927
          [Hojjati et al.(2011)]hoj-etal11 
          Hojjati, A., Pogosian, L. & Zhao, G. 2011, , 1108, 005
          [Hu & Sawicki(2007)]HS07 
          Hu, W., & Sawicki, I. 2007, , 76, 064004
          [Komatsu et al.(2011)]wmap7 
          Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, , 192, 18
          [Lee(2012)]lee12 
          Lee, J. 2012, , 752, 40
          [Lewis et al.(2000)]camb 
          Lewis, A., Challinor, A., & Lasenby, A. 2000, , 538, 473
          [Li & Barrow(2007)]LB07 
          Li, B., & Barrow, J. D. 2007, , 75, 084010 
          [Mainini & Bonometto(2006)]MB06 
          Mainini, R., & Bonometto, S. 2006, , 74, 043504 
          [Macciò et al.(2004)]mac-etal04 
          Macciò, A. V., Quercellini, C., Mainini, R., et al. 2004, , 69, 123516 
          [Maggiore & Riotto(2010a)]MR10a
          Maggiore, M., & Riotto, A. 2010a, , 711, 907
          [Maggiore & Riotto(2010b)]MR10b
          Maggiore, M., & Riotto, A. 2010b, , 717, 515
          [Mangano et al.(2003)]man-etal03 
          Mangano, G., Miele, G., & Pettorino, V. 2003, Modern Physics Letters A, 18, 831 
          [Merten et al.(2019)]mer-etal19
          Merten, J., Giocoli, C., Baldi, M., et al. 2019, , 487, 104
          [Naidoo et al.(2020)]nai-etal20 
          Naidoo, K., Whiteway, L., Massara, E., et al. 2020, , 491, 1709
          [Pace et al.(2010)]pac-etal10 
          Pace, F., Waizmann, J. C., & Bartelmann, M. 2010, , 1865, 1874
          [Peel et al.(2018)]pel-etal18
          Peel, A., Pettorino, V., Giocoli, C., et al. 2018, A&A, 619, A38
          [Peel et al.(2019)]pel-etal19
          Peel, A., Lalande, F., Stark, J.-L., et al. 2018, , 100, 023508
          [Pettorino & Baccigalupi(2008)]PB08 
          Pettorino, V., & Baccigalupi, C. 2008, , 77, 103003 
          [Planck Collaboration et al.(2016)]planck16 
          Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A24
          [Puchwein et al.(2013)]mggadget 
          Puchwein, E., Baldi, M., & Springel, V. 2013, , 436, 348
          [Ryu & Lee(2020a)]RL20a 
          Ryu, S., & Lee, J. 2020a, , 889, 62
          [Ryu & Lee(2020b)]RL20b 
          Ryu, S., & Lee, J. 2020b, , 894, 65
          [Shim & Lee(2013)]SL13 
          Shim, J. & Lee, J. 2013, , 777, 74
          [Springel(2005)]Gadget2
          Springel, V. 2005, , 364, 1105
          [Springel et al.(2001)]spr-etal01
          Springel, V., White, S. D. M., Tormen, G., et al. 2001, , 328, 726
          [Starobinsky(1980)]sta80
          Starobinsky, A. A. 1980, , 91, 99
          [Tyson(2002)]lsst 
          Tyson, J. A. 2002, , 4836, 10
          [Viel et al.(2010)]vie-etal10 
          Viel, M., Haehnelt, M. G., & Springel, V. 2010, , 2010, 015
          [Virtanent et al.(2020)]vir-etal20 
          Virtanen, P. et al. 2020, Nature Methods, 17, 261
          [Wetterich(1995)]wet95 
          Wetterich, C. 1995, , 301, 321 
          [Wintergerst & Pettorino(2010)]WP10 
          Wintergerst, N., & Pettorino, V. 2010, , 82, 103516
          [Wright et al.(2019)]wri-etal19
          Wright, B. S., Koyama, K., Winther, H. A., et al. 2019, , 06, 040
          [Zhao et al.(2009)]zha-etal09 
          Zhao, G., Pogosian, L., Silvestri, A., et al. 2009, , 79, 083513
          [Zucca et al.(2019)]zuc-etal19 
          Zucca, A., Pogosian, L., Silvestri, A., et al. 2019, , 1905, 001
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          
          cccccccc0ptBest-fit Parameters of β(z) for the CoDECS cosmologies.Model    V(ϕ)   s   σ_8   β_A   q_z   z_cΛCDM		    -												    - 			    0.809    -0.16±0.01   0.31±0.04   1.10±0.02
          
          EXP001					   e^-0.08ϕ    0.05 	    0.825    -0.16±0.01   0.31±0.05   1.04±0.02
          
          EXP002					   e^-0.08ϕ    0.10		    0.875    -0.17±0.02   0.35±0.08   1.32±0.04
          
          EXP003					   e^-0.08ϕ    0.15 	    0.967    -0.14±0.01   0.19±0.06   1.44±0.05
          
          EXP008e3				   e^-0.08ϕ    0.40		    0.895    -0.16±0.01   0.27±0.04   1.19±0.03
          
          SUGRA003				   ϕ^-2.15e^ϕ^2/2    -0.15 	    0.806    -0.16±0.01   0.39±0.07   1.35±0.03
          
          ccccccc0ptBest-fit Parameters of β(z) for the DUSTGRAIN-pathfinder cosmologies.Model    |f_R0|   ∑ m_ν [ eV]   σ_8   β_A   q_z   z_cΛCDM		    - 					    0.0		    0.847    -0.11±0.01   0.22±0.06   1.24±0.03
          
          fR4							   10^-4    0.0		    0.967    -0.10±0.01   0.16±0.04   1.39±0.03
          
          fR5							   10^-5    0.0 		    0.903    -0.16±0.02   0.50±0.11   1.40±0.04
          
          fR6							   10^-6    0.0 		    0.861    -0.08±0.01   0.09±0.04   1.24±0.04
          
          fR4+0.3eV				   10^-4    0.3 		    0.893    -0.09±0.01   0.29±0.09   1.52±0.06
          
          fR5+0.15eV			   10^-5    0.15		    0.864    -0.15±0.05   0.85±0.45   1.73±0.14
          
          fR6+0.06eV			   10^-6    0.06		    0.847    -0.08±0.01   0.11±0.04   1.27±0.04
  