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.
- 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
ilc_clean_cij(cij, w)
: return a covariance matrix multiplied by weights,w' * cij * w
, withw
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 inilc_weights()
,cilc_weights()
ormilca_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.
- This would be a variance of the extracted component if
Arguments
cij::Array{<:AbstractFloat,2}
:nν
-by-nν
symmetric covariance matrix, wherenν
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()
andmilca_weights()
are multiple dispatch, which detect the dimension of the inputcij
automatically.
- or,
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 isnν
.b::Array{<:AbstractFloat,1}
: vector of the frequency response, for the component to be nulled. The number of elements isnν
.B::Array{<:AbstractFloat,2}
:nν
-by-nrc
matrix of the frequency response, fornrc
components to be nulled.
Optional arguments
e::Array{<:AbstractFloat,1}
: vector of the frequency response, for the component to be extracted. The default value ise=[1,...,1]
, i.e., CMB. The number of elements isnν
.ℓid::Integer=3
: location of the index for thenℓ
domain, if the dimension of the inputcij
is(nℓ, nν, nν)
,(nν, nℓ, nν)
or(nν, nν, nℓ)
.ℓid=1
ifcij[nℓ,nν,nν]
,ℓid=2
ifcij[nν,nℓ,nν]
, andℓid=3
(the default value) ifcij[nν,nν,nℓ]
.
Parametric Maximum Likelihood Method
loglike_beta(nij, A, d)
orloglike_beta(nij, A, cij)
:: return log(likelihood) for frequency response vectorsA
given a data vectord
or data covariance matrixcij
, based on Equation (9) of Stompor et al., MNRAS, 392, 216 (2009).
Arguments
nij::Array{<:AbstractFloat,2}
:nν
-by-nν
symmetric noise covariance matrix, wherenν
is the number of frequency bands.A::Array{<:AbstractFloat,2}
:nν
-by-nc
matrix of the frequency response, fornc
components in sky.- E.g.,
A = [a B]
wherea = ones(nν)
for CMB andB
is anν
-by-nc-1
matrix for the frequency response of foreground components.
- E.g.,
d::Array{<:AbstractFloat,1}
: data vector for a given pixel, (ℓ,m), or any other appropriate domain. The number of elements isnν
.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, setunits = "rj"
. The default is the CMB units.Tcmb::Real
: present-day temperature of the CMB in units of Kelvin. The default value is2.7255
.Td::Real
: dust temperature. The default value is19.6
(Kelvin).βd::Real
: dust emissivity index. The default value is1.6
.νd::Real
: frequency at which the output dust spectrum is normalized to unity. The default value is353
(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 is23
(GHz).Cs::Real
: curvature of the synchrotron spectrum. The default value is0
.νC::Real
: pivot frequency for curvature of the synchrotron spectrum. The default value is40
(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.
- This code applies
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.
-
- 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:
- Read in a hits map and calculate weights for inhomogeneous noise.
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.- Read in and smooth foreground maps by beams of the telescope.
- Generate Gaussian random realisations of the CMB and noise.
- Calculate covariance matrices
cij
of the foreground, noise, and CMB+foreground+noise maps. - Calculate ILC weights in harmonic domain using
ilc_weights()
. - Calculate power spectra of the clean CMB maps using
ilc_clean_cij()
. - Show results for visual inspection, if
showresults = true
. - Calculate the tensor-to-scalar ratio and its uncertainty from the simulated realisations.
- 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 figureilc_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
- This code applies
-
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 ILCmilca_weights()
. - When calculating the weights, it uses the total covariance matrix
cij
(like for the ILC) rather than the noise covariance matrixnij
(like for the maximum likelihood). This minimises further the foreground contribution that is not modeled.
- The code applies
-
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
- Calculate the best-fitting
βs
andβd
by minimising-loglike_beta()
with respect to them. - Calculate weights using
milca_weights()
with the best-fittingβs
andβd
, and calculate power spectra of the clean CMB maps usingilc_clean_cij()
.
- Calculate the best-fitting
-
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 bysmooth_FWHM
(in units of degrees; the default value issmooth_FWHM = 0.5
). This is equivalent to findingβs
andβd
on smoothed maps with the common resolution.
- The code finds the best-fitting
-
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
- Same as examples/ILCPipelineSOSAT.jl, but add simulated data of the CCAT-prime at 350, 410, and 850 GHz with specifications given in Table 1 of Choi et al., JLTP, 199, 1089 (2020).
-
examples/MILCAPipelineSOSATCCATp.jl
- Same as examples/MILCAPipelineSOSAT.jl, but vary the dust temperature
Td
as the third foreground parameter, and use Gaussian priors forβs
,βd
andTd
. - Also add simulated data of the CCAT-prime at 350, 410, and 850 GHz with specifications given in Table 1 of Choi et al., JLTP, 199, 1089 (2020).
- Same as examples/MILCAPipelineSOSAT.jl, but vary the dust temperature
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)
wherenν
is the number of frequency channels andnbands
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
andTd
with Gaussian priors.
- Note: examples/PerformMILCA.jl varies three parameters,
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.