Joint Bayesian analysis of large angular scale CMB temperature anomalies

TL;DR

This paper presents a comprehensive Bayesian analysis of large-scale anomalies in the CMB, specifically power suppression and hemispherical asymmetry, using Bipolar Spherical Harmonics to jointly analyze isotropy violations and power spectrum features.

Contribution

It introduces a joint Bayesian framework employing Bipolar Spherical Harmonics to analyze CMB anomalies and assesses the scale dependence of the hemispherical asymmetry.

Findings

The asymmetry amplitude at l=16 is approximately 0.064.

The power law index for scale dependence is about -0.92.

The data slightly favor anisotropic models but not conclusively.

Abstract

Cosmic microwave background measurements show an agreement with the concordance cosmology model except for a few notable anomalies: Power Suppression, the lack of large scale power in the temperature data compared to what is expected in the concordance model, and Cosmic Hemispherical Asymmetry, a dipolar breakdown of statistical isotropy. An expansion of the CMB covariance in Bipolar Spherical Harmonics naturally parametrizes both these large-scale anomalies, allowing us to perform an exhaustive, fully Bayesian joint analysis of the power spectrum and violations of statistical isotropy up to the dipole level. Our analysis sheds light on the scale dependence of the Cosmic Hemispherical Asymmetry. Assuming a scale-dependent dipole modulation model with a two-parameter power law form, we explore the posterior pdf of amplitude $A(l = 16)$ and the power law index $\alpha$ and find the maximum a posteriori values $A_*(l = 16) = 0.064 \pm 0.022$ and $\alpha_* = -0.92 \pm 0.22$. The maximum a posteriori direction associated with the Cosmic Hemispherical Asymmetry is $(l,b) = (247.8^o, -19.6^o)$ in Galactic coordinates, consistent with previous analyses. We evaluate the Bayes factor $B_{SI-DM}$ to compare the Cosmic Hemispherical Asymmetry model with the isotropic model. The data prefer but do not substantially favor the anisotropic model ($B_{SI-DM}=0.4$). We consider several priors and find that this evidence ratio is robust to prior choice. The large-scale power suppression does not soften when jointly inferring both the isotropic power spectrum and the parameters of the asymmetric model, indicating no evidence that these anomalies are coupled.

Tables

Table 1. Table 1 : Analysis details and summaries of parameter posterior marginals (Step model)

Details of the map used in analysis: NSIDE = 256, $f_{sky}=84\%$, $l_{cut}=64$
Parameter	$m_{10}$	$m_{11}^r$	$m_{11}^i$	$A_{\ast }$	$l$	$b$
Mean	-0.039	0.035	-0.056	0.050	${238.0}^o$	$-{22.4}^o$
Standard Deviation	0.038	0.029	0.029	0.018	30.4	20.6

Table 2. Table 2 : Analysis details and summaries of parameter posterior marginals (Power law model)

Prior on $A\left(l_p\right)$ is chosen such that (see section 6 for details)
Prior on $\alpha$: Uniform prior over the range $\alpha =-2$ to $\alpha =0$
Details of the map used in analysis: NSIDE = 256, $l_{max}=256$, $f_{sky}=84\%$
Parameter	$m_{10}\left(l_p=16\right)$	$m_{11}^r\left(l_p\right)$	$m_{11}^i\left(l_p\right)$	$\alpha$	$A_{\ast }\left(l_p\right)$	$l$	$b$
Mean	-0.044	0.033	-0.081	-0.92	0.064	${247.8}^o$	$-{19.6}^o$
Standard Deviation	0.043	0.031	0.032	0.22	$0.022$	22.0	18.0

Table 3. Table 3 : Savage-Dickey density ratio (SDDR) for different prior densities

Prior choice	$(w_x,w_y,w_z)$	SDDR
${\Pi }_2$	$(ϵ,ϵ,ϵ),ϵ={10}^{-3}$	0.4
${\Pi }_3$	$(ϵ,ϵ,ϵ),ϵ={10}^{-3}$	0.5

Equations

$$\Delta T(\hat{n})=(1+A\hat{p}\cdot \hat{n})\Delta T^{SI}(\hat{n}).$$

$$a_{lm}=a_{lm}^{SI}+\sum _{N=-1}^1{m}_{1N}(l)\int \Delta T^{SI}(\hat{n})Y_{L=1,N}(\hat{n})Y_{lm}^{\ast }(\hat{n})d^2{\Omega }_{\hat{n}},$$

$$A(l)\hat{p}\cdot \hat{n}=\sum _{N=-1}^1{m}_{1N}(l)Y_{L=1,N}(\hat{n}).$$

$$A(l)=\sqrt{\frac{3}{4\pi}}\sqrt{m_{10}(l)^{2}+2m^{r}_{11}(l)^{2}+2m^{i}_{11}(l)^{2}},\quad\theta_{p}=\cos^{-1}\Big{[}\frac{m_{10}(l)}{A(l)}\sqrt{\frac{3}{4\pi}}\Big{]},\quad\phi_{p}=-\tan^{-1}\Big{[}\frac{m^{i}_{11}(l)}{m^{r}_{11}(l)}\Big{]}.$$

$$A(l)=A(l_{p})\Big{(}\frac{l}{l_{p}}\Big{)}^{\alpha},$$

$$m_{1N}(l)=m_{1N}(l_{p})\Big{(}\frac{l}{l_{p}}\Big{)}^{\alpha},$$

$$S_{l_1{m}_1{l}_2{m}_2}$$

$$S_{l_1{m}_1{l}_2{m}_2}=(-1{)}^{m_2}\sum _{LN}{A}_{l_1{l}_2}^{LN}{C}_{l_1{m}_1{l}_2-m_2}^{LN}.$$

$$A_{ll}^{00}=(-1{)}^l\sqrt{2l+1}{C}_l.$$

$$⟨a_{lm}{a}_{lm}^{\ast }⟩=C_l+\frac{{\Pi }_l^2}{12\pi }\sum _{l^{\prime }{m}^{\prime }}{C}_{l^{\prime }}{\Pi }_{l^{\prime }}^2[C_{l0l^{\prime }0}^{10}{]}^2\sum _{N_1,N_2=-1}^1{m}_{1N_1}(l^{\prime })m_{1N_2}^{\ast }(l^{\prime })C_{lm{l}^{\prime }{m}^{\prime }}^{1N_1}{C}_{lm{l}^{\prime }{m}^{\prime }}^{1N_2}.$$

$${\tilde{C}}_l=C_l+\frac{C_l}{4\pi }\sum _{N=-1}^1|m_{1N}{|}^2,$$

$${\Delta }_l\equiv \frac{{\tilde{C}}_l-C_l}{C_l}=\frac{1}{4\pi }\sum _{N=-1}^1|m_{1N}{|}^2=\frac{A^2(l)}{3}.$$

$$A_{ll+1}^{1N}(l_p,\alpha )=m_{1N}(l_p)G_{ll+1}^1,$$

$$G^{1}_{ll+1}\equiv\frac{\Pi_{ll+1}}{\sqrt{12\pi}}\Big{[}\Big{(}\frac{l}{l_{p}}\Big{)}^{\alpha}C_{l}+\Big{(}\frac{l+1}{l_{p}}\Big{)}^{\alpha}C_{l+1}\Big{]}C^{10}_{l0l+10}.$$

$$\mathcal{P}(\textbf{d}|\textbf{a})=\frac{1}{\sqrt{|\textbf{N}|}(2\pi)^{n/2}}\exp\Big{[}-\frac{1}{2}(\textbf{d}-\textbf{a})^{\dagger}\textbf{N}^{-1}(\textbf{d}-\textbf{a})\Big{]}.$$

$$\mathcal{P}(\textbf{a}|\textbf{S})=\frac{1}{\sqrt{|\textbf{S}|}(2\pi)^{n/2}}\exp\Big{[}-\frac{1}{2}\textbf{a}^{\dagger}\textbf{S}^{-1}\textbf{a}\Big{]},$$

$$\mathcal{P}(\text{S},\text{a}|\text{d})=\frac{\mathcal{P}(\text{d}|\text{a})\mathcal{P}(\text{a}|\text{S})\mathcal{P}(\text{S})}{\mathcal{P}(\text{d})}.$$

$$\mathcal{P}(\textbf{S},\textbf{a}|\textbf{d})=\frac{1}{\sqrt{|\textbf{N}||\textbf{S}|}(2\pi)^{n}}\exp{\Big{\{}-\frac{1}{2}\Big{[}(\textbf{d}-\textbf{a})^{\dagger}\textbf{N}^{-1}(\textbf{d}-\textbf{a})+\textbf{a}^{\dagger}\textbf{S}^{-1}\textbf{a}\Big{]}\Big{\}}}.$$

$$P(\textbf{S}|\textbf{d})=\frac{1}{\sqrt{|\textbf{N}+\textbf{S}|}(2\pi)^{n/2}}\exp{\Big{[}-\frac{1}{2}\textbf{d}^{T}(\textbf{S}+\textbf{N})^{-1}\textbf{d}\Big{]}}.$$

$$P(C_{l},m_{10}|\textbf{d})=\frac{1}{\sqrt{|\textbf{D}|}(2\pi)^{n/2}}\exp{\Big{[}-\frac{1}{2}\textbf{d}^{\dagger}\textbf{D}^{-1}\textbf{d}\Big{]}}\exp{\Big{[}-\frac{(m_{10}-\mu)^{2}}{2\sigma^{2}}\Big{]}}\exp{\Big{[}\frac{\mu^{2}}{2\sigma^{2}}\Big{]}},$$

$$\mu =\frac{{\text{d}}^{†}{\text{D}}^{-1}{\text{O}}_1{\text{D}}^{-1}\text{d}}{Tr[({\text{D}}^{-1}{\text{O}}_1{)}^2]}{M}^{-1}\text{and}{\sigma }^2=\frac{2M^{-1}}{Tr[({\text{D}}^{-1}{\text{O}}_1{)}^2]},$$

$$M=\Big{[}\frac{2\textbf{d}^{\dagger}(\textbf{D}^{-1}\textbf{O}_{1})^{2}\textbf{D}^{-1}\textbf{d}}{Tr[(\textbf{D}^{-1}\textbf{O}_{1})^{2}]}-1\Big{]}.$$

$$Tr[(\textbf{D}^{-1}\textbf{O}_{1})^{2}]=\frac{1}{4\pi}\sum^{l_{max}}_{l=l_{min}}\Big{[}\Big{(}2+\frac{D_{l-1}}{D_{l}}+\frac{D_{l}}{D_{l-1}}\Big{)}l+\Big{(}2+\frac{D_{l+1}}{D_{l}}+\frac{D_{l}}{D_{l+1}}\Big{)}(l+1)\Big{]},$$

$$w_z=m_{10},w_x=-\sqrt{2}{m}_{11}^r,\text{and}{w}_y=\sqrt{2}{m}_{11}^i.$$

$$r\equiv \sqrt{w_z^2+w_x^2+w_y^2}=\sqrt{m_{10}^2+2{m_{11}^r}^2+2{m_{11}^i}^2}.$$

$$\mathcal{F}(r)=\sqrt{\frac{1}{(2\pi\sigma^{2})^{3}}}4\pi r^{2}\exp{\Big{[}-\frac{r^{2}}{2\sigma^{2}}\Big{]}},$$

$$\mathcal{F}(A)=\frac{8\pi}{3}\sqrt{\frac{2}{3}}\frac{A^{2}}{\sigma^{3}}\exp{\Big{[}-\frac{2\pi A^{2}}{3\sigma^{2}}\Big{]}}.$$

$$\mathcal{F}(r)=\frac{r}{\sigma^{2}}\sqrt{\frac{r}{r_{*}}}\exp\Big{[}-\frac{r^{2}+r^{2}_{*}}{2\sigma^{2}}\Big{]}I_{1/2}(\frac{rr_{*}}{\sigma^{2}}),$$

$$I_{1/2}(x)=\sqrt{\frac{x}{2\pi}}\int^{\pi}_{0}d\theta\sin(\theta)\exp\Big{[}x\cos(\theta)\Big{]}.$$

$$\mathcal{F}(A)=\frac{4\pi}{3}\frac{A}{\sigma^{2}}\sqrt{\frac{A}{A_{*}}}\exp\Big{[}-\frac{2\pi(A^{2}+A^{2}_{*})}{3\sigma^{2}}\Big{]}I_{1/2}(\frac{4\pi AA_{*}}{3\sigma^{2}}),$$

Full text

1 Introduction

The statistical homogeneity and isotropy of the universe are fundamental assumptions underlying cosmology. These assumptions have been put to the test using various cosmological probes [1, 2, 3, 4, 5, 6, 7]. Precise and almost full sky measurements of CMB temperature anisotropy provide an unique opportunity to carry out tests of statistical isotropy. As a consequence, soon after WMAP published its first results, several studies that submitted the data to tests of statistical isotropy, uncovered hints of some anomalous features in the CMB data [1, 2, 3]. Prominent among the CMB anomalies is the Power Asymmetry or the Cosmic Hemispherical Asymmetry (CHA), which is characterized by a $\approx 14\%$ excess of power in one hemisphere of CMB temperature fluctuations compared to the other and points to a possible challenge to the assumption of statistical isotropy in cosmology [1, 2]. CHA was later confirmed to also exist in the Planck CMB map at similar amplitude as seen in the WMAP data [5, 4]. Further, detailed analysis enabled in the Bipolar Spherical Harmonic representation by the Planck collaboration confirmed its dipolar nature, frequency independence and the fact that it existed at low multipoles (large angular scales) and died off at high multipoles (small angular scales) [5, 4]. Various physical models have been put forth for the plausible origin of CHA in the CMB, see for example [8, 9, 10, 11, 12, 13].

Apart from the presence of CHA at large-angular scales, temperature data from WMAP and Planck also indicates lower temperature fluctuation on large angular scales than predicted by the concordance Lambda-Cold Dark Matter ($\Lambda$CDM) model with power law form of the power spectrum of initial density fluctuations [14, 15, 16, 17]. An existing body of literature also explores connections between different anomalies. For example, Muir et al. (2018) [18] investigate the covariance between various anomalies. Polastri et al. (2015) [19] study the connection between directional anomaly of low multipole alignment and dipole modulation. Schwarz et al. (2016) [20] investigate different anomalies to find features that can explain more than one observed anomalous signal.

In this paper, we revisit the Planck-2015 temperature data [25] to scrutinize both these anomalies jointly within a common mathematical framework of Bipolar Spherical Harmonics (BipoSH). BipoSH [26, 27] forms a complete basis for the two point function on the sphere and captures the entire structure of the covariance matrix. Hence, it is a natural choice of basis for this analysis. Any similar estimator of the two-point statistics in spherical harmonic space for a random field can be represented in the BipoSH space. The angular power spectrum is represented by the $L=0$ BipoSH coefficients. We explore the probability distribution of $L=1$ (dipole) BipoSH coefficients which is the next order term in the BipoSH expansion of the covariance matrix, hence avoiding a posteriori choice. We perform the Bayesian inference of $L=1$ BipoSH coefficients jointly with the angular power spectrum ($L=0$ BipoSH coefficients). This is the first time a joint analysis of two anomalies in observed CMB data is being done without assuming a particular shape of the angular power spectrum. We extend the Hamiltonian Monte Carlo (HMC) sampling method [28, 29] in the BipoSH framework by Das et al. [30] to incorporate essential observational non-idealities in the Planck data 111While this work was in progress, one of the coauthors (Santanu Das) has posted two arXiv preprints [31, 32] (PDF file of one of those was withdrawn later [31]) which are related to the HMC analysis for statistically anisotropic signal in the presence of partial sky and non-isotropic noise. However, the analysis done in this paper did not use the code made public in [32]. The method and the code used in this analysis are developed independently from the initial code of [30] for full sky and isotropic noise.. We test our algorithm on simulated maps generated using CoNIGS algorithm [33] and then apply it to the Planck-2015 temperature data. Our joint estimation validates the existence of the anomalies in the Planck temperature data and recovers a non-trivial and physically interesting profile for the scale-dependent CHA signal when modeled as a power-law form. The recovered scale-dependent CHA profile is one of the key findings of this analysis and is depicted in figure 1. The details of the analysis and its significance are discussed in section 5. The obtained power-law profile sheds light onto the scale-dependent nature of the signal and goes beyond the simplistic step-model with an arbitrarily chosen cut-off [2, 3, 5, 4] as shown by the blue curve in figure 1. For comparison, we also plot the scale-independent signal due to our local motion [23] which also exhibits a dipolar anisotropy in the CMB data. We use the parameter samples drawn from the probability distribution to perform the model comparison between statistically isotropic (SI) model and power law dipole modulation model of CHA using the ratio of Bayesian evidence obtained under two models. We note that these parameter samples are obtained after marginalizing over the angular power spectrum instead of simply conditioning on a fiducial angular power spectrum. The prior on the CHA parameters is carefully constructed so as to maintain the positive definite nature of the covariance matrix, hence not favoring an artificially large amplitude of CHA signal.

We present the paper in the following manner. In section 2 we provide an introduction to the subject of this paper, Cosmic Hemispherical Asymmetry (CHA), and describe the phenomenological model of CHA and its imprint on the covariance matrix of the CMB temperature fluctuations. Section 3 provides an analytic insight into the probability distributions of the quantities relevant to the problem. These are the probability distributions that are sampled using the HMC algorithm. Details of the HMC sampling method are described in section 4. In appendix C we demonstrate the application of the method to simulated maps. Section 5 and 6 discuss results from the analysis of Planck-2015 SMICA temperature map and the corresponding Bayesian evidence in comparison to the statistically isotropic model respectively. Finally, we conclude in section 7.

2 Cosmic Hemispherical Asymmetry

2.1 Introduction to Dipole Modulation Model

CHA is modeled as dipole modulation of statistically isotropic (SI) CMB sky [34]. Dipole modulated CMB temperature anisotropy field ($\varDelta T(\hat{n})\equiv T(\hat{n})/T_{0}-1$)222$T_{0}=2.7255$ K is the all sky average temperature of CMB [35, 36]. in the direction $\hat{n}$, is given by

[TABLE]

Here, $\varDelta T^{SI}(\hat{n})$ is SI CMB temperature anisotropy field. The direction of the modulation dipole is given by $\hat{p}\equiv(\theta_{p},\phi_{p})$ and $A$ is the amplitude of the dipole modulation. Within this phenomenological model, characterizing the signature of CHA then reduces to estimating these three parameters that capture the departure from statistical isotropy. Various methods have been used to study the CHA observed in the CMB. One approach is to construct an estimator for the amplitude and the direction of the dipole. Akrami et al. (2014) have used local variance estimator to look for such signal in CMB data [37]. This estimator searches for a dipole in the “local variance map” of temperature fluctuations. In Hanson et al. (2009), Quadratic Maximum Likelihood estimator is constructed for these parameters [38]. For analysis based on BipoSH representation, the contribution of CHA is captured in $L=1$ BipoSH coefficients [39, 40, 41]. Another approach is to get the likelihood distribution of these parameters. Hoftuft et al. (2009) sample the likelihood distribution in real space, where the likelihood is defined in terms of real space parameters $A$ and $(\theta_{p},\phi_{p})$. The likelihood is sampled using Markov Chain Monte Carlo method to estimate the parameters of the modulation dipole [3]. Results of these different estimators are also included in the Planck-2013 [5] and Planck-2015 [4] analysis. Different methods have detected the signal of CHA at a significance level around $3\sigma$ with a consistent direction.

It has been now established that the dipole modulation signal is significant only at large angular scales and becomes insignificant at small angular scales [3, 42, 38]. If we allow the amplitude of the dipole to vary with angular scale, it is called scale dependent dipole modulation. Given as it is written, eq. (2.1) can not explicitly express the scale dependence of the dipole amplitude. It is best expressed in harmonic space using following equation

[TABLE]

where $a^{SI}_{lm}$ and $a_{lm}$ are the spherical harmonic coefficients of $\varDelta T^{SI}(\hat{n})$ and $\varDelta T(\hat{n})$, respectively, and $m_{1N}(l)$ (N = -1, 0, 1) are the spherical harmonic coefficients of the dipole field whose amplitude depends on the multipole $l$. From eq. (2.2), it is clear that the SI temperature field $\Delta T^{SI}$ also contributes to the statistically anisotropic (non-SI) part of the signal (second term), and hence a joint estimation of the $a_{lm}$ along with the $m_{1N}$ term is required to adequately capture the signal of CHA.

Here, $m_{11}$ and $m_{1-1}$ are complex numbers, that are related to each other by $m^{*}_{1-1}=-m_{11}$. It is convenient to use $m_{10}$, the real part of $m_{11}$ ($m^{r}_{11}$) and the imaginary part of $m_{11}$ ($m^{i}_{11}$) as the three equivalent independent variables. These three real numbers completely determine the direction $\hat{p}$ and the amplitude $A$ of the dipole. The scale-dependent modulation field using real space variables amplitude $A(l)$ and the dipole direction $\hat{p}$ is expressed as

[TABLE]

The amplitude $A$ and the direction $(\theta_{p},\phi_{p})$ in terms of $m_{10}$, $m^{r}_{11}$ and $m^{i}_{11}$ are given by following equations

[TABLE]

In most of the previous studies, a step model of the scale dependence was assumed, wherein one assumes that the dipole modulation amplitude is nonzero with a constant value below a certain multipole $l_{cut}$ and is zero for $l>l_{cut}$. In our work, we study the power law model for the dipole modulation amplitude $A(l)$

[TABLE]

where $A(l_{p})$ is the amplitude at a chosen pivot multipole $l_{p}$. The power law dependence of the asymmetry results in a power-law form harmonic space modulation field, that can be expressed as

[TABLE]

where $m_{1N}(l_{p})$ are the spherical harmonic coefficients of the dipole $A(l_{p})\hat{p}\cdot\hat{n}$, that are related to each other by eq. (2.4).

The description of a non-SI CMB temperature anisotropy field that captures the power law dipole modulation model requires four additional parameters along with the usual spherical harmonic coefficients $a_{lm}$. In this analysis, we sample the joint probability distribution of $a_{lm}$, $C_{l}$, ${m_{10},m^{r}_{11},m^{i}_{11}}$ and $\alpha$ under the Bayesian framework using HMC method [28, 30]. Previously Aiola et al. [24] studied a power-law scale dependence form to estimate only the CHA signal. However, due to the structural form of eq. (2.2), it is important to perform a joint analysis of both the angular power spectrum and CHA. We provide the first joint analysis in this paper accounting for all the important observational non-idealities such as partial sky coverage, anisotropic noise and instrumental beam response function.

2.2 Covariance matrix of dipole modulated CMB temperature sky

Since the CMB anisotropy is predicted and observed to be consistent with a Gaussian random field [4, 5], its covariance matrix (S) should encode all of its information content. For an SI CMB sky, the covariance matrix in spherical harmonic space is diagonal. However, in the presence of any kind of non-SI signal (a specific form of dipole modulation is given in eq. (2.2)), the covariance matrix of $a_{lm}$ contains non-zero off-diagonal terms. The exact333To all orders in $m_{1N}$ covariance matrix of $a_{lm}$ for the dipole modulation case can be written as

[TABLE]

where $C_{l}$ is angular power spectrum of $a^{SI}_{lm}$. Here, $\Pi_{l_{1}l_{2}\ldots l_{n}}$ denotes $\sqrt{(2l_{1}+1)(2l_{2}+1)\ldots(2l_{n}+1)}$ and $C^{1N}_{l_{1}m_{1}l_{2}m_{2}}$ are well-known Clebsch-Gordon coefficients.

We make use of BipoSH representation of the CMB two-point correlation function to carry out the analysis. The covariance matrix, $S_{l_{1}m_{1}l_{2}m_{2}}$, in terms of BipoSH coefficients can be expressed as [26]

[TABLE]

The BipoSH coefficients $A^{LN}_{l_{1}l_{2}}$ act as equivalent auxiliary variables in our calculations and computations. For $L=0$ these capture the SI correlations and are trivially related to the angular power spectrum, $C_{l}$, through the relation

[TABLE]

The leading order term of $m_{1N}$ in eq. (2.7) contributes only to the off-diagonal terms of the covariance matrix. The off-diagonal elements of the covariance matrix has the dependence on the magnitude and direction of the modulation dipole through $m_{1N}$. The second order term of eq. (2.7) is generally neglected on account of $m_{1N}$ being very small. However, this term has to be treated carefully because it leads to the modification of the diagonal of the covariance matrix that is used in the estimation of the angular power spectrum.

Diagonal terms: Following is the modification to the diagonal element due to the second order term in $m_{1N}$

[TABLE]

Note that $\langle a_{lm}a^{*}_{lm}\rangle$ depends also on $m$ (and not only on $l$) through the Clebsch-Gordan coefficients. For a smooth functional form of $m_{1N}(l)$ ($m_{1N}(l+1)\approx m_{1N}(l)$), we can simplify eq. (2.10) as

[TABLE]

where $\tilde{C}_{l}\equiv\langle a_{lm}a^{*}_{lm}\rangle$. Using eq. (2.3), the relative change in $C_{l}$ in terms of $A(l)$ can be written as

[TABLE]

In our analysis, we neglect the second order term in the covariance matrix, but provide an estimate of the relative correction in Sec. 5. It is shown to be sub-dominant and hence negligible.

Off-diagonal terms: The $L=1$ BipoSH coefficients for the dipole modulation temperature field can be deduced from Eq. 2.7 to be

[TABLE]

where $G^{1}_{ll+1}$ is referred to as the shape factor of the specific non-SI signal. Specifically, here

[TABLE]

eq. (2.13) when substituted in eq. (2.8) provides the off-diagonal terms of the covariance matrix in terms of $m_{1N}(l_{p})$ and $G^{1}_{ll+1}$. Scale dependence is an important aspect of this empirical modeling of the CHA that is captured by the index $\alpha$ of the power law. The shape factor encodes the signature of the dipole modulation through the mixing term between $l$ and $l\pm 1$ modes, which is sufficient to measure the non-SI signal from the map.

3 Probability distributions

In this section, we define the probability distributions of variables involved in the problem. We denote444Bold fonts are used to denote matrices. the diagonal of S by D and the off-diagonal part by O (so that $\textbf{S}=\textbf{D}+\textbf{O}$). We express the set of $a_{lm}$ by the vector $\mathbf{a}$. For the dipole modulation signal, eq. (2.8) and eq. (2.13) jointly express O in terms of $C_{l}$, $m_{1N}(l_{p})$ and $\alpha$.

3.1 Joint Probability distribution of $\textbf{S}(C_{l},m_{1N},\alpha)$ and

$\{a_{lm}\}$

We assume that the noise in the measurement of CMB temperature is Gaussian distributed with the noise covariance matrix N. The probability of the data d given $\textbf{a}\equiv\{a_{lm}\}$ is

[TABLE]

Here, $\textbf{N}^{-1}$ is the inverse of the noise covariance matrix N and $n$ is total number of $a_{lm}$ coefficients, $(n=l^{2}_{max}+2l_{max}-3)$, excluding monopole and dipole. The probability distribution of the zero mean Gaussian ${a_{lm}}$ can be written as

[TABLE]

where $\textbf{S}^{-1}$ denotes the inverse of the covariance matrix S. Using Bayes theorem, the joint probability distribution of S and a given d in terms of above two probability distributions is

[TABLE]

With uniform prior $\mathcal{P}(S)$ and up to a normalization constant $\mathcal{P}(\textbf{d})$,

[TABLE]

In above expression, the dependence of $\mathcal{P}(\textbf{S},\textbf{a}|\textbf{d})$ on $a_{lm}$ is explicit. $\mathcal{P}(\textbf{S},\textbf{a}|\textbf{d})$ depends on angular power spectrum $C_{l}$, $m_{1N}(l_{p})$ and $\alpha$ through S. The angular power spectrum, $C_{l}$, is the diagonal part of S. Dependence of S on $m_{1N}(l_{p})$ and $\alpha$ is given by eq. (2.7).

3.2 Probability distribution of the dipole modulation amplitude

We derive an analytical expression for the probability distribution of $m_{1N}$ and the dipole modulation amplitude $A$ with some approximations. In literature, it has been argued (on heuristic grounds) that $m_{1N}$ have Gaussian distribution (see section 6.3 of [4]). As a result of this, $A$ will have Maxwell-Boltzmann distribution if the distribution of $m_{1N}$ is centered at the origin. The probability distributions of $m_{1N}$ and $A$ play an important role in the interpretation of our results.

We integrate out a from eq. (3.4) and obtain the joint probability distribution of $C_{l}$ and $m_{1N}$. Marginalizing a from eq. (3.4) yields

[TABLE]

The explicit expression for the joint probability distribution of $C_{l}$ and $m_{10}$ can be written as

[TABLE]

where

[TABLE]

with the quantity $M$ given by

[TABLE]

The detailed derivation of the eq. (3.6) is given in appendix A. $M$ is a real number and is approximately equal to 1, which can be checked by replacing $M$ by its ensemble average. The expression for $Tr[(\textbf{D}^{-1}\textbf{O}_{1})^{2}]$ is

[TABLE]

where $D_{l}=C_{l}+N_{l}$ and $(l_{min},l_{max})$ is the range of multipoles that are dipole modulated. Inspection of eq. (3.6) reveals that $m_{10}$ is Gaussian distributed with mean $\mu$ and variance $\sigma^{2}$. It is informative to have a qualitative look at the expression for $\sigma^{2}$. As $C_{l}\approx C_{l+1}$ and CMB measurement are noise sub-dominant at large angular scales ($N_{l}<<C_{l}$), the variance $\sigma^{2}$ is approximately equal to $2\pi/(l^{2}_{max}+2l_{max}-l^{2}_{min}+1)$. This implies (as expected) that the variance is inversely proportional to the total number of independent modulated modes. This compact form is obtained for the scale independent $m_{1N}$. However, the same essence is also preserved for any scale dependent modulation field.

The components $m_{10}$ and $m_{11}$ have same variance, hence $m^{r}_{11}$ and $m^{i}_{11}$ have $\sigma^{2}/2$ as their variance. For the purpose of discussing the probability distribution, it is convenient to deal with the variables that have the same variance. So, we define three variables $w_{x}$, $w_{y}$, and $w_{z}$ as

[TABLE]

Then, $w_{x},w_{y},w_{z}$ are distributed with same variance $\sigma^{2}$. The norm of the signal $r$ can be defined as

[TABLE]

For a SI CMB map, $w_{x}$, $w_{y}$, and $w_{z}$ are Gaussian distributed with zero mean. Hence $r$ is distributed according to Maxwell-Boltzmann distribution

[TABLE]

where $\sigma^{2}$ is given by eq. (3.7). The relation $r=\sqrt{\frac{4\pi}{3}}A$ leads us to the probability distribution of the dipole amplitude $A$ in case of SI map as

[TABLE]

For the dipole modulated CMB sky, $w_{x}$, $w_{y}$, and $w_{z}$ are Gaussian distributed but are not centred at origin. We denote the Maximum Likelihood values of $w_{x}$, $w_{y}$, and $w_{z}$ by $w_{x*}$, $w_{y*}$, and $w_{z*}$ respectively, and the corresponding $r_{*}$ denotes the signal amplitude through eq. (3.11). Then the distribution of $r$ marginalized over $(\theta,\phi)$ is

[TABLE]

where $I_{1/2}(x)$ is a modified Bessel function of first kind of order $1/2$,

[TABLE]

The probability distribution for $A$ is

[TABLE]

where $A_{*}$ is the dipole amplitude corresponding to the radial coordinate $r_{*}$. For an SI map, $A_{*}$ and $r_{*}$ are equal to zero. It is important to note that $A_{*}$ is not same as the maximum probable value of $A$ where $\mathcal{F}(A)$ attains its maximum.

4 Method: Hamiltonian Monte Carlo sampling for non-SI CMB sky in presence of observational non-idealities

For the given problem, we deal with $\approx l^{2}_{max}+3l_{max}$ parameters. Markov Chain Monte Carlo method with the Metropolis-Hastings algorithm is known to have a much smaller acceptance rate in high dimensional parameter space; hence it is not practical for our task. This problem is circumvented in the Gibbs sampling method, which has been used in the CMB analysis for the estimation of the angular power spectrum [43, 44] as well as component separation [45]. The Gibbs sampling uses the conditional distributions of the parameters as proposal distributions. However, in our problem, with the addition of dipole modulation parameters, it is not straightforward to get the conditional distributions and draw samples from those distributions. In comparison, HMC is an efficient method to draw samples from high dimensional distributions. It makes use of Hamiltonian dynamics to propose the next sample, as discussed next in this section.

In this section, we describe the Hamiltonian Monte Carlo (HMC) method employed to sample the posterior probability distribution. We give details of the analytical computations required to implement this Monte Carlo method. The mask, anisotropic noise, and the beam need to be accounted for to analyze the Planck CMB maps. The validation of this method is performed on simulated non-SI maps obtained using CoNIGS [33] as described in Appendix C.

4.1 Review of the basic formalism as a HMC sampling problem

HMC has been applied to various problems in cosmology, notably for cosmological parameter estimation in [46], for CMB power spectrum estimation in [47], for estimation of the large scale structure power spectrum in [48, 49], the inference of non-linear dynamics of large scale structure [50, 51], and lensing potential reconstruction from lensed CMB in [52]. HMC makes use of Hamiltonian Dynamics to sample a given probability distribution $\mathcal{P}(q)$ of random variable $q$ using a Hamiltonian ($\mathcal{H}$) defined in the following way

[TABLE]

where $q$ represents the position of this particle and $p$ is the conjugate momentum corresponding to $q$. The quantity $-\ln(\mathcal{P}(q))$ acts as the potential energy and $\mu$ is “the mass” assigned to “the particle”. In HMC, the choice of the mass parameter, $\mu$ is specific to the problem (akin to the width of the proposal distribution in Markov Chain Monte Carlo method). Hamiltonian dynamics stipulates the time evolution of $q$ and $p$ through the Hamilton’s equations

[TABLE]

For the detailed steps of HMC we follow the algorithm given in [46] and to evolve the sympletic Hamiltonian dynamics equations we use the Forest-Ruth algorithm (see [30] for details of the implementation). HMC algorithm generates the samples of $p$ and $q$ drawn from the following probability distribution

[TABLE]

In the end, we marginalize over $p$ and obtain a fair sample of $q$ drawn from $\mathcal{P}(q)$.

In our problem, the variables of dynamics are the parameters ($a^{r}_{lm}$, $a^{i}_{lm}$, $C_{l}$, $m_{10}$, $m^{r}_{11}$, $m^{i}_{11}$, $\alpha$) and corresponding momentum ($p^{r}_{lm}$, $p^{i}_{lm}$, $p_{C_{l}}$, $\bar{p}_{10}$, $\bar{p^{r}}_{11}$, $\bar{p^{i}}_{11}$, $p_{\alpha}$). For the probability distribution defined in eq. (3.4)

[TABLE]

In above equation, $\textbf{p}_{lm}$ and $\textbf{p}_{C_{l}}$ stand for the set of conjugate momenta for $a_{lm}$ and $C_{l}$ respectively. $\bar{\textbf{p}}_{1N}$ stands for conjugate momenta $\{\bar{p}_{10},\bar{p}^{r}_{11},\bar{p}^{i}_{11}\}$ of ($m_{10},m^{r}_{11},m^{i}_{11}$) and $p_{\alpha}$ is conjugate momentum of $\alpha$. Further, $\mu^{r}_{a_{lm}},\mu^{i}_{a_{lm}},\mu_{C_{l}}$ are the masses assigned to the real part of $a_{lm}$, the imaginary part of $a_{lm}$ and $C_{l}$, respectively. Whereas, $\bar{\mu}_{10}$, $\bar{\mu}^{r}_{11}$, $\bar{\mu}^{i}_{11}$, and $\mu_{\alpha}$ are the masses assigned to $m_{10}$, $m^{r}_{11}$, $m^{i}_{11}$, and $\alpha$, respectively. We follow arguments given in [30] and [47] to choose masses for $a^{r}_{lm}$, $a^{i}_{lm}$, and $C_{l}$. The masses for $a^{r}_{lm}$ and $a^{i}_{lm}$ are chosen equal to inverse of their variance. Therefore,

[TABLE]

The mass for $C_{l}$ is chosen as inverse of the cosmic-variance of $C_{l}$

[TABLE]

In the implementation of the computation, we find it convenient to use scaled parameters $10^{2}m_{1N}$ to ensure that variable takes values of the order unity. We choose unit masses corresponding to $m_{10},m^{r}_{11},m^{i}_{11}$, and $\alpha$.

To perform the Hamiltonian dynamics, it is necessary to compute derivatives of the conjugate momentum. This involves taking derivative of the distribution function with respect to the corresponding parameter. For the probability distribution given in eq. (3.4), it is possible to analytically compute the derivatives with some approximations, bypassing the need for numerical computation of the derivatives. The detailed expressions for the momentum derivatives for all the parameters are given in Appendix B. The formalism applies to any covariance matrix with non-zero off-diagonal terms. Hence, we discuss the methodology for the general modulation case and note explicitly if any result is specifically applicable only to the scenario of dipole modulation.

We simulate the Hamiltonian dynamics for the Hamiltonian given in eq. (4.1) using above equations. Following the HMC algorithm then leads to the samples of ($a^{r}_{lm},a^{i}_{lm},C_{l},m_{10},m^{r}_{11},m^{i}_{11},\alpha$) which are drawn from the multidimensional probability distribution \mathcal{P}\Big{(}S_{lml^{\prime}m^{\prime}},a_{lm}\Big{|}d_{lm}\Big{)}, where $S_{lml^{\prime}m^{\prime}}$ depends on $(C_{l},m_{10},m^{r}_{11},m^{i}_{11},\alpha)$.

4.2 Incorporation of observational non-idealities

The observed sky temperature contains the true signal marred by several observational non-idealities the dominant effects being the residual foreground contamination, anisotropic noise, and instrumental beam. CMB data is masked to reduce the residual foreground contaminations of other emissions from the galactic plane and extragalactic point sources. Hence we deal with a partial CMB sky appropriately masked to minimize residual foreground contamination for the purpose of data analysis. To deal with the masked sky and anisotropic noise, we resort to the real space variables.

Let T represent the vector of CMB temperature signal in the pixel space and d the data vector in the same space. We assume that noise is Gaussian distributed with variance $\sigma^{2}_{i}$ for the $i^{th}$ pixel and the noise between two pixels are not correlated. Then the probability distribution of T given d is

[TABLE]

where $N_{pix}$ is the total number of the pixels. We model the masked portion by setting the noise variance in those pixels to a large (effectively infinite) value. We take into account the finite resolution of CMB data that is convolved with the beam through a circularly symmetric beam transfer function $b_{l}$. Hence, the time derivative of conjugate momentum corresponding to $a_{lm}$ is obtained as

[TABLE]

where $p_{l}$ is the pixel window function in harmonic space. Note, the anisotropic noise is implicitly incorporated in the above equation due to the dependence of the $\sigma^{2}_{i}$ on pixels. We used map2alm and alm2map routines of HEALPix555http://healpix.sourceforge.net [53] to compute the second term in the above equation. Given the difference map in harmonic space, $(d_{l^{\prime}m^{\prime}}-b_{l}p_{l}a_{l^{\prime}m^{\prime}})$, the alm2map subroutine executes the following computation

[TABLE]

With the map $M_{D}(\hat{n}_{i})$ as input to the subroutine map2alm, we get following quantity as the output

[TABLE]

This section concludes all the essential ingredients to perform an HMC analysis of a non-SI CMB sky in the presence of the anisotropic noise and mask.

5 Results: Analysis of Planck-2015 CMB temperature map

In this section, we provide the result of the analysis done on the CMB temperature anisotropy map provided by Planck collaboration666Planck-2015 temperature data are substantially similar to the Planck-2018 release that focused on large angle polarization anisotropies. We use the SMICA CMB temperature map777Name of the file:COM_CMB_IQU-smica_1024_R2.02_full.fits available at https://pla.esac.esa.int/pla/#maps and the corresponding SMICA mask with sky fraction $84\%$ in this analysis. The masked SMICA map is depicted in the figure 2(a). Since we are only interested in the large angular scale anomalies, we choose to carry out our analysis at the resolution of NISDE = 256 to reduce the computational cost. The estimation of the noise variance in every pixel is obtained from $100$ Full Focal Plane (FFP8.1) noise simulations [54] available at NERSC888 http://crd.lbl.gov/departments/computational-science/c3/c3-research/cosmic-microwave-background/cmb-data-at-nersc/ which we first low pass filter at multipole l = 256 before computing noise variance. The noise variance map obtained in this way and used in our analysis is shown in figure 2(b).

We performed the analysis assuming both the models for dipole modulation: (1) step model and (2) power law model. The analysis with the step model is useful for comparison of our results with the previous studies [4, 5]. Detailed results of the analysis with step model are given in section 5.1. Analysis for the power law model is given in section 5.2. A summary of the results for both these cases is provided in Table 1 and 2. We give the mean values and standard deviation of the parameters $m_{10}$, $m^{r}_{11}$, and $m^{i}_{11}$. Values of $A_{*}$, $l$, and $b$ given in the table are obtained from the mean values of $m_{10}$, $m^{r}_{11}$, and $m^{i}_{11}$, whereas their standard deviations are obtained from their respective Monte Carlo samples.

5.1 Step-model of scale dependence

To compare the results of our method with the published literature [5, 4], we present the results of our analysis with a step model. In the step model the dipole modulation parameters take non-zero, constant values over the multipole range $l=2$ to $l=l_{cut}$. For the results presented in this section we take $l_{cut}=64$, a choice motivated by the literature [5, 4]. Figure 3 shows the posterior distribution of $m_{10}$, $m^{r}_{11}$, and $m^{i}_{11}$ for SMICA map obtained from $10^{5}$ Monte Carlo samples after removing the $5\times 10^{4}$ samples as Burn-In. The inferred joint distribution of $m_{10},m^{r}_{11}$, and $m^{i}_{11}$ parameters do not show any significant correlation among these variables as evident from the joint distributions of these parameters given in figure 3. We compare our estimates of the dipole parameters with the estimates of corresponding parameters obtained using the minimum variance estimation method for BipoSH (BipoSH-MVE) given in [4]. We note that the current estimate of the parameter $m^{r}_{11}$ deviates from the BipoSH-MVE estimate by around $1\sigma$. Other two parameters agree with their corresponding value from BipoSH-MVE. It is evident from the mean ($\mu$) and standard deviation ($\sigma$) of these parameters that each parameter is away from zero more than one $\sigma$. In the analysis, assumption about the multipole range of modulation goes in the form of shape factor. For the range of multipoles considered here $l=2$ to $l=64$, $\sigma=0.039$. The standard deviation for $m^{r}_{11}$ and $m^{i}_{11}$ is $\sigma/\sqrt{2}=0.027$. The standard deviations of sampled distributions are close to the values expected from the analytical arguments.

The combined effect of this deviation from the origin, on the dipole amplitude $A$, is reflected in the posterior distribution of $A$ given in the figure 4(a). The posterior distributions of $A,$ in figure 4(a) and the posterior distribution of $\theta_{p},\phi_{p}$ in figure 4(b) are obtained by transforming the Monte Carlo chains of {$m_{10}$, $m^{r}_{11},m^{i}_{11}$} to {$A,\theta_{p},\phi_{p}$} using coordinate transformation equations given in eq. (2.4). We also quote the estimates of {${A,\theta_{p},\phi_{p}}$} from the previous analysis [4] in figure 4(a) and 4(b). This analysis shows good agreement with all the previous analysis, with the maximum difference (about one $\sigma$) with the BipoSH minimum variance estimator.

5.2 Power-law model of scale dependence

The power-law model of dipole modulation is given by eq. (2.5). We perform the HMC analysis for the power-law model by taking the pivot point as $l_{p}=16$ for NSIDE $=256$ SMICA temperature map. In figure 5, we show the posterior of $m_{10}(l_{p}=16)$, $m^{r}_{11}(l_{p}=16)$ and $m^{i}_{11}(l_{p}=16)$ obtained from $10^{5}$ Monte Carlo samples after the Burn-In period of $5\times 10^{4}$ samples. Figure 5 also shows the posterior of the power law index $\alpha$ for the SMICA map. Similar to the analysis with the step model, the $m_{10}(l_{p}),m^{r}_{11}(l_{p})$ and $m^{i}_{11}(l_{p})$ parameters are treated as independent variables and correlation among these three parameters is not expected. This is evident from the joint distributions of these parameters given in the figure 5. However, the choice of pivot multipole $l_{p}$ can lead to correlation between power law index $\alpha$ and ($m_{10}(l_{p}),m^{r}_{11}(l_{p})$, $m^{i}_{11}(l_{p})$). But our choice of pivot multipole $l_{p}=16$ does not lead to any such significant correlation. The distribution of $A(l_{p}=16)$ given in the figure 6(a) reflects the combined effect of harmonic space parameters being away from the origin. The consistency of our result for a different choice of pivot point ($l_{p}=32$) is shown in Appendix D. The posterior distributions of {$A(l_{p}=16),\theta_{p},\phi_{p}$} are obtained by transforming the Monte Carlo chains of {$m_{10}(l_{p}=16)$, $m^{r}_{11}(l_{p}=16),m^{i}_{11}(l_{p}=16)$} to {$A(l_{p}=16),\theta_{p},\phi_{p}$} using coordinate transformation formulas given in eq. (2.4). Significant nonzero value of amplitude $A(l_{p}=16)$ along with the nonzero value of $\alpha$ points to the presence of the dipole modulation with scale dependence. The recovered power-law profile of the CHA signal is depicted in figure 1 along with the profile for the step-model and the Doppler Boost signal which is of the order $10^{-3}$ [21, 22, 23, 55, 56, 57]. The recovered power-law profile indicates more than $0.5\%$ modulation effect up to the multipole of $l=256$ and hence is a contamination to the measurement of Doppler Boost signal at these low multipoles [23]. The recovered profile can also lead to direction dependence in the cosmological parameters [58, 59]. Any theoretical studies in future to understand the origin of CHA need to satisfy the recovered profile shown in this analysis.

The distribution of the dipole direction for the SMICA map is given in the figure 6(b). The peak of the distribution gives our best-fit estimate of the dipole direction, which is $(l,b)=(247.8^{o},-19.6^{o})$ which is mildly deviant from the direction as reported by the previous analysis [5, 4]. We find that the direction of the dipole modulation under step model and power law model are in agreement and hence does not depend on the profile of the dipole modulation. The direction estimated under step model (see figure 4(b)) $l,b=(238.0^{o},-22.4^{o})$ is consistent with that estimated under power law model $l,b=(247.8^{o},-19.6^{o})$ (see figure 6(b)). A more quantitative comparison between the direction recovered in our analysis and the direction reported in the literature is given in figure 7. The estimates of the dipole direction provided by Pixel Based Likelihood (PBL) estimate and BipoSH minimum variance estimate (BipoSH-MVE) are consistent within the 68% confidence level, whereas Quadratic Maximum Likelihood (QML) estimate of the dipole direction lies outside the 68% contour. In particular, we note that the South Ecliptic Pole (SEP) falls within the boundary of the 68% confidence region.

In remaining part of this section, we discuss the probability distribution of $C_{l}$ obtained using HMC analysis of the SMICA map. In figure 8 we show the posterior distribution for $C_{l}$ jointly estimated with the dipole modulation parameters for SMICA map at select multipoles. We fit the analytical probability distribution of $C_{l}$ given in [60] to the histogram of $C_{l}$ samples. We use SciPy routine scipy.optimise.curve_fit to implement the fitting [61]. After fitting the probability distribution of $C_{l}$ to the histogram of the samples of $C_{l}$, we get the best-fit estimate of the angular power spectrum ($C^{HMC}_{l}$), which matches very well with the quadratic estimation of the angular power spectrum obtained from the map. In figure 9 we show two estimates of the angular power spectrum: (1) angular power spectrum estimated jointly with the dipole modulation parameters ($C^{DM}_{l}$) and (2) angular power spectrum estimated without the dipole modulation parameters ($C_{l}$). Both the estimates are compared with the best-fit $\Lambda$CDM angular power spectrum provided by Planck [62]. We note that irrespective of whether we estimate $C_{l}$ with or without dipole modulation parameters, the estimates of $C_{l}$ do not show any significant differences. The suppression in the power spectrum in the multipole range $l\in\{20,30\}$ remains at the same statistical significance even in the joint analysis of $A(l)$ and $C_{l}$. The estimated CHA signal also remains unaltered when the $C_{l}$ are kept fixed at the fiducial $\Lambda$CDM values. The corresponding results are shown in Appendix E.

The recovered power-law of the dipole modulation signal depicted in figure 1 shows a strong modulation field in the low $l$. As a result, the amplitude of the second order terms $A^{2}(l)$ gets stronger at the low $l$ (as mentioned in eq. (2.7)). We estimate the correction to the angular power spectrum due to second order term and depict the relative correction to the angular power spectrum $\Delta_{l}$ in figure 10. For comparison, we also show the $1\sigma$ cosmic variance error-bar by the shaded region. This establishes that the second order term is negligible and does not play any significant role in the analysis.

6 Model comparison between SI and non-SI models of CMB map: Estimation of Bayes factor

Bayesian model comparison:

Bayesian analysis allows comparison of two models (say $M_{0}$ and $M_{1}$). One compares the probability of the model $M_{0}$ given data ($\mathcal{P}(M_{0}|d)$) and the probability of the model $M_{1}$ given data ($\mathcal{P}(M_{1}|d)$) using the following ratio called Bayes factor

[TABLE]

In the above equation, $\mathcal{E}(d|M)$ stands for the probability of the data given model $M$ and $\mathcal{E}(d|M)$ is called Bayesian Evidence for $M$. $\mathcal{P}(M)$ stands for the prior probability of the model $M$. In the situation where two models are considered a priori equally probable ($\mathcal{P}(M_{0})=\mathcal{P}(M_{1})$), the Bayes factor is the ratio of the evidence of the models under the same data. If the Bayes factor is close to unity, the inference is that the data do not discriminate between the two models.

In our work, we compare the SI model ($M_{SI}$) of CMB with the power-law dipole modulation model ($M_{DM}$). The Bayes factor in favor of the SI model, $B_{SI-DM}$ is

[TABLE]

where we have assumed that prior probability $\mathcal{P}(M_{SI})=\mathcal{P}(M_{DM})$. In above equation, $\mathcal{E}(d|M_{SI})$ and $\mathcal{E}(d|M_{DM})$ are evidence under $M_{SI}$ and $M_{DM}$, respectively, which is defined as

[TABLE]

where, $w_{x}$, $w_{y}$, and $w_{z}$ are convenient re-parameterizations of $m^{r}_{11}$, $m^{i}_{11}$, and $m_{10}$ respectively, as defined in eq. (3.10). $\mathcal{L}(w_{x},w_{y},w_{z},\alpha)$ denotes the likelihood and $\Pi(w_{x},w_{y},w_{z},\alpha|M)$ is the prior on parameters.

Bayes factor using Savage-Dickey density ratio:

The power-law dipole modulation model has four additional parameters compared to the SI model, which are $w_{x}(l_{p}),w_{y}(l_{p}),w_{z}(l_{p})$, and $\alpha$. The model $M_{SI}$ is nested within the model $M_{DM}$ in the sense that using $(w_{x}(l_{p}),w_{y}(l_{p}),w_{z}(l_{p}))=(0,0,0)$ reduces $M_{DM}$ to $M_{SI}$ for all values of $\alpha$. The Bayes factor for nested models can be obtained by using Savage-Dickey density ratio (SDDR) [63, 64]. For SI model nested in the power law dipole modulation model, the SDDR is

[TABLE]

where $\mathcal{P}_{C_{l}}(w_{x}=0,w_{y}=0,w_{z}=0,\alpha|d,DM)$ is posterior marginalized over $C_{l}$. $\alpha$ dependence of SDDR in eq. (6.4) is weak because the likelihood conditioned at $(w_{x},w_{y},w_{z})=(0,0,0)$ is independent of $\alpha$. The prior mentioned in (6.4) can be written as

[TABLE]

We choose uniform prior on $\alpha$ over the range $\alpha=-2\text{ to }0$. Hence, the normalized prior density of $\alpha$ is

[TABLE]

The amplitude of the dipole modulation $A(l)$, decides the relative magnitude of the off-diagonal terms in the covariance matrix compared to the diagonal terms ($C_{l}$). In order to ensure the positive definiteness of the covariance matrix

[TABLE]

In spherical polar coordinates, with $r$ as defined in eq. (3.11), the above argument gives an upper bound: $\mathbf{r\leq\sqrt{(\pi/3)}\approx 1}$, which is denoted by $R$. To obey eq. (6.7), we need to make the upper bound on $r$ dependent on $\alpha$. Otherwise, for low values of $\alpha$ close to $-2$, $A(l)$ can become greater than $1/2$ at low multipoles for sufficiently high $A(l_{p})$. Hence, $\alpha$ dependent upper bound $R(\alpha)$ is

[TABLE]

At $\alpha=-2$, with $l_{p}=16$ and $l=2$, $R(\alpha)=0.016$.

We choose a prior on $(w_{x},w_{y},w_{z})$ which has uniform density within the sphere of radius $R(\alpha)$, centred at $(w_{x},w_{y},w_{z})=(0,0,0)$ and zero outside. Normalized form of such prior is

[TABLE]

Note that the prior density is independent of the direction. The form of prior in eq. (6.9) ensures positive definiteness of the covariance matrix at all scales for the range of $\alpha$ considered in the analysis. Choosing an uniform density prior as given in eq. (6.9) implies that the angle marginalized prior density of the amplitude is

[TABLE]

As mentioned in section 3.1, we have sampled the likelihood distribution without augmenting it with any prior. We obtain the posterior by weighting the likelihood samples by the prior density given in eq. (6.9). With this procedure, the SDDR in favour of the SI model is found to be

[TABLE]

In Appendix F, we obtain the SDDR for two more choices of prior densities and show that the variations are negligible.

This implies that the data do not favor the concordance model of cosmology over the non-SI model. The large prior range for the non-SI parameters {${A(l),\alpha}$} is the key reason for the inconclusive Bayes factor, though the likelihood clearly peaks away from zero (the expected value under the concordance model of cosmology) as shown in figure 6(a). In summary, the best available nearly full sky CMB temperature data does not rule out non-SI cosmological model even after carefully including the observational non-idealities (masking, anisotropic noise, instrument beam). Other cosmological probes like CMB polarization, weak lensing, and large scale structure surveys keep the hope alive to provide conclusive evidence for or against this enigmatic anomaly in the future.

7 Discussions and Conclusions

In this paper, we study the combination of two large angular scale CMB anomalies, namely the power suppression and dipole power asymmetry (also called the Cosmic Hemispherical Asymmetry (CHA)) in the Bayesian framework using the Planck-2015 SMICA temperature map. Most of the analysis in the literature assume a step model which have a non-zero amplitude of the dipole signal only below an arbitrarily chosen value of multipole ($l_{cut}$). This value of the $l_{cut}$ is usually chosen as $64$, though a few analysis are also performed setting $l_{cut}=600$ [65]. We model the amplitude of the modulation dipole with a two-parameter power-law, the amplitude at the pivot multipole ($A(l_{p})$) and the spectral index ($\alpha$). The parameters $\theta_{p}$ and $\phi_{p}$ give the direction of the modulation dipole. The joint estimate of these four parameters in harmonic space together with the temperature coefficients of the map $a_{lm}$ and power spectrum $C_{l}$ are performed after considering the observational non-idealities like the partial sky, instrument beam, and anisotropic instrument noise.

The jointly estimated angular power spectrum $C_{l}$ matches well with the value of $C_{l}$ reported by Planck [62]. The low value of the quadrupole, suppression of power between the multipole range $20-30$ holds at the same significance. Hence, we reconfirm the existence of power suppression from this joint study even within the more general statistical model including an anisotropic covariance matrix.

Our estimate reveals a scale-dependent CHA signal with the maximum a posteriori values of the parameters as $A_{*}(l_{p}=16)=0.064\pm 0.022$ and $\alpha=-0.92\pm 0.22$ for a power law parametrization. The direction of the maximum signal in the galactic coordinates is $(l,b)=(247.8^{o},-19.6^{o})$. The spectral index of the power-law differs by about one $\sigma$ from the previously reported value by Aiola et al. [24], when compared with their result for $l_{max}=300$. However, the amplitude at the pivot multipole $l_{p}=16$ ($A(l_{p}=16)=0.071$) is consistent within one sigma of our result. These mild differences can arise due to the choice of different pivot points, the maximum value of the multipoles used ($l_{max}$), use of the different Planck noise simulations (FFP6 instead of FFP8) and most importantly due to the difference in the analysis framework. The direction of the asymmetry is also recovered consistently in both the analysis. We note that, the best fit direction of the dipole is almost the same for both step-model and power-law model analysis. This is expected because both models differ only in the treatment of the dipole amplitude.

We also make a model comparison between power law form of the scale-dependent dipole modulation model and the concordance SI model of CMB. The estimated Bayes factor in favor of the concordance model is $0.4$ and does not rule strongly in favor of SI model, hence remains inconclusive. The current best available temperature data does not stand out loud in favor of the concordance model and shows only a mild inclination in favor of the non-SI model. Future datasets from CMB-polarization, weak lensing and large scale structure are capable to shed light to this enigmatic anomaly.

Though the two-parameter power-law model is simplistic, it has the ability to capture the scaling of scale dependence. A future study with the additional parameters beyond the power-law form can be performed with the availability of the CMB polarization and large scale structure datasets. The extrapolation of the retrieved power-law form to small angular scales leads to a CHA signal less than $0.1\%$ for angular scales below $10$-arcmin. This is sub-dominant (as required) from the well-predicted Doppler Boost signal [55, 56, 57] due to our local motion with a velocity $\beta\equiv v/c=1.23\times 10^{-3}$ [21, 22], which is also measured by Planck using the high resolution temperature data [23]. Hence, the retrieved signal profile obtained in this paper satisfies the observational requirement at both small and large angular scales. The profile of the power-law signal obtained in this analysis will be useful guide in any pursuit of a theoretical understanding of CHA signal and also to make predictions for other cosmological probes including CMB polarization [66, 67, 68, 69, 70, 71]. Bayesian inference of the BipoSH coefficients using HMC discussed in this work can be generalized to CMB polarization. The systematic noise due to the instrument and foreground effects become more important for low multipole polarization data provided by Planck and including these effects accurately requires an in-depth analysis that we defer to future work.

Acknowledgements: S.S. acknowledges University Grants Commission (UGC), India for providing the financial support as Senior Research Fellow. S. M. and B. D. W. acknowledge the support of Simons Foundation for this work. This work is also supported by the Labex ILP (reference ANR-10-LABX-63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the programme Investissements d’avenir under the reference ANR-11-IDEX-0004-02. Authors would like to thank Planck collaboration and NERSC for providing the FFP simulations and other auxiliary files. The present work is carried out using the High Performance Computing facility at IUCAA 999http://hpc.iucaa.in. Some of the results in this paper have been derived using the HEALPix package [53]. We also acknowledge the use of CAMB [72], IPython [73], Python packages Matplotlib [74], NumPy [75], SciPy[61] and corner [76].

Appendix A Derivation of $\mathcal{P}(C_{l},m_{1M}|\textbf{d})$

This appendix provides a detailed derivation of eq. (3.5) and eq. (3.6). To derive eq. (3.5), we marginalize over a from eq. (3.4). We simplify the following expression in eq. (3.4)

[TABLE]

by completing the square [77].

[TABLE]

where $\tilde{\textbf{a}}=\textbf{S}(\textbf{S}+\textbf{N})^{-1}\textbf{d}$ is Wiener filter applied to the data. eq. (3.4) now becomes

[TABLE]

Using the results

[TABLE]

and

[TABLE]

in eq. (3.4) to get $P(\textbf{S}|\textbf{d})$ as given in eq. (3.5). Now we derive the joint probability distribution $C_{l}$ and $m_{10}$ given in eq. (3.6) from eq. (3.5). We write the sum of signal covariance matrix S and noise covariance matrix N in the following manner

[TABLE]

Since we assume N to be diagonal, it contributes only to D. $\textbf{O}_{1}$ is the off-diagonal part of S without $m_{10}$. Taylor series expansion of inverse of the sum $\textbf{S}+\textbf{N}$, upto second order in $m_{10}$ gives

[TABLE]

Determinant $|\textbf{S}+\textbf{N}|$ in eq. (3.5) also has the $m_{10}$ dependence expressed as

[TABLE]

Making use of the matrix identity $\ln|M|=Tr[\ln(M)]$ to compute the logarithm of the determinant, we get

[TABLE]

Taylor series expansion of logarithm gives

[TABLE]

Then, keeping the terms up to second order in $m_{10}$, we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We substitute eq. (A.6) and eq. (A.13) in the following form of eq. (3.5)

[TABLE]

to obtain

[TABLE]

eq. (A.15) can be further rearranged in following manner

[TABLE]

where $M$ represents following quantity

[TABLE]

eq. (3.6) follows from eq. (A.16) given above.

Appendix B Position and Momentum derivatives

This appendix contains the discussion of the position and momentum derivatives required for the Hamiltonian dynamics in HMC.

Following are the Hamilton’s equations for $(a^{r}_{lm},a^{i}_{lm})$

[TABLE]

This involves time derivative of conjugate momentum corresponding to $a_{lm}$

[TABLE]

After taking derivative with respect to $a_{lm}$, above equation reduces to

[TABLE]

To compute $\dot{p^{r}}_{lm}$ and $\dot{p^{i}}_{lm}$, the momentum derivatives corresponding to the real and imaginary part of $a_{lm}$ respectively, we make use of the following relations

[TABLE]

where $\Re$ and $\Im$ denote the real and imaginary part of the complex number, respectively. Hamilton’s equations for $(m_{10},m^{r}_{11},m^{i}_{11})$ are

[TABLE]

[TABLE]

where $\frac{\partial\mathcal{H}}{\partial m_{1N}}$ are given in terms of derivative of $\mathcal{H}$ with respect to $A^{1N}_{ll+1}$ as

[TABLE]

The factor of 2 in above expression appears because $A^{1N}_{ll+1}=A^{1N}_{l+1l}$ for the dipole modulation model. The time derivative of the momentum corresponding to $\alpha$ is obtained as

[TABLE]

From eq. (B.7) and eq. (B.8), we note that to compute derivatives of momentum corresponding to $m_{10},m^{r}_{11},m^{i}_{11}$ and $\alpha$, we need to know $\frac{\partial\mathcal{H}}{\partial A^{1M}_{ll+1}}$.

The time derivative of conjugate momentum corresponding to the BipoSH coefficient $A^{LN}_{ll^{\prime}}$ is

[TABLE]

Following is a detailed derivation of the expression for $\dot{p}^{LN}_{ll^{\prime}}$ starting from eq. (B.9). To explicitly evaluate eq. (B.3) and to simplify eq. (B.9) we use a Taylor series expansion of $\textbf{S}^{-1}$ up to first order in O as

[TABLE]

We make use of two standard results from matrix calculus

[TABLE]

and

[TABLE]

The derivative of element of the covariance matrix with respect to BipoSH coefficient is

[TABLE]

Using the Taylor series expansion, the element of the inverse of the covariance matrix is

[TABLE]

Using the aboe Taylor series expression and the derivative in eq. (B.13), eq. (B.11) can be expressed as

[TABLE]

Using the derivative in eq. (B.13), eq. (B.12) becomes

[TABLE]

The term in the above equation with the Taylor series approximation gives

[TABLE]

In our computation, we replace the second and third term by its ensemble average and keep terms only upto first order in O.

[TABLE]

Eq. (B.16) along with above simplification and eq. (B.15) leads to following expression for $\dot{p}^{LN}_{ll^{\prime}}$

[TABLE]

The time derivative of the conjugate momentum of $C_{l}$ is obtained from eq. (B.19) by substituting $L=0$ and $M=0$ in the above equation (as $A^{00}_{ll}$ are related to the $C_{l}$)

[TABLE]

By using relation between $A^{00}_{ll}$ and $C_{l}$, we get

[TABLE]

Equation B.21 is same as equation (25) in [47]. Note that, in eq. (B.20) and eq. (B.21), $A^{00}_{ll^{\prime}}$ and $C_{l}$ are the variables of dynamics, whereas $\hat{A}^{00}_{ll^{\prime}}$ and $\hat{C}_{l}$ are quantities obtained using other variables of dynamics, namely {$a_{lm}$}.

Appendix C Demonstration of method on simulated map using CoNIGS

We test our algorithm on the simulated maps with the power law dipole modulation signal, in presence of anisotropic noise and masking. The power law dipole modulated CMB map is generated using CoNIGS algorithm [33]. We add the anisotropic noise [54] to this simulated map and do the analysis in presence of SMICA mask. We take $\{m_{10}(l_{p}=32),m^{r}_{11}(l_{p}=32),m^{i}_{11}(l_{p}=32)\}=\{-0.018,0.027,-0.030\}$ and $\alpha=-0.7$. The values of $\{m_{10}(l_{p}),m^{r}_{11}(l_{p}),m^{i}_{11}(l_{p})\}$ are chosen so as to mimic the direction of the observed CHA based on the BipoSH-MVE estimate for the SMICA map [4]. We do the analysis over the multipole range $l=2-256$. In this way we confirm that we recover the signal consistent with that injected in the simulated map. Results of such an HMC run on one map are presented in figure 11, figure 12(a), figure 12(b) and figure 13.

Histograms in figure 11 represent the distributions of harmonic parameters of dipole, $m_{10}(l_{p}),m^{r}_{11}(l_{p})$ and $m^{i}_{11}(l_{p})$ at the pivot multipole $l_{p}=32$, and the power law index $\alpha$. These distributions are obtained using Monte Carlo chain of $10^{5}$ samples after discarding $5\times 10^{4}$ samples as Burn-In. The set of 1D histograms, marginalized distributions of the parameters, in figure 11 show that the sampled distribution is consistent with the input values of the respective parameters for the map. We test the algorithm on different map realizations. We find that the mean values of $m_{10},m^{r}_{11}$ and $m^{i}_{11}$ samples vary fairly randomly around the values assumed in simulation. This suggests that our estimates do not have any bias.

From eq. (3.7) for the variance of $m_{10}$, we know that the variance depends on the range of multipoles over which the dipole modulation signal exists and the power law index. This information is not available prior to the analysis of real data. For the simulated map, $\alpha=-0.7$ and $l=2$ to $l=256$. This leads to the following: $\sigma_{m_{10}}=0.023$ and $\sigma_{m^{r}_{11}}=\sigma_{m^{i}_{11}}=0.016$. The corresponding figures given in figure 11 are in agreement with these estimates. The standard deviations of sampled distributions are close to the values expected from the analytical arguments. Using $\{m_{10}(l_{p}=32),m^{r}_{11}(l_{p}=32),m^{i}_{11}(l_{p}=32)\}$ we obtain the samples of $A(l_{p}=32),\theta_{p},\phi_{p}$ using eq. (2.4). Figure 12(a) shows the distribution of $A(l_{p}=32)$ obtained from the harmonic space variables shown in figure 11. Figure 12(a) indicates that our method does recover the input value present in the simulated map. Figure 12(b) shows the distribution of $\theta_{p},\phi_{p}$ samples. We see that our best-fit direction is close to the dipole modulation direction introduced in the simulated map.

A significant fraction of our computation is expended in incorporating the estimation of the map angular power spectrum. A summary statistics of $C_{l}$ distribution is given in figure 13. Using the Monte Carlo chain of $C_{l}$ samples, we obtain the maximum a posteriori estimate of $C_{l}$ $(C^{HMC}_{l})$ from our samples. We fit the analytical probability distribution of $C_{l}$ given in [60] to the histogram of $C_{l}$ samples. We use SciPy routine scipy.optimise.curve_fit to do the fitting [61]. In figure 13 we compare $C^{HMC}_{l}$ with the angular power spectrum of the map realization, $(C^{map}_{l})$. It is clear that we faithfully recover the estimate of angular power spectrum using our algorithm.

Appendix D Establishing robustness to the choice of pivot multipole

For the power law model analysis of SMICA map presented in the main body of the article, we chose $l=16$ as the pivot multipole. To corroborate the answer we get for $l_{p}=16$, in this appendix we provide the result for a different choice of pivot multipole, $l_{p}=32$. The distributions of parameters $m_{10}(l_{p}=32),m^{r}_{11}(l_{p}=32),m^{i}_{11}(l_{p}=32)$ and $\alpha$ are given in the figure 14. Taking $\alpha=-0.92$ and extrapolating from the results given in figure 5 for $l_{p}=16$, we expect $m_{10}(l=32)=-0.023,m^{r}_{11}(l=32)=0.017,m^{i}_{11}(l=32)=-0.043$. These values are consistent with the estimates of corresponding parameters given in figure 14.

Appendix E HMC analysis of CHA with fixed fiducial $\Lambda$CDM $C_{l}$

From eq. (2.13) and eq. (2.14) we see that the nature of $C_{l}$ affects the magnitude of $m_{1N}$. In the main body of this article we performed the joint analysis of $C_{l}$ and $m_{1N}$ and find that we recover the $C_{l}$ consistent with the one reported in the literature [62]. We also estimate the correction to $C_{l}$ due to second order term in the dipole modulation covariance matrix, which turns out to relatively small except at low multipole $l<5$, below which correction is around 1%. In this section, we provide the result of the sampling of $m_{10}(l_{p}),m^{r}_{11}(l_{p}),m^{i}_{11}(l_{p}),\alpha$ where the angular power spectrum $C_{l}$ is held fixed at the best-fit $\Lambda$CDM power spectrum provided by Planck [62]. The distribution of $A(l_{p})$ and that of $(\theta,\phi)$ given in figure 15(a) and figure 15(b) respectively are similar to those given in figure 6(a) and figure 6(b) with $C_{l}$ as variable.

Appendix F Establishing robustness of the evidence ratio against prior choice

In section 6, we discussed a prior, say, Prior-1, which is spherically symmetric and has a constant density within the sphere of radius $R(\alpha)$. The angle marginalized prior density of the amplitude for Prior-1 goes as $r^{2}$. Here, we explore two other priors, Prior-2 for which the angle marginalized amplitude prior density goes as $r$ and Prior-3 for which angle marginalized amplitude prior density is constant with respect to $r$.

Prior-2 is spherically symmetric, centered at $(0,0,0)$ and is inversely proportional to the separation from $(0,0,0)$ till $R(\alpha)$ and zero afterward. The normalized form of such a prior is

[TABLE]

With Prior-2, the angle marginalized prior density of the amplitude is

[TABLE]

Prior-3 is spherically symmetric, centered at $(0,0,0)$ and is inversely proportional to the square of separation from $(0,0,0)$ till $R(\alpha)$. Prior-3 is again zero outside the sphere of radius $R(\alpha)$. The normalized form of such a prior is

[TABLE]

The prior density $\Pi_{3}$ has the property that it gives uniform probability density to the amplitude of the modulation dipole because the angle marginalized prior density on the amplitude is

[TABLE]

The prior densities $\Pi_{2}$ and $\Pi_{3}$ diverge as $(w_{x},w_{y},w_{z})\rightarrow(0,0,0)$. Hence we can not use the formula for the SDDR given in eq. (6.4) directly to evaluate the Bayes factor. In practice, the numerator and denominator of the SDDR are computed by estimating the densities in a small but finite volume around the nested point. In this section we present an attempt to evaluate the posterior to prior density ratio at a point $(x,y,z)=(\epsilon,\epsilon,\epsilon)$ close to $(x,y,z)=(0,0,0)$, where $\epsilon$ is small enough that it cannot be distinguished from zero at the level of the noise. This essentially amounts to cutting off the divergence in the prior at $r=\epsilon$. Modifying eq. (6.4), the Bayes factor in favour of the SI model is

[TABLE]

For Prior-1, the above expression is evaluated at $(w_{x},w_{y},w_{z})=(0,0,0)$ and is exactly the SDDR in favor of SI model. We estimate the posterior density $\mathcal{P}_{C_{l}}(w_{x}=\epsilon,w_{y}=\epsilon,w_{z}=\epsilon,\alpha|d,DM)$, which appears at the numerator of eq. (F.5), using the Monte-Carlo samples and need a finite bin-width (around $3\times 10^{-2}$) to get the estimate of posterior at any point. As a result of this binning, the posterior estimate does not change over the scale of the bin and is independent of $\epsilon$ as long as $(w_{x},w_{y},w_{z})=(0,0,0)$ and $(w_{x},w_{y},w_{z})=(\epsilon,\epsilon,\epsilon)$ are within the same bin. We obtain the ratio of posterior density and prior density in the grid containing the point $(x,y,z)=(\epsilon,\epsilon,\epsilon)$ with $\epsilon=10^{-3}$. SDDR thus obtained for Prior-2 and Prior-3 is given in Table  3.

From eq. (F.5), SDDR depends on $\alpha$. However, at the nested point, which reduces the power law dipole modulation model to SI model, the likelihood does not depend on $\alpha$ and the SDDR also remains roughly constant over the range of $\alpha$. In Table  3 we quote the value of SDDR averaged over $\alpha$.