Probing the baryon content of a galaxy cluster

The aim of this tutorial is to take you through the basic methods with which XGA can probe the hot gas that makes of the intra-cluster medium (ICM). We’ll focus on measuring the radial gas density profiles of a sample of galaxy clusters here, then move on to using density profiles to measure gas masses.

Please note that I will not really be fully explaining how the methods to measure gas density from clusters work in this tutorial, that is the subject of a section in the ‘Under the Hood’ section, where the maths and concepts are detailed.

from xga.samples import ClusterSample
from xga.sourcetools.density import inv_abel_fitted_model, ann_spectra_apec_norm
from xga.imagetools.psf import rl_psf
from xga.models import KingProfile1D, SimpleVikhlininDensity1D

from astropy.units import Quantity
import numpy as np
import pandas as pd

Defining a sample of galaxy clusters

I’m using the same sample of galaxy clusters that I used in the spectroscopy tutorial for the sake of convenience. I’ve chosen to use a ClusterSample rather than a GalaxyCluster as there are extra, convenient, methods built into ClusterSamples for dealing with samples of gas masses, though I will also show how you interact with gas density profiles/gas masses for an individual cluster:

# Setting up the column names and numpy array that go into the Pandas dataframe
column_names = ['name', 'ra', 'dec', 'z', 'r500', 'r200', 'richness', 'richness_err']
cluster_data = np.array([['XCSSDSS-124', 0.80057775, -6.0918182, 0.251, 1220.11, 1777.06, 109.55, 4.49],
                         ['XCSSDSS-2789', 0.95553986, 2.068019, 0.11, 1039.14, 1519.79, 38.90, 2.83],
                         ['XCSSDSS-290', 2.7226392, 29.161021, 0.338, 935.58, 1359.37, 105.10, 5.99],
                         ['XCSSDSS-134', 4.9083898, 3.6098177, 0.273, 1157.04, 1684.15, 108.60, 4.79]])

# Possibly I'm overcomplicating this by making it into a dataframe, but it is an excellent data structure,
#  and one that is very commonly used in my own analyses.
sample_df = pd.DataFrame(data=cluster_data, columns=column_names)
sample_df[['ra', 'dec', 'z', 'r500', 'r200', 'richness', 'richness_err']] = \
    sample_df[['ra', 'dec', 'z', 'r500', 'r200', 'richness', 'richness_err']].astype(float)

# Defining the sample of four XCS-SDSS galaxy clusters
demo_smp = ClusterSample(sample_df["ra"].values, sample_df["dec"].values, sample_df["z"].values,
                         sample_df["name"].values, r200=Quantity(sample_df["r200"].values, "kpc"),
                         r500=Quantity(sample_df["r500"].values, 'kpc'), richness=sample_df['richness'].values,
Declaring BaseSource Sample: 100%|██████████| 4/4 [00:01<00:00,  2.74it/s]
Generating products of type(s) ccf: 100%|██████████| 4/4 [00:21<00:00,  5.45s/it]
Generating products of type(s) image: 100%|██████████| 4/4 [00:00<00:00,  8.05it/s]
Generating products of type(s) expmap: 100%|██████████| 4/4 [00:00<00:00,  8.64it/s]
Setting up Galaxy Clusters: 100%|██████████| 4/4 [00:04<00:00,  1.21s/it]

Surface brightness based gas density measurement

This first method uses the surface brightness profile of a cluster, as well as assumptions about its geometry and plasma emission, to measure the three-dimensional gas density profile. There are two broad steps to this process; we make a surface brightness profile of the cluster, fit a model to it, then use that to deduce the three-dimensional emissivity and then we use a plasma emission model, combined with an absorption model, to convert that emissivity into a density.

Of course the implementation (and details) of this method are more complicated than that makes it sound, but a convenience function has been added to XGA that goes through the whole process, inv_abel_fitted_model(), for which the full docstring can be found here.

PSF correction

As we are likely to be using small annular bins for the brightness profile, and the XMM-Newton telescope has a non-negligible point spread function (PSF), we have to account for this effect with a correction, bearing in mind that these effects also vary spatially. The PSFs of the three XMM cameras have been well studied, with Read et al. (2011) producing comprehensive models that take into account the complex two-dimensional structure of the effects. In XGA we can use those models to correct images directly, and prior to the generation of a brightness profile, whereas previously a one-dimensional model was used to correct after the profile had been generated.

A full explanation of how this process works can be found in PSF correction of XMM images in the ‘Under the Hood’ section, and the documentation for the correction function is here. Please note that the inv_abel_fitted_model() function will automatically generate and use PSF corrected images, but I’m calling the correction function here manually as a demonstration (you can set custom chunking and iteration factors in the function call but I’m leaving it default). You can also automatically run PSF correction when a ClusterSample is declared by setting psf_corr=True:

Preparing PSF generation commands: 100%|██████████| 4/4 [00:00<00:00, 15.82it/s]
Generating products of type(s) psf: 100%|██████████| 12/12 [00:04<00:00,  2.77it/s]
PSF Correcting Observations - Currently complete: 100%|██████████| 4/4 [00:35<00:00,  8.82s/it]
Generating products of type(s) image: 100%|██████████| 4/4 [00:00<00:00,  7.76it/s]

We can see the effect of this correction on the images and brightness profiles, with a demonstration using the first cluster in the sample (XCSSDSS-124):

rt = demo_smp[0].get_combined_ratemaps(psf_corr=False)

rt = demo_smp[0].get_combined_ratemaps(psf_corr=True)

Running the density calculation function

There are many arguments that can be passed to this function, and I’m not going to explain them all here, please see the function documentation for explanations. However, I will go through the most important arguments you need to consider when running this function:

  • model - This is the model that will be fitted to the surface brightness profile, and then transformed to infer the three-dimensional countrate per volume element radial distribution. This can either be a model name, a list of model names (if dealing with a sample), a model instance (see profile tutorial for context), or a list of model instances (if dealing with a sample).

  • outer_radius - The radius (or radii, if we’re dealing with a sample) out to which the profiles should be generated. This has to be carefully balanced, as you don’t want to choose too small a radius and have too much cluster emission emission in the background annulus, and you don’t want to choose too large a radius and start hitting the edge of the detector.

  • num_dens - If this is true then a number density profile will be calculated, if False then a mass density profile.

  • obs_id and inst - You can pass an ObsID and instrument (or lists of these, if dealing with a sample) if you wish to specify a single image to measure density with.

  • conv_temp - The temperature used for the plasma emission model can be specified manually (for samples as well as individual clusters). The default behaviour is to measure the temperature within the conv_outer_radius and use that, and if that fit fails then to assume the temperature is 3keV.

  • conv_outer_radius - The radius (radii for samples) for which spectra are generated (and fitted if conv_temp isn’t set) to calculate the conversion between count-rate per volume element and density, which involves the use of the ARFs (see the spectroscopic tutorial for an explanation)

Here I generate profiles out to 1.2R\(_{500}\) of each cluster, and use spectra generated within R\(_{500}\) to measure the temperature and calculate the conversion. I also choose to fit a beta profile, and increase the number of steps the MCMC samplers will take from the default (20000) to 40000. This function also returns a list of profiles (if a profile fails to generate for some reason a None will take the profile’s place in the list), as well as storing the profiles within the source objects in the sample:

sb_d_profs = inv_abel_fitted_model(demo_smp, 'beta', outer_radius=1.2*demo_smp.r500, conv_outer_radius='r500',
Generating products of type(s) spectrum: 100%|██████████| 12/12 [38:34<00:00, 192.89s/it]
Running XSPEC Fits: 100%|██████████| 4/4 [00:20<00:00,  5.17s/it]
Running XSPEC Simulations: 100%|██████████| 4/4 [00:02<00:00,  1.50it/s]
Fitting data, inverse Abel transforming, and measuring densities: 100%|██████████| 4/4 [07:42<00:00, 115.59s/it]
[<xga.products.profile.GasDensity3D at 0x7f0890f79970>,
 <xga.products.profile.GasDensity3D at 0x7f0890f79ca0>,
 <xga.products.profile.GasDensity3D at 0x7f0890f79b50>,
 <xga.products.profile.GasDensity3D at 0x7f0890f79e20>]

Annular spectra gas density method

The alternative method is based on a similar idea, but instead we generate annular spectra, fit an absorbed plasma model to each annulus, and use the normalisation of that model (and assumptions about the cluster geometry) to measure the density. The main disadvantage of this method is that you require a lot of X-ray counts in a spectrum to get a good fit, so the annuli must be much larger than the brightness profile method. However, for clusters with high enough quality data this method can be a check of the results measured from the surface brightness profile method.

This method doesn’t really require any PSF correction due to the large annuli required for the annular spectra to have enough counts for a successful fit, and we will set a minimum size to ensure that PSF considerations are unnecessary. Again there are more arguments for this function than I want to go into here, but you can read about all of them in the full function documentation. I will again take you through the most important arguments, however:

  • outer_radius - The radius (or radii, if we’re dealing with a sample) out to which the profiles should be generated. This has to be carefully balanced, as you don’t want to choose too small a radius and have too much cluster emission emission in the background annulus, and you don’t want to choose too large a radius and start hitting the edge of the detector.

  • num_dens - If this is true then a number density profile will be calculated, if False then a mass density profile.

  • annulus_method - The method by which the annuli are designated, the default is a minimimum signal to noise value (min_snr) which is targetted, though if it cannot be attained in all annuli the function will allow the number of annuli to drop to four, use them, and give a warning to the user.

  • min_snr - The minimum signal to noise target for annuli.

  • min_width - The minimum width of each annulus, default is 20 arcseconds to try and remove the need for accounting for PSF effects.

  • use_worst - If set to True then the annuli generation algorithm will use the ‘worst’ observation (defined by global signal to noise of the cluster), to decide how big to make the annuli.

as_d_profs = ann_spectra_apec_norm(demo_smp, 1.2*demo_smp.r500, min_snr=30)
/home/dt237/code/PycharmProjects/XGA/xga/sourcetools/ UserWarning: The requested annuli for XCSSDSS-290 cannot be created, the data quality is too low. As such a set of four annuli will be returned
  warn("The requested annuli for {s} cannot be created, the data quality is too low. As such a set "
Generating products of type(s) spectrum: 100%|██████████| 12/12 [50:47<00:00, 253.96s/it]
Generating products of type(s) annular spectrum set components: 100%|██████████| 129/129 [29:49<00:00, 13.87s/it]
Running XSPEC Fits: 100%|██████████| 43/43 [01:12<00:00,  1.69s/it]
Generating density profiles from annular spectra:  25%|██▌       | 1/4 [00:08<00:26,  8.85s/it]/home/dt237/software/anaconda3/envs/conda4-xcsim/lib/python3.8/site-packages/astropy/units/ RuntimeWarning: invalid value encountered in sqrt
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
Generating density profiles from annular spectra: 100%|██████████| 4/4 [00:36<00:00,  9.06s/it]
[<xga.products.profile.GasDensity3D at 0x7f0890f8fac0>,
 <xga.products.profile.GasDensity3D at 0x7f08913a3700>,
 <xga.products.profile.GasDensity3D at 0x7f0890f8b2b0>,
 <xga.products.profile.GasDensity3D at 0x7f0890f8b130>]
$[0.000784351,~0.00154955,~0.00140047,~0.00120906,~0.00072586,~0.000546847,~0.000442686,~0.000395873,~0.000489743,~0.000352902,~0.000286331,~0.000187971] \; \mathrm{\frac{1}{cm^{5}}}$

Viewing the gas density profiles

The profile tutorial demonstrated how to use the profile view method to make a visualisation of the radial profile data, as well as the fact that we can add two compatible profiles together to compare them. I will use that ability again here to view (and compare) the density profiles from the different methods. At first I will simply draw the profiles from the lists which were outputted by the measurement functions, but then I will demonstrate how to retrieve a profile directly from a source:

for src_ind in range(0, len(demo_smp)):

We can also view the surface brightness profiles that went into the creation of the density profiles from the first method, by cycling through them and using the generation_profile property:

for d_prof in sb_d_profs:

With the same property, generation_profile, we can actually have a look at the APEC normalisation profiles that we generated the annular spectra based density profiles from:

for d_prof in as_d_profs:

Measuring a gas mass from density profiles

For a gas mass to be measured, we first must fit a model to the density profile, again using the profile fit() method. After this process has run through, we will then be able to request a gas mas at any radius we like, and have both an individual gas mass/uncertainty returned, and a gas mass distribution.

I am choosing to fit a simple King function to the density profile, and though I have imported the KingProfile1D class, I’ve only done this to show you the output of the info() method, and I won’t actually be setting up a custom instance:

| King Profile    |                                                                         |
| DESCRIBES       | Gas Density                                                             |
| UNIT            | solMass / Mpc3                                                          |
| PARAMETERS      | beta, r_core, norm                                                      |
| PARAMETER UNITS | , kpc, solMass / Mpc3                                                   |
| AUTHOR          | placeholder                                                             |
| YEAR            | placeholder                                                             |
| PAPER           | placeholder                                                             |
| INFO            | The un-projected version of the beta profile, suitable for a simple fit |
|                 |  to 3D density distributions. Describes a simple isothermal sphere.     |

Now I’m going to demonstrate how you can retrieve a gas density profile from a source object directly, rather than relying on the list output by the density profile generation functions. When you’ve only generated one density profile for a given source, just calling get_density_profiles() will return that for you without any extra information being passed in. If you’ve generated multiple profiles (be it with different methods or simply for different outermost radii) like we’ve done in this demonstration, then it will return a list containing all of the available profiles for the source:

# Just calling it by itself gets you both of the density profiles generated for this particular source
#  Not necessarily what we want!
[<xga.products.profile.GasDensity3D at 0x7f0890f79970>,
 <xga.products.profile.GasDensity3D at 0x7f0890f8fac0>]

As such, we’ll need to pass in some extra information to retrieve the exact profile we want to fit a model to. Seeing as you can generate these profiles in various different ways, the get method can take quite a lot of information to retrieve the exact profiles you want. The outermost radius can be quite a good one to use, but in this case, where we have two different profiles generated to the same outer radius but with two different methods, we can set the method argument. If, for instance, we’ve decided we want to use the profile that was generated by fitting the surface brightness profile then we would pass method='beta', where beta is the name of the model that was fitted to the profile. We, for no particular reason, are going to use the profile generated from annular spectra:

cur_dprof = demo_smp[1].get_density_profiles(method='spec')
<xga.products.profile.GasDensity3D at 0x7f08913a3700>

Now that we’ve successfully retrieved the density profile we want to use, we’re going to fit a model to it. So to start with we’re going to see what models are allowed for this type of profile:

| MODEL NAME            | EXPECTED PARAMETERS                          | DEFAULT START VALUES                            |
| simple_vikhlinin_dens | beta, r_core, alpha, r_s, epsilon, norm      | 1.0, 100.0 kpc, 1.0, 300.0 kpc, 2.0,            |
|                       |                                              | 10000000000000.0 solMass / Mpc3                 |
| king                  | beta, r_core, norm                           | 1.0, 100.0 kpc, 10000000000000.0 solMass / Mpc3 |
| vikhlinin_dens        | beta_one, r_core_one, alpha, r_s, epsilon,   | 1.0, 100.0 kpc, 1.0, 300.0 kpc, 2.0,            |
|                       | gamma, norm_one, beta_two, r_core_two,       | 3.0, 10000000000000.0 solMass / Mpc3,           |
|                       | norm_two                                     | 1.0, 50.0 kpc, 5000000000000.0 solMass / Mpc3   |

I’m going to choose to fit a simple Vikhlinin density model, which I’ve imported at the top of this tutorial so I can call the info() method to give you some extra details about it:

| Simplified Vikhlinin Profile   |                                                                       |
| DESCRIBES                      | Gas Density                                                           |
| UNIT                           | solMass / Mpc3                                                        |
| PARAMETERS                     | beta, r_core, alpha, r_s, epsilon, norm                               |
| PARAMETER UNITS                | , kpc, , kpc, , solMass / Mpc3                                        |
| AUTHOR                         | Ghirardini et al.                                                     |
| YEAR                           | 2019                                                                  |
| PAPER                          |                           |
| INFO                           | A simplified form of Vikhlinin's full density model, a type of broken |
|                                |  power law that deals well with most galaxy cluster density profile.  |

Now we’ll run the actual fit using MCMC, then view the model:

fit_mod ='simple_vikhlinin_dens', num_steps=50000, method='mcmc', show_warn=False)
100%|██████████| 50000/50000 [02:09<00:00, 387.30it/s]

Now that we have successfully fitted a model to the chosen density profile, we can very easily measure a gas mass by calling the gas_mass() method that has been built into the gas density profile class. All we need to tell it is the name of the model which we fitted (this has to be specified as fitting multiple models to a single profile is supported) and the radius within which we would like to know the gas mass. In this case I have chosen to calculate the gas mass within the R\(_{500}\) of this cluster.

Two objects are returned by this method, the first is the median gas mass, with - and + uncertainties at the 68.2% confidence level (~1:math:sigma), and the second is the entire gas mass distribution calculated by the method. If you didn’t wish to use that particular confidence level, you could use the conf_level keyword argument to change it:

gmass, gmass_dist = cur_dprof.gas_mass('simple_vikhlinin_dens', demo_smp[1].r500)
$[4.3807353 \times 10^{13},~6.5393064 \times 10^{11},~6.7858617 \times 10^{11}] \; \mathrm{M_{\odot}}$
$[4.3468629 \times 10^{13},~4.4912317 \times 10^{13},~4.3363341 \times 10^{13},~\dots,~4.3413252 \times 10^{13},~4.4216002 \times 10^{13},~4.3484162 \times 10^{13}] \; \mathrm{M_{\odot}}$

Viewing the gas mass distribution at a given radius

A convenient view method for the gas mass distributions calculated by this object has been implemented, view_gas_mass_dist(). It takes the same arguments as gas_mass() and simply plots the mass distribution output as a histogram. This seems a good time to point out that any gas mass values and distribution that has been calculated by a density profile (at a specific radius) is stored for later use, so if a gas mass at the same radius (from the same model) is requested again later then the relatively expensive calculated doesn’t have to be repeated.

cur_dprof.view_gas_mass_dist('simple_vikhlinin_dens', demo_smp[1].r500)

You can request a gas mass within any radius you like:

cur_dprof.view_gas_mass_dist('simple_vikhlinin_dens', Quantity(400, 'kpc'))

Even if that radius is outside the scope of the data which was used to measure the density profile in the first place (though of course you should always be cautious when extrapolating from models):

cur_dprof.view_gas_mass_dist('simple_vikhlinin_dens', Quantity(2000, 'kpc'))

Generating gas mass profiles

As I just demonstrated, we can measure gas masses at whatever radius we like, which of course means that it is possible to create a profile of how the cumulative gas mass changes with radius. This has been implemented as another method of the GasDensity3D class, and all you need to supply is the name of the model from which the profile should be generated, and the radii where you wish to measure the values.

Here I choose to use the same radii at which the gas densities were measured (this is the default behaviour when you don’t pass custom radii and deg_radii values), but you could use whatever radii you like:

gmass_prof = cur_dprof.gas_mass_profile('simple_vikhlinin_dens')
/home/dt237/code/PycharmProjects/XGA/xga/products/ UserWarning: The outer radius you supplied is greater than or equal to the outer radius covered by the data, so you are effectively extrapolating using the model.
  warn("The outer radius you supplied is greater than or equal to the outer radius covered by the data, so"

You have to pass both the radii and the equivalent radii in degrees, as the profile is unable to convert proper radii for itself, as it has no knowledge of the chosen cosmology and the redshift of the source. This was by design, as all XGA products were meant to be able to be defined completely separate from a source object. If you did want to choose custom radii however, you could set them up like this:

rads = Quantity([100, 200, 300, 400, 500, 600, 700], 'kpc')
deg_rads = demo_smp[1].convert_radius(rads, 'deg')
$[0.013412552,~0.026825104,~0.040237656,~0.053650209,~0.067062761,~0.080475313,~0.093887865] \; \mathrm{{}^{\circ}}$

Here the convert_radius() method of the source we generated the density profile from is used to convert the kpc radii to degrees, as it has knowledge of the cluster redshift and the cosmology that we chose to analyse it.

Now we can take a quick look at the gas mass profile we just made:


Retrieving a sample of gas masses

Finally, I will demonstrate how to retrieve gas masses from a whole sample of density profiles. The ClusterSample class has a method called gas_mass() built in which will iterate through all the clusters in the sample, attempt to retrieve the specified gas density profile, and then calculate the specified gas mass.

You have to tell the method which density model you want to measure the gas mass for, and in this case I’ve chosen the King model to demonstrate that if you haven’t already fitted a model then it will attempt to do so for you (using the default settings for the fit() method of the profile.

demo_smp.gas_mass('r500', 'king', method='beta')
100%|██████████| 30000/30000 [01:14<00:00, 403.91it/s]
100%|██████████| 30000/30000 [01:22<00:00, 362.25it/s]
100%|██████████| 30000/30000 [01:11<00:00, 417.61it/s]
100%|██████████| 30000/30000 [01:13<00:00, 408.02it/s]
$[[1.1816284 \times 10^{14},~8.825163 \times 10^{10},~8.5862588 \times 10^{10}],~ [2.9581939 \times 10^{13},~3.9242556 \times 10^{10},~4.0931691 \times 10^{10}],~ [4.8668179 \times 10^{13},~1.9847375 \times 10^{11},~2.0144428 \times 10^{11}],~ [9.6621655 \times 10^{13},~1.0204632 \times 10^{11},~1.0242727 \times 10^{11}]] \; \mathrm{M_{\odot}}$