## CleanCMB.jl

Algorithms to remove foreground emission and obtain clean maps of the cosmic microwave background
Author komatsu5147
Popularity
5 Stars
Updated Last
8 Months Ago
Started In
June 2020

# CleanCMB

This package contains functions to enable extraction of clean maps of the cosmic microwave background (CMB). Some functions can be used to extract non-CMB astrophysical components as well.

Different algorithms exist for extraction of clean maps of the CMB (as well as of astrophysical components). The package currently supports:

• Internal Linear Combination (ILC) Method
• Parametric Maximum Likelihood Method

See this note for the relationship between them.

## Internal Linear Combination (ILC) Method

• `ilc_weights(cij[, e, ℓid=3])`: return weights for the internal linear combination (ILC) method, following Equation (12) of Tegmark et al., PRD, 68, 123523 (2003).
• `cilc_weights(cij, a, b[, ℓid=3])`: return weights for the constrained internal linear combination (CILC) method for two components, following Equation (20) of Remazeilles, et al., MNRAS, 410, 2481 (2011).
• `milca_weights(cij, a, B[, ℓid=3])`: return weights for the modified internal linear combination algorithm (MILCA) method (CILC for N components), following Equation (15) of Hurier, et al., A&A, 558, A118 (2013).
• Note: These papers define weights in various domain, including harmonic, wavelet, and pixel domain. You can choose to work in any domain by supplying a covariance matrix `cij` in the appropriate domain.
• `ilc_clean_cij(cij, w)`: return a covariance matrix multiplied by weights, `w' * cij * w`, with `w` returned by any of the above functions, e.g., `w = cilc_weights(cij, a, b)`.
• This would be a variance of the extracted component if `cij` were the same as that used in `ilc_weights()`, `cilc_weights()` or `milca_weights()`.
• If `cij` is a noise covariance matrix, this function returns the noise variance of the cleaned map for each multipole, band-power bin, pixel, etc.

### Arguments

• `cij::Array{<:AbstractFloat,2}`: `nν`-by-`nν` symmetric covariance matrix, where `nν` is the number of frequency bands.
• or, `cij::Array{<:AbstractFloat,3}`:symmetric covariance matrix with the dimention of `(nℓ, nν, nν)`, `(nν, nℓ, nν)` or `(nν, nν, nℓ)` (default). Here, `nℓ` is the number of elements in the relevant domain, e.g., multipoles, band-power bins, pixels, etc.
• The functions `ilc_weights()`, `cilc_weights()` and `milca_weights()` are multiple dispatch, which detect the dimension of the input `cij` automatically.
• `a::Array{<:AbstractFloat,1}`: vector of the frequency response, for the component to be extracted. E.g., `a=[1,...,1]` for CMB. The number of elements is `nν`.
• `b::Array{<:AbstractFloat,1}`: vector of the frequency response, for the component to be nulled. The number of elements is `nν`.
• `B::Array{<:AbstractFloat,2}`: `nν`-by-`nrc` matrix of the frequency response, for `nrc` components to be nulled.

### Optional arguments

• `e::Array{<:AbstractFloat,1}`: vector of the frequency response, for the component to be extracted. The default value is `e=[1,...,1]`, i.e., CMB. The number of elements is `nν`.
• `ℓid::Integer=3`: location of the index for the `nℓ` domain, if the dimension of the input `cij` is `(nℓ, nν, nν)`, `(nν, nℓ, nν)` or `(nν, nν, nℓ)`. `ℓid=1` if `cij[nℓ,nν,nν]`, `ℓid=2` if `cij[nν,nℓ,nν]`, and `ℓid=3` (the default value) if `cij[nν,nν,nℓ]`.

## Parametric Maximum Likelihood Method

• `loglike_beta(nij, A, d)` or `loglike_beta(nij, A, cij)`:: return log(likelihood) for frequency response vectors `A` given a data vector `d` or data covariance matrix `cij`, based on Equation (9) of Stompor et al., MNRAS, 392, 216 (2009).

### Arguments

• `nij::Array{<:AbstractFloat,2}`: `nν`-by-`nν` symmetric noise covariance matrix, where `nν` is the number of frequency bands.
• `A::Array{<:AbstractFloat,2}`: `nν`-by-`nc` matrix of the frequency response, for `nc` components in sky.
• E.g., `A = [a B]` where `a = ones(nν)` for CMB and `B` is a `nν`-by-`nc-1` matrix for the frequency response of foreground components.
• `d::Array{<:AbstractFloat,1}`: data vector for a given pixel, (ℓ,m), or any other appropriate domain. The number of elements is `nν`.
• `cij::Array{<:AbstractFloat,2}`: `nν`-by-`nν` symmetric covariance matrix for a given multipole, pixel, or any other appropriate domain.

## Foreground models

The package contains the following functions to return frequency dependence of foreground components:

• `tsz(νGHz; units="cmb", Tcmb=2.7255)`: Spectrum of the thermal Sunyaev-Zeldovich effect given in Equation (V) in Appendix of Zeldovich, Sunyaev, Astrophys. Space Sci. 4, 301 (1969).
• `dust1(νGHz; Td=19.6, βd=1.6, νd=353, units="cmb", Tcmb=2.7255)`: Spectrum of 1-component modified black-body thermal dust emission. The output is normalized to unity at `νGHz = νd`.
• `synch(νGHz; βs=-3, νs=23, Cs=0, νC=40, units="cmb", Tcmb=2.7255)`: Spectrum of synchrotron emission. The output is normalized to unity at `νGHz = νs`.

### Arguments

• `νGHz::Real`: frequency in units of GHz.

### Optional keyword arguments

• `units::String`: units of the spectrum. For Rayleigh-Jeans temperature (brightness temperature) units, set `units = "rj"`. The default is the CMB units.
• `Tcmb::Real`: present-day temperature of the CMB in units of Kelvin. The default value is `2.7255`.
• `Td::Real`: dust temperature. The default value is `19.6` (Kelvin).
• `βd::Real`: dust emissivity index. The default value is `1.6`.
• `νd::Real`: frequency at which the output dust spectrum is normalized to unity. The default value is `353` (GHz).
• `βs::Real`: synchrotron power-law index. The default is `-3`.
• `νs::Real`: frequency at which the output synchrotron spectrum is normalized to unity. The default is `23` (GHz).
• `Cs::Real`: curvature of the synchrotron spectrum. The default value is `0`.
• `νC::Real`: pivot frequency for curvature of the synchrotron spectrum. The default value is `40` (GHz).

## Example Julia Codes

• examples/CleanWMAP.jl
• This code applies `ilc_weights()` in pixel domain to produce a clean map of CMB from the temperature maps of Wilkinson Microwave Anisotropy Probe (WMAP) at five frequency bands.
• The code requires Healpix.jl.

### Pipelines

Here we provide example codes for pipelines, which do everything from generation of simulated maps to estimation of the cosmological parameter (tensor-to-scalar ratio of the primordial gravitational waves) in one go.

• examples/ILCPipelineSOSAT.jl

• This code applies `ilc_weights()` in harmonic domain to produce power spectra of clean polarisation maps of the CMB.
• The code requires Healpix.jl and Libsharp.jl. It also calls python wrappers pymaster and classy via `PyCall`.
• This is a simulation pipeline for the Small Aperture Telescope (SAT) of the Simons Observatory. The code performs the following operations:
1. Read in a hits map and calculate weights for inhomogeneous noise.
2. `include("compute_cl_class.jl")` (see examples/compute_cl_class.jl) to generate theoretical scalar- and tensor-mode power spectra of the CMB using CLASS.
3. Read in and smooth foreground maps by beams of the telescope.
4. Generate Gaussian random realisations of the CMB and noise.
5. Calculate covariance matrices `cij` of the foreground, noise, and CMB+foreground+noise maps.
6. Calculate ILC weights in harmonic domain using `ilc_weights()`.
7. Calculate power spectra of the clean CMB maps using `ilc_clean_cij()`.
8. Show results for visual inspection, if `showresults = true`.
9. Calculate the tensor-to-scalar ratio and its uncertainty from the simulated realisations.
10. Write out the results (tensor-to-scalar ratios with and without residual foreground marginalisation) to a file `ilc_results_sosat.csv` and create a PDF figure `ilc_clbb_sim_sosat.pdf` showing the cleaned B-mode power spectrum, minus noisebias, and minus the residual foreground.
• For your reference, the results from 300 realisations starting with the initial seed of `5147` are given in examples/results/ilc_results_sosat_300sims_seed5147.csv. You can compute the mean and standard deviation of the tensor-to-scalar ratios. You should find, for 30 < ℓ < 260:
• Without FG marginalisation: r = (2.0 ± 1.6) x 10-3
• With marginalisation: r = (-1.0 ± 2.7) x 10-3
• examples/MILCAPipelineSOSAT.jl

• This code performs a hybrid of the maximum likelihood and ILC methods. See this note for the relationship between them.

• The code applies `loglike_beta()` in harmonic domain to find the best-fitting synchrotron and dust spectral indices (`βs` and `βd`), and calculates weights using the N-component constrained ILC `milca_weights()`.
• When calculating the weights, it uses the total covariance matrix `cij` (like for the ILC) rather than the noise covariance matrix `nij` (like for the maximum likelihood). This minimises further the foreground contribution that is not modeled.
• The code performs the same operations as the above ILC pipeline for the Small Aperture Telescope (SAT) of the Simons Observatory, but replaces the steps f and g with

1. Calculate the best-fitting `βs` and `βd` by minimising `-loglike_beta()` with respect to them.
2. Calculate weights using `milca_weights()` with the best-fitting `βs` and `βd`, and calculate power spectra of the clean CMB maps using `ilc_clean_cij()`.
• More detail of the procedure:

• The code finds the best-fitting `βs` and `βd` for each band-power at multipoles below the "switching" multipole, `ℓ ≤ ℓswitch` (the default value is `ℓswitch = 30`). This may better handle complex foreground properties on large angular scales.
• For higher multipoles, the code finds the global `βs` and `βd` using the covariance matrix smoothed to a common resolution specified by `smooth_FWHM` (in units of degrees; the default value is `smooth_FWHM = 0.5`). This is equivalent to finding `βs` and `βd` on smoothed maps with the common resolution.
• For your reference, the results from 300 realisations starting with the initial seed of `5147` are given in examples/results/milca_results_sosat_300sims_seed5147.csv. You can compute the mean and standard deviation of the tensor-to-scalar ratios. You should find, for 30 < ℓ < 260:

• Without FG marginalisation: r = (1.2 ± 2.7) x 10-3
• With marginalisation: r = (-0.8 ± 5.1) x 10-3
• examples/ILCPipelineSOSATCCATp.jl

• examples/MILCAPipelineSOSATCCATp.jl

### Splitting Pipelines

The above codes do everything in one go. However, sometimes it is convenient to split the processes. For example, we do not have to re-generate maps and their covariance matrices everytime we make a small modification to the rest of the pipeline.

Here we provide example codes for splitting the pipelines into two pieces: (1) Generation of simulated maps and their covariance matrices, and (2) Application of foreground cleaning methods to the covariance matrices.

• examples/GenerateCovMatrices.jl: Perform the steps (a)-(e) of the pipeline and write out the covariance matrices to binary files in arrays of `(nν, nν, nbands)` where `nν` is the number of frequency channels and `nbands` is the number of band-powers. It also writes out the binned scalar and tensor power spectra used to generate the simulations to text files.
• examples/PerformILC.jl and examples/PerformMILCA.jl: Perform the steps (f)-(j) of of the pipeline and write out the estimated tensor-to-scalar ratios to a csv file.
• Note: examples/PerformMILCA.jl varies three parameters, `βs`, `βd` and `Td` with Gaussian priors.

## Acknowledgment

Development of the functions provided in this package was supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC-2094 - 390783311 and JSPS KAKENHI Grant Number JP15H05896. The Kavli IPMU is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan.