so_pysm_models Documentation

This is the documentation for so_pysm_models.

Getting started

Installation

Requirements:

Clone the repository:

pip install https://github.com/simonsobs/so_pysm_models/archive/master.zip

Development installation

Clone from Github and install:

git clone https://github.com/simonsobs/so_pysm_models
cd so_pysm_models
pip install -e .

Run unit tests:

python setup.py test -V

Build docs:

python setup.py build_docs -w

Example Usage

This repository implements new models for PySM that can be added as additional components.

For example, create and configure a component:

from so_pysm_models import GaussianSynchrotron
synchrotron = GaussianSynchrotron(target_nside = 16)

Create a PySM sky and add this component:

sky = pysm.Sky({})
sky.add_component("gaussian_synch", gaussian_synch)

Then get a map at a specific frequency in GHz with standard PySM functionalities:

m_synch = sky.gaussian_synch(2.3)

see example notebooks:

Summary of Models

This page contains high-level documentation about the available models, check the classes docstrings, or the online documentation, for the specific arguments.

The input template maps for many models are available at NERSC on the cmb project space at:

/global/project/projectdirs/cmb/www/so_pysm_models_data

they are also published via web at http://portal.nersc.gov/project/cmb/so_pysm_models_data/.

GaussianSynchrotron

This class implements Gaussian simulations for Galactic synchrotron emission. The inputs are a bunch of parameters defining the properties of the synchrotron power spectra, and of synchrotron Spectral Energy Distribution (SED), the output are the stokes IQU maps simulated as Gaussian random fields of the defined spectra. In particular, synchrotron power spectra \(C_{\ell}\) are assumed to follow a power law as a function of \(\ell\): \(C_{\ell}^{TT/TE/EE/BB}\propto\ell^{\alpha}\). Spectra are defined by:

  1. The slope \(\alpha\) (same for all the spectra)
  2. The amplitude of TT and EE spectra at \(\ell=80\),
  3. The ratio between B and E-modes

Stokes Q and U maps are generated as random realization of the polarization spectra. For the temperature map the situation is slightly different as we want the total intensity map to be positive everywhere. The Stokes I map is generated in the following way:

if target Nside<=64:
  1. The TT power spectrum is \(C_\ell \propto \ell^\alpha\) and \(C_\ell[0]=0\)
  2. A first temparature map T is generated as a gaussian realization of this power spectrum
  3. A new map is obtained by adding to T an offset whose value is taken from a reference map
  4. If T+offset is positive everywhere than this is the output temperature map
  5. Otherwise a cut in the TT power spectrum is applied in the following way: \(C_\ell[1:\ell_{cut}] = C_\ell[\ell_{cut}]\)
  6. A new T+offset map is generated. The value of \(\ell_{cut}\) is the minimum one for which T+offset is positive everywhere
if target Nside>64:
  1. a map at Nside=64 is generated following the procedure above and then filtered to retain only large angular scales (ell<30)
  2. a map at the target Nside is generated including only small scales (ell>30) with the same seed as the map at point 1.
  3. the two maps are added together
  4. In case the coadded map still has negative pixels a small offset is added to make it positive everywhere

The default parameters are optimized for SO-SAT observations. Meaning that the amplitudes of power spectra are normalized in the 10% sky region observed by the instrument. In particular:

1. The amplitude of TT spectrum is taken from PySM-s0 model at 23GHz. TT_amplitude = 20 \(\mu K^2\) (for \(D_\ell\) at \(\ell=80\)) 1. The offset for T map is also taken from PySM-s0 model at 23GHz. Toffset = 72 \(\mu K\) 1. The amplitude of EE spectrum is taken from S-PASS at 2.3GHz extrapolated at 23GHz with a powerlaw with \(\beta_s=-3.1\) EE_amplitude = 4.3 math:mu K^2 (for \(D_\ell\) at \(\ell=80\)) 1. ratio between B and E modes from Krachmalnicoff et al. 2018, B_to_E = 0.5 1. spectral tilt from Krachmalnicoff et al 2018, alpha = -1 1. spectral index from Planck IX 2018, beta = -3.1 1. Default value for curvature is zero

GaussianDust

This class implements Gaussian simulations for Galactic thermal dust emission. The inputs are a bunch of parameters defining the properties of dust power spectra, and of dust Spectral Energy Distribution (SED), the output are the stokes IQU maps simulated as Gaussian random fields of the defined spectra. In particular, dust power spectra \(C_{\ell}\) are assumed to follow a power law as a function of \(\ell\): \(C_{\ell}^{TT/TE/EE/BB}\propto\ell^{\alpha}\). Spectra are defined by:

  1. The slope \(\alpha\) (same for all the spectra)
  2. The amplitude of TT and EE spectra at \(\ell=80\),
  3. The ratio between B and E-modes
  4. The degree of correlation between T and E.

Stokes Q and U maps are generated as random realization of the polarization spectra. For the temperature map the situation is slightly different as we want the total intensity map to be positive everywhere. The Stokes I map is generated in the following way:

if target Nside<=64:
  1. The TT power spectrum is \(C_\ell \propto \ell^\alpha\) and \(C_\ell[0]=0\)
  2. A first temparature map T is generated as a gaussian realization of this power spectrum
  3. A new map is obtained by adding to T an offset whose value is taken from a reference map
  4. If T+offset is positive everywhere than this is the output temperature map
  5. Otherwise a cut in the TT power spectrum is applied in the following way: \(C_\ell[1:\ell_{cut}] = C_\ell[\ell_{cut}]\)
  6. A new T+offset map is generated. The value of \(\ell_{cut}\) is the minimum one for which T+offset is positive everywhere
if target Nside>64:
  1. a map at Nside=64 is generated following the procedure above and then filtered to retain only large angular scales (ell<30)
  2. a map at the target Nside is generated including only small scales (ell>30) with the same seed as the map at point 1.
  3. the two maps are added together
  4. In case the coadded map still has negative pixels a small offset is added to make it positive everywhere

Typical values for \(\ell_{cut}\) are between \(\ell=4\) and \(\ell=9\), depending on realization (and also on the Nside of the output map). This implementation removes some power at the very large scales.

The default parameters are optimized for SO-SAT observations. Meaning that the amplitudes of power spectra are normalized in the 10% sky region observed by the instrument. In particular:

  1. The amplitude of TT spectrum is taken from PySM-d0 model at 353GHz. TT_amplitude = 350 \(\mu K^2\) (for \(D_\ell\) at \(\ell=80\))
  2. The offset for T map is also taken from PySM-d0 model at 353GHz. Toffset = 18 \(\mu K\)
  3. The amplitude of EE spectrum is taken from Planck map at 353GHz, EE_amplitude = 100 math:mu K^2 (for \(D_\ell\) at \(\ell=80\))
  4. ratio between B and E modes from Planck IX 2018, B_to_E = 0.5
  5. spectral tilt from Planck IX 2018, alpha = -0.42
  6. spectral index and temperature from Planck IX 2018, beta = 1.53, T=19.6 K

COLines

COLines is not a standard PySM component because PySM does not allow to distinguish between a case where a component is evaluated for the purpose of integrating over the bandpass or evaluated for separate channels. Therefore this class should be instantiated choosing the desired line and summed to the output of PySM. For example:

from so_pysm_models import COLines
co = COLines(target_nside=16, output_units="uK_CMB", line="10")
pysm_map += bandpass_weight * hp.smoothing(co.signal(), fwhm=fwhm)

Where bandpass_weight is the scalar transmission at the line frequency (which is available at co.line_frequency), i.e. if the bandpass is a top-hat between 110 and 120 GHz, the “10” line emission should be multiplied by 0.1.

This class implements simulations for Galactic CO emission involving the first 3 CO rotational lines, i.e. \(J=1-0,2-1,3-2\) whose center frequency is respectively at \(\nu_0 = 115.3, 230.5,345.8\) GHz. The CO emission map templates are the CO Planck maps obtained with MILCA component separation algorithm (See Planck paper). The CO maps have been released at the nominal resolution (10 and 5 arcminutes). However, to reduce noise contamination from template maps (especially at intermediate and high Galactic latitudes), we convolved them with a 1 deg gaussian beam.

The Stokes I map is computed from the template one as it follows:

if target Nside <= 512:

  1. The template map at a nside=512 is downgraded at the target nside

if target Nside > 512 :

  1. The template map at a nside=2048 is downgraded(eventually upgraded) at the target nside

Q and U maps can be computed from the template CO emission map, \(I_{CO}\), assuming a constant fractional polarization, as:

\[ \begin{align}\begin{aligned}Q = f_{pol} I_{CO} g_d \cos( 2 \psi)\\U = f_{pol} I_{CO} g_d \sin( 2 \psi)\end{aligned}\end{align} \]

with \(g_d\) and \(\psi\) being respectively the depolarization and polarization angle maps estimated from a dust map as :

\[ \begin{align}\begin{aligned}g_d = \frac{ \sqrt{Q^2_{d,353} + U^2_{d,353} } }{f_{pol} I_{d,353} }\\\psi = \frac{1}{2} \arctan {\frac{U_{d,353}}{Q_{d,353}}}\end{aligned}\end{align} \]

Most of the CO emission is expected to be confined in the Galactic midplane. However, there are still regions at high Galactic latitudes where the CO emission has been purely assessed (by current surveys) and where the Planck signal-to-noise was not enough to detect any emission.

The PySM user can include the eventuality of molecular emission (both unpolarized and polarized) at High Gal. Latitudes by coadding to the emission maps one realization of CO emission simulated with MCMole3D together with the Planck CO map. The polarization is simulated similarly as above.

The MCMole3D input parameters are are obtained from best fit with the Planck CO 1-0 map (see Puglisi et al. 2017 and the documentation). If include_high_galactic_latitude_clouds=True, a mock CO cloud map is simulated with MCMole3D, encoding high Galactic latitudes clouds at latitudes above and below than 20 degres. The mock emission map is then coadded to the Planck CO emission map. The polarization is simulated similarly as above.

The installation of mcmole3d is not required, HGL clouds can be input to the CO emission by setting run_mcmole3d=False (which is the default). However, if one wants to run several mock CO realizations observing high Galactic latitude patches we encourage to run mcmole3d by changing random_seed in the CO class constructor. The parameter theta_high_galactic_latitude_deg set the latitude above which CO emission from high Galactic latitudes can be included and it has an impact only when run_mcmole3d=True.

The default parameters are set to include CO 1-0 emission and polarization (with 0.1% constant polarization fraction), in particular:

  1. polarization_fraction= 0.001, on average is the expected level on 10% regions of the sky. However, polarization from CO emission have been detected at larger fluxes in Orion and Taurus complexes (Greaves et al.1999 )
  2. theta_high_galactic_latitude_deg = 20, includes CO emission at \(|b|>\theta_{hgl}\) from one realization of mcmole3d maps. Be aware that the larger \(theta_{hgl}\), the farther is the Galactic plane and the more unlikely is to find high Galactic latitude clouds.

PrecomputedAlms

This class generates a PySM component based on a set of precomputed \(a_{\ell,m}\) coefficients stored in a folder in FITS format. This is mostly targeted at simulations of the Cosmic Microwave Background, the input \(a_{\ell,m}\) can be in K_{RJ} or K_{CMB} as defined in the constructor, the unit conversion is performed assuming the CMB black body spectrum. The output unit is specified in the signal method, default is mu K_{RJ}, as expected by PySM. In case the input is in K_{RJ}, it is necessary also to specify input_reference_frequency_GHz.

The transformation between Spherical Harmonics and pixel domain can be performed either during initialization or in the signal method based on precompute_output_map.

See the documentation about mapsims about specific simulated datasets.

InterpolatingComponent

Adds a custom emission to the sky simulated by PySM defined as a set of template maps at pre-defined frequencies to be interpolated at the frequencies requested through PySM.

Inputs

A folder of maps named with their frequency in GHz with the flux in any unit supported by PySM (e.g. Jysr, MJsr, uK_RJ, K_CMB). They don’t need to be equally spaced

For example:

ls `cib_precomputed_maps/`
0010.0.fits 0015.0.fits 0018.0.fits

Usage

Instantiate InterpolatingComponent and point it to the folder, define the unit and the target nside (same used by PySM). It supports all interpolation_kind of scipy.interpolate.interp1d(), e.g. “nearest”, “linear”, “quadratic”, “cubic”:

cib = InterpolatingComponent(path="cib_precomputed_maps", input_units="MJysr", target_nside=nside, interpolation_kind="linear",
                         has_polarization=False, verbose=True)

Full example notebook

High resolution templates

so_pysm_models also provides access to templates with higher resolution and with updated data compared to the models included in PySM.

They can be accessed with the function get_so_models() which works similarly to the models function available in PySM, and you can mix them, for example:

from so_pysm_models import get_so_models
from pysm.nominal import models
from pysm import Sky
sky = Sky({
        "dust" : get_so_models("SO_d0", nside=128),
        "synchrotron" : models("s1", nside=128)
})

so_pysm_models retrieves the templates when needed from NERSC via web accessing: http://portal.nersc.gov/project/cmb/so_pysm_models_data/ Downloaded files are stored in the astropy cache, generally astropy/cache and are accessible using astropy.utils.data, e.g. astropy.utils.data.get_cached_urls() gives the list of downloaded files. If running at NERSC, the module automatically uses the files accessible locally from the /project filesystem.

Low-resolution templates are standard PySM ones at \(N_{side}\) 512, often with updated parameters based on Planck results. High-resolution templates are computed from the low-resolution ones, by extrapolating power spectra considering a simple power law model, and by generating small scales as Gaussian realization of these spectra. High-resolution templates therefore have Gaussian small scales (for \(\ell > ~ 1000\)) modulated with large scale signal for both temperature and polarization.

You can access the high resolution parameters at \(N_{side}\) 4096 appending s (for small scale) at the end of each model name, for example:

from so_pysm_models import get_so_models
from pysm import Sky
sky_highres = Sky({
        "dust" : get_so_models("SO_d0s", nside=4096),
        "synchrotron" : get_so_models("SO_s0s", nside=4096)
})

Whatever the \(N_{side}\) of the input model and the requested \(N_{side}\) in get_so_models(), PySM will automatically use healpy.ud_grade() to adjust the map resolution.

Details about individual models

Append “s” after a model name to access the \(N_{side}\) 4096 template, i.e. SO_f0s.

Dust

  • SO_d0: Thermal dust is modeled as a single-component modified black body, with same templates as in PySM model d1. There is no spatial variation of temperature and emissivity in the sky: \(T=19.6\) K and \(\beta_d=1.53\) (values taken from Planck Collaboration IX 2018).

Synchrotron

  • SO_s0: Templates from PySM model s1. Power law spectral energy distribution, with fixed spectral index \(\beta_s=-3.1\) (from Planck Collaboration IX 2018).

Free Free

  • SO_f0: same model as PySM f1, no spatial variation of spectral index equal to -2.4.

AME

  • SO_a0: sum of two spinning dust populations (as in PySM model a1) with spatially constant peak frequency. No polarization.

Reference/API

so_pysm_models Package

Functions

get_so_models(key, nside[, pixel_indices, …])
test(**kwargs) Run the tests for the package.

Classes

COLines(target_nside, output_units[, …]) Class defining attributes for CO line emission.
GaussianDust(target_nside[, …]) Gaussian dust model
GaussianSynchrotron(target_nside[, …]) Gaussian synchrotron model
InterpolatingComponent(path, input_units, …) PySM component interpolating between precomputed maps
PrecomputedAlms(filename[, input_units, …]) Generic component based on Precomputed Alms
UnsupportedPythonError

Class Inheritance Diagram

Inheritance diagram of so_pysm_models.co_lines.COLines, so_pysm_models.dust.GaussianDust, so_pysm_models.synchrotron.GaussianSynchrotron, so_pysm_models.interpolating.InterpolatingComponent, so_pysm_models.alms.PrecomputedAlms, so_pysm_models.UnsupportedPythonError