\documentclass[10pt,letterpaper]{article}

%\usepackage[utf8]{inputenc}
%\usepackage{times}

\usepackage{amssymb,amsmath}
\usepackage{feynmp}
\usepackage{xfrac}
\usepackage{mathrsfs}
\usepackage{bm}
%\usepackage{color}
\usepackage[usenames, dvipsnames]{color}
\usepackage{listings}
\usepackage{natbib}
\usepackage{graphicx}
\usepackage{float}
\graphicspath{{./Pics/}}
\usepackage{cancel}
\usepackage{amsthm}
\usepackage{cancel}
\usepackage{hyperref}
\usepackage{titlesec}
\usepackage{fancyvrb}
\usepackage{dsfont}
\usepackage{scrextend}
%\usepackage{subfigure}
\usepackage{subcaption}
\usepackage{enumitem}
\usepackage{framed}
\usepackage{wrapfig}
\usepackage{multicol}
\usepackage{empheq}
\usepackage{booktabs}
\usepackage[listings]{tcolorbox}
\tcbuselibrary{theorems,breakable}
\titlelabel{\thetitle.\quad}
\parskip 1ex 
\setlength\parindent{0pt}
\usepackage[margin=0.8in]{geometry}
\usepackage{lastpage}

\def \ul{\underline}
\def \vp{\varphi}
\def \ds{\displaystyle}
\def \pa{\partial}
\def \tcr{\textcolor{red}}
\def \tcb{\textcolor{blue}}
\def \tcg{\textcolor{green}}

\def \dangsym{{\fontencoding{U}\fontfamily{futs}\selectfont\char 66\relax}}



\newcommand{\dbar}{d\hspace*{-0.08em}\bar{}\hspace*{0.1em}}

\allowdisplaybreaks


\pagenumbering{arabic}
\usepackage{fancyhdr}
\setlength{\headheight}{15.2pt}
\pagestyle{fancy}
\rhead{\thepage}
\cfoot{}
\rhead[<even output>]{Page\ \thepage\ of\ \pageref{LastPage}}
\chead[<even output>]{CMB}
\lhead[<even output>]{\textit{Sayan}}

\newcommand{\com}[1]{{\bf {\color{Turquoise} Comment:  #1}}}

\usepackage[auth-lg]{authblk}
\author[1]{Sayan Mandal}
\affil[1]{\small McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213}
\begin{document}
\title{CMB Temperature Spectrum in Bianchi-I Spacetime}
\date{\today\vspace{-4ex}}
\maketitle
\noindent\makebox[\linewidth]{\rule{\textwidth}{1pt}}

\newenvironment{colored}{\par\color{blue}}{\par}
\newenvironment{colored2}{\par\color{magenta}}{\par}
\maketitle


We follow the outline of \textbf{arXiv:1102.5235} and \textbf{arXiv:1007.4539}.\\

\begin{tcolorbox}[colback=cyan!10!white,colframe=black!25!white,coltitle=brown!25!black,breakable,title=\textbf{Mathematical notes}]
To be clear with conventions, let us recall the definitions of the special functions used here. The \textbf{\textit{spherical harmonics}} are defined for $l,m\in\mathbb{Z}\cup\{0\}, |m|\leq l$ as,
\begin{equation}\label{eSphHarm}
Y_{lm}(\theta,\phi)=\sqrt{\frac{2l+1}{4\pi}\cdot\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta)\,e^{im\phi}
\end{equation}
where $P_l^m(x)$ are the \textbf{\textit{associated Legendre polynomials}} defined in terms of the ordinary Legendre polynomials $P_l(x)$,
\begin{equation}\label{eAssocLegen}
P_l^m(x)=(-1)^m(1-x^2)^{\sfrac{m}{2}}\frac{d^m}{dx^m}P_l(x)
\end{equation}
These polynomials are normalized as,
\begin{equation}\label{eAssocLegenNormal}
\int_{-1}^{+1} dx\,P_l^m(x)\,P_{l'}^m(x)=\frac{2(l+m)!}{(2l+1)(l-m)!}\delta_{ll'}
\end{equation}
and the spherical harmonics themselves are normalized as,
\begin{equation}\label{eSphHarmNormal}
\int d\Omega\,Y_{lm}(\theta,\phi)\,Y^*_{l'm'}(\theta,\phi)=\delta_{ll'}\delta_{mm'}
\end{equation}
An useful recurrence relation for a fixed $m$ is,
\begin{equation}\label{eAssocLegenFixedMRecur}
(l-m+1)P_{l+1}^m(x)=(2l+1)xP_l^m(x)-(l+m)P_{l-1}^m(x)
\end{equation}
We also note that an exponential can be decomposed in terms of spherical harmonics and the modified Bessel function,
\begin{equation}\label{ePlaneWaveExpn}
e^{i\mathbf{k}\cdot\mathbf{x}}=4\pi\sum_{L=0}^\infty\sum_{m=-L}^L i^L j_L(kx)Y^*_{LM}(\hat{\mathbf{k}})Y_{LM}(\hat{\mathbf{x}})
\end{equation}
\end{tcolorbox}

We choose, without loss of generalization, that the Bianchi-I spacetime is axisymmetric about the $z$-axis, and the metric can be written as,
\begin{equation}\label{eMetric}
ds^2=-dt^2+a^2(t)\left[e^{-2\sigma}\left(dx^2+dy^2\right)+e^{4\sigma}dz^2\right]
\end{equation}
We denote the scale factor in the $z$ direction as $a_\parallel$ and that in the $x$ and $y$ directions as $a_\perp$,
\begin{equation}\label{eScaleFac}
a_\parallel=ae^{2\sigma},\quad\quad\quad a_\perp=ae^{-\sigma}
\end{equation}
The `average' scale factor is defined as,
\begin{equation}\label{eScaleFacAvg}
a_\mathrm{avg}\equiv \left(a_\parallel a^2_\perp\right)^{\sfrac{1}{3}}=a
\end{equation}
We work in the regime where $\sigma$ is small, i.e., $\sigma\ll 1$, and the difference of the scale factors is also small,
\begin{equation}\label{eScaleFacDiff}
\delta\equiv\frac{a_\parallel-a_\perp}{a}=e^{2\sigma}-e^{-\sigma}= 3\left(\sigma+\frac{\sigma^2}{2}\right)+\mathcal{O}(\sigma^3)
\end{equation}
The four momentum of the photon is written as,
\begin{equation}\label{eFourMom}
p=T_0\left(\alpha,\frac{n_1}{a_\perp},\frac{n_2}{a_\perp},\frac{n_3}{a_\parallel}\right)
\end{equation}
Here, $\hat{\mathbf{n}}=(n_1,n_2,n_3)$ is a unit vector in the direction of propagation of the photon. To determine $\alpha$, we use the condition that $p_\mu p^\mu=0$,
\begin{equation}\label{eAlpha}
\alpha=\frac{\sqrt{e^{2\sigma}\sin^2\theta+e^{-4\sigma}\cos^2\theta}}{a}
\end{equation}
A comoving observer sees the temperature,
\begin{equation}\label{eTempGen}
T^\sigma(t,\theta)=\frac{T_0\sqrt{e^{2\sigma}\sin^2\theta+e^{-4\sigma}\cos^2\theta}}{a}
\end{equation}
where the subscript $\sigma$ denotes the contribution due to the anisotropic spacetime. The direction-averaged temperature is calculated as,
\begin{equation}\label{eTempGenAvg}
\bar{T}^\sigma(t)=\frac{1}{4\pi}\int d\Omega\,T(t,\theta)
\end{equation}
The temperature fluctuation are expanded as usual,
\begin{equation}\label{eTempFlucExp}
\Theta^{\sigma}(\theta,\phi)\equiv\left.\frac{\Delta T}{\bar{T}}\right|^\sigma (\theta,\phi)=\frac{T^\sigma(t,\theta)-\bar{T}^\sigma}{\bar{T}^\sigma}=\sum_{l=0}^\infty\sum_{m=-l}^l \bar{a}_{lm}Y_{lm}(\theta,\phi)
\end{equation}

Let us now consider the case of an isotropic background, i.e., instead of using eq. \eqref{eTempGenAvg} to calculate the direction-averaged temperature, we set it to $T_0/a$,
\begin{equation}\label{eTempIsoAvg}
\bar{T}^\sigma(t)=\frac{T_0}{a}
\end{equation}
In this case we can write,
\begin{equation}\label{eTempFlucIso}
\Theta^{\sigma}(\theta,\phi)=\sqrt{(1+2\sigma)\sin^2\theta+(1-4\sigma)\cos^2\theta}-1=\sqrt{1+2\sigma(1-3\cos^2\theta)}-1\simeq \sigma(1-3\cos^2\theta)
\end{equation}
where the final approximation is for small $\sigma$. We can rewrite this as,
\begin{equation}\label{eTempFlucIso2}
\tcboxmath{\Theta^{\sigma}(\theta,\phi)=-4\sigma\sqrt{\frac{\pi}{5}}Y_{20}(\theta,\phi)+\mathcal{O}(\sigma^2)}
\end{equation}
We can extract the amplitude in a similar way as in eq. \eqref{eTempO1FlucAmp}, giving us back the same expressions,
\begin{equation}\label{eTempFlucIsoAmp}
\bar{a}_{20}=-4\sigma\sqrt{\frac{\pi}{5}} \quad\quad\quad \Rightarrow \quad\quad\quad \bar{C}_2=\frac{16\pi\sigma^2}{25}
\end{equation}


\begin{tcolorbox}[colback=cyan!10!white,colframe=black!25!white,coltitle=brown!25!black,breakable,title=$\sigma$ \textbf{dependent background}]
The use of eq. \eqref{eTempIsoAvg} is justified for small $\sigma$, and to show this, let us consider the case when the background is anisotropic, i.e., the direction-averaged temperature is dependent on $\sigma$. Expanding the exponentials in eq. \eqref{eTempGen} to linear order in $\sigma$, we get,
\begin{equation}\label{eTempO1}
T^\sigma(t,\theta)=\frac{T_0\sqrt{(1+2\sigma)\sin^2\theta+(1-4\sigma)\cos^2\theta}}{a}
\end{equation}
and the averaged temperature is,
\begin{equation}\label{eTempO1Avg}
\bar{T}^\sigma(t)=\frac{T_0}{24 a\sqrt{\sigma(1+2\sigma)}}\left[\sqrt{6}\pi(1+2\sigma)^{\sfrac{3}{2}}+12\sqrt{\sigma(1-2\sigma-8\sigma^2)}-2\sqrt{6}(1+2\sigma)^{\sfrac{3}{2}}\sin^{-1}\sqrt{\frac{1-4\sigma}{1+2\sigma}}\right]
\end{equation}
and it is easily checked that,
\begin{equation}\label{eTempO1AvgLim}
\lim_{\sigma\rightarrow 0}\bar{T}^\sigma(t)=\frac{T_0}{a}
\end{equation}
The temperature fluctuations as expanded like in eq. \eqref{eTempFlucExp}, where $\Delta T/\bar{T}|_\sigma$ is kept up to $\mathcal{O}(\sigma)$,
and the only non-vanishing amplitude for $l=2$ is calculated as,
\begin{equation}\label{eTempO1FlucAmp}
\bar{a}_{20}=\frac{1}{4\pi}\int d\Omega\,\left.\frac{\Delta T}{\bar{T}}\right|^\sigma\,Y^*_{20}(\theta,\phi)=-4\sigma\sqrt{\frac{\pi}{5}}
\end{equation}
giving the quadrupole contribution,
\begin{equation}\label{eTempO1Quad}
\bar{C}_2=\frac{1}{5}\sum_{m=-2}^{2}|\bar{a}_{2m}|^2=\frac{16\pi\sigma^2}{25}
\end{equation}
These agree with eq. \eqref{eTempFlucIsoAmp}. \\~\\

If we keep terms up to $\mathcal{O}(\sigma^2)$ inside the square root in eq. \eqref{eTempGen}, we get,
\begin{equation}\label{eTempO2}
T^\sigma(t,\theta)=\frac{T_0\sqrt{(1+2\sigma+2\sigma^2)\sin^2\theta+(1-4\sigma+8\sigma^2)\cos^2\theta}}{a}
\end{equation}
If we also expand the ensuing expression $\Delta T/\bar{T}|^\sigma$ up to $\mathcal{O}(\sigma^2)$, we end up with the quadrupole contribution,
\begin{equation}\label{eTempO2Quad}
\bar{C}_2=\frac{16\pi(7-5\sigma)^2\sigma^2}{1225}
\end{equation}
which agrees with eq. \eqref{eTempO1Quad} for small $\sigma$.
\end{tcolorbox}



We will now work with the expressions in eq. \eqref{eTempFlucIso2}. The temperature fluctuation for the general case considering the Bianchi metric can be written as,
\begin{equation}\label{eTempFlucFull}
\tcboxmath{\Theta(\theta,\phi)=\left[1-4\sigma\sqrt{\frac{\pi}{5}}Y_{20}(\theta,\phi)\right]\Theta^{(\sigma=0)}(\theta,\phi)+\mathcal{O}(\sigma^2)}
\end{equation}
where $\Theta^{(\sigma=0)}(\theta,\phi)$ denotes the temperature anisotropies in an FLRW universe. These can be expanded as,
 \begin{equation}\label{eTempFlucFRWExp}
\Theta^{(\sigma=0)}(\theta,\phi)=\sum_{l=0}^\infty\sum_{m=-l}^l b_{lm}Y_{lm}(\theta,\phi), \quad\quad\quad \Theta(\theta,\phi)=\sum_{l=0}^\infty\sum_{m=-l}^l a_{lm}Y_{lm}(\theta,\phi)
 \end{equation}
 
 \begin{tcolorbox}[colback=cyan!10!white,colframe=black!25!white,coltitle=brown!25!black,breakable,title=\textbf{The amplitudes}]
 We can thus extract the amplitudes $a_{lm}$ in eq. \eqref{eTempFlucFRWExp} as,
 \begin{equation}\label{eTempFullAmpl}
 \begin{aligned}
 a_{lm}&=\sum_{l'=0}^\infty\sum_{m'=-l'}^{l'} \int_0^{2\pi}d\phi\int_0^\pi d\theta\,\sin\theta \left[1-4\sigma\sqrt{\frac{\pi}{5}}Y_{20}(\theta,\phi)\right] b_{l'm'}Y_{l'm'}(\theta,\phi) Y^*_{lm}(\theta,\phi)+\mathcal{O}(\sigma^2) \\
 &=\sum_{l'=0}^\infty\sum_{m'=-l'}^{l'} b_{l'm'}\left\{\int d\Omega\,Y_{l'm'}(\theta,\phi) Y^*_{lm}(\theta,\phi) -4\sigma\sqrt{\frac{\pi}{5}}\int d\Omega\,Y_{20}(\theta,\phi)Y_{l'm'}(\theta,\phi) Y^*_{lm}(\theta,\phi)\right\}+\mathcal{O}(\sigma^2) \\
 &=b_{lm}-\sigma\sum_{l'=0}^\infty b_{l'm}\left[\frac{(2l+1)(2l'+1)(l-m)!(l'-m)!}{(l+m)!(l'+m)!}\right]^{\sfrac{1}{2}}\int_{-1}^1 dx\,P_2^0(x)\,P_{l'}^m(x)\,P_l^m(x)+\mathcal{O}(\sigma^2)
 \end{aligned}
 \end{equation}
 Let us now tackle the integral above,
 \begin{equation}\label{eThreePlm}
 \begin{aligned}
 \mathscr{I}_{ll'}^m&\equiv\int_{-1}^1 dx\,P_2^0(x)\,P_{l'}^m(x)\,P_l^m(x) \\
 &=\frac{1}{2}\int_{-1}^1 dx\,(3x^2-1)\,P_{l'}^m(x)\,P_l^m(x) \\
 &=\frac{3}{2}\mathscr{A}_{ll'}^m-\frac{(l+m)!}{(2l+1)(l-m)!}\delta_{l'l}
 \end{aligned}
 \end{equation}
 where $\mathscr{A}_{ll'}^m=\ds{\int_{-1}^1 dx\,x^2\,P_{l'}^m(x)\,P_l^m(x)}$. To evaluate this, we make use of the relation in eq. \eqref{eAssocLegenFixedMRecur} to obtain,
 \begin{equation}\label{eAssocLegenFixedMRecur-2}
 xP_l^m(x)=\frac{(l-m+1)P_{l+1}^m(x)+(l+m)P_{l-1}^m(x)}{2l+1}
 \end{equation}
 We plug this into the definition of $\mathscr{A}_{ll'}^m$, and using the expression in eq. \eqref{eAssocLegenNormal}, we finally arrive at,
 \begin{equation}\label{eScri1}
 \begin{aligned}
 \mathscr{I}_{ll'}^m&=\frac{3(l+m)!}{(2l+1)(2l'+1)}\left[\frac{(l'+m)(l+m)}{(2l+3)(l-m)!}\delta_{l',l+2}+\frac{(l'-m+1)}{(2l-1)(l-m-1)!}\delta_{l',l-2}\right. \\
 &+\left.\left\{\frac{(l'-m+1)(l+m)}{(2l+3)(l-m)!}+\frac{(l'+m)}{(2l-1)(l-m-1)!}-\frac{2l'+1}{3(l-m)!}\right\}\delta_{l'l}\right]
 \end{aligned}
 \end{equation}
 Eq. \eqref{eTempFullAmpl} can thus be written in a compact way like so,
 \begin{equation}\label{eTempFullAmpl-2}
 a_{lm}=\sum_{l'=0}^\infty b_{l'm}\Xi^m_{ll'}(\sigma)
 \end{equation}
 where
 \begin{equation}\label{eTempFullAmpl-3}
 \begin{aligned}
 \Xi^m_{ll'}(\sigma)=\delta_{ll'}-\sigma&\bigg[(\Pi^{-2})_{ll'}^m\delta_{l',l-2}+(\Pi^0)_{l'l}^m\delta_{l'l}+(\Pi^{+2})_{ll'}^m\delta_{l',l+2}\bigg]+\mathcal{O}(\sigma^2)
 \end{aligned}
 \end{equation}
 with
 \begin{equation}\label{eTempFullAmpl-3.5}
 \begin{aligned}
(\Pi^{+2})_{ll'}^m&=3\sqrt{\frac{(l-m)!(l'-m)!}{(2l+1)(2l'+1)}\cdot\frac{(l+m)!}{(l'+m)!}}\frac{(l'+m)(l+m)}{(2l+3)(l-m)!}\\
(\Pi^0)_{ll'}^m&=3\sqrt{\frac{(l-m)!(l'-m)!}{(2l+1)(2l'+1)}\cdot\frac{(l+m)!}{(l'+m)!}}\left\{\frac{(l'-m+1)(l+m)}{(2l+3)(l-m)!}+\frac{(l'+m)}{(2l-1)(l-m-1)!}-\frac{2l'+1}{3(l-m)!}\right\}\\
(\Pi^{-2})_{ll'}^m&=3\sqrt{\frac{(l-m)!(l'-m)!}{(2l+1)(2l'+1)}\cdot\frac{(l+m)!}{(l'+m)!}}\frac{(l'-m+1)}{(2l-1)(l-m-1)!}
 \end{aligned}
 \end{equation}
 Putting this in eq.\eqref{eTempFullAmpl-2}, we finally get,
 \begin{equation}\label{eTempFullAmpl-4}
 a_{lm}=b_{lm}-3\sigma\bigg[\Upsilon^{-2}_{lm} b_{l-2,m}+\Upsilon^0_{lm}b_{lm}+\Upsilon^{+2}_{lm} b_{l+2,m}\bigg]+\mathcal{O}(\sigma^2)
 \end{equation}
 where,
 \begin{equation}\label{eTempFullAmpl-5}
 \begin{aligned}
 \Upsilon^0_{lm}&=\frac{l+m}{2l+1}\left(\frac{l-m+1}{2l+3}+\frac{l-m}{2l-1}\right)-\frac{1}{3}\\
 \Upsilon^{+2}_{lm}&=\sqrt{\frac{(l+m)(l-m+1)(l-m+2)^2(l-m+3)}{(2l+1)(2l+3)^2 (2l+5)}}\\
 \Upsilon^{-2}_{lm}&=\sqrt{\frac{(l+m-1)(l+m)(l-m)^2(l-m-1)}{(2l-3)(2l-1)^2(2l+1)}}
 \end{aligned}
 \end{equation}
\end{tcolorbox}
 
 Let us henceforth replace $(\theta,\phi)$ by $\hat{\mathbf{n}}$ to simplify our notation. We will now calculate the ensemble average of the coefficients,
 \begin{equation}\label{eCoeffEnsAvgFull-1}
\left<a^*_{lm}a_{l'm'}\right>\equiv\int d\Omega_1\int d\Omega_2\, Y_{lm}(\hat{\mathbf{n}}_1)Y^*_{l'm'}(\hat{\mathbf{n}}_2)\left<\Theta(\hat{\mathbf{n}}_1)\Theta(\hat{\mathbf{n}}_2)\right>
 \end{equation}
Using eq. \eqref{eTempFlucFull}, the correlator inside the integral in eq. \eqref{eCoeffEnsAvgFull-1} can be written as,
\begin{equation}\label{eEnsAvgFull-1}
\begin{aligned}
\left<\Theta(\hat{\mathbf{n}}_1)\Theta(\hat{\mathbf{n}}_2)\right>&=\left[1-4\sigma\sqrt{\frac{\pi}{5}}Y_{20}(\hat{\mathbf{n}}_1)\right]\left[1-4\sigma\sqrt{\frac{\pi}{5}}Y_{20}(\hat{\mathbf{n}}_2)\right]\left<\Theta^{(\sigma=0)}(\hat{\mathbf{n}}_1)\Theta^{(\sigma=0)}(\hat{\mathbf{n}}_2)\right>+\mathcal{O}(\sigma^2)\\
&=\left[1-4\sigma\sqrt{\frac{\pi}{5}}\bigg\{Y_{20}(\hat{\mathbf{n}}_1)+Y_{20}(\hat{\mathbf{n}}_2)\bigg\}\right]\sum_{L=0}^\infty C_L^{(\sigma=0)}\sum_{m=-L}^L Y^*_{LM}(\hat{\mathbf{n}}_1)Y_{LM}(\hat{\mathbf{n}}_2)+\mathcal{O}(\sigma^2)
\end{aligned}
 \end{equation}
 where we have used,
 \begin{equation}\label{eEnsAvgFull-2}
\left<\Theta^{(\sigma=0)}(\hat{\mathbf{n}}_1)\Theta^{(\sigma=0)}(\hat{\mathbf{n}}_2)\right>=\sum_{L=0}^\infty\frac{2L+1}{4\pi}C_L^{(\sigma=0)}P_L(\hat{\mathbf{n}}_1\cdot \hat{\mathbf{n}}_2)=\sum_{L=0}^\infty C_L^{(\sigma=0)}\sum_{M=-L}^L Y^*_{LM}(\hat{\mathbf{n}}_1)Y_{LM}(\hat{\mathbf{n}}_2)
\end{equation}

Putting eq. \eqref{eEnsAvgFull-1} into eq. \eqref{eCoeffEnsAvgFull-1}, we get,
\begin{equation}\label{eCoeffEnsAvgFull-2}
\begin{aligned}
\left<a^*_{lm}a_{l'm'}\right>&=\sum_{L=0}^\infty\sum_{M=-L}^L C_L^{(\sigma=0)}\int d\Omega_1 \int d\Omega_2\, Y_{lm}(\hat{\mathbf{n}}_1) Y^*_{l'm'}(\hat{\mathbf{n}}_2) Y^*_{LM}(\hat{\mathbf{n}}_1)Y_{LM}(\hat{\mathbf{n}}_2) \\
&-4\sigma\sqrt{\frac{\pi}{5}}\sum_{L=0}^\infty\sum_{M=-L}^L C_L^{(\sigma=0)}\int d\Omega_1 \int d\Omega_2\, Y_{lm}(\hat{\mathbf{n}}_1) Y^*_{l'm'}(\hat{\mathbf{n}}_2) \bigg\{Y_{20}(\hat{\mathbf{n}}_1)+Y_{20}(\hat{\mathbf{n}}_2)\bigg\} Y^*_{LM}(\hat{\mathbf{n}}_1)Y_{LM}(\hat{\mathbf{n}}_2)+\mathcal{O}(\sigma^2) \\
&=\sum_{L=0}^\infty\sum_{M=-L}^L C_L^{(\sigma=0)}\delta_{Ll}\delta_{Mm}\delta_{Ll'}\delta_{Mm'}\\
& -4\sigma\sqrt{\frac{\pi}{5}}\sum_{L=0}^\infty\sum_{M=-L}^L C_L^{(\sigma=0)} \\
&\times\left[\delta_{Ll'}\delta_{Mm'}\int d\Omega_1\,Y_{lm}(\hat{\mathbf{n}}_1) Y_{20}(\hat{\mathbf{n}}_1)Y^*_{LM}(\hat{\mathbf{n}}_1) + \delta_{Ll}\delta_{Mm}\int d\Omega_2\,Y_{LM}(\hat{\mathbf{n}}_2) Y_{20}(\hat{\mathbf{n}}_2)Y^*_{l'm'}(\hat{\mathbf{n}}_2)\right]+\mathcal{O}(\sigma^2) \\
&=C_l^{(\sigma=0)} \delta_{l'l}\delta_{m'm}\\
&-4\sigma\sqrt{\frac{\pi}{5}}\left[C_{l'}^{(\sigma=0)}\int d\Omega_1\,Y_{lm}(\hat{\mathbf{n}}_1) Y_{20}(\hat{\mathbf{n}}_1)Y^*_{l'm'}(\hat{\mathbf{n}}_1) + C_{l}^{(\sigma=0)}\int d\Omega_2\,Y_{lm}(\hat{\mathbf{n}}_2) Y_{20}(\hat{\mathbf{n}}_2)Y^*_{l'm'}(\hat{\mathbf{n}}_2)\right] +\mathcal{O}(\sigma^2)\\
&=C_l^{(\sigma=0)} \delta_{l'l}\delta_{m'm} \\
&-\sigma\left\{\left(C_l^{(\sigma=0)}+C_{l'}^{(\sigma=0)}\right)\left[\frac{(2l+1)(2l'+1)(l-m)!(l'-m')!}{(l+m)!(l'+m')!}\right]^{\sfrac{1}{2}}\int_{-1}^1 dx\,P_2^0(x)\,P_{l'}^{m'}(x)\,P_l^m(x)\right\}\delta_{m'm}+\mathcal{O}(\sigma^2)
\end{aligned}
\end{equation}

Going forward, we can set $m=m'$ in our expressions. The last integral above is then known to us -- it is $\mathscr{I}_{ll'}^m$ of eq. \eqref{eThreePlm}, modulo a $\delta_{m'm}$ factor. We can thus write,
\begin{equation}\label{eCoeffEnsAvgFull-3}
\left<a^*_{lm}a_{l'm}\right>=C_l^{(\sigma=0)} \delta_{l'l}-\sigma \left(C_l^{(\sigma=0)}+C_{l'}^{(\sigma=0)}\right)\left[\frac{(2l+1)(2l'+1)(l-m)!(l'-m)!}{(l+m)!(l'+m)!}\right]^{\sfrac{1}{2}}\mathscr{I}_{ll'}^m+\mathcal{O}(\sigma^2)
\end{equation}
Using the equations in eq. \eqref{eTempFullAmpl-2} -- \eqref{eTempFullAmpl-3.5}, we can finally write this as,
\begin{equation}\label{eCoeffEnsAvgFull-4}
\tcboxmath{\left<a^*_{lm}a_{l'm}\right>=C_l^{(\sigma=0)} \delta_{l'l}-\sigma \left(C_l^{(\sigma=0)}+C_{l'}^{(\sigma=0)}\right) \bigg[(\Pi^{-2})_{ll'}^m\delta_{l',l-2}+(\Pi^0)_{l'l}^m\delta_{l'l}+(\Pi^{+2})_{ll'}^m\delta_{l',l+2}\bigg]+\mathcal{O}(\sigma^2)}
\end{equation}
This tells us,
\begin{equation}\label{eCoeffEnsAvgLimits}
\begin{aligned}
\left<a^*_{lm}a_{lm}\right>&=C_l^{(\sigma=0)} \bigg[1-2\sigma(\Pi^0)_{ll}^m\bigg]+\mathcal{O}(\sigma^2)\\
\left<a^*_{lm}a_{l\pm 2,m}\right>&=-\sigma \left(C_l^{(\sigma=0)}+C_{l\pm 2}^{(\sigma=0)}\right)(\Pi^{\pm 2})_{l,l\pm 2}^m+\mathcal{O}(\sigma^2)
\end{aligned}
\end{equation}
We can now make estimations on the magnitude of $\sigma$. Let us use the notation,
\begin{equation}\label{eDlSymbol}
\bar{D}_l\equiv\frac{1}{2l+1}\sum_{m=-l}^l\left<a^*_{l-1,m}a_{l+1,m}\right>
\end{equation}
Then from the latter of eq. \eqref{eCoeffEnsAvgLimits}, one obtains,
\begin{equation}\label{eSigmaLimOffDiag1}
\tcboxmath{\sigma=-\frac{\bar{D}_l}{\dfrac{1}{2l+1}\left[C_{l-1}^{(\sigma=0)}+C_{l+1}^{(\sigma=0)}\right]\ds{\sum_{m=-l}^l}(\Pi^{+2})_{l-1,l+1}^m}}
\end{equation} 
From Fig. 1 of \textbf{arXiv:astro-ph/0403695}, we estimate that $\bar{D}_{200}\sim 3.8\times 10^{-2}\,\mu\mathrm{K}^2$, while from the WMAP database\footnote{I am taking the WMAP data (\url{https://lambda.gsfc.nasa.gov/data/map/powspec/wmap_lcdm_pl_model_yr1_v1.txt}) rather than \textit{Planck} since \textbf{astro-ph/0403695} is based on it.}, we get $C_{200}^{(\sigma=0)}\sim 8.5\times 10^{-1}\,\mu\mathrm{K}^2$. In addition, $\sum_{m=-l}^l(\Pi^{+2})_{l-1,l+1}^m/(2l+1)\approx 0.5$. This gives an estimate of $\sigma\sim -4.4\times 10^{-2}$.



\end{document}