Cosmological constraints from weak lensing scattering transform using HSC Y1 data
Sihao Cheng Gabriela A. Marques Daniela Grandón Leander Thiele Masato Shirasaki Brice Ménard Jia Liu


firstpage–lastpage
UTF8gkai 
Multivariate Fidelities
    Theshani NuradhaSchool of Electrical and Computer Engineering, Cornell University, Ithaca, NY 14850, USA
  (pt388@cornell.edu, hemant.mishra@cornell.edu, wilde@cornell.edu).
Hemant K. Mishra[1]
Felix LeditzkyDepartment of Mathematics &
Illinois Quantum Information Science and Technology (IQUIST) Center, University of Illinois Urbana-Champaign, Urbana, IL 61801, USA
  (leditzky@illinois.edu)  
Mark M. Wilde[1]

     
=================================================================================================================================================================================================================================================================================================================================================================================================================================




As weak lensing surveys go deeper, there is an increasing need for reliable characterization of non-Gaussian structures at small angular scales. Here we present the first cosmological constraints with weak lensing scattering transform, a statistical estimator that combines efficiency, robustness, and interpretability. With the Hyper Suprime-Cam survey (HSC) year 1 data, we obtain Ω_m=0.29_-0.03^+0.04, S_8≡σ_8(Ω_m/0.3)^0.5=0.83±0.02, and intrinsic alignment strength A_IA=1.0±0.4 through simulation-based forward modeling. 

Our constraints are consistent with those derived from Planck. 

The error bar of Ω_m is 2 times tighter than that obtained from the power spectrum when the same scale range is used. This constraining power is on par with that of convolutional neural networks, suggesting that further investment in spatial information extraction may not yield substantial benefits.

We also point out an internal tension of S_8 estimates linked to a redshift bin around z∼1 in the HSC data. We found that discarding that bin leads to a consistent decrease of S_8 from 0.83 to 0.79, for all statistical estimators.
We argue that 

photometric redshift estimation is now the main limitation in the estimation of S_8 using HSC. This limitation is likely to affect other ground-based weak lensing surveys reaching redshifts greater than one. Alternative redshift estimation techniques, like clustering redshifts, may help alleviate this limitation.





methods: statistical –
gravitational lensing: weak –
cosmological parameters –
large-scale structure of Universe




§ INTRODUCTION
 



Due to gravitational lensing, matter density fluctuations in our universe distort the images of background galaxies. On cosmological scales, this phenomenon is referred to as cosmic shear <cit.>. This effect can be used to probe the geometry and the matter distribution of our Universe. At relatively large angular scales, the lensing field is approximately Gaussian and can be characterized by two-point statistics, such as the power spectrum or correlation function. However, as illustrated in Figure <ref>, ongoing weak lensing surveys can probe scales small enough to reveal rich non-Gaussian structures, requiring additional statistics to extract the full cosmological information in the field.



As opposed to Gaussian fields, non-Gaussianities possess much more degrees of freedom, and different statistics may exhibit sensitivity to diverse features and systematic effects in a map.


The efficacy of non-Gaussian statistics has been explored in the weak lensing context, including on the most recent datasets

<cit.>. However, there has been no simple agreement on which summary statistic to use. Traditional moment based statistics are not robust to outliers in real data. Other simple statistics are usually too focused on particular features, and thus cannot extract the full cosmological information in a field. On the other extreme, inference with deep learning methods provide better constraining power but still lack interpretability.

Recently, a new summary statistic called the scattering transform has been proposed to extract cosmological information from weak lensing maps <cit.>.  Similar statistics using related ideas, such as the phase harmonics <cit.> and scattering spectra <cit.>, have also been proposed and studied on mock data for cosmology.  Originating from the computer vision community <cit.>, the scattering transform is shown to exhibit high information content and robustness in mock cosmological <cit.> and astrophysical <cit.> analysis. 
In this paper, we present the first application of this statistic to observational weak lensing data[Application of the scattering transform has been explored to galaxy clustering with BOSS data <cit.> but not yet on weak lensing data.]. Our analysis uses the Subaru Hyper Suprime-Cam first-year (HSC Y1) data. We employ a forward modeling approach with a large set of N-body simulations and mock catalogs customized for the data, and derive cosmological constraints.




Weak lensing data is mainly sensitive to two cosmological parameters and provides roughly uncorrelated constraints on S_8≡σ_8√(Ω_m/0.3) and Ω_m. Here, Ω_m is the total matter density today, and σ_8 represents the linear matter fluctuation within a sphere of radius 8 h^-1Mpc. The scattering transform
and in general non-Gaussian statistics are expected to significantly tighten the constraint of Ω_m by distinguishing fields with similar amplitude of lensing fluctuations but different extent into the non-linear regime. In addition, in light of the recent tension between S_8 values derived from the cosmic microwave background <cit.> and large-scale structures found in various studies <cit.>, the exploration of non-Gaussian statistics provides a consistency check and may open a new avenue for understanding the underlying cause.












§ STATISTICAL ANALYSIS OF HSC Y1 DATA
 

We perform a statistical analysis of the first year data of the HSC survey[<https://hsc.mtk.nao.ac.jp/ssp/>] <cit.>, which is a wide-area imaging survey conducted using the 8.2m Subaru telescope in Hawaii in five broad-bands, grizy. The HSC is the the deepest large-area weak lensing survey to date, with 5σ point-source depth of i∼ 26 and typically i-band seeing of 0.58". Hence, the HSC data is well suited for studying the non-Gaussian structures in the lensing field.

The first year data release of HSC, referred to as HSC S16A <cit.>, contains data taken from March 2014 to April 2016 and covers an area of 136.9 deg^2 area in six disjoint fields: , , , , , and . Due to the limited 10×10 deg^2 size of our cosmo-varied simulations, we trim off some edge regions in  and  that exceed a width of 10 degrees. 




 §.§ Weak lensing mass map



Our analysis is based on weak lensing convergence maps (also called mass maps) created using the galaxy shapes measured by <cit.> from the i-band coadded images with the re-Gaussianization PSF correction method <cit.>. Figure <ref> shows the resulting mass maps, which clearly reveals non-Gaussian structures such as peaks and ridges. Below, we describe the procedures for creating the mass maps. 

To perform a tomographic lensing analysis, we select source galaxies in 4 redshift bins, with edges [0.3, 0.6, 0.9, 1.2, 1.5]. 
Below, we will refer to those bins as bin 1 to 4, from low to high redshift.
The redshift bin assignment is determined using the “best-fit” point estimate of redshift provided by the  photo-z code <cit.>. The effective number density of galaxies in each photo-z bin is 5.14, 5.23, 3.99, 2.33 arcmin^-2, respectively. The total number of galaxies is about 8.5 million.

We begin by creating shear maps from the galaxy shape catalog. The survey area is discretized into pixels, each with a size of 0.88 arcminutes. Then, we calibrate and average 

the observed galaxy ellipticity e to obtain an unbiased estimator γ̂ of the underlying cosmic shear γ for each pixel,

    γ̂ = 1/1+ m̅_tot(e̅/2 ℛ̅ -c̅) ,

where 1+m_tot is the total multiplicative bias, c is the additive bias for each galaxy, and ℛ≡ 1-e_rms^2 denotes the responsivity. This estimator originates from taking expectation of both sides of the shear distortion equation at small shear limit: e = e^int + 2(1+m)ℛγ + c + e^mea, where e^int and e^mea are the intrinsic galaxy ellipticity and measurement error, both assumed to have zero expectation.

Following <cit.>, we include three terms in the total multiplicative bias m_tot= m + m_sel + m_ℛ, where m is that obtained from image simulations <cit.> and averaged within each of the individual fields and for each of the four redshift bins, m_sel=[0.0086, 0.0099, 0.0091, 0.0091] is the multiplicative bias caused by galaxy size selection for each redshift bin, and m_ℛ=[0.000, 0.000, 0.015, 0.030] is the correction for redshift-dependent responsivity <cit.>.

In eq. <ref>, the bar ·̅ represents averages weighted by w:

    ϕ̅(x,y) ≡∑_in pixel (x,y) w_i ϕ_i/∑_in pixel (x,y) w_i ,

where ϕ(x,y) represent the value of any quantity in {e, ℛ, m_tot, c} at the pixel (x,y), and the weights w_i are read directly from the HSC S16A catalog <cit.>, where it is defined as w=(e_rms^2 + σ_e^2)^-1, where σ_e^2 is the mean square of galaxy ellipticity in each component, estimated from all galaxies with similar color, brightness, and seeing, and σ_e^2 is the measurement error <cit.>. We further smooth our shear map with a Gaussian smoothing kernel of 2 arcmin in size to reduce the effect caused by the non-uniform weight map (galaxy counts per pixel). For Gaussian statistics, an uneven weight map can be analytically incorporated in the definition of the statistics. However, for non-Gaussian statistics except moment based ones, there is no general method to analytically incorporate the weight map, therefore a proper smoothing is usually needed. We note that compared to our two companion papers <cit.>, the smoothing algorithm used in this paper incorporates the weight map in a way closer to the optimal weighting in power spectrum analysis and thus provides lower noise level. The algorithm of smoothing and more discussions are presented in appendix <ref>. 
After smoothing, pixels with galaxy density lower than half of the mean density are masked and set to zero. For visualization purpose, we in-paint the masked regions <cit.> when generating Figure <ref>, but keep the masked regions as zero for statistical analysis. The shear maps are then converted into the mass maps using their arithmetic relation in Fourier space <cit.> and implemented assuming periodic boundary conditions for each field of view.





 §.§ Scattering transform and power spectra


We consider two types of summary statistics of the convergence maps: the power spectrum extracting Gaussian information, and scattering transform coefficients providing substantial non-Gaussian information. 

The scattering transform was recently introduced to cosmology and has desirable properties when performing statistical analyses: it compresses the information in a robust way <cit.>. It achieves the constraining power of a deep convolutional neural network without the need to train millions of neurons. The number of coefficients used is only comparable to that of the power spectrum, which is smaller than most other non-Gaussian statistics. In practice, such an internal compression of information significantly reduces the estimation noise of the covariance matrix and likelihood emulator. Compared to high-order moment based statistics, the scattering coefficients do not amplify outliers, making the estimator stable under cosmic variance and robust to outliers in real observations.

The scattering transform also serves as a bridge between traditional moment-based statistics and convolutional neural networks. In the context of weak lensing maps, it can measure the scale-dependent sparsity of structures in the field, which reflects how non-linearly the fluctuations have evolved and therefore significantly improve the constraint of Ω_m.
For a more intuitive understanding, we refer the reader to our pedagogical review <cit.>. 






From a lensing convergence field κ, the scattering transform first separates scales by wavelet convolutions, and then applies a point-wise modulus as non-linearity:

    I_1^j,θ (x,y) = |κ * ψ^j,θ| ,

where '*' represents convolution, and wavelets ψ are band-pass filters labeled by their size j and orientation θ. As shown in Figure <ref>, we use wavelets that are one-sided in Fourier space  and thus complex-valued in real space. The modulus therefore converts fluctuations into their local amplitude (i.e. envelope) so that the I_1 fields are local intensity maps of fluctuations in a frequency range. 
In particular we use directional Morlet wavelets, whose profiles in frequency space are close to off-centered Gaussians. The mathematical form is given in Appendix <ref>.



The 1st-order scattering coefficients are the average of these I_1 fields over space and orientations

    s_1(j) = ⟨ I_1⟩_x,y,θ .

The s_1 coefficients are qualitatively similar to the power spectrum at the central frequency of the wavelet ψ^j – they just differ in the sense of being L_1-norm of the I_1 fields instead of the square of L_2-norms.

To characterize more non-Gaussian features, each of the I_1 fields are further analyzed following the same procedure,

    I_2^j_1, θ_1, j_2, θ_2(x,y)    = |I_1^j_1,θ_1 * ψ^j_2, θ_2| = ||κ * ψ^j_1, θ_1| * ψ^j_2, θ_2| 
    
        s_2(j_1,j_2)    = ⟨ I_2⟩_x,y,θ_1,θ_2

For the second-order scattering coefficients, as pointed out in <cit.>, most of the cosmological information lies in the orientation-average coefficients. Therefore, we average over both orientations θ_1 and θ_2. It has also been shown that only coefficients with j_1 < j_2 are informative. Intuitively, this is because the modulus converts fluctuations into their local amplitude (envelope), which tends to be smoother than the original fluctuations. The number of scattering coefficients is therefore J+J(J-1)/2, where J is the number of wavelet scales used. 


Now we consider the power spectrum which carries only the Gaussian information. To cover exactly the same frequency range and make a fair comparison between the Gaussian and non-Gaussian information, we use the same wavelets ψ to weigh and bin the power spectrum, thereby obtaining a wavelet power spectrum,

    C_l^wavelet(j) = ⟨ |κ * ψ^j,θ|^2 ⟩_x,y,θ = ⟨ I_1^2 ⟩_x,y,θ

The wavelet power spectrum C_l^wavelet(j) is the averaged power spectrum C_l within the frequency coverage of the wavelets[Strictly speaking, the averaging kernel is ⟨∫(ψ^j,θ(k,θ_k))^2 k dθ_k ⟩_θ, while Figure <ref> shows ψ^j,θ(k,θ_k=0). More discussion can be found in the appendix C of <cit.>.]. For brevity, we shall use C_l to represent C_l^wavelet(j) in figures and tables and specify the j-range in the captions. 
Compared to the normal binning of power spectrum, the wavelet binning guarantees that the scale coverage of the power spectrum is exactly the same as the scattering transform. We also compute the normally binned power spectrum, where 8 bins are set in the range 300<l<3500 with equal spacing in logarithmic of l. They are computed directly from the κ map without further treatment on masks. As we follow a forward modeling approach and compute the statistics for observational and mock data in the same way, the mask should not bias the inference.  
When computing the scattering transform and power spectra, we adopt the flat sky approximation and assume periodic boundary condition.
To reduce the influence of masks, when calculating the scattering transform s_1, s_2 and wavelet power spectrum C_l(j), we exclude the masked regions in the last step of spatial average, i.e., ⟨·⟩_x,y becomes a weighted average with the mask as its weight.




















 §.§ Mock data and emulator


We follow a simulation-based inference approach to infer cosmological parameters from the observational data, as there is no analytical predictions for the scattering transform statistics so far. To build an emulator for the summary statistics, we use two sets of N-body simulations <cit.> to generate multiple realizations of mock shear catalogs customized for our observational data (HSC Y1), and then make mass maps following the same pipeline as on real data. 

Detailed description of the procedures to make the mock catalogs, the cosmological inference framework, and inclusion of systematics are presented in Appendix <ref>, and details about the N-body simulations are presented in Appendix <ref>. Below we only briefly summarize the information. 

In our mock dataset, the N-body simulations assume the ΛCDM model. One set is used to sample cosmic variance at a fiducial cosmology <cit.>, and the other is used to explore the cosmological dependence of Ω_m and S_8 <cit.>, which samples 100 cosmology in the range of 0.1 < Ω_m < 0.7 and 0.23 < S_8 < 1.1 . After the mock mass maps are made, we build a likelihood emulator of the summary statistics to perform cosmological inference. We fit the simulations with Gaussian sampling distribution whose mean vector changes with Ω_m and S_8, but covariance matrix is independent of cosmological parameters.

We also include the intrinsic alignment strength A_IA as a nuisance parameter to be inferred, and we marginalize over two systematic effects: multiplicative bias and PSF residual. For baryonic feedback and different photometric redshift, we only test their influence on different prescriptions but not marginalizing them. We do not include source clustering in our mock data as it is not expected to bias our statistics, which are even functions of the field values <cit.>. 








 §.§ Non-Gaussian tomography



Correlated Gaussian information between different tomographic maps can be exhausted by computing all the 4 auto- and 6 cross-power spectra. However, when considering non-Gaussian statistics, because many non-Gaussian statistics are defined on single maps, the analysis has usually been performed on individual tomographic maps, and there has not been an agreed way to extract correlated non-Gaussian information. In the context of peak count statistics, <cit.> proposed to include the cross information by simply computing non-Gaussian statistics on joint maps,

    κ_joint   = ∑_i N_iκ_i

where N is the effective galaxy number density for a tomographic map and the summation is over subsets of 1,2,3,4. There are 11 such combinations plus the 4 individual maps. 
We point out that this is also similar to the way neural networks deal with multi-channel images, and it can be interesting to consider not only the summation, but more general linear combinations, and further look for the optimal choices. Nevertheless, in this study we just follow the idea of <cit.> and compute the scattering transform on the 15 maps. 
We find that including the cross-map information significantly improves the constraining power 

on Ω_m and A_IA. Figure <ref> shows the mock inference for Ω_m and S_8 with no intrinsic alignment. Compared to the usual strategy of using only 4 tomographic maps or only 1 map with galaxies of all redshifts, a factor of two is earned by including the joint maps. In fact, when assuming A_IA=0, it suffices to add only 3 joint maps, κ_1+2+3+4, κ_2+3+4, and κ_3+4 to the 4 individual maps κ_1, κ_2, κ_3, and κ_4, which provides essentially the same constraining power as using all 15 maps. However, when marginalizing over A_IA, including all the 15 maps improves the constraint of Ω_m by about 20%. A quantitative comparison using real HSC Y1 data is shown in Table <ref>.












§ RESULTS AND DISCUSSIONS


Figure <ref> shows the cosmological constraints on HSC Y1 data. 
As the default setting, we use maps smoothed by a 2 arcmin Gaussian kernel and wavelets with j= 2,3,4 corresponding to central frequencies at l = 600, 1200, 2400 and together covering l≈ 300–3500 (see Figure <ref>).  In total 90 scattering coefficients (45 s_1 and 45 s_2) are computed from 15 tomographic maps and joint maps. The wavelet power spectra is computed using the same wavelets for frequency binning and on 10 tomographic and joint maps, which is enough to cover all Gaussian information.

The posteriors of Ω_m and S_8 are obtained with the intrinsic alignment parameter A_IA and other nuisance parameters marginalized.  The constraints are also listed in Table <ref>. Our tight constraint obtained from the scattering transform and power spectra, Ω_m=0.287_-0.030^+0.042 and S_8≡σ_8(Ω_m/0.3)^0.5=0.832±0.018, is consistent with the value inferred by <cit.> (TT,TE,EE+lowE+lensing) in both parameters. The small error bar of Ω_m and its consistency with Planck demonstrate the power and reliability of the scattering transform statistics.

In Table <ref>, we also show the constraints with several adjustments to the default setting to demonstrate the robustness of our result. In particular, we find that an Illustris-TNG baryonic feedback only shifts the S_8 inference up by 0.004 (0.2σ). Further study of the impact of baryons can be found in <cit.>.  Different photo-z codes or PSF-residual modeling do not shift S_8 significantly, either.  







 §.§ Non-Gaussian estimate improvements


By adding the scattering transform coefficients to the wavelet power spectra, the constraint on Ω_m is improved by a factor of 2[When compared to the power spectra with 6 bins between l = 300–1900, as adopted by <cit.>, the improvement is a factor of 3, but apparently, the additional factor of 1.5 comes from extending the frequency range rather than including non-Gaussian information.].
<cit.> have shown that on a single-redshift convergence map, the constraining power of the scattering transform is on par with that of the convolutional neural networks reported in <cit.>. Here we confirm the same constraining power on tomographic analysis of real data by comparing to <cit.>'s convolutional neural net (CNN) results with the same data set.

Recently, new types of neural-net-based models have also been proposed for cosmological inference <cit.>.  In particular, the multiscale flow model combines the wavelet decomposition idea as used in the scattering transform with the normalizing flow idea to model the likelihood of the whole field, and has been reported to achieve better constraining power than CNN and the scattering transform on simulated weak lensing data <cit.>.  However, its performance in real data with source clustering, survey mask, and other systematics remains to be confirmed.



As for the constraining power of S_8, we find on both mocks and real data no significant improvement from including non-Gaussian information, consistent with previous forecast on mock data <cit.>. There has been some confusion in the literature about the tightening of S_8 constraint from non-Gaussian statistics. We find that almost all the reported tightening of S_8 in the literature can be explained by two reasons. One is the aforementioned mismatch of scales between the non-Gaussian and Gaussian statistics. Many non-Gaussian statistics rely on map smoothing to cut scales, but a Gaussian smoothing kernel does not draw a sharp frequency cut (Figure <ref>). The second reason is leakage from the tightening of Ω_m due to remaining correlation between S_8 and Ω_m, which can be eliminated by optimizing the power index a in the definition of S_8≡σ_8 (Ω_m/0.3)^a or by fixing (conditioning on) Ω_m. Therefore, we conclude that there is no evidence that non-Gaussian statistics can significantly tighten the error bar of S_8. Nevertheless, non-Gaussian statistics do provide consistency checks for the S_8 constraint.











 §.§ Influence of photometric redshift on S_8



In this section we first point out an internal tension of S_8 estimates between different redshift bins, suggesting that photometric redshift is a main limiting factor for an accurate estimate of S_8. We then show that when adopting a similar treatment of photometric redshift, our result is consistent with other major results of HSC data analyses.  However, changing the treatment of photometric redshift will substantially change the inferred value of S_8.

Figure <ref> shows an internal consistency check of photo-z calibration using the inference from different subsets of tomographic maps. For example, the label `12 4' means that all the 6 auto- and cross-power spectra without bin 3, where bin 1–4 correspond to photo-z ranges of [0.3–0.6], [0.6–0.9], [0.9–1.2], and [1.2–1.5]. A systematic trend is found that as long as the 3rd redshift bin is absent, the inferred S_8 value is lower than those with bin 3 by 0.03.  Similar patterns are also found when using the scattering transform or binned power spectra as summary statistics, in addition to the wavelet power spectra inference shown in Figure <ref>, so this pattern should not be caused by particular statistics but rather linked to the galaxy catalog or map itself.

Though such a lowering of S_8 (1.5σ) may be caused by statistical fluctuations, it may also be caused by an underestimate of the density of high redshift sources, i.e. a tail in the redshift distribution of the 3rd photo-z bin.  Indeed, there are multiple challenges to obtaining reliable photometric redshift estimates when going to deeper magnitude and higher redshifts <cit.>.  It requires modeling not only the spectra of galaxies within the redshift range of interest (0.3–1.5 in our case) but all galaxies, and the modeling error for any other galaxies can cause contamination into the redshift range of interest.  
The sensitivity comes from sharp features, mainly the Balmer and Lyman alpha breaks, in galaxy spectra.  
However, the narrow range of optical bands and photometric error lead to degenerate photo-z solutions. As a result, the underlying true redshift distribution of galaxies selected within our photo-z bins tends to have a long tail towards high redshifts, making it difficult to robustly estimate the redshift distribution of the population.  In addition, galaxies above z = 2 in general have high star formation rates and rich emission lines, making it difficult to build reliable templates for galaxy spectra.


In fact, a similar photo-z tension is also found for the 3rd and 4th redshift bins in the HSC Year 3 (Y3) data with 3 times more galaxies <cit.>.  As shown in Figure <ref>, the inferred S_8 values also cluster into two groups: with the raw photo-z used, values of S_8=0.823^+0.032_-0.028, 0.824±0.021, 0.813±0.024 are obtained by the correlation analysis of HSC Y1 data <cit.>,  correlation analysis of HSC Y3 data <cit.>, and pseudo-C_l analysis of HSC Y3 <cit.>. All of the three are consistent with each other and our result[When comparing their results to ours, we note that for photo-z prescription, they all use the COSMOS re-weighted redshift distribution from galaxies, while we use the stacked redshift posterior with the  code. According to the tests presented in table 7 of <cit.> and figure 7 of <cit.>, the inferred S_8 is higher by about 0.02 when switching from their method to ours. Therefore, before any special adjustment on photo-z, most HSC analyses, including our simulation-based inference, provide almost identical high values of S_8 at 0.82. The only analysis giving a low S_8 without photo-z adjustment is the pseudo-C_l analysis of HSC Y1 data reported in <cit.> (S_8 = 0.78). However, they also obtain an unusually low Ω_m=0.16 not shared by the other analyses, which casts doubts on the reliability of their analysis.] at a `high value' of S_8.  


On the other hand, when a photo-z adjustment is performed, a similar lowering of S_8 is also found in the other HSC analyses. Instead of just excluding bin 3 as we did, the larger area of Y3 data allows for a simultaneous inference (self calibration) of photo-z errors. <cit.> and <cit.> assume that the shape of true redshift distribution in bin 3 and 4 are the same as that provided by the photo-z code, but an unknown overall offset is to be fitted.  It turns out that the HSC Y3 data favors an uncomfortably large offset of Δ z ∼ 0.1 and 0.2 in the two tomographic bins, which dramatically lowers the inferred S_8 value: S_8=0.769^+0.031_-0.034 and 0.776±0.032 are obtained from the correlation <cit.> and power spectrum <cit.> analyses, consistent with the lowering of S_8 we found in HSC Y1 data when excluding galaxies with photo-z between 0.9 and 1.2. 
Moreover, this internal tension is also observed with non-Gaussian statistics, including in our scattering transform analysis and the peak count analysis of the HSC Y1 data <cit.>, despite the different choices of tomography and scale cut.





 §.§ Comparison to other S_8 estimates


Figure <ref> shows the comparison of our S_8 inference with others in the literature. Regardless of which statistic is used, our results based on the raw redshift estimate from the HSC pipeline <cit.> have a `high' value of S_8 consistent with Planck. 
As discussed in section <ref>, if we remove background sources with photometric redshifts in the 0.9–1.2 range (based on their "best point estimate" using the   code), the S_8 estimates become relatively lower and more consistent with those obtained with the KiDS[<https://kids.strw.leidenuniv.nl/>] and DES[<https://www.darkenergysurvey.org/>]
weak lensing surveys. 
Because of this internal redshift tension, it is currently difficult to robustly conclude on which value of S_8 is correct. 

For a deep survey like the HSC, a significant fraction of galaxies have relatively high redshifts. Therefore, simply removing them from the analysis causes substantial loss of constraining power.
The self-calibration method adopted by <cit.> and <cit.> is a practical solution but still not ideal, as the uncertainty of the redshift distribution is likely from the tail part which is not well captured by the simple overall shift.  To definitively identify the nature of this internal tension from photo-z estimate, it may be essential to obtain independent redshift estimates, for example using the clustering redshift technique <cit.>.

The above internal tension of S_8 suggests that photo-z calibration is a main limiting factor for S_8 inference in HSC and other ground based surveys reaching a similar depth, including the Rubin Observatory LSST[<https://www.lsst.org/>] <cit.>. Before robust photo-z estimates become available, improvement of S_8 is unlikely to come from optimizing the choice of spatial statistic or machine learning tools [Nevertheless, the inference of Ω_m using non-Gaussian statistics seems robust to the potential photo-z bias.]. As the photo-z bias is likely caused by the ambiguity of spectral structures such as the Balmer and Lyman alpha breaks, space-based surveys such as Euclid[<https://www.cosmos.esa.int/web/euclid>] <cit.>, Xuntian space telescope survey<cit.>, and Roman space telescope survey[<https://roman.gsfc.nasa.gov/>] <cit.> which extend photometry to near infrared may significantly alleviate the photo-z problem and therefore draw conclusion on the S_8 tension.








§ CONCLUSION


We present the first application of the scattering transform to weak lensing observations. The scattering transform is a summary statistic that borrows ideas from convolutional neural nets, but requires no training and is compact and interpretable (see  for more details). To best demonstrate its power with real data, we choose to analyze the HSC Year 1 data <cit.>, which is the deepest wide-field weak lensing survey data available to date.

From a full tomographic analysis of the HSC Y1 data, we conclude:


    
  * Spatial statistics: We find that the scattering transform significantly improves the constraint on Ω_m. It tightens the error bar by a factor of 2 when compared to Gaussian statistics in the same scale range, reaching the constraining power of convolutional neural nets <cit.>. We also find that to improve Ω_m constraint, the non-Gaussian information across tomographic maps cannot be ignored, and it can be captured by measuring non-Gaussian statistics on linear combinations of different tomographic maps.

    
  * Photo-z calibration: We find that the inference of S_8 is more limited by photo-z calibration than choosing spatial statistics. When excluding galaxies with photo-z in [0.9-1.2] (the 3rd bin in HSC) from the analysis, the inferred S_8 changes from 0.82 to 0.79 independent of which statistic is used, suggesting an under-estimation of redshift of those galaxies.  A similar shift is also found in the analyses of HSC Y3 data <cit.>. 
    The S_8 estimation is currently limited by photometric redshift estimation. This photo-z problem is likely shared by all ground-based weak lensing surveys reaching a similar depth. Calibration with clustering redshift and infrared photometry from space-based telescopes may help to overcome this problem and provide final answer to the S_8 tension.














§ ACKNOWLEDGEMENTS

We thank Joachim Harnois-Déraps for providing the simulated intrinsic alignment maps and Ken Osato for providing the kappa-TNG maps.
We thank Alex Drlica-Wagner, Marco Gatti, Kevin M. Huffenberger, Nickolas Kokron, Bhuvnesh Jain, and Masahiro Takada for useful discussions.
S. C. thanks Siyu Yao for her constant inspiration and encouragement.
S. C. acknowledges the support of the Martin A. and Helen Chooljian Member Fund, fund from Zurich Insurance Company and the Fund for Natural Sciences at the Institute for Advanced Study.
This work was supported by JSPS KAKENHI Grants 19K14767 and 20H05861 (to MS) and 23K13095 and 23H00107 (to JL). We acknowledge support from the MEXT Program for Promoting Researches on the Supercomputer Fugaku hp230202 (to JL). 
This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. We thank the Yukawa Institute for Theoretical Physics at Kyoto University, where part of the discussions of this work took place during YITP-T-22-03 on “Cosmology with Weak Lensing: Beyond the Two-point Statistics” and YITP-W-23-02 on “Future Science with CMB x LSS”. Part of this work was performed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. This research used computing resources at Kavli IPMU. This research used resources at the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231. 



§ DATA AVAILABILITY



The derived data generated in this research will be shared on reasonable request to the corresponding author.

mnras







§ SMOOTHING OF MASS MAP
 

With smoothing, the weighted average should be extended from within each pixel to galaxies in its neighborhood, therefore the weighting becomes a combination of the weights of individual galaxies, w_i and a smoothing kernel W(Δ r). Approximating the positions of each galaxy by the center of the pixels where they are located, a smooth ϕ̅ can be obtained by smoothing individually the numerator maps and the denominator map in eq. <ref>,

    ϕ(x,y)_smooth≡W(Δ r) * ∑_in pixel (x,y) w_i ϕ_i/W(Δ r) * ∑_in pixel (x,y) w_i ,
 
where `*' represents convolution on a grid, and the smoothed fields ϕ(x,y)_smooth for the four quantities {e, ℛ, m_tot, c} are used to generate the γ̂ map through eq. <ref>. We use the Gaussian smoothing kernel with standard deviation σ_r = 2 arcmin:

    W(Δ r) = 1/2πσ_r^2exp(-Δ r ^2/2σ_r^2) .

In the literature, as little variation is expected in the map of m_tot and ℛ in eq. <ref>, they are usually treated as constants over a field of view and redshift bin. We have confirmed that replacing m̅_tot and ℛ̅ maps by the global weighted averages over all galaxies only results in negligible change in the cosmological inference. Therefore, our map making method is essentially the same as <cit.>. 

Below we discuss the advantage of this smoothing method compared to a direct convolution of ϕ̅(x,y), as adopted in <cit.>.
Compared to eq. <ref> where both the numerator and denominator are convolved by the smoothing kernel, an alternative is to directly convolve the pixelized map by the smoothing kernel,

    ϕ'(x,y)_smooth≡ W(Δ r) * 
     ∑_in pixel (x,y) w_i ϕ_i/∑_in pixel (x,y) w_i ,
 
However, such a direct convolution is not an optimal estimator because it weighs pixels equally while the effective number of galaxies in each pixel is different. Thus it increases the noise level of the map and reduces the constraining power on cosmological parameters. This is a main difference between the maps used in <cit.> and this paper.
To give a quantitative comparison, on the non-tomographic mass map at fiducial cosmology smoothed by a 1 arcmin kernel, our weighted smoothing method provides a correlation coefficient of 0.38 (an S/N of 0.37) between a noisy mock and the noiseless convergence map, whereas the equal-weight smoothing gives a correlation coefficient of 0.32 (an S/N of 0.22). The error bar of inference increases by 50% in Ω_m and 25% in S_8≡σ_8 (Ω_m/0.3)^0.5 on both mock and real data. In power spectrum and correlation function analyses <cit.>, although there is no explicit smoothing in the map making step, a treatment of weighting equivalent to eq. <ref> is included in the formalism of computing their statistics.





§ MORLET WAVELET



The profile of the Morlet wavelets we use is given by

    ψ^j,θ(x⃗)   =1/σexp( -x^2/2σ^2) [exp(ik⃗_⃗0⃗·x⃗) - β] ,
    ψ^j,θ(k⃗)    = exp[-(k⃗-k⃗_⃗0⃗)^2σ^2/2] - βexp(-k^2σ^2/2) ,

where · represents Fourier transform, and 
    σ=0.8× 2^j ,  |k_0| = 0.75 π× 2^-j , angle(k_0) = πθ / L

specify the dependence on wavelet size index j and orientation index θ. The offset β=exp(-k_0^2σ^2/2) is introduced to guarantee that the wavelets are strictly band-pass instead of low-pass filters: ψ(0)=0. More explanations of the mathematical form can be found in the appendix B of <cit.>.

As a result, the central frequency of the j=0 wavelet kernel is defined to be 3/8 pixel^-1, and when j increases by 1, the wavelet becomes larger in real space by a factor of 2. In this study we use wavelets with j = 2–4. The j=2 wavelet roughly covers l ≈ 1200–3500, 

j=3 wavelet covers l ≈ 600–1800,
and j=4 wavelet covers l ≈ 300–900. 
In principle the sampling of j does not have to be integers, nevertheless we still sample j at integers as we find no significant improvement of cosmological constraint from finer sampling. As for orientations, we follow previous studies and set the total number of orientations to be L = 4, i.e., sampling every 45 degrees.











§ SIMULATION BASED INFERENCE METHOD





 §.§ Mock catalogs


The mock catalogs are generated following <cit.> by first randomly rotating the observed galaxies and statistically separating the measurement noise from intrinsic galaxy shape,

    e^int   = random rotate(e_rms/√(e_rms^2 + σ_e^2)) e^obs
    e^mea   = σ_e N


where e_rms and σ_e are read from the observed galaxy catalog <cit.> and N is a random variable independently sampled for each galaxy from a standard Gaussian distribution with zero mean and unity variance 

Then, the mock ellipticity catalog is obtained by applying cosmic shear distortion δ^sh to the intrinsic galaxy shapes and then adding mock measurement noise <cit.>[With some algebra, one can rewrite these equations in a compact form, which can be found in the literature, in terms of the complex-valued reduced shear g ≡γ/1-κ and complex ellipticity: e = e^int + 2g + g^2 e^int*/1 + |g|^2 + 2Re(g e^int*), by using the relation δ^sh=2g/1+|g|^2 and a=1+|g|^2/2.]

    e_1    = e_1^int + δ^sh_1 + a δ^sh_2 (δ^sh_1 e_2^int - δ^sh_2 e_1^int)/1 + δ^sh·e^int + e_1^mea
    
        e_2    = e_2^int + δ^sh_2 + a δ^sh_1 (δ^sh_2 e_1^int - δ^sh_1 e_2^int)/1 + δ^sh·e^int + e_2^mea
    
        a    ≡1 - √(1-(δ^sh)^2)/(δ^sh)^2 = 1/1+√(1-(δ^sh)^2)
    δ^sh   ≡2(1-κ)/(1-κ)^2+γ^2γ ,

where a is a factor close to 0.5 at the small shear limit, and the shear distortion δ^sh is obtained by ray-tracing the N-body simulations. We also simulate the effect of multiplicative bias m_tot by replacing γ with (1+m_tot)γ.

The fiducial simulation set <cit.> includes 108 quasi-independent full-sky realizations at one fiducial cosmology with the ΛCDM model and WMAP9 parameters <cit.>: Ω_cdm = 0.233, Ω_b = 0.046, Ω_m = Ω_cdm + Ω_b = 0.279, Ω_Λ = 0.721, h = 0.7, σ_8 = 0.82, and n_s = 0.97. The sky footprint of the HSC fields was rotated in 21 angles to take full advantage of the large volume of each simulation, resulting in a total of 2,268 realizations of mock HSC lensing maps. More details are described in Appendix <ref> and the references. This set of mock data at the fiducial cosmology is used to estimate the cosmic variance. 

The cosmic-varying simulation set in <cit.> samples 100 different cosmologies with 0.1<Ω_m<0.7 and 0.23<S_8<1.1. The other parameters are consistent with the Planck 2015 results <cit.>: Ω_b = 0.049, Ω_m+Ω_Λ = 1, h = 0.6727, and n_s = 0.9645. We make 50 realizations per cosmology by random arrangement of the simulation boxes and further enlarge the size to 150 by adding different realizations of galaxy noise. 

To obtain the κ and γ in eq. <ref> to enable generation of the mock galaxy shape catalog, one needs to assign a redshift to each of the galaxies. We follow the same method detailed in <cit.>. In brief, we use the stacked redshift posterior of galaxies provided by the  code, a photo-z code using a self-organizing map, in each tomographic bin as the assumed redshift distribution.




 §.§ Likelihood emulator


We parameterize the likelihood function by assuming that at each cosmology, the summary statistics s⃗ follow a multivariate Gaussian distribution under cosmic variance and galaxy shape noise: 

    p(s⃗|θ) = exp[-1/2 (s⃗-μ⃗(θ))^T  C_fid^-1 (s⃗-μ⃗(θ))]/(2π | C_fid|)^dim(s⃗)/2 .

Usually, this assumption is justified by the fact that the statistical estimator is an average over a large field of view, for which the central limit theorem Gaussianizes its distribution. This assumption is particularly well justified for the scattering transform coefficients because the modulus operation avoids amplification of outliers and thus accelerates the Gaussianization <cit.>.

We estimate the mean μ (θ) as a function of cosmological parameters θ = (Ω_m, S_8) from the set of cosmic varied simulations, by first calculating the sample mean at each cosmology μ̂ (θ_i) and then fit it with 6th-order polynomials of (logΩ_m, log S_8).

    μ⃗ (Ω_m, S_8) = ∑_0 ≤ n+m ≤ 6a⃗_n,m (logΩ_m)^n (log S_8)^m ,

where n and m are natural numbers. The reason to use logarithms is that μ(Ω_m, S_8) or its derivatives may have singularity at Ω_m=0 or S_8=0. As the summary statistics may change quickly when Ω_m approaches zero, we use a logarithm to regularize this behavior and make the polynomial regression converge faster. The regression is implemented through least square method for each component of a⃗_n,m. 
We verify the robustness of our emulator by fitting it with 5th-order and 7th-order polynomials and find no significant difference, which means our choice does not under-fit or over-fit. The covariance matrix C_fid in our model is fixed and does not change with cosmological parameters. We estimate it only with the 2,268 fiducial realizations that originate from 108 quasi-independent simulations and a full-sky setting. In addition, we find small differences between the fiducial mocks and the emulator prediction based on the cosmo-varied simulations, which is likely caused by the small and adaptive box sizes used in the cosmo-varied simulations. We calibrate the emulator using the mean value of fiducial mocks, which is described in detail in appendix <ref>.





 §.§ Systematics


To account for observational and astrophysical systematics, we include the following parameters in our emulator: multiplicative calibration error Δ m, galaxy intrinsic alignment strength A_IA, and baryonic feedback strength A_tng. We further model their influence to the summary statistics  with a simple form,

    μ⃗(Ω_m, S_8, A_IA, Δ m, A_tng) =  b⃗^A_IAc⃗^A_tngd⃗^Δ mμ⃗(Ω_m, S_8, 0,0,0) ,

where the coefficients b⃗, c⃗, d⃗ are constants estimated using simulations described below and the operations are meant as element-wise. That parameterization is equivalent to a linear model for the logarithm of our statistics,

    log   μ⃗(Ω_m, S_8, A_IA, Δ m, A_tng) = 
       A_IAlogb⃗ + A_tnglogc⃗ + Δ m logd⃗ + logμ⃗(Ω_m, S_8, 0,0,0) .

We consider it better than a linear model of the original statistics because both the scattering coefficients and the power spectrum are always positive and related to fluctuation strength. This relation is accurate in the neighborhood of the fiducial cosmology, though its accuracy may gradually degrade as one moves away from the fiducial.




  §.§.§ Intrinsic alignment


To incorporate intrinsic alignment to our emulator for non-Gaussian statistics, we use a mock catalog provide in <cit.> which simulates this effect. This catalog was made following a prescription called δ-NLA, which computes the intrinsic alignment field e^IA according to the non-linear tidal alignment model (NLA) <cit.> but also samples galaxies naturally according to the matter density field δ. Therefore, compared to the NLA model, it realistically up-weights regions with higher matter density. The intrinsic alignment strength is still parametrized by A_IA as in NLA. For simplicity we assume no evolution of A_IA with redshift. 






Previous studies with HSC Y1 data obtains consistent constraints with different error bar sizes, A_IA=0.4±0.7 <cit.> and A_IA=0.9±0.3 <cit.>, with the NLA model.
According to those approximate range of A_IA, we build mocks with A_IA = 0 and 1.5, and determine the coefficient b⃗ in eq. <ref> by b⃗ = ( ⟨s⃗⟩_A_IA=1.5/⟨s⃗⟩_A_IA=0) ^ 1/1.5 .
The average is measured from 400 mocks with 10×10 deg^2 field of view, generated by adding different galaxy shape noise realizations to the 10 mock catalogs provided in <cit.>.




  §.§.§ Error in multiplicative bias


To account for the effect of mis-calibration of multiplicative bias, we model its influence on summary statistics and then marginalize over it. A prior of a Gaussian with 0.01 standard deviation is applied, which is the level of calibration residual estimated from image simulations <cit.>. We determine the coefficient c in eq. <ref> by its value at fiducial cosmology d⃗ = (⟨s⃗⟩_fid, Δ m=+0.01/⟨s⃗⟩_fid, Δ m=-0.01)^1/0.02, which are obtained by adding Δ m = +0.01 and –0.01 to the multiplicative bias of mock galaxies in the 2,268 realizations of fiducial mock catalogs, while keeping all other ingredients including the random seed for galaxy shape noise unchanged. 
We find that the effect of multiplicative bias error basically shifts S_8 by the same fraction, while not changing the inference of Ω_m. Data itself is not able to constrain Δ m and it is mainly driven by the prior. 








  §.§.§ Baryonic feedback

We use the kappa-TNG mocks <cit.> to test baryonic effects. They are a set of convergence maps made from the hydrodynamic IllustrisTNG simulations <cit.>, with 10,000 realizations of 5×5 deg^2 maps ray-traced from the highest resolution of the largest box TNG300-1 and its corresponding dark-matter-only counterpart TNG300-1-Dark. Those simulations adopts a cosmology consistent with <cit.>.
Since these mocks are not tailored for the HSC Y1 field of view, the map statistics are not defined identically to the emulator or observation due to mask and edge effects. 
The coefficient c⃗ in eq. <ref> is then determined by the ratio between hydro and dark matter only simulations: c⃗ = ⟨s⃗⟩_TNG/⟨s⃗⟩_TNG-dark.



  §.§.§ Photometric redshift algorithms


To test the influence of photometric redshift error, we generate 210 new mocks at the fiducial cosmology using photo-z from two other codes,  and  <cit.>. Due to the limited number of additional mocks, we are not able to incorporate Δ z_phot as a nuisance parameter 

and marginalize over it.  Therefore, we only apply the ratio ⟨s⃗⟩_fid, / ⟨s⃗⟩_fid, and ⟨s⃗⟩_fid, /⟨s⃗⟩_fid, to the emulator and examine the shift of cosmological inference.





  §.§.§ PSF residual


We also include PSF systematics in our model. Errors in the measurement of point spread function (PSF) and the PSF leakage from imperfect deconvolution lead to additive bias of the shear measurement. The residual of this calibration acts like correlated noises added to the shear map. To account for its effect, we simulate this noise in our mock maps by adding a Gaussian random field to the mock convergence maps, with the four redshift bins sharing the same PSF systematics. We set the power spectrum of this noise using Figure 7 of <cit.>, which corresponds to setting α_psf=0.03 and β_psf=-0.89 
as measured from data by <cit.>. For simplicity we omit the uncertainty of this power spectrum. Assuming no PSF residual typically shifts S_8 up by 0.004 and Ω_m up by 0.01. A more complicated treatment such as a marginalization of the parametrized power spectrum of the PSF systematics is left for future work.





  §.§.§ Source clustering


We have not modeled the effect of source clustering, but according to an exploration performed by <cit.>, we expect negligible biased for the scattering transform and power spectrum.  Source clustering is the fact that background galaxies are not uniformly distributed but clustered with some correlation of the foreground. Therefore, the observed lensing map is sampled in a biased way and have different statistical properties from uniformly sampled ones as used in our simulations.  Nevertheless, source clustering mostly influences statistics that are odd function of the field (i.e., statistics that change their sign when the field is flipped) <cit.>.  As a result, the power spectrum and the scattering transform are not much influenced but source clustering, whereas third-order moment based statistics, phase harmonic coefficients that are odd, and machine-learning methods may be significantly biased if the source clustering effect is not properly modeled. 






§ N-BODY SIMULATIONS


The fiducial simulation set <cit.> includes 108 quasi-independent full-sky realizations at one fiducial cosmology with the ΛCDM model and WMAP9 parameters <cit.>: Ω_cdm = 0.233, Ω_b = 0.046, Ω_m = Ω_cdm + Ω_b = 0.279, Ω_Λ = 0.721, h = 0.7, σ_8 = 0.82, and n_s = 0.97. 
In total 6 independent runs, each with 14 simulations with different box sizes (450–6300h^-1 Mpc), were implemented to keep the angular resolution roughly the same along the line of sight. Each box has 2048^3 particles corresponding to an angular particle grid size of 3.4–6.8 arcmin at the initial condition. The particles in each box are also randomly shifted and ray-traced to create 108 quasi-independent full sky maps with HEALPix resolution N_side=8192. It has been tested that for lensing convergence with sources at z = 1, the power spectrum in simulation is consistent with halo-fit results within 5% till around l=3000 (see Fig. 23 in <cit.> and eq. 26 in <cit.>).

The cosmic-varying simulation set in <cit.> samples 100 different cosmologies with 0.1<Ω_m<0.7 and 0.23<S_8<1.1. The other parameters are consistent with the Planck 2015 results <cit.>: Ω_b = 0.049, Ω_m+Ω_Λ = 1, h = 0.6727, and n_s = 0.9645.  At each cosmology, four N-body simulations with different box sizes and resolutions were run and connected along the light cone to cover a 10 × 10 deg^2 field of view. The box sizes are designed to match the opening angle at all redshifts, with typical values at Ω_m = 0.3 to be 280–930 Mpc/h. 
The projected particle grid at the initial condition in these simulations has a size of 1.3–2.5 arcmin, which is 2.5 times finer than the fiducial set of simulations.
By randomly shifting each simulation box assuming periodic boundary condition, 50 cosmic shear realizations were generated for each cosmology. From stability tests such as using only a subset of the realizations to build the emulator, we find that 50 realizations are not enough to overcome the numerical noise of the cosmological dependence to mean data vector, so we enlarge the mocks to 150 realizations per cosmology by adding different galaxy noise. 


We find that the μ(θ) predicted by the emulator at the fiducial cosmology has an offset from the sample mean ⟨s⃗⟩ of the fiducial simulations. For all tomographic maps, our emulator based on cosmo-varied simulations predicts a lower power spectrum. The deviation is highest for low redshift and for large scales. The deviation in cosmic shear power spectrum (noise removed) is about 20% for the 1st tomographic bin and 10% for the 4th bin at l = 300. The deficit gradually changes to zero at around l = 2500 and becomes an excess after that. The excess is expected as the cosmo-varied simulations have higher resolution, but the cause of power deficit is unclear. This difference, dominated by the large scale power deficit, results in an offset in parameter inference. When using the mean of fiducial data vector as input for mock inference, we obtain a best fit of Ω_m=0.261 and S_8 = 0.808, 6.5% lower in Ω_m and 2% higher in S_8 than the fiducial values. We suspect that the deficit of power of our cosmo-varied maps based on <cit.> is caused by the adaptive box size of simulations along the light cone, as this is the main difference from another simulation based analysis of HSC year 1 data <cit.>, where no such effect is observed. The adaptive box size results in the fact that a large angular scale always corresponds to roughly the simulation box size, at which the power may be suppressed due to spectral leakage under windowing. For this study, we just correct for this effect by multiplying the emulator prediction of μ⃗ by a factor vector ⟨s⃗⟩ _fid / μ⃗(0.279, 0.791)_emulator. A full understanding of the deficit is left for future study.
















