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)


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)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
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)

/var/jenkins_home/workspace/yt_docs/sandbox/temp/lib64/python3.4/site-packages/yt-3.3.5-py3.4-linux-x86_64.egg/yt/analysis_modules/sunyaev_zeldovich/projection.py in on_axis(self, axis, center, width, nx, source)
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))
193
194         self.ds.field_info.pop(("gas","beta_par"))

/var/jenkins_home/workspace/yt_docs/sandbox/temp/lib64/python3.4/site-packages/yt-3.3.5-py3.4-linux-x86_64.egg/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")

/usr/lib64/python3.4/site-packages/matplotlib/__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

warnings.warn(_use_error_msg)


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)

/usr/lib64/python3.4/site-packages/IPython/kernel/__init__.py:13: ShimWarning: The IPython.kernel package has been deprecated. You should import from ipykernel or jupyter_client instead.
"You should import from ipykernel or jupyter_client instead.", ShimWarning)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-5-877886b6eeef> in <module>()
----> 1 szprj.write_fits("SZ_example.fits", clobber=True)

/var/jenkins_home/workspace/yt_docs/sandbox/temp/lib64/python3.4/site-packages/yt-3.3.5-py3.4-linux-x86_64.egg/yt/utilities/parallel_tools/parallel_analysis_interface.py in root_only(*args, **kwargs)
318     def root_only(*args, **kwargs):
319         if not parallel_capable:
--> 320             return func(*args, **kwargs)
321         comm = _get_comm(args)
322         rv = None

/var/jenkins_home/workspace/yt_docs/sandbox/temp/lib64/python3.4/site-packages/yt-3.3.5-py3.4-linux-x86_64.egg/yt/analysis_modules/sunyaev_zeldovich/projection.py in write_fits(self, filename, sky_scale, sky_center, clobber)
388         w.wcs.ctype = ["LINEAR"]*2
389
--> 390         fib = FITSImageData(self.data, fields=self.data.keys(), wcs=w)
391         if sky_scale is not None and sky_center is not None:
392             fib.create_sky_wcs(sky_center, sky_scale)

/var/jenkins_home/workspace/yt_docs/sandbox/temp/lib64/python3.4/site-packages/yt-3.3.5-py3.4-linux-x86_64.egg/yt/utilities/fits_image.py in __init__(self, data, fields, units, width, wcs)
130                 self.hdulist.append(hdu)
131
--> 132         self.shape = self.hdulist[0].shape
133         self.dimensionality = len(self.shape)
134

/usr/lib64/python3.4/site-packages/astropy/io/fits/hdu/hdulist.py in __getitem__(self, key)
175
176         idx = self.index_of(key)
--> 177         return super(HDUList, self).__getitem__(idx)
178
179     def __contains__(self, item):

IndexError: list index out of range

which would write all of the projections to a single FITS file, including coordinate information in kpc. The optional keyword clobber allows a previous file to be overwritten.