\section{Introduction}
Holography \cite{tHooft,Sussk} has an important role in the study of entanglement entropy \cite{Calab,Holzhey,Vidal} of quantum field theories.  An evidence for \emph{gauge/gravity duality} is the equality of the modular Hamiltonian \cite{Haag-bk,Casini-1} in the field theory side and the asymptotic isometries of the gravity in the bulk. The modular Hamiltonian is defined through the density matrix of the theory. This Hamiltonian is charge of symmetry group of the system that transforms any operator to itself \cite{Basu}.

The modular Hamiltonian has been used to derive the holographic entanglement entropy in the AdS/CFT correspondence \cite{Casini-2}. Another importance of the modular Hamiltonian is in the derivation of the first law of entanglement \cite{Blanco}. Given importance of this Hamiltonian in AdS/CFT correspondence, we would like to study its role in non-AdS holography. So far the vacuum modular flow generator and the corresponding modular Hamiltonian have been determined in a specific class of non-AdS holographies. In the $Flat_3/BMSFT$ holography as a case of holography beyond AdS/CFT, it is considered that gravity in asymptotically flat three-dimensional spacetimes is dual to a BMS-invariant field theory (BMSFT) \cite{Bagchi-1,Bagchi-5}. In this class of the holography, the flow generated by modular Hamiltonian differs from the modular flow in AdS/CFT \cite{Jiang,Apolo}.

In this paper we consider flat holography in the framework of generalized minimal massive gravity model (GMMG) model. The generalized minimal massive gravity model \cite{Setare 2014} that is realized by adding the CS deformation term, the higher derivative
deformation term, and an extra term to pure Einstein gravity with a negative cosmological constant, is considered in the bulk side of the holography.
 Usually the theories including the terms given by the square of the curvatures have the massive spin
2 mode and the massive scalar mode in addition to the massless graviton. Also the theory has ghosts
due to negative energy excitations of the massive tensor. In \cite{Setare 2014} it is discussed that the GMMG is free
of negative-energy bulk modes, and also avoids the aforementioned
\emph{bulk-boundary unitarity clash}. By a Hamiltonian
analysis one can show that the GMMG model
has no Boulware-Deser ghosts and this model propagate
only two physical modes.
 The group of asymptotic symmetries of asymptotically flat spacetimes at the
future null infinity is the BMS group \cite{Bondi-1,Sach-1,Sach-2}. This symmetry is in correspondence with the symmetry group of the BMSFTs.
An evidence for this holography is equality of the gravitational charge that is related to the asymptotic isometries of the asymptotically flat three dimensional spacetime in GMMG in bulk and the modular Hamiltonian on the field theory side. To prove this correspondence, we find the conserved charge associated with asymptotic isometries of asymptotically flat solutions of  GMMG model. To this end, we consider the metric of asymptotically flat three dimensional spacetime in Bondi gauge. We focus on three special cases of this metric. To compute the modular Hamiltonian we consider BMSFTs in thermal and non-thermal states. Using the generalized Rindler transformation, we compute the modular flow generator \cite{Casini-2,Wen,Jiang,Apolo,Castro} (another method for producing the local modular flow without using the explicit form of the Rindler transformation is proposed in \cite{Wen2}). In this general prescription of the Rindler transformations, some conditions are made on the symmetry transformations $x\to \tilde{x}$ of the coordinates in the field theory. This transformation is invariant under a thermal identification $\tilde{x}^i\sim \tilde{x}^i+i\tilde{\beta}^i$. The vacuum is invariant along this translation and the vacuum state is mapped to a thermal state in the Rindler spacetime.

This paper is organized as follows. In section \ref{sec:2} we introduce briefly the GMMG model. Using the Bondi fall-off conditions, we write the Killing vectors that preserve these conditions. In section \ref{sec:3} we describe the modular Hamiltonian and its relation with the vacuum modular flow generator. Then we find the Killing vectors in two special cases of the asymptotically flat spacetime. In the filed theory side we find the modular Hamiltonian in thermal and non-thermal BMSFTs. We find the first law of holographic entanglement in these cases. In section \ref{sec:4} we provide a summary of the paper and give an overview on this class of $Flat_3/BMSFT$ holography.

\section{Generalized minimal massive gravity in asymptotically flat spaces }\label{sec:2}
Generalized minimal massive gravity (GMMG) is an example of the Chern-Simons-like theories of gravity.
By adding the Chern-Simons term, the higher derivative and an extra term to pure Einstein gravity, the Lagrangian \emph{3-form} of GMMG is given by \cite{Setare 2014}
\begin{eqnarray}\label{L-GMMG}
  L_{GMMG} &=& -\sigma e_aR^a+\frac{\Lambda_0}{6}\epsilon^{abc}e_ae_be_c+h_aT^a+\frac{1}{2\mu}[\omega_ad\omega^a+\frac{1}{3}
  \epsilon^{abc}\omega_a\omega_b\omega_c]\nonumber \\
   &-& \frac{1}{\mathrm{m}^2}[f_aR^a+\frac{1}{2}\epsilon^{abc}e_af_bf_c]+\frac{\alpha}{2}\epsilon^{abc}e_ah_bh_c
\end{eqnarray}

There are four fields that form four flavours of \emph{one-form}, $a^r=\{e,~\omega,~h,~f\}$.  $e$ is dreibein, $h$ and $f$ are the auxiliary \emph{one-form} field and $\omega$ is dualised spin-connection. $T( \omega)$ is Lorentz covariant torsion and $R( \omega)$ is curvature \emph{2-form}. $\mathrm{m}$ is mass parameter of the new massive gravity (NMG) term and $\sigma$ is a sign. The last term of the first line of (\ref{L-GMMG})is the \emph{Lorentz Chern-Simons} term that $\mu$ is mass parameter of that term. $\Lambda_0$ that has dimension mass squared is a cosmological parameter and $\alpha $ is a dimensionless parameter. The last term in the Lagrangian of GMMG (\ref{L-GMMG}) is an extra  term that is added to the model.

The following fall-off conditions is considered for the asymptotically $2+1$-dimensional flat spacetime \cite{Setare-Adami Nucl 2018,Detournay 2017}
\begin{eqnarray}\label{B.C.}
  g_{uu} &=& \mathcal{M}(\phi)+\mathcal{O}(r^{-2}) \nonumber\\
  g_{ur} &=& -e^{\mathcal{A}(\phi)}+\mathcal{O}(r^{-2}) \nonumber\\
  g_{u\phi} &=& \mathcal{N}(u,\phi)+\mathcal{O}(r^{-1}) \nonumber\\
  g_{rr} &=& \mathcal{O}(r^{-2}) \nonumber\\
  g_{r\phi} &=& -e^{2\mathcal{A}(\phi)}\mathcal{E}(u,\phi)+\mathcal{O}(r^{-1}) \nonumber\\
  g_{\phi\phi} &=& e^{\mathcal{A}(\phi)}r^2+\mathcal{E}(u,\phi)[2\mathcal{N}(u,\phi)-
  \mathcal{M}(\phi)\mathcal{E}(u,\phi)]+\mathcal{O}(r^{-1})
\end{eqnarray}
where
\begin{eqnarray}\label{N,E}
  \mathcal{N}(u,\phi) &=& \mathcal{L}(\phi)+\frac{u}{2}\partial_{\phi}\mathcal{M}(\phi)\nonumber \\
  \mathcal{E}(u,\phi) &=& \mathcal{B}(\phi)+u\partial_{\phi}\mathcal{A}(\phi)
\end{eqnarray}

In the above equations $\mathcal{A}(\phi)$, $\mathcal{L}(\phi)$, $\mathcal{B}(\phi)$ and $\mathcal{M}(\phi)$ are arbitrary functions and $r$, $u$ and $\phi$ are the coordinates. The conditions admit the Bondi-Sachs  fall-off conditions at future null infinity. The metric transformations generated by vector field $\xi$ are obtained by $\delta_{\xi}g_{\mu\nu}=\mathcal{\L}_{\xi}g_{\mu\nu}$. $\mathcal{\L}_{\xi}$ is the Lie derivative along $\xi$. The Killing vector fields that preserve the boundary conditions invariant are as follows \cite{Setare-Adami Nucl 2018,Detournay 2017}
\begin{eqnarray} \label{Killing GMMG}
   \xi^u & = & \alpha(u,\phi)-\frac{1}{r}e^{-\mathcal{A}(\phi)}\mathcal{E}(u,\phi)\beta(u,\phi)+\mathcal{O}(r^{-2})\nonumber \\
    \xi^r &=& r X(\phi)+e^{-\mathcal{A}(\phi)}[\mathcal{E}(u,\phi)\partial_{\phi}X(\phi)-\partial_{\phi}\beta(u,\phi)]\nonumber \\
     & +& \frac{1}{r}e^{-2\mathcal{A}(\phi)}\beta(u,\phi)[\mathcal{N}(u,\phi)-\mathcal{M}(\phi)\mathcal{E}(u,\phi)]
     +\mathcal{O}(r^{-2}),\nonumber \\
     \xi^{\phi} &=& Y(\phi)+\frac{1}{r}e^{-\mathcal{A}(\phi)}\beta(u,\phi)+\mathcal{O}(r^{-2}),
\end{eqnarray}
where
\begin{eqnarray}\label{alfa-beta}
  \alpha(u,\phi) &=& T(\phi)+u\partial_{\phi}Y(\phi), \nonumber \\
  \beta(u,\phi) &=& Z(\phi)+u\partial_{\phi}X(\phi), \nonumber \\
   \beta(u,\phi) &=& - \partial_{\phi}\alpha(u,\phi)
\end{eqnarray}

$T(\phi)$, $Y(\phi)$, $Z(\phi)$ and $X(\phi)$ are arbitrary functions of $\phi$ that are related by (\ref{alfa-beta}). The conserved charges of asymptotically flat spacetime corresponds to the Killing vector field (\ref{Killing GMMG}) are obtained in \cite{Setare-Adami Nucl 2018}. Using the Fourier modes of the charges, it has been shown the asymptotic symmetry algebra is a semidirect product of a $BMS_3$ algebra and two $U(1)$ current algebras \cite{Setare-Adami Nucl 2018}. These charges are not only for the Killing vector fields are satisfied by spacetime everywhere but also these are admitted for the Killing vectors in asymptotic configuration.
The authors in \cite{Setare-Adami Nucl 2018} have used the quasi-local conserved charge perturbation associated with a field dependent vector field $\xi$ in GMMG model for asymptotically flat spacetimes and they have used an integration over the one-parameter path on the solution space from the charge perturbation to find the following conserved charges
\begin{equation}\label{Total charge}
  Q(T,X,Y,Z)=M(T)+J(X)+L(Y)+P(Z),
\end{equation}
where

\begin{eqnarray}
  M(T) &=& -\frac{1}{16\pi G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\int_{0}^{2\pi}T(\phi)\mathcal{M}(\phi)
  d\phi, \label{M- charge}\\
  J(X) &=& \frac{1}{8\pi G}\int_{0}^{2\pi}X(\phi)\big[(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})
  \partial_{\phi}\mathcal{B}(\phi)-\frac{1}{2\mu}\partial_{\phi}\mathcal{A}(\phi)\big] d\phi,\label{J-charge} \\
  L(Y) &=& -\frac{1}{8\pi G} \int_{0}^{2\pi}Y(\phi)\bigg\{(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\mathcal{L}(\phi)\nonumber    \\
   &-& \frac{1}{4\mu}\big[2\mathcal{M}(\phi)+(\partial_{\phi}\mathcal{A}(\phi))^2-2\partial_{\phi}^2\mathcal{A}(\phi) \big] \bigg\}d\phi, \label{L-charge} \\
   P(Z) &=& \frac{1}{8\pi G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\int_{0}^{2\pi}Z(\phi)\partial_{\phi}\mathcal{A}(\phi)d\phi. \label{P-charge}
\end{eqnarray}

$H$ and $F$ are constant parameters.

\section{Flat holography in the framework of GMMG }\label{sec:3}
The BMS group is the symmetry group of the asymptotically flat spacetime in three dimensions \cite{Bondi-1,Sach-1,Sach-2,Barni-Comp}. This group is an
infinit-dimensional one in $3$ and $4$-dimensions . The BMS algebra in three dimensions is the same as the symmetry algebra of a BMS-invariant field theory (BMSFT) in two dimensions. This motivate us to a $Flat_3/BMSFT$ duality \cite{Bagchi-1,Barni-Troes}.
An evidence for this correspondence is computation of the entanglement entropy \cite{Bagchi-2,Basu}. The entanglement entropy of 2d BMSFT at the boundary is the same as the entropy in the gravity side. Another evidence for the correspondence is reproducing the Bekenstein-Hawking entropy of Flat space cosmologies using a modified Cardy formula for the BMSFTs \cite{Barni-1,Bagchi-3}. Anyone who is interested could see in \cite{Bagchi-4} to a review a number of applications of this holography.

In our previous work \cite{we}, we considered the GMMG model in the asymptotically flat space in three cases: a flat space cosmological geometry (FSC), a null orbifold and a Poincare patch. The 2d field theory on the boundary was considered as a BMSFT. Using the Cardy formula, the entaglement entropy of the field theory side was computed and was compared with the Beckenestein-Hawking entropy of the garvity side. The entropy in two sides of the GMMG/BMSFT correspondence is another evidence for $Flat_3/BMSFT$ holography conjecture. In this section we study the modular Hamiltonian in both sides of the correspondence and derive the first law of the entanglement entropy.

The modular Hamiltonian has an important role in holography and black hole physics. On a subregion $\mathcal{A}$ of theory in the boundary, the density matrix $\rho_{\mathcal{A}}$ is described through its modular Hamiltonian as follows \cite{Haag-bk,Casini-1}
\begin{equation}\label{H mod : int}
  \mathcal{H}_{mod}=-\log\rho_{\mathcal{A}}
\end{equation}
where the modular Hamiltonian is generically independent of some local expressions constructed with the fields on $\mathcal{A}$. This Hamiltonian generically is not a local operator \cite{Casini-2}. The modular Hamiltonian is charge of a symmetry of the system through the unitary operator $U(s)=e^{-i\mathcal{H}_{mod}s}$. This symmetry group transforms any operator in the causal domain of $\mathcal{A}$ into itself as follows \cite{Haag-bk,Casini-2,Boer}
\begin{equation}\label{H mod:cp}
  e^{i\mathcal{H}_{mod}s}\mathcal{O}(\mathcal{A})e^{-i\mathcal{H}_{mod}s}=\mathcal{O}(\mathcal{A})
\end{equation}
where $\mathcal{O}(\mathcal{A})$ is the algebra of operators in $\mathcal{A}$. This internal time transformations on the charge algebra in $\mathcal{A}$ is a one-parameter group. This group is called the modular flow corresponding to $\mathcal{A}$. Based on the generalized Rindler method \cite{Wen,Jiang}, the modular flow generator can be defined as a linear combination of the vacuum symmetry generators $h_i$ as follows
\begin{equation}\label{flow gen:int}
  \zeta=\Sigma_i a_ih_i
\end{equation}
where $a_i$ are some arbitrary constants and the relation (\ref{flow gen:int}) is just true for the local modular flows. The vacuum symmetry generators $h_i$ preserve the vacuum of the field theory invariant. Corresponding to the dictionary of holography, the isometry group of the gravity theory in the bulk $H$ is dual to the isometry group of the vacuum state of the field theory on the boundary. There is the following correspondence between the vacuum symmetry generators $h_i$ and the Killing vectors $H_i$ of diffeomorphism in the bulk
\begin{equation}\label{H=h}
  H_i|_{\partial M}=h_i
\end{equation}
where $H_i|_{\partial M}$ are defined on a surface $\partial M$. According to (\ref{H=h}), we find the following correspondence between the modular flow generator in the bulk and the modular flow generator on the boundary
\begin{equation}\label{xi=zeta}
  \xi|_{\partial M}=\zeta
\end{equation}
where
\begin{equation}\label{xi:gen}
\xi=\Sigma_ia_iH_i.
\end{equation}

The coefficients $a_i$ in (\ref{xi:gen}) are the same as the parameters in (\ref{flow gen:int}). In the Bondi gauge, we can re-parameterize the fall-off conditions (\ref{B.C.}) as follows \cite{Barni-Troes}
\begin{equation}\label{flat space-2}
  ds^2=Mdu^2-2dudr+Jdud\phi+r^2d\phi^2
\end{equation}
where $M$ and $J$ can be regarded as the mass and the angular momentum of the spacetime respectively. Comparing Eqs. (\ref{flat space-2}) and (\ref{B.C.}), we have the following identification
\begin{eqnarray} \label{Ids. 1}
  \mathcal{M}(\phi) &=& M \nonumber\\
  \mathcal{A} &=& 0 \nonumber\\
  2\mathcal{N}(u,\phi) &=& J\nonumber \\
  \mathcal{B} &=& 0\nonumber \\
  \mathcal{E}(u,\phi) &=& 0
\end{eqnarray}

We are especially interested in three cases corresponding to the parameters $M$ and $J$ values:

(i) $M=-1$ and $J=0$: In this case the background spacetime is converted to the Minkowski spacetime in three dimensions (called global Minkowski slotuion). This metric has the following form
\begin{equation}\label{Mink}
  ds^2=-du^2-2dudr+r^2d\phi^2
\end{equation}

(ii) $M=J=0$: The metric has the following form
\begin{equation}\label{Poin}
  ds^2=-2dudr+r^2d\phi^2,
\end{equation}
that is the Poincare patch in the $AdS_3$ spacetime. The manifold is a null orbifold that is defined by a null boost \cite{Horowitz}.

(iii)$M>0$: In this case we have a flat space cosmological geometry (FSC) that is
time-dependent solution of Einstein equation in three dimensions of spacetime \cite{we,FSC-1,FSC-2}.

Using the identifications (\ref{Ids. 1}), the Killing vector fields (\ref{Killing GMMG}) and (\ref{alfa-beta}) that preserve the Bondi boundary conditions invariant are as follows
\begin{eqnarray} \label{Killing GMMG-2}
  \xi^u &=& T(\phi)+u\partial_{\phi}Y(\phi)\nonumber \\
  \xi^r &=& -\frac{J}{2r}\partial_{\phi}T(\phi)-\frac{J}{2r}u\partial_{\phi}^2Y(\phi)+u\partial_{\phi}^3Y(\phi)+\partial_{\phi}^2T(\phi)
  -r\partial_{\phi}Y(\phi)\nonumber \\
  \xi^{\phi} &=& Y(\phi)- \frac{1}{r} \partial_{\phi}T(\phi)-\frac{u}{r}\partial_{\phi}^2Y(\phi),
\end{eqnarray}

 By substituting $\mathcal{A}=0$ and $\mathcal{B}=0$ into the conserved charges GMMG model (\ref{M- charge}), (\ref{J-charge}), (\ref{L-charge}) and (\ref{P-charge}) for asymptotically flat spacetime
 we have the following form for the charges
 \begin{eqnarray}
  M(T) &=& -\frac{1}{16\pi G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\int_{0}^{2\pi}T(\phi)\mathcal{M}(\phi)
  d\phi, \label{M- charge2}\\
  J(X) &=& 0,\label{J-charge2} \\
  L(Y) &=& -\frac{1}{8\pi G} \int_{0}^{2\pi}Y(\phi)\bigg\{(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\mathcal{L}(\phi)-\frac{1}{2\mu}\mathcal{M}(\phi)\bigg\}d\phi \label{L-charge2} \\
   P(Z) &=& 0 \label{P-charge2}
\end{eqnarray}
By adding these results together we obtain following expression
\begin{eqnarray} \label{F-charge GMMG}
  Q &=& M(T)+J(X)+L(Y)+P(Z)\nonumber \\
   &=&  -\frac{1}{16\pi G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\int_{0}^{2\pi}\bigg\{T(\phi)\mathcal{M}(\phi)+2Y(\phi)\mathcal{L}(\phi)\bigg\}d\phi\nonumber \\
   &+& \frac{1}{16\pi G \mu}\int_{0}^{2\pi}Y(\phi)\mathcal{M}(\phi)d\phi
\end{eqnarray}

This is the gravitational charge corresponds to the modular flow generator (\ref{Killing GMMG-2}). The infinitesimal variation corresponding to this charge is as follows
\begin{eqnarray} \label{delta F-charge GMMG}
  \delta Q &=& \delta M(T)+\delta J(X)+\delta L(Y)+\delta P(Z)\nonumber \\
   &=&  -\frac{1}{16\pi G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\int_{0}^{2\pi}\bigg\{T(\phi)\delta \mathcal{M}(\phi)+2Y(\phi)\delta \mathcal{L}(\phi)\bigg\}d\phi\nonumber \\
   &+& \frac{1}{16\pi G \mu}\int_{0}^{2\pi}Y(\phi)\delta \mathcal{M}(\phi)d\phi
\end{eqnarray}

It is interesting and as one expect, the above result reduces to the expression (3.51)  in the paper \cite{Apolo} for Einstein gravity in $2+1$ dimension, when $\mu , \mathrm{m}^2 \to \infty $ and $\sigma=-1$. In these limiting cases, GMMG model reduces to the Einstein gravity. We use this result to show an evidence for the $flat_3/BMSFT$ correspondence. In the field theory side, the story is a little different. To compute the vacuum modular Hamiltonian for a single interval in BMSFTs on the plane, we first find the modular flow generator from the generalized Rindler transformation \cite{Casini-2,Jiang,Apolo,Castro}. In this prescription a symmetry transformation of the vacuum state is used to map an interval $\mathcal{I}$ (and its causal domain $\mathcal{D}$) to a Rindler spacetime. If \textbf{G} is a symmetry group that preserves the field theory invariant, a subgroup $\textbf{H}\subset\textbf{G}$ is a symmetry group that preserves the vacuum state invariant. In the Rindler spacetime, the generator of time-like translations is identified with the modular Hamiltonian \cite{Castro,Apolo}.

In a generalized Rindler prescription \cite{Castro,Apolo}, it could be possible to find the modular flow generator without an explicit expression for the Rindler transformations. For this end, the modular flow generator is considered as a linear combination of the vacuum symmetry generators similar to \ref{flow gen:int}. There is a condition to determine the coefficients $a_i$ and that is the invariance of the causal domain $\partial\mathcal{D}$ under the modular flow generator. To determine the overall normalization of the modular flow generator, it should be assumed that $e^{i \zeta}$ maps any point back to itself.

To apply the generalized Rindler method, we consider an interval $\mathcal{I}$ on the vacuum state as follows
\begin{equation}\label{interval}
  \partial\mathcal{I}=\{(u_-\phi_-),(u_+,\phi_+)\}
\end{equation}
where $u$ and $\phi$ are the coordinates on the plane and the interval $\mathcal{I}$ is identified with its endpoints. The modular flow generator as linear combination of vacuum symmetry generators $l_j$ and $m_j$ can be write as follows
\begin{equation}\label{zeta-1}
  \zeta=\zeta^u\partial_u+\zeta^\phi\partial_\phi=\Sigma_{j=-1}^1(a_jl_j+b_jm_j),
\end{equation}
where $l_j$ and $m_j$ are defined as follows
\begin{eqnarray}
  l_j &=& -\phi^{j+1}\partial_\phi-(j+1)\phi^ju\partial_u,\nonumber \\
  m_j &=& -\phi^{j+1}\partial_u \label{l&m1}
\end{eqnarray}
and satisfy the BMS algebra
\begin{eqnarray}
  [l_j,l_k] &=& (j-k)l_{j+k}, \\ \label{l&l}
  [l_j,m_k] &=& (j-k)m_{j+k}, \\ \label{l&m2}
  [m_j,m_k] &=& 0 \label{m&m}
\end{eqnarray}

To determine the coefficients $a_j$ and $b_j$, it is required to assume the boundary of the causal domain $\partial\mathcal{D}$ is invariant under the modular flow. By this requirement, the modular flow generator vanishes at the endpoints of the interval $\mathcal{I}$ (\ref{interval}) as follows
\begin{eqnarray} \label{zeta right}
  \zeta_{(u_+,\phi_+)} &=& 0 \nonumber \\
  &=& a_+l_++a_0l_0+a_-l_-+b_+m_++m_0b_0+b_-m_-
\end{eqnarray}
 \begin{eqnarray} \label{zeta left}
  \zeta_{(u_-,\phi_-)} &=& 0 \nonumber \\
  &=& a_+l_++a_0l_0+a_-l_-+b_+m_++m_0b_0+b_-m_-
\end{eqnarray}

Substituting $l_j$ and $m_j$ definitions (\ref{l&m1}) into the above equations, we find four relations for $a_j$ and $b_j$ as follows
\begin{eqnarray}\label{coef}
  a_+\phi_+^2+a_0\phi_++a_- &=& 0,  \\
  a_+\phi_-^2+a_0\phi_-+a_- &=& 0, \label{b} \\
  2a_+\phi_+u_++a_0u_++b_+\phi_+^2+b_0\phi_++b_- &=& 0, \label{c} \\
  2a_+\phi_-u_-+a_0u_-+b_+\phi_-^2+b_0\phi_-+b_- &=& 0. \label{d}
\end{eqnarray}
These relations are corresponding to the coefficients of $\partial_u$ and $\partial_\phi$ when the $l_j$ and $m_j$ definitions (\ref{l&m1}) is substituted into the equations (\ref{zeta right}) and (\ref{zeta left}). Using Eqs. (\ref{coef}) and (\ref{b}), we have the following expressions for the coefficients $a_j$
\begin{equation}\label{aj}
  (a_+,a_0,a_-)=a_+\big(1,-(\phi_++\phi_-),\phi_+\phi_-\big).
\end{equation}
The solution of the Eqs. (\ref{c}) and (\ref{d}) for the coefficients $b_j$ is as follows
\begin{equation}\label{bj}
  (b_+,b_0,b_-)=b_+\big(1,-(\phi_++\phi_-),\phi_+\phi_-\big)+a_+\big(0,-(u_-+u_+),u_+\phi_-+u_-\phi_+\big)
\end{equation}
To find $a_+$ and $a_-$, the following differential equations are considered \cite{Apolo}
\begin{equation}\label{traject}
  \partial_su(s)=\zeta^u,~~~~\partial_s\phi(s)=\zeta^\phi,
\end{equation}
where $(\phi(s),u(s))$ is the trajectory under the modular flow $e^{s\zeta}$. By solving the differential equations (\ref{traject}) and applying the conditions $\phi(s)=\phi(s+i)$ and $u(s)=u(s+i)$, $a_+$ and $a_-$ are found as follows  \cite{Apolo}
\begin{eqnarray}\label{ab+}
  a_+ &=& \frac{2\pi}{\phi_+-\phi_-},\nonumber \\
  b_+ &=& -\frac{2\pi(u_+-u_-)}{(\phi_+-\phi_-)^2}.
\end{eqnarray}

Substituting the results for $a_+$ and $a_-$ into eqs. (\ref{aj}) and (\ref{bj}), the coefficients are computed as follows
\begin{eqnarray}\label{aj2}
  (a_+,a_0,a_-) &=& \frac{2\pi}{\phi_+-\phi_-}\big(1,-(\phi_++\phi_-),\phi_+\phi_-\big) \\
  (b_+,b_0,b_-) &=& \frac{2\pi}{(\phi_+-\phi_-)^2}\big(u_--u_+,2u_+\phi_--2u_-\phi_+,u_-\phi_+^2-u_+\phi_-^2\big)\label{bj2}
\end{eqnarray}

By the coefficients (\ref{aj2}) and (\ref{bj2}), we have the modular flow generator $\zeta$. We now can compute the vacuum modular Hamiltonian for an interval $\mathcal{I}$ in BMSFT on the plane. To find this Hamiltonian, the vacuum symmetry generators $l_j$ and $m_j$ should be replaced with the corresponding conserved charges $\mathcal{L}_j$ and $\mathcal{M}_j$. These charges are defined as follows

\begin{eqnarray}
  \mathcal{L}_j &=& \frac{1}{2\pi}\int d\phi \phi^{j+1}\mathcal{J}(\phi),\nonumber \\
  \mathcal{M}_j &=& \frac{1}{2\pi}\int d\phi \phi^{j+1}\mathcal{P}(\phi) \label{chs}
\end{eqnarray}
where $\mathcal{J}(\phi)$ and $\mathcal{P}(\phi)$ are the conserved currents associated with $\mathcal{L}_j$ and $\mathcal{M}_j$, respectively.

Replacing the symmetry generators $l_j$ and $m_j$ with the charges (\ref{chs}) in (\ref{zeta-1}), we find the modular Hamiltonian as follows
\begin{eqnarray}\label{H in FT}
  \mathcal{H}_{\zeta} &=& \Sigma_{j=-1}^1(a_j\mathcal{L}_j+b_j\mathcal{M}_j)\nonumber \\
   &=& -\frac{1}{2\pi}\int_{\phi_-}^{\phi_+}(T(\phi)\mathcal{P}(\phi)+Y(\phi)\mathcal{J}(\phi))
\end{eqnarray}
where by substituting the coefficients (\ref{aj2}) and (\ref{bj2}) into the first line of the above equation, $T(\phi)$ and $Y(\phi)$ are found as follows
\begin{eqnarray}\label{T}
  T(\phi) &=& \frac{2\pi[u_+(\phi-\phi_-)^2-u_-(\phi-\phi_+)^2]}{(\phi_+-\phi_-)^2} \\
  Y(\phi) &=& -\frac{2\pi(\phi-\phi_-)(\phi-\phi_+)}{\phi_+-\phi_-} \label{Y}
\end{eqnarray}
where the functions $T(\phi)$ and $Y(\phi)$ are non-vanishing only within the causal domain $\mathcal{D}$. So far we have calculated the modular Hamiltonian on the BMSFT on a plane in (\ref{H in FT}) and the gravitational charge corresponding to the modular flow generator in (\ref{F-charge GMMG}). A consistency check on these results is that these results satisfy the first law of entanglement entropy. We now show that the modular Hamiltonian (\ref{H in FT}) is the same as the gravitational charge (\ref{F-charge GMMG}). To this end, we compare the charge algebras on two sides of the holography. By introducing the following Fourier modes of the conserved charges (\ref{M- charge2}), (\ref{J-charge2}), (\ref{L-charge2}), (\ref{P-charge2}) and (\ref{F-charge GMMG}) of the asymptotically flat spacetime in GMMG
\begin{eqnarray}\label{fourier1}
  M_n &=& Q(e^{in\phi},0,0,0)=M(e^{in\phi}),\nonumber \\
  J_n &=& Q(0,e^{in\phi},0,0)=J(e^{in\phi}),\nonumber \\
  L_n &=& Q(0,0,e^{in\phi},0)=L(e^{in\phi}),\nonumber \\
  P_n &=& Q(0,0,0,e^{in\phi})=P(e^{in\phi}),
\end{eqnarray}
we find the non-zero Fourier modes as follows
\begin{eqnarray}
  M_n(T) &=& -\frac{1}{16\pi G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2 })\int_{0}^{2\pi}T(\phi)\mathcal{M}(\phi)e^{in\phi}
  d\phi, \label{Mm- charge2}\\
  L_n(Y) &=& -\frac{1}{8\pi G} \int_{0}^{2\pi}Y(\phi)\bigg\{(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\mathcal{L}(\phi)-\frac{1}{2\mu}\mathcal{M}(\phi)\bigg\}e^{in\phi}d\phi \label{Lm-charge2}
\end{eqnarray}
where $m$ indicates the Fourier modes and $\mathrm{m}$ is the mass parameter of NMG.

It could be proved that the gravitational charges (\ref{Mm- charge2}) and (\ref{Lm-charge2}) satisfy three-dimensional BMS algebra as follows \cite{Setare-Adami Nucl 2018}
\begin{eqnarray}\label{BMS algebra}
  \[L_m,L_n\] &=& (m-n)L_{m+n}+\frac{c_L}{12}m^3\delta_{m+n,0}\nonumber \\
  \[M_m,L_n\] &=& (m-n)M_{m+n}+\frac{c_M}{12}m^3\delta_{m+n,0}\nonumber \\
  \[M_m,M_n\] &=& 0
\end{eqnarray}
where the $c_L$ and $c_M$ are the central extensions as follows
\begin{eqnarray}\label{cL,cM}
  c_L &=& \frac{3}{G\mu}\nonumber \\
  c_M &=& -\frac{3}{G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})
\end{eqnarray}
The Fourier modes of the $\mathcal{J}(\phi)$ and $\mathcal{P}(\phi)$ currents (\ref{chs}) are as follows
\begin{eqnarray}\label{fourier2}
  \mathcal{L}_m &=& -\frac{1}{2\pi}\int d\phi e^{im\phi}\mathcal{J}(\phi),\nonumber \\
  \mathcal{M}_m &=& -\frac{1}{2\pi}\int d\phi e^{im\phi}\mathcal{P}(\phi)
\end{eqnarray}
These charges satisfy the same BMS algebra of the asymptotically flat spacetime in GMMG (\ref{BMS algebra}) where
\begin{equation}\label{charge dict1}
  L_m\to \mathcal{L}_m, ~~~~~ M_m\to \mathcal{M}_m
\end{equation}

This correspondence between the bulk and the boundary symmetry algebras is an indication of $flat_3/BMSFT$ holography. By this correspondence, we find the following dictionary between the bulk and the boundary variables
\begin{eqnarray}\label{var.rel}
  \frac{\mathcal{P}}{2\pi} &=& \frac{1}{16\pi G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\mathcal{M}(\phi)\nonumber \\
  \frac{\mathcal{J}}{2\pi} &=& \frac{1}{8\pi G}\bigg\{(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\mathcal{L}(\phi)-\frac{1}{2\mu}\mathcal{M}\bigg\}
\end{eqnarray}

Using this holographic dictionary, the gravitational charge (\ref{F-charge GMMG}) is the same as the modular Hamiltonian (\ref{H in FT}) on the field theory side as follows
\begin{equation}\label{match1}
  \delta Q\[g\]=\delta \langle\mathcal{H}_{\zeta}\rangle=\delta \langle\mathcal{H}_{mod}\rangle
\end{equation}

In \cite{Apolo-2}, the authors have proposed a general prescription for matching the holographic entanglement entropy and the gravitational charge evaluated on an asymptotically flat space. So we have the following result for the first law of entanglement entropy in $flat_3/BMSFT$
\begin{equation}\label{match1}
  \delta S_{EE}=\delta Q\[g\]=\delta \langle\mathcal{H}_{\zeta}\rangle=\delta \langle\mathcal{H}_{mod}\rangle
\end{equation}

To find this result, we have considered the Poincare patch (\ref{Poin}) in the gravity side of the holography. In the field theory side, we have had a BMSFT with zero temperature on the plane (\ref{interval}). We could generalize the result to more general cases in the bulk and boundary of the $flat_3/BMSFT$ correspondence. In the boundary, we consider the thermal BMSFT that is defined on a thermal cylinder that has the following thermal identification of the coordinates
\begin{equation}\label{therm. cyl}
  (u,\phi)\sim (u,,\phi+2\pi)\sim (u+i\beta_u,\phi-i\beta_{\phi})
\end{equation}
where
\begin{equation}\label{beta1}
  \beta_u=\frac{\pi J}{M^{3/2}},~~~~ \beta_{\phi}=\frac{2\pi}{\sqrt{M}}.
\end{equation}

Now, we consider the Flat Space Cosmological solution (FSC) that is third case of general asymptotically flat space time solution of GMMG, this case is considered with the following metric,
\begin{equation}\label{FSC}
  ds^2=Mdu^2-2dudr+Jdud\phi+r^2d\phi^2, ~~~~ M>0.
\end{equation}

To find the modular flow generator, the interval $\mathcal{I}$ at the boundary is parameterized
 in terms of the endpoints as follows
 \begin{equation}\label{therm. cyl.2}
   \partial\mathcal{I}=\{(u_-,\phi_-),(u_+,\phi_+)\},~~~~~ l_u\equiv u_+-u_-,~~l_{\phi}\equiv\phi_+-\phi_-.
 \end{equation}

The generalized Rindler method could be used to find the modular flow generator. The modular flow generator can be written as follows
\begin{eqnarray}\label{zeta-therm}
  \zeta &=& \zeta^u\partial_u+\zeta^\phi\partial_\phi=\Sigma_{j=-1}^1a_jl_j+b_jm_j\nonumber\\
  &=& [T(\phi)+uY'(\phi)]\partial_u+Y(\phi)\partial_\phi
\end{eqnarray}

To write the $2nd$ line of the above equation, we use the following Killing vectors of FSC in $r\to\infty$ limit \cite{Jiang}
\begin{eqnarray}\label{Kil-FSC-2}
  l_1 &=& -e^{-\sqrt{M}\phi}\bigg(\frac{(\sqrt{M}r_c\phi+Mu+r_c)}{M},
  -r,-\frac{1}{\sqrt{M}}\bigg),\nonumber \\
  l_0 &=& (-\frac{r_c}{M},0,\frac{1}{\sqrt{M}}),\nonumber \\
  l_{-1} &=& e^{\sqrt{M}\phi}\bigg(\frac{(\sqrt{M}r_c\phi+Mu-r_c)}{M},
  -r,\frac{1}{\sqrt{M}}\bigg),\nonumber  \\
  m_1 &=& -e^{-\sqrt{M}\phi}\big(\frac{1}{\sqrt{M}},\sqrt{M},0\big),\nonumber \\
  m_0 &=& (-\frac{1}{\sqrt{M}},0,0),\nonumber \\
  m_{-1} &=& e^{\sqrt{M}\phi}\big(-\frac{1}{\sqrt{M}},-\sqrt{M},0\big).
\end{eqnarray}

To find the explicit form of $Y(\phi)$, we expand the first line of (\ref{zeta-therm}) using the Killing vectors (\ref{Kil-FSC-2}), then we obtain
\begin{equation}\label{Y-2}
  Y=b_1(\frac{e^{-\sqrt{M}\phi}}{\sqrt{M}})+b_0(\frac{1}{\sqrt{M}})+b_{-1}(\frac{e^{\sqrt{M}\phi}}{\sqrt{M}})
\end{equation}

Fixing the center point of the interval at $\phi=0$, we set $b_1=b_{-1}$. $b_0$ and $b_1$ could be considered as follows \cite{Jiang}
\begin{eqnarray}
  b_0 &=& \frac{\pi}{\tilde{\beta}_{\phi}}\bigg(\tanh\big(\frac{\pi l_{\phi}}{2\beta_{\phi}}\big)+\coth\big(\frac{\pi l_{\phi}}{2\beta_{\phi}}\big)\bigg)\nonumber \\
  b_1 &=& -\frac{\pi}{\tilde{\beta}_{\phi}}\text{csch}\big(\frac{\pi l_{\phi}}{\beta_{\phi}}\big) \label{bb}
\end{eqnarray}

Substituting the definitions (\ref{bb}) into (\ref{Y-2}) and using some manipulations, we find the following form for $Y(\phi)$
\begin{equation}\label{Y-3}
  Y(\phi)=\frac{2\pi}{\sqrt{M}\sinh\big(\frac{\sqrt{M}l_{\phi}}{2}\big) }
  \big[\cosh\big(\frac{\sqrt{M}l_{\phi}}{2}\big)-\cosh(\sqrt{M}\Delta \phi)\big]
\end{equation}
where we use the following identifications
\begin{equation}\label{ids-2}
  \Delta\phi\equiv \phi-(\phi_++\phi_-)/2.
\end{equation}

A similar calculation could be done to find $T(\phi)$. To this end, we consider $d_1=d_{-1}$ and we have the following definitions \cite{Jiang}
\begin{eqnarray}
  d_1 &=& -\frac{\pi~\text{csch}(\pi l_{\phi}/\beta_{\phi})\big[-\beta^2_{\phi}\tilde{\beta}_u
  +\pi(l_{\phi}\beta_u+l_u\beta_{\phi})\tilde{\beta}_{\phi}
  \coth(\pi l_{\phi}/\beta_{\phi})\big]}{\beta^2_{\phi}\tilde{\beta}^2_{\phi}}, \nonumber \\
  d_0 &=& \frac{\pi~\text{csch}^2(\pi l_{\phi}/\beta_{\phi})\big[-\beta^2_{\phi}\tilde{\beta}_u \sinh(2\pi l_{\phi}/\beta_{\phi})
  +2\pi(l_{\phi}\beta_u+l_u\beta_{\phi})\tilde{\beta}_{\phi}
  \big]}{\beta^2_{\phi}\tilde{\beta}^2_{\phi}} \label{dd}
\end{eqnarray}
We find $T(\phi)$ as follows
\begin{eqnarray}
  T(\phi) &=& \frac{\pi}{2M\sinh\big(\frac{\sqrt{M}l_{}\phi}{2}\big)}\bigg\{(Jl_{\phi}+2l_uM)
  \big[\coth\big(\frac{\sqrt{M}l_{\phi}}{2}\big)\cosh(\sqrt{M}\Delta\phi)-\text{csch}\big(\frac{\sqrt{M}l_{\phi}}{2}\big)
  \big]
  \nonumber \\
   &+& \frac{2J}{\sqrt{M}}\big[\cosh(\sqrt{M}\Delta\phi)-\cosh\big(\frac{\sqrt{M}l_{\phi}}{2}\big)
   \big]-2J\Delta\phi\sinh(\sqrt{M}\Delta\phi)\bigg\}
\end{eqnarray}

A similar computation to that we did to find the modular Hamiltonian (\ref{H in FT}), can be done to find the following result for the modular Hamiltonian for a single interval on a thermal BMSFT
\begin{equation}\label{H in FT-therm}
  \mathcal{H}_{\zeta} = -\frac{1}{2\pi}\int_{\phi_-}^{\phi_+}\big(T(\phi)\mathcal{P}(\phi)+Y(\phi)\mathcal{J}(\phi)\big)
\end{equation}

There is the correspondence (\ref{H=h}) between the vacuum symmetry generators $h_i$ at the boundary and the Killing vectors $H_i$ in the bulk. The Killing vectors in the bulk are given in terms of the $T(\phi)$ and $Y(\phi)$ by (\ref{Killing GMMG-2}). This form of the Killing vectors that we have derived can be used in the FSC geometry. The gravitational charge associated with the Killing vectors is then given by
\begin{eqnarray} \label{F-charge GMMG2}
  Q &=& M(T)+J(X)+L(Y)+P(Z)\nonumber \\
   &=&  -\frac{1}{16\pi G}(\sigma+\frac{\alpha H}{\mu}+\frac{F}{\mathrm{m}^2})\int_{0}^{2\pi}\bigg\{T(\phi)\mathcal{M}(\phi)+2Y(\phi)\mathcal{L}(\phi)\bigg\}d\phi\nonumber \\
   &+& \frac{1}{16\pi G \mu}\int_{0}^{2\pi}Y(\phi)\mathcal{M}(\phi)d\phi
\end{eqnarray}
where $Y(\phi)$, and $T(\phi)$ are given by (3.58), (3.61) respectively.
Using the holographic disctionary (\ref{var.rel}), we find that the gravitational charge (\ref{F-charge GMMG2}) evaluated at the asymptotic boundary agrees with the modular Hamiltonian (\ref{H in FT-therm}) in the field theory side. This is another evidence for the $flat_3/BMSFT$ correspondence. As shown in \cite{Apolo-2}, the gravitational charge along the swing surface equals the entanglement entropy. Thus, we can generalize the first law of entanglement entropy (\ref{match1}) to the thermal BMSFTs.

\section{Conclusion}\label{sec:4}
In this paper we studied a new class of non-AdS holography. This is the $Flat_3/BMSFT$ holography in the framework of GMMG model \cite{Setare-Adami Nucl 2018}. By applying the Bondi fall-off conditions \cite{Bondi-1,Sach-1,Sach-2}, we found the Killing vectors (\ref{Killing GMMG}) of the isometries that preserve the asymptotic form of the metric in (\ref{B.C.}). The conserved charges corresponding to these vectors are dual to the modular Hamiltonian of the field theory at the boundary.

To find the modular Hamiltonian \cite{Haag-bk,Casini-1}, we considered a BMS-invariant field theory at the boundary. The BMSFTs are considered non-thermal on the plane and thermal on the cylinder. To compute the vacuum modular Hamiltonian, the generalized Rindler transformation was used. In this method, the vacuum state was defined on an interval $\mathcal{I}$ that was mapped to a Rindler spacetime under a symmetry transformation. The modular flow generator (\ref{zeta-1}) was considered as a linear combination of the vacuum symmetry generators. Using the conditions on the Rindler transformations, the coefficients (\ref{aj2}) and (\ref{bj2}) of the modular flow generator were computed.

To compute the modular Hamiltonian,  the vacuum symmetry generators $l_j$ and $m_j$ were replaced with the corresponding conserved charges $\mathcal{L}_j$ and $\mathcal{M}_j$. By this replacement the modular flow generator (\ref{zeta-1}) was replaced with the modular Hamiltonian (\ref{H in FT}). Using the equal forms of the charge algebras (\ref{BMS algebra}) of the two sides of the holography, we found the dictionary (\ref{var.rel}) between the variables of the bulk and the boundary. Using this dictionary of the holography, we showed that the gravitational charge (\ref{F-charge GMMG}) is the same as the
modular Hamiltonian (\ref{H in FT}). Finally we were able to show the first law of entanglement in (\ref{match1}).

We first did these steps for a Poincare patch (\ref{Poin}) of GMMG model and a zero temperature BMSFT on a plane (\ref{interval}) on two sides of the holography. Then, we considered a general case in which the metric has a flat space cosmological form (\ref{FSC}). In this case, the holographic dual of the gravity is a finite temperature BMSFT on a cylinder (\ref{therm. cyl}). Using the asymptotic form of the Killing vectors of FSC (\ref{Kil-FSC-2}), the modular Hamiltonian for a single interval on a thermal BMSFT was computed in (\ref{H in FT-therm}). The computations associated with the conserved charges of GMMG model were repeated in (\ref{F-charge GMMG2}). In this general case, the first law of holographic entropy was obtained too.
It is interesting that our results for GMMG model reduces to the result for Einstein gravity in $2+1$ dimension, when $\mu , m^2 \to \infty $, and $\sigma=-1$. For example one can explicitly see that our expression for variation of conserved gravitational charge in Eq(3.17) reduce to the
 expression (3.51)  in the paper \cite{Apolo} for Einstein gravity in mentioned limiting cases. These limiting cases, are a place in parameter space of model, where GMMG model reduces to the Einstein gravity.


\section{Introduction}

The early universe was filled with a hot dense plasma composed primarily of electrons, protons, and helium nuclei in a bath of radiation tightly coupled to the plasma due to the high rate of Thomson scattering. As the universe cooled due to cosmic expansion, the plasma recombined to form a transparent gas of neutral hydrogen and helium. Shortly after recombination, the photons decoupled and have traveled mostly unimpeded ever since. These photons, which continued to redshift due to the cosmic expansion, make up the cosmic microwave background (CMB).

Observations of the CMB have provided a wealth of information about the history, contents, and evolution of the universe. The existence of a thermal bath of radiation provides very strong evidence for the Hot Big Bang model of cosmology, and measurements of the CMB temperature and polarization anisotropies have allowed for percent-level determination of the parameters that define the current $\Lambda$CDM model of cosmology~\cite{Planck2018cosmo}. Future CMB surveys will extract a great deal more information from the temperature and polarization anisotropies.

Looking beyond primary CMB anisotropies, additional processes after recombination leave imprints on the CMB. One such process is Rayleigh scattering of CMB photons by neutral atoms, especially by hydrogen, which achieves maximum density shortly after recombination. Rayleigh scattering is a classical scattering process of long wavelength radiation by the induced dipole of neutral species. Its impact on recombination and CMB anisotropies has been studied in \citet{Lewis2013,Yu2001, Alipour2014}. The spectral function for Rayleigh is well described by a cross-section that scales as $\nu^4$. However, as we will show, observing Rayleigh scattering of the CMB remains challenging due to the presence of astrophysical foregrounds that act as a source of confusion at the relevant frequencies.

Physically, Rayleigh scattering causes a frequency-dependent damping of CMB anisotropies by scattering photons from hot spots out of the line-of-sight and from cold spots into the line of sight. At high frequencies the effect of Rayleigh scattering contributes to the angular power spectrum, and thus it must be consistently included in cosmological modeling.

Observing the effects of Rayleigh scattering of the CMB would provide additional information that can be used to probe the early universe, including the physics of recombination,  the ionization history, the expansion rate of the early universe, and the properties of primordial perturbations. One interesting application of CMB Rayleigh scattering is to improve measurements of the light relic density ($N_\mathrm{eff}$), which further probes physics beyond the Standard Model in the dark sector \citep{Green2019,Dvorkin:2022jyg}. With primary CMB observations alone, measuring $N_\mathrm{eff}$ requires significantly more resources when compared to including Rayleigh scattering information \citep[e.g.,][]{CMBS4,Beringue2020}. However, the value of Rayleigh in constraining physics beyond the Standard Model requires a significant detection of the signal. For example, for experiments such as CCAT-prime~\citep{CCATp}, the Simons Observatory (SO)~\citep{SOforecasts}, and the \textit{Planck} satellite~\citep{Planck2018}, even under ideal circumstance, where foregrounds and other possible systematic effects are neglected, the expected statistical significance of the Rayleigh scattering signal would rise to the level of a detection, but would not significantly improve constraints on cosmological parameters~\citep{Beringue2020}.

The amplitude of the Rayleigh auto-spectrum is very small, both in temperature and in polarization. At around $\nu=500~\mathrm{GHz}$, the temperature Rayleigh scattering auto-spectrum is still around $7$ orders of magnitude lower than the primary CMB spectrum~\citep[][]{Lewis2013}. Hence a detection of the auto-spectrum will remain out of reach for the foreseeable future. On the other hand, cross-correlations between the Rayleigh scattering signal and the primary temperature or polarization anisotropies are more promising~\citep{Lewis2013,Alipour2014,Beringue2020}. The amplitude of the cross signals is around $3$ orders of magnitude larger than the auto-spectrum of Rayleigh scattering, as shown in Fig.~\ref{contribution}. Raw noise levels from current ground-based CMB experiments, the \textit{Planck} satellite and future CMB satellite missions~\citep[e.g., PICO][]{PICO} suggests the Rayleigh scattering cross-correlations should be detectable~\citep[e.g.,][]{Lewis2013,Alipour2014,Beringue2020}. The main challenge is to isolate the Rayleigh scattering signal in the presence of foregrounds and the atmosphere, both of which can have a large impact on detectability. While the role of the atmosphere is only relevant for ground-based experiments, astrophysical foregrounds will affect both satellite missions and ground-based experiments. This paper aims to address the effect of foregrounds on detectability of the Rayleigh scattering signal. Indeed, the frequency dependence of some foreground emissions make their amplitude significant in the frequency range where the Rayleigh scattering is brighter. It is timely to assess how large this impact will be for future CMB experiments.

Previously, \citet{Alipour2014} studied the effect of foregrounds and foregrounds removal on the Rayleigh scattering signal and concluded that the $TE$ cross-power spectra between primary temperature and Rayleigh scattering polarization signals should be relatively unbiased by foregrounds. More recently, \citet{Beringue2020} forecasted a detection of the $TT$ cross-correlation with $4.7\sigma$ for \textit{Planck}, and a whopping $715\sigma$ detection with PICO, both in the absence of foregrounds~\citep{Lewis2013,Beringue2020}. Throughout this paper, we denote the signals as $XY=\mathrm{CMB}\;(T/E)\times \mathrm{RS}\;(T/E)$, where $T$ represents the temperature fluctuations, $E$ refers to the $E$-mode polarization, and RS denotes the Rayleigh scattering component.
These calculations show that $TT$ cross-correlation signal-to-noise ratio is roughly an order of magnitude larger than the $TE$ or $ET$  signal-to-noise. Therefore, given the large difference in signal-to-noise in the absence of foregrounds, a more comprehensive study of all Rayleigh scattering cross-correlations, including foreground removal is warranted.

We use a component separation method to quantify the impact of foreground removal on the Rayleigh scattering signal-to-noise. Several component separation methods and algorithms already exist in the literature, including {\it blind} methods like Spectal Matching Independent Component Analysis (SMICA)~\citep[][]{2003MNRAS.346.1089D}, which was used in \textit{Planck} data analysis \citet{2014A&A...571A..12P}, Fast Independent Component Analysis (FastICA)~\citep[][]{2002MNRAS.334...53M}, the Joint Approximate Diagonalization of Eigenmatrices
method (JADE)~\citep[][]{2006A&A...455..741P}, and Correlated Component Analysis (CCA)~\citep[][]{2005EJASP2005.2400B}.
We chose to be conservative in our approach to component separation and selected the method known as Internal Linear Combination (ILC)~\citep[][]{Delabrouille2009}, which does not assume particular parametrization of the foreground components or templates.

Since we are dealing with the cross-spectra of two signals, we use a variation of this method, the constrained Internal Linear Combination (cILC)~\citep{Remazeilles_2010}. With this method, we avoid residual bias from the CMB primary on the Rayleigh signal, and vice versa in the cross-spectra \citep[][]{RemChluba2018}. Given the conservative nature of our approach, these forecasts should be considered as lower limits on the possible signal-to-noise.

In this paper, we introduce the physical modeling and scientific significance of the Rayleigh scattering signal in Section~\ref{sec:RSPhys}. We set up our modeling of the foregrounds and describe the experimental configurations we consider in Section~\ref{sec:fgmodel} and Section~\ref{sec:Tele}; additional information about the surveys can be found in Appendix~\ref{sec:Telescopes}. The CMB and Rayleigh spectra are generated used a modified version of \texttt{CAMB}\footnote{\url{https://camb.info/}}~\citep{Lewis:1999bs}. Details on the cILC method are presented in Section ~\ref{sec:ILC}. In Section~\ref{sec:Results} we forecast the signal-to-noise for CCAT-prime, SO and \textit{Planck}. In addition, we explore the effect of systematics: bias and gain calibration. We show how foregrounds impact our observation. Then, we perform the same analysis for PICO and explore combinations of data from LiteBIRD with CCAT-prime and \textit{Planck}. We use the extrapolated atmosphere noise model on all frequencies from ACT~\citep[][]{D_nner_2012}. The atmospheric noise for CCAT-prime and SO might differ from these estimates, so our results should be considered with this in mind. We sometimes use $\mathcal{D}_{\ell}$ which is related to ${C}_{\ell}$ as:
\begin{equation}
  \mathcal{D}_{\ell}=\frac{\ell(\ell+1)}{2\pi}C_{\ell}\, .
\end{equation}

\section{Rayleigh scattering Physics}
\label{sec:RSPhys}
Rayleigh scattering describes the classical scattering of electromagnetic waves by particles with size is much smaller than the wavelength of the radiation. To first order, the cross section of this scattering process is proportional to $\nu^4$, where $\nu$ is the photon frequency. Its impact on the recombination history and decoupling of CMB photons has been described in details in Refs.~\citep{Yu2001,Alipour2014,Lewis2013}. We will summarise the cosmological implications of Rayleigh scattering in this section.

In the early universe, photons are tightly coupled with the free electrons in the primordial plasma by Thomson scattering. As the universe expands, the temperature decreases until recombination of free electrons and protons becomes thermally favored. As the fraction of free electrons drops, the Thomson scattering rate falls. Around redshift $z\sim 1100$, photons experience a last scattering event and start free-streaming.

However, the recombination processes also yield the production of neutral species, mainly hydrogen and helium, by which CMB photons can be scattered through Rayleigh scattering. Unlike Thomson scattering, the Rayleigh scattering cross section is frequency dependent:
\begin{equation}
  \sigma_R(\nu)=\sigma_T \sum_{n=1} a_n\left( \frac{\nu}{\nu_0}\right)^{2n+2},
\end{equation}
where $a_1 = 1$, $a_2 = 2.626$, $a_3 = 5.502$ \citep{Lewis2013}, $\sigma_T$ is the cross section of Thomson scattering and $\nu_0=3.102\times 10^6~\mathrm{GHz}$ is the Lyman limit frequency.
After recombination, neutral species quickly dilute due to the cosmic expansion (scaling as $\propto a^{-3}$) and photons are diluted and redshifted ($\propto a^{-4}$), making the probability of Rayleigh scattering $\propto a^{-7}$. As a result, most Rayleigh scattering events occur shortly after recombination. This can be seen in the visibility function $g(z) \equiv \dot{\tau}e^{-\tau}$ where $\tau$ is the optical depth. The visibility function describes the probability that a photon last scattered at a given redshift. The visibility functions are shown in Fig.~\ref{visibility} for both Thomson and Rayleigh scattering.

\begin{figure}[t!]

  \centering

  \includegraphics[width=\linewidth]{other/visibility.pdf}

\caption{Change to the visibility function $g(\tau)=\Dot{\tau}e^{-\tau}$ due to Rayleigh scattering at $\nu=280~\mathrm{GHz}$ and $\nu=500~\mathrm{GHz}$. The dashed lines represent the normalized contribution from Rayleigh scattering alone. The Rayleigh scattering signal shifts the peak of the visibility function towards lower redshifts.}
\label{visibility}
\end{figure}

By increasing the comoving opacity of the plasma in a frequency-dependent way, Rayleigh scattering leaves characteristic imprints on the CMB:
\begin{itemize}
\item {\bf Frequency-dependent shift of the visibility function:}
    As can be seen in Fig.~\ref{visibility}, the visibility function is shifted towards lower redshifts (later time). Crucially, the position of the peak of the visibility function is made frequency dependent with the shift towards later time being larger at higher frequencies.
\item {\bf Increase of diffusion damping:}
    Rayleigh scattering globally increases the comoving opacity of the plasma. The mean free path of photons in the plasma is therefore reduced. The strength of diffusion damping is directly controlled by the mean free path~\citep{visibility}. At lower frequencies, Rayleigh scattering has negligible impact. However, at higher frequencies, the visibility function is shifted towards later time when cosmic expansion has led to a global decrease of the comoving opacity, leading to an overall increase of diffusion damping.
\item {\bf Frequency-dependent sound horizon:}
    As the peak of the visibility function is rendered frequency dependent by Rayleigh scattering, the size of the sound horizon at last scattering similarly inherits the same frequency dependence, being larger at higher frequencies. The size of the sound horizon dictates the location of the acoustic peaks in both the matter and the CMB power spectra, a larger sound horizon shifts the acoustic peaks towards larger scales.
\item {\bf Frequency-dependent amplitude of polarization signal:}
    The shift of the visibility function induced by Rayleigh scattering increases the amplitude of the local temperature quadrupole around the time of recombination, leading to a boost in $E$ modes on large scales.
\end{itemize}

If the shift of the visibility function is small enough, the effect of Rayleigh scattering on the temperature and polarization power spectra is accurately described by a linear, frequency-dependent distortion \citep{Lewis2013}:
\begin{equation}
    a_{\ell m}^{X}(\nu)=a_{\ell m}^{X}+\Big(\frac{\nu}{\nu_0}\Big)^4\Delta a_{\ell m}^{X,4}+ \cdots,
    \label{eq:rs-sht}
\end{equation}
where $a_{\ell m}$ are the coefficients in the spherical harmonic transform of the temperature ($X = T$) and polarization ($X = E,B$) anisotropies. Note that we have neglected higher order corrections to the Rayleigh scattering cross section that would give rise to additional terms on the right hand side of Eq.~\eqref{eq:rs-sht}. The power spectra are then given by
\begin{equation}
  \begin{aligned}
  C_{\ell}^{XY}(\nu_1,\nu_2)&=\left\langle a_{\ell m}^{{X}}(\nu_1)a_{\ell m}^{{Y}}(\nu_2) \right\rangle \\ &=C_{\ell}^{XY}+ \left(\frac{1}{\nu_0}\right)^4 \left(\nu_1^4C_{\ell}^{X\Delta Y_4}+\nu_2^4C_{\ell}^{\Delta X_4Y}\right)\\
  &\quad+\left(\frac{\nu_1\nu_2}{\nu_0^2}\right)^4 C_{\ell}^{\Delta X_4\Delta Y_4}+\cdots .
  \end{aligned}
\end{equation}
The first term in this expansion is the primary CMB, sourced by Thomson scattering only. The following term contains the cross-correlation spectra between the primary CMB and the Rayleigh scattering signal. The last term is proportional to the Rayleigh scattering auto-spectrum. These terms are also ordered with decreasing amplitude (see Fig.~\ref{contribution}) making the cross-correlation the natural target to achieve a first detection.

\begin{figure*}[!ht]
\centering
 \subfloat{\includegraphics[width=0.49\textwidth]{other/280contribution.pdf}}
 \hfill
 \subfloat{\includegraphics[width=0.49\textwidth]{other/500contribution.pdf}}
\caption{Contributions from Rayleigh scattering at $\nu=280~\mathrm{GHz}$ and $\nu=500~\mathrm{GHz}$. Red lines show the cross-spectra between the Rayleigh signal and the primary CMB. The solid lines present positive and the dotted parts are negative correlations. Yellow lines show the auto-spectra of the Rayleigh signal. The cross-spectra are roughly $3$ orders of magnitude larger than the auto-spectra of the Rayleigh signal.
}

\label{contribution}
\end{figure*}

With the low-noise and large-area maps expected from future CMB experiments, the detection of this faint Rayleigh scattering signal should be possible~\citep{Lewis2013,Beringue2020}. More realistic estimates of the expected detection significance require forecasts including foregrounds and testing the robustness of these predictions after applying foreground cleaning techniques.

\section{Foreground Modeling}
\label{sec:fgmodel}
The microwave sky foreground models we considered in this paper include the cosmic infrared background (CIB), the thermal Sunyaev-Zeldovich effect (tSZ), the kinetic Sunyaev-Zeldovich (kSZ), radio sources, and both galactic dust and synchrotron emissions. The CIB is sourced by star-forming dusty galaxies and is modelled according to \citet{Dunkley2013}. The tSZ and kSZ are the two dominant effects in the scattering of CMB photons by free electrons in clusters. We model the tSZ and kSZ auto-spectra using templates from \citet{BBPSS}. We also include the cross-spectra between tSZ and CIB in our foreground modeling, following \citet{addison2012}. In principle there is also a polarized version of the SZ effect (pSZ) and the relativistic SZ effect (rSZ), but these are too small to have an impact on our forecasts, so we do not consider them in this analysis. Radio point sources (radio galaxies) are modeled following \citet{SOforecasts}. The angular power spectrum in temperature and polarization for all of these foregrounds, the CMB, and Rayleigh scattering are shown in Figure~\ref{fg}.

\begin{figure*}[!ht]
\centering
 \subfloat{\includegraphics[width=0.48\textwidth]{noisefig/comp.pdf}}
 \hfill
 \subfloat{\includegraphics[width=0.48\textwidth]{noisefig/comppol.pdf}}
\caption{Signal and foreground power spectra at $\nu=280~\mathrm{GHz}$. The observed temperature includes contributions from the CMB primary signal, Rayleigh scattering auto signal, radio point sources, thermal Sunyaev–Zeldovich (tSZ), kinetic Sunyaev–Zeldovich (kSZ), cosmic infrared background (CIB), and cross correlation of tSZ and CIB; for polarization, the polarized primary CMB and polarized Rayleigh signal, dust, synchrotron, and radio point sources are included.}
\label{fg}
\end{figure*}

The polarized dust and synchrotron foregrounds from our galaxies are modeled following \citet{DUSTSYNC}. For the dust model the emissivity is assumed to follow:
\begin{equation}
  \epsilon_\mathrm{d}(\nu,\beta_\mathrm{d})=\left(\frac{\nu}{\nu_0}\right)^{\beta_\mathrm{d}-2}\frac{B_{\nu}(T_\mathrm{d})}{B_{\nu_0}(T_\mathrm{d})}.
\end{equation}
Here $\beta_\mathrm{d}=1.59$, $T_\mathrm{d}=19.6~\mathrm{K}$, $\nu_0=353~\mathrm{GHz}$, and $B_{\nu}$ is the Planck function.

We can then express the dust and synchrotron power spectrum as:
\begin{align}
  \mathcal{D}^{\mathrm{d},\mathrm{s}}_{\ell}(\nu_1,\nu_2)   =& \left[\alpha_\mathrm{s}^2\left(\frac{\nu_1\nu_2}{\nu_K^2}\right)^{\beta_\mathrm{s}}+\alpha_\mathrm{d}^2\epsilon_\mathrm{d}(\nu_1,\beta_\mathrm{d})\epsilon_\mathrm{d}(\nu_2,\beta_\mathrm{d})\right] \nonumber \\
  &\times
  g(\nu_1)g(\nu_2)\left(\frac{\ell}{80}\right)^{-0.42} \, ,
\end{align}
where $\alpha_\mathrm{s}$ and $\alpha_\mathrm{d}$ are two scaling constants, which we approximate as $\alpha_\mathrm{s}^2=10.3~\mu \mathrm{K}^2$ and $\alpha_\mathrm{d}^2=0.53~\mu \mathrm{K}^2$ in antenna temperature, $\beta_\mathrm{s}=-3.1$ is the synchrotron spectral index. These parameters are measured from the cleanest 45\% of the sky, on angular scales with $\ell=50\mhyphen 110$ \citep[][]{DUSTSYNC}. We ignore the cross-spectra between dust and synchrotron since it is small \citep[][]{DUSTSYNC}. Finally, $g(\nu)$ is the conversion function from antenna to CMB thermodynamic temperature at frequency $\nu$:
\begin{equation}
  g(\nu)=\frac{\left(e^{\frac{h\nu}{k_\mathrm{B} T_\mathrm{CMB}}}-1\right)^2}{e^{\frac{h\nu}{k_\mathrm{B} T_\mathrm{CMB}}}\left(\frac{h\nu}{k_\mathrm{B} T_\mathrm{CMB}}\right)^2},
\end{equation}
where $h$ is the Planck constant and $k_\mathrm{B}$ is the Boltzmann constant.

For radio point sources, we write the auto-spectrum as:
\begin{equation}
  \mathcal{D}_{\ell}^\mathrm{ps}(\nu_1,\nu_2)=A_\mathrm{ps}\left(\frac{\nu_1\nu_2}{\nu_r^2}\right)^{-2.5}\frac{g(\nu_1)g(\nu_2)}{(g(\nu_r))^2}\frac{\ell^2}{\ell_0^2},
  \label{eq:radio-ps}
\end{equation}
where $A_\mathrm{ps}=3.1~\mu \mathrm{K}^2$, $\ell_0=3000$, $\nu_r=150~\mathrm{GHz}$.

For polarized radio point sources, we modify Eq.~\eqref{eq:radio-ps} as:
\begin{equation}
  \mathcal{D}_{\ell}^\mathrm{ps-pol}(\nu_1,\nu_2)=A_\mathrm{ps-pol}\left(\frac{\nu_1\nu_2}{\nu_r^2}\right)^{-2.5}\frac{g(\nu_1)g(\nu_2)}{(g(\nu_r))^2}\frac{\ell^2}{\ell_0^2}
\end{equation}
where $A_\mathrm{ps-pol}=0.015~\mu\mathrm{K}^2$ at $\ell_0=3000$ and $\nu_r=150~\mathrm{GHz}$\citep{SOforecasts}.
We assume that the galactic dust and synchrotron do not depend on the position of the sky observed. We show all auto-spectra at $\nu=280~\mathrm{GHz}$ in Fig.~\ref{fg} for temperature and polarization. We apply the best-fit value from \citet{Dunkley2013} to the extragalactic foreground models.

\section{Instrumental and Atmospheric Noise Models}
\label{sec:Tele}
We consistently model the noise power spectra of the surveys we considered, CCAT-prime, \textit{Planck}, Simons Observatory (SO), PICO and LiteBIRD using a model that consists in a {\it red} and a {\it white} component:
\begin{equation}
  n^\mathrm{ins}_{\ell}=n_\mathrm{red}\left(\frac{\ell}{\ell_\mathrm{knee}}\right)^{\alpha_\mathrm{knee}}+n_\mathrm{white} \, .
\end{equation}
Here $n_\mathrm{white}$ is the white noise due to the instrument alone, $n_\mathrm{red}$ is the red noise due to the atmosphere, where $\ell_\mathrm{knee}$ and $\alpha_\mathrm{knee}$ parameterize the scales affected by atmospheric contamination \citep{SOforecasts}. For satellite experiments such as {\it Planck}, we simply set $n_\mathrm{red}=0$. For ground-based observations, the red noise component is estimated from empirical measurements of the large angular scale noise from ACT at $150 ~{\rm GHz}$, and extrapolated to higher frequencies for SO, CCAT-prime, and CMB-S4 \citep{SOforecasts,CCATp,CMBS4}. Our forecasts will include \textit{Planck} to complement future ground-based experiments on the largest angular scales~\citep{Planck2018}.

The Simons Observatory consists in a set of telescopes under construction in Chile. Its large aperture telescope (LAT) will have a small beam allowing for precision measurement of small angular scales.  We refer to the SO LAT as `SO' throughout \citep{SOforecasts}. The Simons Observatory experiment also includes three small aperture telescopes (SATs), with larger beams. Because of the large beams, these SATs will be limited to $\ell<400$. With polarization modulators, SATs will only produce polarization maps; hence we do not consider temperature anisotropies measurements from the SATs.

The CCAT-prime observatory, currently under construction in Chile, will consist in the Fred Young Sub-millimeter Telescope (FYST) equipped with two instruments: an imaging camera, Prime-CAM \citep{Choi_2020} and a multi-beam submillimeter heterodyne spectrometer, CHAI \citep{CHAI}. We will refer to the FYST + Prime-CAM configuration as `CCAT-prime' for simplicity. Similar to SO, CCAT-prime will have a high sensitivity at small angular scales \citep{Choi_2020}. Compared to the SO LAT, CCAT-prime will focus on higher frequencies. Besides, the higher altitude of the site is expected to yield lower atmospheric noise.

LiteBIRD is a planned satellite mission covering a large frequency range. Due to the relatively small dish the beam will limit its resolution to $\ell < 400$~\citep{LiteBird}.

The Probe of Inflation and Cosmic Origin (PICO) is a recently proposed space satellite \cite{PICO}, with extremely low noise and many frequency channels. PICO will cover a large range of multipoles. If funded, it will be capable of a high significance measurement of Rayleigh scattering \citep{Beringue2020}.

All the noise parameters for the experiments mentioned above are listed in Appendix~\ref{sec:Telescopes}. For the ground-based experiments, we expect some differences between the measured red noise components and extrapolated values we use. Additionally we do not consider that the atmospheric red noise can be removed using proposed analysis techniques \citep{SPTinprep}. Thus, our forecast in this respect are conservative, since if such techniques are viable this would only improve the forecasts we present.

\begin{figure*}[!ht]
\centering
 \subfloat{\includegraphics[width=0.483\linewidth]{noisefig/ins.pdf}}
 \subfloat{\includegraphics[width=0.49\linewidth]{noisefig/pol.pdf}}
\caption{Comparison of instrumental noises of \textit{Planck}, CCAT-prime, SO LAT, and SO SAT at similar frequency channels. CCAT-prime and SO have smaller noise on small scales, while \textit{Planck} and SAT have smaller noise on large scales. In practice, we only use data from SO SAT for $\ell\leq400$.}

\label{fig:instrumental}
\end{figure*}

\section{Component Separation and Reconstructed Noise Spectra}
\label{sec:ILC}

The forecasts carried out in this work follow the constrained Internal Linear Combination (cILC) framework derived in \citet{Remazeilles_2010}\footnote{Technically, our analysis is not an cILC, since we are not constructing the covariance matrix from a map or data, but we will refer to it as such.}. Since we are interested in measuring the cross correlation between the primary CMB and the Rayleigh scattering signals, we schematically need to produce two maps: one that contains only the CMB signal and no Rayleigh scattering and another one that contains only the Rayleigh scattering signal and no primary CMB. This procedure will ensure that no residual CMB contamination in the Rayleigh scattering map (and vice-versa) would bias the cross correlation. We model the total data vector at a pixel $p$ as:
\begin{equation}
\mathbf{d}(p) = \mathbf{a}^{\rm CMB} s^{\rm CMB}(p) + \mathbf{b}^{\rm RS} s^{\rm RS}(p) + \mathbf{n}(p),
\label{eq:signal}
\end{equation}
where $\mathbf{a}$ is the frequency dependence of the CMB (assumed to be unity at all frequencies), $\mathbf{b}$ is the frequency dependence of the Rayleigh scattering signal ($\propto \nu^4$). $\mathbf{n}$ captures the effect of foregrounds, instrumental and atmospheric noises. A pixel here refers to either an actual pixel in a sky map or a set of harmonic coefficients $(\ell,m)$ of the spherical harmonic transform.

The cILC procedure yields a set of weights $\mathbf{w}$ that, for the CMB map, satisfy:
\begin{itemize}
    \item $\mathbf{w}^\mathsf{T}{}\cdot\mathbf{a} = 1$.
    \item $\mathbf{w}^\mathsf{T}{}\cdot\mathbf{b} = 0$.
    \item $\hat{s}^{\rm CMB}(p) \equiv \mathbf{w}^\mathsf{T}{}\cdot\mathbf{d}(p)$ has minimum variance.
\end{itemize}
These weights are given by \citep{Remazeilles_2010}:
\begin{equation}
\mathbf{w_{a=1,b=0}}^\mathsf{T}{} = \frac{ \left(\mathbf{b}^\mathsf{T}{} \mathbf{\Tilde{R}_d}^{-1}   \mathbf{b}\right)\mathbf{a}^\mathsf{T}{}  \mathbf{\Tilde{R}_d}^{-1} - \left(\mathbf{a}^\mathsf{T}{} \mathbf{\Tilde{R}_d}^{-1}   \mathbf{b}\right)\mathbf{b}^\mathsf{T}{}  \mathbf{\Tilde{R}_d}^{-1} }{ \left(\mathbf{a}^\mathsf{T}{} \mathbf{\Tilde{R}_d}^{-1}   \mathbf{a}\right) \left(\mathbf{b}^\mathsf{T}{} \mathbf{\Tilde{R}_d}^{-1}   \mathbf{b}\right) - \left(\mathbf{a}^\mathsf{T}{} \mathbf{\Tilde{R}_d}^{-1}   \mathbf{b}\right)^2}.
\label{eq:cilc-weights}
\end{equation}
Here $\mathbf{\Tilde{R}_d}$ is the covariance of the input maps. If the cILC is carried out in the harmonic domain, the weights will inherit a scale dependence: $\mathbf{w} \longrightarrow \mathbf{w}_\ell$. Finally, we note that the same procedure yields the weights for the Rayleigh scattering map ($\mathbf{w_{a=0,b=1}}$) by simply inverting $\mathbf{a}$ and $\mathbf{b}$ in Eq.~\eqref{eq:cilc-weights}.

From the linear combinations $\hat{s}^{\rm CMB}(p) = \mathbf{w_{a=1,b=0}}^\mathsf{T}{} \cdot \mathbf{d}(p)$ and $\hat{s}^{\rm RS}(p) = \mathbf{w_{a=0,b=1}}^\mathsf{T}{} \cdot \mathbf{d}(p)$ we can compute the primary CMB and Rayleigh scattering auto and cross-spectra, that, given the constraints obeyed by the weights, read:
\begin{align}
    \hat{C}_\ell^{\rm CMB} &= C_\ell^{\rm CMB} + \mathbf{w_{a=1,b=0}}^\mathsf{T}{}  \mathbf{N}_\ell   \mathbf{w_{a=1,b=0}}, \nonumber \\
    &= C_\ell^{\rm CMB} + N_\ell^{\rm CMB}, \nonumber \\
    \hat{C}_\ell^{\rm RS} &= C_\ell^{\rm RS} + \mathbf{w_{a=0,b=1}}^\mathsf{T}{}  \mathbf{N}_\ell   \mathbf{w_{a=0,b=1}}, \nonumber \\
    &= C_\ell^{\rm RS} + N_\ell^{\rm RS}, \nonumber \\
    \hat{C}_\ell^{\rm CMB \times RS} &= C_\ell^{\rm CMB\times RS} + \mathbf{w_{a=0,b=1}}^\mathsf{T}{}  \mathbf{N}_\ell   \mathbf{w_{a=1,b=0}}, \nonumber \\
    &= C_\ell^{\rm CMB\times RS} + N_\ell^{\rm CMB \times RS}.
\end{align}
Here $\mathbf{N}_\ell$ is the covariance of the foregrounds, instrumental and atmospheric noises. We also note that when we consider the cross-correlation between the primary temperature CMB and the Rayleigh scattering $E$-mode polarization signal (or vice-versa), the noise term $N_\ell^{\rm CMB \times RS}$ vanishes due to our modelling of foregrounds and noise.

The variance of the reconstructed cross-spectrum is given by\citep{Knox_1995},
\begin{align}
   \Delta C_{\ell} &= \frac{1}{\sqrt{f_\mathrm{sky}(2\ell+1)}}\Big[ \left(C_{\ell}^\mathrm{RS}+N_{\ell}^\mathrm{RS}\right)
  \left(C_{\ell}^\mathrm{CMB}+N_{\ell}^\mathrm{CMB}\right)\nonumber \\
  & +\left(C_{\ell}^{\mathrm{CMB} \times \mathrm{RS}} + N_\ell^{\rm CMB \times RS}\right)^2\Big]^{1/2}.
\end{align}
The calculations above pertain to a given experiment. However, we are interested in the combination of experiments to maximize our ability to measure the Rayleigh scattering signal. To obtain a combined noise on the cross-spectra for a given set of experiments, we augment the data vector $\mathbf{d}(p)$ by taking the set of frequency channels to be:
\begin{equation}
    f^\mathrm{total}=(f^1,f^2,\cdots,f^i).
\end{equation}

The signal-to-noise of the cross-correlation is given by:
\begin{equation}
  \mathrm{SNR}^2
  =\sum_\ell\frac{\left(C_{\ell}^{\mathrm{CMB} \times \mathrm{RS}}\right)^2}{\Delta C_{\ell}^2}
\label{knox}
\end{equation}

\section{Results}
\label{sec:Results}

The detectability for different experiments and combinations of experiments are summarized in Tab.~\ref{table}. In Fig.~\ref{fghitsn} we focus on the combination of CCAT-prime, SO and \textit{Planck} and show the SNR for the four CMB-Rayleigh scattering cross-spectra with and without foregrounds. For this combined configuration we estimate SNR of $1.1\sigma$ for the $TT$ signal with foregrounds present. For all other combinations, we find signal-to-noise below $1\sigma$. Without foregrounds, the combined configuration of CCAT-prime, SO and \textit{Planck} could reach $4.6\sigma$ in $TT$, $1.1\sigma$ in $EE$, and $1.6\sigma$ in $ET$. This clearly demonstrates the significant impact of foregrounds; for example the $TT$ SNR is reduced by $77\%$ by the inclusion of foregrounds and their removal. Similarly, the signal-to-noise in $ET$  reduces by $81\%$, $EE$ reduces by $45\%$ and $TE$ reduces by $40\%$. The latter is a result of the fact that in $TE$, foregrounds are not very significant.

\begin{table}[tb]
\begin{tabular}{|l|cccc|}
\hline
\rowcolor[HTML]{E6E8FF}
                                                                       & \multicolumn{4}{c|}{\cellcolor[HTML]{E6E8FF}\textbf{SNR}}                                                                                \\ \hline
\rowcolor[HTML]{E6E8FF}
\textbf{Configuration}                                                 & \multicolumn{1}{c|}{\cellcolor[HTML]{E6E8FF}$TT$} & \multicolumn{1}{c|}{\cellcolor[HTML]{E6E8FF}$TE$} & \multicolumn{1}{c|}{\cellcolor[HTML]{E6E8FF}$ET$} & $EE$ \\ \hline
\cellcolor[HTML]{E6E8FF}CCATp+SO+\textit{Planck}                          & \multicolumn{1}{c|}{1.1}                          & \multicolumn{1}{c|}{0.3}                         & \multicolumn{1}{c|}{0.3}                         & 0.6 \\ \hline
\rowcolor[HTML]{EFEFEF}
\cellcolor[HTML]{E6E8FF}CCATp+SO+\textit{Planck} no Foreground                             & \multicolumn{1}{c|}{\cellcolor[HTML]{EFEFEF}4.6}  & \multicolumn{1}{c|}{\cellcolor[HTML]{EFEFEF}0.5} & \multicolumn{1}{c|}{\cellcolor[HTML]{EFEFEF}1.6}  & 1.1  \\ \hline
\cellcolor[HTML]{E6E8FF}LiteBIRD                                       & \multicolumn{1}{c|}{1.9}                          & \multicolumn{1}{c|}{0.1}                         & \multicolumn{1}{c|}{0.9}                         & 0.2 \\ \hline
\rowcolor[HTML]{EFEFEF}
\cellcolor[HTML]{E6E8FF}LiteBIRD+CCATp+\textit{Planck} & \multicolumn{1}{c|}{\cellcolor[HTML]{EFEFEF}2.2}  & \multicolumn{1}{c|}{\cellcolor[HTML]{EFEFEF}0.2} & \multicolumn{1}{c|}{\cellcolor[HTML]{EFEFEF}0.9} & 0.4 \\ \hline
\cellcolor[HTML]{E6E8FF}PICO                                        & \multicolumn{1}{c|}{85}                          & \multicolumn{1}{c|}{17}                           & \multicolumn{1}{c|}{43}                           & 26   \\ \hline
\rowcolor[HTML]{EFEFEF}
\end{tabular}
\caption{Forecasted signal-to-noise ratio for various configurations. The CCATp+SO+\textit{Planck} configuration can potentially achieve a $1.1\sigma$ detection with foreground present. We also show the impact of foregrounds for the combination CCATp+SO+\textit{Planck}. PICO, with its broad frequency coverage, is capable of a high-significance detection of Rayleigh scattering even after foreground removal.}
\label{table}
\end{table}

The SNR without foregrounds is calculated using a modified covariance matrix with only instrumental and atmospheric noise. As a consistency check we compare our SNR without foregrounds to the Fisher forecasts from \citet{Beringue2020}, which did not include any foregrounds. We find good agreements between the two methods, with the SNR from the cILC being slightly smaller due to the increase in variance caused by the deprojections.

In Fig.~\ref{noisecomp} we show the $2\sigma$ upper limit on the detectability of the $TT$ and $EE$ Rayleigh scattering signals for the CCAT-prime {\it Planck}, SO-LAT, and SO-SAT ($EE$ only) as well as their combination. As expected we find {\it Planck} to dominate the constraints on large scales in temperature as ground-based experiments are heavily impacted by the atmosphere. On small scales and in polarization, the improvements in both sensitivity and resolution enabled by ground-based surveys lead to SO and CCAT-prime driving the constraints.

\begin{figure}[tb]

  \centering

  \includegraphics[width=\linewidth]{s2nfig/s2nfgcspnew.pdf}

\caption{Forecast of SNR with and without foregrounds for CMB Rayleigh scattering cross-spectrum using CCATp+SO+\textit{Planck}. The $TT$ SNR is significantly impacted by the foregrounds, being reduced by $77\%$ when foregrounds are present, which demonstrates that foregrounds are an hindrance for the first detection of the Rayleigh scattering signal.}
\label{fghitsn}
\end{figure}

We also estimate how much CCAT-prime is impacted by the atmosphere. The baseline signal-to-noise including atmosphere for $TT$ signal is $\mathrm{SNR}=0.07$, white it is $\mathrm{SNR}=1.79$ without atmosphere. If we combine the atmosphere-less CCAT-prime with SO and \textit{Planck}, we obtain $\mathrm{SNR}=1.88$ for the $TT$ channel. This suggests that the atmosphere severely limits detectability prospects for ground-based observations, thus mitigating atmospheric contributions plays an important role in the detection of Rayleigh scattering from the ground \citep{SPTinprep}.

\begin{figure*}[ht
\centering
 \subfloat{\includegraphics[width=0.49\linewidth]{noisefig/TTnew.pdf}}
 \hfill
 \subfloat{\includegraphics[width=0.49\textwidth]{noisefig/EEnew.pdf}}

\caption{Forecasted $2\sigma$ upper limits on the Rayleigh scattering cross-spectrum with CCAT-prime, Simons Observatory LAT+SAT, \textit{Planck}, and the combination of all of these telescopes. The SO SAT only enters for the $EE$ signal, since it only targets polarized signals. As expected, the combined configuration has the lowest $2\sigma$ upper limits; \textit{Planck} is most helpful on large scales, and CCAT-prime and SO contribute mainly on small scales.}
\label{noisecomp}
\end{figure*}

\begin{figure*}[!ht]
\centering
 \subfloat{\includegraphics[width=.49\linewidth]{noisefig/fgttnew.pdf}}
 \quad
 \subfloat{\includegraphics[width=.49\linewidth]{noisefig/fgtenew.pdf}}
 \quad
 \subfloat{\includegraphics[width=.49\linewidth]{noisefig/fgetnew.pdf}}
 \quad
 \subfloat{\includegraphics[width=.49\linewidth]{noisefig/fgeenew.pdf}}
\caption{Forecasted $2\sigma$ upper limits on the $TT$ and $EE$ Rayleigh scattering cross-spectrum with and without foregrounds, demonstrating how the foregrounds impact detectability. The SO SAT only enters for the $EE$ signal, since it only targets polarized signals.}
\label{fghit}
\end{figure*}

\subsection{Residual Foreground Bias}
\label{sec:bias}

The cILC foreground cleaning method does not perfectly null all foregrounds, only those explicitly chosen to be nulled. As a result, residual signals from foreground not explicitly nulled remain and could possibly contaminate or bias the detection of Rayleigh scattering. For Rayleigh scattering this is especially important to check due to the relative amplitudes of foregrounds compared to Rayleigh cross-spectra. Knowing the amplitude of residuals will also inform whether further foreground cleaning is required. For this study we focus on temperature, since the foregrounds are less bright in polarization and hence less of a concern, as highlighted in Fig.~\ref{fghitsn}.

We define the bias $B_{\ell}$ induced by residual foreground from a component with frequency-covariance matrix $\mathbf{S}_\ell$ to be:
\begin{equation}
  B_{\ell}= \mathbf{w_{a=0,b=1}}^\mathsf{T}{}  \mathbf{S}_{\ell}   \mathbf{w_{a=1,b=0}},
  \label{eq:bias}
\end{equation}
where $\mathbf{w_{a=0,b=1}}$ and $\mathbf{w_{a=1,b=0}}$ are defined in Eq.~\ref{eq:cilc-weights}. We consider tSZ, CIB, and radio point sources for $\mathbf{S}_\ell$. We are especially interested in the CIB contamination since CIB becomes the dominant foreground signal at sub-millimeter wavelengths where the Rayleigh signal is stronger.

In Fig.~\ref{bias} we show the level of residual bias ($B_{\ell}$) resulting from tSZ, CIB, and radio point sources foregrounds calculated for both CCAT-prime and a combination of CCAT-prime, SO and {\it Planck}. Fig.~\ref{bias} illustrates that for CCAT-prime alone the CIB and tSZ are $1$ to $2$ orders of magnitude stronger than the CMB-Rayleigh cross-spectrum at $280~\mathrm{GHz}$, while radio point source residuals are around $1$ to $4$ orders of magnitude lower than the signal. Between $\ell=2000-3000$, we can see that the tSZ bias is slightly lower than the signal at $280~\mathrm{GHz}$.

Combining the experiments adds more frequency coverage, and the CIB residual is more suppressed, but the tSZ residual is less well suppressed. Given the results found here, ground-based detection of Rayleigh scattering is further complicated by CIB and tSZ residuals after the cILC cleaning. This necessitates the need to find methods to better constrain those components in the future or further constrain them in the cILC at the cost of SNR.

\begin{figure}[!ht]
\centering

    \includegraphics[width=\columnwidth]{noisefig/biasnew.pdf}

\caption{CIB, tSZ, and radio point source bias for CCAT-prime and for the combined configuration. With only CCAT-prime, CIB is $1$ to $2$ orders of magnitude stronger than the Rayleigh scattering cross-spectrum at $280~\mathrm{GHz}$, and radio point source bias is around $1$ to $4$ orders of magnitude lower than the signal. The tSZ bias is comparable to the signal itself. With the combined configuration of CCAT-prime, SO and \textit{Planck}, we the CIB is more suppressed, but tSZ is less suppressed. In general, the tSZ and CIB biases are mostly comparable or even larger than the cross Rayleigh scattering cross-spectrum at $280~\mathrm{GHz}$.}
\label{bias}
\end{figure}

\subsection{Experimental Calibration Tolerance}
\label{sec:gain}

Another potential source of hindrance to a first detection of Rayleigh scattering are instrumental systematic effects. In particular, we study the impact of imperfect calibration of a given experiment. We will refer to this overall calibration systematic as the gain. We model this gain systematic following \citet{Abitbol_2021}. Such gain calibration systematics might arise from fabrication variation in the on-chip bandpass filters between both detectors and wafers, time variation from atmospheric fluctuations, and systematic effects in the bandpass calibrations from measurements using Fourier transform spectrometers or other techniques. A formal model for gain calibration systematics would include an effective shift in frequency $\Delta\nu$, and a gain difference $\Delta g$, which describes how much of the sky emission signal the extra gain picks up \citep{Abitbol_2021}. We estimate an approximately $\Delta g=1\%$ calibration gain of the sky emission and do not consider the possible effective shift in the mean frequency to simplify our analysis.  We can express the gain as
\begin{equation}
  G_{\ell}=\Delta g\;\mathbf{w_{a=0,b=1}}^\mathsf{T}{}  \mathbf{E}_{\ell}   \mathbf{w_{a=1,b=0}},
\end{equation}
where the sky emission $\mathbf{E}=\mathbf{S_{fg}}+\mathbf{S}$. Here $\mathbf{S_{fg}}$ represents the foregrounds, and $\mathbf{S}$ represents the signals.

We show the sky emission in temperature for each channel for CCAT-prime in the right panel of Fig.~\ref{gain}. For the CCAT-prime lower frequency channels the sky emission looks like a CMB auto-spectrum, but at higher frequencies the extragalactic foreground signals are more dominant at high-$\ell$, especially the CIB.

We calculate the gain bias for the combined configuration of CCAT-prime, SO, and \textit{Planck} with $\Delta g=1\%$ and $\Delta g=35\%$. These biases are shown in left panel of Fig.~\ref{gain}. For $\Delta g=1\%$ of the sky emission the gain bias is around $3$ to $4$ orders of magnitude smaller than the noise power spectrum and would therefore be a sufficient calibration of instrument for our purpose. This is, in principle, informative in terms of the order of magnitude of the gain. It is helpful to verify that the gain is generally $4$ to $6$ orders of magnitude lower than our noise power spectrum. If we raise the gain to $\Delta g=35\%$, we see that the cross over happens at around $\ell=3000$. This shows that if we want to set our $\ell_\mathrm{max}$ of our observation higher than $\ell=3000$, a calibration around $1\%$ would be sufficient. A more systematic study is required to establish gain calibration thresholds that contain 2 significant figures for CCAT-prime, SO, and Planck.

\begin{figure*}[!ht]
\centering
 \subfloat{\includegraphics[width=0.49\linewidth]{noisefig/gainperch.pdf}}
 \hfill
 \subfloat{\includegraphics[width=0.49\linewidth]{noisefig/gainfull.pdf}}
\caption{The left panel shows the sky emission on each channel of CCAT-prime. We omitted the $850~\mathrm{GHz}$ channel because its noise is so high that it is given a nearly zero weight by our cILC. We can see that at low frequencies, the sky emission is close to the CMB itself, but in higher frequency channels, the foreground terms start to dominate. The right panel is the gain model for the combined configuration of CCAT-prime, SO, and \textit{Planck}. We assumed $\Delta g=1\%$ gain from sky emission, and we can see that it is roughly on the same scale as the Rayleigh signal. With CCAT-prime channels, in lower frequency channels, the gain looks like CMB itself, but in the higher frequency channels, the gain gradually flattens out due to the dominance of radio point sources. The gain for the combined configuration is approximately $3$ to $4$ orders of magnitude lower than the noise power spectrum, meaning a $1\%$ gain calibration would be sufficient for our purpose. We also showed that a $35\%$ gain calibration, in dashed red line, would be comparable to the noise power spectrum at high $\ell$.}
\label{gain}
\end{figure*}

\subsection{The Detection of Rayleigh Scattering by Satellite Missions}

\label{sec:space}
Ground-based observations are heavily impacted by atmospheric noise. Given the large impact of atmosphere on the detectability of the Rayleigh scattering signal, ideally one would attempt such observations in space where there is no atmosphere and where large scale anisotropies can be efficiently recovered. We calculate the detectability of the Rayleigh scattering signal for the planned satellite mission LiteBIRD~\cite{LiteBird} and the proposed mission PICO~\cite{PICO}.

The results are summarized in Table.~\ref{table}. The small aperture of LiteBIRD restricts the range of scales it can probe to $\ell\lesssim 400$. Including realistic foregrounds, LiteBIRD could reach a $1.9\sigma$ detection in $TT$, which is better than the combined configuration of CCAT-prime, SO, and \textit{Planck}, but insufficient for a significant first detection. PICO would lead to a significant detection even in the presence of foregrounds. With its larger aperture, low noise, and multiple frequency channels that span the millimeter to sub-millimeter wavelength range, we forecast that PICO should achieve a $85\sigma$ detection in $TT$ after foregrounds removal. In Fig.~\ref{sp}, we show the error bars for PICO and LiteBIRD. LiteBIRD has much larger error bars than PICO as expected.

\begin{figure*}[!ht]
\centering
 \subfloat{\includegraphics[width=0.49\linewidth]{noisefig/nlttsp.pdf}}
 \hfill
 \subfloat{\includegraphics[width=0.49\linewidth]{noisefig/nleesp.pdf}}

\caption{Noise comparison for PICO and LiteBIRD. PICO has much tighter error bars, especially on small scales, which will enable a high significance measurement of Rayleigh scattering. LiteBIRD has larger error bars, but can still contribute information on large scales.
}
\label{sp}
\end{figure*}

The fact that PICO is forecast to significantly detect the Rayleigh scattering signal after removing foregrounds, naturally called for a careful study of the biases induced by foregrounds residuals. Following the method outlined in  Sec.~\ref{sec:bias} we calculate the residual foreground bias for PICO. The biases from tSZ, CIB and radio point sources for PICO are shown in Fig.~\ref{picobias}. We find that the CMB-Rayleigh cross-spectrum at $280~\mathrm{GHz}$ is $5$ to $6$ orders of magnitude larger than the bias of tSZ and radio point sources. It is also $8$ to $10$ orders of magnitude larger than the bias of the CIB. It is clear that the multitude of frequency channels and low noise levels not only aid in PICO's ability to detect Rayleigh scattering, but also in its ability to evade residual foreground bias.

The highly significant detection of Rayleigh scattering forecast for PICO additionally leads to the question of calibration uncertainty. Following Sec.~\ref{sec:gain} we calculate the gain bias for a $\Delta g=1\%$. The results are shown in Fig.~\ref{picogain}. At the 1\% threshold the gain bias for PICO is $1$ to $4$ orders of magnitude lower than the reconstructed noise at $\ell < 3 <3000$. This means that a 1\% overall calibration is more than sufficient to safely detect Rayleigh scattering without significant biases. Overall, PICO is robust against the systematics studied here for its forecast detection of the Rayleigh scattering signal.

\begin{figure}[htp]

  \centering

  \includegraphics[width=\columnwidth]{noisefig/picobias.pdf}

  \caption{Bias of tSZ, CIB, and radio point sources for PICO. The bias of tSZ, CIB and radio point sources are all lower than the CMB-Rayleigh cross-spectrum at $280~\mathrm{GHz}$ by several orders of magnitude.  }
  \label{picobias}
\end{figure}

\begin{figure}[htp]

  \centering

  \includegraphics[width=\linewidth]{noisefig/picogain.pdf}

  \caption{Gain of PICO. The gain of $\Delta g=1\%$ of sky emission is around $3$ to $4$ orders of magnitude lower than the noise power spectrum. This means that a $1\%$ gain calibration would be sufficient for PICO.
  }
  \label{picogain}
\end{figure}

Given LiteBIRD's low forecast SNR we do not calculate the residual foreground and gain biases for it. However, LiteBIRD's sensitivity on large scale could still provide important information when combined with other measurements. We consider the combination LiteBird and \textit{Planck} with CCAT-prime.

LiteBIRD contributes no information for $\ell\gtrsim 400$. On large scales $\ell<400$, the reconstructed noise converges to LiteBIRD noise. If we limit the sky fraction to the patch observed by CCAT-prime ($f_\mathrm{sky}=0.4$), we will lose essential information on large scales from LiteBIRD ($f_\mathrm{sky}=0.7$). We calculate the signal-to-noise gained by each telescope in different sky fractions and add them together in quadrature. Essentially, we assign $0$ weight to CCAT-prime in $\ell\leq400$, and $0$ weight to LiteBIRD in $\ell>400$, which is consistent with the cILC method, i.e.,
\begin{align}
     (\mathrm{SNR})^2 & =  \sum_{\ell=2}^{400}\left(\frac{\left(C_{\ell}^{\mathrm{CMB} \times \mathrm{RS}}\right)^2}{\Delta C_{\ell}^2}\right)_\mathrm{LiteBIRD} \nonumber \\
    & + \sum_{\ell=401}^{3000}\left(\frac{\left(C_{\ell}^{\mathrm{CMB} \times \mathrm{RS}}\right)^2}{\Delta C_{\ell}^2}\right)_\mathrm{CCATp+\textit{Planck}}
\end{align}

With the addition of CCAT-prime and \textit{Planck}, the SNR is $2.2$ for the $TT$ signal in the presence of foregrounds, which is $12.5\%$ better than the individual LiteBIRD configuration.

\section{Conclusion}

The observation of CMB Rayleigh scattering has the potential to improve constraints on various cosmological parameters, such as the sum of neutrino masses ($\sum m_{\nu}$), the primordial helium fraction ($Y_\mathrm{He}$), and the light relic density ($N_\mathrm{eff}$)~\cite{Alipour2014,Beringue2020}. However, before one can utilize the full potential of the Rayleigh scattering signal for cosmological parameter inference, a first detection is needed. In this paper, we forecast the detectability of the Rayleigh signal in the presence of foregrounds, considering CCAT-prime, Simons Observatory, \textit{Planck}, PICO and LiteBIRD, and applying a cILC to mitigate various foregrounds contamination. We explored combinations of those experiments to minimize the reconstructed noise and determine the effectiveness of these configurations.

Our analysis shows that in the presence of foregrounds a combination of CCAT-prime, SO and \textit{Planck} can reach an SNR of $1.1$ for the $TT$ CMB-Rayleigh cross-spectrum. As expected, this is below previous estimates which were derived in the absence of foregrounds. This is a clear demonstration that foreground removal, as well as atmospheric contamination, have a large impact on detectability. We compared the capability of these ground-based experiments to future satellite missions. PICO is capable of an $85\sigma$ detection of the $TT$ CMB-Rayleigh cross-spectrum even in the presence of foregrounds. LiteBIRD can potentially achieve an SNR of $1.9$ of the same cross signal. When combined with \textit{Planck} and CCAT-prime, LiteBIRD's SNR is improved to $2.2$. We found that the signal-to-noise of $EE$ cross-spectrum is typically the second highest among all four cross-correlations. While the $EE$  spectrum has the weakest signal strength, it also has the lowest cILC reconstruction noise, since polarized foregrounds are weaker. This also explains why the $EE$ signal-to-noise is least impacted by the presence of foregrounds.

Looking beyond the SNR of Rayleigh scattering, we calculated the biases introduced by tSZ, CIB, and radio point source residuals for the combined configuration of CCAT-prime, SO, and \textit{Planck} and for PICO. For CCAT-prime, SO, and \textit{Planck} the biases from tSZ and CIB have a similar magnitude or are larger than the Rayleigh signal at $\nu=280~\mathrm{GHz}$. Thus a detection of Rayleigh scattering requires further foreground cleaning, for example the deprojection of additional foreground components. Similarly, a denser and broader frequency coverage would help. Indeed for PICO, we found that none of those foregrounds introduce a significant bias which would impact a detection.

Furthermore, we estimate the required instrumental calibration to prevent errors in the gain that may impact the detectability of Rayleigh scattering. We found that an overall gain calibration error of $\Delta g=1\%$ is sufficient for CCAT-prime, SO, and \textit{Planck} and  PICO. Nevertheless, a more comprehensive study of the calibration of CCAT-prime, SO, and PICO is still warranted and practical for upcoming analysis.

With CCAT-prime, SO, and \textit{Planck}, a first detection is unlikely to be achievable by relying solely on constrained-ILC methods. Due to the blind nature of such methods, they do not allow us to leverage the scale dependence of the Rayleigh scattering signal. Indeed, for a fixed cosmology, the Rayleigh scattering signal is straightforward to model, which would be an incentive to explore more parametric approaches such as maximum likelihood algorithm for parametric component separation \citep{stompor_maximum_2009, abylkairov_partially_2021}. More realistic analysis of instrumental calibration requirements will also be necessary and will be the topic of future investigations. Telescopes not considered here, such as SPT \citep[][]{SPT3G}, CMB-S4 \citep[][]{CMBS4} and CMB-HD \citep[][]{CMBHD}, are guaranteed to help in an attempt to make a first detection of CMB Rayleigh scattering.

The detection of the Rayleigh signal is an exciting prospect which will eventually allow us to observe a frequency-dependent last-scattering surface. With advancing telescope technology and new data analysis techniques, there is hope to make a first detection of Rayleigh signal this decade. Eventually a significant detection will allow us to use the cosmological information hidden in the Rayleigh signal to improve parameter constraints beyond those achievable from the primary CMB alone.

\section*{Acknowledgements}
We thank Jens Chluba, Tom Crawford, and Blake Sherwin for their useful comments on our manuscript. NB acknowledges support from NSF grant AST-1910021 and NASA grants 21-ADAP21-0114 and 21-ATP21-0129. PDM acknowledges support from the Netherlands organisation for scientific research (NWO) VIDI grant (dossier 639.042.730). JM is supported by the US~Department of Energy under Grant~\mbox{DE-SC0010129}. BB acknowledges postdoctoral support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 849169). SKC acknowledges support from NSF award AST-2001866.