Free Access
Volume 639, July 2020
Article Number A73
Number of page(s) 10
Section Numerical methods and codes
Published online 09 July 2020

© ESO 2020

1. Introduction

Galaxy clusters trace the backbone of the Universe, and their thermodynamic properties are valuable assets, for instance, for probing the physical structure and substructure of the intracluster medium (ICM) to unveil the occurrence of phenomena such as merger events (Mroczkowski et al. 2019). For example, the cluster thermal history is fully captured by the entropy, and large-scale deviations from a power-law behaviour can be used to examine how the ICM is affected by non-gravitational processes, such as heating and radiative cooling of active galactic nuclei (AGN; Voit 2005). The cooling-time profile may be used to distinguish between clusters with and without a cool core (Hudson et al. 2010), while the gas fraction directly relates to the strength of radiative cooling and star formation (Sun et al. 2009). Sharp jumps in the temperature or pressure profiles generally indicate shocks and cold fronts (Markevitch & Vikhlinin 2007). Furthermore, thermodynamic profiles allow us to infer the mass distribution inside the cluster, and from the latter, to improve cosmological constraints (Bocquet et al. 2015; Ruppin et al. 2019), such as the dark energy equation of state (w), the number of neutrino species, the matter density (ΩM), and the amplitude of matter power spectral fluctuations on 8 Mpc h−1 scales (σ8).

Thermodynamic profiles of galaxy clusters can be derived from observations in the optical band, the X-ray band, or in microwaves in the shape of the Sunyaev-Zeldovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972). In particular, SZ measurements became widespread in the past decade (e.g. Birkinshaw & Lancaster 2005; Mroczkowski et al. 2009; Korngut et al. 2011; Sayers et al. 2013; Adam et al. 2015; Romero et al. 2017) because of the huge increment in high-resolution SZ instruments (see Mroczkowski et al. 2019; Castagna & Andreon 2019, for details).

Combining SZ and X-ray observations to gather thermodynamic profile is extremely advantageous because the two wavelengths encode information about the ICM (Ameglio et al. 2007). The highest benefit of this approach applies to high redshift and to the outskirts of the cluster, where the very low surface brightness of the X-ray signal and its low contrast against the background leads to noisy temperature estimates that are often affected by systematics related to the difficult X-ray background spectral modelling (e.g. Eckert et al. 2011). On the other hand, SZ measurements do not experience dimming because they are nearly independent of redshift. However, the joint derivation method is not straightforward, mostly because of cross-correlation among thermodynamic measures: in SZ, the conversion factor from Compton y to surface brightness is temperature dependent (a change within 5% between T = 5 keV and T = 10 keV at 150 GHz, within 15% at 260 GHz), and in X-ray, the conversion from electron density to a soft-band count rate depends on both temperature (by 2% for the same temperature change) and metallicity (20% change from 0.3 solar to 1 solar; Ettori et al. 2013). Furthermore, wavelength-specific temperatures may differ from one another: the SZ temperature is gas-mass weighted, while the X-ray temperature is derived spectroscopically.

Several authors have paved the way for joint X-SZ analyses, mostly considering electron density and SZ data (e.g. Kitayama et al. 2004; Ameglio et al. 2007; Mroczkowski et al. 2009; Eckert et al. 2013; Adam et al. 2015; Shitanishi et al. 2018; Siegel et al. 2018; Ghirardini et al. 2018; Ruppin et al. 2020), where electron density and SZ strength were often derived for a given temperature and metallicity. In a different approach, we replace the electron density fit with a full spatial-spectral X-ray data fit. This provides consistent and improved estimates of all thermodynamic profiles, and the temperature and metallicity dependencies are included in the error budget.

This paper presents JoXSZ, the first publicly available program for performing a multiwavelength joint fit of SZ and X-ray data on galaxy clusters. JoXSZ supports data coming from different X-ray or SZ telescopes, follows a fully Bayesian forward-modelling approach, and supports flexible modelling of the thermodynamic components of a cluster. Users are free to choose which parameter to fit, to decide whether hydrostatic equilibrium should be adopted, whether there are systematics between X-ray and SZ temperatures, and other options. JoXSZ returns the maximised likelihood, model estimates with uncertainties, joint and marginal probability contours, convergence diagnostics, and projected radial profiles, as detailed below. When data are missing in the input file, for instance, measurements are lacking at some radii, JoXSZ is able to deal with them automatically by excluding values marked as missing (nan) from the likelihood computation.

The paper is organised as follows: in Sect. 2 we provide an overview of the software and the technical requirements; in Sect. 3 we describe in detail the method behind each step of the program; in Sect. 4 we present an application of JoXSZ to real data from the galaxy cluster CL J1226.9+3332; and we conclude with the discussion and final remarks in Sect. 5. Appendix A presents the technical implementation of JoXSZ.

2. JoXSZ

2.1. Program flow

As represented in Fig. 1, JoXSZ adopts a descriptive model for a spherically symmetric galaxy cluster with a given centre that parametrises its pressure and electron density profile, and derives the temperature profile as the ratio of these quantities through the ideal gas law. The X-ray and SZ-based temperatures can be asked to be consistent or to differ, for example to study the cluster elongation along the line of sight, gas clumping, or calibration uncertainties. The model assumes a flat metallicity profile for the cluster, whose value can be fitted, or fixed, at user request. Other thermodynamic quantities (entropy, cooling time, and gas mass) are automatically computed, output, and plotted by JoXSZ. JoXSZ fits the surface brightness profile in the SZ domain and across multiple X-ray bands. Under the HE assumption, JoXSZ also computes the mass profile and the gas fraction profile, and derives rΔ and MΔ. The SZ and X-ray modelling structures rely on the PreProFit (Castagna & Andreon 2019) and MBProj2 (Sanders et al. 2018) pipelines, respectively. JoXSZ merges these two processes into a unique joint and consistent model based on a Markov chain Monte Carlo (MCMC) fitting algorithm.

thumbnail Fig. 1.

Block diagram showing the program flow. Radial profiles are plotted in red. Options are shown in green. Analysis pipelines are depicted in yellow. Data enter in the blue box.

Open with DEXTER

2.2. Requirements and installation

JoXSZ was developed and tested with Python 3.6. The following libraries are required: mbproj2, PyAbel, numpy, scipy, astropy, emcee, six, matplotlib, and corner. JoXSZ can be downloaded from GitHub1.

3. Methods

3.1. Intracluster medium modelling

3.1.1. Pressure profile

The pressure profile is described by the generalised Navarro, Frenk & White (gNFW) model proposed by Nagai et al. (2007),


where P0 is a normalising constant and rp is a scale radius. The exponentials b and c describe the logarithmic slopes at r/rp ≫ 1 and r/rp ≪ 1, respectively, while a governs the turnover rate between these two slopes. The five parameters mean that the model is very flexible in fitting current data.

The three-dimensional pressure model is numerically integrated along the line of sight in order to obtain the two-dimensional map of the Compton y parameter, built on a regular grid assuming radial symmetry. As discussed in detail in Castagna & Andreon (2019), the integration is conducted via Abel transform and the upper integration limit Rb is at user choice.

3.1.2. Density profile

To parametrise the electronic density profile, we adopted the model introduced by Vikhlinin et al. (2006),


where ne0 is a normalising constant, α is the logarithmic slope at r/rc ≪ 1, rc represents the core radius, β is the shape parameter for the isothermal β-model (Cavaliere & Fusco-Femiano 1976), rs is the radius at which the density profile steepens with respect to the traditional β-model, ϵ interprets the change in slope near rs, and γ measures the width of the transition region. The second β component, that is, the additive term in Eq. (2), increases model flexibility close to the centre of the cluster by including a small core radius rc2 with shape parameter β2 and additive constant ne02.

3.1.3. Temperature profile

Assuming the ideal gas law, we derive the temperature profile as the ratio between the pressure and density profiles,


where kB is the Boltzmann constant. This temperature calculation is simultaneously influenced by the pressure constraints from SZ observations and the electron density constraints from X-ray observations. Temperature gradients manifest themselves as gradients in the ratio of X-ray surface brightness profiles in different bands. As already mentioned and as detailed in the diagram in Fig. 1, JoXSZ allows users either to consider a unique temperature profile TSZ = TX, or to make a distinction between the gas-mass weighted temperature TSZ and the X-ray temperature TX, introducing the multiplicative parameter log(TX/TSZ). Thus, JoXSZ assumes identical shapes for the two temperature profiles, only allowing different normalisations. The profiles of the entropy K(r), cooling time tcool(r), and gas mass Mgas(< r) are computed as usual (Sanders et al. 2018).

3.1.4. Mass distribution

If requested, the total mass distribution can be derived assuming hydrostatic equilibrium (HE). The mass M(< r) enclosed within a radius r is related to the pressure and electronic density profiles through


where μgas  =  0.61 is the mean molecular gas mass (Anders & Grevesse 1989), mp is the proton mass, and G is Newton’s constant. Because the total mass distribution as defined in Eq. (4) directly depends on the pressure and electron density parameters, a non-monotonically increasing mass profile may occur, and this means that the program would assign negative values of mass at some radii. To avoid this, a positive prior on the radial derivative of the enclosed mass is imposed by default, although users can remove it.

The adopted version of the hydrostatic equilibrium formula exploits the advantages of SZ and X-ray surface brightness profiles over X-ray temperature estimates, especially at large distances from the cluster centre. JoXSZ uses the analytic expression for the pressure derivative to improve program execution time.

By comparing the mass distribution in Eq. (4) with its definition in terms of volume and density,


the overdensity radius rΔ is derived as the radius within which the average density is Δ times the critical density at the cluster redshift, ρc(z). MΔ computation ensues as MΔ = M(rΔ).

The gas fraction profile is straightforwardly deduced as


where Mgas(< r) is obtained by integrating the product between the gas density profile and the volume of the shell.

3.1.5. Conversions from physical to instrumental quantities

In X-ray, the conversion from emissivity to count-rate depends on the gas temperature TX, metallicity, and Galactic absorption. As in Sanders et al. (2018), JoXSZ uses XSPEC (Arnaud 1996) to compute this conversion, which requires users to have computed the response matrix file (RMF) and the ancillary response file (ARF).

In SZ, the conversion from Compton y to surface brightness depends on the temperature TSZ because of the relativistic corrections to the SZ effect (Itoh & Nozawa 2004). The exact value depends on the specific response of the instrument, and thus users are required to provide the instrument-specific conversion factor as an input file. JoXSZ also accounts for the SZ calibration uncertainty, implemented as a Gaussian whose parameters can be set by the user. Both SZ and X-ray conversion factors are depend on the radius and are updated at each step of the chain.

3.2. Processing backbones: PreProFit and MBProj2

JoXSZ relies on PreProFit (Castagna & Andreon 2019) and MBProj2 (Sanders et al. 2018), where full details on the SZ and X-ray processing are presented. In a nutshell, PreProFit takes as input the pressure profile, whose sampling step is at user choice, then projects it onto a two-dimensional map using forward Abel transform, convolves the map with the instrumental beam and the transfer function, and finally derives the surface brightness profile through opportune conversion factors.

About the X-ray analysis, MBProj2 makes full use of the spatial-spectral X-ray data cube and allows users to model one or more thermodynamic profiles, but not the pressure profile. JoXSZ instead parametrises it together with the electron density profile and derives the remaining ones from the ideal gas law (see Appendix A for the technical implementation in Python). As mentioned, with the further assumption of hydrostatic equilibrium, JoXSZ computes the mass profile and the characteristic radius rΔ. The X-ray fitting procedure includes a component that accounts for the background level systematic. As anticipated in Sect. 3.1.3, the user is free to choose whether to include a single temperature profile that is simultaneously fitted from SZ and X-ray data or to consider two separate temperature profiles.

3.3. Model definition

To set up our joint model, we defined the likelihood function as the product of the two distinct functions for the SZ and X-ray analyses: the former is a χ2 function, the latter is a Poisson likelihood. As a result, the log-likelihood function ln(ℒJoXSZ), which is computationally more convenient, is determined as the sum of the wavelength-specific components ln(ℒSZ) and ln(ℒX),


The likelihood function for SZ data is the same as in Castagna & Andreon (2019), where fdataSZ and fmodelSZ are the observed and estimated surface brightness values, respectively, while σdataSZ is the error measure, and nSZ represents the total number of available data points. The likelihood function for the nX X-ray data (Sanders et al. 2018) compares the predicted model values M (including background) with the observed profiles D, assuming that the X-ray counts follow a Poisson distribution. Γ is the usual gamma function.

3.4. Prior, posterior sampling, and diagnostics

As presented in Sect. 3.1, JoXSZ allows extremely flexible modelling of the ICM, which involves a large number of parameters. Users have to select their prior distributions, and then the program estimates the posterior distribution thorugh an MCMC, marginalising over all of them according to the Bayesian approach. For the whole MCMC setup, we relied on the implementation refined by Sanders for MBProj2. The posterior is sampled with an affine-invariant ensemble sampler as proposed by Goodman & Weare (2010) and implemented in emcee (Foreman-Mackey et al. 2013). The user has to specify the list of parameters to be fitted and optionally can change the prior distribution, whether uniform or Gaussian. The user is free to fix the desired number of random walkers, the number of iterations, and the burn-in period extent. The MCMC is initialised with default parameter values that can be also changed upon user request, and is automatically followed by a preliminary optimisation towards higher values of likelihood to facilitate convergence. Multi-threading computation is supported by JoXSZ and is strongly encouraged to minimise the execution time.

Qualitative and quantitative diagnostics (described in detail in Sect. 4) are provided by JoXSZ to evaluate the convergence of the chains to the stationary distribution: traceplot and the cornerplot are automatically displayed, informing the user of the parameter evolution across iterations and of the joint posterior distribution; the acceptance fraction is reported in the program output; and plots with best-fit profiles and uncertainty intervals (at a probability level set at user choice) overplotted on the observed data are automatically produced.

Because the user can choose to fit different combinations of parameters, model selection plays an essential role in performing optimal analyses. When the models to be compared are nested (e.g. a fit with log(TX/TSZ) = 0 compared to an analysis with free ratio), the Savage-Dickey density ratio is recommended (e.g. Trotta 2007). This index represents an approximation of the Bayes factor and can be computed as the ratio between the marginal posterior of the more complex model evaluated at the simpler model parameter value and the prior density of the more complex model evaluated at the same point.

3.5. Execution-time balance

The fit performed with JoXSZ can be slow, largely because of highly time-consuming operations included in the SZ processing, as was discussed in Castagna & Andreon (2019). The precise splitting of the CPU time across the SZ and X-ray parts of the analyses depends on the relative size of the fitted data, on the relative precision adopted for the Abel transform of the SZ and X-ray sides, and on whether the PSF convolution is performed on both sides or only on one. In our example below, we used a widely unbalanced setting: the X-ray side does not use any convolution by the point spread function (PSF), while the SZ convolution adopts an unnecessary small pixel to compute it and the Abel transform adopts again an excessively small pixel in the SZ side. In such unbalanced conditions, the SZ fit accounted for 90% of the CPU time, while the remaining 10% was used for the X-ray fit. The fit can be made faster and less unbalanced between the two sides (e.g. 75% vs. 25%) without loosing accuracy by doubling the integration or convolution pixel size on the SZ side. As always, it remains at user charge to determine the trade-off between execution time and the precision that satisfies their needs and requirements.

3.6. Limitations

To mention the main limitations of the program, JoXSZ analyses one object at a time (it cannot deal with superposed clusters or point sources; flagging should be used instead), assumes a centre for the cluster, spherical symmetry (information on the third dimension is mandatory for integrating along the line of sight), radius-independent metallicity (because joint X-SZ data sets with robustly measurable metallicity gradients are rare), and negligible covariance between errors of the different radial bins.

4. A worked example

To highlight the capabilities of JoXSZ, we present an application of the program to real data. Our analysis purely illustrates the use of the code, it is not meant to be the most accurate astrophysical analysis of the cluster.

We analysed the high-redshift cluster of galaxies CL J1226.9+3332 (z = 0.89), which has been variously studied by several authors in the past years in the X-ray (Maughan et al. 2004, 2007; Donahue et al. 2014) and in SZ (Mroczkowski et al. 2009; Korngut et al. 2011; Sayers et al. 2013; Adam et al. 2015; Romero et al. 2017, 2018). CL J1226.9+3332 is a hot and massive cluster that was discovered in the WARPS survey (Ebeling et al. 2001). The cluster presents a relaxed morphology on a large scale, with some possible evidence of a disturbed core (Maughan et al. 2007; Korngut et al. 2011; Adam et al. 2015; Romero et al. 2017). We flagged the point source detected by Adam et al. (2015) at (RA, Dec) = (12:26:59.855, +33:32:35.21) with an aperture of 46″.

4.1. NIKA and Chandra data

The SZ observation, instrumental beam, transfer function, and the table of temperature-dependent Compton y to mJy conversion factors used in the example come from the publicly available NIKA data release2 (Catalano et al. 2014). We refer to Adam et al. (2015) for details about the cluster observation, whose reduction and filtering differ slightly, however.

We used Chandra X-ray observations (OBSID = 3180, PI. Ebeling, exposure time of 32 ks), reduced following the standard procedures using CIAO 4.1 and CALDB 4.5.2 (e.g. Andreon et al. 2019). First, data were flare filtered. Then, point sources were detected by a wavelet detection algorithm and masked. Events in ten bands were extracted: [0.7−1], [1−1.3], [1.3−1.6], [1.6−2], [2−2.7], [2.7−3.4], [3.4−3.8], [3.8−4.3], [4.3−5.0], and [5−7] keV. Similar results were obtained by merging the central eight bands into just three. We computed energy-dependent exposure maps to calculate the effective exposure time, accounting for dithering, vignetting, CCD defects, gaps, and flagged pixels. We then measured counts and effective exposure time in the ten bands in circular annuli with increasing width with radius to counterbalance the decreasing intensity of the cluster. The minimum width is taken to be 3 arcsec, which is larger than the Chandra PSF. We only considered radii where the exposure time was longer than 50% of the on-axis exposure time and annuli included in the field of view by more than two-thirds. The cluster centre was iteratively computed as the centroid of X-ray emission within the inner 20 kpc, and the same centre was adopted for the SZ data. As background, we used blank field images (using CIAO BLANKSKY, Fruscione et al. 2006), normalised to the count rate in the hard band [9−13] keV, and we derived the background surface brightness in the ten bands accounting for exposure time variations, dithering, vignetting, CCD defects, gaps, and flagged pixels. In the X-ray data, information about the temperature profile is encoded in the ratio of the cluster count-rate profiles, while much of the metallicity information is contained in the [3.4−3.8] keV band.

As a proof of the capabilites of JoXSZ, we used data with different radial sampling between SZ and X-ray, the latter sampled on an irregular grid. Because the model is integrated in the same bins of the observations, results do not depend on binning (except for unreasonable choices, such as having one single bin, which completely loses spatial information).

4.2. Model definition

As defined by the equations in Sect. 3.1, many parameters are involved in the JoXSZ fit: the gNFW pressure profile has five parameters (P0,  rp,  a,  b, and c), the Vikhlinin density profile has 10 (ne0,  rc,  α,  β,  rs,  γ,  ϵ,  ne02,  rc2, and β2), to which need to be added the metallicity Z; the temperature profile ratio log(TX/TSZ); the parameter that accounts for the calibration uncertainty of the SZ measurement, implemented as a multiplicative component on the conversion; and the backscale parameter, which controls the scaling of the X-ray background. This yields a total maximum number of 19 parameters that can be fitted at the same time.

For our example, we cancelled out the additive component of the electron density profile by fixing ne02 = 0, we assumed α = 0 and γ = 3 following the reasoning of Mroczkowski et al. (2009), Comis et al. (2011), Adam et al. (2015), and we fixed c = 0.3081, in accordance with the universal pressure profile from Arnaud et al. (2010). In addition to these parameters, we accounted for the calibration and the backscale. We alternately considered the logarithm of the temperature profile ratio log(TX/TSZ) fixed to 0 or let free to vary, which means that we considered 13 parameters at maximum.

We took a Gaussian prior with σ = 0.07 centred on 1 for the SZ calibration (Adam et al. 2015), a Gaussian prior with σ = 0.1 centred on 1 for the backscale parameter, and uniform priors for all the remaining parameters, with ranges 0 <  P0 <  1, 100 <  rp <  1000, 0.5 <  a <  20, 0.5 <  b <  10, 0.0001 <  ne0 <  1, 0.1 <  rc <  rs <  2.5, 0 <  β <  4, 0 <  ϵ <  10, 0 <  Z <  1, and −1 <  log(TX/TSZ) < 1, where P0 is expressed in keV cm−3, ne0 in cm−3, and rp, rc, rs in kpc.

We ran 70 000 iterations of 30 walkers and considered the first 20 000 as the burn-in period. In order to reduce autocorrelation, we stored N/100 samples for each walker.

We considered an SZ pixel size of 2 arcsec (the PSF is ∼18 arcsec) and a cluster radial extent of 5 Mpc for the Abel integral computation. We fixed the absorbing column at the Galactic value (Kalberla et al. 2005). We adopted a flat ΛCDM cosmology with H0 = 67.32 km s−1 Mpc−1, ΩM = 0.3158 and ΩΛ = 0.6842 (Planck Collaboration VI 2020).

4.3. Results

We fitted the data several times to understand the effect of our assumptions on the derived thermodynamic profiles. We also compared JoXSZ performances with a restricted joint analysis and an X-ray alone analysis.

4.3.1. Disjointed temperatures TX and TSZ with or without HE

Our reference analysis assumes HE, allows X-ray and SZ temperatures to differ from each other, and has 13 free parameters (4 of which for the pressure profile, and 5 for the electron density profile). Figures 2 and B.1, automatically produced by JoXSZ, show the trace plot and joint plus marginal posterior distribution, respectively. They inform about the convergence of the chains and the parameter estimates. The acceptance fraction is found to be 0.11, which is lower than recommended; nevertheless, the wide thinning we adopted allowed us to have sufficiently uncorrelated samples. Some expected degeneracies are evident in the joint distributions: β and ϵ are highly negatively correlated, rc and β are positively correlated, as are rp and b, while ne0 is negatively correlated with rc. The a parameter is not well constrained by the data we considered. Figure 3, which is automatically generated as well, shows the best-fit profiles and their 95% uncertainty (equal-tailed credible interval) on top of the fitted data for each of the X-ray and SZ bands. Figure 4, also automatically produced, shows the main radial thermodynamic profiles of the cluster (median values and 95% intervals). The estimated entropy profile is in line with the Voit et al. (2005) fit to non-radiative simulations, as adapted by Pratt et al. (2010).

thumbnail Fig. 2.

Trace plot automatically produced by JoXSZ. Trace plots for four variables are only shown here, but are automatically produced for all parameters.

Open with DEXTER

thumbnail Fig. 3.

X-ray and SZ projected radial surface brightness profiles automatically produced by JoXSZ. Red lines show the best-fit profile. Shaded yellow areas represent 95% intervals. Blue points show X-ray data. Black points show SZ data. The X-ray profiles flatten off to the background value.

Open with DEXTER

Figure 4 shows the X-ray and SZ temperature profiles. The parameter log(TX/TSZ) is found to be , and its posterior distribution is shown in Fig. 5. The ratio between the X-ray and SZ temperatures has an unusual value, it is particularly far from 1 compared to what is expected from clumping effects as estimated from numerical expectations (Nagai & Lau 2011). This might be related to calibration systematics, such as the long-standing Chandra-XMM T systematics (Schellenberger et al. 2015), or possible inaccuracies of the transfer function (Romero et al. 2020). We emphasise that the conclusions accuracy does not depend exclusively on the correctness of operations that are performed on the data (the fitting code), but also on the lack of systematics in the input quantities (data, PSF, transfer function, instrument calibrations, etc.).

thumbnail Fig. 4.

Deprojected radial profiles for the main thermodynamic properties of the cluster (median with 95% intervals).

Open with DEXTER

thumbnail Fig. 5.

Posterior distribution of the temperature profile ratio parameter log(TX/TSZ). Median value and 95% intervals are highlighted.

Open with DEXTER

In absence of systematics, to estimate whether the data provide evidence in favour of differences between the X-ray and SZ temperatures, users may compute the Bayes factor from the JoXSZ output using the Savage-Dickey ratio, that is, the ratio between the posterior and prior density probabilities at log(TX/TSZ) = 0. In our case, we found a Bayes factor of 3.1 (= 1.53/(1/2)), which is almost inconclusive.

Figure 6 shows the mass profile estimated assuming HE and its 95% uncertainty. The intercept with 500ρc(z)V gives kpc and M, consistent with the measures obtained by Maughan et al. (2007), Mroczkowski et al. (2009), Mantz et al. (2010), and Adam et al. (2015). The profile shape is in accordance with the shape obtained by Adam et al. (2015). Relaxing the assumption of HE (Sect. 3.1.4) has barely any effect on our results: we performed a separate analysis without assuming a monotonically increasing mass profile, and only 0.1% of the samples returned unphysical values, likely because the parameters are sufficiently well constrained to work without this additional prior.

thumbnail Fig. 6.

Mass profile. The blue line shows the median mass profile derived assuming hydrostatic equilibrium (95% intervals shaded). Points with error bars represent r500 and M500 from the literature (references in the text).

Open with DEXTER

4.3.2. Restricted joint SZ+ne and X-ray alone

In a first fit, a restricted joint analysis SZ+ne, we fixed the metallicity to Z = 0.3 solar (Maughan et al. 2007; Arnaud et al. 2010), we discarded all X-ray photons except for those in the [0.7−1] keV band, but kept HE and log(TX/TSZ) = 0. This is meant to mimic some standard SZ+ne fits mentioned in Sect. 1 (for an exception, see Romero et al. 2017). Unlike most analyses, we used the correct Poisson-Gauss expression for the likelihood and a more inclusive error budget, and we fully propagated errors on the derived thermodynamic profiles.

The left panel of Fig. 7 shows our SZ gas-mass weighted posterior temperature profile. Because we discarded the information in the X-ray spectrum, the temperature derived here is between 1.1 and 1.3 times more uncertain than in the fully joint fit (i.e. more uncertain than in Fig. 4) for all r >  200 kpc. Adam et al. (2015) performed a restricted joint analysis of SZ data (of much the same NIKA data, supplemented by Planck) and electron density profile based on Chandra X-ray data. Our restricted joint SZ+ne posterior temperature profile agrees with theirs (points).

thumbnail Fig. 7.

Comparison of CL J1226.9+3332 deprojected TSZ (left panel) and TX (right panel) temperature profiles (median with 95% interval). The restricted joint analysis is shown in yellow-orange. Green shows JoXSZ. Red shows the X-ray analysis alone. The posterior derived by Adam et al. (2015) is marked by points.

Open with DEXTER

The second analysis completely excluded the SZ data and was performed with the original MBProj2 code. In particular, we used the same electron density parametrisation as adopted in JoXSZ and a simplified Vikhlinin et al. (2006) temperature profile, as in McDonald et al. (2014). The right panel of Fig. 7 compares the X-ray temperature profile derived in this way with the profile derived in our reference fully joint fit. The two profiles largely overlap because the free log(TX/TSZ) reduces the information passage between the SZ and X-ray parts of the analysis, and allows the JoXSZ X-ray temperature to agree with the temperature derived from X-ray data alone.

Figure 8 explains why the X-ray temperature is systematically lower than the SZ temperature, as highlighted in Fig. 4. The observed SZ surface brightness profile is higher in absolute value than the profile predicted based on the X-ray temperature and density profiles alone.

thumbnail Fig. 8.

Comparison of CL J1226.9+3332 projected SZ surface brightness profiles (median with 95% interval). Blue shows JoXSZ. Red presents predictions from the X-ray alone analysis. Points with 68% error bars show the SZ data from NIKA.

Open with DEXTER

5. Conclusion

The recent spread of SZ observations, together with the vast availability of X-ray data, caused joint SZ+X analyses of galaxy clusters to flourish. The need for a public tool to perform such analyses was shared among the community of workers within the field. We presented JoXSZ, the first publicly available code that combines the multiwavelength SZ+X approach and the X-ray multiband fitting technique in a full and consistent way. JoXSZ supports data coming from multiple instruments and can even handle missing data for either SZ or X-ray surface brightness, follows a Bayesian forward-modelling approach, and adopts flexible parametrisations of the thermodynamic cluster profiles. It makes full use of the information contained in the observations, it uses the correct Poisson-Gauss expression for the joint likelihood, accounts for beam smearing and transfer function, includes SZ calibration and X-ray background systematics, adopts a consistent temperature in the various parts of the code, and allows differences between SZ and X-ray temperatures, for example if users wish to study cluster elongation or clumping. By also fitting the cluster metallicity, JoXSZ accounts for the metallicity dependence of the X-ray conversion factor. It supports the use of different radial binnings for SZ data and X-ray data and consents either to adopting or to relaxing the assumption of HE on user request. For these reasons, JoXSZ goes beyond simple SZ and electron density fits. When HE holds, JoXSZ uses a physical (positive) prior on the derivative of the mass profile and derives the mass profile and overdensity radii rΔ. As in other approaches, JoXSZ makes the usual assumptions about cluster sphericity, and when requested, HE.

JoXSZ returns convergence diagnostics, parameter estimates with uncertainties, joint and marginal probability distributions, projected radial distributions, and three-dimensional thermodynamic radial profiles.

The code is fully documented, and the users are free to customise their analysis in accordance with their needs and requirements. JoXSZ has been released as an open-source Python project, and its code is publicly available on GitHub.

We also provided an application using real data of the high-redshift galaxy cluster CL J1226.9+3332 to illustrate the features of JoXSZ. When compared to a restricted joint SZ+ne fit or an X-ray alone fit, JoXSZ derives smaller uncertainties by accounting for the whole information of SZ and multiband X-ray data. While the uncertainty reduction is not striking for our case, a more relevant improvement may occur in different circumstances. The easiness with which thermodynamic profiles are derived with JoXSZ allows the users to focus on astronomically interesting features.

We plan to further develop the code in the near future to also fit the shear (i.e. weak-lensing information). The execution time of the SZ part of the analysis will also be improved.


We warmly thank Charles Romero, Roberto Scaramella and Luca Di Mascolo for their useful comments on the manuscript. We are also thankful to the anonymous referee for the detailed revision that helped improve the quality of the paper. F.C. acknowledges financial contribution from the agreement ASI-INAF n.2017-14-H.0 and PRIN MIUR 2015 Cosmology and Fundamental Physics: Illuminating the Dark Universe with Euclid.


  1. Adam, R., Comis, B., Macías-Pérez, J. F., et al. 2015, A&A, 576, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Ameglio, S., Borgani, S., Pierpaoli, E., & Dolag, K. 2007, MNRAS, 382, 397 [NASA ADS] [CrossRef] [Google Scholar]
  3. Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
  4. Andreon, S., Moretti, A., Trinchieri, G., & Ishwara-Chandra, C. H. 2019, A&A, 630, A78 [CrossRef] [EDP Sciences] [Google Scholar]
  5. Arnaud, K. A. 1996, in XSPEC: The First Ten Years, eds. G. H. Jacoby, & J. Barnes, ASP Conf. Ser., 101, 17 [Google Scholar]
  6. Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Birkinshaw, M., & Lancaster, K. 2005, in Background Microwave Radiation and Intracluster Cosmology, eds. F. Melchiorri, & Y. Rephaeli, 127 [Google Scholar]
  8. Bocquet, S., Saro, A., Mohr, J. J., et al. 2015, ApJ, 799, 214 [NASA ADS] [CrossRef] [Google Scholar]
  9. Castagna, F., & Andreon, S. 2019, A&A, 632, A22 [CrossRef] [EDP Sciences] [Google Scholar]
  10. Catalano, A., Calvo, M., Ponthieu, N., et al. 2014, A&A, 569, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Cavaliere, A., & Fusco-Femiano, R. 1976, A&A, 500, 95 [NASA ADS] [Google Scholar]
  12. Comis, B., de Petris, M., Conte, A., Lamagna, L., & de Gregori, S. 2011, MNRAS, 418, 1089 [NASA ADS] [CrossRef] [Google Scholar]
  13. Donahue, M., Voit, G. M., Mahdavi, A., et al. 2014, ApJ, 794, 136 [NASA ADS] [CrossRef] [Google Scholar]
  14. Ebeling, H., Jones, L. R., Fairley, B. W., et al. 2001, ApJ, 548, L23 [NASA ADS] [CrossRef] [Google Scholar]
  15. Eckert, D., Molendi, S., Gastaldello, F., & Rossetti, M. 2011, A&A, 529, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Eckert, D., Molendi, S., Vazza, F., Ettori, S., & Paltani, S. 2013, A&A, 551, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Ettori, S., Donnarumma, A., Pointecouteau, E., et al. 2013, Space Sci. Rev., 177, 119 [NASA ADS] [CrossRef] [Google Scholar]
  18. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [NASA ADS] [CrossRef] [Google Scholar]
  19. Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in CIAO: Chandra’s Data Analysis System, SPIE Conf. Ser., 6270, 62701V [Google Scholar]
  20. Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
  22. Hudson, D. S., Mittal, R., Reiprich, T. H., et al. 2010, A&A, 513, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Itoh, N., & Nozawa, S. 2004, A&A, 417, 827 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Kitayama, T., Komatsu, E., Ota, N., et al. 2004, PASJ, 56, 17 [NASA ADS] [Google Scholar]
  26. Korngut, P. M., Dicker, S. R., Reese, E. D., et al. 2011, ApJ, 734, 10 [NASA ADS] [CrossRef] [Google Scholar]
  27. Mantz, A., Allen, S. W., Ebeling, H., Rapetti, D., & Drlica-Wagner, A. 2010, MNRAS, 406, 1773 [NASA ADS] [Google Scholar]
  28. Markevitch, M., & Vikhlinin, A. 2007, Phys. Rep., 443, 1 [NASA ADS] [CrossRef] [Google Scholar]
  29. Maughan, B. J., Jones, L. R., Ebeling, H., & Scharf, C. 2004, MNRAS, 351, 1193 [NASA ADS] [CrossRef] [Google Scholar]
  30. Maughan, B. J., Jones, C., Jones, L. R., & Van Speybroeck, L. 2007, ApJ, 659, 1125 [NASA ADS] [CrossRef] [Google Scholar]
  31. McDonald, M., Benson, B. A., Vikhlinin, A., et al. 2014, ApJ, 794, 67 [NASA ADS] [CrossRef] [Google Scholar]
  32. Mroczkowski, T., Bonamente, M., Carlstrom, J. E., et al. 2009, ApJ, 694, 1034 [NASA ADS] [CrossRef] [Google Scholar]
  33. Mroczkowski, T., Nagai, D., Basu, K., et al. 2019, Space Sci. Rev., 215, 17 [NASA ADS] [CrossRef] [Google Scholar]
  34. Nagai, D., & Lau, E. T. 2011, ApJ, 731, L10 [NASA ADS] [CrossRef] [Google Scholar]
  35. Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1 [NASA ADS] [CrossRef] [Google Scholar]
  36. Planck Collaboration VI. 2020, A&A, in press, [Google Scholar]
  37. Pratt, G. W., Arnaud, M., Piffaretti, R., et al. 2010, A&A, 511, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Romero, C. E., Mason, B. S., Sayers, J., et al. 2017, ApJ, 838, 86 [NASA ADS] [CrossRef] [Google Scholar]
  39. Romero, C., McWilliam, M., Macías-Pérez, J. F., et al. 2018, A&A, 612, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Romero, C. E., Sievers, J., Ghirardini, V., et al. 2020, ApJ, 891, 90 [CrossRef] [Google Scholar]
  41. Ruppin, F., Mayet, F., Macías-Pérez, J. F., & Perotto, L. 2019, MNRAS, 490, 784 [NASA ADS] [CrossRef] [Google Scholar]
  42. Ruppin, F., McDonald, M., Brodwin, M., et al. 2020, ApJ, 893, 74 [CrossRef] [Google Scholar]
  43. Sanders, J. S., Fabian, A. C., Russell, H. R., & Walker, S. A. 2018, MNRAS, 474, 1065 [NASA ADS] [CrossRef] [Google Scholar]
  44. Sayers, J., Czakon, N. G., Mantz, A., et al. 2013, ApJ, 768, 177 [NASA ADS] [CrossRef] [Google Scholar]
  45. Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Shitanishi, J. A., Pierpaoli, E., Sayers, J., et al. 2018, MNRAS, 481, 749 [NASA ADS] [CrossRef] [Google Scholar]
  47. Siegel, S. R., Sayers, J., Mahdavi, A., et al. 2018, ApJ, 861, 71 [CrossRef] [Google Scholar]
  48. Sun, M., Voit, G. M., Donahue, M., et al. 2009, ApJ, 693, 1142 [NASA ADS] [CrossRef] [Google Scholar]
  49. Sunyaev, R. A., & Zeldovich, Y. B. 1970, Ap&SS, 7, 3 [NASA ADS] [Google Scholar]
  50. Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comm. Astrophys. Space Phys., 4, 173 [Google Scholar]
  51. Trotta, R. 2007, MNRAS, 378, 72 [NASA ADS] [CrossRef] [Google Scholar]
  52. Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691 [NASA ADS] [CrossRef] [Google Scholar]
  53. Voit, G. M. 2005, Rev. Mod. Phys., 77, 207 [NASA ADS] [CrossRef] [Google Scholar]
  54. Voit, G. M., Kay, S. T., & Bryan, G. L. 2005, MNRAS, 364, 909 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Creation of Python classes pressure and SZ data and joint likelihood definition

The JoXSZ code is grafted onto the structure designed for the MBProj2 code by Sanders et al. (2018). To fit X-ray surface brightness profiles alone, MBProj2 allows users to model one or more thermodynamic profiles except for the pressure profile, and it does not support any fit to SZ data.

The novelty we introduced in JoXSZ is the inclusion and treatment of the SZ data and the possibility of parametrising the pressure profile by creating a specific class for it. The functions within the class replicate the standard functions for the other thermodynamic profiles: defPars defines the default parameter values, press_fun computes the analytic expression of the pressure profile, and press_derivative returns the analytically computed first derivative of it. Similarly, a new class for the temperature is defined as the ratio between pressure and electron density.

We also created a specific Python class to recollect all the elements required to perform the fit on SZ data that do not change in the iterations: the sampling step and the corresponding radius vector, the matrices of beam and transfer function built from observed or approximated data, the conversion factors from Compton y to surface brightness, and the observed SZ data.

As outlined in Fig. 1 and described in Sect. 3.3, the likelihood function ℒJoXSZ is defined as the product of ℒX from MBProj2 and ℒSZ from PreProFit, the former updated for ignoring missing values and the latter updated to allow the temperature dependence of the Compton y to surface brightness conversion. Prior to the likelihood calculation, if the unphysical mass exclusion is set, the mass profile is evaluated according to Eq. (4) and ℒJoXSZ is properly set to −∞.

Appendix B: Joint and marginal posterior probability distributions

thumbnail Fig. B.1.

Joint and marginal posterior distributions that are automatically produced by JoXSZ. Dashed lines in the diagonal panels indicate median values.

Open with DEXTER

All Figures

thumbnail Fig. 1.

Block diagram showing the program flow. Radial profiles are plotted in red. Options are shown in green. Analysis pipelines are depicted in yellow. Data enter in the blue box.

Open with DEXTER
In the text
thumbnail Fig. 2.

Trace plot automatically produced by JoXSZ. Trace plots for four variables are only shown here, but are automatically produced for all parameters.

Open with DEXTER
In the text
thumbnail Fig. 3.

X-ray and SZ projected radial surface brightness profiles automatically produced by JoXSZ. Red lines show the best-fit profile. Shaded yellow areas represent 95% intervals. Blue points show X-ray data. Black points show SZ data. The X-ray profiles flatten off to the background value.

Open with DEXTER
In the text
thumbnail Fig. 4.

Deprojected radial profiles for the main thermodynamic properties of the cluster (median with 95% intervals).

Open with DEXTER
In the text
thumbnail Fig. 5.

Posterior distribution of the temperature profile ratio parameter log(TX/TSZ). Median value and 95% intervals are highlighted.

Open with DEXTER
In the text
thumbnail Fig. 6.

Mass profile. The blue line shows the median mass profile derived assuming hydrostatic equilibrium (95% intervals shaded). Points with error bars represent r500 and M500 from the literature (references in the text).

Open with DEXTER
In the text
thumbnail Fig. 7.

Comparison of CL J1226.9+3332 deprojected TSZ (left panel) and TX (right panel) temperature profiles (median with 95% interval). The restricted joint analysis is shown in yellow-orange. Green shows JoXSZ. Red shows the X-ray analysis alone. The posterior derived by Adam et al. (2015) is marked by points.

Open with DEXTER
In the text
thumbnail Fig. 8.

Comparison of CL J1226.9+3332 projected SZ surface brightness profiles (median with 95% interval). Blue shows JoXSZ. Red presents predictions from the X-ray alone analysis. Points with 68% error bars show the SZ data from NIKA.

Open with DEXTER
In the text
thumbnail Fig. B.1.

Joint and marginal posterior distributions that are automatically produced by JoXSZ. Dashed lines in the diagonal panels indicate median values.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.