# Mock Observations of the Sunyaev-Zeldovich Effect¶

Notebook

The change in the CMB intensity due to Compton scattering of CMB photons off of thermal electrons in galaxy clusters, otherwise known as the Sunyaev-Zeldovich (S-Z) effect, can to a reasonable approximation be represented by a projection of the pressure field of a cluster. However, the full S-Z signal is a combination of thermal and kinetic contributions, and for large frequencies and high temperatures relativistic effects are important. For computing the full S-Z signal incorporating all of these effects, there is a library: SZpack (Chluba et al 2012).

The sunyaev_zeldovich analysis module in yt makes it possible to make projections of the full S-Z signal given the properties of the thermal gas in the simulation using SZpack. SZpack has several different options for computing the S-Z signal, from full integrations to very good approximations. Since a full or even a partial integration of the signal for each cell in the projection would be prohibitively expensive, we use the method outlined in Chluba et al 2013 to expand the total S-Z signal in terms of moments of the projected optical depth $\tau$, projected electron temperature $T_e$, and velocities $\beta_{c,\parallel}$ and $\beta_{c,\perp}$ (their equation 18):

$$S(\tau, T_{e},\beta_{c,\parallel},\beta_{\rm c,\perp}) \approx S_{\rm iso}^{(0)} + S_{\rm iso}^{(2)}\omega^{(1)} + C_{\rm iso}^{(1)}\sigma^{(1)} + D_{\rm iso}^{(2)}\kappa^{(1)} + E_{\rm iso}^{(2)}\beta_{\rm c,\perp,SZ}^2 +~...$$

yt makes projections of the various moments needed for the calculation, and then the resulting projected fields are used to compute the S-Z signal. In our implementation, the expansion is carried out to first-order terms in $T_e$ and zeroth-order terms in $\beta_{c,\parallel}$ by default, but terms up to second-order in can be optionally included.

## Installing SZpack¶

SZpack can be downloaded here. Make sure you install a version later than v1.1.1. For computing the S-Z integrals, SZpack requires the GNU Scientific Library. For compiling the Python module, you need to have a recent version of swig installed. After running make in the top-level SZpack directory, you'll need to run it in the python subdirectory, which is the location of the SZpack module. You may have to include this location in the PYTHONPATH environment variable.

**NOTE**: Currently, use of the SZpack library to create S-Z projections in yt is limited to Python 2.x.

## Creating S-Z Projections¶

Once you have SZpack installed, making S-Z projections from yt datasets is fairly straightforward:

In [1]:
%matplotlib inline
import yt
from yt.analysis_modules.sunyaev_zeldovich.api import SZProjection

freqs = [90.,180.,240.]
szprj = SZProjection(ds, freqs)
<ipython-input-1-90c2c756cdc6>:3: VisibleDeprecationWarning: Development of the SZProjection module has been moved to the yt_astro_analysis package. This version is deprecated and will be removed from yt in a future release. See https://github.com/yt-project/yt_astro_analysis for further information.
from yt.analysis_modules.sunyaev_zeldovich.api import SZProjection

freqs is a list or array of frequencies in GHz at which the signal is to be computed. The SZProjection constructor also accepts the optional keywords, mue (mean molecular weight for computing the electron number density, 1.143 is the default) and high_order (set to True to compute terms in the S-Z signal expansion up to second-order in $T_{e,SZ}$ and $\beta$).

Once you have created the SZProjection object, you can use it to make on-axis and off-axis projections:

In [2]:
# An on-axis projection along the z-axis with width 10 Mpc, centered on the gas density maximum
szprj.on_axis("z", center="max", width=(10.0, "Mpc"), nx=400)
/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py:167: UserWarning: Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py:45: UserWarning: Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
ds.add_field(("gas", "t_squared"), function = _t_squared,
/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py:50: UserWarning: Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
ds.add_field(("gas","beta_par_squared"), function = _beta_par_squared,
/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py:55: UserWarning: Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
ds.add_field(("gas","beta_perp_squared"), function = _beta_perp_squared,
/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py:60: UserWarning: Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
ds.add_field(("gas","t_beta_par"), function = _t_beta_par,
/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py:65: UserWarning: Because 'sampling_type' not specified, yt will assume a cell 'sampling_type'
ds.add_field(("gas","t_sz"), function = _t_sz,
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-a68982c006a8> in <module>
1 # An on-axis projection along the z-axis with width 10 Mpc, centered on the gas density maximum
----> 2 szprj.on_axis("z", center="max", width=(10.0, "Mpc"), nx=400)

/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py in on_axis(self, axis, center, width, nx, source)
188         self.nx = nx
189
--> 190         self._compute_intensity(np.array(tau), np.array(Te), np.array(bpar),
191                                 np.array(omega1), np.array(sigma1),
192                                 np.array(kappa1), np.array(bperp2))

/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py in _compute_intensity(self, tau, Te, bpar, omega1, sigma1, kappa1, bperp2)
331             for j in range(ny):
332                 xo[:] = self.xinit[:]
--> 333                 SZpack.compute_combo_means(xo, tau[i,j], Te[i,j],
334                                            bpar[i,j], omega1[i,j],
335                                            sigma1[i,j], kappa1[i,j], bperp2[i,j])

NameError: name 'SZpack' is not defined

To make an off-axis projection, szprj.off_axis is called in the same way, except that the first argument is a three-component normal vector.

Currently, only one projection can be in memory at once. These methods create images of the projected S-Z signal at each requested frequency, which can be accessed dict-like from the projection object (e.g., szprj["90_GHz"]). Projections of other quantities may also be accessed; to see what fields are available call szprj.keys(). The methods also accept standard yt keywords for projections such as center, width, and source. The image buffer size can be controlled by setting nx.

## Writing out the S-Z Projections¶

You may want to output the S-Z images to figures suitable for inclusion in a paper, or save them to disk for later use. There are a few methods included for this purpose. For PNG figures with a colorbar and axes, use write_png:

In [3]:
szprj.write_png("SZ_example")

For simple output of the image data to disk, call write_hdf5:

In [4]:
szprj.write_hdf5("SZ_example.h5")

Finally, for output to FITS files which can be opened or analyzed using other programs (such as ds9), call export_fits.

In [5]:
szprj.write_fits("SZ_example.fits", clobber=True)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-5-5f51e413cf3b> in <module>
----> 1 szprj.write_fits("SZ_example.fits", clobber=True)

/tmp/yt/yt/utilities/parallel_tools/parallel_analysis_interface.py in root_only(*args, **kwargs)
317     def root_only(*args, **kwargs):
318         if not parallel_capable:
--> 319             return func(*args, **kwargs)
320         comm = _get_comm(args)
321         rv = None

/tmp/yt/yt/analysis_modules/sunyaev_zeldovich/projection.py in write_fits(self, filename, sky_scale, sky_center, overwrite, **kwargs)
392         w.wcs.ctype = ["LINEAR"]*2
393
--> 394         fib = FITSImageData(self.data, fields=self.data.keys(), wcs=w)
395         if sky_scale is not None and sky_center is not None:
396             fib.create_sky_wcs(sky_center, sky_scale)

/tmp/yt/yt/visualization/fits_image.py in __init__(self, data, fields, units, width, wcs)
173                 self.hdulist.append(hdu)
174
--> 175         self.shape = self.hdulist[0].shape
176         self.dimensionality = len(self.shape)
177

/tmp/yt/venv/lib/python3.8/site-packages/astropy/io/fits/hdu/hdulist.py in __getitem__(self, key)
302             # Raise a more helpful IndexError if the file was not fully read.
--> 304                 raise e
305             else:
306                 raise IndexError('HDU not found, possibly because the index '

/tmp/yt/venv/lib/python3.8/site-packages/astropy/io/fits/hdu/hdulist.py in __getitem__(self, key)
298         try:
--> 299             return self._try_while_unread_hdus(super().__getitem__,
300                                                self._positive_index_of(key))
301         except IndexError as e:

/tmp/yt/venv/lib/python3.8/site-packages/astropy/io/fits/hdu/hdulist.py in _try_while_unread_hdus(self, func, *args, **kwargs)
1094         while True:
1095             try:
-> 1096                 return func(*args, **kwargs)
1097             except Exception: