This is the fourth and final report of the course Data Science Laboratory. Here I’m summarizing everything I’ve done during this course in detail and give some outlook on what other tasks can be done regarding this subject as a project work.
During the course I’ve worked on a project, where I studied the Cosmic Microwave Background radiation. This phenomenon is an actively researched topic in astronomy, because it can provide us an enormous amount of information about cosmological parameters and the early state of the universe. By measuring the angular power spectrum of the CMB temperature anisotropy, we can tune parameters of the currently most widely accepted (standard) model of cosmology, the $\Lambda$CDM model. Measuring the CMB radiation can give us insight about numerous astrophysical objects too, as they leave traces and “shadows” on its temperature map. These objects can be eg. distant and dense galaxy clusters, Active Galactic Nuclei (AGNs), Dust Star Forming Galaxies (DSFGs), etc.
In my second report I gave a short overview on the history of the CMB research. Went into detail how the scientific consensus established slowly around the idea of the CMB being the remnant of the Big Bang and how it inspired the technological advancements in space exploration. In the past two and the half decades, three notable space observatories (COBE[1], WMAP[2] and Planck[3]) were launched to gather unprecedented amount data about the CMB radiation and map the temperature anisotropy of the whole sky. The most recent one was the Planck space telescope, which provided us the most detailed and highest resolution datasets of the CMB radiation. In my project I also used these datasets.
In the third report I’ve already summarized in detail the programming aspects and progress of this project work. This time I will give a general overlook on the subject, talk about the tasks I’ve worked on, and outline some inspirational ideas about what other questions can be examined and researched in the context of this project. The main goal for me was to learn some insight about how to simulate and how to analyze CMB temperature anisotropy maps. In the first case I had to research and learn about the methods, which are usually used to simulate realistic images of these temperature maps and finally implement them. Similarly, in the second task I’ve been ought to process real observational data and I had to look for and learn the methods that were needed to study true CMB temperature maps and then analyze them.
I started off the work by reading through the initial supplementary materials, which included a gentle introduction to the theoretical basics of the CMB radiation[4], the presentation slides from a lecture about a short introduction to the CMB radiation, held on the University of Oslo[5] and finally the most recent paper about the research of the Planck space observatory. An extra source given were the presentation materials from the event called CMB Summer School 2019, which was organized by The McMahon Cosmology Laboratory[6] in the form of Jupyter Notebooks, written and organized by Jeff McMahon, Renée Hložek, Mat Madhavacheril, Sigurd Naess, and Alex Van Engelen. While the first three materials helped me to understand the theoretical basics of the CMB radiation, the latter one served as a perfect base to learn the methods of creating CMB simulations. It also gave an overview of the different astrophysical and cosmological components, which form the layers of a real image depicting the CMB radiation.
During the semester I’ve written two completely different pipelines for the simulation of CMB temperature anisotropy maps and one to load, process and analyze - true or generated - CMB data, which was encoded using the HEALPix spherical indexing standard[7]. One of the simulation methods was based on the materials from the mentioned CMB Summer School, while the other method utilized the built-in subroutines of the HEALPix library.
The codes were written in Python 3.9.0, using the latest release of all the utilized Python packages. The finished material features some very small Python “libraries” alongside Jupyter Notebooks for demonstration purposes, which showcases the usage of these libraries. All other online or other tool I’ve used are mentioned throughout this report at the appropriate places.
It is a well-known fact for science today that the Cosmic Microwave Background radiation is the remnant of the Big Bang in the form of heavily redshifted photons. The intensity of this black-body radiation is uniform throughout all directions on the sky and as our most recent measurements indicates, its average colour temperature is approximately $2.725$ K.
According to our current cosmological knowledge, the photons consisting the CMB radiation came to existence at the time of the Big Bang. The early universe was very hot and dense for hundreds of thousands years long since the Big Bang. Because of its excessive temperature, electrons existed only in an unbounded, free gas state and completely filled the early universe. This caused the existing photons to constantly scatter on particles in this electron gas and thus their mean free path was immensely short. Because of the constant scattering, photons and electron were also in thermal equilibrium. Since photons couldn’t travel to long distances, the universe was opaque to electromagnetic radiation at this point. As the universe expanded it slowly became more spares and colder. Approximately after $370\,000$ years, the universe cooled down enough (to $\approx 3000$ K), where the formation of hydrogen atoms was energetically favorable.
This period in the history of the universe is called as the “recombination era” or “decoupling”. The electrons combined (“recombined”) with the existing hydrogen nuclei and thus the photons “decoupled” from the matter, they stopped being in thermal equilibrium and was able to travel freely at virtually infinite distances without being constantly scattered and thus being “binded/coupled to matter”. The term “recombined” refers to the now-defeated cosmological hypothesis that this event had already happened before, maybe numerous times. However physics still kept this original name to describe this era because of historical reasons, but now it’s more appropriate to refer to this period as the “decoupling”.
At this point the dense electron cloud filling the universe cleared up and photons scattered for the last time on free electrons. These photons are those, which we observe today as the CMB radiation itself. At the time of the last scattering the photons was the same temperature as the dissipated electron gas because of being on thermal equilibrium, but after that the CMB photons traveled practically without anything to cool them down. However the universe constantly expanded since then, and now because of this stretching, the CMB photons redshifted from their initial $3000$ K temperature the the currently measured $2.725$ K magnitude, which is equivalent to a redshift of $z \approx 1200$.
As it was already mentioned, the CMB radiation is completely uniform in intensity in all directions on the sky. Because of this nature, the CMB photons may be considered to coming towards us from a very distant surface - where we are in the center -, which is usually referred to as the “Last Scattering Surface” (LSS). Photons travel at a finite speed, which means as time goes forward, this surface constantly expanding. First those photons reach our telescopes, which are the closest to us. But after some time $t$, farther photons will now had enough time to travel the distance between their positions at $t=0$ (the period of decoupling) and us. As time goes on, CMB photons reach us from farther and farther distances. This distance is - according to our best knowledge about the structure of the universe - is equal in every direction. It thus defines a spherical surface at any time $t>0$, which we can consider to be the origin of the currently observed CMB radiation, or as it called, the LSS.
The mathematical description of the CMB must have to build on this property. The general idea is to describe the observable quantities as a function over a spherical surface. Fortunately this problem arises in a lot of topics in physics and mathematics and thus it was already described by harmonic analysis and by studying spherical harmonics way before humanity reached to the point of considering the existence of the CMB.
In my first report I mentioned that the CMB temperature anisotropy is just a function over a spherical surface (this is the LSS) and can be denoted as
\begin{equation} \label{eq:1} f \left( \vartheta, \varphi \right) \equiv T \left( \vartheta, \varphi \right). \end{equation}
Our goal is to measure this temperature function and interpret the results of our observations. We hope to extract some hidden or deeper information from this observed function, which gives us insight about other, directly unobservable quantities or parameters. To construct the necessary conditions for any analysis like this, we have to break down the function $T \left( \vartheta, \varphi \right)$ and formulate it in an appropriate way for this. A function over a spherical surface can be expanded into a spherical harmonic series as
\begin{equation} \label{eq:2} f \left( \vartheta, \varphi \right) = \sum_{\ell\,=\,0}^{\infty} \sum_{m\,=\,-\ell}^{\ell} a_{\ell m} Y_{\ell m} \left( \vartheta, \varphi \right), \end{equation} where $\vartheta$ and $\varphi$ denote the spherical coordinates with polar angle $\vartheta \in \left[ 0, \pi \right]$ and azimuth $\varphi \in \left[0, 2 \pi \right)$. The term $Y_{\ell m} \left( \vartheta, \varphi \right)$ denotes the well-known formula of spherical harmonics
\begin{equation} \label{eq:3} Y_{\ell m} \left( \vartheta, \varphi \right) = \sqrt{ \frac{2 \ell + 1}{4 \pi} \frac{\left( \ell - m \right)!}{\left( \ell + m \right)!} } P_{\ell m} \left( \cos \left( \vartheta \right) \right) e^{i m \varphi}. \end{equation} As seen from equation \eqref{eq:2}, the real information about an observed CMB temperature map $T \left( \vartheta, \varphi \right)$ will be encoded in the $a_{\ell m}$ coefficients, so we need to determine these values. By reformulating the equation \eqref{eq:2} and take into account the orthonormal properties of the $Y_{\ell m}$ spherical harmonics, we can express $a_{\ell m}$ for any $f \left( \vartheta, \varphi \right)$ as
\begin{equation} \label{eq:4} a_{\ell m} = \int_{0}^{2 \pi} \int_{0}^{\pi} Y_{\ell m}^{\ast} \left( \vartheta, \varphi \right) f \left( \vartheta, \varphi \right) \sin \left( \vartheta \right)\,\mathrm{d}\vartheta\,\mathrm{d}\varphi \end{equation}
As it was also mentioned in my previous reports, our measurements and observations in physics are not perfect and definitely not continuous. On the macroscopic level we can sample only discrete chunks of reality or making observation only at finite time steps, so we need to work with discrete data all the time, regardless of the measured quantity. However this couldn’t (and shouldn’t) stop us to approximate continuous quantities with discrete measurements.
In the case of the CMB radiation this discretization arises from the finite resolution and finite beam size of our instruments. This means we can’t measure on every frequency ranges, can’t resolve too small and too distant objects or separate two objects on an image that are too close to each other. Our final observations of the sky will be pixelated images, which means that our measured $T \left( \vartheta, \varphi \right)$ temperature function will be discrete. As a result we can and should reformulate equation \eqref{eq:4} for discrete values of $f \left( \vartheta, \varphi \right)$. That’s what I’ve already presented in my second report talking about the HEALPix standard.
HEALPix is on itself is an indexing convention and pixelization standard of the sphere and a library for efficient discrete spherical harmonic analysis. This framework utilizes the discrete form of the formula for $a_{\ell m}$ values, which is
\begin{equation} \label{eq:5} a_{\ell m} = \frac{4 \pi}{N_{\mathrm{pix}}} \sum_{p\,=\,0}^{N_{\mathrm{pix}} - 1} Y_{\ell m}^{\ast} \left( \theta_{p}, \varphi_{p} \right) f_{p} \left( \theta_{p}, \varphi_{p} \right), \end{equation} where $N_{\mathrm{pix}}$ is the number of pixels of the image and $p$ denotes the index of pixels. HEALPix helps us by standardizing the order of these pixels and offers us complete subroutines to work with data encoded in this format. Using the discretized $a_{\ell m}$ coefficients we can calculate an approximation of the angular power spectrum as
\begin{equation} \label{eq:6} C_{\ell} = \frac{1}{2 \ell + 1} \sum_{m\,=\,-\ell}^{\ell} \left| a_{\ell m} \right|^{2}. \end{equation}
This quantity derived from the measured $T \left( \vartheta, \varphi \right)$ map can be finally used for further research, as it can be utilized to fit cosmological parameters and thus can be used to make our cosmological models more precise. The available methods for this parameter search were not discussed in this project.
To maintain the consistency of my previous reports, I’m referring as “Generation method 1.” to the simulation method involving the use of the HEALPix library, while referring as “Generation method 2.” to the other one.
The HEALPix library has been already implemented in Python under the package name healpy. It contains all subroutines from its original Fortran90 implementation with almost very small to no changes[8]. This method waits for an arbitrary (true, generated or completely random) power spectrum as an input. Reversing the operation in equations \eqref{eq:5} and \eqref{eq:6}, we’re able to randomly generate an $a_{\ell m}$ and a corresponding $T \left( \vartheta, \varphi \right)$ 1D HEALPix array from the input $C_{\ell}$ dataset. The HEALPix routine synfast does all the work here and it generates completely random maps with $C_{\ell}$ angular power spectrum for every rerun.
Independently from the creation of a random HEALPix array, healpy offers us a routine that transforms an arbitrary input HEALPix dataset to any geographical projections that are implemented in the healpy package. This routine creates a rectangular matrix, which then can be simply visualized using eg. matplotlib or other graphic packages. This projector method healpy.projector creates rectangular matrices already in the correct, chosen projections. This way all non-rectangular maps will be encompassed between the “borders” of the matrix, where elements outside the projection will be marked as -inf values. I created a custom colormap to make the images created with it as similar as possible to the CMB maps of the Planck space telescope. All the color values for a seamless gradient was generated from an input file. The colormap was also configured to correctly display these special, -inf values, as well as the actual CMB map itself.
Some sample images can be seen on Figure 2. above. The synfast function is capable to take into account the effective beam size of the simulated instrument and create images accordingly. The comparison of two different beam sizes can be also seen on this figure.
The only “disadvantage” of using this method is that it’s only capable to generate full-sky images. The output HEALPix array will always cover the whole sky and it makes the routine a bit slower and memory-intensive. Also generating custom sized cut-outs or tiles of the CMB cannot be efficiently done using this algorithm.
The second generation method is completely different in every aspects. This method is the one, which was based on the presentation materials of the CMB Summer School of The McMahon Cosmology Laboratory. The original material consist of $10$ different Jupyter Notebooks containing really instructive, introductory materials to the CMB research, but they lack of organization, modern Python syntax and compliance with clear code conduct. I couldn’t use these codes to create such figures and results as I wanted to have at the end of the project, so I had to rebuild a new Python “library” (and visualization pipeline) on top of them using the methods detailed in a written and programmed form in these materials.
At the end of the semester I’ve created a small Python library along with three Jupyter Notebooks to demonstrate the usage and function of my completed routines. The base idea of this method is to create each different layers of the CMB temperature map individually. This simulation is limited however in its capabilities and creates more like just “mock images” of the CMB temperature anisotropy with close-to-correct power spectrum, than completely realistic CMB maps that include all of astrophysical and instrumental effects.
The naïve method contains routines for simulating the pure CMB temperature map, which serves as a base layer. It also contains routines to randomly generate different foreground effects, namely point sources and the Sunyaev-Zeldovic effect. At last but not least it is able to simulate the effects arising from the finite size of the instrumental beam, as well as different noise components.
The full simulation starts off by creating the pure CMB temperature map using an arbitrary angular power spectrum $C_{\ell}$. This steps begins by creating a noise map with Gaussian distribution in Fourier space, which parameters are based on the input power spectrum. After that the map is converted to image space and it results us a randomized temperature noise map with angular power spectrum close-to-identical to the original input spectrum. This can simply be considered to a randomly generated CMB temperature map, which serves as our base layer for the simulation.
A short summary was already given in my third report that what effects the foreground map intends to simulate. As I also noted, in reality the CMB radiation is completely obscured by the radiation of objects inside and outside of our galaxy. Proper filtering pipelines needs to be utilized to detach the foreground effects from the CMB. Because we’re speaking about billions of objects in reality, I’ve only simulated this effect in a much more modest way by only replicating the effect of just a couple thousands of them.
In this simulation three types of foreground objects are generated:
Point sources can be arise from a number of (mostly) bright astronomical objects, like Active Galactic Nuclei (AGNs), Dust Star Forming Galaxies (DSFGs), and even the bright tail of lensed DSFGs, as they are listed as examples in the original Jupyter Notebooks of Jeff McMahon and Renée Hložek. To simulate a mock image of these sources, we’re approximating their population by combining a layer of abundant, but faint sources with a layer of a few, but bright ones. This is what the two different distribution for point sources represent.
On the other hand the Sunyaev-Zeldovich “sources” are not refers to specific types of objects, but rather to any object, which causes the phenomenon itself. When CMB photons propagate through galaxy clusters filled with high-energy electrons, they receive a small energy boost by interacting with these electrons through inverse Compton scattering. This changes the photon temperatures and polarizations as well. Since we’re only working with temperature maps, only the temperature component of the SZ-effect is considered here. Speaking on technical terms, observing the Sunyaev-Zeldovich effect can be used to discover distant and dense galaxy clusters and thus is a field of great interest in cosmology to study.
The effects of finite beam sizes were already detailed above. Because of its finite size and FWHM, our instrument simply don’t have infinite resolution. While creating an image it can’t resolve too faint objects neither objects that are seen too close to each other on the sky. This effect results in a bit blurry image, as we can see it on Figure 1. The naive method however able to generate images where the effects of the beam size are much more trivial.
Here we assume that the beam profile is Gaussian. To simulate its effect on the measured temperature map, we convolve a Gaussian distribution with a given FWHM over the generated CMB image. The wider the beam, the more blurry our image will be.
In the framework of the naïve method, three different types of noise are simulated. Two of them occurs everywhere in CMB measurements, while one of them (atmospheric noise) only affects ground-based instruments. To generalize the CMB generation to these observatories too, it is mandatory to implement this noise component too. The origin of this noise is trivial, as it arises from the atmospheric fluctuations of Earth, which are well-known to every human being on this planet in the form of wind and other phenomenons.
White noise and $1/f$ noise (or pink noise) are literally everywhere. They both occurs in every measurement and affects every scientific observations. White noise refers generally to random noise and this one is also well-known to most human beings but maybe on another name. The main characteristics of white noise is to have a completely uniform power spectrum on every frequency ranges. In contrast, pink noise is probably much less famous, however it can be found in just as many natural processes or measurements as white noise does. It’s other name, “$1/f$ noise” stands for the behavior of its power spectrum density to be inversely proportional to the frequency. In the context of electronic devices its effect is called as “flicker noise” and can affect all measurements made with instruments with electronic parts.
A comparison of all the three generated noises can be seen on Figure 6.
After the generation of every layer in the model, they can be assembled by simply adding the generated images together. The beam effect should be applied always before the noise map added to the batch, either on a map with or without foreground effects.
In this part of the project my goal was to load a real, observed and processed CMB temperature map and try to extract the CMB’s power spectrum from the image. I used the most recent measurements of the CMB, made by the Planck space observatory between 2009 and 2013. Since these measurement currently are the most detailed, available datasets, the choice to use them was obvious. All the datasets made by the Planck telescope can be accessed at the website of the NASA/IPAC Infrared Science Archive[9]. For the project and for demonstration I used the output of one of Planck’s preprocessing pipelines, called “Commander”. Another main advantage of using the data from Planck is that they store every dataset encoded using the HEALPix standard, so processing them is quiet easy with healpy.
Using the mathematical conventions described in Section III.2, we only need two steps to reconstruct the power spectrum from an available CMB dataset. After reading in the $T \left( \vartheta, \varphi \right)$ temperature data using the HEALPix library, we have to determine the $a_{\ell m}$ coefficients and then we can build the power spectrum using them immediately. Fortunately the healpy package includes a function called anastat[8]. This subroutine calculates the $a_{\ell m}$ coefficients, as well as the $C_{\ell}$ power spectrum from an input $f_{p} \left( \theta_{p}, \varphi_{p} \right)$ HEALPix vector. After acquiring a HEALPix dataset of one arbitrary Planck CMB map, the angular power spectrum $C_{\ell}$ can be immediately reconstructed. However the $C_{\ell}$ quantity is what we call as “angular power spectrum”, it is traditional to plot another quantity against the multipoles $\ell$, which is only a function of $C_{\ell}$. This is the quantity
\begin{equation} D_{\ell} = \frac{\ell \left( \ell + 1 \right)}{2 \pi} C_{\ell}. \end{equation}
The reconstructed and calculated $D_{\ell}$ spectrum of the true CMB radiation can be seen on Figure 8., with the effect of instrumental noise seen on higher multipoles. The actual CMB temperature map used for the power spectrum extraction can be seen on Figure 1.
On the last figure I compared the obtained $D_{\ell}$ curve to the theoretical one. The latter was created with the CAMB software (Code for Anisotropies in the Microwave Background), which can be accessed from the website of NASA’s Goddard Space Flight Center, which institution maintains the (LAMBA) Legacy Archive for Microwave Background Data Analysis[10] project.
During the course I successfully accomplished every task of the project work, while learning a complete set of new methods in data analysis and visualization. I also got better insight about the behaviour and research of the Cosmic Microwave Background radiation.
The main reason of science to research the CMB radiation is to get more knowledge from it about the early stage and the structure of the universe. Currently the $\Lambda$CDM parametrization is considered to be (or referred to) the standard model of cosmology. Our main goal is to make the parameters in this model more accurate. Bring them closer to reality. More precise a physical model is the more better predictions it can give and the more accurately it can be used to describe and interpret the events and phenomenons in the universe around us.
As an interesting extra task I’ve tried to extract some cosmological parameters from an input power spectrum. For this I’ve used an MCMC algorithm and tried to converge the model parameters into a stable and meaningful position. Until the end of the semester I didn’t had enough time to get any meaningful results from this analysis, neither my simulations converged to stable values. The wandering of some example parameters during the MCMC iteration can be seen on Figure 9. Nonetheless this analysis can be created and rerun using the data that I’ve worked with and even meaningful cosmological parameters can be obtained. (But for now it’s up to someone else and not me…)