The Python Sky Model 3 software

The Python Sky Model (PySM) is a Python package used by Cosmic Microwave Background (CMB) experiments to simulate maps, in HEALPix pixelization, of the various diffuse astrophysical components of Galactic emission relevant at CMB frequencies (i.e. dust, synchrotron, free-free and Anomalous Microwave Emission), as well as the CMB itself. These maps may be integrated over a given instrument bandpass and smoothed with a given instrument beam. PySM 2, released in 2016, has become the de-facto standard for simulating Galactic emission, for example it is used by CMB-S4, Simons Observatory, LiteBird, PICO, CLASS, POLARBEAR and other CMB experiments, as shown by the 80+ citations of the PySM 2 publication. As the resolution of upcoming experiments increases, the PySM 2 software has started to show some limitations, the solution to these issues was to reimplement PySM from scratch focusing on these features: reimplement all the models with the numba Just-In-Time compiler for Python to reduce memory overhead and optimize performance; use MPI through mpi4py to coordinate execution of PySM 3 across multiple nodes and rely on libsharp, for distributed spherical harmonic transforms; employ the data utilities infrastructure provided by astropy to download the input templates and cache them when requested. At this stage we strive to maintain full compatibility with PySM 2, therefore we implement the exact same astrophysical emission models with the same naming scheme. In the extensive test suite we compare the output of each PySM 3 model with the results obtained by PySM 2.


Statement of Need
The Cosmic Microwave Background (CMB) radiation, emitted just 370 thousand years after the Big Bang, is a pristine probe of the Early Universe. After being emitted at high temperatures, the CMB was redshifted by the subsequent 13.8 billion years of cosmic expansion, such that it is brightest at microwave frequencies today. However, our own Milky Way galaxy also emits in the microwave portion of the spectrum, obscuring our view of the CMB. Examples of this emission are thermal radiation by interstellar dust grains, and synchrotron emission by relativistic electrons spiraling in magnetic fields. Cosmologists need to create synthetic maps of the CMB and of the galactic emission based on available data and on physical models that extrapolate observations to different frequencies. The resulting maps are useful to test data reduction algorithms, to understand residual systematics, to forecast maps produced by future instruments, to run Monte Carlo analysis for noise estimation, and more.

Summary
The Python Sky Model (PySM) is a Python package used by Cosmic Microwave Background (CMB) experiments to simulate maps, in HEALPix (Górski et al., 2005;Zonca et al., 2019) pixelization, of the various diffuse astrophysical components of Galactic emission relevant at CMB frequencies (i.e. dust, synchrotron, free-free and Anomalous Microwave Emission), as well as the CMB itself. These maps may be integrated over a given instrument bandpass and smoothed with a given instrument beam. The template emission maps used by PySM are based on Planck (Planck Collaboration, 2018) and WMAP (Bennett et al., 2013) data and are noise-dominated at small scales. Therefore, PySM simulation templates are smoothed to retain the large-scale information, and then supplemented with modulated Gaussian realizations at smaller scales. This strategy allows one to simulate data at higher resolution than the input maps.
PySM 2 (Thorne et al., 2017), released in 2016, has become the de-facto standard for simulating Galactic emission, for example it is used by CMB-S4, Simons Observatory, LiteBird, PICO, CLASS, POLARBEAR and other CMB experiments, as shown by the 80+ citations of the PySM 2 publication. As the resolution of upcoming experiments increases, the PySM 2 software has started to show some limitations: • Emission templates are provided at 7.9 arcminutes resolution (HEALPix N side = 512), while the next generation of CMB experiments will require sub-arcminute resolution. • The software is implemented in pure numpy, meaning that it has significant memory overhead and is not multi-threaded, precluding simply replacing the current templates with higher-resolution versions • Emission templates are included in the PySM 2 Python package, this is still practical when each of the roughly 40 input maps is~10 Megabytes, but will not be if they are over 1 Gigabyte.
The solution to these issues was to reimplement PySM from scratch focusing on these features: • Reimplement all the models with the numba (Lam et al., 2015) Just-In-Time compiler for Python to reduce memory overhead and optimize performance: the whole integration loop of a template map over the frequency response of an instrument is performed in a single pass in automatically compiled and multi-threaded Python code. • Use MPI through mpi4py to coordinate execution of PySM 3 across multiple nodes, this allows to support template maps at a resolution up to 0.4 arcminutes (HEALPix N side = 8192). • Rely on libsharp (Reinecke & Seljebotn, 2013), a distributed implementation of spherical harmonic transforms, to smooth the maps with the instrument beam when maps are distributed over multiple nodes with MPI. • Employ the data utilities infrastructure provided by astropy (Astropy Collaboration et al., 2018 to download the input templates and cache them when requested.
At this stage we strive to maintain full compatibility with PySM 2, therefore we implement the exact same astrophysical emission models with the same naming scheme. In the extensive test suite we compare the output of each PySM 3 model with the results obtained by PySM 2.

Performance
As an example of the performance improvements achieved with PySM 3 over PySM 2, we run the following configuration: • An instrument with 3 channels, with different beams, and a top-hat bandpass defined numerically at 10 frequency samples. • A sky model with the simplest models of dust, synchrotron, free-free and AME [a1,d1, s1,f1 in PySM terms]. • Execute on a 12-core Intel processor with 12 GB of RAM.
The following tables shows the walltime and peak memory usage of this simulation executed at the native PySM 2 resolution of N side = 512 and at two higher resolutions: At the moment it is not very useful to run at resolutions higher than N side = 512 because there is no actual template signal at smaller scales. However, it demonstrates the performance improvements that will make working with higher resolution templates possible.

Future work
PySM 3 opens the way to implement a new category of models at much higher resolution. However, instead of just upgrading the current models to smaller scales we want to also update them with the latest knowledge of Galactic emission and gather feedback from each of the numerous CMB experiments. For this reason we are collaborating with the Panexperiment Galactic Science group to lead the development of the new class of models to be included in PySM 3.

How to cite
If you are using PySM 3 for your work, please cite this paper for the software itself; for the actual emission modeling please also cite the original PySM 2 paper (Thorne et al., 2017). There will be a future paper on the generation of new PySM 3 astrophysical models.