Constraining $f(R)$ Gravity with a $k$-cut Cosmic Shear Analysis of the Hyper Suprime-Cam First-Year Data

Using Subaru Hyper Suprime-Cam (HSC) year 1 data, we perform the first $k$-cut cosmic shear analysis constraining both $\Lambda$CDM and $f(R)$ Hu-Sawicki modified gravity. To generate the $f(R)$ cosmic shear theory vector, we use the matter power spectrum emulator trained on COLA (COmoving Lagrangian Acceleration) simulations. The $k$-cut method is used to significantly down-weight sensitivity to small scale ($k>1 \ h {\rm Mpc }^{-1}$) modes in the matter power spectrum where the emulator is less accurate, while simultaneously ensuring our results are robust to baryonic feedback model uncertainty. We have also developed a test to ensure that the effects of poorly modeled small scales are nulled as intended. For $\Lambda$CDM we find $S_8 = \sigma_8 (\Omega_m / 0.3) ^ {0.5} = 0.789 ^{+0.039}_{-0.022}$, while the constraints on the $f(R)$ modified gravity parameters are prior dominated. In the future, the $k$-cut method could be used to constrain a large number of theories of gravity where computational limitations make it infeasible to model the matter power spectrum down to extremely small scales.


I. INTRODUCTION
Weak gravitational lensing offers a unique test of gravity on cosmological scales. The precision of these cosmological tests of gravity will increase in the coming decade thanks to large amounts of precise data coming from Stage IV experiments including Euclid 1 [2], the Nancy Grace Roman Space Telescope 2 [3] and the Rubin Observatory 3 .
With the improved statistical precision of these surveys, extreme care must be taken to account for and remove systematic errors which could introduce biases in the parameter inference. Modeling biases are of particular concern as weak lensing is sensitive to scales down to k ∼ 10 hMpc −1 in the matter power spectrum, deep into the nonlinear regime [4]. While it is not the primary focus of this work, the impact of baryonic feedback on these scales is also highly uncertain [5].
To perform parameter inference without loss of information, one must be able to model small scales over a large volume in cosmological parameter space. It is not possible to do this analytically [6], but one can obtain the nonlinear power spectrum from an emulator [7][8][9] (or halo model [10]) trained (calibrated) on a suite of O(100) [8] high resolution N-body simulations, run over a large volume of cosmological parameter space.
Simulations must be run with more than a trillion particles to meet the percent-level matter power spectrum accuracy requirements of upcoming surveys [11]. This makes a search for deviations from general relativity (GR) challenging as it is likely infeasible to run a suite of N-body simulations at this precision for every modified theory of gravity which we would like to test, particularly since we must also account for baryonic feedback. Thus, we must cut scales from the cosmic shear data vector.
Care must be taken to retain useful information while removing poorly modeled small scales to avoid bias. The k-cut [12] cosmic shear estimator and its generalization to the combination of cosmic shear and galaxy clustering [13] (k-cut 3 × 2 point statistics [14]) were developed to meet these criteria. 4 The k-cut cosmic shear method works by applying the Bernardeau-Nishimichi-Taruya (BNT) transformation [16] to the tomographic weak lensing power spectra C ij ( ). This reorganizes the information so that it is possible to relate angular scales, , to physical scales, k, in the matter power spectrum. Then angular scale cuts correspond to physical scale cuts, making it easier to remove small-scale data.
In this paper we use this method to constrain Hu-Sawicki f (R) gravity [17] using the Hyper Suprime-Cam Year 1 cosmic shear data [16]. To generate the f (R) matter power spectrum, P (k, z), we use the emulator developed in [1,18]. As this emulator was trained on fast COLA (COmoving Lagrangian Acceleration) simulations [19,20], it becomes inaccurate above k = 1 hMpc −1 . Therefore we take a k-cut at this scale removing the smaller-scale data while preserving useful information at larger scales.
The primary aim of this work is to provide a roadmap for future tests of modified gravity with Stage IV photometric surveys. We envisage a three step process: • Step 1: For each theory of gravity, develop a model for the nonlinear power spectrum using an emulator trained on simulations as in this work, or using a halo model approach [21][22][23][24].
• Step 2: Determine the k-mode, k max , where the emulator fails i.e. the scale at which the accuracy of the power spectrum model fails to meet the requirements of the survey. • Step 3: Use k-cut cosmic shear to perform the inference removing sensitivity to scales with k > k max .
In this paper, we first review k-cut cosmic shear in Sec. II. We discuss the Subaru Hyper Suprime-Cam Year 1 Survey (hereafter HSCY1) data and covariance matrix in Sec. III. Finally, in Sec. IV, we present the k-cut cosmic shear parameter constraints for both f (R) gravity and ΛCDM before concluding in Sec. V. Unless explicitly stated otherwise, we take the Planck 2018 cosmology [25] as our baseline cosmology.

A. Lensing power spectrum
We use the lensing power spectrum, C( ), to extract the two-point information of the shear field. Under the Limber [26,27], spatially-flat universe [28], flat-sky [27], reduced shear [29,30] and Zel'dovich approximations [27] it is given by Here χ is the radial comoving distance, P is the matter power spectrum, χ H is the distance to the horizon and q(χ) is the lensing efficiency kernel, defined as: In this expression, H 0 is the Hubble constant, Ω m is the fractional matter density parameter today, c is the speed of light, a is the scale factor and n i (χ ) is the probability distribution of the effective number density of galaxies as a function of comoving distance. We also account for intrinsic alignments due to galaxy alignments with the local tidal field. The total angular power spectrum is given by where C ij IG ( ) accounts for correlations between tidally aligned foreground galaxies with gravitationally sheared background galaxies and C ij II ( ) accounts for the autocorrelation between local tidally aligned galaxies. These are given by and We use the extended nonlinear alignment model (eNLA) [31] used in the HSCY1 analysis [32] and Dark Energy Survey (DES) Year 1 analysis [33]. This is similar to the standard nonlinear alignment model (NLA) [34], but the amplitude of the intrinsic alignments, A IA (z) is allowed to vary with redshift following where A IA and α are free nuisance parameters and we fix z 0 = 0.62 [35].
In all that follows we use CosmoSIS [36] to compute the lensing spectrum using CAMB [37] to generate the linear power spectrum and Halofit [38] to generate the ΛCDM nonlinear power spectrum. We describe how the nonlinear f (R) power spectrum is generated in Sec. II E.
B. The BNT transform and k-cut cosmic shear Weakly lensed galaxies probe structure along the entire line-of-sight. Lensing structures of the same physical scale at different redshifts subtends different angles on the sky, breaking the correspondence between angular and physical scales. Hence, taking an angular scale cut does not correspond to taking a physical scale cut.
In order to find a more precise correspondence between physical and angular scales, in the k-cut method [12,14] we make use of the Bernardeau-Nishimichi-Taruya (BNT) transform [15]. Specifically we take a linear combination of lensing kernels, q i , so thatq  where we sum over repeated indices, M ij is the BNT matrix 5 andq i (z) are the BNT lensing kernels. 5 The BNT matrix is implicitly a function of cosmology, but we compute it once at the fiducial cosmology and leave it fixed for The BNT matrix is found by first setting M ii = 1 for all i and M ij = 0 for i < j, the remaining BNT matrix elements are found by solving the system and z max is the maximum redshift of the survey.
We use the publicly available code at: https://github.com/pltaylor16/x-cut to compute the BNT matrix.
The top row of Fig. 1 shows the four photometric redshift bins used in the HSCY1 cosmic shear analysis [35]. The corresponding lensing kernels are plotted in the second row. The kernels are broad in z, which as we have argued makes it difficult to relate physical and angular scales. Meanwhile the BNT kernels are shown in the bottom row of Fig. 1. These are much narrower in z, which will allow us to relate physical and angular scales at the two-point level.
Since the lensing power spectrum depends on two sets of kernels labeled by tomographic bin numbers i and j, the correct transform at the two-point level is where we sum over repeated indices. We refer to C ij ( ) as the k-cut cosmic shear power spectrum. The power spectra (top right) and BNT-transformed power spectra (bottom right) are shown in Fig. 2. The linear and nonlinear contributions are shown in blue and orange respectively and the combination of the two is shown in black. We expect a bigger change between the BNT-transformed power spectrum and the standard power spectrum for higher redshift correlation, which we indeed find. Crucially the linear and non-linear power spectra contributions are comparable at higher-in the BNT-transformed case. This implies that we can take scale cuts at higher -values in the BNT-transformed case, while still removing sensitivity to poorly modeled nonlinear scales.
Formally, since each kernel is narrow in z, we can define a "typical distance," χ i , which is given by comoving distance at the redshift where kernel q i (z) obtains its maximum value 6 , assuming a fiducial cosmology. Then the remainder of this work. It was found in [13] that the choice of fiducial cosmology has a negligible impact. 6 One could alternatively take χ i as the co-moving distance at the average redshift over the kernel. Top left: tomographic theoretical angular power spectra, C( ) for the best-fit cosmology of HSCY1. The tomographic bin combination is indicated in the top left of each subplot. Black lines are used for the total angular power spectra, while blue and orange lines are used to indicate the linear and nonlinear contributions from the nonlinear matter power spectrum respectively. Bottom right: BNT-transformed angular power spectra, C( ). Due to the sorting of physical scales, the nonlinear contributions become important at higher-modes than for standard tomographic angular power spectra. This implies we can take a cut at a higher-mode and extract more information while remaining robust to modeling uncertainties at small physical scales.
if we want to remove sensitivity to scales smaller than some k-mode in the matter power spectrum, we can use the Limber relation and cut all > k χ ij , where χ ij = min( χ i , χ j ). This procedure is referred to as kcut cosmic shear.
The effectiveness of the method at removing sensitivity to small scales comes with several caveats. The BNT-transformed kernels, q j (z), shown in Fig. 1 still have some width so the -k correspondence is not exact. This will become less problematic in future surveys as the number of tomographic bins is increased and photometric redshift errors improve. Since the BNT transform depends on n j (z) the effectiveness of the transform also relies on the accuracy of the photometric redshift estimates. 7 In light of these issues, we perform a test in Sec. IV A to ensure that k-cut method cuts sensitivity to small scales as intended.
C. k-cut cosmic shear likelihood and covariance We develop our own module to evaluate the k-cut cosmic shear likelihood in CosmoSIS [39], a modular opensource package for likelihood-based cosmological parameter inference.
We assume that the k-cut cosmic shear likelihood is Gaussian. This will need to be verified in a follow-up study following the method in [40] or [41]. The data and theory vectors are found by BNT transforming the theory and data vectors of the standard angular power spectra, C ij ( ).
The covariance matrix of the BNT transformed data vector, C ab,cd BNT ( , ), can be computed from an estimate of the covariance of the tomographic angular power spectra. The elements are given by where C ef,gh ( , ) denotes the covariance between C ef ( ) and C gh ( ). We list analogous expressions for the covariance elements in the k-cut 3 × 2 point formalism [14] in the Appendix A. One can also use the likelihood sampling method to compute the covariance of the BNT-transformed data vector (see [42] for more details).
In this paper we used the covariance matrix, C ef,gh ( , ), computed in [35] at the best-fit cosmology therein. This covariance includes all Gaussian, non-Gaussian and super-sample covariance terms. We refer the reader to the Appendix of [32] for more details.
The HSCY1 and k-cut cosmic shear correlation matrices 8 are shown in Fig. 3. The matrices are ordered into block matrices with tomographic bins ordering {i, j} = 11, 12, 13, 14, 21, 22, 23, 24, 33, 34, 44 and increasing inside each block. The k-cut covariance is much sparser. This is expected as by construction the BNT transform reduces the redshift overlap between the lensing efficiency kernels (see Fig. 1) weakening the covariance between different tomographic bins.

D. f (R) gravity
In f (R) gravity, the Einstein-Hilbert action S is modified by including an extra function of the Ricci curvature, f (R) [43] so that, where L M gives the Lagrangian of matter in the Universe, G the gravitational constant and the speed of light, c, is set to unity. In principle the extra degrees of freedom in this model can drive cosmic acceleration. 8 For a covariance matrix, C, the correlation matrix is Top: the angular power spectra correlation matrix from [35], computed at the best fit cosmology from the HSCY1 analysis. The matrices are ordered into block matrices with tomographic bin ordering {i, j} = 11, 12, 13, 14, 21, 22, 23, 24, 33, 34, 44 and increasing inside each block. Bottom: the k-cut cosmic shear correlation matrix. As expected the k-cut covariance is much sparser as the BNT transform reduces the redshift overlap between the lensing efficiency kernels (see Fig. 1) weakening the covariance between tomographic bins.
In this paper we consider the Hu-Sawicki parametrization [17], which to date has passed all observational tests in the GR limit [44] and is well-studied in the literature, making it an ideal choice for the first application of the k-cut framework to a test of modified gravity. In this parametrization where m = H 0 √ Ω m and n, c 1 and c 2 are free parameters. Matching the background expansion to that of a ΛCDM universe imposes as additional constraint on the derivative of f (R) with respect to background Ricci scalar, ∂f (R) ∂R , evaluated today .
The Hu-Sawicki model is typically parameterized in terms of two free parameters: |f R0 | and n, 9 and the theory reduces to GR in the limit f R0 → 0.

E. f (R) power spectrum emulator
We use the f (R) Hu-Sawicki modified gravity Gaussian process emulator developed in [1] and available at https://github.com/LSSTDESC/mgemu to compute the nonlinear f (R) power spectrum. This emulator was trained inside a five-dimensional parameter space using a suite of f (R) comoving Lagrangian acceleration (COLA) simulations [20]. The parameters and their prior ranges are listed in Table I. The Hubble parameter is fixed as h 0 = 0.67 as assumed in [1]. Note in particular that we restrict log(f R0 ) to a flat prior with f R0 ∈ [10 −4 , 10 −8 ] as this is the range in which the f (R) emulator is accurate [1].
COLA uses a combination of second order Lagrangian perturbation theory (2LPT) and an N-body component to retain most of the accuracy of full N-body simulations at a fraction of the computational cost. The emulator is accurate at the 5% level up to k = 1 hMpc −1 , achieving 1% accuracy for f R0 < 10 −5 . We refer the reader to [1] for more details.
The emulator computes the ratio of the f (R) and ΛCDM nonlinear power spectra i.e. P f (R) (k, z)/P ΛCDM (k, z).
We use Halofit [38] to compute P ΛCDM (k, z) and we have developed our own CosmoSIS module to call the emulator and compute P f (R) (k, z).
Although baryonic feedback represents a subdominant effect on scales larger than 1 hMpc −1 [5], it is noted in Appendix A of [35] that inaccuracies of the Halofit dark matter power spectrum model can result in 1σ biases on σ 8 for a cut at = 2000. We expect this to serve as an upper bound on the biases in this k-cut study since [10] found that the Halofit model is accurate to within ∼ 4% below 1 hMpc −1 and typically accurate to ∼ 2% below 1 hMpc −1 over most of this k-range. As we will show in Sec. IV B this is smaller than our measurement error on S 8 . Nevertheless future studies can sidestep this issue by using a more accurate emulator such as [45] or [46].

A. Hyper Suprime-Cam year 1 data
We use data from the first year of the Hyper Suprime-Cam (HSC) survey [16]. The data covers 136.9 deg 2 over 6 disjoint fields with an effective source number density of 17 galaxies per arcmin 2 [35]. We refer the reader to [16] for a detailed description of the shear estimation.
Galaxies were binned into 4 redshift bins inside the range 0.3 < z < 1.5 with a magnitude cut of i < 24.5 [35] and redshifts were estimated from the COS-MOS 30-band photo-z catalog (see [32,[47][48][49][50] for more details). The resulting redshift probability density functions, n j (z), for each bin are plotted in Fig. 1. Power spectra were estimated in 15 logarithmically spaced bandpowers in the range 30 < < 6500 using the pseudo-C pipeline described in Sec. II B . We cut all < 300 as in [35] and take additional -cuts after applying the BNT transform following the k-cut cosmic shear procedure. In this paper we remove sensitivity to all scales corresponding to k larger than k = 1 hMpc −1 which corresponds to cutting bandpowers with effective > 1205, 1814, 2329 and 2781 for the four tomographic bins respectively. In contrast the HSCY1 analysis took a global multipole cut at = 2000 [35].

B. Residual systematics
As in [32], we consider two sources of residual bias: multiplicative shear bias and linear photometric redshift error.
The theoretical angular power spectra are rescaled by a multiplicative factor where Here m i sel and m i R are fixed parameters which account for the selection bias and a residual bias in the shear response in tomographic bin i (see [32] Sec. 5.7 for more details), while ∆m i is a residual multiplicative left as a free parameter in this analysis. We use the same values and prior ranges as in [32] 10 , but they are outlined in Table II for convenience. 10 Unlike in [32] each tomographic bin is assigned its own residual multiplicative bias as there is no a priori reason to assume that they should be the same for each bin We also allow for a linear bias in the redshift probability distribution function n i (z) → n i (z + ∆z i ). (17) We use the same prior ranges for ∆z i as in the fiducial HSCY1 [32] analysis. They are displayed in Table II. We do not marginalize over any free parameters to account for point spread function residuals, as this was found to have negligible impact in [35].

A. Verification of k-cut method
In this section we perform a test to ensure the kcut method removes sensitivity to the desired scales. We fix the cosmological parameters choosing a Planck 2018 cosmology and constrain the amplitude, A P , of the power spectrum above the target k cut = 1 hMpc −1 using: 1) a k-cut analysis and 2) a standard C( ) analysis. In both cases the largest bandpower has an effective = 3080. The amplitude, A P , is taken relative to the Planck 2018 prediction so that A P = 1 corresponds to the matter power spectrum for a Planck cosmology. Since the k-cut data vector should be insensitive to modes k cut > 1 hMpc −1 by construction, if the errors on A P are significantly larger in the k-cut case, this implies that small scales are nulled as intended.
Formally we rescale the power spectrum where Θ(k cut ) = 0, if k < k cut 1, if k ≥ k cut (19) and use the Emcee sampler [51] in CosmoSIS to constrain A P . The results of this analysis are shown in Fig. 4. The error on the amplitude, σ(A P ) increases from 0.06 to 0.31 or a factor of 5.5 validating that sensitivity to small scales is significantly down-weighted. We recommend this validation step is performed in all future k-cut analyses.
A similar rescaling of the power spectrum was used to search for deviations from ΛCDM in [52] using The Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS) data. We find no evidence that our constraints on A p are inconsistent with ΛCDM. Although we find the majority of the posterior lies below A p < 1 in both the C( ) and k-cut analysis, this is the expected behavior for two reasons. First, we did not include baryonic physics which is known to suppress the power spectrum above k = 1 hMpc. Second, we use the Planck 2018 cosmology. Weak lensing constraints are known to favor lower values for σ 8 , i.e. the amplitude of the power spectrum.

B. ΛCDM comparison with HSCY1 analysis
Using our CosmoSIS pipeline and the Multinest sampler [53], we perform the first k-cut likelihood analysis using real data with a target k cut = 1 hMpc −1 . This scale cut was chosen as it is the point where the f (R) power spectrum becomes inaccurate [1] to maintain consistency with our f (R) results. Conveniently this choice also removes sensitivity to scales where baryonic corrections to the nonlinear matter power spectrum [5,13] become important. Thus we use the dark matter-only power spectrum without marginalizing over any additional nuisance parameters to account for baryonic feedback.
The resulting HSCY1 k-cut parameter constraints in the Ω m − σ 8 and Ω m − S 8 planes are shown Fig. 5 and the full constraints are shown in Fig. 6. We find S 8 = 0.789 +0.039 −0.022 in the k-cut analysis. Despite taking a very conservative choice of scale cut, the symmetric error on S 8 increases by just 3% compared to the HSCY1 constraints.

C. f (R) results
We perform a k-cut analysis to constrain Hu-Sawicki f (R) gravity using the nonlinear emulator developed in [1] to generate the nonlinear power spectrum. As in the ΛCDM case we take a target k cut = 1 hMpc −1 corresponding to the scale where the emulator's accuracy exceeds the 5% level [1]. We take our priors to cover the region of parameter space where the emulator was trained (see Table I), while for the nuisance and intrinsic alignment parameters we take the same priors as in the ΛCDM case. We show the results of the 2 Hu-Sawicki parameters in Fig. 7. Our constraints, log(f R0 ) = −6.38 +0.94 −1.41 and n = 1.8 +1.1 −1.5 , are almost completely prior dominated. We also repeat the likelihood analysis fixing n = 1 as in the analysis of [55]. In this case we find log(f R0 ) = −6.67 +1. 37 −0.88 . These results are consistent with those in [55] which found that log(f R0 ) could not be constrained within a uniform prior −8 < log(f R0 ) < −2 using Kilo-Degree Survey (KiDS) data. External con- Gauss(0, 0.0135) ∆z 3 Gauss(0, 0.0383) ∆z 4 Gauss(0, 0.0376) TABLE II. Prior range used in the fiducial k-cut ΛCDM analysis. The parameters are arranged into cosmological parameters (top), intrinsic alignment parameters (middle) and nuisance parameters (bottom). The prior ranges are the same as in the fidicial HSCY1 analysis [35], except we assign a multiplicative bias to each tomographic bin as there is no a priori reason to expect the multiplicative bias in each bin to be the same.
Although the region of parameter space containing ΛCDM is not contained in the prior, i.e., (f R0 = 0), there is a very weak preference for small values of log(f R0 ) consistent with ΛCDM.

V. CONCLUSION
In this work, we have performed the first ever kcut cosmic shear analysis to constrain ΛCDM and Hu-Sawicki f (R) gravity. The k-cut method is ideally suited to test theories of modified gravity for which we do not have accurate models for the nonlinear power spectrum as deep into the nonlinear regime as for the ΛCDM case. We choose a target k-cut of k cut = 1 hMpc −1 for both analyses as this is the scale where the nonlinear f (R) emulator becomes inaccurate.
By taking a simple linear transformation of the original data vector, the BNT transform tightens the correspondence between angular and physical scales. An-  6. The 68% and 95% confidence regions for the k-cut analysis with a target k-cut of 1 hMpc −1 . It is worth comparing with the HSCY1 fiducial analysis shown in Fig. 21 of [35]. Despite the extremely conservative choice of scale cuts, the majority of the information is preserved (see also Fig. 4). All contour plots appearing in this work were generated using Chainconsumer [54] and the width of the kernel density estimator (KDE) bandpass is set to 1.5 unless explicitly stated otherwise.
gular scale cuts can then be used to cut sensitivity to poorly modeled small scales. The method comes with virtually no additional computational overhead and the strength of the − k relationship will only tighten as the number of tomographic bins is increased and photometric redshift error is reduced with future datasets.
The k-cut cosmic shear method comes with several caveats. The BNT-transformed kernels although nar-row in z-space still have some width, so the relationship between and k is not exact. Furthermore the BNT transformation itself depends on the tomographic selection function n j (z), so the efficacy of the transformation at separating physical scales could in principle be compromised by photometric redshift errors.
To ensure that the method removes sensitivity to small scales as intended, we fix the cosmological param- eters and constrain a single free parameter, A p , which controls the amplitude of the matter power spectrum above the target scale cut in k-space. We find the error on this amplitude increase by a factor of 5.5 in the k-cut case demonstrating that the sensitivity to these scales is significantly reduced as intended. We recommend this test is performed in all future k-cut cosmic shear and k-cut 3 × 2 point [14] studies.
In the ΛCDM case, we find S 8 = 0.789 +0.039 −0.022 . Our choice of scale cut ensures the constraints are robust to baryonic physics model uncertainties. Despite a conservative choice of scale cut, the symmetric error on S 8 increases by just 3% relative to the HSCY1 angular power spectrum analysis [35].
To constrain f (R) gravity, we used the nonlinear power spectrum emulator developed in [1]. After marginalizing over all other parameters we find log(f R0 ) = −6.38 +0.94 −1.41 and n = 1.8 +1.1 −1.5 . Although we find weak lensing constraints on f (R) gravity are prior dominated, this will soon change as data from Euclid, Roman and the Rubin Observatory arrives over the next decade. The k-cut method, implemented for the first time in this paper, is a promising approach to test modified gravity into the nonlinear regime while avoiding model bias. We would also like to Alex Hall for pointing out a more streamlined approach to compute the BNT covariance matrix and Vincenzo Cardone for useful discussions.