{'fraction_non_alphanumeric': 0.03720068017729211, 'fraction_numerical': 0.02000111504474117, 'mean_word_length': 4.383582201545734, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Probing Cosmology and Cluster Astrophysics with Multi-Wavelength Surveys I. Correlation Statistics

Upcoming multi-wavelength astronomical surveys will soon discover all massive galaxy clusters and provide unprecedented constraints on cosmology and cluster astrophysics. In this paper, we investigate the constraining power of the multi-band cluster surveys, through a joint analysis of three observables associated with clusters of galaxies, including thermal Sunyaev-Zel'dovich (tSZ) effect in cosmic microwave background (CMB), X-ray emission of ionized gas, and gravitational weak lensing effect of background galaxies by the cluster's gravitational potential. We develop a theoretical framework to predict and interpret two-point correlation statistics among the three observables using a semi-analytic model of intracluster medium (ICM) and halo-based approach. In this work, we show that the auto- and cross-angular power spectra in tSZ, X-ray and lensing statistics from upcoming missions (eROSITA, CMB-S4, and LSST) can help break the degeneracy between cosmology and ICM physics. These correlation statistics are less sensitive to selection biases, and are able to probe ICM physics in distant, faint and small clusters that are otherwise difficult to be detected individually. We show that the correlation statistics are able to provide cosmological constraints comparable to the conventional cluster abundance measurements, while constraining cluster astrophysics at the same time. Our results indicate that the correlation statistics can significantly enhance the scientific returns of upcoming multi-wavelength cluster surveys.


INTRODUCTION
Upcoming multi-wavelength cosmological surveys will provide unprecedented images of the distribution of dark matter, gas and stars in the universe. Clusters of galaxies are the most massive gravitationally bound objects in the universe. They form at the highest density peaks in the large-scale structure of the universe. Their abundance traces the growth of structure of the universe, providing unique constraints on cosmology (e.g. Allen et al. 2011, for review).
With upcoming large-scale galaxy cluster surveys, astrophysical systematics such as feedback from SNe and AGN, non-thermal pressure support, and gas clumping, need to be controlled for the large samples of clusters in upcoming X-ray and microwave. In the X-ray band, the upcoming eROSITA mission will potentially improve cluster cosmological constraints by detecting very large number (∼ 10 5 ) of galaxy clusters out to redshifts z > 1 ( Merloni et al. 2012). In the microwave, ongoing ground-based experiments by the Atacama Cosmology Telescope (ACT) and the South Pole Telescope (SPT) will perform comprehensive studies of galaxy clusters through the tSZ effect (e.g. Hilton et al. 2018;Chown et al. 2018;Bocquet et al. 2019), while the nextgeneration ground-based experiment of CMB, the CMB-S4, will obtain precise measurements of the properties of ICM via tSZ effect and the mass distribution via CMB lensing (Abazajian et al. 2019). In the optical regime, ongoing galaxy imaging surveys will measure the lensing effect of large-scale structures with a few percent level accuracy (e.g. Hildebrandt et al. 2017;Troxel et al. 2018;Hikage et al. 2019), which will be further improved with future surveys with the Wide Field Infrared Survey Telescope (WFIRST 1 ), the Large Synoptic Survey Telescope (LSST 2 ), and the Euclid satellite 3 .
Multi-wavelength correlation statistics (i.e., correlation function or its Fourier-space equivalent, the power spectrum) are promising in providing robust statistical constraints on cluster astrophysics (e.g. Shaw et al. 2010;Battaglia et al. 2010) and cosmology (e.g. Komatsu & Seljak 2002;Shirasaki et al. 2016;Bolliet et al. 2018;Makiya et al. 2019). Correlation statistics are complementary to the conventional cluster abundance measurements, as they are less sensitive to selection biases, and are able to probe distant, faint and small clusters which will not be resolved on an individual-by-individual basis. In particular, given the large amount datasets coming online, cross-correlation or clusters across different wavebands will allow us to probe the baryonic and dark matter content in galaxy clusters that will improve our understanding of cluster astrophysics and cosmology (e.g. Van Waerbeke et al. 2014;Hill & Spergel 2014;Battaglia et al. 2015b;Ma et al. 2015;Hojjati et al. 2017;Osato et al. 2018;Lakey & Huffenberger 2019).
This paper is the first part of a series of papers where we investigate the use of multi-wavelength correlation statistics in constraining cosmology and cluster astrophysics. In this paper, we introduce a framework for jointly constraining cluster astrophysics and cosmology using multi-wavelength angular power spectra of galaxy clusters. Specifically, we use a semi-analytic model of ICM based from Shaw et al. (2010); Flender et al. (2017), which provide physical prescriptions for ICM profiles in dark matter halos, including gas physics relevant to the tSZ and X-ray observables, such as cluster cool-cores, non-thermal pressure. In this paper, we extend the previous model by incorporating the effects of gas clumping in cluster outskirts. We then provide Fisher forecasts on cluster astrophysics and cosmology from the multi-wavelength angular power spectrum analysis for upcoming surveys, in particular eROSITA in X-ray, CMB-S4 in microwave, and LSST in optical. Our forthcoming projects include validating our modeling of the correlation statistics with hydrodynamical simulations and applying our model to observational datasets. This paper is organized as follows. In §2, we provide an overview of galaxy cluster observables in X-ray, microwave, and lensing. We present our semi-analytic model of ICM in §3 and then describe a theoretical framework for computing two-point auto-and cross-correlations of the galaxy cluster observables in §4. In §5, we summarize the statistical uncertainties of two-point correlations from our Fisher analysis to forecast the parameter constraints of cosmology and ICM physics in upcoming surveys. §6 presents the main results in this paper, which includes (i) the expected parameter constraints in a hypothetical set of future multi-wavelength surveys with a large sky coverage of 20,000 square degrees, (ii) the dependence of survey configurations and analysis methods on our Fisher forecast, and (iii) possible systematic uncertainties induced by the statistical fluctuation in reconstruction method of tSZ effects and mis-estimation of lensing systematics. We discuss the current issues in our modeling of two-point correlation analysis in §7, and conclude this paper in §8.

OBSERVABLES
Throughout this paper, we assume General Relativity under a Friedmann Robertson Walker (FRW) metric with a spatially flat curvature. Within this framework, the governing equation of late-time expansion of the universe is given by where a is the scale factor, H(a) is known as the Hubble parameter and H 0 = 100h km s −1 Mpc −1 is the present value of H(a). Note that we work with the usual normalization a = 1 at present. The Ω m0 a −3 term in Eq.
(1) represents the contribution from cosmic matter density, while the second term in the right-hand side arises from dark energy causing the accelerated expansion of the universe at a > ∼ 0.5. The equation of state of dark energy is specified as w DE (a) = P DE (a)/ρ DE (a), where P DE and ρ DE denote pressure and density of dark energy, respectively. Ongoing and future surveys of large-scale structure of the universe aim at measuring w DE (a) at different epoch with high precision and exploring a possible deviation from the cosmological constant with w DE (a) = −1. In this sense, we adopt the simplest model of dark energy with constant value of w 0 , i.e. w DE (a) = w 0 , and study the constraining power of w 0 using several observables accessible to upcoming multi-wavelength astronomical surveys.

Thermal Sunyaev-Zel'dovich effect
The thermal Sunyaev-Zel'dovich (tSZ) effect probes the thermal pressure of hot electrons in galaxy clusters. At frequency ν, the change in CMB temperature by the tSZ effect is expressed as where T 0 = 2.725 K is the CMB temperature (Fixsen 2009), g(x) = xcoth(x/2) − 4 with x = h P ν/k B T 0 , h P and k B are the Planck constant and the Boltzmann constant, respectively 4 .
Compton parameter y is obtained as the integral of the electron pressure P e along a line of sight: where σ T is the Thomson cross section, χ is the comoving radial distance to redshift z, f K ( χ) is the comoving angular diameter distance, and χ H is the comoving distance up to z → ∞.

X-ray surface brightness
Another probe of ICM is X-ray surface brightness, defined as the number of X-ray photons per unit time per unit area per solid angle in a given range of energy E min to E max . The cumulative X-ray surface brightness from ionized gas along a line of sight is given by where X represents the X-ray emissivity for the energy band of interest, n e and n H ≈ n e /1.2 are the number density of electrons and protons, respectively. The X-ray emissivity depends on the temperature of plasma T and gas metalicity Z.
To be specific, X is computed as where Λ(T, Z, E) is the volume emissivity (the number of photons emitted per unit volume, time, and energy range) at energy E for a plasma of temperature T. The receiver's frame energy range is set as [E min , E max ] = [0.5, 2.0] keV. In this paper, we use the atomic database AtomDB 5 to compute Λ(T, Z, E) with Z = 0.3Z , where Z is the solar metallicity. We include Galactic absorption exp[−σ(E)N H ] where σ(E) is the effective photoelectric absorption cross section of hydrogen atom, taken from Morrison & McCammon (1983). The hydrogen column density is set to be N H = 2.5 × 10 20 cm −2 , which is estimated as the averaged N H for the full sky with the Galaxy masked 6 . Note that we define the surface brightness in units of photons s −1 cm −2 arcsec −2 throughout this paper.

Weak gravitational lensing in galaxy's image
Weak gravitational lensing offers a direct probe of the mass distribution of galaxy clusters. Relevant quantity of weak lensing effect in galaxy imaging surveys is known as convergence field κ, which is related to the second derivative of the gravitational potential Φ (Bartelmann & Schneider 2001). Using the Poisson equation and the Born approximation, one can express the weak lensing convergence field as the weighted integral of matter overdensity field δ m (x): where q( χ) is called lensing kernel. For a given redshift distribution of source galaxies, the lensing kernel is expressed as where p( χ) represents the redshift distribution of source galaxies normalized to ∫ d χ p( χ) = 1.

A SEMI-ANALYTIC MODEL OF INTRACLUSTER MEDIUM
To predict the statistical property of three fields defined by Eqs. (3), (4) and (6), we adopt a physically motivated semi-analytic ICM model as developed in Shaw et al. (2010);Flender et al. (2017). This ICM model is the modified version of the one originally proposed in Ostriker et al. (2005). The updated features from the original model include a radially-dependent non-thermal pressure component (e.g., turbulence or bulk flows) (Nelson et al. 2014b) and gascooling effect in cluster cores to provide a better description of recent X-ray data (Flender et al. 2017).

Dark matter distribution
The ICM model assumes that the gas inside cluster-sized dark matter halos initially follows the underlying dark matter density distributions which is modeled as Navarro-Frenk-White (NFW) profile (Navarro et al. 1996): where r s and ρ s are a scale radius and scale density, respectively. These parameters r s and ρ s are related to the concentration c ∆ = r ∆ /r s , through the definitions of the halo mass where ρ crit (z) is the critical density in the universe at redshift z and ∆ is the overdensity parameter. Eqs (9) and (10) reduce to Various overdensity parameters have been adopted in the literature. In this paper, we define the virial overdensity parameter as in Bryan & Norman (1998): where . Throughout this paper, we set the fiducial mass to be the virial mass with Eq. (12). For the concentration parameter, we adopt the model developed in Diemer & Kravtsov (2015, DK15 hereafter), which has been calibrated with dark-matter-only N-body simulations. The halo concentration can be affected by baryonic processes such as the star formation and additional AGN feedback (e.g. Gnedin et al. 2004;Duffy et al. 2010;Fedeli 2012;Shirasaki et al. 2018). We include baryonic effects in halo concentration by varying the normalization of halo concentration: where c vir = r vir /r s and c vir,DK15 is the predicted concentration by the DK15 model. We use the fitting formula from Hu & Kravtsov (2003) for the conversion between masses of different spherical overdensities.

Gas distribution
Given the gravitational potential determined by NFW profile, we then assume the gas inside every halo re-distributes so that its state can be described by hydrostatic equilibrium. The governing equation of hydrostatic equilibrium is given by where P tot is the total (thermal+non-thermal) pressure of gas, ρ g is the gas density, and Φ represents the NFW gravitational potential. We further assume a polytrope solution of Eq. (14), which is written as where n is the polytropic index and θ(r) is the polytropic variable In Eq. (17), Φ 0 is the NFW potential at the halo centre, and Γ = 1 + 1/n is the polytropic exponent. Cosmological hydrodynamical simulations suggest that the ICM follows Γ 1.2 except in high-density core (Ostriker et al. 2005;Shaw et al. 2010;Battaglia et al. 2012a). Gas cooling due to some strong feedback at the cluster core can change the polytropic exponent. The gas density in the inner region, relative to the critical density of the universe, is shown to be evolving with time (McDonald et al. 2013). Hence we adopt a two-value polytropic model, one for the outskirts, and one for the core: where r 500 is the spherical over-density radius with ∆ = 500 in Eq. (9) andΓ, γ are free parameters in our ICM model, and x break is fixed to be 0.195, which defines the cluster core region that is dominated by cooling and feedback physics. The two unknowns, i.e. the normalizations in gas pressure profile P 0 and gas density profile ρ 0 , are solved by two conditions. The first condition describes the energy balance between initial and final states of gas, given as where the left hand side of the equation is the final energy of the ICM, and the right hand side includes the initial energy of the ICM (denoted as E g,i ) and the work done by gas as it expands relative to its initial state (∆E p ). The two terms of DM |E DM | and f M * c 2 in Eq. (19) control the efficiency of two different feedback processes with two free parameters DM and f . The former is the energy introduced into the ICM during major halo mergers via dynamical friction heating, where E DM is the total energy in the dark matter halo (i.e., the sum of kinetic and potential energy). The latter is the energy injected into the ICM due to feedback from SNe and AGN, where M * is the stellar mass in halo. Throughout this paper, we parametrize the stellar-to-halo mass relation as where f * and S * are free parameters in the model. The second condition to set P 0 and ρ 0 is the boundary condition of gas pressure. The boundary condition is determined by setting the total gas pressure P tot at the boundary radius r f to the effective DM "pressure" P DM at the virial radius of the NFW halo scaled by the cosmic baryon fraction: where Ω b0 is the cosmic baryon density parameter and P DM is given by the velocity dispersion inside the NFW halo times the NFW density profile (e.g. see Shaw et al. 2010, for details). Once r f is set, we then determine ρ 0 and P 0 by solving together with the energy constraint equation in Eq. (19). We refer the reader to Ostriker et al. (2005) for details on setting P 0 and ρ 0 .

Non-thermal pressure
As in §2, the thermal pressure of ICM is closely related to the tSZ effect. The ICM model in this paper assumes a radiallydependent non-thermal pressure calibrated in hydrodynamic simulation (Nelson et al. 2014b) and thus the thermal pressure can be modeled by the residual between the total and non-thermal pressure profiles: where r 200m is the spherical over-density radius with respect to 200 times the mean matter density of the universe. In this paper, we assume the outermost radius of non-thermal pressure component to be r max = 4r 500 . The non-thermal pressure profile is independent of redshift when the radial scale of the profile is scaled in units of r 200m (z) (Nelson et al. 2014b). The non-thermal pressure in this model consists of only three free parameters, A nt , B nt , and γ nt . Note that the thermal electron pressure probed by the tSZ effect is given by P e = (2X + 2)/(5X + 3)P th , where X = 0.76 is the hydrogen mass fraction.

Gas density clumping
The X-ray emissions from the hot ICM (> 10 7 K) are dominated by thermal bremsstrahlung emission due to Couloumb interactions between electron and proton. The X-ray emissivity is proportional to the product of electron number density n e and proton number density n p : where Λ(T, Z) is the bolometric volume emissivity of the ICM, i.e. the number of X-ray photons emitted per volume dV for given gas temperature T and metallicity Z. In a fully ionized plasma with metallicty with 0.3 solar value, the electron density is related to the proton density as n e = 1.17n p . Thus, the X-ray surface brightness is proportional to the line-of-sight integral over X-ray emissivity or the gas density squared: x ∝ X ∝ n 2 e ∝ ρ 2 g . Gas density is expected to be highly inhomogeneous due to freshly accreted gas inside infalling galaxies and accretion from penetrating filaments (Nagai & Lau 2011;Vazza et al. 2013;Battaglia et al. 2015a). Failing to account for density inhomogeneity will lead to overestimates in gas density and mass by √ C, where C ≡ ρ 2 g / ρ g 2 1 is the clumping factor that quantifies the degree of gas density inhomogeneities. Gas density inhomogeneity is likely the physical reason for the apparent flattening of measured ICM entropy (e.g. Walker et al. 2013), as well as for the measured gas mass fraction that exceeded the expected cosmic baryon value (Simionescu et al. 2011).
We account of the effect of gas density clumping in our model by including the following generalized NFW form of gas density clumping as where C 0 , α C , β C and γ C are free parameters in the model. This model is motivated from the clumping factor profile computed from the Omega500 cosmological hydrodynamical simulations (Nelson et al. 2014a). Specifically, the parametrization in Eq. (26) provides a flexible fit to the clumping factor profile at r < ∼ 3r 200m , independent of the implemented baryonic physics, such as gas cooling, star formation and feedback from SNe and AGN. Scaling the radius by r 200m removes mass and redshift dependencies in the clumping profile, as r 200m is shown to be closely tracking the mass accretion rate of the halo (Lau et al. 2015), which governs the level of gas clumping in cluster outskirts. This leads to fewer parameters of our gas clumping model, compared to previous work (Battaglia et al. 2015a). Note that the gas clumping model includes gas inhomogeneities due to infalling substructures, such as merging halos and accreting filaments, as well as smaller scale gas density inhomogeneities due to turbulent motions in the ICM. We refer the reader to the companion paper (Lau et al. in prep.) for more details about our gas clumping model. The X-ray surface brightness profile in a single halo with mass M at redshift z is then given by where X,h , n e,h , and n H,h are the volume emissivity profile, the electron number density profile, and the proton density profile of an individual halo, respectively. We use this X-ray surface brightness profile to take into account the effect of gas clumping on the correlation statistics in §4. Table 1 summarizes the 15 free parameters of the ICM model used in the paper and their fiducial values. For the normalization of halo concentration, we take A C = 1 as a baseline model which is consistent with the dark-matter-only N-body prediction by DK15. For the polytropic exponent in the cluster core we refer the best-fit value to recent X-ray measurements of gas density profiles of clusters, and gas masses of groups and clusters as in Flender et al. (2017). We also adopt the value of parameters for feedback processes, stellar mass fraction as constrained in Flender et al. (2017) and nonthermal pressure as discussed in Nelson et al. (2014b). The parameters for gas clumping are derived from X-ray angular power spectrum measurements from ROSAT All-Sky Survey (Lau et al. in prep.). We fix a number of parameters in our parameter forecast. We fix x break = 0.195 because this value is set by the cooling radius where the gas cooling time becomes shorter than the Hubble time. Given the ICM profiles, there is little ambiguity in the location of the cool core radius. Furthermore, we fix the shape parameters of of the profiles of non-thermal pressure fraction and the gas clumping model: B nt , γ nt , α C , β C , and γ C , as they show little dependence in halo mass, redshift, and baryonic physics in cosmological hydrodynamical simulations.

Summary of model parameters
As a result, we vary 9 parameters for ICM physics in our fiducial analysis, shown in Table 1.

POWER SPECTRUM FORMALISM
In this section, we describe the formalism to compute the auto-and cross-power spectra between the tSZ compton parameter y, the X-ray surface brightness x, and lensing covergence κ, as shown in Eqs. (3), (4) and (6), respectively. These three fields are expected to be sensitive to the ICM physics and the dark matter potential well, allowing us to probe the connection between the ICM and dark matter content of galaxy clusters, and the growth of large-scale structure traced by galaxy clusters.
The cross-power spectrum between any two fields is given by: where · · · indicates the operation of ensemble average, δ (n) D (r) represents the Dirac delta function in n-dimensional space, A and B are the two-dimensional (i.e. projected) fields of interest.

Halo Model Approach
The auto and cross power spectra for any two fields C A B , under the flat-sky approximation, can be decomposed into two components within the halo-model framework (e.g. Cooray & Sheth 2002) where the first term on the r.h.s represents the two-point clustering in a single halo (i.e. the 1-halo term), and the second corresponds to the clustering term between a pair of halos (i.e. the 2-halo term where z min = 10 −3 , z max = 3, M min = 10 13 h −1 M and M max = 10 16 h −1 M , P L (k, z) is the linear matter power spectrum, dn/dM is the halo mass function, and b is the linear halo bias. We define the halo mass M by virial overdensity (Bryan & Norman 1998). We adopt the simulationcalibrated halo mass function presented in Tinker et al. (2008) and linear bias in Tinker et al. (2010). In Eqs. (30) and (31), A (z, M) and B (z, M) represent the Fourier transforms of profiles of fields A and B of a single halo with mass of M at redshift z, respectively. The fields A and B can be combinations of the Compton-y signal y, X-ray emissivity x, and the lensing convergence κ. In this paper, we consider the combinations: AB = x x, yy, x y, xκ, yκ. Note that we use a flat-sky approximation in Eqs. (30) and (31), The exact expression for the curved sky can be found in Appendix A of Hill & Pajer (2013). The Fourier transforms of the thermal electron pressure profile, y (M, z), the X-ray emissivity profile x (M, z), and the lensing convergence profile κ (M, z) of the halo with mass M and redshift z are expressed as where r 500 is cluster radius enclosing spherical region centered on the cluster center with average mass density 500 times critical density, P e,h is 3D electron pressure profile in single halo, E h is the X-ray emissivity profile defined in Eq. (27), ρ m is the mass density profile in a single halo and ρ m is the mean matter density in the universe, and q is the lensing kernel defined in Eq. (7). We define u = ar/r 500 and 500 = a χ/r 500 . The auto power spectrum of lensing convergence κ can be expressed, under the Limber approximation (Limber 1954), as where P m (k, z) is the non-linear matter power spectrum for wave number k at redshift z. We adopt the fitting formula of P m developed in Takahashi et al. (2012), which has been calibrated using a set of dark-matter-only N-body simulations for a variety of cosmological models. Note that this model does not capture baryonic effects on the cosmic matter density distribution (Jing et al. 2006;Rudd et al. 2008;Duffy et al. 2010;Zentner et al. 2013;Schneider & Teyssier 2015;Shirasaki et al. 2018;Chisari et al. 2018). Given the importance of this effect, we decide not to include C κκ in the our forecast of ICM physics and cosmology. We defer this to a future work where we account for baryonic effects on the matter power spectrum (e.g. Zentner et al. 2013;Eifler et al. 2015;Schneider & Teyssier 2015;Mohammed & Gnedin 2018;Chisari et al. 2018). Note that we still use this DM-only based C κκ for computing the covariance between the cross spectra between X-ray and lensing, C xκ and between tSZ and lensing, C yκ (see §5.1).

Effective halo redshift and mass probed by power spectra
To show the halo redshift and mass dependence in the auto and cross-power spectra, we compute the derivatives of the 1-halo term of the power spectrum C 1h with respect to redshift and mass. For any given fields AB, they are We skip the 2-halo term since its contribution is subdominant at angular scales > 100 for all the cases presented in this paper. Figure 1 shows the redshift-and mass-derivatives of tSZ and X-ray auto power spectra. In each panel, black, red and blue lines represent the derivative at = 100, 1000, and 10000, respectively. At larger (i.e. corresponding to the fluctuations at smaller angular separation), the power is dominated by lower mass halos at higher redshifts. Within the halo-model framework, the fluctuations in tSZ and X-ray maps for a given are determined by the halos whose angular size is of −1 . The typical angular size of a single halo is given by r vir (z, M)/ f K (z). Hence, at larger value, halos with smaller mass and at higher redshift (i.e. with smaller virial radii r vir and larger f K ) dominates the 1-halo power spectrum. Figure 2 summarizes the halo mass and redshift dependence in the cross power spectra. The general trend behaves similarly to the auto power spectra in Figure 1. However, the lensing cross-power spectra C xκ and C yκ are sensitive to a narrower range of halo redshifts since the lensing does not probe local halos with z ∼ 0 and halos with redshift close to the source z ∼ z source . For simplicity, we assume the source redshift z source = 1, with the effective redshift in lensing kernel (Eq. 7) ranging from 0.2 − 0.6. To probe the redshift evolution of ICM properties with lensing cross spectra, one need to adopt a "tomographic" approach, i.e., by using lensing sources with multiple redshifts.

STATISTICAL METHOD
In this section, we describe our methods for (1) estimating the statistical uncertainties in the model auto-and cross-power spectra presented in the §4 and (2) forecasting ICM physics and cosmological parameter constraints us-ing the auto-and cross-power spectra from upcoming multiwavelength cluster surveys.

Covariance modeling
We first estimate the covariance of a given power spectrum measurements. Assuming that power spectra follow Gaussian statistics, we can write the covariance matrix between two spectra C A B and C X Y as: where δ i j is the Kronecker delta, f sky is the fraction of sky used in the analysis, and we assume the power spectrum is measured with a finite bin width of ∆ . In Eq. (38), C A X = C A X + N A X represents the observed power spectrum including the observational noise, N A X . In general, the noise term is non-zero for the auto power spectrum. For instance, the auto power spectrum of tSZ has the noise term associated with the reconstruction method of tSZ map, while the cross power spectrum of tSZ and lensing is noiseless 7 . In the following, we estimate the noise terms in the auto power spectra of C yy , C xx and C κκ .

Noise power spectrum in tSZ
The noise in tSZ map arises from the reconstruction of a Compton-y map from observed temperature brightness maps at multiple frequencies in the microwave sky. Among several reconstruction methods in the literature (e.g. Bobin et al. 2008;Remazeilles et al. 2011Remazeilles et al. , 2013Hurier et al. 2013Hurier et al. , 2017Khatri 2015), we adopt the minimum-variance method as developed in Hill & Pajer (2013). This method finds the weight of multi-frequency maps in microwave so as to minimize the variance of the tSZ auto power spectrum.
In the framework of minimum-variance method, the noise power spectrum in tSZ can be expressed as the weighted sum of cross-correlation at different frequencies of foreground components: where w mv is the minimum variance weight and C i j is the cross-correlation at different frequencies. Note that the weight of w mv can be set from the observed C i j alone. The detail in the modeling of C i j is found in Appendix A.

X-ray point sources
In X-ray observations, there are various contributions to the cumulative X-ray surface brightness apart from the ICM. Among them, AGNs are believed to be a dominant component of the extragalactic X-ray emission (e.g.  The left panels shows the result for one-halo term of C y y , while the right is for one-halo term of C x x . In left and right panels, the upper panel represents the derivative of one-halo power spectra in redshift and the lower is the derivative in virial halo mass. In each panel, different colored lines show the difference in multipole of (black, red and blue correspond to = 100, 1000, and 10000, respectively).
Figure 2. Similar to Figure 1, but we show the differential contribution to cross power spectra with respect to redshift or halo mass. Left, medium, and right panels show the case of C x y , C x κ and C y κ , respectively. Note that we assume every source galaxy locates at z = 1 for lensing analyses in this figure.
X-ray flux count of normal galaxies is found to be close to one of the AGN in the faint end (Lehmer et al. 2012). Therefore, we consider two populations of X-ray point sources, AGNs and normal galaxies, as a noise in the cumulative Xray brightness of ICM (Eq. 4).
In this paper, we assume that all resolved point sources are masked. Assuming that unresolved X-ray point sources follow Poisson statistics, we can express the noise power spectrum of X-ray surface brightness N xx as where I X,α is the X-ray intensity from unresolved point sources α, n photon is the observed angular number density of photons in the energy range of E min to E max , and B X, is the Fourier transform of the point spread function in the X-ray survey (which we assume to be a Gaussian). In Appendix B, we summarize a phenomenological model of I X,α .

Shape noises in WL measurement
In modern galaxy imaging surveys, galaxy ellipticities are commonly used as estimators of the gravitational lensing field. Since the estimator is dominated by the intrinsic ellipticities of individual galaxies, the observed power spectrum of lensing convergence has the shot noise term, given by where we assume two source galaxy populations "a" and "b" with different redshift distributions, σ int,a is the root-mean-   square of the intrinsic ellipticity in the source catalogue of "a", and n gal,a represents the angular number density in the catalogue "a".

Survey specifications
We here define the specifications of future multi-wavelength cluster surveys in the microwave, X-ray and optical frequencies with a sky coverage of 20, 000 square degrees. In this work, as an illustration, we will focus on the three major upcoming cluster surveys: CMB-S4 in the microwave, eROSITA in X-ray, and LSST in the optical.

CMB Stage 4 (CMB-S4)
First, we consider the next-generation ground-based CMB experiment, referred to as CMB Stage 4 (CMB-S4; Abazajian et al. 2019). Following Madhavacheril et al. (2017), we adopt the experimental configuration consisting of a single large telescope with seven band-passes. Table 2 summarizes the different experimental configurations of CMB-S4 experiment setups considered in this paper. The baseline config-uration has a white noise level of s(ν) = 1.5 µK arcmin at the 150 GHz and 90 GHz channels. We examine different beam full-width half-maximum (FWHM) θ b in the 150 GHz channel from 1 arcmin to 3 arcmin, while we scale the beam FWHM in the other channels at a given ν as θ b ∝ 1/ν. For our baseline model, we use the configuration "S4:2arcmin" in Table 2.
Using the method in §5.1.1, we compute the tSZ noise power spectrum for a given experimental configuration. Figure 3 shows the noise tSZ power spectrum for the configurations described in Table 2. In the figure, the black line is the fiducial model of tSZ power spectrum by the halomodel approach (see Table 3 for the fiducial parameters of cosmology and ICM) and different colored lines represent the noise power spectrum for various experimental configurations. The red solid line is our baseline model of the noise power spectrum, while the red dashed and dotted lines are for the model with different knee . As shown in the figure, the beam size controls the amplitude of noise power spectrum at high multipoles at > ∼ 2000, while the angular parameter in atmospheric noise knee determines the amplitude of noise at < ∼ 1000. It is worth noting that the noise power spectrum dominates the observed power spectrum at < ∼ 1000 for all the cases considered in this paper. For simplicity, we assume perfect subtraction of the noise power spectrum from the observed one in the fisher analysis, whereas we discuss the possible biases in parameter estimation due to the imperfect source substraction in §6.4.1.

eROSITA
For X-ray, we consider the upcoming eROSITA X-ray allsky survey (Merloni et al. 2012). We adopt a beam size of 30 arcsec FWHM, an energy range of 0.5 − 2 keV in the observer's frame, and an effective area of 1365 cm 2 . The total observation time is assumed to be four years. The flux limit in the energy range of 0.5 − 2 keV is set to be F lim = 4.4 × 10 −14 erg cm −2 s −1 , and the mean X-ray intensity from unresolved point sources is estimated to be 14.7 cm −2 s −1 str −1 . We assume all the point sources with the flux greater then F lim are removed.

LSST
For lensing, we consider a future galaxy imaging survey performed with the Large Synoptic Survey Telescope (LSST; LSST Science Collaboration et al. 2009). We assumed the source redshift distribution to follow where α = 1.27, β = 1.02, z 0 = 0.5. The normalization is set so that the total number density of the source galaxies is 26 arcmin −2 , which represents the largest number density of sources achievable by ground-based imaging surveys. For the shape noise, we adopt the RMS of source ellipticity to be 0.367, the typical value found in the ground-based imaging surveys (e.g., see Mandelbaum et al. 2018).
To probe the redshift evolution of the ICM and the growth of clusters, we divide the source galaxies into N tomographic subsamples, where we take into account the un- certainties in the photometric redshift z ph with the following probability distribution for a given true redshift z: We set the scatter of z ph to be σ z = 0.05(1 + z) −1 , a typical estimate in the ongoing imaging survey (Tanaka et al. 2018). The bias parameter z b is assumed to be zero. We discuss possible biases in parameter estimation due to the imperfect knowledge of z b in §6.4.2. Given Eqs. (42) and (43), we define the redshift distribution of i-th source sample as and we set the bin width of z ph so that each subsample of source galaxies has the same number density. Figure 4 shows the source redshift distributions when we adopt 10 tomographic bins (N = 10) and its corresponding kernel for lensing convergence (Eq. 7).
In addition, we also consider the multiplicative uncertainties in the lensing convergence of i-th tomographic bin, κ obs → (1 + m i ) κ true , where κ true is given by Eq. (6) and m i is the so-called multiplicative bias. As in our fiducial setup, we set m i = 0 for each tomographic bin. We then examine the impact of non-zero multiplicative biases in parameter estimation in §6.4.2. Table 3. The list of parameters in our fisher analysis. Assuming N tomographic bins in lensing survey, we have 15+2N parameters in total, consisting of 6 cosmological parameters, 9 ICM parameters (see Table 1  Lensing systematics m i (i = 1, · · · , N ) 0.00 0.10 z bi (i = 1, · · · , N ) 0.00 0.05

Parameter forecast: Fisher analysis framework
We here summarize a Fisher analysis and the fiducial parameters for computing power spectra. For a multivariate Gaussian likelihood, the Fisher matrix F i j can be written as where Cov is the data covariance matrix, µ represents the assumed model, and p i describes parameters of interest. The Fisher matrix provides an estimate of the error covariance for two parameters as where ∆p α represents the statistical uncertainty of parameter p α . Since the data covariance is expected to scale inversely with the survey area, the second term in Eq. (45) becomes dominant for a large area survey. We therefore consider only the second term. We vary 6 cosmological and 9 ICM parameters, and fix the rest as shown in Table 3. Note that the cosmological parameters consist of physical dark matter density Ω cdm h 2 , physical baryon density Ω b0 h 2 , the amplitude of initial curvature fluctuations A s (at the pivot scale k = 0.05 Mpc −1 ), the spectrum index parameter of linear power spectrum n s , the parameter of equation of state of dark energy w 0 , and the dimensionless Hubble parameter h. In addition, we introduce 2N parameters for lensing systematics when using N tomographic bins in lensing field κ (see §5.2.3).
We then construct the data vector D from a set of binned spectra C yy , C xx , C xy , C xκ,i and C yκ,i , where the index i is for the tomographic redshift bins of the lensing source galaxies. To compute the power spectra, we perform the logarithmic binning in the multipole range from min = 30 to max = 3000. The number of bins in each power spectrum is set to be 10. Hence, we have 10+10+10+10N +10N = 20N +30 elements in the data vector for N tomographic bins in total.
Owing to the strong degeneracy among cosmological and ICM parameters in the power spectra (especially for C yy ), we need to introduce loose priors for realistic forecast of parameter constraints. In this paper, we compute the Fisher matrix as follows: where F model is computed by Eq. (45) with the power spectra, and F prior represents the prior information. We assume that the loose prior on the i-th parameter is given by Gaussian with the variance, σ prior,i . In this case, the prior term is expressed as F prior,i j = σ −2 prior,i δ i j . The fiducial value of σ prior,i is found in Table 3. We set the value of σ prior to be 30 − 450% of fiducial value of parameters.

Parameter dependence
We begin by studying the dependence of parameters on various power spectra. Figure 5 shows some representative examples of predicted power spectra for a range of cosmological models and ICM parameters. In this figure, we consider six cases as follows: (i) The best-fit cosmological model by Planck Collaboration et al. (2016a) (denoted as Planck15) and fiducial ICM parameters; (ii) The best-fit cosmological model by Hinshaw et al. (2013) (denoted WMAP9) and fiducial ICM parameters; (iii) The Planck15 cosmology and the ICM model with a parameter of non-thermal pressure profile being A nt = 0.60 (A nt = 0.451 for fiducial case); (iv) The Planck15 cosmology and the ICM model with a parameter in stellar-to-halo mass relation being S * = 0.62 (S * = 0.12 for fiducial case); (v) The Planck15 cosmology and the ICM model with a parameter of stellar feedbacks being f = 10 −5 ( f = 4 × 10 −6 for fiducial case); (vi) The Planck15 cosmology and the ICM model with a parameter of gas clumping being C 0 = 2.0 (C 0 = 0.9 for fiducial case).
Here, we set the weak lensing source redshift to be 1 and source number density to be 10 arcmin −2 for simplicity. The black points in the figure show the model (i), while the blue, green, red, cyan and purple lines indicate the model (ii), (iii), (iv), (v), and (vi), respectively. The gray error bars represent the statistical uncertainty evaluated in §5.1 for a survey specification in §5.2.
The difference in the black points and blue lines highlights the dependence of cosmological parameters on the power spectra. Note that the main difference between the Planck15 and WMAP9 cosmological models lies in the cosmic matter density Ω m0 and the present-day normalization of linear matter power spectrum σ 8 (both being lower in the WMAP9). Since our power spectra will be mainly determined by the 1-halo term, the dependence on cosmological parameters stems primarily from that of the halo mass function. In particular, the tSZ power spectrum has a significant dependence on cosmology compared to other spectra, because it contains the information of halo mass function at massive cluster scales (see Figure 1). Furthermore, the cosmological dependence among power spectra is different, which allows the breaking of the degeneracy among the parameters, thus in turn improving the cosmological constraints. Interestingly, the X-ray auto power spectrum depends on cosmology in an opposite way as the tSZ power spectrum does. We expect that the X-ray auto power spectrum to be sensitive to Ω m0 , since it depends on the gas density squared in a single halo. Given that the cosmic baryon density Ω b0 is relatively unchanged, as one reduces Ω m0 , the gas mass in a single halo increases, thus increasing the overall amplitude of the X-ray auto power spectrum. This explains why the X-ray auto power spectrum amplitude is higher in the WMAP9 cosmology (with lower Ω m0 = 0.27) than Planck15 (with Ω m0 = 0.30). While the tSZ power spectrum is more sensitive to σ 8 (comparing the black points and red lines in Figure 5), the X-ray power spectrum is a more efficient probe of Ω m0 . Uncertainties in ICM model will degrade the cosmological constraints from C yy and C xx . The comparison among the black points, green, red, and cyan lines in Figure 5 demonstrates how ICM physics can hide the cosmological dependence on the power spectra. Reducing the non-thermal pressure fraction (green lines) can increase the tSZ power spectrum, while the X-ray auto power spectrum is less sensitive to the change in A nt . On the other hand, increase in the slope in the stellar-to-halo mass relations S * (red lines) can affect the X-ray observables C xx and C xκ . The parameter degeneracy in ICM can be seen in the comparison between red and cyan lines in Figure 5. Increasing the amplitude of the stellar+AGN feedback can have a similar effect as increasing the non-thermal pressure in C yy and C yκ , as both reduce the thermal pressure in less massive group-sized halos. Importantly, this degeneracy can be broken if one includes X-ray information, since the X-ray observable is less sensitive to the non-thermal pressure, but still dependent on stellar and AGN feedback. However, the X-ray observable is dependent on gas clumping, while the tSZ is insensitive to this effect (purple dash-dotted lines). The cross correlation with lensing will provide additional constraints on cosmology and ICM physics, but the amplitude in lensing observables is subject to the imperfect shape measurement and photo-z accuracy.

Fisher forecast for future wide-area surveys
We present the results of Fisher forecast in our fiducial setup as in §5.3. Figure 6 summarizes the forecast for a survey with a sky coverage of 20,000 square degrees when 10 tomographic bins in lening observables are available. In this  Table 1 for the ICM parameters), while the gray squared box represents the expected statistical uncertainty in hypothetical surveys with a sky coverage of 20,000 square degrees (see also §5.2). The black dotted line in the bottom panels shows the 2-halo term which is always sub-dominant. Different colored lines in each panel show the parameter dependence on the model power spectra. Blue lines correspond to the model with the best-fit cosmology to the WMAP CMB measurement (Hinshaw et al. 2013) and fiducial ICM, while green, red, cyan and purple are for the model with the Planck cosmology and higher A nt , higher S * , higher f , and higher C 0 , respectively. Note that A nt represents the parameter about overall amplitude in non-thermal pressure, S * is the power-law index in the stellar-to-halo mass relations in cluster-sized halos, and f controls the strength of stellar+AGN feedback inside galaxy clusters. The parameter of C 0 determines amplitude of the the clumping factor in gas density profile in the cluster outskirts. The y-values in the upper and lower panels are multiplied by 1.5 and 1.3 respectively to better show the dependence on the parameters. figure, each panel shows a set of relevant physical effects in the joint analysis of tSZ, X-ray and lensing observables. In each panel, the orange shaded region is derived by the Fisher matrix analysis when we limit the information of C yy and C yκ as examined in Battaglia et al. (2015b) and Osato et al. (2018). The purple region represents the constraints if we can combine the information from C yy , C yκ , C xx , C xy and C xκ .

Cluster Physics Constraints
In addition to cosmolology, the joint analysis will also provide an important insight to the ICM physics inside virial radius and in the cluster outskirts. In our ICM model, we vary only the normalization parameter of the non-thermal pressure, A nt , while keeping the shape of its profile fixed. We define the 68% confidence level of the fraction of non-thermal pressure by using the marginalized fisher matrix of A nt . The top left panel of Figure 6 shows the non-thermal fraction of ICM in cluster-sized halos as a function of the normalized cluster-centric radius, r/r 200m . The joint analysis allows us to constrain the fraction of non-thermal pressure and gas clumping at the level of 22/50% (1σ error), respectively at r = r 200m . This represents a factor of two improvements compared to previous studies (Battaglia et al. 2015b;Osato et al. 2018). Note that there is currently no other method that can observationally constrain these ICM physics for a statistical sample of galaxy clusters. The bottom right panel shows the constraints on the astrophysical feedback parameters in the ICM model. The ICM model without feedback processes can be ruled out with a ∼ 1σ level. We also find the stellar mass fraction in massive galaxy clusters can be constrained to within ∼ 5%. Feedback processes Figure 6. The summary of parameter forecasts with a joint analysis of tSZ, X-ray and lensing observables. The orange error contours assume the analysis with the power spectra related with tSZ and lensing (C y y and C y κ ), while the purple one shows the results by the joint analysis of tSZ, X-ray and lensing (C y y , C x x , C x y , C x κ and C y κ ). In this figure, we assume our fiducial set up of multi-wavelength surveys covering 20,000 square degrees. All the error contours are for 68% confidence level. Note that the gas clumpiness can be constrained only if the X-ray information would be available.

Cosmological Constraints
The top right panel of Figure 6 shows the constraint on the cosmic expansion model summarized by two parameters of w 0 and H 0 . We find that the dark energy equation-of-state parameter can be constrained at the level of 1σ = 0.080 from the joint analysis of tSZ, X-ray and lensing observables. That is comparable to the current constraints (e.g. Huterer & Shafer 2018, for review), but we do not use the tightest prior information of cosmology from CMB and SNe measurements. Therefore, our joint analysis approach is best to be applied at lower redshifts z < ∼ 1, where the dark energy is expected to dominate the cosmic energy density.

Breaking Parameter Degeneracy with X-ray
As shown in Figure 6, the X-ray information can improve the constraints of ICM model, and help break some of the degeneracy among ICM model and cosmological parameters. Note that we include the uncertainty of clumping effects in X-ray modeling and marginalize it in the Fisher forecast. Although some nuisance parameters can increase the uncertainty in the amplitude of power spectra, the joint analysis is still able to improve the parameter constraints in ICM and cosmology. Figures 7 and 8 show the marginalized error circles in different two-parameter planes as we combine various information of the power spectra, highlighting that the Xray information can play a crucial role in breaking parameter degeneracies. The marginalized and un-marginalized errors of each parameter in our fiducial Fisher analysis are also provided in Table 4.

Dependence of analysis and survey setups
We then discuss how we can tighten the constraints furthermore by changing survey configurations and analysis methods. In the following, we examine to what extent the constraint on non-thermal presssure A nt will be affected by various choices in prior information, number of lensing tomographic bins, smallest scales in the analysis, and the detailed configuration in CMB-S4. The results are summarized in Figure 9. Similar tests on constraint for the equation-ofstate parameter of dark energy are given in Appendix C.

Impact due to priors
Since our power spectra have a complex parameter dependence and most of them degenerates each other, it would be important to set appropriate prior information for improving parameter constraints.
We first examine the prior information from the CMB C yy C yy + C yκ C yy + C xx + C xy C yy + C xx + C xy + C xκ + C yκ   Figure 7. The corner plot of 6 cosmological parameters by different combinations of power spectra. In this figure, different colored lines show the parameter constraints with a 68% confidence level by various combination of tSZ, X-ray and lensing observables (as shown in the legend). Note that the pink circles are based on the combined analysis with tSZ and lenisng as proposed in the literature (Battaglia et al. 2015b;Osato et al. 2018), while the blue one shows the significant improvement in our joint-analysis method.
measurements by Planck satellite (Planck Collaboration et al. 2016a). To do so, we use the approximate covariance of parameter estimates as shown in Makiya et al. (2018), denoted as C P15 . After marginalizing the uncertainty of optical depth, we set the prior term in Fisher matrix (Eq. 47) to be F prior =C −1 P15 , whereC P15 represents the prior infor-mation of covariance over five cosmological parameters of Ω cdm h 2 , Ω b0 h 2 , A s , n s and h. We find that using the prior information from the Planck CMB measurments can tighten the constraint of A nt by a factor of 3.8.
Another possible informative prior information will be included in the systematic effects of lensing measurements.  Figure 8. Similar to Figure 7, but focusing on 9 parameters of cluster astrophysics. The same legend is adopted as in Figure 7. It would be worth noting that the constraints without X-ray information are mostly determined by the priors.
We change the prior Gaussian uncertainty of multiplicative bias m and photometric redshift bias z b in each tomographic bin to be 0.01, but fixing other prior scatters same as in Table 3. When applying a 1% Gaussian scatter in m and z b , we find that the expected constraint on A nt is almost unaffected (it differs from the fiducial value by < ∼ 0.5% level). This indicates that the joint analysis of tSZ, X-ray and lenisng power spectra will be insensitive to the choice of prior information of lensing systematics.
Detailed observations of individual galaxy clusters can further improve our understanding of ICM physics in the cluster cores and outskirts. To see the impact of prior information on the ICM parameters, we simply examine the case with the perfect calibration of ICM parameters for the cooling in core regime (Γ and γ) and gas clumping effects Table 4. The 1σ uncertainty evaluated by fisher analysis with a joint analysis of tSZ, X-ray, and lensing observables assuming a sky coverage of 20,000 square degrees. Note that we assume 10 tomographic bins are available in lensing information. See Table 1 for details of the ICM model.
(C 0 ). Even if we can determine those parameters with an unprecedented accuracy, the joint analysis of tSZ, X-ray and lensing power spectra can constrain A nt at a level of 0.078 (i.e. ∼ 10% improvement). In addition, varying the shape of universal profile of non-thermal pressure and gas clumping (i.e. B nt , γ nt , α C , β C and γ C ), the expected 1σ error on A nt will be degraded by 58%. This highlights the importance of external priors on the ICM model, especially in the X-ray band.

Efficiency of lensing tomography
Tomographic lensing analysis is useful for probing redshift evolution in ICM physics. Since there is no standard procedure to set the number of tomographic bins, we first examine the different number of tomographic bins from 2 to 12 with a fixed total source number density and check how our fiducial analysis is affected by the number of tomographic bins. We find that the analysis with smaller tomographic bins slightly degrades the constraint on A nt (by 6% and 3% from the fiducial value when adopting the 2-bin and 4-bin tomography, respectively), but the effect is generally small. Increasing the number of tomographic bins does not always provide better constraints on the redshift evolution of the parameters, since there remains finite cross covariances between different tomographic bins (i.e., the correlation statistics between two different tomographic bins are not independent). It is difficult to collect independent information in redshift direction from lening observable alone. As the lensing effect is caused by the intervening gravitational potential between source and observer, the lensing kernels share the redshift information among different tomographic bins. From the practical point of view, it is better to set the elements of data vector as small as possible, because the inversion of covariance matrix will be inaccurate for large dimensions.

Importance of small-scale information
An aggressive analysis of the power spectra with smallerscale information has the potential to improve the constraints on cluster astrophysics. It is, therefore, interesting to examine how much additional information can be gained by increasing the maximum multipole max in the Fisher analysis. We examine four different values of max = 1000, 2000, 4000 and 5000 in order to assess the effect of max on our fiducial analysis. We find the constraint of A nt is degraded by 104% for max = 1000, whereas it is improved by 12% for max = 5000.

Dependence of noise power spectrum in tSZ measurements
Finally, we examine the impact of the noise property in tSZ measurements on the Fisher analysis. We consider four different noise power spectra in the tSZ auto power spectrum N yy , indicated using blue solid, yellow solid, red dashed and red dotted lines in Figure 3. We also vary the covariamce matrix of the power spectra with different N yy , while keeping other survey parameters fixed. We find that the beam size in CMB-S4 is less important for the cosmological constraints, whereas the amplitude of atmospheric noise can change the constraint of A nt at a level of ∼ 1%. This indicates that precise measurements of the tSZ auto power spectrum at large scales (smaller ) will be critical for constraining the nonthermal pressure in cluster outskirts.

Possible biases in parameter estimation
Next, we investigate possible sources of systematic effects in parameter estimation. In particular, we would like to understand how large parameter biases will be induced within our current modeling framework. The current framework assumes the perfect subtraction of noise in tSZ and X-ray auto power spectra and the systematic effects in lensing measurements are under control. In practice, astrophysical and cosmological constraints are affected by observational biases. Here, we consider two main sources of observational biases originating from (i) imperfect subtraction of noise power spectrum in tSZ auto power spectrum 8 , and (ii) misestimation of multiplicative and photometric redshift biases in lensing measurements. To do so, we adopt the method to estimate the bias in parameter estimation (Huterer et al. 2006): where p is the parameter of interest, D fid is the data vector for fiducial model, D obs is the data vector affected by the statistical fluke, and Cov represents the data covariance.

Impact of astrophysical sources
We first study the impact of imperfect subtraction of noise power spectrum in tSZ measurements. The expected signalto-noise ratio of tSZ noise power spectrum is given by where Var[N yy ( i )] represents the variance of N yy at i-th multipole i and we assume the Gaussian variance in the second line in Eq. (49). In our fiducial setup, we find (S/N) = , where Ω survey is the survey coverage. Therefore, the amplitude of N yy is estimated to be ∼ 1/1470 = 0.068% in the survey covering 20,000 square degrees. Note that N yy can be estimated with the observational information alone (see Eq. 39).
In other words, the imperfect subtraction of N yy will lead to an uncertainty of 0.068% in the tSZ measurements. 8 The imperfect subtraction of X-ray shot noise is found to have a negligible impact of our Fisher forecast, because the shot noise term is sufficiently small compared to the expected C x x in the range of < ∼ 3000.  This small residual may induce the biased parameter estimations. We now model the residual of N yy by a simple form of AN yy , where A is the constant value in multipole and should be an order of 0.068%. Assuming there remains a residual, we can express the difference of D obs − D fid in Eq. (48) as AN yy for the term of C yy and zero otherwise. Figure 10 shows the possible parameter bias as a function of A. The left panel represents the bias of cosmological parameters with respect to the 1σ marginalized error, while the right is for some ICM parameters. We find the residual of N yy with a level of 0.1% can induce < ∼ 0.1σ-level bias in the estimated parameter.

Photometric redshift estimates and Imperfect shape measurements
We next study the impact of systematic effects in lensing observables. We assume the systematic parameter in lensing (m and z b ) can be biased in actual observations. Specifically, we consider the bias in m or z b in each tomographic bin with a level of 0.01 and investigate how large bias δ p can be induced by this small bias in m or z b . Figure 11 shows the bias of parameter estimation when the photo-z estimate in a single tomographic bin is biased by 0.01. According to this figure, the 1% error in photo-z bias z b can induce the bias in parameter estimation by ∼ 0.005σ. We also confirm that the 1% mis-estimation of multiplicative bias m causes ∼ 0.01σ-level parameter biases. Hence, we expect the systematic error due to multiplicative and photo-z biases will be subdominant in the joint analysis of tSZ, X-ray and lensing observables, as long as we calibrate the systematic bias in lensing at 1% level.

LIMITATIONS
We summarize the major limitations in our correlation statistics for multi-wavelength cluster surveys. All of the following issues will be addressed in forthcoming studies.

Uncertainties in the ICM model
Although our semi-analytic ICM model is an important step toward an effective utilization of multi-wavelength observables for precision cosmology and understanding ICM physics, it is still insufficient to apply to actual observation data sets. For instance, our model assumes that there are no contributions to the observables from diffuse gas in large-scale filaments. We also ignore the dependence of mass assembly history on gas profiles (e.g. Nelson et al. 2014b;Shi et al. 2015;Lau et al. 2015) and dark matter density profile (e.g. Wechsler et al. 2002;Diemer & Kravtsov 2014) throughout this paper, which might induce the bias in estimation of the parameter of interest in real data analysis. To study this issue, we require a synthetic observational data generated from cosmological hydrodynamic simulations to fully take into account relevant ICM physics and non-linear gravitational growth.

Additional correlation from astrophysical sources
Our model assumes that the observed power spectra are dominated by the ICM. In actual observations, correlations can be induced by other faint astrophysical sources, including the extragalctic radio source and CIB in tSZ, and faint unresolved AGNs and galaxies in X-rays. Since these sources also correlate with each other and with the lensing observables, it may induce non-negligible levels of parameter biases. The modeling of astrophysical sources in correlation statistics is still developing (e.g. Cappelluti et al. 2013;Zentner et al. 2013;Hurier 2015;Planck Collaboration et al. 2016b;Li et al. 2018;Shirasaki 2019) and require appropriate parametrizations to be taken into account in future studies.

Intrinsic alignment in cluster-sized halos
The intrinsic alignment (IA) effect of observed galaxies in galaxy imaging surveys is one of the major contaminants in modern lensing analyses (e.g. Troxel & Ishak 2015, for review). This effect is likely more prominent in redder galaxies, indicating that member galaxies in cluster regimes will be more affected by the IA (e.g. see Singh et al. 2015;Samuroff et al. 2018;Hikage et al. 2019;Johnston et al. 2019, for recent observational studies). Unfortunately, it is still difficult to include the IA effect in our halo-model approach in a selfconsistent way. Although there exists a simple halo model of the IA effect (e.g. Schneider & Bridle 2010), the halo-model parameters have been poorly constrained so far. Further advances in theoretical modeling and observational studies are needed to implement IA effects in the halo model.

Non-Gaussian covariance
The statistical error in correlation statistics is assumed to be Gaussian throughout this paper, whereas the non-Gaussian covariance may play an important role in actual observations. The non-Gaussian covariance should arise from the non-linear gravitational growth of structures and it is computed as the four-point correlation function among observables of interest. For instance, the non-Gaussian covariance in tSZ auto power spectrum C yy dominates the Gaussian error at < ∼ 100 in the recent measurement (e.g. Horowitz & Seljak 2017;Bolliet et al. 2018;Osato et al. 2018). Nevertheless, the impact of non-Gaussian covariance on C yy is less critical for our cases, because we work with groundbased experiments and the noise at small is dominated by the atomospheric noises (see Figure 3). On the other hand, the noise covariance in X-ray auto power spectrum C xx in eROSITA-like survey will be dominated by sample variance, indicating that non-Gaussian error may become important.
Analogously to the C yy case (e.g. Komatsu & Seljak 2002;Hill & Pajer 2013), we expect that our Gaussian covariance in C xx may be underestimated by at most a factor of ∼ 10. It is still uncertain if the non-Gaussian covariance in C xx can induce a large degradation of parameter constraints, since we have 30 − 40 model parameters in total. We will address this issue with a set of synthetic observations to evaluate precise covariance and perform a realistic likelihood analysis.

CONCLUSIONS
Upcoming multi-wavelength galaxy cluster surveys are poised to significantly improve our understanding of astrophysics and cosmology of the large-scale structures. In this paper, we present a theoretical framework to predict the correlation statistics among the different galaxy cluster observables: thermal Sunyaev-Zeldovich (tSZ) effect in cosmic microwave background (CMB), X-ray emission from ionized intracluster gas, and weak gravitational lensing effect of background galaxies due to the cluster's gravitational field.
Our model framework is based on a physicallymotivated semi-analytic model of ICM as presented in §3.
Using the model, we compute the auto and cross power spectra of the observables and perform a Fisher analysis to forecast the constraints on cosmology and ICM physics in future CMB ground-based experiment, X-ray satellite mission, and galaxy imaging survey. Our findings are summarized as follows: • With joint tSZ, X-ray, and lensing auto and cross angular power spectra from multiwavelength surveys with 20,000 square degrees sky area coverage, we will be able to constrain the poorly understood cluster astrophysics: non-thermal pressure, gas clumping, and feedback from supernovae (SNe) and active galactic nuclei (AGN), at a level of 4.4σ, 2.0σ, and 1σ respectively. Although the improvements on individual cluster astrophysics parameters are modest due to the degeneracy among them, our gas model is sufficiently flexible in marginalising over the poorly understood cluster astrophysics parameters in a physically motivated way.
• Marginalising over the aforementioned cluster astrophysics, the joint tSZ, X-ray, and lensing auto-and crossangular power spectra will allow us to constrain the equation-of-state parameter of Dark Energy, w 0 , at a precision of 8%, which is comparable to the precision achieved by cluster abundance measurements. Further improvements of the ICM model are required for its applications to real data sets ( §7).
• Our analyses indicate that tSZ noises will become the major systematic uncertainties, compared to lensing systematics, in affecting the level of bias in the ICM physics con-straints. In particular, imperfect subtraction of tSZ noises at a 1σ level of the signal can induce ∼ 10% bias in the ICM physics parameters with respect to the marginalized uncertainties. Whereas 1%-level biases in lensing systematics leads only to 0.5 − 1% bias.
With data from a plethora of multi-wavelength largescale astronomical surveys coming online in the next decade, our framework will be useful for extracting valuable information about cosmology and ICM physics (such as feedback from SNe and AGN, non-thermal pressure, gas clumpng in cluster outskirts) simultaneously. Since these ICM physics are the main source of systematic uncertainties for cluster cosmology, their constraints from the joint analysis can provide meaningful priors for constraining cosmology with measurements of cluster abundance and spatial clustering.
Aside from ICM physics, our physics-based framework can also be extended to smaller-mass haloes, such as groups and galaxies, as well as different cosmological models. In the future works, we plan to extend our model in order to interpret the cross-correlation measurements involving galaxy properties (e.g. Planck Collaboration et al. 2013;Greco et al. 2015;Vikram et al. 2017;Hill et al. 2018;Makiya et al. 2018) as well as test non-standard cosmologies involving modified gravity, massive neutrinos, and variants of non-cold dark matter models. and their frequency scaling, (iii) the tSZ-CIB correlation, (iv) the extragalactic radio source, (v) and thermal dust emission from Galactic cirrus. We use the best-fit parameters of the model fitted to the data from Atacama Cosmology Telescope (ACT, see Dunkley et al. 2013, for details). For the C N term, we adopt the form presented in Madhavacheril et al. (2017) and given by where s(ν) is the instrumental white noise at frequency ν and the remaining term in the right hand side expresses the effect of atmospheric noise for a ground-based experiment. Throughout this paper, we fix α = −4.5, but examine three different values of knee = 2500, 3500 and 4500. We adopt knee = 3500 unless otherwise stated.

APPENDIX B: A MODEL OF X-RAY POINT SOURCES
In this appendix, we summarize a phenomenological model of X-ray point sources. In the following, we consider two major populations at X-ray bands, AGN and normal galaxies. Suppose that all resolved point sources can be masked, the remaining X-ray intensity from an unresolved point source α is given by (also see Zandanel et al. 2015) W X,α (E, z) = 1 4π ∫ L X,max (z) L X,min dL X Φ X,α (L X , z) L X,α (E, z), (B2) where α = {AGN, Galaxy}, Φ X,α is the X-ray luminosity function of the source α, and L X,α represents the differential luminosity at an energy E and redshift z. Note that L X,α is defined as a number of photon emitted per unit time and per unit energy range. Hence, it depends on X-ray energy spectrum of the source α and the definition of L X . Also, the upper bound of luminosity L X,max (z) is set so as to compute all contributions with X-ray flux less than the flux limit of the survey. To compute Eq. (B1), we adopt the model of Φ X,AGN developed by Aird et al. (2010) in the rest-frame energy range of 2 − 10 keV, while we use the model in Ptak et al. (2007) for Φ X,Galaxy in the rest-frame energy range of 0.5 − 2 keV. The differential luminosity is then computed as where E min,ref and E max,ref are the minimum and maximum rest-frame energy in the definition of L X , respectively. We assume the AGN spectrum index to be Γ α = 1.7 (e.g. Tozzi et al. 2006) and Γ α = 2 for normal galaxies (e.g. Young et al. 2012) in this paper. For the lower and upper limit of luminosity integration in Eq. (B2), we set L X,min = 10 41 erg s −1 and L X,max (z) = 4πd 2 L (z) F lim K(z, E min , E max ), where d L (z) is the luminosity function to redshift z, F lim is the X-ray flux limit in the survey, and K is the K-correction given by −1.0 −0.8 −0.6 Fisher forecast on w 0 Fiducial set up Planck Prior 1% calibration in photo-z bias 1% calibration in multiplicative bias Varied the shape of universal profile w/o lensing systematics Fixed cooling core and clumping physics 2-bin tomography 4-bin tomography 6-bin tomography 8-bin tomography 12-bin tomography max = 1000 max = 2000 max = 4000 max = 5000 1 arcmin beam in CMBS4 3 arcmin beam in CMBS4 knee = 2500 in CMBS4 atmospheric noise knee = 4500 in CMBS4 atmospheric noise LSST × CMB − S4 × eROSITA (20, 000 deg 2 ) Figure B1. Similar to Figure 9, but we consider a representative example of cosmological parameters. The panel shows the dependence of analysis and survey setups on the parameter of dark energy w 0 .
for AGN with Γ = 1.7 and K(z, E min , E max ) = ln(E max /E min ) ln(E max,ref /E min,ref ) .
for normal galaxies.

APPENDIX C: DEPENDENCE OF ANALYSIS AND SURVEY SETUPS ON COSMOLOGICAL CONSTRAINTS
In this appendix, we provide the results of the expected constraints of the parameter of dark energy when varying the analysis and survey setups as examined in §6.3. Figure B1 summarizes the results. To improve the constraint of w 0 , we need a prior cosmological information from primary CMB fluctuations to reduce the statistical uncertainty in cosmic baryon and mass density, and the spectral index of primordial curvature fluctuations. Also, a small-scale information would be helpful to break the parameter degeneracy among ICM physics.
This paper has been typeset from a T E X/L A T E X file prepared by the author.