Free-form modelling of galaxy clusters: a Bayesian and data-driven approach

TL;DR

This paper introduces a Bayesian, data-driven free-form modelling method for galaxy clusters that automatically determines the complexity of physical property reconstructions from observational data, improving flexibility over traditional parametric models.

Contribution

The authors develop a novel free-form, Bayesian approach using control points to model galaxy cluster properties without assuming specific functional forms, adaptable to data complexity.

Findings

Successfully modeled electron pressure profiles from simulated and real SZ data.

Constrained total Comptonisation parameter and gas extent in galaxy clusters.

Demonstrated potential for detecting clusters in blind SZ surveys.

Abstract

A new method is presented for modelling the physical properties of galaxy clusters. Our technique moves away from the traditional approach of assuming specific parameterised functional forms for the variation of physical quantities within the cluster, and instead allows for a 'free-form' reconstruction, but one for which the level of complexity is determined automatically by the observational data and may depend on position within the cluster. This is achieved by representing each independent cluster property as some interpolating or approximating function that is specified by a set of control points, or 'nodes', for which the number of nodes, together with their positions and amplitudes, are allowed to vary and are inferred in a Bayesian manner from the data. We illustrate our nodal approach in the case of a spherical cluster by modelling the electron pressure profile Pe(r) in analyses both of simulated Sunyaev-Zel'dovich (SZ) data from the Arcminute MicroKelvin Imager (AMI) and of real AMI observations of the cluster MACS J0744+3927 in the CLASH sample. We demonstrate that one may indeed determine the complexity supported by the data in the reconstructed Pe(r), and that one may constrain two very important quantities in such an analysis: the cluster total volume integrated Comptonisation parameter (Ytot) and the extent of the gas distribution in the cluster (rmax). The approach is also well-suited to detecting clusters in blind SZ surveys.

Tables

Table 1. Table 1: Jeffreys’ scale for interpreting PORs (or Bayes factors). As ln ⁡ 𝒫 i ​ j = − ln ⁡ 𝒫 j ​ i subscript 𝒫 𝑖 𝑗 subscript 𝒫 𝑗 𝑖 \ln\mathcal{P}_{ij}=-\ln\mathcal{P}_{ji} negative values imply reversed model favouring.

$\ln {𝒫}_{ij}$	Favouring of ${ℳ}_j$ over ${ℳ}_i$
	None
	Slight
	Significant
$5.0\le \ln {𝒫}_{ij}$	Decisive

Equations

$$\mathrm{Pr}(\mathbf{\Theta }|\mathbf{D},\mathcal{M})=\frac{\mathrm{Pr}(\mathbf{D}|\mathbf{\Theta },\mathcal{M})\mathrm{Pr}(\mathbf{\Theta }|\mathcal{M})}{\mathrm{Pr}(\mathbf{D}|\mathcal{M})}\equiv \frac{\mathcal{L}(\mathbf{\Theta })\pi (\mathbf{\Theta })}{\mathcal{Z}},$$

$$\mathrm{Pr}(\mathcal{M}|\mathbf{D})=\frac{\mathrm{Pr}(\mathbf{D}|\mathcal{M})\mathrm{Pr}(\mathcal{M})}{\mathrm{Pr}(\mathbf{D})}\equiv \frac{\mathcal{Z}{\pi }_{\mathcal{M}}}{\mathrm{Pr}(\mathbf{D})},$$

$${\mathcal{P}}_{ij}\equiv \frac{\mathrm{Pr}({\mathcal{M}}_j|\mathbf{D})}{\mathrm{Pr}({\mathcal{M}}_i|\mathbf{D})}=\frac{Z_j{\pi }_{{\mathcal{M}}_j}}{Z_i{\pi }_{{\mathcal{M}}_i}}={\mathcal{B}}_{ij}\frac{{\pi }_{{\mathcal{M}}_j}}{{\pi }_{{\mathcal{M}}_i}},$$

$$\mathcal{Z}=\int \mathcal{L}(\mathbf{\Theta })\pi (\mathbf{\Theta }){\mathrm{d}}^D\mathbf{\Theta },$$

$$\delta I_{\nu}=T_{\rm CMB}yf(\nu)\frac{\partial B_{\nu}}{\partial T}\Big{|}_{T=T_{\rm CMB}}\,.$$

$$y=\frac{{\sigma }_T}{m_{\mathrm{e}}{c}^2}{\int }_{-\infty }^{+\infty }{P}_{\mathrm{e}}(r)\mathrm{d}l,$$

$$Y_{\mathrm{cyl}}(R)$$

$$Y_{\mathrm{sph}}(r)$$

$${\tilde{I}}_{\nu }(\mathbf{u})=\int A_{\nu }(\mathbf{x})I_{\nu }(\mathbf{x})\exp (2\pi i\mathbf{u}\cdot \mathbf{x})\mathrm{d}\mathbf{x},$$

$$V_{\nu }(\mathbf{u})={\tilde{I}}_{\nu }(\mathbf{u})+N_{\nu }(\mathbf{u}),$$

$$\mathcal{L}(\mbox{\boldmath$\Theta$})\propto\exp\left(-{\textstyle\frac{1}{2}}\chi^{2}\right)\,,$$

$$\chi^{2}=\sum_{\nu,\nu^{\prime}}(\mbox{\boldmath$D$}_{\nu}-\mbox{\boldmath$D$}_{\nu}^{\rm p})^{T}(\mbox{\boldmath$C$}_{\nu,\nu^{\prime}})^{-1}(\mbox{\boldmath$D$}_{\nu^{\prime}}-\mbox{\boldmath$D$}_{\nu^{\prime}}^{\rm p})\,,$$

$$\mbox{\boldmath$\Theta$}_{\rm c}\equiv(x,y,r_{0},P_{0},\ldots,r_{N-1},P_{N-1})\,,$$

$$\pi(\mbox{\boldmath$\Theta$}_{\rm c})=\pi(x)\,\pi(y)\,\pi(r_{0},r_{1},\ldots,r_{N-1})\,\prod_{n=0}^{N-1}\pi(P_{n})\,.$$

Full text

††pagerange: LABEL:firstpage–Free-form modelling of galaxy clusters: a Bayesian and data-driven approach††pubyear: 2017

1 Introduction

Determining the properties of clusters of galaxies, such as their total and baryonic mass, has the potential to provide an independent tool for constraining the parameters of the $\Lambda\rm{CDM}$ model, since cluster population properties are sensitive to several cosmological parameters, most notably $\sigma_{8}$ (see e.g. Sievers et al. 2013, Planck Collaboration et al. 2016, and de Haan et al. 2016). The challenge lies, however, in obtaining a robust estimate of the cluster masses as these are not directly observable. The mass distribution within a cluster is usually measured using a variety of observational methods, including X-ray (see e.g. Vikhlinin et al. 2009, Ettori 2013, Mantz et al. 2014 and Olamaie et al. 2015), the Sunyaev–Zel’dovich (SZ) effect (see e.g. Planck Collaboration et al. 2013, Schammel et al. 2013, Hasselfield et al. 2013, Perrott et al. 2015, Rumsey et al. 2016, and Rodríguez-Gonzálvez et al. 2017), and gravitational lensing (see e.g. Rozo et al. 2010, Hoekstra et al. 2013, Rozo et al. 2014 and Battaglia et al. 2016).

Each of these approaches relies on developing some method for determining the cluster’s mass (distribution) from its observable properties, namely its distributions of X-ray surface brightness, SZ Comptonisation parameter, and weak-lensing shear distribution. This is usually achieved by modelling the physical properties through the cluster in terms of some specific parameterised functional forms. Typical examples include assuming an NFW profile (Navarro, Frenk, & White, 1996; Navarro, Frenk & White, 1997) for the dark matter density distribution, a $\beta$-model (Cavaliere & Fusco-Femiano, 1976, 1978) for the gas density, or a generalised NFW (GNFW) profile (Nagai, Kravtsov & Vikhlinin, 2007) for the gas pressure. The cluster mass (distribution) is then usually calculated under the standard assumption of a spherical cluster model obeying hydrostatic equilibrium (HSE) and/or some scaling relationships.

Even for physically-based cluster models (see e.g. Olamaie et al. 2012, 2013), there still remains considerable uncertainty regarding the appropriate form one should assume for the radial variation of cluster properties, and this can lead to different cluster mass estimates (see e.g. Olamaie et al. 2012, Giodini et al. 2013, and Köhlinger, Hoekstra & Eriksen 2015, De Martino & Atrio-Barandela 2016). This may result from either adopting inappropriate functional forms, often by extrapolating their use to cluster masses and/or redshifts that are not well sampled by observations or simulations, or by fitting models that depend on parameters to which the data are insensitive.

In this paper, we therefore move away from the traditional approach of assuming specific parameterised forms for the variation of cluster properties, such as the pressure, density and/or temperature distribution, and instead allow for a ‘free-form’ reconstruction, but one for which the level of complexity is determined automatically by the observational data and may depend on position within the cluster. This is achieved by representing each independent cluster property as some interpolating or approximating function (for example, a piecewise linear interpolation or a spline) that is specified by a set of control points, or ‘nodes’. The positions and amplitudes of these nodes and, most importantly, the number of nodes used are allowed to vary and constitute the set of parameters to be inferred from the data in a Bayesian manner.

We note that we have already successfully applied such a Bayesian nodal modelling approach to a number of cosmological analysis problems, including the reconstruction of the primordial anisotropy power spectrum and the variation of the dark energy equation-of-state parameter with redshift (see e.g. Vázquez et al. 2012a,b and Hee et al. 2016, 2017 ).

The structure of this paper is as follows. In Section 2, we briefly summarise Bayesian inference, in particular parameter estimation and model selection. Our nodal approach to modelling galaxy clusters is presented in Section 3, and in Section 4 we apply it to the particular case of modelling the pressure profile $P_{\rm e}(r)$ in a spherical cluster observed via its SZ effect. Section 5 outlines our Bayesian methodology for inferring the cluster parameters from interferometric SZ observations, and summarises our simulated cluster observations and real observations of the cluster MACS J0744+3927 using AMI. In Section 6, we present the results of our Bayesian nodal analysis of these simulations and real observations, and we conclude in Section 7.

2 Bayesian inference

For the analysis of some data $\mathbf{D}$ in the context of a model (or hypothesis) $\mathcal{M}$ that depends on some set of parameters $\mathbf{\Theta}$, Bayes’ theorem states that

[TABLE]

which is usually interpreted as the prior probability distribution $\Pr({\mathbf{\Theta}}|\mathcal{M})\equiv\pi(\mathbf{\Theta})$ of the parameters being updated by the likelihood $\Pr({\mathbf{D}}|{\mathbf{\Theta}},\mathcal{M})\equiv\mathcal{L}(\mathbf{\Theta})$ of obtaining the observed data given some set of parameter values to yield the posterior probability distribution $\Pr({\mathbf{\Theta}}|{\mathbf{D}},\mathcal{M})$ of the parameters, which is normalised by the Bayesian evidence $\Pr({\mathbf{D}}|\mathcal{M})\equiv\mathcal{Z}$ (which does not depend on the parameters $\mathbf{\Theta}$). The joint (unnormalised) posterior provides the complete inference in Bayesian parameter estimation, and can be subsequently marginalised over each parameter to obtain individual parameter constraints.

Similarly, in Bayesian model selection, one can calculate the probability of a model given the data as

[TABLE]

Taking the ratio of the probabilities of two models signifies our degree of belief in one model over another, and is given by the posterior odds ratio (POR)

[TABLE]

where $\mathcal{B}_{ij}\equiv\mathcal{Z}_{j}/\mathcal{Z}_{i}$ is the Bayes factor (Jeffreys, 1961). Clearly, if $\pi_{\mathcal{M}_{j}}=\pi_{\mathcal{M}_{i}}$, so that the two models are considered a priori equally probable, then $\mathcal{P}_{ij}=\mathcal{B}_{ij}$. Table 1 lists a modern version of Jeffreys’ criteria, which are used to give meaning to this quantification (Kass & Raftery, 1995; Feroz, 2013).

Thus, for Bayesian model selection, in contrast to parameter estimation, the evidence takes a central role. Typically, the evidence for each model is calculated separately and their ratios evaluated. From (1), the evidence is the normalisation constant for the posterior, which is given by

[TABLE]

where $D$ is the dimensionality of the parameter space. As the average of the likelihood over the prior, the evidence embodies the notion of Occam’s razor (see e.g, Jaynes 1986 and Sivia 2005): a simple theory with a compact parameter space will have a larger evidence than a more complicated one, unless the latter is significantly better at explaining the data.

The evaluation of the multidimensional integral (4) over the whole parameter space is a challenging numerical task. We perform this calculation here using the nested sampling algorithm Multinest (Feroz & Hobson, 2008; Feroz, Hobson & Bridges, 2009; Feroz et al., 2013). This Monte-Carlo method is targeted at the efficient calculation of the evidence, but as a by-product also produces posterior inferences for parameter estimation; it is also very efficient at exploring posteriors that contain multiple modes and/or large (curving) degeneracies.

We note that elsewhere (Hee et al., 2016, 2017), we have presented an alternative method for performing Bayesian model selection, without explicitly computing evidences, which uses a combined likelihood and introduces an integer model selection parameter $n$. If the total number of models under consideration is specified a priori, the full joint parameter space of the models is of fixed dimensionality and can be explored using standard sampling methods, without the need for trans-dimensional techniques, although the posterior is usually highly multimodal and hence nested sampling is again an obvious choice. Bayes factors, or more generally posterior odds ratios, may then be read off directly from the posterior on $n$, which is obtained by straightforward marginalization. To keep our discussion simple, however, we will not use this method here, but plan to apply it to nodal modelling of galaxy clusters in a forthcoming publication.

In closing this section, it should be mentioned that, in general, a gain in information via a Bayesian analysis may be achieved in several ways: it can occur because of a tightening of the parameter constraints, a shift in position of the peak(s) of the distribution from prior to posterior, or an increase in the evidence (see e.g. Trotta et al. 2008, Seehars et al. 2014 and Seehars et al. 2016).

3 Nodal model for a galaxy cluster

As described briefly earlier, in our nodal approach to free-form modelling of a galaxy cluster, each independent physical property of the cluster is represented by some interpolating or approximating function that is determined by a set of control points, or nodes. In principle, these functions can be fully three-dimensional to allow modelling of arbitrary structure in each property of interest in the cluster. To illustrate the method simply, however, we will consider here the special case of a spherical cluster, such that each property is a function only of radius $r$ from the cluster centre. Moreover, we will specialise still further to the case where one constructs a nodal model of just a single property of interest, described by some function $f(r)$. It is a straightforward matter to extend the following analysis to multiple properties of interest.

The basic idea is to represent $f(r)$ not with some standard parameterised functional form, as in most current cluster analyses, but in terms of a number $N$ of nodes in $(r,f)$-space, together with their corresponding positions $r_{n}$ and amplitudes $f_{n}$ $(n=0,1,2,\ldots,N-1)$ (thus each node can ‘move’ both horizontally and vertically– see figure 1; although we will in fact restrict the movement of the ‘end’ nodes, as described below). The $N$ nodes act as control points for the continuous function $f(r)$. In this way, one obtains a continuous free-form reconstruction of the profile $f(r)$, for which the complexity is regularised by the data under analysis.

One is free to choose from a wide range of possible interpolation or approximation methods, such as polynomials, rational functions, splines, Bézier curves or even Gaussian processes. It should be noted, however, that some forms of smooth interpolating or approximating function, such as cubic splines, have continuity requirements on the function and its derivatives which can significantly reduce the ability of the resulting $f(r)$ to reproduce abrupt features in the cluster profile (Vázquez et al., 2012). In principle, the nature of the interpolating or approximating function could be determined by performing a straightforward Bayesian model selection between the various options, although we will not consider that further here. Instead, for illustration, we choose here the simplest approach of linearly interpolating between the nodes; since we are performing an interpolation, rather than an approximation, one has $f_{n}=f(r_{n})$.

As mentioned above, given the nature of the one-dimensional function $f(r)$ that we wish to model in the case of a galaxy cluster, we typically restrict the movement of the first and last nodes (or end nodes) as follows. The first node has a fixed position, at the origin, such that $r_{0}=0$ always. Consequently, its corresponding amplitude parameter $f_{0}$ represents the central value $f(0)$, which is often of interest in galaxy cluster analyses. In contrast, the last node has a fixed amplitude of zero, such that $f_{N-1}=0$. Thus, its corresponding position parameter $r_{N-1}$ can sometimes be interpreted as an extent of the cluster $r_{\rm max}$, which is again often of interest (although this interpretation does depend on the nature of the quantity $f(r)$ being modelled). In our demonstration of the method presented in Section 4, $f(r)$ represents the electron pressure profile $P_{\rm e}(r)$ of the cluster and so $r_{N-1}$ represents the extent of the cluster gas. Figure 1 illustrates our linearly-interpolated nodal model of $P_{\rm e}(r)$ for $N=4$ nodes.

In the context of Bayesian inference, we denote the model consisting of $N$ nodes by $\mathcal{M}_{N}$, which has the parameters ${\mathbf{\Theta}}=\{r_{n},f_{n}\}$ for $n=0,1,2,\ldots,N-1$ (with delta function priors on $r_{0}$ and $P_{N-1}$), plus any further ‘global’ parameters, such as the cluster position on the sky or parameters describing any contaminating signals or noise in the data. The priors on the parameters may be chosen to accommodate whatever information is available a priori. In general, however, we typically also impose a ‘sorting condition’ on the node positions, such that $r_{n}<r_{n+1}$. This can be done straightforwardly, without the need to reject any samples. Indeed, if each node position is considered to be drawn from a uniform distribution in some given range, as is usually the case, the corresponding (non-separable) joint prior $\pi(r_{0},r_{1},r_{2},\ldots,r_{N-1})$ on the node positions, which incorporates the sorting condition $r_{n}<r_{n+1}$, has an analytic form in terms of the beta distribution (Handley, Hobson & Lasenby 2015). It is also worth mentioning that, although we will assume here that each model $\mathcal{M}_{N}$ is equally likely a priori, such that $\pi_{\mathcal{M}_{i}}=\pi_{\mathcal{M}_{j}}$ (and hence PORs are equivalent to Bayes factors), this is not necessary. One may view this assumption as imposing a uniform prior (within some range) on the number $N$ of nodes, but one could equally well impose, for example, a Poisson prior on $N$ with some given mean, from which one can read off the corresponding values of $\pi_{\mathcal{M}_{N}}$.

In the straightforward model selection approach used here, one begins by analysing the $N=2$ model, using Multinest to calculate the evidence $\mathcal{Z}_{2}$ and obtain samples from the posterior $\Pr({\mathbf{\Theta}}|{\mathbf{D}},\mathcal{M}_{2})$. This process is then repeated separately for $N=3,4,\ldots$ until the evidence $\mathcal{Z}_{N}$ has decreased well below that of the most favoured model. One may then proceed either by conditioning on the most favoured model $\mathcal{M}_{\hat{N}}$ and simply infer parameters and the corresponding constraints on $f(r)$ from the posterior $\Pr({\mathbf{\Theta}}|{\mathbf{D}},\mathcal{M}_{\hat{N}})$ or perform a ‘multi-model’ inference (see e.g. Parkinson & Liddle 2013, Hee et al. 2016, 2017, and Planck Collaboration et al. 2016) in which the constraints on $f(r)$ are determined by averaging over all the models $\mathcal{M}_{N}$ considered, weighted by their PORs. In the interests of simplicity we will adopt the former approach here.

4 Application to SZ observations

Although our nodal approach may be used to model observations of galaxy clusters in any waveband, we illustrate the method here by applying it to SZ observations (see e.g. Sunyaev & Zeldovich 1970, Birkinshaw 1999, and Carlstrom, Holder & Reese 2002). Since the SZ surface brightness is proportional to the line-of-sight integral of the electron pressure, as we show below, it is natural to use the nodal approach to model the radial electron pressure profile $P_{\rm e}(r)$ of the ionised gas within the cluster.

The observed SZ surface brightness $\delta I_{\nu}$ in the direction of an electron reservoir is given by

[TABLE]

Here $B_{\nu}$ is the blackbody spectrum, $T_{\rm CMB}=2.73$ K (Fixsen et al. 1996) is the temperature of the CMB radiation, $f(\nu)=\left(x\coth(x/2)-4\right)(1+\delta(x,T_{\rm e}))$ is the frequency dependence of thermal SZ signal, $x=\frac{h_{\rm p}\nu}{k_{\rm B}T_{\rm CMB}}$, $h_{\rm p}$ is Planck’s constant, $\nu$ is the frequency and $\rm{k_{\rm B}}$ is Boltzmann’s constant. The function $\delta(x,T_{\rm e})$ takes into account the relativistic corrections in the study of the thermal SZ effect which is due to the presence of thermal weakly relativistic electrons in the ICM and is derived by solving the Kompaneets equation up to the higher orders (Rephaeli 1995, Itoh, Kohyama, & Nozawa 1998, Nozawa, Itoh, & Kohyama 1998, Pointecouteau, Giard, & Barret 1998, Challinor & Lasenby 1998). It should be noted that at 15 GHz (AMI observing frequency) $x=0.3$ and therefore the relativistic correction, as shown by Rephaeli (1995), is negligible for $k_{\rm B}T_{\rm e}\leq 15\,\rm{keV}$. The dimensionless parameter $y$, known as the Comptonization parameter, is the integral of the number of collisions multiplied by the mean fractional energy change of photons per collision, along the line of sight

[TABLE]

where $P_{\rm e}(r)$ is the electron pressure at radius $r$ respectively, $\sigma_{\rm T}$ is Thomson scattering cross-section, $m_{\rm e}$ is the electron rest mass, $c$ is the speed of light and $dl$ is the line element along the line of sight. It should be noted that (6) assumes an ideal gas equation of state.

The integral ($Y_{\rm SZ}$) of the Comptonisation $y$ parameter over the solid angle $\Omega$ subtended by the cluster is proportional to the volume integral of the gas pressure. It is thus a good estimate for the total thermal energy content of the cluster and hence its mass (see e.g. Bartlett & Silk 1994). The $Y_{\rm SZ}$ parameter in cylindrical and spherical geometries, respectively, may be described as

[TABLE]

where $R$ is the projected radius of the cluster on the sky. Thus, by using our nodal approach to model $P_{\rm e}(r)$, one may constrain $Y_{\rm SZ}$ in either geometry.

5 Analysis of interferometric SZ data

In order to verify that our proposed model, with its corresponding assumptions, can describe profiles of cluster physical properties accurately, we carry out a Bayesian analysis of simulated SZ observations using AMI (AMI Consortium: Zwart et al., 2008) of a set of three clusters, as well as real AMI SZ observations of the cluster MACS J0744+3927 (Rumsey et al., 2016).

5.1 Bayesian methodology

An interferometer like AMI operating at a frequency $\nu$ measures samples from the complex visibility plane $\widetilde{I}_{\nu}({\bf u})$. These are given by a weighted Fourier transform of the surface brightness $I_{\nu}({\bf x})$, namely

[TABLE]

where ${\bf x}$ is the position on the sky relative to the phase centre, $A_{\nu}({\bf x})$ is the (power) primary beam of the antennas at observing frequency $\nu$ (normalised to unity at its peak) and ${\bf u}$ is the baseline vector in units of wavelength.

Details of our Bayesian methodology for modelling interferometric SZ data, primordial CMB anisotropies, and resolved and unresolved radio point-source are given in Hobson & Maisinger (2002); Feroz & Hobson (2008); Feroz et al. (2009); AMI Consortium: Davies et al. (2011) and AMI Consortium: Olamaie M. et al. (2012).

In short, the measured interferometer visibilities in our model are assumed to have the form

[TABLE]

where the signal part, $\widetilde{I}_{\nu}({\bf u})$, contains the contributions from the SZ cluster and identified radio point sources, whereas the generalised noise part, $N_{\nu}({\bf u})$, contains contributions from the background of unsubtracted radio point sources, primary CMB anisotropies and instrumental noise. The last two contributions to the generalised noise are well described by Gaussian processes, whereas the background of unsubtracted radio sources is strictly a Poisson process. Nonetheless, in the limit of a large number of faint unsubtracted sources, this contribution can also be well approximated as Gaussian (Feroz et al. 2009). Thus, the generalised noise is assumed to be Gaussian distributed, so that the likelihood function has the form

[TABLE]

where $\chi^{2}$ is the standard statistic quantifying the misfit between the observed data $D$ and the predicted data \mbox{\boldmathD}^{\rm p}(\mbox{\boldmath\Theta}):

[TABLE]

in which \mbox{\boldmathC}_{\nu,\nu^{\prime}} is the generalised noise matrix relating the frequency channels $\nu$ and $\nu^{\prime}$. Under the assumption that the three contributions to the generalized noise discussed above are independent, this matrix is simply the sum of the covariance matrices for each component. These matrices are described in Feroz et al. (2009) and are assumed to be independent of the parameters to be fitted in the analysis. In particular, the primary CMB anisotropies are assumed to be consistent with the concordance cosmology.

In the Bayesian analysis, the point sources are modelled using delta-function priors on position and Gaussian priors on flux and spectral index, usually determined from higher-resolution observations from the AMI Large Array, as discussed in Feroz et al. (2009). We focus here, however, on the inference of the cluster parameters. For the model $\mathcal{M}_{N}$, having $N$ nodes, the cluster parameters are

[TABLE]

where $x$ and $y$ are the cluster projected position on the sky. We carry out the analysis for $N=2,3,4,\ldots,8$. We assume that the priors on these sampling parameters are separable, apart from imposing the ‘sorting condition’ on the nodes positions, as discussed in Section 3, so that

[TABLE]

The position of the first node is fixed to $r_{0}=0$, whereas the position of the last node has a uniform prior in the range $0.5\leq r_{N-1}/\mbox{Mpc}\leq 2.5$. The positions of the intervening nodes have uniform priors in the range $0\leq r_{n}/\mbox{Mpc}\leq 1$. The prior on the pressure $P_{0}$ of the first node is a truncated exponential distribution with mean $\lambda=100$ MeVm*-3* in the range $10\leq P_{0}/\rm{MeVm}^{-3}\leq 5\times 10^{3}$, whereas the pressure $P_{N-1}$ of the last node is fixed to zero. The pressures $P_{n}$ of the intervening nodes have uniform priors in the range $0\leq P_{n}/\rm{MeVm}^{-3}\leq 500$. Finally, we assume Gaussian priors on the cluster position parameters $x$ and $y$, centred on the origin with a standard deviation of 1 arcmin. It is also worth mentioning that we assume here that each model $\mathcal{M}_{N}$ is equally likely a priori, such that $\pi_{\mathcal{M}_{i}}=\pi_{\mathcal{M}_{j}}$ (and hence PORs are equivalent to Bayes factors), which is equivalent to imposing a uniform prior on $N$ in the range $2\leq N\leq 8$. We note that the above priors are chosen to be consistent with the results of $N$-body simulations and real cluster observations of galaxy clusters (see e.g. Olamaie et al. 2012, Olamaie et al. 2013 and references therein). In all analyses, the redshift of the cluster is assumed known.

5.2 Simulated AMI observations

To generate simulated SZ skies and observe them with a model AMI instrument, we use the methods outlined in Hobson & Maisinger (2002), Grainge et al. (2002), Feroz et al. (2009), and AMI Consortium: Olamaie M. et al. (2012). In particular, we generate simulated observations of three different model clusters, which we call $\rm{SIM}_{1}$, $\rm{SIM}_{2}$ and $\rm{SIM}_{3}$. These differ in the input pressure profile used to generate them, and all include primordial CMB anisotropies, receiver thermal noise, and simulated residual points sources typical of AMI pointed observations of clusters.

For $\rm{SIM}_{1}$ and $\rm{SIM}_{2}$, we assume an input pressure profile that matches our nodal model precisely, in order to investigate the ability of our approach to recover the true parameter values in the simplest case. In particular, the pressure profiles for $\rm{SIM}_{1}$ and $\rm{SIM}_{2}$ are generated by linear interpolation between $N=3$ and $N=4$ nodes, respectively, and are plotted in figure 2. For both simulations, we assume the cluster lies at a redshift of $z=0.5$.

For SIM3, we use a more realistic cluster model to test the ability of our nodal approach to recover a pressure profile that is not of the form assumed in the analysis. In this case, the cluster is simulated using the model described in Olamaie et al. (2012, 2013), which assumes that the dark matter density follows an NFW profile and the ICM plasma pressure is described by the generalised NFW (GNFW) profile. The model also assumes that hydrostatic equilibrium is satisfied and that the local gas fraction is small throughout the cluster. This cluster model is fully specified by just three parameters, for which we assume the values $M_{\rm{tot}}(r_{\rm 200})=5\times 10^{14}\,\rm{M_{\odot}}$, $z=0.54$ and $f_{\rm{gas}}(r_{200})=0.13$. The resulting pressure profile is shown in figure 2, and is formally singular at the origin, decreasing sharply with radius in the central regions of the cluster.

5.3 AMI observations of MACS J0744+3927

We also analyse real AMI observations of MACS J0744+3927, one of the clusters in the CLASH (Cluster Lensing And Supernova survey with Hubble) sample (Postman et al., 2012). MACS J0744+3927 is a rich cluster at redshift z = 0.689 and has been studied through its X-ray emission, strong lensing, weak lensing and SZ effect (Schmidt & Allen, 2007; Ettori & Balestra, 2009; Umetsu et al., 2016; Rumsey et al., 2016). The SZ signal (decrement) on the AMI map appears circular, (figure 7 in Rumsey et al. 2016) in agreement with the X-ray surface brightness from the Chandra archive data (figure 6 in Postman et al. 2012). Details of AMI pointed observation towards the cluster, data reduction pipeline and mapping are described in Rumsey et al. (2016). In particular, we note that the Bayesian analysis includes 23 radio point sources in the AMI field. We focus here, however, on the determination of the parameters defining our nodal model of the cluster.

6 Results and Discussion

We now present the results of our Bayesian nodal analysis applied to the three simulated data sets $\rm{SIM}_{1}$, $\rm{SIM}_{2}$, and $\rm{SIM}_{3}$, and real AMI observations of the cluster MACS J0744+3927.

6.1 Simulation SIM1

The results of our evidence-based Bayesian model selection analysis to determine the number $N$ of nodes in the $P_{\rm e}(r)$ reconstruction for $\rm{SIM}_{1}$ are given in figure 3. This shows the histogram of Bayes factors $\mathcal{B}_{2j}$, i.e. relative to the $N=2$ ‘straight line’ model, as a function of the number of nodes $N$ in the reconstruction of the pressure profile $P_{\rm e}(r)$ for simulation SIM1. Given our prior choice $\pi_{\mathcal{M}_{i}}=\pi_{\mathcal{M}_{j}}$ on the models, the Bayes factors are equal to the PORs $\mathcal{P}_{2j}$, and so the plotted histogram is equivalent to the marginalised posterior on $N$. Recalling that the input pressure profile for SIM1 is constructed by linear interpolation between $N=3$ nodes, one sees that our analysis has recovered the true value of $N$ as the most favoured. In particular, it is worth noting that the log Bayes factor $\ln\mathcal{B}_{23}=4.0\pm 0.12$, indicating strong evidence for the $N=3$ model over the $N=2$ (straight line) model, according to Table 1. The Bayes factors then gradually decline for $N>3$, ultimately reaching the value $\log\mathcal{B}_{28}=0.0\pm 0.12$, which indicates no preference for $N=8$ over the $N=2$ model, demonstrating that the ability of the $N=8$ model to (over)fit the data is offset by the penalty of its increased complexity.

As mentioned earlier, we will adopt here the straightforward approach of determining the constraints on parameters by conditioning on the most favoured model $\hat{N}=3$, rather than performing model averaging according to their Bayes factors. Since it is of interest to understand any biases or constraints on the parameters imposed by our choice of priors, we first consider the ‘posterior’ $\Pr({\mathbf{\Theta}}_{\rm c}|{\mathbf{D}}={\mathbf{0}},\mathcal{M}_{3})$ on the cluster parameters obtained in the absence of any data. This is calculated simply by setting the likelihood to a constant value, so that the sampler explores just the prior $\pi({\mathbf{\Theta}}_{\rm c})$. The resulting 1-D and 2-D marginalised distributions for the sampling parameters ${\mathbf{\Theta}}_{\rm c}$ are shown in figure 4 (left panel); in addition we plot the derived parameter $Y_{\rm tot}$ defined in equation (7). These plots show that we correctly recover the assumed prior distributions, and also reveal the constraints that our choice of priors has placed on the derived parameter $Y_{\rm tot}$. The plots are produced using the open source Python library corner.py (Foreman–Mackey, 2016).

The corresponding plot obtained after analysing the simulation SIM1 is shown in figure 4 (right panel), and shows the effect that the data have in updating the prior, via the likelihood function, to produce the posterior, in the spirit of Bayes’ theorem. In particular, one sees that the cluster position ($x$ and $y$) on the sky is firmly constrained and the true values lie within few arcseconds of the means of the posterior probability distributions for $x$ and $y$. Turning to the parameters defining the nodal model of the pressure profile, one sees that the amplitude $P_{0}$ of the first node is not well constrained, and we are essentially recovering just the prior distribution. This is to be expected, since an interferometric SZ observation is insensitive to lengths scales on the sky that correspond to Fourier modes in the visibility plane lying well below the shortest baseline (in units of wavelengths) of the interferometer. Consequently, for this simulation, the observations cannot probe the cluster inner core and thus provide no information on the pressure $P_{0}$ at the centre. By contrast, the position and amplitude of the second node, ($r_{1}$, $P_{1}$), are both constrained relative to their prior distributions, although their 2D marginal reveals a clear degeneracy between them. The position $r_{2}$ of the third (and final) node is very well constrained. Since this node corresponds to the point at which the gas pressure drops to zero, it is a valuable quantity for defining the extent of the gas distribution in the cluster. Indeed, obtaining a robust constraint on this quantity can provide insight to the dynamical state of the cluster. Finally, we note that the important derived parameter $Y_{\rm tot}$ is also very well constrained. From their 2-D marginal, however, one sees that there is some degeneracy between the parameters $r_{2}$ and $Y_{\rm tot}$.

Rather than viewing the posterior constraints on the individual parameters $r_{1}$, $r_{2}$, $P_{0}$ and $P_{1}$, as in figure 4, it can be more intuitive and instructive to plot the corresponding inference on the reconstructed pressure profile $P_{\rm e}(r)$ directly. This may be performed in a number of ways. For example, one may plot the posterior probability $\Pr(P|r,{\mathbf{D}},\mathcal{M}_{3})$, in normalised slices at constant $r$ (Hee et al., 2016). We will not pursue that method here, however, and instead adopt the simpler approach of plotting the mean position and amplitude of each node, together with horizontal and vertical error bars representing the 68 percent marginalised Bayesian credible interval in each direction; these nodes are then linearly interpolated to obtain the reconstructed pressure profile, as shown in figure 5. The input pressure profile used to generate the simulation is also plotted. As one might expect from figure 4 (right panel), the reconstructed pressure profile is consistent with the input profile, with the true node locations in $(r,P)$-space all lying within the $68\%$ Bayesian credible intervals of the corresponding inferred node locations. It is again clear, however, that the central pressure $P_{0}$ is poorly constrained, as discussed above, but that the remaining node parameters are reasonably well determined.

In the remainder of the paper, we will plot the reconstructed $P_{\rm e}(r)$ directly, as in figure 5, together with the posterior distributions on the position ($x$ and $y$) of the cluster and the important physical parameters $r_{N-1}\equiv r_{\rm max}$ and $Y_{\rm tot}$.

6.2 Simulation SIM2

The histogram of Bayes factors $\mathcal{B}_{2j}$ as a function of $N$ obtained in the analysis of simulation SIM2 is shown in figure 6. Recalling that the input pressure profile for SIM2 is constructed by linear interpolation between $N=4$ nodes, one sees that our analysis has again recovered the true value of $N$ as the most favoured. In this case, one sees the $N=2$ (straight-line) model is again strongly disfavoured. In particular, one finds $\ln\mathcal{B}_{24}=4.0\pm 0.12$, indicating strong evidence for the most favoured $N=4$ model over the $N=2$ model, according to Table 1. The Bayes factors then gradually decline for $N>4$ in a similar way to that found for SIM1, ultimately reaching the value $\log\mathcal{B}_{28}=0.0\pm 0.12$; this again indicates no preference for $N=8$ over the $N=2$ model, for the reasons discussed above.

Conditioning on $N=4$ nodes, figure 7 shows the reconstructed pressure profile $P_{\rm e}(r)$ (top panel), and the 1-D and 2-D marginal posterior distributions on the cluster position parameters $x$ and $y$ (middle panel) and on the gas extent $r_{3}$ and total Comptonisation parameter $Y_{\rm tot}$ (bottom panel). One again sees that the reconstructed pressure profile is consistent with the input profile used to generate the simulation, but that cluster gas pressure $P_{0}$ is poorly constrained, for the reasons we discussed above in the context of SIM1. The remaining node parameters are again reasonably well constrained, especially the first internal node parameters $(r_{1},P_{1})$. Moreover, one again obtains tight constraints on the cluster position, consistent with the input values. The cluster extent $r_{3}$ and total Comptonisation parameter are also both well constrained and in agreement with the input values. As in SIM1, however, the 2-D marginal distribution of $r_{2}$ and $Y_{\rm tot}$ reveals a mild degeneracy.

6.3 Simulation SIM3

The histogram of Bayes factors $\mathcal{B}_{2j}$ as a function of $N$ obtained in the analysis of simulation SIM3 is shown in figure 8. The input pressure profile for SIM2 is not constructed from a linear interpolation between nodes, but instead from the cluster model of Olamaie et al. (2012, 2013), which assumes that the pressure is described by the generalised NFW (GNFW) profile (Nagai, Kravtsov & Vikhlinin, 2007). Hence, in this case, there is no ‘correct’ number of nodes to recover. Instead, the most favoured value of $\hat{N}=3$ nodes gives an indication of the level of complexity in the pressure profile reconstruction that is supported by the data. Thus, for interferometric SZ observations of the type simulated, the data support only a very simple reconstruction of the pressure profile, favouring a representation consisting of just two straight-line segments. Indeed, the variation of the Bayes factor with $N$ is very similar to that shown in figure 3 for SIM1, for which the input pressure profile had precisely this simple form.

Conditioning on $N=3$ nodes, figure 9 shows the reconstructed pressure profile $P_{\rm e}(r)$ (top panel), and the 1-D and 2-D marginal posterior distributions on the cluster position parameters $x$ and $y$ (middle), and the gas extent $r_{3}$ and total Comptonisation parameter $Y_{\rm tot}$ (bottom). Although in this simulation there are no ‘correct’ locations for the nodes in $(r,P)$-space, one sees that the reconstructed node locations yield reconstructed pressure profiles that are consistent with the input one, although the plotted ‘best-fit’ profile is somewhat shallower than the GNFW profile assumed in the simulation. One can understand this behaviour by recalling that the SZ Comptonisation parameter given in (7) is proportional to the integral of $P_{\rm e}(r)$ along the line-of-sight through the cluster. Thus, the larger amplitude of the first (only) internal node relative to the input profile can offset the smaller amplitude of the node at the origin. Once again our analysis serves to highlight the rather coarse level of detail in the reconstructed pressure profile that is achievable with SZ observations of the type simulated. Nonetheless, one still obtains tight constraints on the cluster position, which are consistent with the input values, and also on the cluster extent $r_{3}$ and total Comptonisation parameter. We again see, however, that there is a slight degeneracy in the 2-D marginal distribution of $r_{2}$ and $Y_{\rm tot}$.

6.4 Cluster MACS J0744+3927

The histogram of Bayes factors $\mathcal{B}_{2j}$ as a function of $N$ obtained in the analysis of real AMI observations of the cluster MACS J0744+3927 is shown in figure 10. For this real cluster, one sees that the variation of the Bayes factor $\mathcal{B}_{2j}$ is broadly similar to that shown in figure 8, obtained from the analysis of simulation SIM3, for which an input GNFW pressure profile was assumed. In particular, the most favoured model again has $N=3$ nodes, after which the Bayes factors gradually decline with increasing $N$. This indicates that the data support a model for the pressure profile that is no more complex than two straight-line segments. It is worth noting, however, that in comparing the favoured $N=3$ model with the base $N=2$ (straight-line) model, one obtains $\ln\mathcal{B}_{23}=6.0\pm 0.2$, which corresponds to a decisive favouring of the former model, according to Table 1. Thus, one may deduce at high confidence from our analysis, in a model independent manner, that the pressure profile is not simply a linear function of $r$.

Conditioned on $N=3$ nodes, figure 11 shows the reconstructed pressure profile $P_{\rm e}(r)$ (top panel), and the 1-D and 2-D marginal posterior distributions on the cluster position parameters $x$ and $y$ (middle), and the gas extent $r_{2}$ and total Comptonisation parameter $Y_{\rm tot}$ (bottom). Once again the central pressure $P_{0}$ is poorly constrained, but the remaining node locations in $(r,P)$-space are better determined. It is interesting that the resulting reconstructed pressure profile is convex, rather than concave such as those obtained from the simulated observations, although the uncertainty in the location of the single internal node also allows for concave profiles. It is also worth recalling from our analysis of simulated data how SZ observations of this type allow for only a very coarse reconstruction of the cluster pressure profile, owing to the SZ effect being proportional to the line-of-sight integral of the pressure. Indeed, from figure 9, we recall that in the analysis of simulation SIM3 the reconstructed pressure profile was somewhat flatter than the input GNFW profile, although still consistent with it to within the error-bars, and a similar effect could be occurring in figure 11.

Nonetheless, as we found in our analysis of simulated data, the cluster position ($x$ and $y$) is well constrained, and the corresponding 2-D marginal posterior distribution shows no sign of degeneracy. The gas extent $r_{3}$ and total Comptonisation parameter $Y_{\rm tot}$ are also both very well determined, although their 2-D marginal posterior distribution shows that same slight degeneracy as seen in the analyses of the simulation observations.

7 Conclusions

Almost all current approaches to modelling observations of galaxy clusters rely on assuming some parameterised functional form for the properties of the cluster, such as gas density, dark matter density or temperature. A generic weakness of this approach is that these functional forms have usually been arrived at through empirical means, via the analysis of $N$-body simulations or observations, and are often chosen to have simple analytic expressions, rather than being fundamental or physically well motivated.

In this paper, we have moved away from this approach and presented a free-form model for the physical properties of galaxy clusters. Previous attempts to model clusters in this way have typically relied simply on dividing the cluster into a predefined number of cells, or concentric shells for spherical clusters, and determining the value of each physical quantity of interest within these subregions (Tchernin et al., 2015). Such approaches typically lead to under-determined inverse problems that therefore need to be regularised in some way. There is considerable freedom in how to choose the level or nature of the regularisation to apply, and the results can vary significantly depending on how this choice is made. We have therefore presented an alternative approach to free-form reconstruction in which the complexity of the model is determined directly from the data. This is achieved by representing each independent cluster property as some interpolating or approximating function that is specified by a set of control points, or ‘nodes’, for which the number of nodes, together with their positions and amplitudes, are allowed to vary and are inferred from the data in a Bayesian manner, employing both model selection and parameter estimation.

To demonstrate our approach in a simple setting, we have applied it to the particular case of modelling interferometric SZ observations of spherical galaxy clusters. In this context, the free-form part of the cluster model is simply a nodal representation of the electron pressure profile $P_{\rm e}(r)$. We have performed Bayesian analyses of simulated observations with the Arcminute Microkelvin Imager (AMI) of three separate model clusters.

In the first two simulations, the input pressure profile has the same form as that assumed in the analysis, namely a linear interpolation between a set of $N$ nodes (with $N=3$ and $N=4$, respectively). We showed that, in both cases, our Bayesian model selection analysis returned the true value of $N$ as the most favoured. Moreover the resulting reconstructed pressure profiles were consistent with those used as input. In our third simulation, in which the input pressure was assumed to follow a GNFW, the most favoured model again had $N=3$ nodes, and the resulting reconstructed pressure profile was consistent with the input one. In all cases, we found that the central pressure of the cluster is not well determined, since interferometric observations of the type simulated do not probe length scales corresponding to the inner core. In the analysis of our third simulation, we also noted that the reconstructed pressure profile was somewhat shallower that the singular GNFW profile used to generate the simulation (although still consistent with it), which results from the SZ effect being proportional only to the line-of-sight integral of the pressure in the cluster. A general feature of our results is that SZ interferometric observations of this type allow for only a very coarse reconstruction of the cluster pressure profile. Nonetheless, we also find that in all cases one obtains tight constraints on the cluster position, and that the cluster extent and total Comptonisation parameter are also both well determined.

We also applied our approach to real AMI observations of the cluster MACS J0744+3927. We found that the most favoured model has $N=3$ nodes, and that the corresponding best-fit reconstructed pressure profile is convex, although the uncertainties of the node locations in $(r,P)$-space also allow for concave profiles. As we found in the analysis of simulations, the central pressure is poorly determined but the remaining node parameters are reasonably well constrained. Once again, we found that cluster position, cluster gas extent and total Comptonisation parameter are all very well constrained.

In closing, some further general points and avenues for future research are worth discussing. First, the tight constraints obtained on the cluster position and on the two very important cluster parameters $r_{\rm max}$ and $Y_{\rm tot}$ demonstrate the robustness of our approach. Moreover, with only minor modification, the method may prove very useful in cluster detection. Although, for the sake of illustration, we assumed the cluster redshift in our analyses presented here, this is not necessary. One can easily re-perform the analysis by instead constructing a nodal model for the pressure profile $P_{\rm e}(\theta)$, where $\theta$ is the projected angle on the sky from the centre of the cluster. In this way, the approach does not depend on the redshift, but will still produce the tight constraints on the cluster position, angular extent and $Y_{\rm tot}$. Finally, since our Bayesian approach to the inference produces posterior weighted samples in the parameter space, further directions for future development include defining other derived parameters that capture particular features of interest in the pressure profile. One example, motivated by our analysis of the cluster MACS J0744+3927, would be a statistic that embodies the concavity or convexity of the pressure profile. Others might include a parameter that quantifies the cuspy versus core nature of the central region of the cluster. In any case, one may easily use the posterior samples to determine the full (joint) posterior distribution of such derived parameters.

Acknowledgments

The authors thank William Handley for illuminating discussions regarding Bayesian inference. This work was performed using both the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/), and COSMOS Shared Memory system at DAMTP, University of Cambridge operated on behalf of the STFC DiRAC HPC Facility. Darwin Supercomputer is 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. COSMOS Shared Memory system is funded by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1, ST/K00333X/1. We are grateful to Stuart Rankin and COSMOS management team for their computing assistance. MO thanks Astro Hack Week 2016 for valuable discussions and insights on Bayesian inference and statistics. YCP acknowledges support from a Trinity College Junior Research Fellowship.