Joint Bayesian analysis of large angular scale CMB temperature anomalies
Shabbir Shaikh Suvodip Mukherjee Santanu Das Benjamin D. Wandelt Tarun Souradeep


Ground-state energy estimation of the water molecule on a trapped ion quantum computer
    Jungsang Kim
    Received ... 2018/ Accepted ?
======================================================================================





§ 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 <cit.>. 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 <cit.>. Prominent among the CMB anomalies is the Power Asymmetry or the Cosmic Hemispherical  Asymmetry (CHA), which is characterized by a ≈ 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 <cit.>. CHA was later confirmed to also exist in the Planck CMB map at similar amplitude as seen in the WMAP data <cit.>. 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) <cit.>. Various physical models have been put forth 
for the plausible origin of CHA in the CMB, see for example <cit.>. 

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 (ΛCDM) model with power law form of the power spectrum of initial density fluctuations <cit.>. An existing body of literature also explores connections between different anomalies. For example, Muir et al. (2018) <cit.> investigate the covariance between various anomalies. Polastri et al. (2015) <cit.> study the connection between directional anomaly of low multipole alignment and dipole modulation. Schwarz et al. (2016) <cit.> 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 <cit.> to scrutinize both these anomalies jointly within a common mathematical framework of Bipolar Spherical Harmonics (BipoSH). BipoSH <cit.> 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 <cit.> in the BipoSH framework by Das et al. <cit.> to incorporate essential observational non-idealities in the Planck data 
[While this work was in progress, one of the coauthors (Santanu Das) has posted two arXiv preprints <cit.> (PDF file of one of those was withdrawn later <cit.>) 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 <cit.>. The method and the code used in this analysis are developed independently from the initial code of <cit.> for full sky and isotropic noise.].
We test our algorithm on simulated maps generated using  algorithm <cit.> 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 <ref>. The details of the analysis and its significance are discussed in section <ref>. 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 <cit.> as shown by the blue curve in figure <ref>. For comparison, we also plot the scale-independent signal due to our local motion <cit.> 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 <ref> 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 <ref> 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 <ref>. In appendix <ref> we demonstrate the application of the method to simulated maps. Section <ref> and <ref> 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 <ref>.



§ COSMIC HEMISPHERICAL ASYMMETRY



 §.§ Introduction to Dipole Modulation Model

CHA is modeled as dipole modulation of statistically isotropic (SI) CMB 
sky <cit.>. Dipole modulated CMB temperature anisotropy field 
(T(n̂) ≡ T(n̂)/T_0-1)[T_0=2.7255 K is the all sky average temperature of CMB <cit.>.]
in the direction n̂, is given by 

    T(n̂) = (1 + Ap̂·n̂)  T^SI(n̂).

Here, T^SI(n̂) is SI CMB temperature anisotropy field. The direction of the modulation dipole is given by p̂≡ (θ_p, ϕ_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 
<cit.>. 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 <cit.>. For analysis based on BipoSH representation, the contribution of CHA is captured in L = 1 BipoSH coefficients 
<cit.>. 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 (θ_p, ϕ_p). The likelihood is sampled using Markov 
Chain Monte Carlo method to estimate the parameters of the modulation dipole 
<cit.>. Results of these different estimators are also included in the 
Planck-2013 <cit.> and Planck-2015 <cit.> analysis. Different methods have detected the signal of CHA at a significance level around 3σ 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 <cit.>. 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, T_tilde can not explicitly express the 
scale dependence of the dipole amplitude. It is best expressed in harmonic space using following equation

    a_lm= a^SI_lm + ∑^1_N = -1 m_1N(l) ∫ T^SI() Y_L = 1, N
    (n̂) Y^*_lm(n̂) d^2Ω_n̂,

where a^SI_lm and a_lm are the spherical harmonic coefficients of T^SI() and T(), 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_T_tilde_l, it is clear that the SI temperature field Δ 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 p̂ and the amplitude A of the dipole. 
The scale-dependent modulation field using real space variables amplitude A(l) and the dipole direction p̂ is expressed as

    A(l) p̂·n̂ = ∑^1_N = -1 m_1N(l) Y_L=1,N(n̂).

The  amplitude A and the direction (θ_p, ϕ_p) in terms of m_10, m^r_11 and m^i_11 are given by following equations

    A(l) = √(3/4 π)√(m_10(l)^2 + 2m^r_11(l)^2 + 2m^i_11(l)^2),   θ_p = cos^-1 [m_10(l)/A(l)√(3/4 π) ],   ϕ_p = -tan^-1[ m^i_11(l)/m^r_11(l)  ].

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)

    A(l) = A(l_p)( l/l_p)^α,

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

    m_1N(l) = m_1N(l_p) ( l/l_p)^α,

where m_1N(l_p) are the spherical harmonic coefficients of the dipole A(l_p)·, that are related to each other by m1M_to_Athephi. 


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 α under the Bayesian framework using HMC method 
<cit.>. Previously Aiola et al. <cit.> studied a power-law scale dependence form to estimate only the CHA signal. However, due to the structural form of eq_T_tilde_l, 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.



 §.§ Covariance matrix of dipole modulated CMB temperature sky

Since the CMB anisotropy is predicted and observed to be consistent with a Gaussian random field <cit.>, its covariance matrix () 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_T_tilde_l), the covariance matrix of a_lm contains non-zero off-diagonal terms.  
The exact[To all orders in m_1N] covariance matrix of a_lm
for the dipole modulation case can be written as

    S_l_1m_1l_2m_2   ≡   ⟨ a_l_1 m_1 a^*_l_2 m_2⟩ = 
    C_l_1δ_l_1 l_2δ_m_1 m_2
       +   Π_l_1l_2/√( 12π)∑^1_N = -1 (m_1N(l_1) C_l_1 (-1)^l_1 + l_2 + 1 
    + m_1N(l_2) C_l_2 ) C^1 0_l_1 0 l_2 0 C^1 N_l_1 m_1 l_2 m_2
       +   Π_l_1l_2/12 π∑_l'_1 m'_1 C_l'_1Π^2_l'_1 
    C^1 0_l_1 0 l'_1 0 C^1 0_l_2 0 l'_1 0∑^1_N_1,N_2 = -1
    m_1N_1(l'_1) m^*_1N_2(l'_1) C^1 N_1_l_1 m_1 l'_1 m'_1 C^1 N_2_l_2 m_2 l'_1 m'_1,

where C_l is angular power spectrum of a^SI_lm. Here, Π_l_1l_2 l_n denotes √((2l_1 + 1)(2l_2+1) (2l_n+1)) and C^1N_l_1m_1l_2m_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_1m_1l_2m_2, in terms of BipoSH coefficients can be expressed as <cit.>

    S_l_1m_1l_2m_2 = (-1)^m_2∑_LN A^LN_l_1l_2 C^LN_l_1m_1l_2-m_2.

The BipoSH coefficients A^LN_l_1l_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 

    A^00_ll = (-1)^l √(2l + 1) C_l.



The leading order term of m_1N in cov_mat 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 cov_mat 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

    ⟨ a_l m a^*_l m⟩ = C_l + Π^2_l/12 π∑_l' m' C_l'Π^2_l' [C^1 0_l 0 l' 0]^2
    ∑^1_N_1,N_2 = -1 m_1N_1(l') m^*_1N_2(l') C^1 N_1_l m l' m' C^1 N_2_l m l' m'.

Note that ⟨ a_l m a^*_l m⟩ 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) ≈ m_1N(l)), we can simplify cov_mat_corr_to_diag as

    C̃_l = C_l + C_l/4π∑^1_N = -1 |m_1N|^2,

where C̃_l ≡⟨ a_l m a^*_l m⟩. Using dipole, the relative change in C_l in terms of A(l) can be written as

    Δ_l ≡C̃_l - C_l/C_l = 1/4π∑^1_N = -1 |m_1N|^2 = A^2(l)/3.

In our analysis, we neglect the second order term in the covariance matrix, but provide an estimate of the relative correction in Sec. <ref>. 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. <ref> to be

    A^1N_ll+1(l_p, α) = m_1N(l_p) G^1_ll+1,

where G^1_ll+1 is referred to as the shape factor of the specific non-SI signal. Specifically, here

    G^1_ll+1≡Π_ll+1/√(12π)[ ( l/l_p)^α C_l + 
     ( l+1/l_p)^α C_l+1]  
    C^10_l0l+10.

DM_biposh when substituted in cov_in_term_of_biposh provides the off-diagonal terms of the covariance matrix in terms of m_1N(l_p) and G^1_l l+1. Scale dependence is an important aspect of this empirical modeling of the CHA that is captured by the index α of the power law.
The shape factor encodes the signature of the dipole modulation through the mixing term between l and l± 1 modes, which is sufficient to measure the non-SI signal from the map.



§ PROBABILITY DISTRIBUTIONS

In this section, we define the probability distributions of variables involved in the problem. We denote[Bold fonts are used to denote matrices.] the diagonal of  by  and the off-diagonal part by  (so that =  +). We express the set of a_lm by the vector 𝐚. For the dipole modulation signal, cov_in_term_of_biposh and DM_biposh jointly express  in terms of C_l, m_1N(l_p) and α. 



 §.§ Joint Probability distribution of Lg and 
Lg 

We assume that the noise in the measurement of CMB temperature is Gaussian distributed with the 
noise covariance matrix . The probability of the data  given ≡{a_lm} is

    (d|a) = 1/√(|N|) (2 π)^n/2exp[ -1/2 
    (d - a)^†N^-1 (d - a)  ].

Here, N^-1 is the inverse of the noise covariance matrix  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

    (a|S) = 1/√(|S|) (2 π)^n/2exp[ -1/2a^†S^-1a ],

where ^-1 denotes the inverse of the covariance matrix .
Using Bayes theorem, the joint probability distribution of  and  given  in terms of 
above two probability distributions is

    (S, a | d) = (d | a) 
    (a|S) (S)/(d).

With uniform prior (S) and up to a normalization constant (d),

    (,  | ) = 1/√(|| || ) (2 π)^n exp{ -1/2[ ( - )^†^-1  ( - ) + ^†^-1 ] }.

In above expression, the dependence of (,  | ) on a_lm is explicit. (, 
 | ) depends on angular power spectrum C_l, m_1N(l_p) and α through . The angular power spectrum, C_l, is the diagonal part of . Dependence of  on m_1N(l_p) 
and α is given by cov_mat.



 §.§ 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 <cit.>). 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  from Prob_a_S and obtain the joint probability distribution of 
C_l and m_1N. Marginalizing  from Prob_a_S yields 

    P(|) = 1/√(| + |) (2π)^n/2exp[-1/2^T ( + )^-1 ] .

The explicit expression for the joint probability distribution of C_l and m_10 can be written as

    P(C_l, m_10| ) = 1/√(||) (2 π)^n/2exp[-1/2^†^-1 ]exp[-(m_10 - μ)^2/2 σ^2]exp[μ^2/2 σ^2],

where

    μ = ^†^-1_1^-1/ Tr[(^-1_1)^2]M^-1  and  σ^2 = 2M^-1/ Tr[(^-1_1)^2],

with the quantity M given by

    M = [ 2 ^† (^-1_1)^2 ^-1/Tr[(^-1_1)^2] - 1 ].

The detailed derivation of the joint_P_cl_m10 is given in appendix <ref>. 
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[(^-1_1)^2] is 

    Tr[(^-1_1)^2] = 1/4 π∑^l_max_l = l_min [ (2 + 
    D_l-1/D_l + D_l/D_l-1)l + (2 + D_l+1/D_l + 
    D_l/D_l+1)(l + 1) ],

where D_l = C_l + N_l and (l_min, l_max) is the range of multipoles that are dipole 
modulated. Inspection of joint_P_cl_m10 reveals that m_10 is Gaussian distributed with 
mean μ and variance σ^2.  It is informative to have a qualitative look at the expression for σ^2. As C_l ≈ C_l+1 and CMB measurement are noise sub-dominant at large angular scales (N_l << C_l), the variance σ^2 is approximately equal to 2 π/(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 σ^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 

    w_z = m_10,  w_x = -√(2)m^r_11, and w_y = √(2) m^i_11.

Then, w_x, w_y, w_z are distributed with same variance σ^2. The norm of the signal r can be defined as 

    r ≡√(w^2_z + w^2_x + w^2_y ) = √(m_10^2 + 2m^r_11^2 + 2m^i_11^2).

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

    (r) = √(1/(2 πσ^2)^3) 4 π r^2 
     exp[-r^2/2 σ^2 ],

where σ^2 is given by var_m10. The relation r = √(4 π/3) A leads us to 
the probability distribution of the dipole amplitude A in case of SI map as 

    (A) = 8π/3√(2/3)A^2/σ^3exp[-2 π A^2/3 σ^2 ].

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 def_r.
Then the distribution of r marginalized over (θ, ϕ) is

    (r) = r/σ^2√(r/r_*)exp[ -r^2 + r^2_*/2 σ^2] I_1/2(r r_*/σ^2),

where I_1/2(x) is a modified Bessel function of first kind of order 1/2,

    I_1/2(x)=√(x/2π)∫^π_0 dθsin(θ) exp[ x cos(θ) ].

The probability distribution for A is

    (A) = 4 π/3A/σ^2√(A/A_* )exp[ -2π(A^2 + A^2_* )/3 σ^2] I_1/2(4 π A A_* /3σ^2),

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 (A) attains its maximum.



§ METHOD: HAMILTONIAN MONTE CARLO SAMPLING FOR NON-SI CMB SKY IN PRESENCE OF OBSERVATIONAL NON-IDEALITIES

For the given problem, we deal with ≈ l^2_max + 3 l_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 <cit.> as well as component separation <cit.>. 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 <cit.> as described in Appendix <ref>.



 §.§ 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 <cit.>, for CMB power spectrum estimation in 
<cit.>, for estimation of the large scale structure power spectrum in <cit.>, the inference of non-linear dynamics of large scale structure <cit.>, and lensing potential reconstruction from lensed CMB in <cit.>. HMC makes use of Hamiltonian Dynamics to sample 
a given probability distribution (q) of random variable q using a Hamiltonian 
() defined in the following way 

    = p^2/2μ - ln ((q)),

where q represents the position of this particle and p is the conjugate momentum corresponding 
to q. The quantity - ln ((q)) acts as the potential energy and μ is “the mass” 
assigned to “the particle”. In HMC, the choice of the mass parameter, μ 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

    q̇≡dq/dt = / p  and  ṗ≡dp/dt = -/ q.

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

    exp(-) = exp[ -p^2/2μ] (q).

In the end, we marginalize over p and obtain a fair sample of q drawn from (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, α) and corresponding momentum (p^r_lm, p^i_lm, 
p_C_l, p̅_10, p̅^̅r̅_11, p̅^̅i̅_11, p_α).
For the probability distribution defined in Prob_a_S


    (_lm, _C_l, _1N, p_α, ^r, ^i, {C_l}, m, α) 
       =   ∑^l_max_l = 2{p^2_C_l/2μ_C_l + ∑^l_m = 0[p^r^2_lm/2μ^r_a_lm + p^i^2_lm/2μ^i_a_lm] } + 
    p^2_α/2μ_α
       +   p̅^2_10/2μ̅_10 + p̅^̅r̅^2_11/2μ̅^r_11 + 
    p̅^̅i̅_11^2 /2 μ̅^i_11 
    -ln [ ( S_lml'm', a_lm | d_lm) ].

In above equation, p_lm and _C_l stand for the set of conjugate momenta for 
a_lm and C_l respectively. p̅_1N stands for conjugate momenta 
{p̅_10, p̅^r_11, p̅^i_11} of (m_10, m^r_11, m^i_11) and 
p_α is conjugate momentum of α. Further, μ^r_a_lm, μ^i_a_lm, μ_C_l are the 
masses assigned to the real part of a_lm, the imaginary part of a_lm and C_l, respectively. Whereas, μ̅_10, μ̅^r_11, μ̅^i_11, and μ_α are the masses assigned to m_10, m^r_11, m^i_11, and α, respectively.
We follow arguments given in <cit.> and <cit.> 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,

    μ^r_a_lm = μ^i_a_lm = (2/C_l + 2/N_l)   for   m ≠ 0 
       and  μ_a_l0 = (1/C_l + 1/N_l).

The mass for C_l is chosen as inverse of the cosmic-variance of C_l

    μ_C_l = 2l + 1/2C^2_l.

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 α.

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 Prob_a_S, 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 <ref>. 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 def_H 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, α) which are drawn from the multidimensional probability distribution 
(S_lml'm', a_lm | d_lm), where S_lml'm' depends on (C_l, 
m_10, m^r_11, m^i_11, α).



 §.§ 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  represent the vector of CMB temperature 
signal in the pixel space and  the data vector in the same space. We assume that noise is 
Gaussian distributed with variance σ^2_i for the i^th pixel and the noise between two 
pixels are not correlated. Then the probability distribution of  given  is

    ( | ) = 1/(2 π)^N_pix/2 ||^1/2exp[ -1/2∑^N_pix_i = 1(d_i - T_i)^2/σ^2_i],

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

    ṗ_lm = -1/2∑_l_1 m_1 S^-1_l_1m_1lm a^*_l_1m_1 
    + 1/2∑^N_pix_i = 11/σ^2_i[ ∑_l'm' 
    (d_l'm' - b_l p_l a_l'm') Y_l'm'(n̂_̂î) ]Y_lm(n̂_̂î),

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 σ^2_i on pixels.
We used  and  routines of [<http://healpix.sourceforge.net>]  <cit.>
to compute the second term in the above equation. Given the difference map in harmonic space, (d_l'm' - b_l p_l 
a_l'm'), the  subroutine executes the following computation

    M_D(n̂_i) = ∑_l'm'[ (d_l'm' - b_l p_l a_l'm') Y_l'm'(n̂_̂î) ].

With the map M_D(n̂_i) as input to the subroutine , we get following quantity as the output

    ∑^N_pix_i = 1M_D(n̂_i)/σ^2_i Y_lm(n̂_̂î) 4 π/N_pix.

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.



§ 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 collaboration[Planck-2015 temperature data are substantially similar to the Planck-2018 release that focused on large angle polarization anisotropies]. We use the SMICA CMB temperature map[Name of the 
file: 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 <ref>.
Since we are only interested in the large angular scale anomalies, we choose to carry out our analysis at the resolution of  = 
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 <cit.> available at NERSC[
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 <ref>.



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 <cit.>. Detailed results of the analysis with step model are given in section <ref>. Analysis for the power law model is given in section <ref>. A summary of the results for both these cases is provided in Table <ref> and <ref>.
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.







 §.§ Step-model of scale dependence

To compare the results of our method with the published literature <cit.>,
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 
<cit.>. Figure <ref> 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 × 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 <ref>. 
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 
<cit.>. We note that the current estimate of the parameter m^r_11 deviates from the BipoSH-MVE estimate by around 
1 σ. Other two parameters agree with their corresponding value from BipoSH-MVE.
It is evident from the mean (μ) and standard deviation (σ) of these parameters that each 
parameter is away from zero more than one σ. 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, σ = 0.039. The standard deviation for m^r_11 and m^i_11 is 
σ/√(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 <ref>. The posterior distributions of A, in 
figure <ref> and the posterior distribution of θ_p, ϕ_p in figure 
<ref> are obtained by transforming the Monte Carlo chains of 
{m_10, m^r_11, m^i_11} to {A, θ_p, ϕ_p} using coordinate transformation equations 
given in m1M_to_Athephi. We also quote the estimates of {A, θ_p, ϕ_p} from the previous analysis <cit.> in figure <ref> and <ref>. This analysis shows good agreement with all the previous analysis, with the maximum difference (about  one σ) with the BipoSH minimum variance estimator.





 §.§ Power-law model of scale dependence

The power-law model of dipole modulation is given by eq_power_law_Al. 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 <ref>, 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 × 10^4 samples. Figure <ref> also shows the posterior of the power law index α 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 <ref>. However, the choice of pivot multipole l_p can lead to  correlation between power law index α 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 <ref> 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 <ref>. 
The posterior distributions of {A(l_p = 16) , θ_p, ϕ_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) , θ_p, ϕ_p} using coordinate transformation formulas given in m1M_to_Athephi. Significant nonzero value of amplitude A(l_p = 16) along with the nonzero value of α points to the presence of the dipole modulation with scale dependence.  The recovered power-law profile of the CHA signal is depicted in figure <ref> along with the profile for the step-model and the Doppler Boost signal which is of the order 10^-3 <cit.>. 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 <cit.>. The recovered profile can also lead to direction dependence in the cosmological parameters <cit.>. 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 <ref>. 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 <cit.>. 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 <ref>) 
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 <ref>).
A more quantitative comparison between the direction recovered in our analysis and the direction reported in the literature is given in figure <ref>. 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 <ref> 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 <cit.> to the histogram of C_l samples. We use  
routine  to implement the fitting <cit.>. 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 <ref> 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 ΛCDM angular power spectrum provided by Planck <cit.>. 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∈{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 ΛCDM values. The corresponding results are shown in Appendix <ref>.  

The recovered power-law of the dipole modulation signal depicted in figure <ref> 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 cov_mat). 
We estimate the correction to the angular power spectrum due to second order term and depict the relative correction to the angular power spectrum Δ_l in figure <ref>. For comparison, we also show the 1σ 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.
















§ 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 ((M_0|d)) and the probability of the model M_1 
given data ((M_1|d)) using the following ratio called Bayes factor 

    (M_0|d)/(M_1|d) = ℰ(d|M_0) (M_0)/ℰ(d|M_1) (M_1).

In the above equation, ℰ(d|M) stands for the probability of the data given model M and ℰ(d|M) is 
called Bayesian Evidence for M. (M) stands for the prior probability of the model M.
In the situation where two models are considered a priori equally probable ((M_0) = (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 

    B_SI-DM≡ℰ(d| M_SI) /ℰ(d | M_DM) ,

where we have assumed that prior probability (M_SI) = (M_DM).
In above equation, ℰ(d| M_SI) and ℰ(d| M_DM) are evidence under  M_SI and M_DM, respectively, which is defined as

    ℰ(d|M) = ∫(w_x, w_y, w_z, α ) Π(w_x, w_y, w_z, α | M) dw_xdw_y dw_z dα,

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 def_wxyz. (w_x, w_y, w_z, α ) denotes the likelihood and Π(w_x, w_y, w_z, α | 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 α. 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 α. The Bayes factor for nested models can be obtained by using Savage-Dickey density ratio (SDDR) <cit.>. For SI model nested in the power law dipole modulation model, the SDDR is

    B_SI-DM = _C_l(w_x = 0, w_y = 0, w_z = 0,α |d,DM)/Π(w_x = 0, w_y = 0, w_z = 0, α | DM),

where _C_l(w_x = 0, w_y = 0, w_z = 0,α |d,DM) is posterior marginalized over C_l. α dependence of SDDR in SD_ratio_1 is weak because the likelihood conditioned at (w_x, w_y, w_z) = (0,0,0) is independent of α. The prior mentioned in (<ref>) can be written as 
    Π(w_x, w_y, w_z, α | M_DM) = Π(w_x, w_y, w_z |α; M_DM) Π(α | M_DM).

We choose uniform prior on α  over the range α = -2  to  0. Hence, the normalized prior density of α is

    Π(α | M_DM) = 1/2  for -2 ≤α≤ 0.

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

    A(l) < 1/2    ∀   l.

In spherical polar coordinates, with r as defined in def_r, the above argument gives an upper bound: 𝐫 ≤√((π/3))≈ 1, which is denoted by R. To obey prior_Al, we need to make the upper bound on r dependent on α. 
Otherwise, for low values of α close to -2, A(l) can become greater than 1/2 at low multipoles for sufficiently high A(l_p). Hence, α dependent upper bound R(α) is

    R(α) = 1/2( l/l_p)^-α√(4 π/3) = 
     1.02 ( l/l_p)^-α.

At α = -2, with l_p = 16 and l = 2, R(α) = 0.016.   

We choose a prior on (w_x, w_y, w_z) which has uniform density within the sphere of radius R(α), centred at (w_x, w_y, w_z) = (0,0,0) and zero outside. Normalized form of such prior is

    Π_1(w_x, w_y, w_z|α, DM) = 
    3/4 π R^3(α)   if √(w_x^2 + w_y^2 + w_z^2)≤ R(α) 
    
        0                 otherwise.

Note that the prior density is independent of the direction. The form of prior in Prior_xyz ensures positive definiteness of the covariance matrix at all scales for the range of α considered in the analysis.
Choosing an uniform density prior as given in Prior_xyz implies that the angle marginalized prior density of the amplitude is

    Π_1(r | α, DM) = 3r^2/R^3(α)  for   r ≤ R(α).

As mentioned in section <ref>, 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 Prior_xyz. With this procedure, the SDDR in favour of the SI model is found to be

    B_SI-DM≈ 0.4.

In Appendix <ref>, 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), α} 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 <ref>. 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.



§ 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 <cit.>.
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 (α). The parameters θ_p and ϕ_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 <cit.>. 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 ± 0.022 and α = -0.92 ± 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 σ from the previously reported value by Aiola et al. <cit.>, 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 <cit.> due to our local motion with a velocity β≡ v/c= 1.23 × 10^-3 <cit.>, which is also measured by Planck using the high resolution temperature data <cit.>. 
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 <cit.>. 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.

20pt
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 [<http://hpc.iucaa.in>]. Some of the results in this paper have been derived using the HEALPix package 
<cit.>. We also acknowledge the use of  <cit.>, 
 <cit.>, Python packages  <cit.>, 
 <cit.>, <cit.> and  <cit.>.



§ DERIVATION OF LG


This appendix provides a detailed derivation of P_of_S_given_d and joint_P_cl_m10. To derive 
P_of_S_given_d, we marginalize over  from Prob_a_S. We simplify the following expression in Prob_a_S

    E = (d - a)^†N^-1  (d -a) 
    + a^†S^-1a,

by completing the square <cit.>.

    E = ^†( + )^-1 + 
    ( - )^†(^-1 + ^-1) ( - ),

where = ( + )^-1 is Wiener filter applied to the data.
Prob_a_S now becomes

    (,  | ) = 1/√(|| || ) (2 π)^nexp[-1/2( - )^†(^-1 + ^-1)( - ) ] exp[ -1/2^†( + )^-1 ].

Using the results

    ∫exp[ -1/2 ( - )^† (^-1 + ^-1) 
     ( - ) ] d = (2 π)^n/2√(|(^-1 + ^-1)^-1|),

and 

    |(^-1 + ^-1)| = |( + )|/||||,

in Prob_a_S to get P(S| d) as given in P_of_S_given_d. Now we derive 
the joint probability distribution C_l and m_10 given in joint_P_cl_m10 from P_of_S_given_d. We write the sum of signal covariance matrix  and noise covariance matrix  in the following manner

    +  =  + m_10_1.

Since we assume  to be diagonal, it contributes only to . _1 is the off-diagonal 
part of  without m_10. Taylor series expansion of inverse of the sum +, upto 
second order in m_10 gives

    ( + )^-1 = ( + m_10_1)^-1 = 
     ^-1 - m_10^-1_1 ^-1 + m^2_10 (^-1_1)^2 ^-1.

Determinant | + | in P_of_S_given_d also has the m_10 dependence expressed as 

    | + | = || (| + m_10^-1_1|).

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

    ln(| + m_10^-1_1|) = Tr[ln( + m_10^-1_1)].

Taylor series expansion of logarithm gives

    ln( + m_10^-1_1) = ∑^∞_k = 1(-1)^k+1/k (m_10^-1_1)^k.

Then, keeping the terms up to second order in m_10, we get 

    ln( + m_10^-1_1) = m_10^-1_1 -m^2_10/2 (^-1_1)^2.


    Tr[ln( + m_10^-1_1)]  = 
     m_10 Tr[^-1_1] - m^2_10/2 Tr[(^-1_1)^2].


    ln(| + m_10^-1_1|)  = 
    - m^2_10/2 Tr[(^-1_1)^2]         ∵ Tr[^-1_1] = 0.


    ln(| + |) = ln(||) - m^2_10/2 Tr[(^-1_1)^2].

We substitute matrix_sum_inverse and matrix_log_det in the following form of 
P_of_S_given_d

    2 ln((|)) = 2 ln((C_l, m_10 |)) = 
     - n ln(2π) - ln| + | - ^†( + )^-1,

to obtain

    -    2 ln((C_l, m_10 |)) = 
     n ln(2π) + ln|| + ^†^-1
       -    m_10^†^-1_1 ^-1 - 
     m^2_10/2 Tr[(^-1_1)^2] 
     + m^2_10^† (^-1_1)^2 ^-1.

P_Cl_m10_semifinal_eq can be further rearranged in following manner

    - ln((C_l, m_10 |)) 
        =   n/2ln(2π) + 1/2ln|| + 1/2^†^-1
       +   Tr[(^-1_1)^2]/4 M { m^2_10 - 2 m_10^†^-1_1 ^-1/Tr[(^-1_1)^2] M^-1},

where M represents following quantity

    M = [ 2 ^† (^-1_1)^2 ^-1/Tr[(^-1_1)^2] - 1 ].

 joint_P_cl_m10 follows from P_Cl_m10_final_eq given above.



§ 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)

    ȧ^̇ṙ_lm = p^r_lm/μ^r_lm,   ṗ^̇ṙ_lm = -∂/∂ 
    a^r_lm,  ȧ^̇i̇_lm = p^i_lm/μ^i_lm,   and  ṗ^̇i̇_lm = -∂/∂ 
    a^i_lm.

This involves time derivative of conjugate momentum corresponding to a_lm

    ṗ_lm = -1/2^†^-1/ a_lm - 
     1/2 [( - )^†^-1 ( - )]/ a_lm.

After taking derivative with respect to a_lm, above equation reduces to

    ṗ_lm = -1/2∑_l_1 m_1 S^-1_l_1m_1lm a^*_l_1m_1 + 
     1/2∑_l_1 m_1N^-1_l_1m_1lm (d^*_l_1m_1 - a^*_l_1m_1).

To compute ṗ^̇ṙ_lm and ṗ^̇i̇_lm, the momentum derivatives corresponding to the real and imaginary part of a_lm respectively, we make use of the following relations

    / a^r_lm = 2 [ / a_lm] 
       and  / a^i_lm = -2 [ / a_lm],

where  and  denote the real and imaginary part of the complex number, respectively.
Hamilton's equations for (m_10, m^r_11, m^i_11) are

    ṁ_10 = p̅_10/μ_10,  ̅̇p_10 = 
    		-∂/∂ m_10,
      ṁ^̇ṙ_11 = p̅^̅r̅_11/μ^r_11,  ̅̇p^r_11 = 
                    -∂/∂ m^r_11,
      ṁ^̇i̇_11 = p̅^̅i̅_11/μ^i_11,  ̅̇p^i_11 = 
                    -∂/∂ m^i_11.


    / m^r_11 = 2[ / m_11]
       and  / m^i_11 = -2[ / m_11],

where / m_1N are given in terms of derivative of  with respect to 
A^1N_ll+1 as

    / m_1N = 
    2 ∑_l A^1N_l l+1/ m_1N/ A^1N_ll+1.

The factor of 2 in above expression appears because A^1N_l l+1 = A^1N_l+1 l for the dipole modulation model.
The time derivative of the momentum corresponding to α is obtained as

    ṗ_α = /α = ∑_l  [ / 
    A^10_ll+1(l_p, α)  A^10_ll+1(l_p, α)/α + 
    / A^11_ll+1(l_p, α)  A^11_ll+1(l_p, α)/α ].

From Pm1M_dot and Palpha_dot, we note that to compute derivatives of momentum 
corresponding to m_10, m^r_11, m^i_11 and α, we need to know / A^1M_ll+1. 

The time derivative of conjugate momentum corresponding to the BipoSH coefficient A^LN_ll' is

    ṗ^LN_ll' = - / A^LN_ll' = 
     - 1/2 ln ||/ A^LN_ll' 
     - 1/2^†^-1/ A^LN_ll'.


Following is a detailed derivation of the expression for ṗ^LN_ll' starting from PALMll_dot.
To explicitly evaluate palm_dot_2 and to simplify PALMll_dot we use a Taylor series expansion of ^-1 up to first order in  as 
 
    ^-1 =  - ^-1^-1.

We make use of two standard results from matrix calculus

    ln ||/ A^LN_ll' = Tr[ ^-1/ A^LN_ll'],

and 

    ^†^-1/ A^LN_ll' = 
     -Tr[ (^-1) (^†^-1) / A^LN_ll'].

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

    S_l_1 m_1 l_2 m_2/ A^LN_ll' = (-1)^m_2 C^LN_l_1 m_1 l_2 -m_2δ_l_1lδ_l_2l'.

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

    S^-1_l' m' l m = D_l' m' l m - D^-1_l' m' l' m'O_l' m' l m D^-1_l m l m.

Using the aboe Taylor series expression and the derivative in der_Sl1m1l2m2,
der_ln_matrix can be expressed as

    ln ||/ A^LN_ll' = ∑_m'(-1)^m'/D_l'm'l'm' C^LN_lm'l'-m'δ_ll' - ∑_m', m(-1)^m' O_l'm'lm/D_l'm'l'm' D_lmlm C^LN_lml'-m'.

Using the derivative in der_Sl1m1l2m2, der_inv_matrix becomes

    ^†^-1/ A^LN_ll' = 
     - ∑_m', m (^-1)_l'm' (^†^-1)_lm (-1)^m_1 C^LN_lml'-m'.

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

    (^-1)_l_1m_1 (^†^-1)_l_2m_2   =   a_l_1 m_1 a^*_l_2 m_2/D_ l_1 m_1 l_1 m_1 D_l_2 m_2 l_2 m_2 
    - a_l_1 m_1/D_ l_1 m_1 l_1 m_1 D_l_2 m_2 l_2 m_2∑_l'_2 m'_2O_l'_2 m'_2 l_2 m_2 a^*_l'_2 m'_2/D_l'_2 m'_2 l'_2 m'_2
       -   a_l_2 m_2/D_ l_1 m_1 l_1 m_1 D_l_2 m_2 l_2 m_2∑_l'_1 m'_1O_l'_1 m'_1 l_1 m_1 a^*_l'_1 m'_1/D_l'_1 m'_1 l'_1 m'_1
       +   ∑_l'_1 m'_1∑_ l'_2 m'_2O_l_1 m_1 l'_1 m'_1 O_l'_2 m'_2 l_2 m_2/D_l'_1 m'_1 l'_1 m'_1 D_l'_2 m'_2 l'_2 m'_2a_l'_1 m'_1 a^*_l'_2 m'_2/D_ l_1 m_1 l_1 m_1 D_l_2 m_2 l_2 m_2.

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

    (^-1)_l_1m_1 (^†^-1)_l_2m_2 
    = a_l_1 m_1 a^*_l_2 m_2/D_ l_1 m_1 l_1 m_1 D_l_2 m_2 l_2 m_2 
    - O_l_1 m_1 l_2 m_2/D_ l_1 m_1 l_1 m_1 D_l_2 m_2 l_2 m_2
    - O_l_2 m_2 l_1 m_1/D_ l_1 m_1 l_1 m_1 D_l_2 m_2 l_2 m_2.

Eq. (<ref>) along with above simplification and der_ln_matrix_2 leads to following 
expression 
for ṗ^LN_ll'

    ṗ^LN_ll' =    -   1/2∑_m'(-1)^m'/D_l'm'l'm' C^LN_lm'l'-m'δ_ll'
       -   1/2∑_m', m(-1)^m' O_l'm'lm/D_l'm'l'm' D_lmlm C^LN_lml'-m'
     + 1/2∑_m', m(-1)^m' a_l'm' a^*_lm/D_l'm'l'm' D_lmlm 
    C^LN_lml'-m'.

The time derivative of the conjugate momentum of C_l is obtained from PALMll_dot_2 by 
substituting L = 0 and M = 0 in the above equation (as A^00_ll are related to the C_l)

    ṗ^00_ll' = 
     2l + 1/2A^00_ll (Â^00_ll'/A^00_ll - 1),
       where  Â^00_ll'  = ∑_m m' a_lm a^*_l'm' C^00_lml'm'.

By using relation between A^00_ll and C_l, we get

    ṗ_C_l = 2l + 1/2 C_l(Ĉ_l/C_l - 1),   where  Ĉ_l = 1/2l + 1∑_m |a_lm|^2.

Equation <ref> is same as equation (25) in <cit.>. Note that, in 
eq_PA00ll_dot and eq_PCl_dot, A^00_ll' and C_l are the variables of dynamics, 
whereas Â^00_ll' and Ĉ_l are quantities obtained using other variables of 
dynamics, namely {a_lm}.





§ 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  algorithm <cit.>. We add the anisotropic noise <cit.> 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 α = -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 <cit.>. 
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 <ref>, figure <ref>, figure <ref> and figure <ref>.

Histograms in figure <ref> 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 α. These distributions are obtained using Monte Carlo chain of 10^5 samples after discarding 5 × 10^4 samples as Burn-In. The set of 1D histograms, marginalized 
distributions of the parameters, in figure <ref> 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 var_m10 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, α = -0.7 and l = 2 to l = 256. This leads to the following: σ_m_10 = 0.023 and σ_m^r_11 = σ_m^i_11 = 0.016. The corresponding figures given in figure 
<ref> 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), θ_p, ϕ_p using m1M_to_Athephi. Figure <ref> shows the distribution of A(l_p = 32) obtained from the harmonic space variables shown in figure <ref>. Figure <ref> indicates that our method does recover the input value present in the simulated map. Figure <ref> shows the distribution of θ_p, ϕ_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 <ref>. 
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 <cit.> to the histogram of C_l samples. We use  
routine  to do the fitting <cit.>. In figure 
<ref> 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. 









§ 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 α are given in the figure <ref>. Taking α = -0.92 and extrapolating from the results given in figure <ref> 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 <ref>.  




§ HMC ANALYSIS OF CHA WITH FIXED FIDUCIAL LGCDM LG

From DM_biposh and shape_factor 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 <cit.>. 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), α where the angular power spectrum C_l is held fixed at the best-fit ΛCDM power spectrum provided by Planck <cit.>. The distribution of A(l_p) and that of (θ, ϕ) given in figure <ref> and figure <ref> respectively are similar to those given in figure <ref> and figure <ref> with C_l as variable.





§ ESTABLISHING ROBUSTNESS OF THE EVIDENCE RATIO AGAINST PRIOR CHOICE

In section <ref>, we discussed a prior, say, Prior-1, which is spherically symmetric and has a constant density within the sphere of radius R(α). 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(α) and zero afterward. The normalized form of such a prior is

    Π_2(w_x, w_y, w_z|α, DM) = 
     1/2 π R^2(α) √(w^2_x + w^2_y + w^2_z)   if √(w^2_x + w^2_y + w^2_z)≤ 
    R(α) 
    
        0                 otherwise.

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

    Π_2(r|α, DM)  = 2 r/R(α)^2  for  r ≤ R(α).

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(α). Prior-3 is again zero outside the sphere of radius R(α). The normalized form of such a prior is

    Π_3(w_x, w_y, w_z|α, DM) = 
     1/4 π R(α) (w_x^2 + w_y^2 + w_z^2)   if √(w_x^2 + w_y^2 + w_z^2)≤ 
    R(α) 
    
        0                 otherwise.

The prior density Π_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

    Π_3(r|α, DM) = 1/R(α)  for  r ≤ R(α).


The prior densities Π_2 and Π_3 diverge as (w_x,w_y,w_z) → (0,0,0). Hence we can not use the formula for the SDDR given in SD_ratio_1 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) = (ϵ, ϵ, ϵ) close to (x,y,z) = (0,0,0), where ϵ 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=ϵ. Modifying SD_ratio_1, the Bayes factor in favour of the SI model is

    B_SI-DM = 
    _C_l(w_x = ϵ, w_y = ϵ, w_z = ϵ,α |d,DM)/Π(w_x = ϵ, w_y = ϵ, w_z = ϵ, α | DM).

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 _C_l(w_x = ϵ, w_y = ϵ, w_z = ϵ,α |d,DM), which appears at the numerator of SD_ratio_2, using the Monte-Carlo samples and need a finite bin-width (around 3 × 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 ϵ as long as (w_x, w_y, w_z) = (0,0,0) and (w_x, w_y, w_z) = (ϵ, ϵ, ϵ) are within the same bin. 
We obtain the ratio of posterior density and prior density in the grid containing the point (x, y, z) = (ϵ, ϵ, ϵ) with ϵ = 10^-3. SDDR thus obtained for Prior-2 and Prior-3 is given in Table  <ref>. 

From SD_ratio_2, SDDR depends on α. However, at the nested point, which reduces the power law dipole modulation model to SI model, the likelihood does not depend on α and the SDDR also remains roughly constant over the range of α. In Table  <ref> we quote the value of SDDR averaged over α.











JHEP

