CleanCMB
This package contains functions to enable extraction of clean maps of the cosmic microwave background (CMB). Some functions can be used to extract nonCMB 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, bandpower bin, pixel, etc.
 This would be a variance of the extracted component if
Arguments
cij::Array{<:AbstractFloat,2}
:nν
bynν
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, bandpower 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ν
bynrc
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ν
bynν
symmetric noise covariance matrix, wherenν
is the number of frequency bands.A::Array{<:AbstractFloat,2}
:nν
bync
matrix of the frequency response, fornc
components in sky. E.g.,
A = [a B]
wherea = ones(nν)
for CMB andB
is anν
bync1
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ν
bynν
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 SunyaevZeldovich 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 1component modified blackbody 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 RayleighJeans temperature (brightness temperature) units, setunits = "rj"
. The default is the CMB units.Tcmb::Real
: presentday 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 powerlaw index. The default is3
.ν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 (tensortoscalar 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 tensormode 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 tensortoscalar ratio and its uncertainty from the simulated realisations.
 Write out the results (tensortoscalar 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 Bmode 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 tensortoscalar 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 bestfitting synchrotron and dust spectral indices (βs
andβd
), and calculates weights using the Ncomponent 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 bestfitting
βs
andβd
by minimisingloglike_beta()
with respect to them.  Calculate weights using
milca_weights()
with the bestfittingβs
andβd
, and calculate power spectra of the clean CMB maps usingilc_clean_cij()
.
 Calculate the bestfitting

More detail of the procedure:
 The code finds the bestfitting
βs
andβd
for each bandpower 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 bestfitting

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 tensortoscalar 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 CCATprime 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 CCATprime 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 regenerate 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 bandpowers. 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 tensortoscalar 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  EXC2094  390783311 and JSPS KAKENHI Grant Number JP15H05896. The Kavli IPMU is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan.