Sunyaev--Zel'dovich profile fitting with joint AMI-Planck analysis

TL;DR

This paper introduces a Bayesian joint analysis method for Sunyaev-Zel'dovich data from AMI and Planck, improving cluster parameter constraints and resolving discrepancies in measurements.

Contribution

A novel Bayesian approach for joint SZ data analysis from interferometers and satellites, enabling unbiased cluster parameter estimation and profile shape constraints.

Findings

Joint dataset yields precise cluster parameters when pressure profile is known.

Allowing profile shape variation removes bias and constrains shape parameters.

Method resolves discrepancies in Planck and AMI measurements for a real cluster.

Abstract

We develop a Bayesian method of analysing Sunyaev-Zel'dovich measurements of galaxy clusters obtained from the Arcminute Microkelvin Imager (AMI) radio interferometer system and from the Planck satellite, using a joint likelihood function for the data from both instruments. Our method is applicable to any combination of Planck data with interferometric data from one or more arrays. We apply the analysis to simulated clusters and find that when the cluster pressure profile is known a-priori, the joint dataset provides precise and accurate constraints on the cluster parameters, removing the need for external information to reduce the parameter degeneracy. When the pressure profile deviates from that assumed for the fit, the constraints become biased. Allowing the pressure profile shape parameters to vary in the analysis allows an unbiased recovery of the integrated cluster signal and produces constraints on some shape parameters, depending on the angular size of the cluster. When applied to real data from Planck-detected cluster PSZ2 G063.80+11.42, our method resolves the discrepancy between the AMI and Planck $Y$-estimates and usefully constrains the gas pressure profile shape parameters at intermediate and large radii.

Figures

Tables

Table 1. Table 1: Cluster parameter prior distributions. 𝒩 𝒩 \mathcal{N} denotes a Gaussian distribution parameterised by its mean and standard deviation, 𝒰 𝒰 \mathcal{U} denotes a uniform distribution, and δ 𝛿 \delta is a Dirac delta function. In the cases where the latter is used, the values used for the function’s argument will be stated when the analyses are carried out.

Parameter	Prior distribution(s)
$x_0$	$𝒩(0^{\prime \prime },{300}^{\prime \prime })$
$y_0$	$𝒩(0^{\prime \prime },{300}^{\prime \prime })$
$Y_{\mathrm{tot}}$	$𝒰[0.00{\mathrm{arcmin}}^2,0.02{\mathrm{arcmin}}^2]$
${\theta }_s$	$𝒰[{1.3}^{\prime },{15}^{\prime }]$
$\alpha$	$\delta \left({\alpha }_{\mathrm{model}}\right)$ or $𝒰[0.1,3.5]$
$\beta$	$\delta \left({\beta }_{\mathrm{model}}\right)$ or $𝒰[3.5,7.5]$
$\gamma$	$\delta \left({\gamma }_{\mathrm{model}}\right)$

Table 2. Table 2: Cluster simulation inputs using an observational model and three different pressure profile shapes. x 0 subscript 𝑥 0 x_{0}  and y 0 subscript 𝑦 0 y_{0}  are always 0, i.e. the cluster is at the simulated map centre.

		$Y_{\mathrm{tot}}$	${\theta }_s$	$\gamma$	$\alpha$	$\beta$
		/ arcmin2	/ arcmin
UPP	low SNR	0.001	2	0.3081	1.0510	5.4905
	high SNR	0.007	8
PIP	low SNR	0.001	2	0.31	1.33	4.13
	high SNR	0.007	8
REP	low SNR	0.001	2	0.065	0.33	5.49
	high SNR	0.007	8

Table 3. Table 3: Summary of parameter constraints for PSZ2 G063.80+11.42, varying α 𝛼 \alpha  and β 𝛽 \beta . The errors given are the 68% limits from the one-dimensional marginalised parameter constraints.

	AMI-only	Planck-only	AMI-Planck
$x_0$/ arcsec	$-{59}_{-18}^{+24}$	$-{38}_{-26}^{+31}$	$-55\pm 14$
$y_0$/ arcsec	$-{70}_{-14}^{+26}$	$-45\pm 31$	$-65\pm 15$
$Y_{\mathrm{tot}}$/ (${10}^3$ arcmin2)	${6.9}_{-5.7}^{+2.3}$	${1.32}_{-0.38}^{+0.17}$	${1.40}_{-0.34}^{+0.23}$
${\theta }_s$/ arcmin	${8.9}_{-3.0}^{+3.9}$		${4.27}_{-1.4}^{+0.95}$
$\alpha$	${0.77}_{-0.63}^{+0.10}$		${1.66}_{-0.74}^{+0.55}$
$\beta$	${4.79}_{-0.74}^{+0.42}$	$>5.94$	${5.78}_{-0.87}^{+1.0}$

Equations

$$Pr\left(\mathbf{\Theta }|\mathcal{D},\mathcal{M}\right)=\frac{Pr\left(\mathcal{D}|\mathbf{\Theta },\mathcal{M}\right)Pr\left(\mathbf{\Theta }|\mathcal{M}\right)}{Pr\left(\mathcal{D}|\mathcal{M}\right)},$$

$$\mathcal{L}(\mathbf{\Theta })=\mathcal{L}\left({\mathbf{d}}_{\mathrm{AMI}},{\mathbf{d}}_{\mathrm{Pl}}|\mathbf{\Theta },\mathcal{M}\right).$$

$$\mathcal{L}(\mathbf{\Theta })={\mathcal{L}}_{\mathrm{AMI}}\left({\mathbf{d}}_{\mathrm{AMI}}|\mathbf{\Theta },\mathcal{M}\right){\mathcal{L}}_{\mathrm{Pl}}\left({\mathbf{d}}_{\mathrm{Pl}}|{\mathbf{\Theta }}_{\mathrm{cl}},\mathcal{M}\right).$$

$${\mathcal{L}}_{\mathrm{AMI}}(\mathbf{\Theta })=\frac{1}{Z_D}{e}^{-\frac{1}{2}{\chi }_{\mathrm{AMI}}^2}.$$

$$\begin{array}{r}{\chi }_{\mathrm{AMI}}^2=\sum _{\nu ,{\nu }^{\prime }}{\left[{\mathbf{d}}_{\mathrm{AMI},\nu }-{\mathbf{d}}_{\mathrm{AMI},\nu }^{\mathrm{p}}(\mathbf{\Theta })\right]}^{\mathrm{T}}{\text{C}}_{\mathrm{AMI},\nu ,{\nu }^{\prime }}^{-1}\times \\ \left[{\mathbf{d}}_{\mathrm{AMI},{\nu }^{\prime }}-{\mathbf{d}}_{\mathrm{AMI},{\nu }^{\prime }}^{\mathrm{p}}(\mathbf{\Theta })\right].\end{array}$$

$${\mathbf{d}}_{\mathrm{Pl},\nu }^{\mathrm{p}}({\mathbf{\Theta }}_{\mathrm{cl}})={\mathbf{s}}_{\nu }({\mathbf{\Theta }}_{\mathrm{cl}})+{\mathbf{b}}_{\nu }+{\mathbf{n}}_{\nu }.$$

$$\begin{array}{rl}\ln \left[\frac{{\mathcal{L}}_{\mathrm{Pl}}\left({\mathbf{\Theta }}_{\mathrm{cl}}\right)}{{\mathcal{L}}_{\mathrm{Pl}}\left({\mathbf{\Theta }}_{\mathrm{cl},0}\right)}\right]=& \sum _{\nu ,{\nu }^{\prime }}{\tilde{\mathbf{d}}}_{\mathrm{Pl},\nu }^{\mathrm{p}}({\mathbf{\Theta }}_{\mathrm{cl}}{)}^{\mathrm{T}}{\text{C}}_{\mathrm{Pl},\nu ,{\nu }^{\prime }}^{-1}{\tilde{\mathbf{d}}}_{\mathrm{Pl},\nu }({\mathbf{\Theta }}_{\mathrm{cl}})\\ & -\frac{1}{2}{\tilde{\mathbf{d}}}_{\mathrm{Pl},\nu }^{\mathrm{p}}({\mathbf{\Theta }}_{\mathrm{cl}}{)}^{\mathrm{T}}{\text{C}}_{\mathrm{Pl},\nu ,{\nu }^{\prime }}^{-1}{\tilde{\mathbf{d}}}_{\mathrm{Pl},\nu }^{\mathrm{p}}({\mathbf{\Theta }}_{\mathrm{cl}}),\end{array}$$

$$y=\frac{{\sigma }_{\mathrm{T}}{k}_{\mathrm{B}}}{m_{\mathrm{e}}{c}^2}\int T_{\mathrm{e}}(r)n_{\mathrm{e}}(r)\mathrm{d}l,$$

$$y=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}\int P_{\mathrm{e}}(r)\mathrm{d}l.$$

$$P_{\mathrm{e}}(r)=\frac{P_{\mathrm{ei}}}{{\left(\frac{r}{r_{\mathrm{p}}}\right)}^{\gamma }{\left(1+{\left(\frac{r}{r_{\mathrm{p}}}\right)}^{\alpha }\right)}^{(\beta -\gamma )/\alpha }}.$$

$$Y_{\mathrm{tot}}=\frac{4\pi {\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}{P}_{\mathrm{ei}}{D}_A{\theta }_{\mathrm{s}}^3\frac{\Gamma \left(\frac{3-\gamma }{\alpha }\right)\Gamma \left(\frac{\beta -3}{\alpha }\right)}{\alpha \Gamma \left(\frac{\beta -\gamma }{\alpha }\right)}$$

Full text

††pubyear: 2019††pagerange: Sunyaev–Zel’dovich profile fitting with joint AMI-Planck analysis–Sunyaev–Zel’dovich profile fitting with joint AMI-Planck analysis

1 Introduction

With the advent of large Sunyaev-Zel’dovich (SZ) effect surveys carried out by instruments such as Planck (Planck Collaboration et al., 2016), the Atacama Cosmology Telescope (Hilton et al., 2018), and the South Pole Telescope (Bleem et al., 2015), SZ observations have the potential to become a powerful tool for constraining, for example, cosmological properties via cluster number counts. Numerical simulations show a tight, low-scatter correlation between the SZ observable, the Compton-$y$ parameter and cluster mass (e.g. da Silva et al. 2004, Nagai 2006), but recent attempts to use SZ cluster number counts for cosmological analysis have produced results in tension with other, more mature methods such as the Cosmic Microwave Background (CMB) primary anisotropies (Planck Collaboration et al., 2016). One issue is the uncertain mass-observable calibration, although another issue is the modelling of the observable itself.

Perrott et al. (2015) (hereafter P15) compared properties of 99 galaxy clusters observed in SZ with Planck and the Arcminute Microkelvin Imager (AMI) radio interferometer system. They showed that the discrepancies between the cluster parameters as constrained by AMI and Planck could be explained by the cluster gas pressure profile deviating from the profile assumed for analysis. The AMI observations were shown to be particularly sensitive to this effect when attempting to constrain the total integrated Compton-$y$ parameter due to missing angular scales. It was noted in P15 that the combination of the two instruments would be powerful for investigating the gas pressure profiles of the clusters due to the complementary angular scales measured. In this paper we explore this idea further by developing a joint Bayesian analysis pipeline which combines the data from the two instruments. We note that this pipeline could also be used with other interferometric data, for example from the Atacama Large Millimeter/submillimeter Array (ALMA) which has recently been used for SZ analysis (e.g. Kitayama et al. 2016) and in the future for the Square Kilometre Array (SKA), which will be able to observe the SZ effect in its highest frequency band (Grainge et al., 2015). The pipeline would also be easily extended to include data from other single-dish instruments such as NIKA(2) (e.g. Adam et al. 2014). This work joins a growing body of analysis combining SZ data from different instruments, sensitive to different angular scales. Recent work such as Sayers et al. (2016), Ruppin et al. (2017), and Di Mascolo et al. (2019) all combine Planck Modified Internal Linear Combination Algorithm (MILCA; Planck Collaboration et al. 2016) $y$-maps with SZ data from other instruments to jointly fit the gas pressure profile, while Romero et al. (2018) use a Planck-derived prior on the integrated Compton-$y$ parameter; however to our knowledge this is the first time that Planck frequency maps (rather than $y$-maps) have been jointly analysed with other SZ data.

The paper is arranged as follows. In Section 2 we describe the telescopes and data used in our analysis and in Section 3 we describe our analysis method. In Section 4 we describe our cluster model. In Sections 5 and  6 we verify our method using simulated clusters and in Section 7 we apply the method to a test case using real data. We anticipate further work in Section 8 and conclude in Section 9.

2 Planck and AMI telescopes

2.1 Planck satellite

The combination of the Planck satellite’s low-frequency and high-frequency instruments (LFI and HFI) provides nine frequency channels in the range 37 GHz – 857 GHz. The HFI, used for cluster analysis, has angular resolutions of $10$, $7.1$, and $5.5$ arcminutes at $100$, $143$, and $217$ GHz and $5.0$ arcminutes at each of $353$, $545$, and $857$ GHz. The Planck frequency bands correspond to two decrements, the null, and three increments in the SZ spectrum, making it particularly effective at the blind identification of galaxy clusters despite its relatively low angular resolution. See e.g. Planck Collaboration et al. 2016 for further details.

2.2 AMI

AMI (Zwart et al., 2008) is a dual-array interferometer designed for SZ studies, which is situated near Cambridge, UK. AMI consists of two arrays: the Small Array (SA), optimised for viewing arcminute-scale features, having an angular resolution of $\approx$ 3 arcmin and sensitivity to structures up to $\approx$ 10 arcmin in scale; and the Large Array (LA), with angular resolution of $\approx$ 30 arcsec, which is insensitive to SA angular scales and is used to characterise and subtract confusing radio-sources. Both arrays operate over the same frequency band with a central frequency of $\approx$ 15.5 GHz and a bandwidth of $\approx$ 5 GHz; in P15 this was divided into 6 channels but after a recent correlator upgrade the band is now divided into 4096 channels (Hickish et al., 2018), and binned down to 8 channels for analysis after radio-frequency-interference excision and calibration. The simulated data in this paper have the properties of the new correlator.

3 Joint likelihood analysis

3.1 Bayesian parameter estimation

For a model, $\mathcal{M}$ and a data vector, $\boldsymbol{\mathcal{D}}$, we can obtain the probability distributions of model parameters (also known as input or sampling parameters) $\boldsymbol{\Theta}$ conditioned on $\mathcal{M}$ and $\boldsymbol{\mathcal{D}}$ using Bayes’ theorem:

[TABLE]

where $Pr(\boldsymbol{\Theta}|\boldsymbol{\mathcal{D}},\mathcal{M})\equiv\mathcal{P}(\boldsymbol{\Theta})$ is the posterior distribution of the model parameter set, $Pr(\boldsymbol{\mathcal{D}}|\boldsymbol{\Theta},\mathcal{M})\equiv\mathcal{L}(\boldsymbol{\Theta})$ is the likelihood function for the data, $Pr(\boldsymbol{\Theta}|\mathcal{M})\equiv\pi(\boldsymbol{\Theta})$ is the prior probability distribution for the model parameter set, and $Pr(\boldsymbol{\mathcal{D}}|\mathcal{M})\equiv\mathcal{Z}(\boldsymbol{\mathcal{D}})$ is the Bayesian evidence of the data given a model $\mathcal{M}$. In this paper we will be interested in the posterior distributions of the sampling parameters rather than the evidence, which would be used for model comparison. We use the nested sampling algorithm MultiNest (Feroz et al., 2009) to calculate our posteriors.

3.2 Model parameters

The model parameters can be split into two subsets (which are assumed to be independent of one another): cluster parameters $\boldsymbol{\Theta}_{\rm cl}$ and radio-source or ‘nuisance’ parameters $\boldsymbol{\Theta}_{\rm rs}$. The cluster model parameters are relevant to both AMI and Planck data, and are detailed with their associated priors $\pi(\boldsymbol{\Theta}_{\rm cl})$ in Section 4. $\boldsymbol{\Theta}_{\rm rs}$ are relevant only for AMI data, since they are used to model the radio-source contamination of the SZ cluster signal recorded by the SA, based on values measured with the LA. More information on the prior distributions used for $\boldsymbol{\Theta}_{\rm rs}$ can be found in Section 4.3 of P15.

3.3 Joint likelihood function

If one has an AMI dataset $\boldsymbol{d}_{\mathrm{AMI}}$ and a Planck dataset $\boldsymbol{d}_{\mathrm{Pl}}$, then the joint likelihood function for the data is given by

[TABLE]

In this analysis we treat $\boldsymbol{d}_{\mathrm{AMI}}$ and $\boldsymbol{d}_{\mathrm{Pl}}$ as being independent (see Section 5.3 for justification), and since the Planck predicted data only rely on the cluster parameters we can write

[TABLE]

3.3.1 AMI likelihood function

The AMI likelihood calculation is detailed in Feroz et al. (2009) (herafter F09). Briefly, the AMI likelihood function $\mathcal{L}_{\mathrm{AMI}}\left(\boldsymbol{d}_{\mathrm{AMI}}|\boldsymbol{\Theta},\mathcal{M}\right)\equiv\mathcal{L}_{\mathrm{AMI}}(\boldsymbol{\Theta})$ is given by

[TABLE]

Here $\chi_{\mathrm{AMI}}^{2}$ is a measure of the goodness-of-fit between the real and modelled data and can be expressed as

[TABLE]

In this expression $\boldsymbol{d}_{\mathrm{AMI},\nu}$ are the data observed by AMI at frequency $\nu$, and $\boldsymbol{d}_{\mathrm{AMI},\nu}^{\rm p}(\boldsymbol{\Theta})$ are the predicted data generated by the model at the same frequency. $\textbf{{C}}_{\mathrm{AMI},\nu,\nu^{\prime}}$ is the theoretical covariance matrix for the AMI likelihood, which includes primordial CMB and source confusion noise as described in Hobson and Maisinger (2002) and F09 (Section 5.3). Source confusion noise allows for the remaining radio-sources with flux densities below the flux limit $S_{\rm lim}$ that the LA can subtract down to. The instrumental noise is estimated from the scatter of the visibility measurements within an observation. Referring back to equation 4, $Z_{D}$ is a normalisation constant given by $(2\pi)^{D/2}|\textbf{{C}}_{\mathrm{AMI}}|^{1/2}$ where $D$ is the length of $\boldsymbol{d}_{\mathrm{AMI}}$ (i.e. the combined vector of data from all frequencies).

3.3.2 Planck likelihood function

To calculate the Planck likelihood, we use a version of the PowellSnakes (PwS; Carvalho et al. 2009 and Carvalho et al. 2012) Bayesian detection algorithm designed for detecting galaxy clusters in Planck data, adapted to operate on a previously determined position rather than to conduct a blind search. PwS treats the data observed by Planck as a superposition of background sky emission (including foreground emission and primordial CMB) $\boldsymbol{b}_{\nu}$, instrumental noise $\boldsymbol{n}_{\nu}$, and signal from the (cluster) source $\boldsymbol{s}_{\nu}$. The model for the predicted data vector is thus

[TABLE]

PwS works with patches of sky sufficiently small such that it can be assumed the noise contributions are statistically homogeneous. In this limit it is more convenient to work in Fourier space, as the Fourier modes are uncorrelated. It also assumes that the noise contributions are Gaussian, which is very accurate in the case of instrumental noise and the primordial CMB but may be more questionable for Galactic emission – the deviations from Gaussianity of $\boldsymbol{b}_{\nu}$ are discussed in Section 4.3 of Carvalho et al. 2012. Since PwS is a detection algorithm, it calculates the ratio of the likelihood of detecting a cluster parameterised by $\boldsymbol{\Theta}_{\rm cl}$ and the likelihood of the data with no cluster signal ($\boldsymbol{s}_{\nu}(\boldsymbol{\Theta}_{\rm cl,0}$) = 0). Thus the log-likelihood ratio of the Fourier transformed quantities is

[TABLE]

where tildes denote the Fourier transform of a quantity, and $\textbf{{C}}_{\mathrm{Pl},\nu,\nu^{\prime}}$ is the covariance matrix of the data in Fourier space.

3.3.3 Joint likelihood analysis hyperparameters

Ideally when making a joint inference from two different datasets, likelihood hyperparameters would be used (see Lahav et al. 2000, Hobson et al. 2002, and Ma and Berndsen 2014) so that the relative weighting of the two likelihoods is treated in a Bayesian way. This allows meaningful results to be extracted when the datasets are not in good agreement, which could be due to, for example, systematic bias or an incorrect model. However the log-ratio given by equation 7 is not a probability density due to the fact that it is missing a normalisation factor proportional to $\tilde{\boldsymbol{d}}_{\mathrm{Pl},\nu}(\boldsymbol{\Theta}_{\rm cl})^{\rm T}\textbf{{C}}_{\mathrm{Pl},\nu,\nu^{\prime}}^{-1}\tilde{\boldsymbol{d}}_{\mathrm{Pl},\nu}(\boldsymbol{\Theta}_{\rm cl})$. These hyperparameters affect the normalisation factor of the likelihood (it becomes a function of them), thus the PwS algorithm is incompatible with their use as it does not include the normalisation term in the likelihood calculations (a more detailed explanation is given in Chapter 8 of Javid 2019). We therefore do not include likelihood hyperparameters in our analysis. This should not be problematic when applying our method unless the data have been incorrectly analysed or the cluster model used is not flexible enough to describe both of the datasets.

4 Cluster model

Since AMI and Planck observe the SZ effect caused by the electron gas in clusters (see e.g. Birkinshaw 1999), both telescopes measure signals proportional to the Comptonisation parameter $y$,

[TABLE]

where $k_{\rm B}$ is the Boltzmann constant, $m_{\rm e}$ is the rest mass of an electron, $c$ is the speed of light, and $\sigma_{\rm T}$ is the Thomson scattering cross-section. $T_{\rm e}(r)$ and $n_{\rm e}(r)$ are respectively the electron temperature and number density in the intra-cluster medium, as a function of radius from the centre of the cluster ($r$), and the integral is over the line of sight. If an ideal gas equation of state is assumed for the electron gas then in terms of the electron pressure $P_{\mathrm{e}}(r)$, the Comptonisation parameter is given by

[TABLE]

The cluster model considered here is used to calculate a ‘map’ of $y$ by evaluating equation 9 at different points on the plane of the sky. It assumes a spherically symmetric, generalised-Navarro-Frenk-White (GNFW, Nagai et al. 2007) profile to model the electron pressure

[TABLE]

$P_{\rm ei}$ is an overall pressure normalisation factor and $r_{\rm p}$ is a characteristic radius.

The parameters $\alpha$, $\beta$ and $\gamma$ describe the slope of the pressure profile at $r\approx r_{\rm p}$, $r\gg r_{\rm p}$ and $r\ll r_{\rm p}$ respectively. Our model parameterises a cluster in terms of observational (rather than physical) quantities: $\theta_{\rm s}$, the characteristic angular scale corresponding to $r_{\rm p}$ ($\theta_{\rm s}=r_{\rm p}/D_{A}$ where $D_{A}$ is angular diameter distance) and $Y_{\rm tot}$, the total integrated Comptonisation parameter of the cluster, given by

[TABLE]

where $\Gamma$ is the gamma function.

Our model therefore has input parameters $\boldsymbol{\Theta}_{\rm cl}=\left(x_{0},y_{0},Y_{\rm tot},\theta_{\rm s},\alpha,\beta,\gamma\right)$, where $x_{0}$ and $y_{0}$ are the cluster centre offsets from the designated central sky coordinate.

A set of ‘universal’ pressure profile (UPP) GNFW shape parameter values were derived in Arnaud et al. (2010) as the best fit to a sample of clusters from REXCESS (observed with XMM-Newton, Böhringer et al. 2007). These are $(\gamma,\alpha,\beta,c_{500})$ = (0.3081, 1.0510, 5.4905, 1.177) and are often used as a fixed standard cluster profile. In this analysis we will not restrict our model profiles to the UPP case.

The GNFW profile extends to infinity; in practice some cut-off radius for the $y$-map must therefore be defined when implementing this model. A frequent choice, used for example in the analysis of Planck data, is to cut off at $\theta=5\theta_{500}=5\theta_{\rm s}c_{500}$. For the case of the UPP, this implies that a spherical integral to $5\theta_{500}$ gives $Y_{\mathrm{sph},5\theta_{500}}=0.96\,Y_{\rm tot}$, or a cylindrical integral gives $Y_{\mathrm{cyl},5\theta_{500}}=0.97\,Y_{\rm tot}$ (with the line-of-sight integral extending to infinity within a radius of $5\theta_{500}$ on the sky). In the case of arbitrary values of $(\gamma,\alpha,\beta)$, this fraction can change significantly; plus in the general case $c_{500}$ is not necessarily known so a different cut-off radius must be defined which we denote $\theta_{\mathrm{lim}}$. We choose to define $\theta_{\mathrm{lim}}$ for an arbitrary profile by the radius at which $Y_{\mathrm{sph},\theta_{\mathrm{lim}}}=0.95\,Y_{\rm tot}$ (found via a numerical root-finder) by analogy with the UPP. For some profiles which fall off more slowly with radius this becomes prohibitively large and we therefore impose a maximum radius of $10\,\theta_{\rm s}$. We have verified that in these cases the Comptonization parameter integrated over the line-of-sight at $10\,\theta_{\rm s}$ is $<0.1$% of the value at the centre, i.e. we are not cutting off substantial cluster signal in the outskirts.

We then translate the $y$-map to signal on the sky at each frequency using the non-relativistic approximation (Zeldovich and Sunyaev, 1969). Recently it has been shown that relativistic corrections may be significant for high-temperature clusters at Planck frequencies (e.g. Remazeilles et al. 2019); we note that relativistic effects reduce the Planck SZ signal, so if corrections were applied the effect would be to increase the discrepancy between AMI and Planck parameter constraints shown in P15.

4.1 Priors

In P15 we used a position prior based on the Planck position and error, and a joint ellisoidal prior on $Y_{\rm tot}$ and $\theta_{\rm s}$ based on the population of clusters detected by Planck. Here we assign wide, non-informative, independent priors to $x_{0}$, $y_{0}$, $Y_{\rm tot}$ and $\theta_{\rm s}$ (see Table 1), to explore how much the combination of the two datasets can constrain the parameters.

In the standard Planck analysis and for the AMI data analysed in P15, the GNFW shape parameter values were fixed to the UPP values. In this analysis we will both simulate and analyse clusters with non-UPP profiles and explore the constraints on $\alpha$ and $\beta$ produced by our joint dataset. The priors used for $\alpha$ and $\beta$ vary throughout our analysis (Table 1); they are either fixed at some specific value (delta function priors) or allowed to vary uniformly (uniform priors on a bounded domain).

5 Cluster simulations

For all simulations the $y$ map of a single cluster is generated with ten distinct noise realisations as follows. Firstly, ten CMB realisations are created by sampling primordial CMB noise from an empirical power-law distribution (Hinshaw et al., 2013) and distributing at random positions on the sky. For each instrument, we add further realisations of the relevant sources of noise, as follows.

5.1 Planck cluster simulations

We construct a Planck all-sky foreground and thermal noise map by adding one noise realisation to the CO, FIRB, free-free, spinning dust, synchrotron and thermal dust emission maps simulated using the Planck sky model, all taken from the Planck Legacy Archive111https://pla.esac.esa.int/#home using the FFP8 release. We note that we do not include point sources since as of this data release they were not separated into strong and weak point source maps, and for cluster analysis on real data the strong point sources would be masked. We also choose not to add the thermal and kinetic SZ emission maps as we wish to see how much information can be extracted from the data in an ideal situation of a single, isolated cluster.

We then randomly select ten patch centres on the sky, with the constraint that the patch centres are above $\delta=0^{\circ}$ (to satisfy AMI’s observing limits) and that all of the 20*∘* square patch is outside of the Planck 20% Galactic plane mask (in which 20% of the sky is masked). We cut patches from the all-sky foreground $+$ noise map at these coordinates and add the resulting patch maps to the CMB and cluster maps to produce the final Planck simulations. Each patch map therefore contains different thermal noise, foreground emission and CMB but the same cluster.

5.2 AMI cluster simulations

The AMI simulations are constructed by adding the cluster and CMB maps to confusion noise realisations, created using the 10C source counts given in AMI Consortium et al. (2011) up to a maximum flux density of 360 $\upmu$Jy (i.e. assuming sources above $4\times$ a typical AMI-LA noise limit of 90 $\upmu$Jy/beam have been detected and removed). A ten-hour mock AMI-SA observation is performed at the ten patch centres used for the Planck simulations, using the in-house package Profile (see e.g. Grainge et al. 2002). Instrumental noise is also added to the mock observations, giving a total noise level of $\approx$ 120 $\upmu$Jy/beam on the map.

Similarly to the Planck simulations, each AMI simulation contains different thermal noise, confusion noise and CMB and the same cluster. Each AMI simulation corresponds directly to a Planck simulation which has the same observing centre and CMB realisation.

5.3 Testing the independence of the AMI and Planck datasets

In Section 3.3 we made the assumption that $\boldsymbol{d}_{\mathrm{AMI},\nu}^{\rm p}$ and $\boldsymbol{d}_{\mathrm{Pl},\nu}^{\rm p}$ are not correlated with each other, so that the likelihoods for the two datasets can be separated. The instrumental noises associated with each telescope can safely be assumed to be independent. Due to the telescopes operating at different angular scales and frequencies, the confusion noise present in AMI data and foreground emission present in Planck data are assumed to be independent of one another. A similar argument can be applied for primordial CMB noise, nevertheless we carried out a simple test to see if this is the case. For a given set of cluster parameters, we ran the joint analysis on Planck and AMI datasets which had different CMB realisations to one another. We found that the resultant parameter constraints were not significantly different to the results obtained using AMI and Planck data which had the same CMB realisations as one another (Figure 1). We thus concluded that the covariance between the datasets introduced by the common CMB background was negligible.

6 Cluster simulation results

In the following we apply the joint analysis to clusters simulated as described in Section 5, and compare results with analyses which use (the same) AMI or Planck data alone.

We consider simulations generated using three different sets of profile shape parameters: the UPP, the parameters fitted to a stacked Planck dataset in Planck Collaboration et al. (2013) (‘Planck Intermediate Profile’, PIP), and the parameters fitted to the cluster RXC J2319.6-7313 in Arnaud et al. (2010) which are the most different in the sample to the UPP (‘REXCESS Extreme Profile’, REP). These values are listed in Table 2. We consider a ‘low’ and a ‘high’ signal-to-noise ratio (SNR) cluster, which correspond to input values of $(Y_{\rm tot},\theta_{\rm s})=(0.001\,\mathrm{arcmin}^{2},2\,\mathrm{arcmin})$ and $(Y_{\rm tot},\theta_{\rm s})=(0.007\,\mathrm{arcmin}^{2},8\,\mathrm{arcmin})$ respectively. We note that ‘low’ and ‘high’ SNR are in relation to the Planck simulations rather than the AMI simulations, where they are both well detected. We then analyse the ten different noise realisations for each cluster using the priors given in Table 1 and plot the resulting posterior distributions using GetDist222http://getdist.readthedocs.io/en/latest/..

6.1 Analysis with fixed profile parameters

We firstly analyse the high- and low-SNR UPP simulations using $(\gamma,\alpha,\beta)$ fixed to the correct, input values. Figs. 2 and 3 show the results of the simulation sets using AMI data only, Planck data only, and the two datasets combined. Each two-dimensional marginalised posterior distribution plot shows the 68% confidence contours for the 10 different noise realisations, with the input values marked with black stars.

In the case of the positional parameters, it is clear that the higher angular resolution of the AMI data means that it drives the posterior inferences, although the addition of the Planck data does improve the constraints slightly in the high-SNR case. This result is consistent for all following analyses, and we do not consider the $x_{0}$ and $y_{0}$ constraints any further.

Both AMI and Planck constraints are degenerate in the $\theta_{\rm s}$/$Y_{\rm tot}$ plane and previous results based on AMI and Planck data have relied on ancillary data to reduce this degeneracy. For AMI, a correlated prior in $\theta_{\rm s}$ and $Y_{\rm tot}$ based on a simulated population of clusters injected into and recovered from Planck data was used (see P15). For Planck, a prior on $\theta_{\rm s}$ has been used to ‘slice’ the $\theta_{500}$/$Y(r_{500})$ posterior constraint based on either an X-ray measurement or a mass-observable scaling relationship (see, e.g. Planck Collaboration et al. 2016), both relying on the assumption of the ‘universal’ $c_{500}$. The combination of AMI and Planck data removes the need for these ancillary priors and produces a much tighter constraint on both $\theta_{\rm s}$ and $Y_{\rm tot}$ since the degeneracy directions are different. This is most striking in the case of the low SNR cluster, but is also evident in the case of the high SNR cluster.

We next analyse the PIP and REP simulations keeping $(\gamma,\alpha,\beta)$ fixed to the UPP values in the analysis, i.e. we now have a mismatch between the cluster profiles used to produce and analyse the simulations. The two-dimensional posterior constraints on $Y_{\rm tot}$ and $\theta_{\rm s}$ in this case are shown in Figs. 4 and 5. In the case of the PIP simulations, all of the constraints are offset from the true position; although the direction of the offset is mostly in $\theta_{\rm s}$ so that the one-dimensional marginal constraint on $Y_{\rm tot}$ is not too badly offset it is clear that any method to reduce the degeneracy by slicing the posterior will be problematic. In the case of the REP simulations, the low-SNR cluster constraints are not too badly affected by the profile mismatch; this is because the change is to the inner part of the profile which is not well-resolved by either instrument. The high-SNR constraints however are significantly biased both in the two-dimensional plane and in the one-dimensional $Y_{\rm tot}$ plane.

6.2 Variable shape parameter analysis

We next consider the same set of simulations described in Section 6, but now allowing the GNFW shape parameters $\alpha$ and $\beta$ to vary in the analysis. We assign the uniform priors stated in Table 1 to $\alpha$ and $\beta$.

Figs. 6 and 7 show two-dimensional posterior constraints for a selection of parameter pairs for the UPP low- and high-SNR cluster simulations respectively. The two sets of simulations share some common features, as follows. The constraints on $Y_{\rm tot}$ are mostly driven by the Planck data, since the lower-resolution Planck data are most suited to measuring the total cluster signal while the interferometric AMI data relies on extrapolations to larger angular scales than are measured; the joint constraints are generally tighter than the Planck-only constraints and appear unbiased. There is a large degeneracy between $\theta_{\rm s}$ and $\beta$ which is inherent to the GNFW model, since decreasing $\beta$ makes the cluster amplitude fall off more slowly with radius which can also be achieved by increasing $\theta_{\rm s}$.

In the case of the low SNR cluster, the joint analysis successfully constrains $Y_{\rm tot}$ and puts an upper limit on $\theta_{\rm s}$ and a lower limit on $\beta$. $\alpha$ is fairly unconstrained since it affects the profile on the scale of $\theta_{\rm s}$ = 2 arcmin which is not well-resolved by either instrument.

In the case of the high SNR cluster, $\alpha$ can be constrained by AMI alone. The joint analysis places a lower limit on $\theta_{\rm s}$ but $\beta$ is unconstrained due to the $\theta_{\rm s}$/$\beta$ degeneracy.

We do not show all of the constraints produced for the PIP and REP simulations as they are qualitatively similar, with the following exceptions. Since $\beta$ is lower for the PIP it effectively makes the cluster much more extended. The significance of the high SNR AMI detection becomes much lower and with $\alpha$ and $\beta$ allowed to vary all of the parameters are essentially unconstrained, so the joint constraints become driven by the Planck data. $Y_{\rm tot}$ is well-constrained; a lower limit can be put on $\theta_{\rm s}$ and an upper limit on $\beta$, while $\alpha$ is unconstrained. In the low SNR case, only $Y_{\rm tot}$ is well-constrained since although the cluster is more extended it is also much less bright so there is not enough signal-to-noise to constrain $\alpha$ and $\beta$. We show the joint constraints for this case in Fig. 8 as a ‘worst-case’ scenario.

In the high SNR, REP case, $\alpha$ is very tightly constrained and $\beta$ is better-constrained by the AMI data alone. This is because the lower $\alpha$ value of the profile causes it to fall off more sharply with radius, putting more signal on AMI scales. In the low SNR, REP case we obtain tight constraints on $Y_{\rm tot}$ and $\alpha$, a tight upper limit on $\theta_{\rm s}$ and a lower limit on $\beta$. We show the joint constraints for this case in Fig. 8 as a ‘best-case’ scenario. We note that for the REP analysis we are fixing $\gamma$ to the incorrect, UPP value; we also varied $\alpha$ and $\beta$ while fixing $\gamma$ to the correct, REP value which had little impact on the parameter constraints.

Overall we conclude that where nothing is known a-priori about the pressure profile of the cluster (other than that it follows a GNFW shape), we can use the combination of AMI and Planck data to successfully constrain $Y_{\rm tot}$; often constrain $\alpha$; and sometimes put limits on $\beta$ and $\theta_{\rm s}$ depending on how well-resolved the cluster is (which depends both on $\theta_{\rm s}$ and on $\alpha$ and $\beta$).

7 Application of joint analysis to real cluster data

As a test case, we apply the joint analysis to PSZ2 G063.80+11.42 ($z=0.426$; $M_{500,\mathrm{SZ}}=6.2\times 10^{14}M_{\sun}$) from the P15 sample. The cluster is well-detected by both Planck (PwS signal-to-noise-ratio = 6.5) and AMI (Bayesian detection significance = 3.8) and there is a significant offset between the AMI and Planck posterior constraints on $Y_{\rm tot}$ and $\theta_{\rm s}$. The AMI radio-source environment is relatively clean. We reobserved the cluster with AMI to benefit from the improved performance of the new correlator, and use the DX11d data release with strong point sources masked (Planck Collaboration, priv. comm.) for the Planck data.

We firstly run the AMI and Planck analyses separately using the priors given in Table 1 (assigning delta priors to the GNFW shape parameters at UPP values), to confirm the discrepancy. For comparison we also ran the AMI analysis with the P15 priors; these three sets of posterior chains are shown in Fig. 9. We note that the tighter positional constraint in the latter case is due to a degeneracy with a radio source flux near one edge of the cluster; when the wider positional prior is used the source flux is allowed to increase, broadening the cluster decrement and shifting the position of the cluster. This has little effect on the $\theta_{\rm s}$/$Y_{\rm tot}$ constraint. We see that the discrepancy is confirmed with the newer AMI data (and the AMI Bayesian detection significance increases to 6.0); the Planck $\theta_{\rm s}$/$Y_{\rm tot}$ constraint lies significantly above that produced by AMI.

We now analyse the PSZ2 G063.80+11.42 data while allowing $\alpha$ and $\beta$ to vary. Figure 10 shows the resulting posterior distributions for the AMI-only, Planck-only and joint analysis methods. Similarly to the simulations, we see that $Y_{\rm tot}$ is well constrained by the joint analysis and the constraint is completely driven by the Planck data. The addition of AMI data improves the constraint on $\theta_{\rm s}$ from an upper limit to a true constraint, and both $\alpha$ and $\beta$ are constrained, although not tightly. The one-dimensional parameter constraints are summarized in Table 3; $\beta$ agrees with the UPP value, while a higher value of $\alpha$ is favoured but only at the $\approx$ $1\sigma$ level. We perform one further analysis on the AMI data, fixing $\alpha$ to the fitted value from the joint analysis, and leaving $\beta$ fixed to the UPP value. The posterior constraints on $\theta_{\rm s}$ and $Y_{\rm tot}$ in this case are shown in comparison with the AMI-only UPP analysis and the Planck analysis varying $\alpha$ and $\beta$; it can be seen that the AMI posterior shifts in the correct direction to overlap with Planck, confirming that this value of $\alpha$ gives better agreement between the two instruments (Fig. 10, upper right-hand corner). No X-ray observations of this cluster are available to compare X-ray-derived constraints on the profile parameters; for our future AMI-Planck sample (see Section 8) comparison with X-ray-derived parameters will be informative and complementary across a different range of angular scales.

8 Future work

Along with the sample of 99 clusters from P15, clusters which were previously excluded from the AMI sample due to difficult radio source environments are currently being reobserved with the new correlator; the superior dynamic range of the new instrument allows us to cope better with these environments and successfully extract cluster parameters. We will analyse all AMI detections with the joint pipeline, giving us a large cluster sample to probe deviations from the UPP and consider the impact this may have on the Planck cluster number counts.

We note that Javid et al. (2019) found that realistic radio source environments could bias the recovery of cluster parameters from AMI data. This issue and its effect on the recovery of the pressure profile parameters will be investigated in conjunction with the analysis of the larger AMI-Planck sample.

With the larger sample we will be able to compare our SZ-derived profile parameters to X-ray-derived parameters. This will allow us to test for any systematic differences and combine information across a broader range of angular scales.

Our pipeline is also simply extensible to the use of any physical model which uses a GNFW profile for the gas pressure, e.g. the models proposed in Olamaie et al. (2012), Javid et al. (2018) and Javid et al. (2019). We also plan to implement a non-parametric model such as that proposed in Olamaie et al. (2018).

9 Conclusions

We have developed a joint likelihood function for SZ data obtained from the Planck satellite and the Arcminute Microkelvin Imager (AMI) radio interferometer system, in order to compare inferences obtained using it with those from the individual likelihoods. The method could apply to any combination of Planck and interferometric data from one or more telescopes. 2. 2.

We generated simulations of clusters using an observational model similar to the one used in Perrott et al. (2015), using gas pressure profile shape parameter values taken from either the UPP (Arnaud et al. 2010, ‘universal’ pressure profile) or two other realistic variations. We considered both a smaller angular size, fainter cluster and a larger angular size, brighter cluster. From looking at the resulting posterior distributions we found the following:

• •

When simulating and analysing the clusters using the model with UPP parameters the joint analysis greatly reduced the degeneracy in $Y_{\rm tot}$-$\theta_{\rm s}$ shown in the individual AMI and Planck analyses, due to the different degeneracy directions for the individual datasets. The improvement on the parameter constraints for the joint analysis is particularly prominent in the small angular size cases. Thus when the profile shape of a cluster is known a-priori, the combination of the two datasets provides accurate, precise constraints on the cluster parameters with no need for external information to reduce the $Y_{\rm tot}$-$\theta_{\rm s}$ degeneracy.

• •

When simulated clusters created using non-UPP profiles are analysed with the profile shape assumed to be UPP, the $Y_{\rm tot}$-$\theta_{\rm s}$ constraints are biased away from the true value. This occurs in the individual AMI and Planck datasets and is particularly problematic in the joint analysis, where a tight, significantly biased constraint is produced.

• •

When allowing the shape parameters to vary in the Bayesian analysis, we generally found that for all the clusters $Y_{\rm tot}$ was well-constrained and unbiased; Planck drove the constraint and the joint analyses improved the constraint slightly.

• •

Furthermore dependent on how well resolved the clusters are, the shape parameter $\alpha$ can often be constrained. Due to the strong $\beta$-$\theta_{\rm s}$ degeneracy these parameters are more difficult to constrain and it is usually only possible to place limits on them. 3. 3.

Finally, we applied the joint analysis to real data for the cluster PSZ2 G063.80+11.42 which is part of the sample of 99 clusters considered in P15. We confirmed the discrepancy in $Y_{\rm tot}$ and $\theta_{\rm s}$ estimates when using updated AMI data and resolved it by allowing $\alpha$ and $\beta$ to vary. Using the joint analysis we could constrain $Y_{\rm tot}$ and $\theta_{\rm s}$ well and place loose constraints on $\alpha$ and $\beta$, finding that a slightly higher value of $\alpha$ than the UPP value was preferred, while the constraint on $\beta$ was consistent with the UPP value. 4. 4.

We plan to apply our method to all the 99 clusters of the P15 sample, plus clusters currently being reobserved by AMI, to investigate deviations from the UPP and possible impact on the Planck cluster number counts. We also plan to extend our method to incorporate different physical and non-parametric cluster models.

Acknowledgements

We thank the staff of the Mullard Radio Astronomy Observatory for their invaluable assistance in the commissioning and operation of AMI, which is supported by Cambridge and Oxford Universities. This work is based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. This work was performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing (HPC) Service (http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England and funding from the Science and Technology Facilities Council. The authors would like to thank Stuart Rankin from HPC and Greg Willatt and David Titterington from Cavendish Astrophysics for computing assistance. KJ and PJE acknowledge Science and Technology Facilities Council studentships. YCP acknowledges support from a Trinity College JRF and Rutherford Discovery Fellowship.